ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLXXInterface.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLXXInterface.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 #include <cmath>
39 
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4INCLXXInterface.hh"
43 #include "G4GenericIon.hh"
44 #include "G4INCLCascade.hh"
46 #include "G4ReactionProduct.hh"
49 #include "G4String.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
53 #include "G4INCLVersion.hh"
54 #include "G4VEvaporation.hh"
55 #include "G4VEvaporationChannel.hh"
56 #include "G4CompetitiveFission.hh"
58 
61  theINCLModel(NULL),
62  thePreCompoundModel(aPreCompound),
63  theInterfaceStore(G4INCLXXInterfaceStore::GetInstance()),
64  theTally(NULL),
65  complainedAboutBackupModel(false),
66  complainedAboutPreCompound(false),
67  theIonTable(G4IonTable::GetIonTable()),
68  theINCLXXLevelDensity(NULL),
69  theINCLXXFissionProbability(NULL)
70 {
71  if(!thePreCompoundModel) {
74  thePreCompoundModel = static_cast<G4VPreCompoundModel*>(p);
76  }
77 
78  // Use the environment variable G4INCLXX_NO_DE_EXCITATION to disable de-excitation
79  if(std::getenv("G4INCLXX_NO_DE_EXCITATION")) {
80  G4String message = "de-excitation is completely disabled!";
82  theDeExcitation = 0;
83  } else {
86  theDeExcitation = static_cast<G4VPreCompoundModel*>(p);
88 
89  // set the fission parameters for G4ExcitationHandler
90  G4VEvaporationChannel * const theFissionChannel =
92  G4CompetitiveFission * const theFissionChannelCast = dynamic_cast<G4CompetitiveFission *>(theFissionChannel);
93  if(theFissionChannelCast) {
95  theFissionChannelCast->SetLevelDensityParameter(theINCLXXLevelDensity);
98  theFissionChannelCast->SetEmissionStrategy(theINCLXXFissionProbability);
99  theInterfaceStore->EmitBigWarning("INCL++/G4ExcitationHandler uses its own level-density parameter for fission");
100  } else {
101  theInterfaceStore->EmitBigWarning("INCL++/G4ExcitationHandler could not use its own level-density parameter for fission");
102  }
103  }
104 
105  // use the envvar G4INCLXX_DUMP_REMNANT to dump information about the
106  // remnants on stdout
107  if(std::getenv("G4INCLXX_DUMP_REMNANT"))
108  dumpRemnantInfo = true;
109  else
110  dumpRemnantInfo = false;
111 
114 }
115 
117 {
118  delete theINCLXXLevelDensity;
120 }
121 
123  // Use direct kinematics if the projectile is a nucleon or a pion
124  const G4ParticleDefinition *projectileDef = aTrack.GetDefinition();
125  if(std::abs(projectileDef->GetBaryonNumber()) < 2) // Every non-composite particle. abs()-> in case of anti-nucleus projectile
126  return false;
127 
128  // Here all projectiles should be nuclei
129  const G4int pA = projectileDef->GetAtomicMass();
130  if(pA<=0) {
131  std::stringstream ss;
132  ss << "the model does not know how to handle a collision between a "
133  << projectileDef->GetParticleName() << " projectile and a Z="
134  << theNucleus.GetZ_asInt() << ", A=" << theNucleus.GetA_asInt();
136  return true;
137  }
138 
139  // If either nucleus is a LCP (A<=4), run the collision as light on heavy
140  const G4int tA = theNucleus.GetA_asInt();
141  if(tA<=4 || pA<=4) {
142  if(pA<tA)
143  return false;
144  else
145  return true;
146  }
147 
148  // If one of the nuclei is heavier than theMaxProjMassINCL, run the collision
149  // as light on heavy.
150  // Note that here we are sure that either the projectile or the target is
151  // smaller than theMaxProjMass; otherwise theBackupModel would have been
152  // called.
153  const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
154  if(pA > theMaxProjMassINCL)
155  return true;
156  else if(tA > theMaxProjMassINCL)
157  return false;
158  else
159  // In all other cases, use the global setting
161 }
162 
164 {
165  G4ParticleDefinition const * const trackDefinition = aTrack.GetDefinition();
166  const G4bool isIonTrack = trackDefinition->GetParticleType()==G4GenericIon::GenericIon()->GetParticleType();
167  const G4int trackA = trackDefinition->GetAtomicMass();
168  const G4int trackZ = (G4int) trackDefinition->GetPDGCharge();
169  const G4int nucleusA = theNucleus.GetA_asInt();
170  const G4int nucleusZ = theNucleus.GetZ_asInt();
171 
172  // For reactions induced by weird projectiles (e.g. He2), bail out
173  if((isIonTrack && (trackZ<=0 || trackA<=trackZ)) || (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
174  theResult.Clear();
178  return &theResult;
179  }
180 
181  // For reactions on nucleons, use the backup model (without complaining)
182  if(trackA<=1 && nucleusA<=1) {
183  return theBackupModelNucleon->ApplyYourself(aTrack, theNucleus);
184  }
185 
186  // For systems heavier than theMaxProjMassINCL, use another model (typically
187  // BIC)
188  const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
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.";
199  }
200  return theBackupModel->ApplyYourself(aTrack, theNucleus);
201  }
202 
203  // For energies lower than cascadeMinEnergyPerNucleon, use PreCompound
204  const G4double cascadeMinEnergyPerNucleon = theInterfaceStore->GetCascadeMinEnergyPerNucleon();
205  const G4double trackKinE = aTrack.GetKineticEnergy();
206  if((trackDefinition==G4Neutron::NeutronDefinition() || trackDefinition==G4Proton::ProtonDefinition())
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.";
217  }
218  return thePreCompoundModel->ApplyYourself(aTrack, theNucleus);
219  }
220 
221  // Calculate the total four-momentum in the entrance channel
222  const G4double theNucleusMass = theIonTable->GetIonMass(nucleusZ, nucleusA);
223  const G4double theTrackMass = trackDefinition->GetPDGMass();
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);
227  const G4ThreeVector theTrackMomentum = aTrack.Get4Momentum().getV().unit() * theTrackMomentumAbs;
228 
229  G4LorentzVector goodTrack4Momentum(theTrackMomentum, theTrackEnergy);
230  G4LorentzVector fourMomentumIn;
231  fourMomentumIn.setE(theTrackEnergy + theNucleusMass);
232  fourMomentumIn.setVect(theTrackMomentum);
233 
234  // Check if inverse kinematics should be used
235  const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus);
236 
237  // If we are running in inverse kinematics, redefine aTrack and theNucleus
238  G4LorentzRotation *toInverseKinematics = NULL;
239  G4LorentzRotation *toDirectKinematics = NULL;
240  G4HadProjectile const *aProjectileTrack = &aTrack;
241  G4Nucleus *theTargetNucleus = &theNucleus;
242  if(inverseKinematics) {
243  G4ParticleDefinition *oldTargetDef = theIonTable->GetIon(nucleusZ, nucleusA, 0);
244  const G4ParticleDefinition *oldProjectileDef = trackDefinition;
245 
246  if(oldProjectileDef != 0 && oldTargetDef != 0) {
247  const G4int newTargetA = oldProjectileDef->GetAtomicMass();
248  const G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
249 
250  if(newTargetA > 0 && newTargetZ > 0) {
251  // This should give us the same energy per nucleon
252  theTargetNucleus = new G4Nucleus(newTargetA, newTargetZ);
253  toInverseKinematics = new G4LorentzRotation(goodTrack4Momentum.boostVector());
254  G4LorentzVector theProjectile4Momentum(0.0, 0.0, 0.0, theNucleusMass);
255  G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
256  aProjectileTrack = new G4HadProjectile(swappedProjectileParticle);
257  } else {
258  G4String message = "badly defined target after swapping. Falling back to normal (non-swapped) mode.";
259  theInterfaceStore->EmitWarning(message);
260  toInverseKinematics = new G4LorentzRotation;
261  }
262  } else {
263  G4String message = "oldProjectileDef or oldTargetDef was null";
264  theInterfaceStore->EmitWarning(message);
265  toInverseKinematics = new G4LorentzRotation;
266  }
267  }
268 
269  // INCL assumes the projectile particle is going in the direction of
270  // the Z-axis. Here we construct proper rotation to convert the
271  // momentum vectors of the outcoming particles to the original
272  // coordinate system.
273  G4LorentzVector projectileMomentum = aProjectileTrack->Get4Momentum();
274 
275  // INCL++ assumes that the projectile is going in the direction of
276  // the z-axis. In principle, if the coordinate system used by G4
277  // hadronic framework is defined differently we need a rotation to
278  // transform the INCL++ reaction products to the appropriate
279  // frame. Please note that it isn't necessary to apply this
280  // transform to the projectile because when creating the INCL++
281  // projectile particle; \see{toINCLKineticEnergy} needs to use only the
282  // projectile energy (direction is simply assumed to be along z-axis).
283  G4RotationMatrix toZ;
284  toZ.rotateZ(-projectileMomentum.phi());
285  toZ.rotateY(-projectileMomentum.theta());
286  G4RotationMatrix toLabFrame3 = toZ.inverse();
287  G4LorentzRotation toLabFrame(toLabFrame3);
288  // However, it turns out that the projectile given to us by G4
289  // hadronic framework is already going in the direction of the
290  // z-axis so this rotation is actually unnecessary. Both toZ and
291  // toLabFrame turn out to be unit matrices as can be seen by
292  // uncommenting the folowing two lines:
293  // G4cout <<"toZ = " << toZ << G4endl;
294  // G4cout <<"toLabFrame = " << toLabFrame << G4endl;
295 
296  theResult.Clear(); // Make sure the output data structure is clean.
298 
299  std::list<G4Fragment> remnants;
300 
301  const G4int maxTries = 200;
302  G4int nTries = 0;
303  // INCL can generate transparent events. However, this is meaningful
304  // only in the standalone code. In Geant4 we must "force" INCL to
305  // produce a valid cascade.
306  G4bool eventIsOK = false;
307  do {
308  const G4INCL::ParticleSpecies theSpecies = toINCLParticleSpecies(*aProjectileTrack);
309  const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack);
310 
311  // The INCL model will be created at the first use
313 
314  const G4INCL::EventInfo eventInfo = theINCLModel->processEvent(theSpecies, kineticEnergy, theTargetNucleus->GetA_asInt(), theTargetNucleus->GetZ_asInt(),0);
315  // eventIsOK = !eventInfo.transparent && nTries < maxTries;
316  eventIsOK = !eventInfo.transparent;
317  if(eventIsOK) {
318 
319  // If the collision was run in inverse kinematics, prepare to boost back
320  // all the reaction products
321  if(inverseKinematics) {
322  toDirectKinematics = new G4LorentzRotation(toInverseKinematics->inverse());
323  }
324 
325  G4LorentzVector fourMomentumOut;
326 
327  for(G4int i = 0; i < eventInfo.nParticles; ++i) {
328  G4int A = eventInfo.A[i];
329  G4int Z = eventInfo.Z[i];
330  G4int PDGCode = eventInfo.PDGCode[i];
331  // G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl;
332  G4double kinE = eventInfo.EKin[i];
333  G4double px = eventInfo.px[i];
334  G4double py = eventInfo.py[i];
335  G4double pz = eventInfo.pz[i];
336  G4DynamicParticle *p = toG4Particle(A, Z, PDGCode, kinE, px, py, pz);
337  if(p != 0) {
339 
340  // Apply the toLabFrame rotation
341  momentum *= toLabFrame;
342  // Apply the inverse-kinematics boost
343  if(inverseKinematics) {
344  momentum *= *toDirectKinematics;
345  momentum.setVect(-momentum.vect());
346  }
347 
348  // Set the four-momentum of the reaction products
349  p->Set4Momentum(momentum);
350  fourMomentumOut += momentum;
352 
353  } else {
354  G4String message = "the model produced a particle that couldn't be converted to Geant4 particle.";
355  theInterfaceStore->EmitWarning(message);
356  }
357  }
358 
359  for(G4int i = 0; i < eventInfo.nRemnants; ++i) {
360  const G4int A = eventInfo.ARem[i];
361  const G4int Z = eventInfo.ZRem[i];
362  // G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl;
363  const G4double kinE = eventInfo.EKinRem[i];
364  const G4double px = eventInfo.pxRem[i];
365  const G4double py = eventInfo.pyRem[i];
366  const G4double pz = eventInfo.pzRem[i];
367  G4ThreeVector spin(
368  eventInfo.jxRem[i]*hbar_Planck,
369  eventInfo.jyRem[i]*hbar_Planck,
370  eventInfo.jzRem[i]*hbar_Planck
371  );
372  const G4double excitationE = eventInfo.EStarRem[i];
373  const G4double nuclearMass = G4NucleiProperties::GetNuclearMass(A, Z) + excitationE;
374  const G4double scaling = remnant4MomentumScaling(nuclearMass,
375  kinE,
376  px, py, pz);
377  G4LorentzVector fourMomentum(scaling * px, scaling * py, scaling * pz,
378  nuclearMass + kinE);
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
386  << "\n Reaction was: " << aTrack.GetKineticEnergy()/MeV
387  << "-MeV " << trackDefinition->GetParticleName() << " + "
388  << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
389  << ", in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
390  theInterfaceStore->EmitWarning(ss.str());
391  }
392 
393  // Apply the toLabFrame rotation
394  fourMomentum *= toLabFrame;
395  spin *= toLabFrame3;
396  // Apply the inverse-kinematics boost
397  if(inverseKinematics) {
398  fourMomentum *= *toDirectKinematics;
399  fourMomentum.setVect(-fourMomentum.vect());
400  }
401 
402  fourMomentumOut += fourMomentum;
403  G4Fragment remnant(A, Z, fourMomentum);
404  remnant.SetAngularMomentum(spin);
405  if(dumpRemnantInfo) {
406  G4cerr << "G4INCLXX_DUMP_REMNANT: " << remnant << " spin: " << spin << G4endl;
407  }
408  remnants.push_back(remnant);
409  }
410 
411  // Check four-momentum conservation
412  const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
413  const G4double energyViolation = std::abs(violation4Momentum.e());
414  const G4double momentumViolation = violation4Momentum.rho();
415  if(energyViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
416  std::stringstream ss;
417  ss << "energy conservation violated by " << energyViolation/MeV << " MeV in "
418  << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
419  << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
420  << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
421  theInterfaceStore->EmitWarning(ss.str());
422  eventIsOK = false;
423  const G4int nSecondaries = theResult.GetNumberOfSecondaries();
424  for(G4int j=0; j<nSecondaries; ++j)
425  delete theResult.GetSecondary(j)->GetParticle();
426  theResult.Clear();
428  remnants.clear();
429  } else if(momentumViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
430  std::stringstream ss;
431  ss << "momentum conservation violated by " << momentumViolation/MeV << " MeV in "
432  << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
433  << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
434  << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
435  theInterfaceStore->EmitWarning(ss.str());
436  eventIsOK = false;
437  const G4int nSecondaries = theResult.GetNumberOfSecondaries();
438  for(G4int j=0; j<nSecondaries; ++j)
439  delete theResult.GetSecondary(j)->GetParticle();
440  theResult.Clear();
442  remnants.clear();
443  }
444  }
445  nTries++;
446  } while(!eventIsOK && nTries < maxTries); /* Loop checking, 10.07.2015, D.Mancusi */
447 
448  // Clean up the objects that we created for the inverse kinematics
449  if(inverseKinematics) {
450  delete aProjectileTrack;
451  delete theTargetNucleus;
452  delete toInverseKinematics;
453  delete toDirectKinematics;
454  }
455 
456  if(!eventIsOK) {
457  std::stringstream ss;
458  ss << "maximum number of tries exceeded for the proposed "
459  << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
460  << " + " << theIonTable->GetIonName(nucleusZ, nucleusA, 0)
461  << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
462  theInterfaceStore->EmitWarning(ss.str());
466  return &theResult;
467  }
468 
469  // De-excitation:
470 
471  if(theDeExcitation != 0) {
472  for(std::list<G4Fragment>::iterator i = remnants.begin();
473  i != remnants.end(); ++i) {
474  G4ReactionProductVector *deExcitationResult = theDeExcitation->DeExcite((*i));
475 
476  for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
477  fragment != deExcitationResult->end(); ++fragment) {
478  const G4ParticleDefinition *def = (*fragment)->GetDefinition();
479  if(def != 0) {
480  G4DynamicParticle *theFragment = new G4DynamicParticle(def, (*fragment)->GetMomentum());
481  theResult.AddSecondary(theFragment);
482  }
483  }
484 
485  for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
486  fragment != deExcitationResult->end(); ++fragment) {
487  delete (*fragment);
488  }
489  deExcitationResult->clear();
490  delete deExcitationResult;
491  }
492  }
493 
494  remnants.clear();
495 
497  theTally->Tally(aTrack, theNucleus, theResult);
498 
499  return &theResult;
500 }
501 
503  return 0;
504 }
505 
507  if( pdef == G4Proton::Proton()) return G4INCL::Proton;
508  else if(pdef == G4Neutron::Neutron()) return G4INCL::Neutron;
509  else if(pdef == G4PionPlus::PionPlus()) return G4INCL::PiPlus;
510  else if(pdef == G4PionMinus::PionMinus()) return G4INCL::PiMinus;
511  else if(pdef == G4PionZero::PionZero()) return G4INCL::PiZero;
512  else if(pdef == G4KaonPlus::KaonPlus()) return G4INCL::KPlus;
513  else if(pdef == G4KaonMinus::KaonMinus()) return G4INCL::KMinus;
514  else if(pdef == G4Deuteron::Deuteron()) return G4INCL::Composite;
515  else if(pdef == G4Triton::Triton()) return G4INCL::Composite;
516  else if(pdef == G4He3::He3()) return G4INCL::Composite;
517  else if(pdef == G4Alpha::Alpha()) return G4INCL::Composite;
519  else return G4INCL::UnknownParticle;
520 }
521 
523  const G4ParticleDefinition *pdef = aTrack.GetDefinition();
524  const G4INCL::ParticleType theType = toINCLParticleType(pdef);
525  if(theType!=G4INCL::Composite)
526  return G4INCL::ParticleSpecies(theType);
527  else {
528  G4INCL::ParticleSpecies theSpecies;
529  theSpecies.theType=theType;
530  theSpecies.theA=pdef->GetAtomicMass();
531  theSpecies.theZ=pdef->GetAtomicNumber();
532  return theSpecies;
533  }
534 }
535 
537  return aTrack.GetKineticEnergy();
538 }
539 
541  if(PDGCode == 2212) return G4Proton::Proton();
542  else if(PDGCode == 2112) return G4Neutron::Neutron();
543  else if(PDGCode == 211) return G4PionPlus::PionPlus();
544  else if(PDGCode == 111) return G4PionZero::PionZero();
545  else if(PDGCode == -211) return G4PionMinus::PionMinus();
546 
547  else if(PDGCode == 221) return G4Eta::Eta();
548  else if(PDGCode == 22) return G4Gamma::Gamma();
549 
550  else if(PDGCode == 3122) return G4Lambda::Lambda();
551  else if(PDGCode == 3222) return G4SigmaPlus::SigmaPlus();
552  else if(PDGCode == 3212) return G4SigmaZero::SigmaZero();
553  else if(PDGCode == 3112) return G4SigmaMinus::SigmaMinus();
554  else if(PDGCode == 321) return G4KaonPlus::KaonPlus();
555  else if(PDGCode == -321) return G4KaonMinus::KaonMinus();
556  else if(PDGCode == 130) return G4KaonZeroLong::KaonZeroLong();
557  else if(PDGCode == 310) return G4KaonZeroShort::KaonZeroShort();
558 
559  else if(PDGCode == 1002) return G4Deuteron::Deuteron();
560  else if(PDGCode == 1003) return G4Triton::Triton();
561  else if(PDGCode == 2003) return G4He3::He3();
562  else if(PDGCode == 2004) return G4Alpha::Alpha();
563  else if(A > 0 && Z > 0 && A > Z) { // Returns ground state ion definition. No hyper-nucleus allows in Geant4
564  return theIonTable->GetIon(Z, A, 0);
565  } else { // Error, unrecognized particle
566  return 0;
567  }
568 }
569 
571  G4double kinE,
572  G4double px,
573  G4double py, G4double pz) const {
574  const G4ParticleDefinition *def = toG4ParticleDefinition(A, Z, PDGCode);
575  if(def == 0) { // Check if we have a valid particle definition
576  return 0;
577  }
578  const G4double energy = kinE * MeV;
579  const G4ThreeVector momentum(px, py, pz);
580  const G4ThreeVector momentumDirection = momentum.unit();
581  G4DynamicParticle *p = new G4DynamicParticle(def, momentumDirection, energy);
582  return p;
583 }
584 
586  G4double kineticE,
587  G4double px, G4double py,
588  G4double pz) const {
589  const G4double p2 = px*px + py*py + pz*pz;
590  if(p2 > 0.0) {
591  const G4double pnew2 = kineticE*kineticE + 2.0*kineticE*mass;
592  return std::sqrt(pnew2)/std::sqrt(p2);
593  } else {
594  return 1.0;
595  }
596 }
597 
598 void G4INCLXXInterface::ModelDescription(std::ostream& outFile) const {
599  outFile
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";
611 }
612 
614  return theDeExcitation->GetModelName();
615 }