34 #define INCLXX_IN_GEANT4_MODE 1
45 #ifndef G4INCLNucleus_hh
46 #define G4INCLNucleus_hh 1
70 theInitialZ(charge), theInitialA(mass), theInitialS(strangess),
71 theNpInitial(0), theNnInitial(0),
72 theNpionplusInitial(0), theNpionminusInitial(0),
73 theNkaonplusInitial(0), theNkaonminusInitial(0),
74 initialInternalEnergy(0.),
75 incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
76 initialCenterOfMass(0.,0.,0.),
80 theUniverseRadius(universeRadius),
81 isNucleusNucleus(
false),
82 theProjectileRemnant(NULL),
146 for(
ParticleIter iter=created.begin(),
e=created.end(); iter!=
e; ++iter) {
148 if(!(*iter)->isOutOfWell()) {
149 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
154 for(
ParticleIter iter=deleted.begin(),
e=deleted.end(); iter!=
e; ++iter) {
159 for(
ParticleIter iter=modified.begin(),
e=modified.end(); iter!=
e; ++iter) {
161 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
165 for(
ParticleIter iter=out.begin(),
e=out.end(); iter!=
e; ++iter) {
166 if((*iter)->isCluster()) {
169 #ifdef INCLXX_IN_GEANT4_MODE
179 totalEnergy += (*iter)->getEnergy();
180 theA -= (*iter)->getA();
181 theZ -= (*iter)->getZ();
182 theS -= (*iter)->getS();
188 for(
ParticleIter iter=entering.begin(),
e=entering.end(); iter!=
e; ++iter) {
190 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
196 INCL_DEBUG(
"A Particle is entering below the Fermi sea:" <<
'\n' << finalstate->
print() <<
'\n');
199 for(
ParticleIter iter=entering.begin(),
e=entering.end(); iter!=
e; ++iter) {
206 INCL_ERROR(
"Energy nonconservation! Energy at the beginning of the event = "
208 <<
" and after interaction = "
209 << totalEnergy <<
'\n'
210 << finalstate->
print());
215 INCL_WARN(
"Useless Nucleus::propagateParticles -method called." <<
'\n');
222 if((*p)->isNucleon())
223 totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
224 else if((*p)->isResonance())
226 else if((*p)->isHyperon())
229 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
252 theSpin -= (*p)->getAngularMomentum();
274 cm += (*p)->getPosition() *
mass;
290 std::stringstream ss;
291 ss <<
"Particles in the nucleus:" <<
'\n'
292 <<
"Inside:" <<
'\n';
296 ss <<
"index = " << counter <<
'\n'
300 ss <<
"Outgoing:" <<
'\n';
312 if((*i)->isDelta()) deltas.push_back((*i));
314 if(deltas.empty())
return false;
317 INCL_DEBUG(
"Decay outgoing delta particle:" <<
'\n'
318 << (*i)->print() <<
'\n');
320 const G4double deltaMass = (*i)->getMass();
325 (*i)->setEnergy((*i)->getMass());
338 newMomentum *= decayMomentum / newMomentum.
mag();
350 nucleon->
boost(beta);
369 const G4bool unphysicalRemnant = (theZ<0 || theZ>
theA);
377 if((*i)->isDelta()) deltas.push_back((*i));
381 INCL_DEBUG(
"Decay inside delta particle:" <<
'\n'
382 << (*i)->print() <<
'\n');
388 if(unphysicalRemnant) {
389 INCL_WARN(
"Forcing delta decay inside an unphysical remnant (A=" <<
theA
390 <<
", Z=" <<
theZ <<
"). Might lead to energy-violation warnings."
409 if(unphysicalRemnant) {
410 INCL_DEBUG(
"Remnant is unphysical: Z=" <<
theZ <<
", A=" <<
theA <<
", emitting all the pions" <<
'\n');
422 const G4bool unphysicalRemnant = (theZ<0 || theZ>
theA);
423 if(unphysicalRemnant){
425 INCL_WARN(
"Remnant is unphysical: Z=" <<
theZ <<
", A=" <<
theA <<
", too much strange particles? -> all emit" <<
'\n');
437 if((*i)->isSigma() || (*i)->isAntiKaon()) stranges.push_back((*i));
438 else if((*i)->isNucleon() && (*i)->getZ() == 1) protons.push_back((*i));
439 else if((*i)->isNucleon() && (*i)->getZ() == 0) neutrons.push_back((*i));
442 if((stranges.size() > protons.size()) || (stranges.size() > neutrons.size())){
443 INCL_WARN(
"Remnant is unphysical: Nproton=" << protons.size() <<
", Nneutron=" << neutrons.size() <<
", Strange particles : " << stranges.size() <<
'\n');
451 for(
ParticleIter i=stranges.begin(),
e=stranges.end(); i!=
e; ++i) {
452 INCL_DEBUG(
"Absorbe inside strange particles:" <<
'\n'
453 << (*i)->print() <<
'\n');
456 decay =
new DecayAvatar((*protonIter), (*i), 0.0,
this,
true);
460 decay =
new DecayAvatar((*neutronIter), (*i), 0.0,
this,
true);
463 else if(
Random::shoot()*(protons.size() + neutrons.size()) < protons.size()){
464 decay =
new DecayAvatar((*protonIter), (*i), 0.0,
this,
true);
468 decay =
new DecayAvatar((*neutronIter), (*i), 0.0,
this,
true);
487 if(pionResonances.empty())
return false;
489 for(
ParticleIter i=pionResonances.begin(),
e=pionResonances.end(); i!=
e; ++i) {
490 INCL_DEBUG(
"Decay outgoing pionResonances particle:" <<
'\n'
491 << (*i)->print() <<
'\n');
493 const G4double pionResonanceMass = (*i)->getMass();
498 (*i)->setEnergy((*i)->getMass());
506 Particle *
const theCreatedParticle1 = created.front();
508 if (created.size() == 1) {
513 newMomentum *= decayMomentum / newMomentum.
mag();
519 theCreatedParticle1->
boost(beta);
525 theModifiedParticle->
boost(beta);
529 else if (created.size() == 2) {
530 Particle *
const theCreatedParticle2 = created.back();
532 theCreatedParticle1->
boost(beta);
534 theCreatedParticle2->
boost(beta);
536 theModifiedParticle->
boost(beta);
542 INCL_ERROR(
"Wrong number (< 2) of created particles during the decay of a pion resonance");
557 if(neutralsigma.empty())
return false;
559 for(
ParticleIter i=neutralsigma.begin(),
e=neutralsigma.end(); i!=
e; ++i) {
560 INCL_DEBUG(
"Decay outgoing neutral sigma:" <<
'\n'
561 << (*i)->print() <<
'\n');
563 const G4double neutralsigmaMass = (*i)->getMass();
568 (*i)->setEnergy((*i)->getMass());
576 Particle *
const theCreatedParticle = created.front();
578 if (created.size() == 1) {
583 newMomentum *= decayMomentum / newMomentum.
mag();
588 theCreatedParticle->
boost(beta);
594 theModifiedParticle->
boost(beta);
599 INCL_ERROR(
"Wrong number (!= 1) of created particles during the decay of a sigma zero");
612 if((*i)->getType() ==
KZero || (*i)->getType() ==
KZeroBar) neutralkaon.push_back((*i));
614 if(neutralkaon.empty())
return false;
616 for(
ParticleIter i=neutralkaon.begin(),
e=neutralkaon.end(); i!=
e; ++i) {
617 INCL_DEBUG(
"Transform outgoing neutral kaon:" <<
'\n'
618 << (*i)->print() <<
'\n');
635 if((*i)->isCluster()) clusters.push_back((*i));
637 if(clusters.empty())
return false;
639 for(
ParticleIter i=clusters.begin(),
e=clusters.end(); i!=
e; ++i) {
642 #ifdef INCLXX_IN_GEANT4_MODE
648 for(
ParticleIter j=decayProducts.begin(), end=decayProducts.end(); j!=end; ++j){
662 for(
ParticleIter j=decayProducts.begin(),
e=decayProducts.end(); j!=
e; ++j){
675 INCL_WARN(
"Forcing emissions of all pions in the nucleus." <<
'\n');
678 const G4double tinyPionEnergy = 0.1;
686 INCL_DEBUG(
"Forcing emission of the following particle: "
687 << thePion->
print() <<
'\n');
693 if(kineticEnergyOutside > 0.0)
700 toEject.push_back(thePion);
716 INCL_DEBUG(
"Forcing emissions of all strange particles in the nucleus." <<
'\n');
725 if((*i)->isSigma() || (*i)->isAntiKaon()) {
727 INCL_DEBUG(
"Forcing emission of the following particle: "
728 << theParticle->
print() <<
'\n');
734 if(kineticEnergyOutside > 0.0)
743 toEject.push_back(theParticle);
759 INCL_DEBUG(
"Forcing emissions of all Lambda in the nucleus." <<
'\n');
768 if((*i)->isLambda()) {
770 INCL_DEBUG(
"Forcing emission of the following particle: "
771 << theLambda->
print() <<
'\n');
777 if(kineticEnergyOutside > 0.0)
785 toEject.push_back(theLambda);
793 return toEject.size();
802 INCL_DEBUG(
"Forcing emissions of all Kaon in the nucleus." <<
'\n');
813 INCL_DEBUG(
"Forcing emission of the following particle: "
814 << theKaon->
print() <<
'\n');
820 if(kineticEnergyOutside > 0.0)
828 toEject.push_back(theKaon);
837 return toEject.size() != 0;
846 if(nEventCollisions==0 && nEventDecays==0 && nEventClusters==0)
883 if(outgoing.size() == 2) {
885 INCL_DEBUG(
"Two particles in the outgoing channel, applying exact two-body kinematics" <<
'\n');
889 Particle *p1 = outgoing.front(), *p2 = outgoing.back();
892 p1->
boost(aBoostVector);
900 p2->adjustEnergyFromMomentum();
902 p1->
boost(-aBoostVector);
903 p2->boost(-aBoostVector);
907 INCL_DEBUG(
"Trying to adjust final-state momenta to achieve energy and momentum conservation" <<
'\n');
909 const G4int maxIterations=8;
911 G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
913 std::vector<ThreeVector> minMomenta;
916 minMomenta.reserve(outgoing.size());
919 totalMomentum.
setX(0.0);
920 totalMomentum.
setY(0.0);
921 totalMomentum.
setZ(0.0);
923 totalMomentum += (*i)->getMomentum();
928 totalEnergy += (*i)->getEnergy();
931 for(
G4int iterations=0; iterations < maxIterations; ++iterations) {
934 oldOldOldVal = oldOldVal;
938 if(iterations%2 == 0) {
944 pOldTot += (*i)->getMomentum().
mag();
945 for(
ParticleIter i=outgoing.begin(),
e=outgoing.end(); i!=
e; ++i) {
947 (*i)->setMomentum(mom + deltaP*mom.
mag()/pOldTot);
948 (*i)->adjustEnergyFromMomentum();
954 for(
ParticleIter i=outgoing.begin(),
e=outgoing.end(); i!=
e; ++i) {
956 G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.
mag2();
958 (*i)->setEnergy((*i)->getEnergy()*energyScale);
959 (*i)->adjustMomentumFromEnergy();
965 totalMomentum.
setX(0.0);
966 totalMomentum.
setY(0.0);
967 totalMomentum.
setZ(0.0);
969 for(
ParticleIter i=outgoing.begin(),
e=outgoing.end(); i!=
e; ++i) {
970 totalMomentum += (*i)->getMomentum();
971 totalEnergy += (*i)->getEnergy();
977 INCL_DEBUG(
"Merit function: val=" << val <<
", oldVal=" << oldVal <<
", oldOldVal=" << oldOldVal <<
", oldOldOldVal=" << oldOldOldVal <<
'\n');
981 INCL_DEBUG(
"New minimum found, storing the particle momenta" <<
'\n');
984 minMomenta.push_back((*i)->getMomentum());
988 if(val > oldOldVal && oldVal > oldOldOldVal) {
989 INCL_DEBUG(
"Search is diverging, breaking out of the iteration loop: val=" << val <<
", oldVal=" << oldVal <<
", oldOldVal=" << oldOldVal <<
", oldOldOldVal=" << oldOldOldVal <<
'\n');
999 std::vector<ThreeVector>::const_iterator
v = minMomenta.begin();
1000 for(
ParticleIter i=outgoing.begin(),
e=outgoing.end(); i!=
e; ++i, ++
v) {
1001 (*i)->setMomentum(*v);
1002 (*i)->adjustEnergyFromMomentum();
1012 G4bool isNucleonAbsorption =
false;
1014 G4bool isPionAbsorption =
false;
1020 isPionAbsorption =
true;
1031 if(outgoingParticles.size() == 0 &&
1034 isNucleonAbsorption =
true;
1041 for(
ParticleIter i=outgoingParticles.begin(),
e=outgoingParticles.end(); i!=
e; ++i ) {
1044 if(isPionAbsorption) {
1045 if((*i)->isPion()) {
1046 isPionAbsorption =
false;
1052 eventInfo->
A[eventInfo->
nParticles] = (*i)->getA();
1053 eventInfo->
Z[eventInfo->
nParticles] = (*i)->getZ();
1054 eventInfo->
S[eventInfo->
nParticles] = (*i)->getS();
1056 eventInfo->
EKin[eventInfo->
nParticles] = (*i)->getKineticEnergy();
1064 eventInfo->
history.push_back(
"");
1089 INCL_WARN(
"Negative excitation energy in projectile-like remnant! EStarRem = " << eventInfo->
EStarRem[eventInfo->
nRemnants] <<
'\n');
1158 theBalance.
Z = theEventInfo.
Zp + theEventInfo.
Zt;
1159 theBalance.
A = theEventInfo.
Ap + theEventInfo.
At;
1160 theBalance.
S = theEventInfo.
Sp + theEventInfo.
St;
1167 for(
ParticleIter i=outgoingParticles.begin(),
e=outgoingParticles.end(); i!=
e; ++i ) {
1168 theBalance.
Z -= (*i)->getZ();
1169 theBalance.
A -= (*i)->getA();
1170 theBalance.
S -= (*i)->getS();
1173 theBalance.
energy -= (*i)->getEnergy();
1174 theBalance.
momentum -= (*i)->getMomentum();
1190 theBalance.
Z -=
getZ();
1191 theBalance.
A -=
getA();
1192 theBalance.
S -=
getS();
1220 const G4double anExcitationEnergy = aMass