Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCL::INCL Class Reference

#include <G4INCLCascade.hh>

Public Member Functions

 INCL (Config const *const config)
 
 ~INCL ()
 
 INCL (const INCL &rhs)
 Dummy copy constructor to silence Coverity warning.
 
INCLoperator= (const INCL &rhs)
 Dummy assignment operator to silence Coverity warning.
 
G4bool prepareReaction (const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S)
 
G4bool initializeTarget (const G4int A, const G4int Z, const G4int S, AnnihilationType theAType)
 
G4double initUniverseRadiusForAntiprotonAtRest (const G4int A, const G4int Z)
 
const EventInfoprocessEvent ()
 
const EventInfoprocessEvent (ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4int targetA, const G4int targetZ, const G4int targetS)
 
void finalizeGlobalInfo (Random::SeedVector const &initialSeeds)
 
const GlobalInfogetGlobalInfo () const
 

Detailed Description

Definition at line 56 of file G4INCLCascade.hh.

Constructor & Destructor Documentation

◆ INCL() [1/2]

G4INCL::INCL::INCL ( Config const *const config)

Definition at line 82 of file G4INCLCascade.cc.

83 :propagationModel(0), theA(208), theZ(82), theS(0),
84 targetInitSuccess(false),
85 maxImpactParameter(0.),
86 maxUniverseRadius(0.),
87 maxInteractionDistance(0.),
88 fixedImpactParameter(0.),
89 theConfig(config),
90 nucleus(NULL),
91 forceTransparent(false),
92 minRemnantSize(4)
93 {
94 // Set the logger object.
95#ifdef INCLXX_IN_GEANT4_MODE
97#else // INCLXX_IN_GEANT4_MODE
98 Logger::initialize(theConfig);
99#endif // INCLXX_IN_GEANT4_MODE
100
101 // Set the random number generator algorithm. The system can support
102 // multiple different generator algorithms in a completely
103 // transparent way.
104 Random::initialize(theConfig);
105
106 // Select the Pauli and CDPP blocking algorithms
107 Pauli::initialize(theConfig);
108
109 // Set the cross-section set
110 CrossSections::initialize(theConfig);
111
112 // Set the phase-space generator
114
115 // Select the Coulomb-distortion algorithm:
117
118 // Select the clustering algorithm:
119 Clustering::initialize(theConfig);
120
121 // Initialize the INCL particle table:
122 ParticleTable::initialize(theConfig);
123
124 // Initialize the value of cutNN in BinaryCollisionAvatar
126
127 // Initialize the value of strange cross section bias
129
130 // Propagation model is responsible for finding avatars and
131 // transporting the particles. In principle this step is "hidden"
132 // behind an abstract interface and the rest of the system does not
133 // care how the transportation and avatar finding is done. This
134 // should allow us to "easily" experiment with different avatar
135 // finding schemes and even to support things like curved
136 // trajectories in the future.
137 propagationModel = new StandardPropagationModel(theConfig->getLocalEnergyBBType(),theConfig->getLocalEnergyPiType(),theConfig->getHadronizationTime());
139 cascadeAction = new AvatarDumpAction();
140 else
141 cascadeAction = new CascadeAction();
142 cascadeAction->beforeRunAction(theConfig);
143
144 theGlobalInfo.cascadeModel = theConfig->getVersionString();
145 theGlobalInfo.deexcitationModel = theConfig->getDeExcitationString();
146#ifdef INCL_ROOT_USE
147 theGlobalInfo.rootSelection = theConfig->getROOTSelectionString();
148#endif
149
150#ifndef INCLXX_IN_GEANT4_MODE
151 // Fill in the global information
152 theGlobalInfo.At = theConfig->getTargetA();
153 theGlobalInfo.Zt = theConfig->getTargetZ();
154 theGlobalInfo.St = theConfig->getTargetS();
155 const ParticleSpecies theSpecies = theConfig->getProjectileSpecies();
156 theGlobalInfo.Ap = theSpecies.theA;
157 theGlobalInfo.Zp = theSpecies.theZ;
158 theGlobalInfo.Sp = theSpecies.theS;
159 theGlobalInfo.Ep = theConfig->getProjectileKineticEnergy();
160 theGlobalInfo.biasFactor = theConfig->getBias();
161#endif
162
163 fixedImpactParameter = theConfig->getImpactParameter();
164 }
static void setBias(const G4double b)
Set the global bias factor.
static void setCutNN(const G4double c)
static std::string const getVersionString()
Get the INCL version string.
G4double getImpactParameter() const
CascadeActionType getCascadeActionType() const
Get the cascade-action type.
G4int getTargetA() const
Get the target mass number.
G4double getProjectileKineticEnergy() const
Get the projectile kinetic energy.
G4double getHadronizationTime() const
Get the hadronization time.
ParticleSpecies getProjectileSpecies() const
Get the projectile species.
LocalEnergyType getLocalEnergyPiType() const
Get the type of local energy for pi-N and decay avatars.
G4int getTargetS() const
Get the target strangess number.
G4double getCutNN() const
LocalEnergyType getLocalEnergyBBType() const
Get the type of local energy for N-N avatars.
G4double getBias() const
Get the bias.
G4int getTargetZ() const
Get the target charge number.
std::string getDeExcitationString() const
Get the de-excitation string.
void initialize(Config const *const theConfig)
Initialize the clustering model based on the Config object.
void initialize(Config const *const theConfig)
Initialize the Coulomb-distortion algorithm.
void initialize(Config const *const theConfig)
void initVerbosityLevelFromEnvvar()
void initialize(Config const *const theConfig=0)
Initialize the particle table.
void initialize(Config const *const aConfig)
Initialise blockers according to a Config object.
void initialize(Config const *const theConfig)
void initialize(Config const *const)
Initialize generator according to a Config object.
std::string deexcitationModel
Name of the de-excitation model.
Short_t St
Target strangeness number given as input.
Float_t Ep
Projectile kinetic energy given as input.
Float_t biasFactor
Bias factor.
Short_t At
Target mass number given as input.
Short_t Zt
Target charge number given as input.
Short_t Ap
Projectile mass number given as input.
Short_t Sp
Projectile strangeness number given as input.
Short_t Zp
Projectile charge number given as input.
std::string cascadeModel
Name of the cascade model.

