33 #define INCLXX_IN_GEANT4_MODE 1
47 #include "G4String.hh"
56 thePreCompoundModel(aPreCompound),
58 complainedAboutBackupModel(false),
59 complainedAboutPreCompound(false),
70 if(getenv(
"G4INCLXX_NO_DE_EXCITATION")) {
71 G4String message =
"de-excitation is completely disabled!";
83 if(getenv(
"G4INCLXX_DUMP_REMNANT"))
109 std::stringstream ss;
110 ss <<
"the model does not know how to handle a collision between a "
132 if(pA > theMaxProjMassINCL)
134 else if(tA > theMaxProjMassINCL)
152 && theNucleus.
GetA_asInt() > theMaxProjMassINCL) {
155 std::stringstream ss;
156 ss <<
"INCL++ refuses to handle reactions between nuclei with A>"
157 << theMaxProjMassINCL
158 <<
". A backup model ("
160 <<
") will be used instead.";
171 && trackKinE < cascadeMinEnergyPerNucleon) {
174 std::stringstream ss;
175 ss <<
"INCL++ refuses to handle nucleon-induced reactions below "
176 << cascadeMinEnergyPerNucleon /
MeV
177 <<
" MeV. A PreCoumpound model ("
179 <<
") will be used instead.";
190 const G4double theTrackEnergy = trackKinE + theTrackMass;
191 const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
192 const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
197 fourMomentumIn.setE(theTrackEnergy + theNucleusMass);
198 fourMomentumIn.setVect(theTrackMomentum);
207 G4Nucleus *theTargetNucleus = &theNucleus;
208 if(inverseKinematics) {
212 if(oldProjectileDef != 0 && oldTargetDef != 0) {
216 if(newTargetA > 0 && newTargetZ > 0) {
218 theTargetNucleus =
new G4Nucleus(newTargetA, newTargetZ);
221 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
224 G4String message =
"badly defined target after swapping. Falling back to normal (non-swapped) mode.";
229 G4String message =
"oldProjectileDef or oldTargetDef was null";
250 toZ.rotateZ(-projectileMomentum.phi());
251 toZ.rotateY(-projectileMomentum.theta());
265 std::list<G4Fragment> remnants;
267 const G4int maxTries = 200;
287 if(inverseKinematics) {
306 momentum *= toLabFrame;
308 if(inverseKinematics) {
309 momentum *= *toDirectKinematics;
310 momentum.setVect(-momentum.vect());
315 fourMomentumOut += momentum;
319 G4String message =
"the model produced a particle that couldn't be converted to Geant4 particle.";
333 eventInfo.
jxRem[i]*hbar_Planck,
334 eventInfo.
jyRem[i]*hbar_Planck,
335 eventInfo.
jzRem[i]*hbar_Planck
344 if(std::abs(scaling - 1.0) > 0.01) {
345 std::stringstream ss;
346 ss <<
"momentum scaling = " << scaling
347 <<
"\n Lorentz vector = " << fourMomentum
348 <<
")\n A = " << A <<
", Z = " << Z
349 <<
"\n E* = " << excitationE <<
", nuclearMass = " << nuclearMass
350 <<
"\n remnant i=" << i <<
", nRemnants=" << eventInfo.
nRemnants
354 <<
", in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
359 fourMomentum *= toLabFrame;
362 if(inverseKinematics) {
363 fourMomentum *= *toDirectKinematics;
364 fourMomentum.setVect(-fourMomentum.vect());
367 fourMomentumOut += fourMomentum;
371 G4cerr <<
"G4INCLXX_DUMP_REMNANT: " << remnant <<
" spin: " << spin <<
G4endl;
373 remnants.push_back(remnant);
377 const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
378 const G4double energyViolation = std::abs(violation4Momentum.e());
379 const G4double momentumViolation = violation4Momentum.rho();
381 std::stringstream ss;
382 ss <<
"energy conservation violated by " << energyViolation/
MeV <<
" MeV in "
385 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
389 for(
G4int j=0; j<nSecondaries; ++j)
395 std::stringstream ss;
396 ss <<
"momentum conservation violated by " << momentumViolation/
MeV <<
" MeV in "
399 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
403 for(
G4int j=0; j<nSecondaries; ++j)
411 }
while(!eventIsOK && nTries < maxTries);
414 if(inverseKinematics) {
415 delete aProjectileTrack;
416 delete theTargetNucleus;
417 delete toInverseKinematics;
418 delete toDirectKinematics;
422 std::stringstream ss;
423 ss <<
"maximum number of tries exceeded for the proposed "
426 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
437 for(std::list<G4Fragment>::iterator i = remnants.begin();
438 i != remnants.end(); i++) {
441 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
442 fragment != deExcitationResult->end(); ++fragment) {
450 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
451 fragment != deExcitationResult->end(); ++fragment) {
454 deExcitationResult->clear();
455 delete deExcitationResult;
508 else if(A == 3 && Z == 2)
return G4He3::He3();
510 else if(A > 0 && Z > 0 && A > Z) {
536 const G4double p2 = px*px + py*py + pz*pz;
538 const G4double pnew2 = kineticE*kineticE + 2.0*kineticE*mass;
539 return std::sqrt(pnew2)/std::sqrt(p2);
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetMaxProjMassINCL() const
Getter for theMaxProjMassINCL.
G4double GetCascadeMinEnergyPerNucleon() const
Getter for cascadeMinEnergyPerNucleon.
G4HadronicInteraction * theBackupModel
G4HadSecondary * GetSecondary(size_t i)
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
Float_t py[maxSizeParticles]
Particle momentum, y component [MeV/c].
void SetAngularMomentum(const G4ThreeVector &value)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4Proton * ProtonDefinition()
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
G4DynamicParticle * toG4Particle(G4int A, G4int Z, G4double kinE, G4double px, G4double py, G4double pz) const
Convert an INCL particle to a G4DynamicParticle.
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
G4INCLXXInterface(G4VPreCompoundModel *const aPreCompound=0)
const G4String & GetModelName() const
void EmitBigWarning(const G4String &message) const
Emit a BIG warning to G4cout.
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
G4INCL::ParticleType toINCLParticleType(G4ParticleDefinition const *const) const
Convert G4ParticleDefinition to corresponding INCL particle type.
G4bool complainedAboutPreCompound
Float_t pyRem[maxSizeRemnants]
Remnant momentum, y component [MeV/c].
Header file for the G4INCLXXInterfaceStore class.
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
G4HadronicInteraction * theBackupModelNucleon
G4int GetAtomicNumber() const
Short_t ZRem[maxSizeRemnants]
Remnant charge number.
G4VPreCompoundModel * thePreCompoundModel
G4bool complainedAboutBackupModel
void SetStatusChange(G4HadFinalStateStatus aS)
Short_t ARem[maxSizeRemnants]
Remnant mass number.
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4double remnant4MomentumScaling(G4double mass, G4double kineticE, G4double px, G4double py, G4double pz) const
Rescale remnant momentum if necessary.
static G4INCLXXInterfaceStore * GetInstance()
Get the singleton instance.
Singleton class for configuring the INCL++ Geant4 interface.
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
Main method to apply the INCL physics model.
Short_t Z[maxSizeParticles]
Particle charge number.
const G4ParticleDefinition * GetDefinition() const
Short_t nParticles
Number of particles in the final state.
G4HadFinalState theResult
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].
static G4Triton * Triton()
static G4Proton * Proton()
static G4PionPlus * PionPlus()
const G4String & GetParticleType() const
G4bool AccurateProjectile(const G4HadProjectile &aTrack, const G4Nucleus &theTargetNucleus) const
static G4Neutron * Neutron()
static const G4double A[nN]
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
static G4Deuteron * Deuteron()
G4int GetAtomicMass() const
G4ParticleDefinition * toG4ParticleDefinition(G4int A, G4int Z) const
Convert A and Z to a G4ParticleDefinition.
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
G4LorentzVector Get4Momentum() const
G4INCL::ParticleSpecies toINCLParticleSpecies(G4HadProjectile const &) const
Convert G4HadProjectile to corresponding INCL particle species.
Int_t nRemnants
Number of remnants.
Float_t pxRem[maxSizeRemnants]
Remnant momentum, x component [MeV/c].
G4HadronicInteraction * FindModel(const G4String &name)
Short_t A[maxSizeParticles]
Particle mass number.
void Set4Momentum(const G4LorentzVector &momentum)
void SetEnergyChange(G4double anEnergy)
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)=0
static G4GenericIon * GenericIon()
G4double GetPDGMass() const
G4double toINCLKineticEnergy(G4HadProjectile const &) const
Convert G4HadProjectile to corresponding INCL particle kinetic energy.
G4double energy(const ThreeVector &p, const G4double m)
G4INCL::INCL * theINCLModel
static G4HadronicInteractionRegistry * Instance()
G4DynamicParticle * GetParticle()
static G4PionMinus * PionMinus()
G4bool GetAccurateProjectile() const
Getter for accurateProjectile.
G4VPreCompoundModel * theDeExcitation
std::map< G4String, G4AttDef > * GetInstance(const G4String &storeKey, G4bool &isNew)
G4IonTable *const theIonTable
Bool_t transparent
True if the event is transparent.
void EmitWarning(const G4String &message)
Emit a warning to G4cout.
G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
G4double GetPDGCharge() const
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
void AddSecondary(G4DynamicParticle *aP)
G4GLOB_DLL std::ostream G4cerr
G4int GetBaryonNumber() const
CLHEP::HepLorentzVector G4LorentzVector
G4INCLXXInterfaceStore *const theInterfaceStore