34 #define INCLXX_IN_GEANT4_MODE 1
80 :propagationModel(0), theA(208), theZ(82), theS(0),
81 targetInitSuccess(
false),
83 maxUniverseRadius(0.),
84 maxInteractionDistance(0.),
85 fixedImpactParameter(0.),
88 forceTransparent(
false),
92 #ifdef INCLXX_IN_GEANT4_MODE
94 #else // INCLXX_IN_GEANT4_MODE
96 #endif // INCLXX_IN_GEANT4_MODE
147 #ifndef INCLXX_IN_GEANT4_MODE
162 #ifndef INCLXX_IN_GEANT4_MODE
163 NuclearMassTable::deleteTable();
171 #ifndef INCLXX_IN_GEANT4_MODE
172 Logger::deleteLoggerSlave();
183 if(A < 0 || A > 300 || Z < 1 || Z > 200) {
184 INCL_ERROR(
"Unsupported target: A = " << A <<
" Z = " << Z <<
" S = " << S <<
'\n'
185 <<
"Target configuration rejected." <<
'\n');
189 (projectileSpecies.
theZ==projectileSpecies.
theA || projectileSpecies.
theZ==0)) {
190 INCL_ERROR(
"Unsupported projectile: A = " << projectileSpecies.
theA <<
" Z = " << projectileSpecies.
theZ <<
" S = " << projectileSpecies.
theS <<
'\n'
191 <<
"Projectile configuration rejected." <<
'\n');
222 if(projectileSpecies.
theA > 0)
257 INCL_WARN(
"Target initialisation failed for A=" << targetA <<
", Z=" << targetZ <<
", S=" << targetS <<
'\n');
308 INCL_DEBUG(
"Selected impact parameter: " << impactParameter <<
'\n');
314 if(effectiveImpactParameter < 0.) {
331 unsigned long loopCounter = 0;
332 const unsigned long maxLoopCounter = 10000000;
346 if(avatar == 0)
break;
383 INCL_DEBUG(
"Trying compound nucleus" <<
'\n');
387 #ifndef INCLXX_IN_GEANT4_MODE
397 if(projectileRemnant) {
419 #ifdef INCLXX_IN_GEANT4_MODE
421 #endif // INCLXX_IN_GEANT4_MODE
449 INCL_DEBUG(
"Cascade resulted in complete fusion, using realistic fusion kinematics" <<
'\n');
455 INCL_WARN(
"Complete-fusion kinematics yields negative excitation energy, returning a transparent!" <<
'\n');
471 INCL_ERROR(
"Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." <<
'\n');
475 #ifndef INCLXX_IN_GEANT4_MODE
477 globalConservationChecks(
false);
489 #ifndef INCLXX_IN_GEANT4_MODE
491 globalConservationChecks(
true);
530 std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
532 std::shuffle(shuffledComponents.begin(), shuffledComponents.end(),
Random::getAdapter());
535 G4bool atLeastOneNucleonEntering =
false;
536 for(std::vector<Particle*>::const_iterator
p=shuffledComponents.begin(),
e=shuffledComponents.end();
p!=
e; ++
p) {
540 (*p)->getPropagationVelocity(),
542 if(!intersectionInteractionDistance.
exists)
546 atLeastOneNucleonEntering =
true;
559 theCNZ += (*p)->getZ();
560 theCNS += (*p)->getS();
570 if(!success || !atLeastOneNucleonEntering) {
571 INCL_DEBUG(
"No nucleon entering in forced CN, forcing a transparent" <<
'\n');
582 theCNEnergy -= theProjectileRemnant->
getEnergy();
583 theCNMomentum -= theProjectileRemnant->
getMomentum();
594 const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.
mag2();
595 if(theCNInvariantMassSquared<0.) {
600 const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
601 if(theCNExcitationEnergy<0.) {
603 INCL_DEBUG(
"CN excitation energy is negative, forcing a transparent" <<
'\n'
604 <<
" theCNA = " << theCNA <<
'\n'
605 <<
" theCNZ = " << theCNZ <<
'\n'
606 <<
" theCNS = " << theCNS <<
'\n'
607 <<
" theCNEnergy = " << theCNEnergy <<
'\n'
608 <<
" theCNMomentum = (" << theCNMomentum.
getX() <<
", "<< theCNMomentum.
getY() <<
", " << theCNMomentum.
getZ() <<
")" <<
'\n'
609 <<
" theCNExcitationEnergy = " << theCNExcitationEnergy <<
'\n'
610 <<
" theCNSpin = (" << theCNSpin.
getX() <<
", "<< theCNSpin.
getY() <<
", " << theCNSpin.
getZ() <<
")" <<
'\n'
616 INCL_DEBUG(
"CN excitation energy is positive, forcing a CN" <<
'\n'
617 <<
" theCNA = " << theCNA <<
'\n'
618 <<
" theCNZ = " << theCNZ <<
'\n'
619 <<
" theCNS = " << theCNS <<
'\n'
620 <<
" theCNEnergy = " << theCNEnergy <<
'\n'
621 <<
" theCNMomentum = (" << theCNMomentum.
getX() <<
", "<< theCNMomentum.
getY() <<
", " << theCNMomentum.
getZ() <<
")" <<
'\n'
622 <<
" theCNExcitationEnergy = " << theCNExcitationEnergy <<
'\n'
623 <<
" theCNSpin = (" << theCNSpin.
getX() <<
", "<< theCNSpin.
getY() <<
", " << theCNSpin.
getZ() <<
")" <<
'\n'
658 theRecoilFunctor(theSolution.
x);
660 INCL_WARN(
"Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." <<
'\n');
665 #ifndef INCLXX_IN_GEANT4_MODE
666 void INCL::globalConservationChecks(
G4bool afterRecoil) {
672 if(theBalance.
Z != 0) {
675 if(theBalance.
A != 0) {
678 if(theBalance.
S != 0) {
681 G4double EThreshold, pLongThreshold, pTransThreshold;
686 pTransThreshold = 1.;
690 pLongThreshold = 0.1;
691 pTransThreshold = 0.1;
696 if(
std::abs(pLongBalance)>pLongThreshold) {
697 INCL_WARN(
"Violation of longitudinal momentum conservation > " << pLongThreshold <<
" MeV/c. pLongBalance = " << pLongBalance <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber <<
'\n');
699 if(
std::abs(pTransBalance)>pTransThreshold) {
700 INCL_WARN(
"Violation of transverse momentum conservation > " << pTransThreshold <<
" MeV/c. pTransBalance = " << pTransBalance <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber <<
'\n');
715 <<
"), stopping cascade" <<
'\n');
721 INCL_DEBUG(
"No participants in the nucleus and no incoming particles left, stopping cascade" <<
'\n');
728 <<
"), stopping cascade" <<
'\n');
734 INCL_DEBUG(
"Trying to make a compound nucleus, stopping cascade" <<
'\n');
778 G4int nUnmergedSpectators = 0;
781 if(dynSpectators.empty() && geomSpectators.empty()) {
783 }
else if(dynSpectators.size()==1 && geomSpectators.empty()) {
794 nUnmergedSpectators = rejected.size();
802 return nUnmergedSpectators;
816 INCL_DEBUG(
"Initialised interaction distance: r0 = " << r0 <<
'\n'
817 <<
" theNNDistance = " << theNNDistance <<
'\n'
827 for(
IsotopeIter i=theIsotopes.begin(),
e=theIsotopes.end(); i!=
e; ++i) {
831 rMax =
std::max(maximumRadius, rMax);
837 rMax =
std::max(maximumRadius, rMax);