◆ ~INCL()

G4INCL::INCL::~INCL ( )

Definition at line 166 of file G4INCLCascade.cc.

166 {
168#ifndef INCLXX_IN_GEANT4_MODE
169 NuclearMassTable::deleteTable();
170#endif
177#ifndef INCLXX_IN_GEANT4_MODE
178 Logger::deleteLoggerSlave();
179#endif
182 cascadeAction->afterRunAction();
183 delete cascadeAction;
184 delete propagationModel;
185 delete theConfig;
186 }
static void deleteBackupParticles()
Release the memory allocated for the backup particles.
void deleteClusteringModel()
Delete the clustering model.
void deleteCoulomb()
Delete the Coulomb-distortion object.
void clearCache()
Clear the INuclearPotential cache.
void deleteBlockers()
Delete blockers.

◆ INCL() [2/2]

G4INCL::INCL::INCL ( const INCL & rhs)

Dummy copy constructor to silence Coverity warning.

Member Function Documentation

◆ finalizeGlobalInfo()

void G4INCL::INCL::finalizeGlobalInfo ( Random::SeedVector const & initialSeeds)

Definition at line 1036 of file G4INCLCascade.cc.

1036 {
1037 const G4double normalisationFactor = theGlobalInfo.geometricCrossSection /
1038 ((G4double) theGlobalInfo.nShots);
1039 theGlobalInfo.nucleonAbsorptionCrossSection = normalisationFactor *
1040 ((G4double) theGlobalInfo.nNucleonAbsorptions);
1041 theGlobalInfo.pionAbsorptionCrossSection = normalisationFactor *
1042 ((G4double) theGlobalInfo.nPionAbsorptions);
1043 theGlobalInfo.reactionCrossSection = normalisationFactor *
1044 ((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents));
1045 theGlobalInfo.errorReactionCrossSection = normalisationFactor *
1046 std::sqrt((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents));
1047 theGlobalInfo.forcedCNCrossSection = normalisationFactor *
1048 ((G4double) theGlobalInfo.nForcedCompoundNucleus);
1049 theGlobalInfo.errorForcedCNCrossSection = normalisationFactor *
1050 std::sqrt((G4double) (theGlobalInfo.nForcedCompoundNucleus));
1051 theGlobalInfo.completeFusionCrossSection = normalisationFactor *
1052 ((G4double) theGlobalInfo.nCompleteFusion);
1053 theGlobalInfo.errorCompleteFusionCrossSection = normalisationFactor *
1054 std::sqrt((G4double) (theGlobalInfo.nCompleteFusion));
1055 theGlobalInfo.energyViolationInteractionCrossSection = normalisationFactor *
1056 ((G4double) theGlobalInfo.nEnergyViolationInteraction);
1057
1058 theGlobalInfo.initialRandomSeeds.assign(initialSeeds.begin(), initialSeeds.end());
1059
1060 Random::SeedVector theSeeds = Random::getSeeds();
1061 theGlobalInfo.finalRandomSeeds.assign(theSeeds.begin(), theSeeds.end());
1062 }
double G4double
Definition G4Types.hh:83
SeedVector getSeeds()
Int_t nNucleonAbsorptions
Number of nucleon absorptions (no outcoming particles)
Float_t forcedCNCrossSection
Calculated forced-compound-nucleus cross section.
Float_t pionAbsorptionCrossSection
Pion absorption cross section.
Float_t errorCompleteFusionCrossSection
Error on the calculated complete-fusion cross section (nParticles==0)
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
std::vector< Int_t > initialRandomSeeds
Initial seeds for the pseudo-random-number generator.
Int_t nShots
Number of shots.
Float_t completeFusionCrossSection
Calculated complete-fusion cross section (nParticles==0)
Float_t nucleonAbsorptionCrossSection
Nucleon absorption cross section.
Float_t errorReactionCrossSection
Error on the calculated reaction cross section.
std::vector< Int_t > finalRandomSeeds
Final seeds for the pseudo-random-number generator.
Int_t nPionAbsorptions
Number of nucleon absorptions (no outcoming pions)
Float_t energyViolationInteractionCrossSection
Cross section for attempted collisions/decays for which the energy-conservation algorithm failed to f...
Int_t nCompleteFusion
Number of complete-fusion events (nParticles==0)
Int_t nTransparents
Number of transparent shots.
Int_t nForcedCompoundNucleus
Number of forced compound-nucleus events.
Float_t reactionCrossSection
Calculated reaction cross section.
Float_t errorForcedCNCrossSection
Error on the calculated forced-compound-nucleus cross section.
Float_t geometricCrossSection
Geometric cross section.

