Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4INCLXXInterface Class Reference

INCL++ intra-nuclear cascade. More...

#include <G4INCLXXInterface.hh>

Inheritance diagram for G4INCLXXInterface:
Collaboration diagram for G4INCLXXInterface:

Public Member Functions

 G4INCLXXInterface (G4VPreCompoundModel *const aPreCompound=0)
 
 ~G4INCLXXInterface ()
 
G4int operator== (G4INCLXXInterface &right)
 
G4int operator!= (G4INCLXXInterface &right)
 
G4ReactionProductVectorPropagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
 
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
 
void DeleteModel ()
 
virtual void ModelDescription (std::ostream &outFile) const
 
G4String const & GetDeExcitationModelName () const
 
- Public Member Functions inherited from G4VIntraNuclearTransportModel
 G4VIntraNuclearTransportModel (const G4String &mName="CascadeModel", G4VPreCompoundModel *ptr=nullptr)
 
virtual ~G4VIntraNuclearTransportModel ()
 
virtual G4ReactionProductVectorPropagateNuclNucl (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4V3DNucleus *theProjectileNucleus)
 
void SetDeExcitation (G4VPreCompoundModel *ptr)
 
void Set3DNucleus (G4V3DNucleus *const value)
 
void SetPrimaryProjectile (const G4HadProjectile &aPrimary)
 
const G4StringGetModelName () const
 
virtual void PropagateModelDescription (std::ostream &outFile) const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4StringGetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
virtual const std::pair
< G4double, G4double
GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double,
G4double
GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void InitialiseModel ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4VIntraNuclearTransportModel
G4V3DNucleusGet3DNucleus () const
 
G4VPreCompoundModelGetDeExcitation () const
 
const G4HadProjectileGetPrimaryProjectile () const
 
- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4VIntraNuclearTransportModel
G4String theTransportModelName
 
G4V3DNucleusthe3DNucleus
 
G4VPreCompoundModeltheDeExcitation
 
const G4HadProjectilethePrimaryProjectile
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

INCL++ intra-nuclear cascade.

Interface for INCL++. This interface handles basic hadron bullet particles (protons, neutrons, pions), as well as light ions.

Example usage in case of protons:

* inclModel -> SetMinEnergy(0.0 * MeV); // Set the energy limits
* inclModel -> SetMaxEnergy(3.0 * GeV);
*
* G4ProtonInelasticProcess* protonInelasticProcess = new G4ProtonInelasticProcess();
* G4ProtonInelasticCrossSection* protonInelasticCrossSection = new G4ProtonInelasticCrossSection();
*
* protonInelasticProcess -> RegisterMe(inclModel);
* protonInelasticProcess -> AddDataSet(protonInelasticCrossSection);
*
* particle = G4Proton::Proton();
* processManager = particle -> GetProcessManager();
* processManager -> AddDiscreteProcess(protonInelasticProcess);
*

The same setup procedure is needed for neutron, pion and generic-ion inelastic processes as well.

Definition at line 103 of file G4INCLXXInterface.hh.

Constructor & Destructor Documentation

G4INCLXXInterface::G4INCLXXInterface ( G4VPreCompoundModel *const  aPreCompound = 0)

Definition at line 59 of file G4INCLXXInterface.cc.

