00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 #include "G4CascadeCheckBalance.hh"
00053 #include "globals.hh"
00054 #include "G4CascadParticle.hh"
00055 #include "G4InuclElementaryParticle.hh"
00056 #include "G4InuclNuclei.hh"
00057 #include "G4InuclParticle.hh"
00058 #include "G4CollisionOutput.hh"
00059 #include "G4LorentzVector.hh"
00060 #include <vector>
00061
00062
00063
00064
00065 const G4double G4CascadeCheckBalance::tolerance = 1e-6;
00066
00067 G4CascadeCheckBalance::G4CascadeCheckBalance(const char* owner)
00068 : G4VCascadeCollider(owner), relativeLimit(G4CascadeCheckBalance::tolerance),
00069 absoluteLimit(G4CascadeCheckBalance::tolerance), initialBaryon(0),
00070 finalBaryon(0), initialCharge(0), finalCharge(0), initialStrange(0),
00071 finalStrange(0) {}
00072
00073 G4CascadeCheckBalance::G4CascadeCheckBalance(G4double relative,
00074 G4double absolute,
00075 const char* owner)
00076 : G4VCascadeCollider(owner), relativeLimit(relative),
00077 absoluteLimit(absolute), initialBaryon(0), finalBaryon(0),
00078 initialCharge(0), finalCharge(0), initialStrange(0),
00079 finalStrange(0) {}
00080
00081
00082
00083
00084 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet,
00085 G4InuclParticle* target,
00086 G4CollisionOutput& output) {
00087 if (verboseLevel)
00088 G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide"
00089 << G4endl;
00090
00091 initial *= 0.;
00092 if (bullet) initial += bullet->getMomentum();
00093 if (target) initial += target->getMomentum();
00094
00095
00096 initialCharge = 0;
00097 if (bullet) initialCharge += G4int(bullet->getCharge());
00098 if (target) initialCharge += G4int(target->getCharge());
00099
00100 G4InuclElementaryParticle* pbullet =
00101 dynamic_cast<G4InuclElementaryParticle*>(bullet);
00102 G4InuclElementaryParticle* ptarget =
00103 dynamic_cast<G4InuclElementaryParticle*>(target);
00104
00105 G4InuclNuclei* nbullet = dynamic_cast<G4InuclNuclei*>(bullet);
00106 G4InuclNuclei* ntarget = dynamic_cast<G4InuclNuclei*>(target);
00107
00108 initialBaryon =
00109 ((pbullet ? pbullet->baryon() : nbullet ? nbullet->getA() : 0) +
00110 (ptarget ? ptarget->baryon() : ntarget ? ntarget->getA() : 0) );
00111
00112
00113 initialStrange = 0;
00114 if (pbullet) initialStrange += pbullet->getStrangeness();
00115 if (ptarget) initialStrange += ptarget->getStrangeness();
00116
00117
00118 final = output.getTotalOutputMomentum();
00119 finalBaryon = output.getTotalBaryonNumber();
00120 finalCharge = output.getTotalCharge();
00121 finalStrange = output.getTotalStrangeness();
00122
00123
00124 if (verboseLevel) {
00125 G4cout << " initial px " << initial.px() << " py " << initial.py()
00126 << " pz " << initial.pz() << " E " << initial.e()
00127 << " baryon " << initialBaryon << " charge " << initialCharge
00128 << " strange " << initialStrange << G4endl
00129 << " final px " << final.px() << " py " << final.py()
00130 << " pz " << final.pz() << " E " << final.e()
00131 << " baryon " << finalBaryon << " charge " << finalCharge
00132 << " strange " << finalStrange << G4endl;
00133 }
00134 }
00135
00136
00137
00138 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet,
00139 G4InuclParticle* target,
00140 const std::vector<G4InuclElementaryParticle>& particles) {
00141 if (verboseLevel)
00142 G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<vector>)"
00143 << G4endl;
00144
00145 tempOutput.reset();
00146 tempOutput.addOutgoingParticles(particles);
00147 collide(bullet, target, tempOutput);
00148 }
00149
00150
00151
00152
00153 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet,
00154 G4InuclParticle* target,
00155 const std::vector<G4InuclNuclei>& fragments) {
00156 if (verboseLevel)
00157 G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<vector>)"
00158 << G4endl;
00159
00160 tempOutput.reset();
00161 tempOutput.addOutgoingNuclei(fragments);
00162 collide(bullet, target, tempOutput);
00163 }
00164
00165
00166
00167
00168 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet,
00169 G4InuclParticle* target,
00170 const std::vector<G4CascadParticle>& particles) {
00171 if (verboseLevel)
00172 G4cout << " >>> G4CascadeCheckBalance(" << theName
00173 << ")::collide(<cparticles>)" << G4endl;
00174
00175 tempOutput.reset();
00176 tempOutput.addOutgoingParticles(particles);
00177 collide(bullet, target, tempOutput);
00178 }
00179
00180
00181
00182
00183 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet,
00184 G4InuclParticle* target,
00185 G4CollisionOutput& output,
00186 const std::vector<G4CascadParticle>& cparticles) {
00187 if (verboseLevel)
00188 G4cout << " >>> G4CascadeCheckBalance(" << theName
00189 << ")::collide(<EP>,<CP>)" << G4endl;
00190
00191 tempOutput.reset();
00192 tempOutput.add(output);
00193 tempOutput.addOutgoingParticles(cparticles);
00194 collide(bullet, target, tempOutput);
00195 }
00196
00197
00198
00199
00200 G4bool G4CascadeCheckBalance::energyOkay() const {
00201 G4bool relokay = (std::abs(relativeE()) < relativeLimit);
00202 G4bool absokay = (std::abs(deltaE()) < absoluteLimit);
00203
00204 if (verboseLevel && !(relokay || absokay)) {
00205 G4cerr << theName << ": Energy conservation: relative " << relativeE()
00206 << (relokay ? " conserved" : " VIOLATED")
00207 << " absolute " << deltaE()
00208 << (absokay ? " conserved" : " VIOLATED") << G4endl;
00209 } else if (verboseLevel > 1) {
00210 G4cout << theName << ": Energy conservation: relative " << relativeE()
00211 << " conserved absolute " << deltaE() << " conserved" << G4endl;
00212 }
00213
00214 return (relokay && absokay);
00215 }
00216
00217 G4bool G4CascadeCheckBalance::ekinOkay() const {
00218 G4bool relokay = (std::abs(relativeKE()) < relativeLimit);
00219 G4bool absokay = (std::abs(deltaKE()) < absoluteLimit);
00220
00221 if (verboseLevel && !(relokay || absokay)) {
00222 G4cerr << theName << ": Kinetic energy balance: relative "
00223 << relativeKE() << (relokay ? " conserved" : " VIOLATED")
00224 << " absolute " << deltaKE()
00225 << (absokay ? " conserved" : " VIOLATED") << G4endl;
00226 } else if (verboseLevel > 1) {
00227 G4cout << theName << ": Kinetic energy balance: relative "
00228 << relativeKE() << " conserved absolute " << deltaKE()
00229 << " conserved" << G4endl;
00230 }
00231
00232 return (relokay && absokay);
00233 }
00234
00235 G4bool G4CascadeCheckBalance::momentumOkay() const {
00236 G4bool relokay = (std::abs(relativeP()) < relativeLimit);
00237 G4bool absokay = (std::abs(deltaP()) < absoluteLimit);
00238
00239 if (verboseLevel && !(relokay || absokay)) {
00240 G4cerr << theName << ": Momentum conservation: relative " << relativeP()
00241 << (relokay ? " conserved" : " VIOLATED")
00242 << " absolute " << deltaP()
00243 << (absokay ? " conserved" : " VIOLATED") << G4endl;
00244 } else if (verboseLevel > 1) {
00245 G4cout << theName << ": Momentum conservation: relative " << relativeP()
00246 << " conserved absolute " << deltaP() << " conserved" << G4endl;
00247 }
00248
00249 return (relokay && absokay);
00250 }
00251
00252 G4bool G4CascadeCheckBalance::baryonOkay() const {
00253 G4bool bokay = (deltaB() == 0);
00254
00255 if (verboseLevel && !bokay)
00256 G4cerr << theName << ": Baryon number VIOLATED " << deltaB() << G4endl;
00257
00258 return bokay;
00259 }
00260
00261 G4bool G4CascadeCheckBalance::chargeOkay() const {
00262 G4bool qokay = (deltaQ() == 0);
00263
00264 if (verboseLevel && !qokay)
00265 G4cerr << theName << ": Charge conservation VIOLATED " << deltaQ()
00266 << G4endl;
00267
00268 return qokay;
00269 }
00270
00271 G4bool G4CascadeCheckBalance::strangeOkay() const {
00272 G4bool sokay = (deltaS() == 0);
00273
00274 if (verboseLevel && !sokay)
00275 G4cerr << theName << ": Strangeness conservation VIOLATED " << deltaS()
00276 << G4endl;
00277
00278 return sokay;
00279 }