33 #define INCLXX_IN_GEANT4_MODE 1
47 #include "G4String.hh"
56 thePreCompoundModel(aPreCompound),
58 complainedAboutBackupModel(false),
59 complainedAboutPreCompound(false),
62 if(!thePreCompoundModel) {
70 if(getenv(
"G4INCLXX_NO_DE_EXCITATION")) {
71 G4String message =
"de-excitation is completely disabled!";
83 if(getenv(
"G4INCLXX_DUMP_REMNANT"))
84 dumpRemnantInfo =
true;
86 dumpRemnantInfo =
false;
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)
145 return theBackupModelNucleon->
ApplyYourself(aTrack, theNucleus);
152 && theNucleus.
GetA_asInt() > theMaxProjMassINCL) {
153 if(!complainedAboutBackupModel) {
154 complainedAboutBackupModel =
true;
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) {
172 if(!complainedAboutPreCompound) {
173 complainedAboutPreCompound =
true;
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.";
182 return thePreCompoundModel->
ApplyYourself(aTrack, theNucleus);
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);
201 const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus);
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";
265 std::list<G4Fragment> remnants;
267 const G4int maxTries = 200;
275 const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack);
287 if(inverseKinematics) {
306 momentum *= toLabFrame;
308 if(inverseKinematics) {
309 momentum *= *toDirectKinematics;
315 fourMomentumOut += momentum;
319 G4String message =
"the model produced a particle that couldn't be converted to Geant4 particle.";
334 eventInfo.
jyRem[i]*hbar_Planck,
335 eventInfo.
jzRem[i]*hbar_Planck
339 const G4double scaling = remnant4MomentumScaling(nuclearMass,
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;
367 fourMomentumOut += fourMomentum;
370 if(dumpRemnantInfo) {
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) {
511 return theIonTable->
GetIon(Z, A, 0);
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);
Hep3Vector boostVector() const
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetMaxProjMassINCL() const
Getter for theMaxProjMassINCL.
G4double GetCascadeMinEnergyPerNucleon() const
Getter for cascadeMinEnergyPerNucleon.
G4HadSecondary * GetSecondary(size_t i)
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 [ ].
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
HepRotation & rotateY(double delta)
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].
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
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 precision function energy(A, Z)
Singleton class for configuring the INCL++ Geant4 interface.
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
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].
static G4Triton * Triton()
static G4Proton * Proton()
static G4PionPlus * PionPlus()
const G4String & GetParticleType() const
static G4Neutron * Neutron()
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
static G4Deuteron * Deuteron()
G4int GetAtomicMass() const
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
G4LorentzVector Get4Momentum() const
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
static G4HadronicInteractionRegistry * Instance()
G4DynamicParticle * GetParticle()
static G4PionMinus * PionMinus()
G4bool GetAccurateProjectile() const
Getter for accurateProjectile.
G4VPreCompoundModel * theDeExcitation
HepRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
std::map< G4String, G4AttDef > * GetInstance(const G4String &storeKey, G4bool &isNew)
void setVect(const Hep3Vector &)
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