59  :
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);
75  if(!thePreCompoundModel) { thePreCompoundModel = new G4PreCompoundModel; }
76  }
77 
78  // Use the environment variable G4INCLXX_NO_DE_EXCITATION to disable de-excitation
79  if(getenv("G4INCLXX_NO_DE_EXCITATION")) {
80  G4String message = "de-excitation is completely disabled!";
81  theInterfaceStore->EmitWarning(message);
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) {
94  theINCLXXLevelDensity = new G4FissionLevelDensityParameterINCLXX;
95  theFissionChannelCast->SetLevelDensityParameter(theINCLXXLevelDensity);
96  theINCLXXFissionProbability = new G4FissionProbability;
97  theINCLXXFissionProbability->SetFissionLevelDensityParameter(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(getenv("G4INCLXX_DUMP_REMNANT"))
108  dumpRemnantInfo = true;
109  else
110  dumpRemnantInfo = false;
111 
112  theBackupModel = new G4BinaryLightIonReaction;
113  theBackupModelNucleon = new G4BinaryCascade;
114 }
G4VEvaporationChannel * GetFissionChannel()
const char * p
Definition: xmltok.h:285
G4VIntraNuclearTransportModel(const G4String &mName="CascadeModel", G4VPreCompoundModel *ptr=nullptr)
G4VEvaporation * GetEvaporation()
void EmitBigWarning(const G4String &message) const
Emit a BIG warning to G4cout.
void SetEmissionStrategy(G4VEmissionProbability *aFissionProb)
void SetFissionLevelDensityParameter(G4VLevelDensityParameter *aLevelDensity)
static G4INCLXXInterfaceStore * GetInstance()
Get the singleton instance.
G4ExcitationHandler * GetExcitationHandler() const
Revised level-density parameter for fission after INCL++.
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
G4HadronicInteraction * FindModel(const G4String &name)
static G4HadronicInteractionRegistry * Instance()
void EmitWarning(const G4String &message)
Emit a warning to G4cout.
void SetLevelDensityParameter(G4VLevelDensityParameter *aLevelDensity)

Here is the call graph for this function:

G4INCLXXInterface::~G4INCLXXInterface ( )

Definition at line 116 of file G4INCLXXInterface.cc.

117 {
118  delete theINCLXXLevelDensity;
119  delete theINCLXXFissionProbability;
120 }

Member Function Documentation

G4HadFinalState * G4INCLXXInterface::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus theNucleus 
)
virtual

Main method to apply the INCL physics model.

Parameters
aTrackthe projectile particle
theNucleustarget nucleus
Returns
the output of the INCL physics model

Implements G4HadronicInteraction.

Definition at line 167 of file G4INCLXXInterface.cc.