◆ getGlobalInfo()

const GlobalInfo & G4INCL::INCL::getGlobalInfo ( ) const
inline

Definition at line 90 of file G4INCLCascade.hh.

90{ return theGlobalInfo; }

◆ initializeTarget()

G4bool G4INCL::INCL::initializeTarget ( const G4int A,
const G4int Z,
const G4int S,
AnnihilationType theAType )

Definition at line 296 of file G4INCLCascade.cc.

296 {
297 delete nucleus;
298
299 if (theAType==PType || theAType==NType) {
300 G4double newmaxUniverseRadius=0.;
301 if (theAType==PType) newmaxUniverseRadius=initUniverseRadiusForAntiprotonAtRest(A+1, Z+1);
302 else newmaxUniverseRadius=initUniverseRadiusForAntiprotonAtRest(A+1, Z);
303 nucleus = new Nucleus(A, Z, S, theConfig, newmaxUniverseRadius, theAType);
304 }
305 else{
306 nucleus = new Nucleus(A, Z, S, theConfig, maxUniverseRadius, theAType);
307 }
308 nucleus->getStore()->getBook().reset();
309 nucleus->initializeParticles();
310 propagationModel->setNucleus(nucleus);
311 return true;
312 }
G4double S(G4double temp)
const G4double A[17]
G4double initUniverseRadiusForAntiprotonAtRest(const G4int A, const G4int Z)
virtual void setNucleus(G4INCL::Nucleus *nucleus)=0
Store * getStore() const
Book & getBook()

Referenced by prepareReaction().

◆ initUniverseRadiusForAntiprotonAtRest()

G4double G4INCL::INCL::initUniverseRadiusForAntiprotonAtRest ( const G4int A,
const G4int Z )

Definition at line 1164 of file G4INCLCascade.cc.

1164 {
1165 G4double rMax = 0.0;
1166 if(A==0) {
1167 IsotopicDistribution const &anIsotopicDistribution =
1169 IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
1170 for(IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
1171 const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, i->theA, Z);
1172 const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, i->theA, Z);
1173 const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
1174 rMax = std::max(maximumRadius, rMax);
1175 }
1176 } else {
1177 const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
1178 const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
1179 const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
1180 rMax = std::max(maximumRadius, rMax);
1181 }
1182 return rMax;
1183 }
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
IsotopeVector::iterator IsotopeIter
std::vector< Isotope > IsotopeVector

Referenced by initializeTarget().

◆ operator=()

INCL & G4INCL::INCL::operator= ( const INCL & rhs)

Dummy assignment operator to silence Coverity warning.

◆ prepareReaction()

G4bool G4INCL::INCL::prepareReaction ( const ParticleSpecies & projectileSpecies,
const G4double kineticEnergy,
const G4int A,
const G4int Z,
const G4int S )

Definition at line 188 of file G4INCLCascade.cc.

188 {
189 if(A < 0 || A > 300 || Z < 1 || Z > 200) {
190 INCL_ERROR("Unsupported target: A = " << A << " Z = " << Z << " S = " << S << '\n'
191 << "Target configuration rejected." << '\n');
192 return false;
193 }
194 if(projectileSpecies.theType==Composite &&
195 (projectileSpecies.theZ==projectileSpecies.theA || projectileSpecies.theZ==0)) {
196 INCL_ERROR("Unsupported projectile: A = " << projectileSpecies.theA << " Z = " << projectileSpecies.theZ << " S = " << projectileSpecies.theS << '\n'
197 << "Projectile configuration rejected." << '\n');
198 return false;
199 }
200
201 // Reset the forced-transparent flag
202 forceTransparent = false;
203
204 // Initialise the maximum universe radius
205 initUniverseRadius(projectileSpecies, kineticEnergy, A, Z);
206 // Initialise the nucleus
207
208//D
209 //reset
210 G4bool ProtonIsTheVictim = false;
211 G4bool NeutronIsTheVictim = false;
212 theEventInfo.annihilationP = false;
213 theEventInfo.annihilationN = false;
214
215 //G4double AnnihilationBarrier = kineticEnergy;
216 if(projectileSpecies.theType == antiProton && kineticEnergy <= theConfig->getAtrestThreshold()){
217 G4double SpOverSn = 1.331;//from experiments with deuteron (E.Klempt)
218 //INCL_WARN("theA number set to A-1 from " << A <<'\n');
219
220 G4double neutronprob;
221 if(theConfig->isNaturalTarget()){ // A = 0 in this case
222 theA = ParticleTable::drawRandomNaturalIsotope(Z) - 1; //43 and 61 are ok (Technetium and Promethium)
223 neutronprob = (theA + 1 - Z)/(theA + 1 - Z + SpOverSn*Z);
224 }
225 else{
226 theA = A - 1;
227 neutronprob = (A - Z)/(A - Z + SpOverSn*Z); //from experiments with deuteron (E.Klempt)
228 }
229
230 theS = S;
231
232 G4double rndm = Random::shoot();
233 if(rndm >= neutronprob){ //proton is annihilated
234 theEventInfo.annihilationP = true;
235 theZ = Z - 1;
236 ProtonIsTheVictim = true;
237 //INCL_WARN("theZ number set to Z-1 from " << Z << '\n');
238 }
239 else{ //neutron is annihilated
240 theEventInfo.annihilationN = true;
241 theZ = Z;
242 NeutronIsTheVictim = true;
243 }
244 }
245 else{ // not annihilation of pbar
246 theZ = Z;
247 theS = S;
248 if(theConfig->isNaturalTarget())
249 theA = ParticleTable::drawRandomNaturalIsotope(Z); //change order
250 else
251 theA = A;
252 }
253
254 AnnihilationType theAType = Def;
255 if(ProtonIsTheVictim == true && NeutronIsTheVictim == false)
256 theAType = PType;
257 if(NeutronIsTheVictim == true && ProtonIsTheVictim == false)
258 theAType = NType;
259
260//D
261
262 initializeTarget(theA, theZ, theS, theAType);
263
264 // Set the maximum impact parameter
265 maxImpactParameter = CoulombDistortion::maxImpactParameter(projectileSpecies, kineticEnergy, nucleus);
266 INCL_DEBUG("Maximum impact parameter initialised: " << maxImpactParameter << '\n');
267
268 // For forced CN events
269 initMaxInteractionDistance(projectileSpecies, kineticEnergy);
270// Set the geometric cross sectiony section
271 if(projectileSpecies.theType == antiProton && kineticEnergy <= theConfig->getAtrestThreshold()){
272 G4int currentA = A;
273 if(theConfig->isNaturalTarget()){
275 }
276 G4double kineticEnergy2=kineticEnergy;
277 if (kineticEnergy2 <= 0.) kineticEnergy2=0.001;
278 theGlobalInfo.geometricCrossSection = 9.7* //normalization factor from Corradini
279 Math::pi*std::pow((1.840 + 1.120*std::pow(currentA,(1./3.))),2)*
280 (1. + (Z*G4INCL::PhysicalConstants::eSquared*(currentA+1))/(currentA*kineticEnergy2*(1.840 + 1.120*std::pow(currentA,(1./3.)))));
281 //xsection formula was borrowed from Corradini et al. https://doi.org/10.1016/j.physletb.2011.09.069
282 }
283 else{
284 theGlobalInfo.geometricCrossSection =
285 Math::tenPi*std::pow(maxImpactParameter,2);
286 }
287
288 // Set the minimum remnant size
289 if(projectileSpecies.theA > 0)
290 minRemnantSize = std::min(theA, 4);
291 else
292 minRemnantSize = std::min(theA-1, 4);
293 return true;
294 }
#define INCL_ERROR(x)
#define INCL_DEBUG(x)
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
G4bool isNaturalTarget() const
Natural targets.
G4bool initializeTarget(const G4int A, const G4int Z, const G4int S, AnnihilationType theAType)
G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n)
Return the maximum impact parameter for Coulomb-distorted trajectories.
const G4double pi
const G4double tenPi
G4int drawRandomNaturalIsotope(const G4int Z)
const G4double eSquared
Coulomb conversion factor [MeV*fm].
G4double shoot()
Bool_t annihilationP
True if annihilation at rest on a proton.
Bool_t annihilationN
True if annihilation at rest on a neutron.

