Main method to apply the INCL physics model.
177 if((isIonTrack && (trackZ<=0 || trackA<=trackZ)) || (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
186 if(trackA<=1 && nucleusA<=1) {
187 return theBackupModelNucleon->
ApplyYourself(aTrack, theNucleus);
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 ("
201 <<
") will be used instead.";
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 ("
219 <<
") will be used instead.";
222 return thePreCompoundModel->
ApplyYourself(aTrack, theNucleus);
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);
235 fourMomentumIn.
setE(theTrackEnergy + theNucleusMass);
236 fourMomentumIn.
setVect(theTrackMomentum);
239 const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus);
245 G4Nucleus *theTargetNucleus = &theNucleus;
246 if(inverseKinematics) {
250 if(oldProjectileDef != 0 && oldTargetDef != 0) {
254 if(newTargetA > 0 && newTargetZ > 0) {
256 theTargetNucleus =
new G4Nucleus(newTargetA, newTargetZ);
259 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
262 G4String message =
"badly defined target after swapping. Falling back to normal (non-swapped) mode.";
267 G4String message =
"oldProjectileDef or oldTargetDef was null";
303 std::list<G4Fragment> remnants;
305 const G4int maxTries = 200;
313 const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack);
325 if(inverseKinematics) {
345 momentum *= toLabFrame;
347 if(inverseKinematics) {
348 momentum *= *toDirectKinematics;
354 fourMomentumOut += momentum;
358 G4String message =
"the model produced a particle that couldn't be converted to Geant4 particle.";
373 eventInfo.
jyRem[i]*hbar_Planck,
374 eventInfo.
jzRem[i]*hbar_Planck
378 const G4double scaling = remnant4MomentumScaling(nuclearMass,
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
393 <<
", in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
398 fourMomentum *= toLabFrame;
401 if(inverseKinematics) {
402 fourMomentum *= *toDirectKinematics;
403 fourMomentum.setVect(-fourMomentum.vect());
406 fourMomentumOut += fourMomentum;
408 remnant.SetAngularMomentum(spin);
409 if(dumpRemnantInfo) {
410 G4cerr <<
"G4INCLXX_DUMP_REMNANT: " << remnant <<
" spin: " << spin <<
G4endl;
412 remnants.push_back(remnant);
416 const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
417 const G4double energyViolation = std::abs(violation4Momentum.
e());
418 const G4double momentumViolation = violation4Momentum.
rho();
420 std::stringstream ss;
421 ss <<
"energy conservation violated by " << energyViolation/
MeV <<
" MeV in "
424 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
428 for(
G4int j=0; j<nSecondaries; ++j)
434 std::stringstream ss;
435 ss <<
"momentum conservation violated by " << momentumViolation/
MeV <<
" MeV in "
438 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
442 for(
G4int j=0; j<nSecondaries; ++j)
450 }
while(!eventIsOK && nTries < maxTries);
453 if(inverseKinematics) {
454 delete aProjectileTrack;
455 delete theTargetNucleus;
456 delete toInverseKinematics;
457 delete toDirectKinematics;
461 std::stringstream ss;
462 ss <<
"maximum number of tries exceeded for the proposed "
464 <<
" + " << theIonTable->
GetIonName(nucleusZ, nucleusA, 0)
465 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
476 for(std::list<G4Fragment>::iterator i = remnants.begin();
477 i != remnants.end(); ++i) {
480 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
481 fragment != deExcitationResult->end(); ++fragment) {
489 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
490 fragment != deExcitationResult->end(); ++fragment) {
493 deExcitationResult->clear();
494 delete deExcitationResult;
500 if((theTally = theInterfaceStore->
GetTally()))
501 theTally->
Tally(aTrack, theNucleus, theResult);
Short_t nRemnants
Number of remnants.
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].
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4Proton * ProtonDefinition()
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
HepRotation & rotateY(double delta)
const G4String & GetModelName() const
void EmitBigWarning(const G4String &message) const
Emit a BIG warning to G4cout.
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
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
HepRotation inverse() const
G4int GetAtomicNumber() 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.
double A(double temperature)
Short_t Z[maxSizeParticles]
Particle charge number.
const G4ParticleDefinition * GetDefinition() 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 [ ].
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
const EventInfo & processEvent()
G4double GetKineticEnergy() const
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
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()
G4double GetPDGMass() const
G4DynamicParticle * GetParticle()
G4VPreCompoundModel * theDeExcitation
HepRotation & rotateZ(double delta)
static constexpr double MeV
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.
G4double GetPDGCharge() const
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()
void SetMomentumChange(const G4ThreeVector &aV)
G4int GetNumberOfSecondaries() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4GLOB_DLL std::ostream G4cerr