168 {
169  G4ParticleDefinition const * const trackDefinition = aTrack.GetDefinition();
170  const G4bool isIonTrack = trackDefinition->GetParticleType()==G4GenericIon::GenericIon()->GetParticleType();
171  const G4int trackA = trackDefinition->GetAtomicMass();
172  const G4int trackZ = (G4int) trackDefinition->GetPDGCharge();
173  const G4int nucleusA = theNucleus.GetA_asInt();
174  const G4int nucleusZ = theNucleus.GetZ_asInt();
175 
176  // For reactions induced by weird projectiles (e.g. He2), bail out
177  if((isIonTrack && (trackZ<=0 || trackA<=trackZ)) || (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
178  theResult.Clear();
179  theResult.SetStatusChange(isAlive);
180  theResult.SetEnergyChange(aTrack.GetKineticEnergy());
181  theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
182  return &theResult;
183  }
184 
185  // For reactions on nucleons, use the backup model (without complaining)
186  if(trackA<=1 && nucleusA<=1) {
187  return theBackupModelNucleon->ApplyYourself(aTrack, theNucleus);
188  }
189 
190  // For systems heavier than theMaxProjMassINCL, use another model (typically
191  // BIC)
192  const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
193  if(trackA > theMaxProjMassINCL && nucleusA > theMaxProjMassINCL) {
194  if(!complainedAboutBackupModel) {
195  complainedAboutBackupModel = true;
196  std::stringstream ss;
197  ss << "INCL++ refuses to handle reactions between nuclei with A>"
198  << theMaxProjMassINCL
199  << ". A backup model ("
200  << theBackupModel->GetModelName()
201  << ") will be used instead.";
202  theInterfaceStore->EmitBigWarning(ss.str());
203  }
204  return theBackupModel->ApplyYourself(aTrack, theNucleus);
205  }
206 
207  // For energies lower than cascadeMinEnergyPerNucleon, use PreCompound
208  const G4double cascadeMinEnergyPerNucleon = theInterfaceStore->GetCascadeMinEnergyPerNucleon();
209  const G4double trackKinE = aTrack.GetKineticEnergy();
210  if((trackDefinition==G4Neutron::NeutronDefinition() || trackDefinition==G4Proton::ProtonDefinition())
211  && trackKinE < cascadeMinEnergyPerNucleon) {
212  if(!complainedAboutPreCompound) {
213  complainedAboutPreCompound = true;
214  std::stringstream ss;
215  ss << "INCL++ refuses to handle nucleon-induced reactions below "
216  << cascadeMinEnergyPerNucleon / MeV
217  << " MeV. A PreCoumpound model ("
218  << thePreCompoundModel->GetModelName()
219  << ") will be used instead.";
220  theInterfaceStore->EmitBigWarning(ss.str());
221  }
222  return thePreCompoundModel->ApplyYourself(aTrack, theNucleus);
223  }
224 
225  // Calculate the total four-momentum in the entrance channel
226  const G4double theNucleusMass = theIonTable->GetIonMass(nucleusZ, nucleusA);
227  const G4double theTrackMass = trackDefinition->GetPDGMass();
228  const G4double theTrackEnergy = trackKinE + theTrackMass;
229  const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
230  const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
231  const G4ThreeVector theTrackMomentum = aTrack.Get4Momentum().getV().unit() * theTrackMomentumAbs;
232 
233  G4LorentzVector goodTrack4Momentum(theTrackMomentum, theTrackEnergy);
234  G4LorentzVector fourMomentumIn;
235  fourMomentumIn.setE(theTrackEnergy + theNucleusMass);
236  fourMomentumIn.setVect(theTrackMomentum);
237 
238  // Check if inverse kinematics should be used
239  const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus);
240 
241  // If we are running in inverse kinematics, redefine aTrack and theNucleus
242  G4LorentzRotation *toInverseKinematics = NULL;
243  G4LorentzRotation *toDirectKinematics = NULL;
244  G4HadProjectile const *aProjectileTrack = &aTrack;
245  G4Nucleus *theTargetNucleus = &theNucleus;
246  if(inverseKinematics) {
247  G4ParticleDefinition *oldTargetDef = theIonTable->GetIon(nucleusZ, nucleusA, 0);
248  const G4ParticleDefinition *oldProjectileDef = trackDefinition;
249 
250  if(oldProjectileDef != 0 && oldTargetDef != 0) {
251  const G4int newTargetA = oldProjectileDef->GetAtomicMass();
252  const G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
253 
254  if(newTargetA > 0 && newTargetZ > 0) {
255  // This should give us the same energy per nucleon
256  theTargetNucleus = new G4Nucleus(newTargetA, newTargetZ);
257  toInverseKinematics = new G4LorentzRotation(goodTrack4Momentum.boostVector());
258  G4LorentzVector theProjectile4Momentum(0.0, 0.0, 0.0, theNucleusMass);
259  G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
260  aProjectileTrack = new G4HadProjectile(swappedProjectileParticle);
261  } else {
262  G4String message = "badly defined target after swapping. Falling back to normal (non-swapped) mode.";
263  theInterfaceStore->EmitWarning(message);
264  toInverseKinematics = new G4LorentzRotation;
265  }
266  } else {
267  G4String message = "oldProjectileDef or oldTargetDef was null";
268  theInterfaceStore->EmitWarning(message);
269  toInverseKinematics = new G4LorentzRotation;
270  }
271  }
272 
273  // INCL assumes the projectile particle is going in the direction of
274  // the Z-axis. Here we construct proper rotation to convert the
275  // momentum vectors of the outcoming particles to the original
276  // coordinate system.
277  G4LorentzVector projectileMomentum = aProjectileTrack->Get4Momentum();
278 
279  // INCL++ assumes that the projectile is going in the direction of
280  // the z-axis. In principle, if the coordinate system used by G4
281  // hadronic framework is defined differently we need a rotation to
282  // transform the INCL++ reaction products to the appropriate
283  // frame. Please note that it isn't necessary to apply this
284  // transform to the projectile because when creating the INCL++
285  // projectile particle; \see{toINCLKineticEnergy} needs to use only the
286  // projectile energy (direction is simply assumed to be along z-axis).
287  G4RotationMatrix toZ;
288  toZ.rotateZ(-projectileMomentum.phi());
289  toZ.rotateY(-projectileMomentum.theta());
290  G4RotationMatrix toLabFrame3 = toZ.inverse();
291  G4LorentzRotation toLabFrame(toLabFrame3);
292  // However, it turns out that the projectile given to us by G4
293  // hadronic framework is already going in the direction of the
294  // z-axis so this rotation is actually unnecessary. Both toZ and
295  // toLabFrame turn out to be unit matrices as can be seen by
296  // uncommenting the folowing two lines:
297  // G4cout <<"toZ = " << toZ << G4endl;
298  // G4cout <<"toLabFrame = " << toLabFrame << G4endl;
299 
300  theResult.Clear(); // Make sure the output data structure is clean.
301  theResult.SetStatusChange(stopAndKill);
302 
303  std::list<G4Fragment> remnants;
304 
305  const G4int maxTries = 200;
306  G4int nTries = 0;
307  // INCL can generate transparent events. However, this is meaningful
308  // only in the standalone code. In Geant4 we must "force" INCL to
309  // produce a valid cascade.
310  G4bool eventIsOK = false;
311  do {
312  const G4INCL::ParticleSpecies theSpecies = toINCLParticleSpecies(*aProjectileTrack);
313  const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack);
314 
315  // The INCL model will be created at the first use
317 
318  const G4INCL::EventInfo eventInfo = theINCLModel->processEvent(theSpecies, kineticEnergy, theTargetNucleus->GetA_asInt(), theTargetNucleus->GetZ_asInt());
319  // eventIsOK = !eventInfo.transparent && nTries < maxTries;
320  eventIsOK = !eventInfo.transparent;
321  if(eventIsOK) {
322 
323  // If the collision was run in inverse kinematics, prepare to boost back
324  // all the reaction products
325  if(inverseKinematics) {
326  toDirectKinematics = new G4LorentzRotation(toInverseKinematics->inverse());
327  }
328 
329  G4LorentzVector fourMomentumOut;
330 
331  for(G4int i = 0; i < eventInfo.nParticles; ++i) {
332  G4int A = eventInfo.A[i];
333  G4int Z = eventInfo.Z[i];
334  G4int PDGCode = eventInfo.PDGCode[i];
335  // G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl;
336  G4double kinE = eventInfo.EKin[i];
337  G4double px = eventInfo.px[i];
338  G4double py = eventInfo.py[i];
339  G4double pz = eventInfo.pz[i];
340  G4DynamicParticle *p = toG4Particle(A, Z, PDGCode, kinE, px, py, pz);
341  if(p != 0) {
342  G4LorentzVector momentum = p->Get4Momentum();
343 
344  // Apply the toLabFrame rotation
345  momentum *= toLabFrame;
346  // Apply the inverse-kinematics boost
347  if(inverseKinematics) {
348  momentum *= *toDirectKinematics;
349  momentum.setVect(-momentum.vect());
350  }
351 
352  // Set the four-momentum of the reaction products
353  p->Set4Momentum(momentum);
354  fourMomentumOut += momentum;
355  theResult.AddSecondary(p);
356 
357  } else {
358  G4String message = "the model produced a particle that couldn't be converted to Geant4 particle.";
359  theInterfaceStore->EmitWarning(message);
360  }
361  }
362 
363  for(G4int i = 0; i < eventInfo.nRemnants; ++i) {
364  const G4int A = eventInfo.ARem[i];
365  const G4int Z = eventInfo.ZRem[i];
366  // G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl;
367  const G4double kinE = eventInfo.EKinRem[i];
368  const G4double px = eventInfo.pxRem[i];
369  const G4double py = eventInfo.pyRem[i];
370  const G4double pz = eventInfo.pzRem[i];
371  G4ThreeVector spin(
372  eventInfo.jxRem[i]*hbar_Planck,
373  eventInfo.jyRem[i]*hbar_Planck,
374  eventInfo.jzRem[i]*hbar_Planck
375  );
376  const G4double excitationE = eventInfo.EStarRem[i];
377  const G4double nuclearMass = G4NucleiProperties::GetNuclearMass(A, Z) + excitationE;
378  const G4double scaling = remnant4MomentumScaling(nuclearMass,
379  kinE,
380  px, py, pz);
381  G4LorentzVector fourMomentum(scaling * px, scaling * py, scaling * pz,
382  nuclearMass + kinE);
383  if(std::abs(scaling - 1.0) > 0.01) {
384  std::stringstream ss;
385  ss << "momentum scaling = " << scaling
386  << "\n Lorentz vector = " << fourMomentum
387  << ")\n A = " << A << ", Z = " << Z
388  << "\n E* = " << excitationE << ", nuclearMass = " << nuclearMass
389  << "\n remnant i=" << i << ", nRemnants=" << eventInfo.nRemnants
390  << "\n Reaction was: " << aTrack.GetKineticEnergy()/MeV
391  << "-MeV " << trackDefinition->GetParticleName() << " + "
392  << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
393  << ", in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
394  theInterfaceStore->EmitWarning(ss.str());
395  }
396 
397  // Apply the toLabFrame rotation
398  fourMomentum *= toLabFrame;
399  spin *= toLabFrame3;
400  // Apply the inverse-kinematics boost
401  if(inverseKinematics) {
402  fourMomentum *= *toDirectKinematics;
403  fourMomentum.setVect(-fourMomentum.vect());
404  }
405 
406  fourMomentumOut += fourMomentum;
407  G4Fragment remnant(A, Z, fourMomentum);
408  remnant.SetAngularMomentum(spin);
409  if(dumpRemnantInfo) {
410  G4cerr << "G4INCLXX_DUMP_REMNANT: " << remnant << " spin: " << spin << G4endl;
411  }
412  remnants.push_back(remnant);
413  }
414 
415  // Check four-momentum conservation
416  const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
417  const G4double energyViolation = std::abs(violation4Momentum.e());
418  const G4double momentumViolation = violation4Momentum.rho();
419  if(energyViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
420  std::stringstream ss;
421  ss << "energy conservation violated by " << energyViolation/MeV << " MeV in "
422  << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
423  << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
424  << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
425  theInterfaceStore->EmitWarning(ss.str());
426  eventIsOK = false;
427  const G4int nSecondaries = theResult.GetNumberOfSecondaries();
428  for(G4int j=0; j<nSecondaries; ++j)
429  delete theResult.GetSecondary(j)->GetParticle();
430  theResult.Clear();
431  theResult.SetStatusChange(stopAndKill);
432  remnants.clear();
433  } else if(momentumViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
434  std::stringstream ss;
435  ss << "momentum conservation violated by " << momentumViolation/MeV << " MeV in "
436  << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
437  << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
438  << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
439  theInterfaceStore->EmitWarning(ss.str());
440  eventIsOK = false;
441  const G4int nSecondaries = theResult.GetNumberOfSecondaries();
442  for(G4int j=0; j<nSecondaries; ++j)
443  delete theResult.GetSecondary(j)->GetParticle();
444  theResult.Clear();
445  theResult.SetStatusChange(stopAndKill);
446  remnants.clear();
447  }
448  }
449  nTries++;
450  } while(!eventIsOK && nTries < maxTries); /* Loop checking, 10.07.2015, D.Mancusi */
451 
452  // Clean up the objects that we created for the inverse kinematics
453  if(inverseKinematics) {
454  delete aProjectileTrack;
455  delete theTargetNucleus;
456  delete toInverseKinematics;
457  delete toDirectKinematics;
458  }
459 
460  if(!eventIsOK) {
461  std::stringstream ss;
462  ss << "maximum number of tries exceeded for the proposed "
463  << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
464  << " + " << theIonTable->GetIonName(nucleusZ, nucleusA, 0)
465  << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
466  theInterfaceStore->EmitWarning(ss.str());
467  theResult.SetStatusChange(isAlive);
468  theResult.SetEnergyChange(aTrack.GetKineticEnergy());
469  theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
470  return &theResult;
471  }
472 
473  // De-excitation:
474 
475  if(theDeExcitation != 0) {
476  for(std::list<G4Fragment>::iterator i = remnants.begin();
477  i != remnants.end(); ++i) {
478  G4ReactionProductVector *deExcitationResult = theDeExcitation->DeExcite((*i));
479 
480  for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
481  fragment != deExcitationResult->end(); ++fragment) {
482  const G4ParticleDefinition *def = (*fragment)->GetDefinition();
483  if(def != 0) {
484  G4DynamicParticle *theFragment = new G4DynamicParticle(def, (*fragment)->GetMomentum());
485  theResult.AddSecondary(theFragment);
486  }
487  }
488 
489  for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
490  fragment != deExcitationResult->end(); ++fragment) {
491  delete (*fragment);
492  }
493  deExcitationResult->clear();
494  delete deExcitationResult;
495  }
496  }
497 
498  remnants.clear();
499 
500  if((theTally = theInterfaceStore->GetTally()))
501  theTally->Tally(aTrack, theNucleus, theResult);
502 
503  return &theResult;
504 }
Short_t nRemnants
Number of remnants.
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetMaxProjMassINCL() const
Getter for theMaxProjMassINCL.
G4double GetCascadeMinEnergyPerNucleon() const
Getter for cascadeMinEnergyPerNucleon.
Int_t PDGCode[maxSizeParticles]
PDG numbering of the particles.
G4HadSecondary * GetSecondary(size_t i)
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
Float_t py[maxSizeParticles]
Particle momentum, y component [MeV/c].
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
HepRotation & rotateY(double delta)
Definition: Rotation.cc:79
const G4String & GetModelName() const
void EmitBigWarning(const G4String &message) const
Emit a BIG warning to G4cout.
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
int G4int
Definition: G4Types.hh:78
Float_t pyRem[maxSizeRemnants]
Remnant momentum, y component [MeV/c].
const G4String & GetParticleName() const
Float_t pz[maxSizeParticles]
Particle momentum, z component [MeV/c].
const G4String & GetIonName(G4int Z, G4int A, G4int lvl=0) const
Definition: G4IonTable.cc:1094
HepRotation inverse() const
G4int GetAtomicNumber() const
double phi() const
Short_t ZRem[maxSizeRemnants]
Remnant charge number.
void SetStatusChange(G4HadFinalStateStatus aS)
Short_t ARem[maxSizeRemnants]
Remnant mass number.
std::vector< G4ReactionProduct * > G4ReactionProductVector
static G4INCLXXInterfaceStore * GetInstance()
Get the singleton instance.
Hep3Vector vect() const
double A(double temperature)
Short_t Z[maxSizeParticles]
Particle charge number.
const G4ParticleDefinition * GetDefinition() const
double theta() const
Short_t nParticles
Number of particles in the final state.
G4INCL::INCL * GetINCLModel()
Get the cached INCL model engine.
Float_t jxRem[maxSizeRemnants]
Remnant angular momentum, x component [ ].
bool G4bool
Definition: G4Types.hh:79
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1335
const EventInfo & processEvent()
G4double GetKineticEnergy() const
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
double rho() const
Float_t EStarRem[maxSizeRemnants]
Remnant excitation energy [MeV].
const G4String & GetParticleType() const
const G4LorentzVector & Get4Momentum() const
G4int GetAtomicMass() const
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
G4LorentzVector Get4Momentum() const
Float_t pxRem[maxSizeRemnants]
Remnant momentum, x component [MeV/c].
Short_t A[maxSizeParticles]
Particle mass number.
void Set4Momentum(const G4LorentzVector &momentum)
void SetEnergyChange(G4double anEnergy)
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:93
G4double GetPDGMass() const
Hep3Vector unit() const
G4DynamicParticle * GetParticle()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:92
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
HepLorentzRotation inverse() const
void setVect(const Hep3Vector &)
Bool_t transparent
True if the event is transparent.
void EmitWarning(const G4String &message)
Emit a warning to G4cout.
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
static constexpr double hbar_Planck
virtual void Tally(G4HadProjectile const &aTrack, G4Nucleus const &theNucleus, G4HadFinalState const &result)=0
G4INCLXXVInterfaceTally * GetTally() const
Getter for the interface tally.
Float_t px[maxSizeParticles]
Particle momentum, x component [MeV/c].
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
void SetMomentumChange(const G4ThreeVector &aV)
G4int GetNumberOfSecondaries() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4GLOB_DLL std::ostream G4cerr
Hep3Vector getV() const
void G4INCLXXInterface::DeleteModel ( )
inline

