| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /** \file G4INCLCascade.cc
39 *
40 * INCL Cascade
41 */
42 #include "G4INCLCascade.hh"
43 #include "G4INCLRandom.hh"
44 #include "G4INCLStandardPropagationModel.hh"
45 #include "G4INCLParticleTable.hh"
46 #include "G4INCLParticle.hh"
47 #include "G4INCLNuclearMassTable.hh"
48 #include "G4INCLGlobalInfo.hh"
49 #include "G4INCLNucleus.hh"
50
51 #include "G4INCLPauliBlocking.hh"
52
53 #include "G4INCLCrossSections.hh"
54
55 #include "G4INCLPhaseSpaceGenerator.hh"
56
57 #include "G4INCLLogger.hh"
58 #include "G4INCLGlobals.hh"
59 #include "G4INCLNuclearDensityFactory.hh"
60
61 #include "G4INCLINuclearPotential.hh"
62
63 #include "G4INCLCoulombDistortion.hh"
64
65 #include "G4INCLClustering.hh"
66
67 #include "G4INCLIntersection.hh"
68
69 #include "G4INCLBinaryCollisionAvatar.hh"
70
71 #include "G4INCLCascadeAction.hh"
72 #include "G4INCLAvatarDumpAction.hh"
73
74 #include <cstring>
75 #include <cstdlib>
76 #include <numeric>
77
78 #include "G4INCLPbarAtrestEntryChannel.hh"
79
80 namespace G4INCL {
81
82 INCL::INCL(Config const * const config)
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
96 Logger::initVerbosityLevelFromEnvvar();
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
113 PhaseSpaceGenerator::initialize(theConfig);
114
115 // Select the Coulomb-distortion algorithm:
116 CoulombDistortion::initialize(theConfig);
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
125 BinaryCollisionAvatar::setCutNN(theConfig->getCutNN());
126
127 // Initialize the value of strange cross section bias
128 BinaryCollisionAvatar::setBias(theConfig->getBias());
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());
138 if(theConfig->getCascadeActionType() == AvatarDumpActionType)
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 }
165
166 INCL::~INCL() {
167 InteractionAvatar::deleteBackupParticles();
168 #ifndef INCLXX_IN_GEANT4_MODE
169 NuclearMassTable::deleteTable();
170 #endif
171 PhaseSpaceGenerator::deletePhaseSpaceGenerator();
172 CrossSections::deleteCrossSections();
173 Pauli::deleteBlockers();
174 CoulombDistortion::deleteCoulomb();
175 Random::deleteGenerator();
176 Clustering::deleteClusteringModel();
177 #ifndef INCLXX_IN_GEANT4_MODE
178 Logger::deleteLoggerSlave();
179 #endif
180 NuclearDensityFactory::clearCache();
181 NuclearPotential::clearCache();
182 cascadeAction->afterRunAction();
183 delete cascadeAction;
184 delete propagationModel;
185 delete theConfig;
186 }
187
188 G4bool INCL::prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S) {
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()){
274 currentA = ParticleTable::drawRandomNaturalIsotope(Z);
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 }
295
296 G4bool INCL::initializeTarget(const G4int A, const G4int Z, const G4int S, AnnihilationType theAType) {
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 }
313
314 const EventInfo &INCL::processEvent(
315 ParticleSpecies const &projectileSpecies,
316 const G4double kineticEnergy,
317 const G4int targetA,
318 const G4int targetZ,
319 const G4int targetS
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") ) {
353 G4ExceptionDescription ed;
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 const G4String& dataPath0(G4FindDataDir("G4INCLDATA"));
360 const G4String& dataPathppbar(dataPath0 + "/rawppbarFS.dat");
361 // const G4String& dataPathnpbar(dataPath0 + "/rawnpbarFS.dat"); // NOT used!
362 const G4String& dataPathppbark(dataPath0 + "/rawppbarFSkaonic.dat");
363 // const G4String& dataPathnpbark(dataPath0 + "/rawnpbarFSkaonic.dat"); // NOT used!
364 #else
365 std::string path;
366 if (theConfig) path = theConfig->getINCLXXDataFilePath();
367 const std::string& dataPathppbar(path + "/rawppbarFS.dat");
368 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
369 const std::string& dataPathnpbar(path + "/rawnpbarFS.dat");
370 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
371 const std::string& 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 const std::string& 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, std::move(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
496 G4double energyOfMesonStar=ParticleTable::getRealMass(Proton)+ParticleTable::getRealMass(antiProton);
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
523 Particle::INCLBiasVector.clear();
524 //Particle::INCLBiasVector.Clear();
525 Particle::nextBiasedCollisionID = 0;
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 }
547
548 G4bool INCL::preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
549 // Reset theEventInfo
550 theEventInfo.reset();
551
552 EventInfo::eventNumber++;
553
554 // Fill in the event information
555 theEventInfo.projectileType = projectileSpecies.theType;
556 theEventInfo.Ap = (Short_t)projectileSpecies.theA;
557 theEventInfo.Zp = (Short_t)projectileSpecies.theZ;
558 theEventInfo.Sp = (Short_t)projectileSpecies.theS;
559 theEventInfo.Ep = kineticEnergy;
560 theEventInfo.St = (Short_t)nucleus->getS();
561
562 if(nucleus->getAnnihilationType()==PType){
563 theEventInfo.annihilationP = true;
564 theEventInfo.At = (Short_t)nucleus->getA()+1;
565 theEventInfo.Zt = (Short_t)nucleus->getZ()+1;
566 }
567 else if(nucleus->getAnnihilationType()==NType){
568 theEventInfo.annihilationN = true;
569 theEventInfo.At = (Short_t)nucleus->getA()+1;
570 theEventInfo.Zt = (Short_t)nucleus->getZ();
571 }
572 else {
573 theEventInfo.At = (Short_t)nucleus->getA();
574 theEventInfo.Zt = (Short_t)nucleus->getZ();
575 }
576 // Do nothing below the Coulomb barrier
577 if(maxImpactParameter<=0.) {
578 // Fill in the event information
579 //Particle *pbar = new Particle;
580 //PbarAtrestEntryChannel *obj = new PbarAtrestEntryChannel(nucleus, pbar);
581 if(projectileSpecies.theType == antiProton && kineticEnergy <= theConfig->getAtrestThreshold()){ //D
582 INCL_DEBUG("at rest annihilation" << '\n');
583 //theEventInfo.transparent = false;
584 } else {
585 theEventInfo.transparent = true;
586 return false;
587 }
588 }
589
590
591 // Randomly draw an impact parameter or use a fixed value, depending on the
592 // Config option
593 G4double impactParameter, phi;
594 if(fixedImpactParameter<0.) {
595 impactParameter = maxImpactParameter * std::sqrt(Random::shoot0());
596 phi = Random::shoot() * Math::twoPi;
597 } else {
598 impactParameter = fixedImpactParameter;
599 phi = 0.;
600 }
601 INCL_DEBUG("Selected impact parameter: " << impactParameter << '\n');
602
603 // Fill in the event information
604 theEventInfo.impactParameter = impactParameter;
605
606 const G4double effectiveImpactParameter = propagationModel->shoot(projectileSpecies, kineticEnergy, impactParameter, phi);
607 if(effectiveImpactParameter < 0.) {
608 // Fill in the event information
609 theEventInfo.transparent = true;
610 return false;
611 }
612
613 // Fill in the event information
614 theEventInfo.transparent = false;
615 theEventInfo.effectiveImpactParameter = effectiveImpactParameter;
616
617 return true;
618 }
619
620 void INCL::cascade() {
621 FinalState *finalState = new FinalState;
622
623 unsigned long loopCounter = 0;
624 const unsigned long maxLoopCounter = 10000000;
625 do {
626 // Run book keeping actions that should take place before propagation:
627 cascadeAction->beforePropagationAction(propagationModel);
628
629 // Get the avatar with the smallest time and propagate particles
630 // to that point in time.
631 IAvatar *avatar = propagationModel->propagate(finalState);
632
633 finalState->reset();
634
635 // Run book keeping actions that should take place after propagation:
636 cascadeAction->afterPropagationAction(propagationModel, avatar);
637
638 if(avatar == 0) break; // No more avatars in the avatar list.
639
640 // Run book keeping actions that should take place before avatar:
641 cascadeAction->beforeAvatarAction(avatar, nucleus);
642
643 // Channel is responsible for calculating the outcome of the
644 // selected avatar. There are different kinds of channels. The
645 // class IChannel is, again, an abstract interface that defines
646 // the externally observable behavior of all interaction
647 // channels.
648 // The handling of the channel is transparent to the API.
649 // Final state tells what changed...
650 avatar->fillFinalState(finalState);
651 // Run book keeping actions that should take place after avatar:
652 cascadeAction->afterAvatarAction(avatar, nucleus, finalState);
653
654 // So now we must give this information to the nucleus
655 nucleus->applyFinalState(finalState);
656 // and now we are ready to process the next avatar!
657
658 delete avatar;
659
660 ++loopCounter;
661 } while(continueCascade() && loopCounter<maxLoopCounter); /* Loop checking, 10.07.2015, D.Mancusi */
662
663 delete finalState;
664 }
665
666 void INCL::postCascade(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy) {
667 // Fill in the event information
668 theEventInfo.stoppingTime = propagationModel->getCurrentTime();
669
670 // The event bias
671 theEventInfo.eventBias = (Double_t) Particle::getTotalBias();
672
673 // Forced CN?
674 if(!(projectileSpecies.theType==antiProton && kineticEnergy<=theConfig->getAtrestThreshold())){
675 if(nucleus->getTryCompoundNucleus()) {
676 INCL_DEBUG("Trying compound nucleus" << '\n');
677 makeCompoundNucleus();
678 theEventInfo.transparent = forceTransparent;
679 // Global checks of conservation laws
680 #ifndef INCLXX_IN_GEANT4_MODE
681 if(!theEventInfo.transparent) globalConservationChecks(true);
682 #endif
683 return;
684 }
685 }
686
687 if(!(projectileSpecies.theType==antiProton && kineticEnergy<=theConfig->getAtrestThreshold())){
688 theEventInfo.transparent = forceTransparent || nucleus->isEventTransparent();
689 }
690
691 if(theEventInfo.transparent) {
692 ProjectileRemnant * const projectileRemnant = nucleus->getProjectileRemnant();
693 if(projectileRemnant) {
694 // Clear the incoming list (particles will be deleted by the ProjectileRemnant)
695 nucleus->getStore()->clearIncoming();
696 } else {
697 // Delete particles in the incoming list
698 nucleus->getStore()->deleteIncoming();
699 }
700 } else {
701
702 // Check if the nucleus contains strange particles
703 theEventInfo.sigmasInside = nucleus->containsSigma();
704 theEventInfo.antikaonsInside = nucleus->containsAntiKaon();
705 theEventInfo.lambdasInside = nucleus->containsLambda();
706 theEventInfo.kaonsInside = nucleus->containsKaon();
707
708 // Capture antiKaons and Sigmas and produce Lambda instead
709 theEventInfo.absorbedStrangeParticle = nucleus->decayInsideStrangeParticles();
710
711 // Emit strange particles still inside the nucleus
712 nucleus->emitInsideStrangeParticles();
713 theEventInfo.emitKaon = nucleus->emitInsideKaon();
714
715 #ifdef INCLXX_IN_GEANT4_MODE
716 theEventInfo.emitLambda = nucleus->emitInsideLambda();
717 #endif // INCLXX_IN_GEANT4_MODE
718
719 // Check if the nucleus contains deltas
720 theEventInfo.deltasInside = nucleus->containsDeltas();
721
722 // Take care of any remaining deltas
723 theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
724 theEventInfo.forcedDeltasInside = nucleus->decayInsideDeltas();
725
726 // Take care of any remaining etas, omegas, neutral Sigmas and/or neutral kaons
727 G4double timeThreshold=theConfig->getDecayTimeThreshold();
728 theEventInfo.forcedPionResonancesOutside = nucleus->decayOutgoingPionResonances(timeThreshold);
729 nucleus->decayOutgoingSigmaZero(timeThreshold);
730 nucleus->decayOutgoingNeutralKaon();
731
732 // Apply Coulomb distortion, if appropriate
733 // Note that this will apply Coulomb distortion also on pions emitted by
734 // unphysical remnants (see decayInsideDeltas). This is at variance with
735 // what INCL4.6 does, but these events are (should be!) so rare that
736 // whatever we do doesn't (shouldn't!) make any noticeable difference.
737 CoulombDistortion::distortOut(nucleus->getStore()->getOutgoingParticles(), nucleus);
738
739 // If the normal cascade predicted complete fusion, use the tabulated
740 // masses to compute the excitation energy, the recoil, etc.
741 if(nucleus->getStore()->getOutgoingParticles().size()==0
742 && (!nucleus->getProjectileRemnant()
743 || nucleus->getProjectileRemnant()->getParticles().size()==0)) {
744
745 INCL_DEBUG("Cascade resulted in complete fusion, using realistic fusion kinematics" << '\n');
746
747 nucleus->useFusionKinematics();
748
749 if(nucleus->getExcitationEnergy()<0.) {
750 // Complete fusion is energetically impossible, return a transparent
751 INCL_WARN("Complete-fusion kinematics yields negative excitation energy, returning a transparent!" << '\n');
752 theEventInfo.transparent = true;
753 return;
754 }
755
756 } else { // Normal cascade here
757
758 // Set the excitation energy
759 nucleus->setExcitationEnergy(nucleus->computeExcitationEnergy());
760
761 // Make a projectile pre-fragment out of the geometrical and dynamical
762 // spectators
763 theEventInfo.nUnmergedSpectators = makeProjectileRemnant();
764
765 // Compute recoil momentum, energy and spin of the nucleus
766 if(nucleus->getA()==1 && minRemnantSize>1) {
767 INCL_ERROR("Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << '\n');
768 }
769 nucleus->computeRecoilKinematics();
770
771 #ifndef INCLXX_IN_GEANT4_MODE
772 // Global checks of conservation laws
773 globalConservationChecks(false);
774 #endif
775
776 // Make room for the remnant recoil by rescaling the energies of the
777 // outgoing particles.
778 if(nucleus->hasRemnant()) rescaleOutgoingForRecoil();
779
780 }
781
782 // Cluster decay
783 theEventInfo.clusterDecay = nucleus->decayOutgoingClusters() || nucleus->decayMe(); //D
784
785 #ifndef INCLXX_IN_GEANT4_MODE
786 // Global checks of conservation laws
787 globalConservationChecks(true);
788 #endif
789
790 // Fill the EventInfo structure
791 nucleus->fillEventInfo(&theEventInfo);
792
793 }
794 }
795
796 void INCL::makeCompoundNucleus() {
797 // If this is not a nucleus-nucleus collision, don't attempt to make a
798 // compound nucleus.
799 //
800 // Yes, even nucleon-nucleus collisions can lead to particles entering
801 // below the Fermi level. Take e.g. 1-MeV p + He4.
802 if(!nucleus->isNucleusNucleusCollision()) {
803 forceTransparent = true;
804 return;
805 }
806
807 // Reset the internal Nucleus variables
808 nucleus->getStore()->clearIncoming();
809 nucleus->getStore()->clearOutgoing();
810 nucleus->getProjectileRemnant()->reset();
811 nucleus->setA(theEventInfo.At);
812 nucleus->setZ(theEventInfo.Zt);
813
814 // CN kinematical variables
815 // Note: the CN orbital angular momentum is neglected in what follows. We
816 // should actually take it into account!
817 ThreeVector theCNMomentum = nucleus->getIncomingMomentum();
818 ThreeVector theCNSpin = nucleus->getIncomingAngularMomentum();
819 const G4double theTargetMass = ParticleTable::getTableMass(theEventInfo.At, theEventInfo.Zt, theEventInfo.St);
820 G4int theCNA=theEventInfo.At, theCNZ=theEventInfo.Zt, theCNS=theEventInfo.St;
821 Cluster * const theProjectileRemnant = nucleus->getProjectileRemnant();
822 G4double theCNEnergy = theTargetMass + theProjectileRemnant->getEnergy();
823
824 // Loop over the potential participants
825 ParticleList const &initialProjectileComponents = theProjectileRemnant->getParticles();
826 std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
827 // Shuffle the list of potential participants
828 std::shuffle(shuffledComponents.begin(), shuffledComponents.end(), Random::getAdapter());
829
830 G4bool success = true;
831 G4bool atLeastOneNucleonEntering = false;
832 for(std::vector<Particle*>::const_iterator p=shuffledComponents.begin(), e=shuffledComponents.end(); p!=e; ++p) {
833 // Skip particles that miss the interaction distance
834 Intersection intersectionInteractionDistance(IntersectionFactory::getEarlierTrajectoryIntersection(
835 (*p)->getPosition(),
836 (*p)->getPropagationVelocity(),
837 maxInteractionDistance));
838 if(!intersectionInteractionDistance.exists)
839 continue;
840
841 // Build an entry avatar for this nucleon
842 atLeastOneNucleonEntering = true;
843 ParticleEntryAvatar *theAvatar = new ParticleEntryAvatar(0.0, nucleus, *p);
844 nucleus->getStore()->addParticleEntryAvatar(theAvatar);
845 FinalState *fs = theAvatar->getFinalState();
846 nucleus->applyFinalState(fs);
847 FinalStateValidity validity = fs->getValidity();
848 delete fs;
849 switch(validity) {
850 case ValidFS:
851 case ParticleBelowFermiFS:
852 case ParticleBelowZeroFS:
853 // Add the particle to the CN
854 theCNA++;
855 theCNZ += (*p)->getZ();
856 theCNS += (*p)->getS();
857 break;
858 case PauliBlockedFS:
859 case NoEnergyConservationFS:
860 default:
861 success = false;
862 break;
863 }
864 }
865
866 if(!success || !atLeastOneNucleonEntering) {
867 INCL_DEBUG("No nucleon entering in forced CN, forcing a transparent" << '\n');
868 forceTransparent = true;
869 return;
870 }
871
872 // assert(theCNA==nucleus->getA());
873 // assert(theCNA<=theEventInfo.At+theEventInfo.Ap);
874 // assert(theCNZ<=theEventInfo.Zt+theEventInfo.Zp);
875 // assert(theCNS>=theEventInfo.St+theEventInfo.Sp);
876
877 // Update the kinematics of the CN
878 theCNEnergy -= theProjectileRemnant->getEnergy();
879 theCNMomentum -= theProjectileRemnant->getMomentum();
880
881 // Deal with the projectile remnant
882 nucleus->finalizeProjectileRemnant(propagationModel->getCurrentTime());
883
884 // Subtract the angular momentum of the projectile remnant
885 // assert(nucleus->getStore()->getOutgoingParticles().empty());
886 theCNSpin -= theProjectileRemnant->getAngularMomentum();
887
888 // Compute the excitation energy of the CN
889 const G4double theCNMass = ParticleTable::getTableMass(theCNA,theCNZ,theCNS);
890 const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.mag2();
891 if(theCNInvariantMassSquared<0.) {
892 // Negative invariant mass squared, return a transparent
893 forceTransparent = true;
894 return;
895 }
896 const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
897 if(theCNExcitationEnergy<0.) {
898 // Negative excitation energy, return a transparent
899 INCL_DEBUG("CN excitation energy is negative, forcing a transparent" << '\n'
900 << " theCNA = " << theCNA << '\n'
901 << " theCNZ = " << theCNZ << '\n'
902 << " theCNS = " << theCNS << '\n'
903 << " theCNEnergy = " << theCNEnergy << '\n'
904 << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
905 << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
906 << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
907 );
908 forceTransparent = true;
909 return;
910 } else {
911 // Positive excitation energy, can make a CN
912 INCL_DEBUG("CN excitation energy is positive, forcing a CN" << '\n'
913 << " theCNA = " << theCNA << '\n'
914 << " theCNZ = " << theCNZ << '\n'
915 << " theCNS = " << theCNS << '\n'
916 << " theCNEnergy = " << theCNEnergy << '\n'
917 << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
918 << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
919 << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
920 );
921 nucleus->setA(theCNA);
922 nucleus->setZ(theCNZ);
923 nucleus->setS(theCNS);
924 nucleus->setMomentum(theCNMomentum);
925 nucleus->setEnergy(theCNEnergy);
926 nucleus->setExcitationEnergy(theCNExcitationEnergy);
927 nucleus->setMass(theCNMass+theCNExcitationEnergy);
928 nucleus->setSpin(theCNSpin); // neglects any orbital angular momentum of the CN
929
930 // Take care of any remaining deltas
931 theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
932
933 // Take care of any remaining etas and/or omegas
934 G4double timeThreshold=theConfig->getDecayTimeThreshold();
935 theEventInfo.forcedPionResonancesOutside = nucleus->decayOutgoingPionResonances(timeThreshold);
936
937 // Take care of any remaining Kaons
938 theEventInfo.emitKaon = nucleus->emitInsideKaon();
939
940 // Cluster decay
941 theEventInfo.clusterDecay = nucleus->decayOutgoingClusters() || nucleus->decayMe(); //D
942
943 // Fill the EventInfo structure
944 nucleus->fillEventInfo(&theEventInfo);
945 }
946 }
947
948 void INCL::rescaleOutgoingForRecoil() {
949 RecoilCMFunctor theRecoilFunctor(nucleus, theEventInfo);
950
951 // Apply the root-finding algorithm
952 const RootFinder::Solution theSolution = RootFinder::solve(&theRecoilFunctor, 1.0);
953 if(theSolution.success) {
954 theRecoilFunctor(theSolution.x); // Apply the solution
955 } else {
956 INCL_WARN("Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." << '\n');
957 }
958 }
959
960 #ifndef INCLXX_IN_GEANT4_MODE
961 void INCL::globalConservationChecks(G4bool afterRecoil) {
962 Nucleus::ConservationBalance theBalance = nucleus->getConservationBalance(theEventInfo,afterRecoil);
963
964 // Global conservation checks
965 const G4double pLongBalance = theBalance.momentum.getZ();
966 const G4double pTransBalance = theBalance.momentum.perp();
967 if(theBalance.Z != 0) {
968 INCL_ERROR("Violation of charge conservation! ZBalance = " << theBalance.Z << " eventNumber=" << theEventInfo.eventNumber << '\n');
969 }
970 if(theBalance.A != 0) {
971 INCL_ERROR("Violation of baryon-number conservation! ABalance = " << theBalance.A << " Emit Lambda=" << theEventInfo.emitLambda << " eventNumber=" << theEventInfo.eventNumber << '\n');
972 }
973 if(theBalance.S != 0) {
974 INCL_ERROR("Violation of strange-number conservation! SBalance = " << theBalance.S << " eventNumber=" << theEventInfo.eventNumber << '\n');
975 }
976 G4double EThreshold, pLongThreshold, pTransThreshold;
977 if(afterRecoil) {
978 // Less stringent checks after accommodating recoil
979 EThreshold = 10.; // MeV
980 pLongThreshold = 1.; // MeV/c
981 pTransThreshold = 1.; // MeV/c
982 } else {
983 // More stringent checks before accommodating recoil
984 EThreshold = 0.1; // MeV
985 pLongThreshold = 0.1; // MeV/c
986 pTransThreshold = 0.1; // MeV/c
987 }
988 if(std::abs(theBalance.energy)>EThreshold) {
989 INCL_WARN("Violation of energy conservation > " << EThreshold << " MeV. EBalance = " << theBalance.energy << " Emit Lambda=" << theEventInfo.emitLambda << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
990 }
991 if(std::abs(pLongBalance)>pLongThreshold) {
992 INCL_WARN("Violation of longitudinal momentum conservation > " << pLongThreshold << " MeV/c. pLongBalance = " << pLongBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
993 }
994 if(std::abs(pTransBalance)>pTransThreshold) {
995 INCL_WARN("Violation of transverse momentum conservation > " << pTransThreshold << " MeV/c. pTransBalance = " << pTransBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
996 }
997
998 // Feed the EventInfo variables
999 theEventInfo.EBalance = theBalance.energy;
1000 theEventInfo.pLongBalance = pLongBalance;
1001 theEventInfo.pTransBalance = pTransBalance;
1002 }
1003 #endif
1004
1005 G4bool INCL::continueCascade() {
1006 // Stop if we have passed the stopping time
1007 if(propagationModel->getCurrentTime() > propagationModel->getStoppingTime()) {
1008 INCL_DEBUG("Cascade time (" << propagationModel->getCurrentTime()
1009 << ") exceeded stopping time (" << propagationModel->getStoppingTime()
1010 << "), stopping cascade" << '\n');
1011 return false;
1012 }
1013 // Stop if there are no participants and no pions inside the nucleus
1014 if(nucleus->getStore()->getBook().getCascading()==0 &&
1015 nucleus->getStore()->getIncomingParticles().empty()) {
1016 INCL_DEBUG("No participants in the nucleus and no incoming particles left, stopping cascade" << '\n');
1017 return false;
1018 }
1019 // Stop if the remnant is smaller than minRemnantSize
1020 if(nucleus->getA() <= minRemnantSize) {
1021 INCL_DEBUG("Remnant size (" << nucleus->getA()
1022 << ") smaller than or equal to minimum (" << minRemnantSize
1023 << "), stopping cascade" << '\n');
1024 return false;
1025 }
1026 // Stop if we have to try and make a compound nucleus or if we have to
1027 // force a transparent
1028 if(nucleus->getTryCompoundNucleus()) {
1029 INCL_DEBUG("Trying to make a compound nucleus, stopping cascade" << '\n');
1030 return false;
1031 }
1032
1033 return true;
1034 }
1035
1036 void INCL::finalizeGlobalInfo(Random::SeedVector const &initialSeeds) {
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 }
1063
1064 G4int INCL::makeProjectileRemnant() {
1065 // Do nothing if this is not a nucleus-nucleus reaction
1066 if(!nucleus->getProjectileRemnant())
1067 return 0;
1068
1069 // Get the spectators (geometrical+dynamical) from the Store
1070 ParticleList geomSpectators(nucleus->getProjectileRemnant()->getParticles());
1071 ParticleList dynSpectators(nucleus->getStore()->extractDynamicalSpectators());
1072
1073 G4int nUnmergedSpectators = 0;
1074
1075 // If there are no spectators, do nothing
1076 if(dynSpectators.empty() && geomSpectators.empty()) {
1077 return 0;
1078 } else if(dynSpectators.size()==1 && geomSpectators.empty()) {
1079 // No geometrical spectators, one dynamical spectator
1080 // Just put it back in the outgoing list
1081 nucleus->getStore()->addToOutgoing(dynSpectators.front());
1082 } else {
1083 // Make a cluster out of the geometrical spectators
1084 ProjectileRemnant *theProjectileRemnant = nucleus->getProjectileRemnant();
1085
1086 // Add the dynamical spectators to the bunch
1087 ParticleList rejected = theProjectileRemnant->addAllDynamicalSpectators(dynSpectators);
1088 // Put back the rejected spectators into the outgoing list
1089 nUnmergedSpectators = (G4int)rejected.size();
1090 nucleus->getStore()->addToOutgoing(rejected);
1091
1092 // Deal with the projectile remnant
1093 nucleus->finalizeProjectileRemnant(propagationModel->getCurrentTime());
1094
1095 }
1096
1097 return nUnmergedSpectators;
1098 }
1099
1100 void INCL::initMaxInteractionDistance(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
1101 if(projectileSpecies.theType != Composite) {
1102 maxInteractionDistance = 0.;
1103 return;
1104 }
1105
1106 const G4double r0 = std::max(ParticleTable::getNuclearRadius(Proton, theA, theZ),
1107 ParticleTable::getNuclearRadius(Neutron, theA, theZ));
1108
1109 const G4double theNNDistance = CrossSections::interactionDistanceNN(projectileSpecies, kineticEnergy);
1110 maxInteractionDistance = r0 + theNNDistance;
1111 INCL_DEBUG("Initialised interaction distance: r0 = " << r0 << '\n'
1112 << " theNNDistance = " << theNNDistance << '\n'
1113 << " maxInteractionDistance = " << maxInteractionDistance << '\n');
1114 }
1115
1116 void INCL::initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z) {
1117 G4double rMax = 0.0;
1118 if(A==0) {
1119 IsotopicDistribution const &anIsotopicDistribution =
1120 ParticleTable::getNaturalIsotopicDistribution(Z);
1121 IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
1122 for(IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
1123 const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, i->theA, Z);
1124 const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, i->theA, Z);
1125 const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
1126 rMax = std::max(maximumRadius, rMax);
1127 }
1128 } else {
1129 const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
1130 const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
1131 const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
1132 rMax = std::max(maximumRadius, rMax);
1133 }
1134 if(p.theType==Composite || p.theType==Proton || p.theType==Neutron) {
1135 const G4double interactionDistanceNN = CrossSections::interactionDistanceNN(p, kineticEnergy);
1136 maxUniverseRadius = rMax + interactionDistanceNN;
1137 } else if(p.theType==PiPlus
1138 || p.theType==PiZero
1139 || p.theType==PiMinus) {
1140 const G4double interactionDistancePiN = CrossSections::interactionDistancePiN(kineticEnergy);
1141 maxUniverseRadius = rMax + interactionDistancePiN;
1142 } else if(p.theType==KPlus
1143 || p.theType==KZero) {
1144 const G4double interactionDistanceKN = CrossSections::interactionDistanceKN(kineticEnergy);
1145 maxUniverseRadius = rMax + interactionDistanceKN;
1146 } else if(p.theType==KZeroBar
1147 || p.theType==KMinus) {
1148 const G4double interactionDistanceKbarN = CrossSections::interactionDistanceKbarN(kineticEnergy);
1149 maxUniverseRadius = rMax + interactionDistanceKbarN;
1150 } else if(p.theType==Lambda
1151 ||p.theType==SigmaPlus
1152 || p.theType==SigmaZero
1153 || p.theType==SigmaMinus) {
1154 const G4double interactionDistanceYN = CrossSections::interactionDistanceYN(kineticEnergy);
1155 maxUniverseRadius = rMax + interactionDistanceYN;
1156 }
1157 else if(p.theType==antiProton) {
1158 maxUniverseRadius = rMax; //check interaction distance!!!
1159 }
1160 INCL_DEBUG("Initialised universe radius: " << maxUniverseRadius << '\n');
1161 }
1162
1163
1164 G4double INCL::initUniverseRadiusForAntiprotonAtRest(const G4int A, const G4int Z) {
1165 G4double rMax = 0.0;
1166 if(A==0) {
1167 IsotopicDistribution const &anIsotopicDistribution =
1168 ParticleTable::getNaturalIsotopicDistribution(Z);
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 }
1184
1185
1186 void INCL::updateGlobalInfo() {
1187 // Increment the global counter for the number of shots
1188 theGlobalInfo.nShots++;
1189
1190 if(theEventInfo.transparent) {
1191 // Increment the global counter for the number of transparents
1192 theGlobalInfo.nTransparents++;
1193 // Increment the global counter for the number of forced transparents
1194 if(forceTransparent)
1195 theGlobalInfo.nForcedTransparents++;
1196 return;
1197 }
1198
1199 // Check if we have an absorption:
1200 if(theEventInfo.nucleonAbsorption) theGlobalInfo.nNucleonAbsorptions++;
1201 if(theEventInfo.pionAbsorption) theGlobalInfo.nPionAbsorptions++;
1202
1203 // Count complete-fusion events
1204 if(theEventInfo.nCascadeParticles==0) theGlobalInfo.nCompleteFusion++;
1205
1206 if(nucleus->getTryCompoundNucleus())
1207 theGlobalInfo.nForcedCompoundNucleus++;
1208
1209 // Counters for the number of violations of energy conservation in
1210 // collisions
1211 theGlobalInfo.nEnergyViolationInteraction += theEventInfo.nEnergyViolationInteraction;
1212 }
1213
1214 G4double INCL::read_file(std::string filename, std::vector<G4double>& probabilities,
1215 std::vector<std::vector<G4String>>& particle_types) {
1216 std::ifstream file(filename);
1217 G4double sum_probs = 0.0;
1218 if (file.is_open()) {
1219 G4String line;
1220 while (getline(file, line)) {
1221 std::istringstream iss(line);
1222 G4double prob;
1223 iss >> prob;
1224 sum_probs += prob;
1225 probabilities.push_back(prob);
1226 std::vector<G4String> types;
1227 G4String type;
1228 while (iss >> type) {
1229 types.push_back(type);
1230 }
1231 particle_types.push_back(std::move(types));
1232 }
1233 } else {
1234 #ifdef INCLXX_IN_GEANT4_MODE
1235 G4cout << "ERROR no fread_file " << filename << G4endl;
1236 #else
1237 std::cout << "ERROR no fread_file " << filename << std::endl;
1238 #endif
1239 }
1240 return sum_probs;
1241 }
1242
1243
1244 G4int INCL::findStringNumber(G4double rdm, std::vector<G4double> yields) {
1245 G4int stringNumber = -1;
1246 G4double smallestsum = 0.0;
1247 G4double biggestsum = yields[0];
1248 //G4cout << "initial input " << rdm << G4endl;
1249 for (G4int i = 0; i < static_cast<G4int>(yields.size()-1); i++) {
1250 if (rdm >= smallestsum && rdm <= biggestsum) {
1251 //G4cout << smallestsum << " and " << biggestsum << G4endl;
1252 stringNumber = i+1;
1253 }
1254 smallestsum += yields[i];
1255 biggestsum += yields[i+1];
1256 }
1257 if(stringNumber==-1) stringNumber = static_cast<G4int>(yields.size());
1258 if(stringNumber==-1){
1259 INCL_ERROR("ERROR in findStringNumber (stringNumber=-1)");
1260 #ifdef INCLXX_IN_GEANT4_MODE
1261 G4cout << "ERROR in findStringNumber" << G4endl;
1262 #else
1263 std::cout << "ERROR in findStringNumber" << std::endl;
1264 #endif
1265 }
1266 return stringNumber;
1267 }
1268
1269
1270 void INCL::preCascade_pbarH1(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
1271 // Reset theEventInfo
1272 theEventInfo.reset();
1273
1274 EventInfo::eventNumber++;
1275
1276 // Fill in the event information
1277 theEventInfo.projectileType = projectileSpecies.theType;
1278 theEventInfo.Ap = -1;
1279 theEventInfo.Zp = -1;
1280 theEventInfo.Sp = 0;
1281 theEventInfo.Ep = kineticEnergy;
1282 theEventInfo.St = 0;
1283 theEventInfo.At = 1;
1284 theEventInfo.Zt = 1;
1285 }
1286
1287
1288 void INCL::postCascade_pbarH1(ParticleList const &outgoingParticles) {
1289 theEventInfo.nParticles = 0;
1290
1291 // Reset the remnant counter
1292 theEventInfo.nRemnants = 0;
1293 theEventInfo.history.clear();
1294
1295 for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
1296 theEventInfo.A[theEventInfo.nParticles] = (Short_t)(*i)->getA();
1297 theEventInfo.Z[theEventInfo.nParticles] = (Short_t)(*i)->getZ();
1298 theEventInfo.S[theEventInfo.nParticles] = (Short_t)(*i)->getS();
1299 theEventInfo.EKin[theEventInfo.nParticles] = (*i)->getKineticEnergy();
1300 ThreeVector mom = (*i)->getMomentum();
1301 theEventInfo.px[theEventInfo.nParticles] = mom.getX();
1302 theEventInfo.py[theEventInfo.nParticles] = mom.getY();
1303 theEventInfo.pz[theEventInfo.nParticles] = mom.getZ();
1304 theEventInfo.theta[theEventInfo.nParticles] = Math::toDegrees(mom.theta());
1305 theEventInfo.phi[theEventInfo.nParticles] = Math::toDegrees(mom.phi());
1306 theEventInfo.origin[theEventInfo.nParticles] = -1;
1307 #ifdef INCLXX_IN_GEANT4_MODE
1308 theEventInfo.parentResonancePDGCode[theEventInfo.nParticles] = (*i)->getParentResonancePDGCode();
1309 theEventInfo.parentResonanceID[theEventInfo.nParticles] = (*i)->getParentResonanceID();
1310 #endif
1311 theEventInfo.history.push_back("");
1312 ParticleSpecies pt((*i)->getType());
1313 theEventInfo.PDGCode[theEventInfo.nParticles] = pt.getPDGCode();
1314 theEventInfo.nParticles++;
1315 }
1316 theEventInfo.nCascadeParticles = theEventInfo.nParticles;
1317 }
1318
1319
1320 void INCL::preCascade_pbarH2(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
1321 // Reset theEventInfo
1322 theEventInfo.reset();
1323
1324 EventInfo::eventNumber++;
1325
1326 // Fill in the event information
1327 theEventInfo.projectileType = projectileSpecies.theType;
1328 theEventInfo.Ap = -1;
1329 theEventInfo.Zp = -1;
1330 theEventInfo.Sp = 0;
1331 theEventInfo.Ep = kineticEnergy;
1332 theEventInfo.St = 0;
1333 theEventInfo.At = 2;
1334 theEventInfo.Zt = 1;
1335 }
1336
1337
1338 void INCL::postCascade_pbarH2(ParticleList const &outgoingParticles, ParticleList const &H2Particles) {
1339 theEventInfo.nParticles = 0;
1340
1341 // Reset the remnant counter
1342 theEventInfo.nRemnants = 0;
1343 theEventInfo.history.clear();
1344
1345 for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
1346 theEventInfo.A[theEventInfo.nParticles] = (Short_t)(*i)->getA();
1347 theEventInfo.Z[theEventInfo.nParticles] = (Short_t)(*i)->getZ();
1348 theEventInfo.S[theEventInfo.nParticles] = (Short_t)(*i)->getS();
1349 theEventInfo.EKin[theEventInfo.nParticles] = (*i)->getKineticEnergy();
1350 ThreeVector mom = (*i)->getMomentum();
1351 theEventInfo.px[theEventInfo.nParticles] = mom.getX();
1352 theEventInfo.py[theEventInfo.nParticles] = mom.getY();
1353 theEventInfo.pz[theEventInfo.nParticles] = mom.getZ();
1354 theEventInfo.theta[theEventInfo.nParticles] = Math::toDegrees(mom.theta());
1355 theEventInfo.phi[theEventInfo.nParticles] = Math::toDegrees(mom.phi());
1356 theEventInfo.origin[theEventInfo.nParticles] = -1;
1357 #ifdef INCLXX_IN_GEANT4_MODE
1358 theEventInfo.parentResonancePDGCode[theEventInfo.nParticles] = (*i)->getParentResonancePDGCode();
1359 theEventInfo.parentResonanceID[theEventInfo.nParticles] = (*i)->getParentResonanceID();
1360 #endif
1361 theEventInfo.history.push_back("");
1362 ParticleSpecies pt((*i)->getType());
1363 theEventInfo.PDGCode[theEventInfo.nParticles] = pt.getPDGCode();
1364 theEventInfo.nParticles++;
1365 }
1366
1367 for(ParticleIter i=H2Particles.begin(), e=H2Particles.end(); i!=e; ++i ) {
1368 theEventInfo.A[theEventInfo.nParticles] = (Short_t)(*i)->getA();
1369 theEventInfo.Z[theEventInfo.nParticles] = (Short_t)(*i)->getZ();
1370 theEventInfo.S[theEventInfo.nParticles] = (Short_t)(*i)->getS();
1371 theEventInfo.EKin[theEventInfo.nParticles] = (*i)->getKineticEnergy();
1372 ThreeVector mom = (*i)->getMomentum();
1373 theEventInfo.px[theEventInfo.nParticles] = mom.getX();
1374 theEventInfo.py[theEventInfo.nParticles] = mom.getY();
1375 theEventInfo.pz[theEventInfo.nParticles] = mom.getZ();
1376 theEventInfo.theta[theEventInfo.nParticles] = Math::toDegrees(mom.theta());
1377 theEventInfo.phi[theEventInfo.nParticles] = Math::toDegrees(mom.phi());
1378 theEventInfo.origin[theEventInfo.nParticles] = -1;
1379 #ifdef INCLXX_IN_GEANT4_MODE
1380 theEventInfo.parentResonancePDGCode[theEventInfo.nParticles] = (*i)->getParentResonancePDGCode();
1381 theEventInfo.parentResonanceID[theEventInfo.nParticles] = (*i)->getParentResonanceID();
1382 #endif
1383 theEventInfo.history.push_back("");
1384 ParticleSpecies pt((*i)->getType());
1385 theEventInfo.PDGCode[theEventInfo.nParticles] = pt.getPDGCode();
1386 theEventInfo.nParticles++;
1387 }
1388 theEventInfo.nCascadeParticles = theEventInfo.nParticles;
1389 }
1390
1391 }
1392