34 #define INCLXX_IN_GEANT4_MODE 1
49 #include "G4String.hh"
62 thePreCompoundModel(aPreCompound),
65 complainedAboutBackupModel(
false),
66 complainedAboutPreCompound(
false),
68 theINCLXXLevelDensity(NULL),
69 theINCLXXFissionProbability(NULL)
79 if(std::getenv(
"G4INCLXX_NO_DE_EXCITATION")) {
93 if(theFissionChannelCast) {
107 if(std::getenv(
"G4INCLXX_DUMP_REMNANT"))
131 std::stringstream ss;
132 ss <<
"the model does not know how to handle a collision between a "
154 if(pA > theMaxProjMassINCL)
156 else if(tA > theMaxProjMassINCL)
173 if((isIonTrack && (trackZ<=0 || trackA<=trackZ)) || (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
182 if(trackA<=1 && nucleusA<=1) {
189 if(trackA > theMaxProjMassINCL && nucleusA > theMaxProjMassINCL) {
192 std::stringstream ss;
193 ss <<
"INCL++ refuses to handle reactions between nuclei with A>"
194 << theMaxProjMassINCL
195 <<
". A backup model ("
197 <<
") will be used instead.";
207 && trackKinE < cascadeMinEnergyPerNucleon) {
210 std::stringstream ss;
211 ss <<
"INCL++ refuses to handle nucleon-induced reactions below "
212 << cascadeMinEnergyPerNucleon /
MeV
213 <<
" MeV. A PreCoumpound model ("
215 <<
") will be used instead.";
224 const G4double theTrackEnergy = trackKinE + theTrackMass;
225 const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
226 const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
231 fourMomentumIn.
setE(theTrackEnergy + theNucleusMass);
232 fourMomentumIn.
setVect(theTrackMomentum);
241 G4Nucleus *theTargetNucleus = &theNucleus;
242 if(inverseKinematics) {
246 if(oldProjectileDef != 0 && oldTargetDef != 0) {
250 if(newTargetA > 0 && newTargetZ > 0) {
252 theTargetNucleus =
new G4Nucleus(newTargetA, newTargetZ);
255 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
258 G4String message =
"badly defined target after swapping. Falling back to normal (non-swapped) mode.";
299 std::list<G4Fragment> remnants;
301 const G4int maxTries = 200;
321 if(inverseKinematics) {
341 momentum *= toLabFrame;
343 if(inverseKinematics) {
344 momentum *= *toDirectKinematics;
354 G4String message =
"the model produced a particle that couldn't be converted to Geant4 particle.";
369 eventInfo.
jyRem[i]*hbar_Planck,
370 eventInfo.
jzRem[i]*hbar_Planck
379 if(
std::abs(scaling - 1.0) > 0.01) {
380 std::stringstream ss;
381 ss <<
"momentum scaling = " << scaling
382 <<
"\n Lorentz vector = " << fourMomentum
383 <<
")\n A = " << A <<
", Z = " << Z
384 <<
"\n E* = " << excitationE <<
", nuclearMass = " << nuclearMass
385 <<
"\n remnant i=" << i <<
", nRemnants=" << eventInfo.
nRemnants
389 <<
", in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
394 fourMomentum *= toLabFrame;
397 if(inverseKinematics) {
398 fourMomentum *= *toDirectKinematics;
402 fourMomentumOut += fourMomentum;
406 G4cerr <<
"G4INCLXX_DUMP_REMNANT: " << remnant <<
" spin: " << spin <<
G4endl;
408 remnants.push_back(remnant);
412 const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
414 const G4double momentumViolation = violation4Momentum.
rho();
416 std::stringstream ss;
417 ss <<
"energy conservation violated by " << energyViolation/
MeV <<
" MeV in "
420 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
424 for(
G4int j=0; j<nSecondaries; ++j)
430 std::stringstream ss;
431 ss <<
"momentum conservation violated by " << momentumViolation/
MeV <<
" MeV in "
434 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
438 for(
G4int j=0; j<nSecondaries; ++j)
446 }
while(!eventIsOK && nTries < maxTries);
449 if(inverseKinematics) {
450 delete aProjectileTrack;
451 delete theTargetNucleus;
452 delete toInverseKinematics;
453 delete toDirectKinematics;
457 std::stringstream ss;
458 ss <<
"maximum number of tries exceeded for the proposed "
461 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
472 for(std::list<G4Fragment>::iterator i = remnants.begin();
473 i != remnants.end(); ++i) {
476 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
477 fragment != deExcitationResult->end(); ++fragment) {
485 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
486 fragment != deExcitationResult->end(); ++fragment) {
489 deExcitationResult->clear();
490 delete deExcitationResult;
563 else if(A > 0 && Z > 0 && A > Z) {
589 const G4double p2 = px*px + py*py + pz*pz;
591 const G4double pnew2 = kineticE*kineticE + 2.0*kineticE*
mass;
592 return std::sqrt(pnew2)/std::sqrt(p2);
600 <<
"The Li�ge Intranuclear Cascade (INCL++) is a model for reactions induced\n"
601 <<
"by nucleons, pions and light ion on any nucleus. The reaction is\n"
602 <<
"described as an avalanche of binary nucleon-nucleon collisions, which can\n"
603 <<
"lead to the emission of energetic particles and to the formation of an\n"
604 <<
"excited thermalised nucleus (remnant). The de-excitation of the remnant is\n"
605 <<
"outside the scope of INCL++ and is typically described by another model.\n\n"
606 <<
"INCL++ has been reasonably well tested for nucleon (~50 MeV to ~15 GeV),\n"
607 <<
"pion (idem) and light-ion projectiles (up to A=18, ~10A MeV to 1A GeV).\n"
608 <<
"Most tests involved target nuclei close to the stability valley, with\n"
609 <<
"numbers between 4 and 250.\n\n"
610 <<
"Reference: D. Mancusi et al., Phys. Rev. C90 (2014) 054602\n\n";