Definition at line 128 of file G4INCLXXInterface.hh.

128  {
129  delete theINCLModel;
130  theINCLModel = NULL;
131  }
G4String const & G4INCLXXInterface::GetDeExcitationModelName ( ) const

Definition at line 606 of file G4INCLXXInterface.cc.

606  {
607  return theDeExcitation->GetModelName();
608 }
const G4String & GetModelName() const

Here is the call graph for this function:

void G4INCLXXInterface::ModelDescription ( std::ostream &  outFile) const
virtual

Reimplemented from G4VIntraNuclearTransportModel.

Definition at line 591 of file G4INCLXXInterface.cc.

591  {
592  outFile
593  << "The Liège Intranuclear Cascade (INCL++) is a model for reactions induced\n"
594  << "by nucleons, pions and light ion on any nucleus. The reaction is\n"
595  << "described as an avalanche of binary nucleon-nucleon collisions, which can\n"
596  << "lead to the emission of energetic particles and to the formation of an\n"
597  << "excited thermalised nucleus (remnant). The de-excitation of the remnant is\n"
598  << "outside the scope of INCL++ and is typically described by another model.\n\n"
599  << "INCL++ has been reasonably well tested for nucleon (~50 MeV to ~15 GeV),\n"
600  << "pion (idem) and light-ion projectiles (up to A=18, ~10A MeV to 1A GeV).\n"
601  << "Most tests involved target nuclei close to the stability valley, with\n"
602  << "numbers between 4 and 250.\n\n"
603  << "Reference: D. Mancusi et al., Phys. Rev. C90 (2014) 054602\n\n";
604 }
G4int G4INCLXXInterface::operator!= ( G4INCLXXInterface right)
inline

Definition at line 112 of file G4INCLXXInterface.hh.

112  {
113  return (this != &right);
114  }
G4int G4INCLXXInterface::operator== ( G4INCLXXInterface right)
inline

Definition at line 108 of file G4INCLXXInterface.hh.

108  {
109  return (this == &right);
110  }
G4ReactionProductVector * G4INCLXXInterface::Propagate ( G4KineticTrackVector theSecondaries,
G4V3DNucleus theNucleus 
)
virtual

Implements G4VIntraNuclearTransportModel.

Definition at line 506 of file G4INCLXXInterface.cc.

506  {
507  return 0;
508 }

The documentation for this class was generated from the following files: