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
00053
00054
00055
00056
00057
00058
00059
00060 #include <algorithm>
00061
00062 #include "G4CollisionOutput.hh"
00063 #include "G4SystemOfUnits.hh"
00064 #include "G4CascadParticle.hh"
00065 #include "G4ParticleLargerEkin.hh"
00066 #include "G4LorentzConvertor.hh"
00067 #include "G4LorentzRotation.hh"
00068 #include "G4LorentzVector.hh"
00069 #include "G4ReactionProductVector.hh"
00070 #include "G4ReactionProduct.hh"
00071 #include "G4ThreeVector.hh"
00072
00073 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
00074 typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
00075
00076
00077 G4CollisionOutput::G4CollisionOutput()
00078 : verboseLevel(0), eex_rest(0), on_shell(false) {
00079 if (verboseLevel > 1)
00080 G4cout << " >>> G4CollisionOutput::G4CollisionOutput" << G4endl;
00081 }
00082
00083
00084 G4CollisionOutput& G4CollisionOutput::operator=(const G4CollisionOutput& right)
00085 {
00086 if (this != &right) {
00087 verboseLevel = right.verboseLevel;
00088 outgoingParticles = right.outgoingParticles;
00089 outgoingNuclei = right.outgoingNuclei;
00090 theRecoilFragment = right.theRecoilFragment;
00091 eex_rest = right.eex_rest;
00092 on_shell = right.on_shell;
00093 }
00094 return *this;
00095 }
00096
00097 void G4CollisionOutput::reset() {
00098 outgoingNuclei.clear();
00099 outgoingParticles.clear();
00100 removeRecoilFragment();
00101 }
00102
00103
00104
00105
00106 void G4CollisionOutput::add(const G4CollisionOutput& right) {
00107 addOutgoingParticles(right.outgoingParticles);
00108 addOutgoingNuclei(right.outgoingNuclei);
00109 theRecoilFragment = right.theRecoilFragment;
00110 }
00111
00112
00113
00114
00115 void G4CollisionOutput::addOutgoingParticles(const std::vector<G4InuclElementaryParticle>& particles) {
00116 outgoingParticles.insert(outgoingParticles.end(),
00117 particles.begin(), particles.end());
00118 }
00119
00120 void G4CollisionOutput::addOutgoingNuclei(const std::vector<G4InuclNuclei>& nuclea) {
00121 outgoingNuclei.insert(outgoingNuclei.end(), nuclea.begin(), nuclea.end());
00122 }
00123
00124
00125
00126 void G4CollisionOutput::addOutgoingParticle(const G4CascadParticle& cparticle) {
00127 addOutgoingParticle(cparticle.getParticle());
00128 }
00129
00130 void G4CollisionOutput::addOutgoingParticles(const std::vector<G4CascadParticle>& cparticles) {
00131 for (unsigned i=0; i<cparticles.size(); i++)
00132 addOutgoingParticle(cparticles[i]);
00133 }
00134
00135
00136
00137 void G4CollisionOutput::addOutgoingParticles(const G4ReactionProductVector* rproducts) {
00138 if (!rproducts) return;
00139
00140 if (verboseLevel) {
00141 G4cout << " >>> G4CollisionOutput::addOutgoingParticles(G4RPVector)"
00142 << G4endl;
00143 }
00144
00145 G4ReactionProductVector::const_iterator j;
00146 for (j=rproducts->begin(); j!=rproducts->end(); ++j) {
00147 G4ParticleDefinition* pd = (*j)->GetDefinition();
00148 G4int type = G4InuclElementaryParticle::type(pd);
00149
00150
00151 G4LorentzVector mom((*j)->GetMomentum(), (*j)->GetTotalEnergy());
00152 mom /= GeV;
00153
00154 if (verboseLevel>1)
00155 G4cout << " Processing " << pd->GetParticleName() << " (" << type
00156 << "), momentum " << mom << " GeV" << G4endl;
00157
00158
00159
00160 if (type) {
00161 outgoingParticles.resize(numberOfOutgoingParticles()+1);
00162 outgoingParticles.back().fill(mom, pd, G4InuclParticle::PreCompound);
00163
00164 if (verboseLevel>1) G4cout << outgoingParticles.back() << G4endl;
00165 } else {
00166 outgoingNuclei.resize(numberOfOutgoingNuclei()+1);
00167 outgoingNuclei.back().fill(mom,pd->GetAtomicMass(),pd->GetAtomicNumber(),
00168 0.,G4InuclParticle::PreCompound);
00169
00170 if (verboseLevel>1) G4cout << outgoingNuclei.back() << G4endl;
00171 }
00172 }
00173 }
00174
00175
00176
00177
00178 void G4CollisionOutput::removeOutgoingParticle(G4int index) {
00179 if (index >= 0 && index < numberOfOutgoingParticles())
00180 outgoingParticles.erase(outgoingParticles.begin()+(size_t)index);
00181 }
00182
00183 void G4CollisionOutput::removeOutgoingNucleus(G4int index) {
00184 if (index >= 0 && index < numberOfOutgoingNuclei())
00185 outgoingNuclei.erase(outgoingNuclei.begin()+(size_t)index);
00186 }
00187
00188
00189
00190 void G4CollisionOutput::removeOutgoingParticle(const G4InuclElementaryParticle& particle) {
00191 particleIterator pos =
00192 std::find(outgoingParticles.begin(), outgoingParticles.end(), particle);
00193 if (pos != outgoingParticles.end()) outgoingParticles.erase(pos);
00194 }
00195
00196 void G4CollisionOutput::removeOutgoingNucleus(const G4InuclNuclei& nuclei) {
00197 nucleiIterator pos =
00198 std::find(outgoingNuclei.begin(), outgoingNuclei.end(), nuclei);
00199 if (pos != outgoingNuclei.end()) outgoingNuclei.erase(pos);
00200 }
00201
00202
00203
00204 void G4CollisionOutput::removeRecoilFragment() {
00205 static const G4Fragment emptyFragment;
00206 theRecoilFragment = emptyFragment;
00207 }
00208
00209
00210
00211 G4LorentzVector G4CollisionOutput::getTotalOutputMomentum() const {
00212 if (verboseLevel > 1)
00213 G4cout << " >>> G4CollisionOutput::getTotalOutputMomentum" << G4endl;
00214
00215 G4LorentzVector tot_mom;
00216 G4int i(0);
00217 for(i=0; i < G4int(outgoingParticles.size()); i++) {
00218 tot_mom += outgoingParticles[i].getMomentum();
00219 }
00220 for(i=0; i < G4int(outgoingNuclei.size()); i++) {
00221 tot_mom += outgoingNuclei[i].getMomentum();
00222 }
00223 tot_mom += theRecoilFragment.GetMomentum()/GeV;
00224
00225 return tot_mom;
00226 }
00227
00228 G4int G4CollisionOutput::getTotalCharge() const {
00229 if (verboseLevel > 1)
00230 G4cout << " >>> G4CollisionOutput::getTotalCharge" << G4endl;
00231
00232 G4int charge = 0;
00233 G4int i(0);
00234 for(i=0; i < G4int(outgoingParticles.size()); i++) {
00235 charge += G4int(outgoingParticles[i].getCharge());
00236 }
00237 for(i=0; i < G4int(outgoingNuclei.size()); i++) {
00238 charge += G4int(outgoingNuclei[i].getCharge());
00239 }
00240 charge += theRecoilFragment.GetZ_asInt();
00241
00242 return charge;
00243 }
00244
00245 G4int G4CollisionOutput::getTotalBaryonNumber() const {
00246 if (verboseLevel > 1)
00247 G4cout << " >>> G4CollisionOutput::getTotalBaryonNumber" << G4endl;
00248
00249 G4int baryon = 0;
00250 G4int i(0);
00251 for(i=0; i < G4int(outgoingParticles.size()); i++) {
00252 baryon += outgoingParticles[i].baryon();
00253 }
00254 for(i=0; i < G4int(outgoingNuclei.size()); i++) {
00255 baryon += G4int(outgoingNuclei[i].getA());
00256 }
00257 baryon += theRecoilFragment.GetA_asInt();
00258
00259 return baryon;
00260 }
00261
00262 G4int G4CollisionOutput::getTotalStrangeness() const {
00263 if (verboseLevel > 1)
00264 G4cout << " >>> G4CollisionOutput::getTotalStrangeness" << G4endl;
00265
00266 G4int strange = 0;
00267 G4int i(0);
00268 for(i=0; i < G4int(outgoingParticles.size()); i++) {
00269 strange += outgoingParticles[i].getStrangeness();
00270 }
00271
00272 return strange;
00273 }
00274
00275
00276 void G4CollisionOutput::printCollisionOutput(std::ostream& os) const {
00277 os << " Output: " << G4endl
00278 << " Outgoing Particles: " << outgoingParticles.size() << G4endl;
00279
00280 G4int i(0);
00281 for(i=0; i < G4int(outgoingParticles.size()); i++)
00282 os << outgoingParticles[i] << G4endl;
00283
00284 os << " Outgoing Nuclei: " << outgoingNuclei.size() << G4endl;
00285 for(i=0; i < G4int(outgoingNuclei.size()); i++)
00286 os << outgoingNuclei[i] << G4endl;
00287
00288 if (theRecoilFragment.GetA() > 0) os << theRecoilFragment << G4endl;
00289 }
00290
00291
00292
00293
00294 void G4CollisionOutput::boostToLabFrame(const G4LorentzConvertor& convertor) {
00295 if (verboseLevel > 1)
00296 G4cout << " >>> G4CollisionOutput::boostToLabFrame" << G4endl;
00297
00298 if (!outgoingParticles.empty()) {
00299 particleIterator ipart = outgoingParticles.begin();
00300 for(; ipart != outgoingParticles.end(); ipart++) {
00301 ipart->setMomentum(boostToLabFrame(ipart->getMomentum(), convertor));
00302 }
00303
00304 std::sort(outgoingParticles.begin(), outgoingParticles.end(),
00305 G4ParticleLargerEkin());
00306 }
00307
00308 if (!outgoingNuclei.empty()) {
00309 nucleiIterator inuc = outgoingNuclei.begin();
00310 for (; inuc != outgoingNuclei.end(); inuc++) {
00311 inuc->setMomentum(boostToLabFrame(inuc->getMomentum(), convertor));
00312 }
00313 }
00314
00315 if (theRecoilFragment.GetA() > 0.) {
00316 theRecoilFragment.SetMomentum(boostToLabFrame(theRecoilFragment.GetMomentum(), convertor));
00317 }
00318 }
00319
00320
00321 G4LorentzVector
00322 G4CollisionOutput::boostToLabFrame(G4LorentzVector mom,
00323 const G4LorentzConvertor& convertor) const {
00324 if (convertor.reflectionNeeded()) mom.setZ(-mom.z());
00325 mom = convertor.rotate(mom);
00326 mom = convertor.backToTheLab(mom);
00327
00328 return mom;
00329 }
00330
00331
00332
00333 void G4CollisionOutput::rotateEvent(const G4LorentzRotation& rotate) {
00334 if (verboseLevel > 1)
00335 G4cout << " >>> G4CollisionOutput::rotateEvent" << G4endl;
00336
00337 particleIterator ipart = outgoingParticles.begin();
00338 for(; ipart != outgoingParticles.end(); ipart++)
00339 ipart->setMomentum(ipart->getMomentum()*=rotate);
00340
00341 nucleiIterator inuc = outgoingNuclei.begin();
00342 for (; inuc != outgoingNuclei.end(); inuc++)
00343 inuc->setMomentum(inuc->getMomentum()*=rotate);
00344
00345 if (theRecoilFragment.GetA() > 0.) {
00346 G4LorentzVector mom = theRecoilFragment.GetMomentum();
00347 theRecoilFragment.SetMomentum(mom*=rotate);
00348 }
00349
00350 }
00351
00352
00353 void G4CollisionOutput::trivialise(G4InuclParticle* bullet,
00354 G4InuclParticle* target) {
00355 if (verboseLevel > 1)
00356 G4cout << " >>> G4CollisionOutput::trivialize" << G4endl;
00357
00358 reset();
00359
00360 if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
00361 outgoingNuclei.push_back(*nuclei_target);
00362 } else {
00363 G4InuclElementaryParticle* particle =
00364 dynamic_cast<G4InuclElementaryParticle*>(target);
00365 outgoingParticles.push_back(*particle);
00366 }
00367
00368 if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
00369 outgoingNuclei.push_back(*nuclei_bullet);
00370 } else {
00371 G4InuclElementaryParticle* particle =
00372 dynamic_cast<G4InuclElementaryParticle*>(bullet);
00373 outgoingParticles.push_back(*particle);
00374 }
00375 }
00376
00377
00378 void G4CollisionOutput::setOnShell(G4InuclParticle* bullet,
00379 G4InuclParticle* target) {
00380 if (verboseLevel > 1)
00381 G4cout << " >>> G4CollisionOutput::setOnShell" << G4endl;
00382
00383 const G4double accuracy = 0.00001;
00384
00385 on_shell = false;
00386
00387 G4LorentzVector ini_mom = bullet->getMomentum();
00388 G4LorentzVector momt = target->getMomentum();
00389 G4LorentzVector out_mom = getTotalOutputMomentum();
00390 if(verboseLevel > 2){
00391 G4cout << " bullet momentum = " << ini_mom.e() <<", "<< ini_mom.x() <<", "<< ini_mom.y()<<", "<< ini_mom.z()<<G4endl;
00392 G4cout << " target momentum = " << momt.e()<<", "<< momt.x()<<", "<< momt.y()<<", "<< momt.z()<<G4endl;
00393 G4cout << " Fstate momentum = " << out_mom.e()<<", "<< out_mom.x()<<", "<< out_mom.y()<<", "<< out_mom.z()<<G4endl;
00394 }
00395
00396 ini_mom += momt;
00397 G4LorentzVector mom_non_cons = ini_mom - out_mom;
00398
00399 G4double pnc = mom_non_cons.rho();
00400 G4double enc = mom_non_cons.e();
00401
00402 setRemainingExitationEnergy();
00403
00404 if(verboseLevel > 2) {
00405 printCollisionOutput();
00406 G4cout << " momentum non conservation: " << G4endl
00407 << " e " << enc << " p " << pnc << G4endl
00408 << " remaining exitation " << eex_rest << G4endl;
00409 }
00410
00411 if (std::fabs(enc) <= accuracy && pnc <= accuracy) {
00412 on_shell = true;
00413 return;
00414 }
00415
00416
00417
00418
00419 if (verboseLevel > 2) G4cout << " re-balancing four-momenta" << G4endl;
00420
00421 G4int npart = outgoingParticles.size();
00422 G4int nnuc = outgoingNuclei.size();
00423
00424 if (npart > 0) {
00425 for (G4int ip=npart-1; ip>=0; ip--) {
00426 if (outgoingParticles[ip].getKineticEnergy()+enc > 0.) {
00427 G4LorentzVector last_mom = outgoingParticles[ip].getMomentum();
00428 last_mom += mom_non_cons;
00429 outgoingParticles[ip].setMomentum(last_mom);
00430 break;
00431 }
00432 }
00433 } else if (nnuc > 0) {
00434 for (G4int in=nnuc-1; in>=0; in--) {
00435 if (outgoingNuclei[in].getKineticEnergy()+enc > 0.) {
00436 G4LorentzVector last_mom = outgoingNuclei[in].getMomentum();
00437 last_mom += mom_non_cons;
00438 outgoingNuclei[in].setMomentum(last_mom);
00439 break;
00440 }
00441 }
00442 } else if (theRecoilFragment.GetA() > 0.) {
00443
00444 G4LorentzVector last_mom = theRecoilFragment.GetMomentum()/GeV;
00445 if ((last_mom.e()-last_mom.m())+enc > 0.) {
00446 last_mom += mom_non_cons;
00447 theRecoilFragment.SetMomentum(last_mom*GeV);
00448 }
00449 }
00450
00451
00452 out_mom = getTotalOutputMomentum();
00453 mom_non_cons = ini_mom - out_mom;
00454 pnc = mom_non_cons.rho();
00455 enc = mom_non_cons.e();
00456
00457 if(verboseLevel > 2){
00458 printCollisionOutput();
00459 G4cout << " momentum non conservation after (1): " << G4endl
00460 << " e " << enc << " p " << pnc << G4endl;
00461 }
00462
00463
00464 G4bool need_hard_tuning = true;
00465
00466 G4double encMeV = mom_non_cons.e() / GeV;
00467 if (theRecoilFragment.GetA() > 0.) {
00468 G4double eex = theRecoilFragment.GetExcitationEnergy();
00469 if (eex > 0.0 && eex + encMeV >= 0.0) {
00470
00471
00472
00473 G4LorentzVector fragMom = theRecoilFragment.GetMomentum();
00474 G4double newMass = fragMom.m() + encMeV;
00475 fragMom.setVectM(fragMom.vect(), newMass);
00476 need_hard_tuning = false;
00477 }
00478 } else if (outgoingNuclei.size() > 0) {
00479 for (G4int i=0; i < G4int(outgoingNuclei.size()); i++) {
00480 G4double eex = outgoingNuclei[i].getExitationEnergy();
00481
00482 if(eex > 0.0 && eex + encMeV >= 0.0) {
00483 outgoingNuclei[i].setExitationEnergy(eex+encMeV);
00484 need_hard_tuning = false;
00485 break;
00486 }
00487 }
00488 if (need_hard_tuning && encMeV > 0.) {
00489 outgoingNuclei[0].setExitationEnergy(encMeV);
00490 need_hard_tuning = false;
00491 }
00492 }
00493
00494 if (!need_hard_tuning) {
00495 on_shell = true;
00496 return;
00497 }
00498
00499
00500 if (verboseLevel > 2)
00501 G4cout << " trying hard (particle-pair) tuning" << G4endl;
00502
00503 std::pair<std::pair<G4int, G4int>, G4int> tune_par = selectPairToTune(mom_non_cons.e());
00504 std::pair<G4int, G4int> tune_particles = tune_par.first;
00505 G4int mom_ind = tune_par.second;
00506
00507 if(verboseLevel > 2) {
00508 G4cout << " p1 " << tune_particles.first << " p2 " << tune_particles.second
00509 << " ind " << mom_ind << G4endl;
00510 }
00511
00512 G4bool tuning_possible =
00513 (tune_particles.first >= 0 && tune_particles.second >= 0 &&
00514 mom_ind >= G4LorentzVector::X);
00515
00516 if (!tuning_possible) {
00517 if (verboseLevel > 2) G4cout << " tuning impossible " << G4endl;
00518 return;
00519 }
00520
00521 G4LorentzVector mom1 = outgoingParticles[tune_particles.first].getMomentum();
00522 G4LorentzVector mom2 = outgoingParticles[tune_particles.second].getMomentum();
00523 G4double newE12 = mom1.e() + mom2.e() + mom_non_cons.e();
00524 G4double R = 0.5 * (newE12 * newE12 + mom2.e() * mom2.e() - mom1.e() * mom1.e()) / newE12;
00525 G4double Q = -(mom1[mom_ind] + mom2[mom_ind]) / newE12;
00526 G4double UDQ = 1.0 / (Q * Q - 1.0);
00527 G4double W = (R * Q + mom2[mom_ind]) * UDQ;
00528 G4double V = (mom2.e() * mom2.e() - R * R) * UDQ;
00529 G4double DET = W * W + V;
00530
00531 if (DET < 0.0) {
00532 if (verboseLevel > 2) G4cout << " DET < 0 " << G4endl;
00533 return;
00534 }
00535
00536
00537 G4double x1 = -(W + std::sqrt(DET));
00538 G4double x2 = -(W - std::sqrt(DET));
00539
00540
00541 G4bool xset = false;
00542 G4double x = 0.0;
00543
00544 if(mom_non_cons.e() > 0.0) {
00545 if(x1 > 0.0) {
00546 if(R + Q * x1 >= 0.0) {
00547 x = x1;
00548 xset = true;
00549 };
00550 };
00551 if(!xset && x2 > 0.0) {
00552 if(R + Q * x2 >= 0.0) {
00553 x = x2;
00554 xset = true;
00555 };
00556 };
00557 } else {
00558 if(x1 < 0.0) {
00559 if(R + Q * x1 >= 0.) {
00560 x = x1;
00561 xset = true;
00562 };
00563 };
00564 if(!xset && x2 < 0.0) {
00565 if(R + Q * x2 >= 0.0) {
00566 x = x2;
00567 xset = true;
00568 };
00569 };
00570 }
00571
00572 if(!xset) {
00573 if(verboseLevel > 2)
00574 G4cout << " no appropriate solution found " << G4endl;
00575 return;
00576 }
00577
00578
00579 mom1[mom_ind] += x;
00580 mom2[mom_ind] -= x;
00581 outgoingParticles[tune_particles.first ].setMomentum(mom1);
00582 outgoingParticles[tune_particles.second].setMomentum(mom2);
00583 out_mom = getTotalOutputMomentum();
00584 std::sort(outgoingParticles.begin(), outgoingParticles.end(), G4ParticleLargerEkin());
00585 mom_non_cons = ini_mom - out_mom;
00586 pnc = mom_non_cons.rho();
00587 enc = mom_non_cons.e();
00588
00589 on_shell = (std::fabs(enc) < accuracy || pnc < accuracy);
00590
00591 if(verboseLevel > 2) {
00592 G4cout << " momentum non conservation tuning: " << G4endl
00593 << " e " << enc << " p " << pnc
00594 << (on_shell?" success":" FAILURE") << G4endl;
00595 }
00596 }
00597
00598
00599
00600
00601 void G4CollisionOutput::setRemainingExitationEnergy() {
00602 eex_rest = theRecoilFragment.GetExcitationEnergy() / GeV;
00603
00604 for(G4int i=0; i < G4int(outgoingNuclei.size()); i++)
00605 eex_rest += outgoingNuclei[i].getExitationEnergyInGeV();
00606 }
00607
00608
00609 std::pair<std::pair<G4int, G4int>, G4int>
00610 G4CollisionOutput::selectPairToTune(G4double de) const {
00611 if (verboseLevel > 2)
00612 G4cout << " >>> G4CollisionOutput::selectPairToTune" << G4endl;
00613
00614 std::pair<G4int, G4int> tup(-1, -1);
00615 G4int i3 = -1;
00616 std::pair<std::pair<G4int, G4int>, G4int> tuner(tup, i3);
00617
00618 if (outgoingParticles.size() < 2) return tuner;
00619
00620 G4int ibest1 = -1;
00621 G4int ibest2 = -1;
00622 G4double pbest = 0.0;
00623 G4double pcut = 0.3 * std::sqrt(1.88 * std::fabs(de));
00624 G4double p1 = 0.0;
00625
00626 for (G4int i = 0; i < G4int(outgoingParticles.size())-1; i++) {
00627 G4LorentzVector mom1 = outgoingParticles[i].getMomentum();
00628
00629 for (G4int j = i+1; j < G4int(outgoingParticles.size()); j++) {
00630 G4LorentzVector mom2 = outgoingParticles[j].getMomentum();
00631
00632 for (G4int l = G4LorentzVector::X; l<=G4LorentzVector::Z; l++) {
00633 if (mom1[l]*mom2[l]<0.0) {
00634 if (std::fabs(mom1[l])>pcut && std::fabs(mom2[l])>pcut) {
00635 G4double psum = std::fabs(mom1[l]) + std::fabs(mom2[l]);
00636
00637 if(psum > pbest) {
00638 ibest1 = i;
00639 ibest2 = j;
00640 i3 = l;
00641 p1 = mom1[l];
00642 pbest = psum;
00643 }
00644 }
00645 }
00646 }
00647 }
00648 }
00649
00650 if (i3 < 0) return tuner;
00651
00652 tuner.second = i3;
00653
00654
00655 if (de > 0.0) {
00656 tuner.first.first = (p1>0.) ? ibest1 : ibest2;
00657 tuner.first.second = (p1>0.) ? ibest2 : ibest1;
00658 } else {
00659 tuner.first.first = (p1<0.) ? ibest2 : ibest1;
00660 tuner.first.second = (p1<0.) ? ibest1 : ibest2;
00661 }
00662
00663 return tuner;
00664 }