Referenced by processEvent().

◆ processEvent() [1/2]

const EventInfo & G4INCL::INCL::processEvent ( )
inline

Definition at line 72 of file G4INCLCascade.hh.

72 {
73 return processEvent(
74 theConfig->getProjectileSpecies(),
75 theConfig->getProjectileKineticEnergy(),
76 theConfig->getTargetA(),
77 theConfig->getTargetZ(),
78 theConfig->getTargetS()
79 );
80 }
const EventInfo & processEvent()

Referenced by G4INCLXXInterface::ApplyYourself(), and processEvent().

◆ processEvent() [2/2]

const EventInfo & G4INCL::INCL::processEvent ( ParticleSpecies const & projectileSpecies,
const G4double kineticEnergy,
const G4int targetA,
const G4int targetZ,
const G4int targetS )

Definition at line 314 of file G4INCLCascade.cc.

320 {
321
322 ParticleList starlistH2;
323
324 if (projectileSpecies.theType==antiProton && (targetA==1 || targetA==2) && targetZ==1 && targetS==0) {
325
326 if (targetA==1) {
327 preCascade_pbarH1(projectileSpecies, kineticEnergy);
328 } else {
329 preCascade_pbarH2(projectileSpecies, kineticEnergy);
330 theEventInfo.annihilationP = false;
331 theEventInfo.annihilationN = false;
332
333 G4double SpOverSn = 1.331; //from experiments with deuteron (E.Klempt)
334
335 ThreeVector dummy(0.,0.,0.);
336 G4double rndm = Random::shoot()*(SpOverSn+1);
337 if (rndm <= SpOverSn) { //proton is annihilated
338 theEventInfo.annihilationP = true;
339 Particle *p2 = new Particle(Neutron, dummy, dummy);
340 starlistH2.push_back(p2);
341 //delete p2;
342 } else { //neutron is annihilated
343 theEventInfo.annihilationN = true;
344 Particle *p2 = new Particle(Proton, dummy, dummy);
345 starlistH2.push_back(p2);
346 //delete p2;
347 }
348 }
349
350 // File names
351#ifdef INCLXX_IN_GEANT4_MODE
352 if (!G4FindDataDir("G4INCLDATA") ) {
354 ed << " Data missing: set environment variable G4INCLDATA\n"
355 << " to point to the directory containing data files needed\n"
356 << " by the INCL++ model" << G4endl;
357 G4Exception("G4INCLDataFile::readData()","rawppbarFS.dat, ...", FatalException, ed);
358 }
359 G4String dataPath0(G4FindDataDir("G4INCLDATA"));
360 G4String dataPathppbar(dataPath0 + "/rawppbarFS.dat");
361 G4String dataPathnpbar(dataPath0 + "/rawnpbarFS.dat");
362 G4String dataPathppbark(dataPath0 + "/rawppbarFSkaonic.dat");
363 G4String dataPathnpbark(dataPath0 + "/rawnpbarFSkaonic.dat");
364#else
365 G4String path;
366 if (theConfig) path = theConfig->getINCLXXDataFilePath();
367 G4String dataPathppbar(path + "/rawppbarFS.dat");
368 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
369 G4String dataPathnpbar(path + "/rawnpbarFS.dat");
370 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
371 G4String dataPathppbark(path + "/rawppbarFSkaonic.dat");
372 INCL_DEBUG("Reading https://doi.org/10.1016/j.physrep.2005.03.002 ppbar kaonic final states" << dataPathppbark << '\n');
373 G4String dataPathnpbark(path + "/rawnpbarFSkaonic.dat");
374 INCL_DEBUG("Reading https://doi.org/10.1007/BF02818764 and https://link.springer.com/article/10.1007/BF02754930 npbar kaonic final states" << dataPathnpbark << '\n');
375 #endif
376
377 //read probabilities and particle types from file
378 std::vector<G4double> probabilities; //will store each FS yield
379 std::vector<std::vector<G4String>> particle_types; //will store particle names
380 G4double sum = 0.0; //will contain a sum of probabilities of all FS in the file
381 G4double kaonicFSprob=0.05; //probability to kave kaonic FS
382
383 ParticleList starlist;
384 ThreeVector mommy; //momentum to be assigned later
385
386 G4double rdm = Random::shoot();
387 ThreeVector annihilationPosition(0.,0.,0.);
388 if (rdm < (1.-kaonicFSprob)) { // pionic FS was chosen
389 INCL_DEBUG("pionic pp final state chosen" << '\n');
390 sum = read_file(dataPathppbar, probabilities, particle_types);
391 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.88 normalize by the sum of probabilities in the file
392 //now get the line number in the file where the FS particles are stored:
393 G4int n = findStringNumber(rdm, probabilities)-1;
394 if ( n < 0 ) return theEventInfo;
395 for (G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++) {
396 if (particle_types[n][j] == "pi0") {
397 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
398 starlist.push_back(p);
399 } else if (particle_types[n][j] == "pi-") {
400 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
401 starlist.push_back(p);
402 } else if (particle_types[n][j] == "pi+") {
403 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
404 starlist.push_back(p);
405 } else if (particle_types[n][j] == "omega") {
406 Particle *p = new Particle(Omega, mommy, annihilationPosition);
407 starlist.push_back(p);
408 } else if (particle_types[n][j] == "eta") {
409 Particle *p = new Particle(Eta, mommy, annihilationPosition);
410 starlist.push_back(p);
411 } else if (particle_types[n][j] == "rho-") {
412 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
413 starlist.push_back(p);
414 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
415 starlist.push_back(pp);
416 } else if (particle_types[n][j] == "rho+") {
417 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
418 starlist.push_back(p);
419 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
420 starlist.push_back(pp);
421 } else if (particle_types[n][j] == "rho0") {
422 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
423 starlist.push_back(p);
424 Particle *pp = new Particle(PiPlus, mommy, annihilationPosition);
425 starlist.push_back(pp);
426 } else {
427 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
428 for (G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++) {
429#ifdef INCLXX_IN_GEANT4_MODE
430 G4cout << "gotcha! " << particle_types[n][jj] << G4endl;
431#else
432 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
433#endif
434 }
435#ifdef INCLXX_IN_GEANT4_MODE
436 G4cout << "Some non-existing FS particle detected when reading pbar FS files" << G4endl;
437#else
438 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
439#endif
440 }
441 }
442 } else {
443 INCL_DEBUG("kaonic pp final state chosen" << '\n');
444 sum = read_file(dataPathppbark, probabilities, particle_types);
445 rdm = ((1.-rdm)/kaonicFSprob)*sum; //2670 normalize by the sum of probabilities in the file
446 //now get the line number in the file where the FS particles are stored:
447 G4int n = findStringNumber(rdm, probabilities)-1;
448 if ( n < 0 ) return theEventInfo;
449 for (G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++) {
450 if (particle_types[n][j] == "pi0") {
451 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
452 starlist.push_back(p);
453 } else if (particle_types[n][j] == "pi-") {
454 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
455 starlist.push_back(p);
456 } else if (particle_types[n][j] == "pi+") {
457 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
458 starlist.push_back(p);
459 } else if (particle_types[n][j] == "omega") {
460 Particle *p = new Particle(Omega, mommy, annihilationPosition);
461 starlist.push_back(p);
462 } else if (particle_types[n][j] == "eta") {
463 Particle *p = new Particle(Eta, mommy, annihilationPosition);
464 starlist.push_back(p);
465 } else if (particle_types[n][j] == "K-") {
466 Particle *p = new Particle(KMinus, mommy, annihilationPosition);
467 starlist.push_back(p);
468 } else if (particle_types[n][j] == "K+") {
469 Particle *p = new Particle(KPlus, mommy, annihilationPosition);
470 starlist.push_back(p);
471 } else if (particle_types[n][j] == "K0") {
472 Particle *p = new Particle(KZero, mommy, annihilationPosition);
473 starlist.push_back(p);
474 } else if (particle_types[n][j] == "K0b") {
475 Particle *p = new Particle(KZeroBar, mommy, annihilationPosition);
476 starlist.push_back(p);
477 } else {
478 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
479 for (G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++) {
480#ifdef INCLXX_IN_GEANT4_MODE
481 G4cout << "gotcha! " << particle_types[n][jj] << G4endl;
482#else
483 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
484#endif
485 }
486#ifdef INCLXX_IN_GEANT4_MODE
487 G4cout << "Some non-existing FS particle detected when reading pbar FS files" << G4endl;
488#else
489 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
490#endif
491 }
492 }
493 }
494
495 //compute energies of mesons with a phase-space model
497 if (starlist.size() < 2) {
498 INCL_ERROR("should never happen, at least 2 final state particles!" << '\n');
499 } else if (starlist.size() == 2) {
500 ParticleIter first = starlist.begin();
501 ParticleIter last = std::next(first, 1);
502 G4double m1 = (*first)->getMass();
503 G4double m2 = (*last)->getMass();
504 G4double s = energyOfMesonStar*energyOfMesonStar;
505 G4double mom1 = std::sqrt(s/4. - (std::pow(m1,2) + std::pow(m2,2))/2. - std::pow(m1,2)*std::pow(m2,2)/s + (std::pow(m1,4) + 2.*std::pow(m1*m2,2) + std::pow(m2,4))/(4.*s));
506 ThreeVector momentello = Random::normVector(mom1);
507 (*first)->setMomentum(momentello);
508 (*first)->adjustEnergyFromMomentum();
509 (*last)->setMomentum(-momentello);
510 (*last)->adjustEnergyFromMomentum();
511 } else {
512 PhaseSpaceGenerator::generate(energyOfMesonStar, starlist);
513 }
514
515 if (targetA==1) postCascade_pbarH1(starlist);
516 else postCascade_pbarH2(starlist,starlistH2);
517
518 theGlobalInfo.nShots++;
519 return theEventInfo;
520 } // pbar on H1
521
522 // ReInitialize the bias vector
524 //Particle::INCLBiasVector.Clear();
526
527 // Set the target and the projectile
528 targetInitSuccess = prepareReaction(projectileSpecies, kineticEnergy, targetA, targetZ, targetS);
529
530 if(!targetInitSuccess) {
531 INCL_WARN("Target initialisation failed for A=" << targetA << ", Z=" << targetZ << ", S=" << targetS << '\n');
532 theEventInfo.transparent=true;
533 return theEventInfo;
534 }
535
536 cascadeAction->beforeCascadeAction(propagationModel);
537
538 const G4bool canRunCascade = preCascade(projectileSpecies, kineticEnergy);
539 if(canRunCascade) {
540 cascade();
541 postCascade(projectileSpecies, kineticEnergy);
542 cascadeAction->afterCascadeAction(nucleus);
543 }
544 updateGlobalInfo();
545 return theEventInfo;
546 }
const char * G4FindDataDir(const char *)
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
#define INCL_WARN(x)
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
std::string const & getINCLXXDataFilePath() const
Set the ABLAXX datafile path.
G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S)
static std::vector< G4double > INCLBiasVector
Time ordered vector of all bias applied.
static G4ThreadLocal G4int nextBiasedCollisionID
G4double getRealMass(const G4INCL::ParticleType t)
Get particle mass (in MeV/c^2)
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
ThreeVector normVector(G4double norm=1.)
ParticleList::const_iterator ParticleIter
Bool_t transparent
True if the event is transparent.

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