83 : verboseLevel(0), eex_rest(0), on_shell(false) {
85 G4cout <<
" >>> G4CollisionOutput::G4CollisionOutput" <<
G4endl;
134 particles.begin(), particles.end());
148 for (
unsigned i=0; i<cparticles.size(); i++)
155 if (!rproducts)
return;
158 G4cout <<
" >>> G4CollisionOutput::addOutgoingParticles(G4RPVector)"
162 G4ReactionProductVector::const_iterator j;
163 for (j=rproducts->begin(); j!=rproducts->end(); ++j) {
173 <<
"), momentum " << mom <<
" GeV" <<
G4endl;
232 G4cout <<
" >>> G4CollisionOutput::getTotalOutputMomentum" <<
G4endl;
251 G4cout <<
" >>> G4CollisionOutput::getTotalCharge" <<
G4endl;
270 G4cout <<
" >>> G4CollisionOutput::getTotalBaryonNumber" <<
G4endl;
289 G4cout <<
" >>> G4CollisionOutput::getTotalStrangeness" <<
G4endl;
302 os <<
" Output: " <<
G4endl
321 G4cout <<
" >>> G4CollisionOutput::boostToLabFrame" <<
G4endl;
340 fmom = ifrag->GetMomentum() /
GeV;
350 mom = convertor.
rotate(mom);
360 G4cout <<
" >>> G4CollisionOutput::rotateEvent" <<
G4endl;
364 ipart->setMomentum(ipart->getMomentum()*=rotate);
368 inuc->setMomentum(inuc->getMomentum()*=rotate);
373 ifrag->SetMomentum(mom*=rotate);
381 G4cout <<
" >>> G4CollisionOutput::trivialize" <<
G4endl;
385 if (
G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
393 if (
G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
406 G4cout <<
" >>> G4CollisionOutput::setOnShell" <<
G4endl;
416 G4cout <<
" bullet momentum = " << ini_mom.e() <<
", "<< ini_mom.x() <<
", "<< ini_mom.y()<<
", "<< ini_mom.z()<<
G4endl;
417 G4cout <<
" target momentum = " << momt.e()<<
", "<< momt.x()<<
", "<< momt.y()<<
", "<< momt.z()<<
G4endl;
418 G4cout <<
" Fstate momentum = " << out_mom.e()<<
", "<< out_mom.x()<<
", "<< out_mom.y()<<
", "<< out_mom.z()<<
G4endl;
431 G4cout <<
" momentum non conservation: " << G4endl
432 <<
" e " << enc <<
" p " << pnc << G4endl
436 if (std::fabs(enc) <= accuracy && pnc <= accuracy) {
453 for (
G4int ip=npart-1; ip>=0; ip--) {
461 }
else if (nnuc > 0) {
462 for (
G4int in=nnuc-1; in>=0; in--) {
470 }
else if (nfrag > 0) {
471 for (
G4int ifr=nfrag-1; ifr>=0; ifr--) {
474 if ((last_mom.e()-last_mom.m())+enc > 0.) {
490 G4cout <<
" momentum non conservation after (1): " << G4endl
491 <<
" e " << enc <<
" p " << pnc <<
G4endl;
495 G4bool need_hard_tuning =
true;
499 for (
G4int i=0; i < nfrag; i++) {
501 if (eex > 0.0 && eex + encMeV >= 0.0) {
506 G4double newMass = fragMom.m() + encMeV;
507 fragMom.setVectM(fragMom.vect(), newMass);
508 need_hard_tuning =
false;
512 }
else if (nnuc > 0) {
513 for (
G4int i=0; i < nnuc; i++) {
515 if (eex > 0.0 && eex + encMeV >= 0.0) {
517 need_hard_tuning =
false;
521 if (need_hard_tuning && encMeV > 0.) {
523 need_hard_tuning =
false;
527 if (!need_hard_tuning) {
534 G4cout <<
" trying hard (particle-pair) tuning" <<
G4endl;
565 std::pair<G4int, G4int> tune_particles = tune_par.first;
566 G4int mom_ind = tune_par.second;
569 (tune_particles.first >= 0 && tune_particles.second >= 0 &&
570 mom_ind >= G4LorentzVector::X);
572 if (!tuning_possible) {
578 G4cout <<
" p1 " << tune_particles.first <<
" p2 " << tune_particles.second
579 <<
" ind " << mom_ind <<
G4endl;
595 pnc = mom_non_cons.rho();
596 enc = mom_non_cons.e();
598 on_shell = (std::fabs(enc) < accuracy || pnc < accuracy);
601 G4cout <<
" momentum non conservation tuning: " << G4endl
602 <<
" e " << enc <<
" p " << pnc
603 << (
on_shell?
" success":
" FAILURE") << G4endl;
620 std::pair<std::pair<G4int, G4int>,
G4int>
623 G4cout <<
" >>> G4CollisionOutput::selectPairToTune" <<
G4endl;
625 std::pair<G4int, G4int> tup(-1, -1);
627 std::pair<std::pair<G4int, G4int>,
G4int> tuner(tup, i3);
634 G4double pcut = 0.3 * std::sqrt(1.88 * std::fabs(de));
643 for (
G4int l = G4LorentzVector::X; l<=G4LorentzVector::Z; l++) {
644 if (mom1[l]*mom2[l]<0.0) {
645 if (std::fabs(mom1[l])>pcut && std::fabs(mom2[l])>pcut) {
646 G4double psum = std::fabs(mom1[l]) + std::fabs(mom2[l]);
661 if (i3 < 0)
return tuner;
667 tuner.first.first = (p1>0.) ? ibest1 : ibest2;
668 tuner.first.second = (p1>0.) ? ibest2 : ibest1;
670 tuner.first.first = (p1<0.) ? ibest2 : ibest1;
671 tuner.first.second = (p1<0.) ? ibest1 : ibest2;
680 G4int mom_ind)
const {
682 G4cout <<
" >>> G4CollisionOutput::tuneSelectedPair" <<
G4endl;
685 G4double R = 0.5 * (newE12*newE12 + mom2.e()*mom2.e() - mom1.e()*mom1.e()) / newE12;
686 G4double Q = -(mom1[mom_ind] + mom2[mom_ind]) / newE12;
688 G4double W = (R * Q + mom2[mom_ind]) * UDQ;
689 G4double V = (mom2.e()*mom2.e() - R*R) * UDQ;
698 G4double x1 = -(W + std::sqrt(DET));
699 G4double x2 = -(W - std::sqrt(DET));
707 if (R + Q * x1 >= 0.0) {
712 if (!xset && x2 > 0.0) {
713 if (R + Q * x2 >= 0.0) {
720 if (R + Q * x1 >= 0.) {
725 if (!xset && x2 < 0.0) {
726 if (R + Q * x2 >= 0.0) {
void trivialise(G4InuclParticle *bullet, G4InuclParticle *target)
G4CollisionOutput & operator=(const G4CollisionOutput &right)
G4bool tuneSelectedPair(G4LorentzVector &mom1, G4LorentzVector &mom2, G4int mom_index) const
std::vector< G4InuclNuclei >::iterator nucleiIterator
void addOutgoingNuclei(const std::vector< G4InuclNuclei > &nuclea)
std::vector< G4Fragment > recoilFragments
CLHEP::HepLorentzRotation G4LorentzRotation
G4LorentzVector getMomentum() const
void removeRecoilFragment(G4int index=-1)
void printCollisionOutput(std::ostream &os=G4cout) const
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
G4LorentzVector rotate(const G4LorentzVector &mom) const
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
G4LorentzVector getTotalOutputMomentum() const
const G4InuclElementaryParticle & getParticle() const
const G4String & GetParticleName() const
G4int GetAtomicNumber() const
G4LorentzVector mom_non_cons
std::vector< G4ReactionProduct * > G4ReactionProductVector
void add(const G4CollisionOutput &right)
G4GLOB_DLL std::ostream G4cout
G4int getTotalStrangeness() const
G4int numberOfOutgoingParticles() const
void boostToLabFrame(const G4LorentzConvertor &convertor)
void removeOutgoingParticle(G4int index)
std::pair< std::pair< G4int, G4int >, G4int > selectPairToTune(G4double de) const
void rotateEvent(const G4LorentzRotation &rotate)
G4int GetAtomicMass() const
std::vector< G4InuclElementaryParticle > outgoingParticles
G4int numberOfOutgoingNuclei() const
std::vector< G4InuclElementaryParticle >::iterator particleIterator
G4bool reflectionNeeded() const
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
G4int getTotalBaryonNumber() const
const G4double x[NPOINTSGL]
static const G4Fragment emptyFragment
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
G4int getTotalCharge() const
std::vector< G4InuclElementaryParticle >::iterator particleIterator
void setRemainingExitationEnergy()
std::vector< G4Fragment >::iterator fragmentIterator
G4int numberOfFragments() const
void removeOutgoingNucleus(G4int index)
const G4Fragment & getRecoilFragment(G4int index=0) const
static const G4double pos
void setOnShell(G4InuclParticle *bullet, G4InuclParticle *target)
std::vector< G4InuclNuclei > outgoingNuclei
CLHEP::HepLorentzVector G4LorentzVector