| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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),
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 algorit
102 // multiple different generator algorithms
103 // transparent way.
104 Random::initialize(theConfig);
105
106 // Select the Pauli and CDPP blocking algo
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 Binary
125 BinaryCollisionAvatar::setCutNN(theConfig-
126
127 // Initialize the value of strange cross s
128 BinaryCollisionAvatar::setBias(theConfig->
129
130 // Propagation model is responsible for fi
131 // transporting the particles. In principl
132 // behind an abstract interface and the re
133 // care how the transportation and avatar
134 // should allow us to "easily" experiment
135 // finding schemes and even to support thi
136 // trajectories in the future.
137 propagationModel = new StandardPropagation
138 if(theConfig->getCascadeActionType() == Av
139 cascadeAction = new AvatarDumpAction();
140 else
141 cascadeAction = new CascadeAction();
142 cascadeAction->beforeRunAction(theConfig);
143
144 theGlobalInfo.cascadeModel = theConfig->ge
145 theGlobalInfo.deexcitationModel = theConfi
146 #ifdef INCL_ROOT_USE
147 theGlobalInfo.rootSelection = theConfig->g
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 = theConf
156 theGlobalInfo.Ap = theSpecies.theA;
157 theGlobalInfo.Zp = theSpecies.theZ;
158 theGlobalInfo.Sp = theSpecies.theS;
159 theGlobalInfo.Ep = theConfig->getProjectil
160 theGlobalInfo.biasFactor = theConfig->getB
161 #endif
162
163 fixedImpactParameter = theConfig->getImpac
164 }
165
166 INCL::~INCL() {
167 InteractionAvatar::deleteBackupParticles()
168 #ifndef INCLXX_IN_GEANT4_MODE
169 NuclearMassTable::deleteTable();
170 #endif
171 PhaseSpaceGenerator::deletePhaseSpaceGener
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 ParticleS
189 if(A < 0 || A > 300 || Z < 1 || Z > 200) {
190 INCL_ERROR("Unsupported target: A = " <<
191 << "Target configuration reje
192 return false;
193 }
194 if(projectileSpecies.theType==Composite &&
195 (projectileSpecies.theZ==projectileSpec
196 INCL_ERROR("Unsupported projectile: A =
197 << "Projectile configuration
198 return false;
199 }
200
201 // Reset the forced-transparent flag
202 forceTransparent = false;
203
204 // Initialise the maximum universe radius
205 initUniverseRadius(projectileSpecies, kine
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 = kineticEn
216 if(projectileSpecies.theType == antiProton
217 G4double SpOverSn = 1.331;//from experim
218 //INCL_WARN("theA number set to A-1 from
219
220 G4double neutronprob;
221 if(theConfig->isNaturalTarget()){ // A =
222 theA = ParticleTable::drawRandomNatura
223 neutronprob = (theA + 1 - Z)/(theA + 1
224 }
225 else{
226 theA = A - 1;
227 neutronprob = (A - Z)/(A - Z + SpOverS
228 }
229
230 theS = S;
231
232 G4double rndm = Random::shoot();
233 if(rndm >= neutronprob){ //proton is
234 theEventInfo.annihilationP = true;
235 theZ = Z - 1;
236 ProtonIsTheVictim = true;
237 //INCL_WARN("theZ number set to Z-1 fr
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::drawRandomNatura
250 else
251 theA = A;
252 }
253
254 AnnihilationType theAType = Def;
255 if(ProtonIsTheVictim == true && NeutronIsT
256 theAType = PType;
257 if(NeutronIsTheVictim == true && ProtonIsT
258 theAType = NType;
259
260 //D
261
262 initializeTarget(theA, theZ, theS, theATyp
263
264 // Set the maximum impact parameter
265 maxImpactParameter = CoulombDistortion::ma
266 INCL_DEBUG("Maximum impact parameter initi
267
268 // For forced CN events
269 initMaxInteractionDistance(projectileSpeci
270 // Set the geometric cross sectiony section
271 if(projectileSpecies.theType == antiProton
272 G4int currentA = A;
273 if(theConfig->isNaturalTarget()){
274 currentA = ParticleTable::drawRandomNa
275 }
276 G4double kineticEnergy2=kineticEnergy;
277 if (kineticEnergy2 <= 0.) kineticEnergy2
278 theGlobalInfo.geometricCrossSection = 9.
279 Math::pi*std::pow((1.840 + 1.120*std::
280 (1. + (Z*G4INCL::PhysicalConstants::eS
281 //xsection formula was borrowed from
282 }
283 else{
284 theGlobalInfo.geometricCrossSection =
285 Math::tenPi*std::pow(maxImpactParamete
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,
297 delete nucleus;
298
299 if (theAType==PType || theAType==NType) {
300 G4double newmaxUniverseRadius=0.;
301 if (theAType==PType) newmaxUniverseRadiu
302 else newmaxUniverseRadius=initUniverseRa
303 nucleus = new Nucleus(A, Z, S, theConfig
304 }
305 else{
306 nucleus = new Nucleus(A, Z, S, theConfig
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
325
326 if (targetA==1) {
327 preCascade_pbarH1(projectileSpecies, k
328 } else {
329 preCascade_pbarH2(projectileSpecies, k
330 theEventInfo.annihilationP = false;
331 theEventInfo.annihilationN = false;
332
333 G4double SpOverSn = 1.331; //from exp
334
335 ThreeVector dummy(0.,0.,0.);
336 G4double rndm = Random::shoot()*(SpOve
337 if (rndm <= SpOverSn) { //proton is a
338 theEventInfo.annihilationP = true;
339 Particle *p2 = new Particle(Neutron,
340 starlistH2.push_back(p2);
341 //delete p2;
342 } else { //neutron is
343 theEventInfo.annihilationN = true;
344 Particle *p2 = new Particle(Proton,
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
355 << " to point to the directory cont
356 << " by the INCL++ model" << G4endl
357 G4Exception("G4INCLDataFile::readData(
358 }
359 const G4String& dataPath0(G4FindDataDir(
360 const G4String& dataPathppbar(dataPath0
361 // const G4String& dataPathnpbar(dataPath
362 const G4String& dataPathppbark(dataPath0
363 // const G4String& dataPathnpbark(dataPat
364 #else
365 std::string path;
366 if (theConfig) path = theConfig->getINCL
367 const std::string& dataPathppbar(path +
368 INCL_DEBUG("Reading https://doi.org/10.1
369 const std::string& dataPathnpbar(path +
370 INCL_DEBUG("Reading https://doi.org/10.1
371 const std::string& dataPathppbark(path +
372 INCL_DEBUG("Reading https://doi.org/10.1
373 const std::string& dataPathnpbark(path +
374 INCL_DEBUG("Reading https://doi.org/10.1
375 #endif
376
377 //read probabilities and particle types
378 std::vector<G4double> probabilities; //
379 std::vector<std::vector<G4String>> parti
380 G4double sum = 0.0; //will contain a su
381 G4double kaonicFSprob=0.05; //probabili
382
383 ParticleList starlist;
384 ThreeVector mommy; //momentum to be ass
385
386 G4double rdm = Random::shoot();
387 ThreeVector annihilationPosition(0.,0.,0
388 if (rdm < (1.-kaonicFSprob)) { // pioni
389 INCL_DEBUG("pionic pp final state chos
390 sum = read_file(dataPathppbar, probabi
391 rdm = (rdm/(1.-kaonicFSprob))*sum; //
392 //now get the line number in the file
393 G4int n = findStringNumber(rdm, std::m
394 if ( n < 0 ) return theEventInfo;
395 for (G4int j = 0; j < static_cast<G4in
396 if (particle_types[n][j] == "pi0") {
397 Particle *p = new Particle(PiZero,
398 starlist.push_back(p);
399 } else if (particle_types[n][j] == "
400 Particle *p = new Particle(PiMinus
401 starlist.push_back(p);
402 } else if (particle_types[n][j] == "
403 Particle *p = new Particle(PiPlus,
404 starlist.push_back(p);
405 } else if (particle_types[n][j] == "
406 Particle *p = new Particle(Omega,
407 starlist.push_back(p);
408 } else if (particle_types[n][j] == "
409 Particle *p = new Particle(Eta, mo
410 starlist.push_back(p);
411 } else if (particle_types[n][j] == "
412 Particle *p = new Particle(PiMinus
413 starlist.push_back(p);
414 Particle *pp = new Particle(PiZero
415 starlist.push_back(pp);
416 } else if (particle_types[n][j] == "
417 Particle *p = new Particle(PiPlus,
418 starlist.push_back(p);
419 Particle *pp = new Particle(PiZero
420 starlist.push_back(pp);
421 } else if (particle_types[n][j] == "
422 Particle *p = new Particle(PiMinus
423 starlist.push_back(p);
424 Particle *pp = new Particle(PiPlus
425 starlist.push_back(pp);
426 } else {
427 INCL_ERROR("Some non-existing FS p
428 for (G4int jj = 0; jj < static_cas
429 #ifdef INCLXX_IN_GEANT4_MODE
430 G4cout << "gotcha! " << particle
431 #else
432 std::cout << "gotcha! " << parti
433 #endif
434 }
435 #ifdef INCLXX_IN_GEANT4_MODE
436 G4cout << "Some non-existing FS pa
437 #else
438 std::cout << "Some non-existing FS
439 #endif
440 }
441 }
442 } else {
443 INCL_DEBUG("kaonic pp final state chos
444 sum = read_file(dataPathppbark, probab
445 rdm = ((1.-rdm)/kaonicFSprob)*sum; //
446 //now get the line number in the file
447 G4int n = findStringNumber(rdm, probab
448 if ( n < 0 ) return theEventInfo;
449 for (G4int j = 0; j < static_cast<G4in
450 if (particle_types[n][j] == "pi0") {
451 Particle *p = new Particle(PiZero,
452 starlist.push_back(p);
453 } else if (particle_types[n][j] == "
454 Particle *p = new Particle(PiMinus
455 starlist.push_back(p);
456 } else if (particle_types[n][j] == "
457 Particle *p = new Particle(PiPlus,
458 starlist.push_back(p);
459 } else if (particle_types[n][j] == "
460 Particle *p = new Particle(Omega,
461 starlist.push_back(p);
462 } else if (particle_types[n][j] == "
463 Particle *p = new Particle(Eta, mo
464 starlist.push_back(p);
465 } else if (particle_types[n][j] == "
466 Particle *p = new Particle(KMinus,
467 starlist.push_back(p);
468 } else if (particle_types[n][j] == "
469 Particle *p = new Particle(KPlus,
470 starlist.push_back(p);
471 } else if (particle_types[n][j] == "
472 Particle *p = new Particle(KZero,
473 starlist.push_back(p);
474 } else if (particle_types[n][j] == "
475 Particle *p = new Particle(KZeroBa
476 starlist.push_back(p);
477 } else {
478 INCL_ERROR("Some non-existing FS p
479 for (G4int jj = 0; jj < static_cas
480 #ifdef INCLXX_IN_GEANT4_MODE
481 G4cout << "gotcha! " << particle
482 #else
483 std::cout << "gotcha! " << parti
484 #endif
485 }
486 #ifdef INCLXX_IN_GEANT4_MODE
487 G4cout << "Some non-existing FS pa
488 #else
489 std::cout << "Some non-existing FS
490 #endif
491 }
492 }
493 }
494
495 //compute energies of mesons with a phas
496 G4double energyOfMesonStar=ParticleTable
497 if (starlist.size() < 2) {
498 INCL_ERROR("should never happen, at le
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*energyO
505 G4double mom1 = std::sqrt(s/4. - (std:
506 ThreeVector momentello = Random::normV
507 (*first)->setMomentum(momentello);
508 (*first)->adjustEnergyFromMomentum();
509 (*last)->setMomentum(-momentello);
510 (*last)->adjustEnergyFromMomentum();
511 } else {
512 PhaseSpaceGenerator::generate(energyOf
513 }
514
515 if (targetA==1) postCascade_pbarH1(starl
516 else postCascade_pbarH2(starl
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(projec
529
530 if(!targetInitSuccess) {
531 INCL_WARN("Target initialisation failed
532 theEventInfo.transparent=true;
533 return theEventInfo;
534 }
535
536 cascadeAction->beforeCascadeAction(propaga
537
538 const G4bool canRunCascade = preCascade(pr
539 if(canRunCascade) {
540 cascade();
541 postCascade(projectileSpecies, kineticEn
542 cascadeAction->afterCascadeAction(nucleu
543 }
544 updateGlobalInfo();
545 return theEventInfo;
546 }
547
548 G4bool INCL::preCascade(ParticleSpecies cons
549 // Reset theEventInfo
550 theEventInfo.reset();
551
552 EventInfo::eventNumber++;
553
554 // Fill in the event information
555 theEventInfo.projectileType = projectileSp
556 theEventInfo.Ap = (Short_t)projectileSpeci
557 theEventInfo.Zp = (Short_t)projectileSpeci
558 theEventInfo.Sp = (Short_t)projectileSpeci
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
565 theEventInfo.Zt = (Short_t)nucleus->getZ
566 }
567 else if(nucleus->getAnnihilationType()==NT
568 theEventInfo.annihilationN = true;
569 theEventInfo.At = (Short_t)nucleus->getA
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 PbarAt
581 if(projectileSpecies.theType == antiProt
582 INCL_DEBUG("at rest annihilation" << '
583 //theEventInfo.transparent = false;
584 } else {
585 theEventInfo.transparent = true;
586 return false;
587 }
588 }
589
590
591 // Randomly draw an impact parameter or us
592 // Config option
593 G4double impactParameter, phi;
594 if(fixedImpactParameter<0.) {
595 impactParameter = maxImpactParameter * s
596 phi = Random::shoot() * Math::twoPi;
597 } else {
598 impactParameter = fixedImpactParameter;
599 phi = 0.;
600 }
601 INCL_DEBUG("Selected impact parameter: " <
602
603 // Fill in the event information
604 theEventInfo.impactParameter = impactParam
605
606 const G4double effectiveImpactParameter =
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 = ef
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 = 10000
625 do {
626 // Run book keeping actions that should
627 cascadeAction->beforePropagationAction(p
628
629 // Get the avatar with the smallest time
630 // to that point in time.
631 IAvatar *avatar = propagationModel->prop
632
633 finalState->reset();
634
635 // Run book keeping actions that should
636 cascadeAction->afterPropagationAction(pr
637
638 if(avatar == 0) break; // No more avatar
639
640 // Run book keeping actions that should
641 cascadeAction->beforeAvatarAction(avatar
642
643 // Channel is responsible for calculatin
644 // selected avatar. There are different
645 // class IChannel is, again, an abstract
646 // the externally observable behavior of
647 // channels.
648 // The handling of the channel is transp
649 // Final state tells what changed...
650 avatar->fillFinalState(finalState);
651 // Run book keeping actions that should
652 cascadeAction->afterAvatarAction(avatar,
653
654 // So now we must give this information
655 nucleus->applyFinalState(finalState);
656 // and now we are ready to process the n
657
658 delete avatar;
659
660 ++loopCounter;
661 } while(continueCascade() && loopCounter<m
662
663 delete finalState;
664 }
665
666 void INCL::postCascade(const ParticleSpecies
667 // Fill in the event information
668 theEventInfo.stoppingTime = propagationMod
669
670 // The event bias
671 theEventInfo.eventBias = (Double_t) Partic
672
673 // Forced CN?
674 if(!(projectileSpecies.theType==antiProton
675 if(nucleus->getTryCompoundNucleus()) {
676 INCL_DEBUG("Trying compound nucleus" <
677 makeCompoundNucleus();
678 theEventInfo.transparent = forceTransp
679 // Global checks of conservation laws
680 #ifndef INCLXX_IN_GEANT4_MODE
681 if(!theEventInfo.transparent) globalCons
682 #endif
683 return;
684 }
685 }
686
687 if(!(projectileSpecies.theType==antiProton
688 theEventInfo.transparent = forceTranspar
689 }
690
691 if(theEventInfo.transparent) {
692 ProjectileRemnant * const projectileRemn
693 if(projectileRemnant) {
694 // Clear the incoming list (particles
695 nucleus->getStore()->clearIncoming();
696 } else {
697 // Delete particles in the incoming li
698 nucleus->getStore()->deleteIncoming();
699 }
700 } else {
701
702 // Check if the nucleus contains strange
703 theEventInfo.sigmasInside = nucleus->con
704 theEventInfo.antikaonsInside = nucleus->
705 theEventInfo.lambdasInside = nucleus->co
706 theEventInfo.kaonsInside = nucleus->cont
707
708 // Capture antiKaons and Sigmas and prod
709 theEventInfo.absorbedStrangeParticle = n
710
711 // Emit strange particles still inside t
712 nucleus->emitInsideStrangeParticles();
713 theEventInfo.emitKaon = nucleus->emitIns
714
715 #ifdef INCLXX_IN_GEANT4_MODE
716 theEventInfo.emitLambda = nucleus->emitI
717 #endif // INCLXX_IN_GEANT4_MODE
718
719 // Check if the nucleus contains deltas
720 theEventInfo.deltasInside = nucleus->con
721
722 // Take care of any remaining deltas
723 theEventInfo.forcedDeltasOutside = nucle
724 theEventInfo.forcedDeltasInside = nucleu
725
726 // Take care of any remaining etas, omeg
727 G4double timeThreshold=theConfig->getDec
728 theEventInfo.forcedPionResonancesOutside
729 nucleus->decayOutgoingSigmaZero(timeThre
730 nucleus->decayOutgoingNeutralKaon();
731
732 // Apply Coulomb distortion, if appropri
733 // Note that this will apply Coulomb dis
734 // unphysical remnants (see decayInsideD
735 // what INCL4.6 does, but these events a
736 // whatever we do doesn't (shouldn't!) m
737 CoulombDistortion::distortOut(nucleus->g
738
739 // If the normal cascade predicted compl
740 // masses to compute the excitation ener
741 if(nucleus->getStore()->getOutgoingParti
742 && (!nucleus->getProjectileRemnant()
743 || nucleus->getProjectileRemnant(
744
745 INCL_DEBUG("Cascade resulted in comple
746
747 nucleus->useFusionKinematics();
748
749 if(nucleus->getExcitationEnergy()<0.)
750 // Complete fusion is energetically
751 INCL_WARN("Complete-fusion kinematic
752 theEventInfo.transparent = true;
753 return;
754 }
755
756 } else { // Normal cascade here
757
758 // Set the excitation energy
759 nucleus->setExcitationEnergy(nucleus->
760
761 // Make a projectile pre-fragment out
762 // spectators
763 theEventInfo.nUnmergedSpectators = mak
764
765 // Compute recoil momentum, energy and
766 if(nucleus->getA()==1 && minRemnantSiz
767 INCL_ERROR("Computing one-nucleon re
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
777 // outgoing particles.
778 if(nucleus->hasRemnant()) rescaleOutgo
779
780 }
781
782 // Cluster decay
783 theEventInfo.clusterDecay = nucleus->dec
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 collis
798 // compound nucleus.
799 //
800 // Yes, even nucleon-nucleus collisions ca
801 // below the Fermi level. Take e.g. 1-MeV
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 i
816 // should actually take it into account!
817 ThreeVector theCNMomentum = nucleus->getIn
818 ThreeVector theCNSpin = nucleus->getIncomi
819 const G4double theTargetMass = ParticleTab
820 G4int theCNA=theEventInfo.At, theCNZ=theEv
821 Cluster * const theProjectileRemnant = nuc
822 G4double theCNEnergy = theTargetMass + the
823
824 // Loop over the potential participants
825 ParticleList const &initialProjectileCompo
826 std::vector<Particle *> shuffledComponents
827 // Shuffle the list of potential participa
828 std::shuffle(shuffledComponents.begin(), s
829
830 G4bool success = true;
831 G4bool atLeastOneNucleonEntering = false;
832 for(std::vector<Particle*>::const_iterator
833 // Skip particles that miss the interact
834 Intersection intersectionInteractionDist
835 (*p)->getPosition(),
836 (*p)->getPropagationVelocity(),
837 maxInteractionDistance));
838 if(!intersectionInteractionDistance.exis
839 continue;
840
841 // Build an entry avatar for this nucleo
842 atLeastOneNucleonEntering = true;
843 ParticleEntryAvatar *theAvatar = new Par
844 nucleus->getStore()->addParticleEntryAva
845 FinalState *fs = theAvatar->getFinalStat
846 nucleus->applyFinalState(fs);
847 FinalStateValidity validity = fs->getVal
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 force
868 forceTransparent = true;
869 return;
870 }
871
872 // assert(theCNA==nucleus->getA());
873 // assert(theCNA<=theEventInfo.At+theEventInfo
874 // assert(theCNZ<=theEventInfo.Zt+theEventInfo
875 // assert(theCNS>=theEventInfo.St+theEventInfo
876
877 // Update the kinematics of the CN
878 theCNEnergy -= theProjectileRemnant->getEn
879 theCNMomentum -= theProjectileRemnant->get
880
881 // Deal with the projectile remnant
882 nucleus->finalizeProjectileRemnant(propaga
883
884 // Subtract the angular momentum of the pr
885 // assert(nucleus->getStore()->getOutgoingPart
886 theCNSpin -= theProjectileRemnant->getAngu
887
888 // Compute the excitation energy of the CN
889 const G4double theCNMass = ParticleTable::
890 const G4double theCNInvariantMassSquared =
891 if(theCNInvariantMassSquared<0.) {
892 // Negative invariant mass squared, retu
893 forceTransparent = true;
894 return;
895 }
896 const G4double theCNExcitationEnergy = std
897 if(theCNExcitationEnergy<0.) {
898 // Negative excitation energy, return a
899 INCL_DEBUG("CN excitation energy is nega
900 << " theCNA = " << theCNA << '\n'
901 << " theCNZ = " << theCNZ << '\n'
902 << " theCNS = " << theCNS << '\n'
903 << " theCNEnergy = " << theCNEner
904 << " theCNMomentum = (" << theCNM
905 << " theCNExcitationEnergy = " <<
906 << " theCNSpin = (" << theCNSpin.
907 );
908 forceTransparent = true;
909 return;
910 } else {
911 // Positive excitation energy, can make
912 INCL_DEBUG("CN excitation energy is posi
913 << " theCNA = " << theCNA << '\n'
914 << " theCNZ = " << theCNZ << '\n'
915 << " theCNS = " << theCNS << '\n'
916 << " theCNEnergy = " << theCNEner
917 << " theCNMomentum = (" << theCNM
918 << " theCNExcitationEnergy = " <<
919 << " theCNSpin = (" << theCNSpin.
920 );
921 nucleus->setA(theCNA);
922 nucleus->setZ(theCNZ);
923 nucleus->setS(theCNS);
924 nucleus->setMomentum(theCNMomentum);
925 nucleus->setEnergy(theCNEnergy);
926 nucleus->setExcitationEnergy(theCNExcita
927 nucleus->setMass(theCNMass+theCNExcitati
928 nucleus->setSpin(theCNSpin); // neglects
929
930 // Take care of any remaining deltas
931 theEventInfo.forcedDeltasOutside = nucle
932
933 // Take care of any remaining etas and/o
934 G4double timeThreshold=theConfig->getDec
935 theEventInfo.forcedPionResonancesOutside
936
937 // Take care of any remaining Kaons
938 theEventInfo.emitKaon = nucleus->emitIns
939
940 // Cluster decay
941 theEventInfo.clusterDecay = nucleus->dec
942
943 // Fill the EventInfo structure
944 nucleus->fillEventInfo(&theEventInfo);
945 }
946 }
947
948 void INCL::rescaleOutgoingForRecoil() {
949 RecoilCMFunctor theRecoilFunctor(nucleus,
950
951 // Apply the root-finding algorithm
952 const RootFinder::Solution theSolution = R
953 if(theSolution.success) {
954 theRecoilFunctor(theSolution.x); // Appl
955 } else {
956 INCL_WARN("Couldn't accommodate remnant
957 }
958 }
959
960 #ifndef INCLXX_IN_GEANT4_MODE
961 void INCL::globalConservationChecks(G4bool a
962 Nucleus::ConservationBalance theBalance =
963
964 // Global conservation checks
965 const G4double pLongBalance = theBalance.m
966 const G4double pTransBalance = theBalance.
967 if(theBalance.Z != 0) {
968 INCL_ERROR("Violation of charge conserva
969 }
970 if(theBalance.A != 0) {
971 INCL_ERROR("Violation of baryon-number c
972 }
973 if(theBalance.S != 0) {
974 INCL_ERROR("Violation of strange-number
975 }
976 G4double EThreshold, pLongThreshold, pTran
977 if(afterRecoil) {
978 // Less stringent checks after accommoda
979 EThreshold = 10.; // MeV
980 pLongThreshold = 1.; // MeV/c
981 pTransThreshold = 1.; // MeV/c
982 } else {
983 // More stringent checks before accommod
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 conservat
990 }
991 if(std::abs(pLongBalance)>pLongThreshold)
992 INCL_WARN("Violation of longitudinal mom
993 }
994 if(std::abs(pTransBalance)>pTransThreshold
995 INCL_WARN("Violation of transverse momen
996 }
997
998 // Feed the EventInfo variables
999 theEventInfo.EBalance = theBalance.energy;
1000 theEventInfo.pLongBalance = pLongBalance;
1001 theEventInfo.pTransBalance = pTransBalanc
1002 }
1003 #endif
1004
1005 G4bool INCL::continueCascade() {
1006 // Stop if we have passed the stopping ti
1007 if(propagationModel->getCurrentTime() > p
1008 INCL_DEBUG("Cascade time (" << propagat
1009 << ") exceeded stopping time (" <<
1010 << "), stopping cascade" << '\n');
1011 return false;
1012 }
1013 // Stop if there are no participants and
1014 if(nucleus->getStore()->getBook().getCasc
1015 nucleus->getStore()->getIncomingParti
1016 INCL_DEBUG("No participants in the nucl
1017 return false;
1018 }
1019 // Stop if the remnant is smaller than mi
1020 if(nucleus->getA() <= minRemnantSize) {
1021 INCL_DEBUG("Remnant size (" << nucleus-
1022 << ") smaller than or equal to mini
1023 << "), stopping cascade" << '\n');
1024 return false;
1025 }
1026 // Stop if we have to try and make a comp
1027 // force a transparent
1028 if(nucleus->getTryCompoundNucleus()) {
1029 INCL_DEBUG("Trying to make a compound n
1030 return false;
1031 }
1032
1033 return true;
1034 }
1035
1036 void INCL::finalizeGlobalInfo(Random::SeedV
1037 const G4double normalisationFactor = theG
1038 ((G4double) theGlobalInfo.nShots);
1039 theGlobalInfo.nucleonAbsorptionCrossSecti
1040 ((G4double) theGlobalInfo.nNucleonAbsor
1041 theGlobalInfo.pionAbsorptionCrossSection
1042 ((G4double) theGlobalInfo.nPionAbsorpti
1043 theGlobalInfo.reactionCrossSection = norm
1044 ((G4double) (theGlobalInfo.nShots - the
1045 theGlobalInfo.errorReactionCrossSection =
1046 std::sqrt((G4double) (theGlobalInfo.nSh
1047 theGlobalInfo.forcedCNCrossSection = norm
1048 ((G4double) theGlobalInfo.nForcedCompou
1049 theGlobalInfo.errorForcedCNCrossSection =
1050 std::sqrt((G4double) (theGlobalInfo.nFo
1051 theGlobalInfo.completeFusionCrossSection
1052 ((G4double) theGlobalInfo.nCompleteFusi
1053 theGlobalInfo.errorCompleteFusionCrossSec
1054 std::sqrt((G4double) (theGlobalInfo.nCo
1055 theGlobalInfo.energyViolationInteractionC
1056 ((G4double) theGlobalInfo.nEnergyViolat
1057
1058 theGlobalInfo.initialRandomSeeds.assign(i
1059
1060 Random::SeedVector theSeeds = Random::get
1061 theGlobalInfo.finalRandomSeeds.assign(the
1062 }
1063
1064 G4int INCL::makeProjectileRemnant() {
1065 // Do nothing if this is not a nucleus-nu
1066 if(!nucleus->getProjectileRemnant())
1067 return 0;
1068
1069 // Get the spectators (geometrical+dynami
1070 ParticleList geomSpectators(nucleus->getP
1071 ParticleList dynSpectators(nucleus->getSt
1072
1073 G4int nUnmergedSpectators = 0;
1074
1075 // If there are no spectators, do nothing
1076 if(dynSpectators.empty() && geomSpectator
1077 return 0;
1078 } else if(dynSpectators.size()==1 && geom
1079 // No geometrical spectators, one dynam
1080 // Just put it back in the outgoing lis
1081 nucleus->getStore()->addToOutgoing(dynS
1082 } else {
1083 // Make a cluster out of the geometrica
1084 ProjectileRemnant *theProjectileRemnant
1085
1086 // Add the dynamical spectators to the
1087 ParticleList rejected = theProjectileRe
1088 // Put back the rejected spectators int
1089 nUnmergedSpectators = (G4int)rejected.s
1090 nucleus->getStore()->addToOutgoing(reje
1091
1092 // Deal with the projectile remnant
1093 nucleus->finalizeProjectileRemnant(prop
1094
1095 }
1096
1097 return nUnmergedSpectators;
1098 }
1099
1100 void INCL::initMaxInteractionDistance(Parti
1101 if(projectileSpecies.theType != Composite
1102 maxInteractionDistance = 0.;
1103 return;
1104 }
1105
1106 const G4double r0 = std::max(ParticleTabl
1107 ParticleTable:
1108
1109 const G4double theNNDistance = CrossSecti
1110 maxInteractionDistance = r0 + theNNDistan
1111 INCL_DEBUG("Initialised interaction dista
1112 << " theNNDistance = " << theNND
1113 << " maxInteractionDistance = "
1114 }
1115
1116 void INCL::initUniverseRadius(ParticleSpeci
1117 G4double rMax = 0.0;
1118 if(A==0) {
1119 IsotopicDistribution const &anIsotopicD
1120 ParticleTable::getNaturalIsotopicDist
1121 IsotopeVector theIsotopes = anIsotopicD
1122 for(IsotopeIter i=theIsotopes.begin(),
1123 const G4double pMaximumRadius = Parti
1124 const G4double nMaximumRadius = Parti
1125 const G4double maximumRadius = std::m
1126 rMax = std::max(maximumRadius, rMax);
1127 }
1128 } else {
1129 const G4double pMaximumRadius = Particl
1130 const G4double nMaximumRadius = Particl
1131 const G4double maximumRadius = std::max
1132 rMax = std::max(maximumRadius, rMax);
1133 }
1134 if(p.theType==Composite || p.theType==Pro
1135 const G4double interactionDistanceNN =
1136 maxUniverseRadius = rMax + interactionD
1137 } else if(p.theType==PiPlus
1138 || p.theType==PiZero
1139 || p.theType==PiMinus) {
1140 const G4double interactionDistancePiN =
1141 maxUniverseRadius = rMax + interactionD
1142 } else if(p.theType==KPlus
1143 || p.theType==KZero) {
1144 const G4double interactionDistanceKN =
1145 maxUniverseRadius = rMax + interactionD
1146 } else if(p.theType==KZeroBar
1147 || p.theType==KMinus) {
1148 const G4double interactionDistanceKbarN
1149 maxUniverseRadius = rMax + interactionD
1150 } else if(p.theType==Lambda
1151 ||p.theType==SigmaPlus
1152 || p.theType==SigmaZero
1153 || p.theType==SigmaMinus) {
1154 const G4double interactionDistanceYN =
1155 maxUniverseRadius = rMax + interactionD
1156 }
1157 else if(p.theType==antiProton) {
1158 maxUniverseRadius = rMax;
1159 }
1160 INCL_DEBUG("Initialised universe radius:
1161 }
1162
1163
1164 G4double INCL::initUniverseRadiusForAntipro
1165 G4double rMax = 0.0;
1166 if(A==0) {
1167 IsotopicDistribution const &anIsotopicD
1168 ParticleTable::getNaturalIsotopicDist
1169 IsotopeVector theIsotopes = anIsotopicD
1170 for(IsotopeIter i=theIsotopes.begin(),
1171 const G4double pMaximumRadius = Parti
1172 const G4double nMaximumRadius = Parti
1173 const G4double maximumRadius = std::m
1174 rMax = std::max(maximumRadius, rMax);
1175 }
1176 } else {
1177 const G4double pMaximumRadius = Particl
1178 const G4double nMaximumRadius = Particl
1179 const G4double maximumRadius = std::max
1180 rMax = std::max(maximumRadius, rMax);
1181 }
1182 return rMax;
1183 }
1184
1185
1186 void INCL::updateGlobalInfo() {
1187 // Increment the global counter for the n
1188 theGlobalInfo.nShots++;
1189
1190 if(theEventInfo.transparent) {
1191 // Increment the global counter for the
1192 theGlobalInfo.nTransparents++;
1193 // Increment the global counter for the
1194 if(forceTransparent)
1195 theGlobalInfo.nForcedTransparents++;
1196 return;
1197 }
1198
1199 // Check if we have an absorption:
1200 if(theEventInfo.nucleonAbsorption) theGlo
1201 if(theEventInfo.pionAbsorption) theGlobal
1202
1203 // Count complete-fusion events
1204 if(theEventInfo.nCascadeParticles==0) the
1205
1206 if(nucleus->getTryCompoundNucleus())
1207 theGlobalInfo.nForcedCompoundNucleus++;
1208
1209 // Counters for the number of violations
1210 // collisions
1211 theGlobalInfo.nEnergyViolationInteraction
1212 }
1213
1214 G4double INCL::read_file(std::string filena
1215 std::vector<std::v
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(ty
1232 }
1233 } else {
1234 #ifdef INCLXX_IN_GEANT4_MODE
1235 G4cout << "ERROR no fread_file " << fil
1236 #else
1237 std::cout << "ERROR no fread_file " <<
1238 #endif
1239 }
1240 return sum_probs;
1241 }
1242
1243
1244 G4int INCL::findStringNumber(G4double rdm,
1245 G4int stringNumber = -1;
1246 G4double smallestsum = 0.0;
1247 G4double biggestsum = yields[0];
1248 //G4cout << "initial input " << rdm << G4
1249 for (G4int i = 0; i < static_cast<G4int>(
1250 if (rdm >= smallestsum && rdm <= bigges
1251 //G4cout << smallestsum << " and " <<
1252 stringNumber = i+1;
1253 }
1254 smallestsum += yields[i];
1255 biggestsum += yields[i+1];
1256 }
1257 if(stringNumber==-1) stringNumber = stati
1258 if(stringNumber==-1){
1259 INCL_ERROR("ERROR in findStringNumber (
1260 #ifdef INCLXX_IN_GEANT4_MODE
1261 G4cout << "ERROR in findStringNumber" <
1262 #else
1263 std::cout << "ERROR in findStringNumber
1264 #endif
1265 }
1266 return stringNumber;
1267 }
1268
1269
1270 void INCL::preCascade_pbarH1(ParticleSpecie
1271 // Reset theEventInfo
1272 theEventInfo.reset();
1273
1274 EventInfo::eventNumber++;
1275
1276 // Fill in the event information
1277 theEventInfo.projectileType = projectileS
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
1289 theEventInfo.nParticles = 0;
1290
1291 // Reset the remnant counter
1292 theEventInfo.nRemnants = 0;
1293 theEventInfo.history.clear();
1294
1295 for(ParticleIter i=outgoingParticles.begi
1296 theEventInfo.A[theEventInfo.nParticles]
1297 theEventInfo.Z[theEventInfo.nParticles]
1298 theEventInfo.S[theEventInfo.nParticles]
1299 theEventInfo.EKin[theEventInfo.nParticl
1300 ThreeVector mom = (*i)->getMomentum();
1301 theEventInfo.px[theEventInfo.nParticles
1302 theEventInfo.py[theEventInfo.nParticles
1303 theEventInfo.pz[theEventInfo.nParticles
1304 theEventInfo.theta[theEventInfo.nPartic
1305 theEventInfo.phi[theEventInfo.nParticle
1306 theEventInfo.origin[theEventInfo.nParti
1307 #ifdef INCLXX_IN_GEANT4_MODE
1308 theEventInfo.parentResonancePDGCode[the
1309 theEventInfo.parentResonanceID[theEvent
1310 #endif
1311 theEventInfo.history.push_back("");
1312 ParticleSpecies pt((*i)->getType());
1313 theEventInfo.PDGCode[theEventInfo.nPart
1314 theEventInfo.nParticles++;
1315 }
1316 theEventInfo.nCascadeParticles = theEvent
1317 }
1318
1319
1320 void INCL::preCascade_pbarH2(ParticleSpecie
1321 // Reset theEventInfo
1322 theEventInfo.reset();
1323
1324 EventInfo::eventNumber++;
1325
1326 // Fill in the event information
1327 theEventInfo.projectileType = projectileS
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
1339 theEventInfo.nParticles = 0;
1340
1341 // Reset the remnant counter
1342 theEventInfo.nRemnants = 0;
1343 theEventInfo.history.clear();
1344
1345 for(ParticleIter i=outgoingParticles.begi
1346 theEventInfo.A[theEventInfo.nParticles]
1347 theEventInfo.Z[theEventInfo.nParticles]
1348 theEventInfo.S[theEventInfo.nParticles]
1349 theEventInfo.EKin[theEventInfo.nParticl
1350 ThreeVector mom = (*i)->getMomentum();
1351 theEventInfo.px[theEventInfo.nParticles
1352 theEventInfo.py[theEventInfo.nParticles
1353 theEventInfo.pz[theEventInfo.nParticles
1354 theEventInfo.theta[theEventInfo.nPartic
1355 theEventInfo.phi[theEventInfo.nParticle
1356 theEventInfo.origin[theEventInfo.nParti
1357 #ifdef INCLXX_IN_GEANT4_MODE
1358 theEventInfo.parentResonancePDGCode[the
1359 theEventInfo.parentResonanceID[theEvent
1360 #endif
1361 theEventInfo.history.push_back("");
1362 ParticleSpecies pt((*i)->getType());
1363 theEventInfo.PDGCode[theEventInfo.nPart
1364 theEventInfo.nParticles++;
1365 }
1366
1367 for(ParticleIter i=H2Particles.begin(), e
1368 theEventInfo.A[theEventInfo.nParticles]
1369 theEventInfo.Z[theEventInfo.nParticles]
1370 theEventInfo.S[theEventInfo.nParticles]
1371 theEventInfo.EKin[theEventInfo.nParticl
1372 ThreeVector mom = (*i)->getMomentum();
1373 theEventInfo.px[theEventInfo.nParticles
1374 theEventInfo.py[theEventInfo.nParticles
1375 theEventInfo.pz[theEventInfo.nParticles
1376 theEventInfo.theta[theEventInfo.nPartic
1377 theEventInfo.phi[theEventInfo.nParticle
1378 theEventInfo.origin[theEventInfo.nParti
1379 #ifdef INCLXX_IN_GEANT4_MODE
1380 theEventInfo.parentResonancePDGCode[the
1381 theEventInfo.parentResonanceID[theEvent
1382 #endif
1383 theEventInfo.history.push_back("");
1384 ParticleSpecies pt((*i)->getType());
1385 theEventInfo.PDGCode[theEventInfo.nPart
1386 theEventInfo.nParticles++;
1387 }
1388 theEventInfo.nCascadeParticles = theEvent
1389 }
1390
1391 }
1392