| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA
29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA
31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege
32 // Davide Mancusi, CEA-Saclay, France << 32 //
>> 33 // INCL++ revision: v5.0_rc3
33 // 34 //
34 #define INCLXX_IN_GEANT4_MODE 1 35 #define INCLXX_IN_GEANT4_MODE 1
35 36
36 #include "globals.hh" 37 #include "globals.hh"
37 38
38 /** \file G4INCLCascade.cc <<
39 * <<
40 * INCL Cascade <<
41 */ <<
42 #include "G4INCLCascade.hh" 39 #include "G4INCLCascade.hh"
43 #include "G4INCLRandom.hh" 40 #include "G4INCLRandom.hh"
>> 41 #include "G4INCLRanecu.hh"
>> 42 #include "G4INCLGeant4Random.hh"
44 #include "G4INCLStandardPropagationModel.hh" 43 #include "G4INCLStandardPropagationModel.hh"
45 #include "G4INCLParticleTable.hh" 44 #include "G4INCLParticleTable.hh"
46 #include "G4INCLParticle.hh" <<
47 #include "G4INCLNuclearMassTable.hh" <<
48 #include "G4INCLGlobalInfo.hh" 45 #include "G4INCLGlobalInfo.hh"
49 #include "G4INCLNucleus.hh" <<
50 46
51 #include "G4INCLPauliBlocking.hh" 47 #include "G4INCLPauliBlocking.hh"
52 << 48 #include "G4INCLIPauli.hh"
53 #include "G4INCLCrossSections.hh" << 49 #include "G4INCLPauliStrict.hh"
54 << 50 #include "G4INCLPauliStandard.hh"
55 #include "G4INCLPhaseSpaceGenerator.hh" << 51 #include "G4INCLPauliStrictStandard.hh"
>> 52 #include "G4INCLPauliGlobal.hh"
>> 53 #include "G4INCLCDPP.hh"
56 54
57 #include "G4INCLLogger.hh" 55 #include "G4INCLLogger.hh"
58 #include "G4INCLGlobals.hh" 56 #include "G4INCLGlobals.hh"
59 #include "G4INCLNuclearDensityFactory.hh" 57 #include "G4INCLNuclearDensityFactory.hh"
60 58
61 #include "G4INCLINuclearPotential.hh" <<
62 <<
63 #include "G4INCLCoulombDistortion.hh" 59 #include "G4INCLCoulombDistortion.hh"
>> 60 #include "G4INCLICoulomb.hh"
>> 61 #include "G4INCLCoulombNone.hh"
>> 62 #include "G4INCLCoulombNonRelativistic.hh"
64 63
65 #include "G4INCLClustering.hh" 64 #include "G4INCLClustering.hh"
>> 65 #include "G4INCLClusteringModelIntercomparison.hh"
>> 66 #include "G4INCLClusteringModelNone.hh"
66 67
67 #include "G4INCLIntersection.hh" << 68 #include <cstring>
68 <<
69 #include "G4INCLBinaryCollisionAvatar.hh" <<
70 <<
71 #include "G4INCLCascadeAction.hh" <<
72 #include "G4INCLAvatarDumpAction.hh" <<
73 <<
74 #include <cstring> <<
75 #include <cstdlib> 69 #include <cstdlib>
76 #include <numeric> <<
77 <<
78 #include "G4INCLPbarAtrestEntryChannel.hh" <<
79 70
80 namespace G4INCL { 71 namespace G4INCL {
81 << 72
82 INCL::INCL(Config const * const config) << 73 INCL::INCL(G4INCL::Config const * const config)
83 :propagationModel(0), theA(208), theZ(82), << 74 :propagationModel(0), theA(208), theZ(82), maxImpactParameter(0.),
84 targetInitSuccess(false), << 75 theConfig(config)
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 { 76 {
94 // Set the logger object. 77 // Set the logger object.
95 #ifdef INCLXX_IN_GEANT4_MODE << 78 G4INCL::Logger::setLoggerSlave(new G4INCL::LoggerSlave(theConfig->getLogFileName()));
96 Logger::initVerbosityLevelFromEnvvar(); << 79 G4INCL::Logger::setVerbosityLevel(theConfig->getVerbosity());
97 #else // INCLXX_IN_GEANT4_MODE <<
98 Logger::initialize(theConfig); <<
99 #endif // INCLXX_IN_GEANT4_MODE <<
100 80
101 // Set the random number generator algorit 81 // Set the random number generator algorithm. The system can support
102 // multiple different generator algorithms 82 // multiple different generator algorithms in a completely
103 // transparent way. 83 // transparent way.
104 Random::initialize(theConfig); << 84 #ifdef INCLXX_IN_GEANT4_MODE
105 << 85 G4INCL::Random::setGenerator(new G4INCL::Geant4RandomGenerator());
106 // Select the Pauli and CDPP blocking algo << 86 #else
107 Pauli::initialize(theConfig); << 87 G4INCL::Random::setGenerator(new G4INCL::Ranecu(theConfig->getRandomSeeds()));
>> 88 #endif // INCLXX_IN_GEANT4_MODE
108 89
109 // Set the cross-section set << 90 // Select the Pauli blocking algorithm:
110 CrossSections::initialize(theConfig); << 91 G4INCL::PauliType pauli = theConfig->getPauliType();
>> 92 if(pauli == G4INCL::StrictStatisticalPauli)
>> 93 G4INCL::Pauli::setBlocker(new G4INCL::PauliStrictStandard);
>> 94 else if(pauli == G4INCL::StatisticalPauli)
>> 95 G4INCL::Pauli::setBlocker(new G4INCL::PauliStandard);
>> 96 else if(pauli == G4INCL::StrictPauli)
>> 97 G4INCL::Pauli::setBlocker(new G4INCL::PauliStrict);
>> 98 else if(pauli == G4INCL::GlobalPauli)
>> 99 G4INCL::Pauli::setBlocker(new G4INCL::PauliGlobal);
>> 100 else if(pauli == G4INCL::NoPauli)
>> 101 G4INCL::Pauli::setBlocker(NULL);
111 102
112 // Set the phase-space generator << 103 if(theConfig->getCDPP())
113 PhaseSpaceGenerator::initialize(theConfig) << 104 G4INCL::Pauli::setCDPP(new G4INCL::CDPP);
>> 105 else
>> 106 G4INCL::Pauli::setCDPP(NULL);
114 107
115 // Select the Coulomb-distortion algorithm 108 // Select the Coulomb-distortion algorithm:
116 CoulombDistortion::initialize(theConfig); << 109 G4INCL::CoulombType coulombType = theConfig->getCoulombType();
>> 110 if(coulombType == G4INCL::NonRelativisticCoulomb)
>> 111 G4INCL::CoulombDistortion::setCoulomb(new G4INCL::CoulombNonRelativistic);
>> 112 else // if(coulombType == G4INCL::NoCoulomb)
>> 113 G4INCL::CoulombDistortion::setCoulomb(new G4INCL::CoulombNone);
117 114
118 // Select the clustering algorithm: 115 // Select the clustering algorithm:
119 Clustering::initialize(theConfig); << 116 G4INCL::ClusterAlgorithmType clusterAlgorithm = theConfig->getClusterAlgorithm();
>> 117 if(clusterAlgorithm == G4INCL::IntercomparisonClusterAlgorithm) {
>> 118 G4INCL::Clustering::setClusteringModel(new G4INCL::ClusteringModelIntercomparison);
>> 119 // Set the maximum mass for the clustering algorithm
>> 120 G4INCL::IClusteringModel::maxClusterAlgorithmMass = theConfig->getClusterMaxMass();
>> 121 }
>> 122 else // if(clusterAlgorithm == G4INCL::NoClusterAlgorithm)
>> 123 G4INCL::Clustering::setClusteringModel(new G4INCL::ClusteringModelNone);
120 124
121 // Initialize the INCL particle table: 125 // Initialize the INCL particle table:
122 ParticleTable::initialize(theConfig); << 126 G4INCL::ParticleTable::initialize();
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 127
130 // Propagation model is responsible for fi 128 // Propagation model is responsible for finding avatars and
131 // transporting the particles. In principl 129 // transporting the particles. In principle this step is "hidden"
132 // behind an abstract interface and the re << 130 // behind an abstract G4interface and the rest of the system does not
133 // care how the transportation and avatar 131 // care how the transportation and avatar finding is done. This
134 // should allow us to "easily" experiment 132 // should allow us to "easily" experiment with different avatar
135 // finding schemes and even to support thi 133 // finding schemes and even to support things like curved
136 // trajectories in the future. 134 // trajectories in the future.
137 propagationModel = new StandardPropagation << 135 propagationModel = new G4INCL::StandardPropagationModel(theConfig->getLocalEnergyBBType(),theConfig->getLocalEnergyPiType());
138 if(theConfig->getCascadeActionType() == Av << 136 eventAction = new EventAction();
139 cascadeAction = new AvatarDumpAction(); << 137 propagationAction = new PropagationAction();
140 else << 138 avatarAction = new AvatarAction();
141 cascadeAction = new CascadeAction(); << 139
142 cascadeAction->beforeRunAction(theConfig); << 140 std::strcpy(theGlobalInfo.cascadeModel, theConfig->getVersionID().c_str());
143 << 141 std::strcpy(theGlobalInfo.deexcitationModel, "none");
144 theGlobalInfo.cascadeModel = theConfig->ge << 142
145 theGlobalInfo.deexcitationModel = theConfi << 143 // Set the target
146 #ifdef INCL_ROOT_USE << 144 if(!theConfig->isNaturalTarget()) {
147 theGlobalInfo.rootSelection = theConfig->g << 145 setTarget(theConfig->getTargetA(), theConfig->getTargetZ());
148 #endif << 146 // Fill in the global information
>> 147 theGlobalInfo.At = theConfig->getTargetA();
>> 148 theGlobalInfo.Zt = theConfig->getTargetZ();
>> 149 } else {
>> 150 // TODO: support for natural targets
>> 151 FATAL("Fatal: natural targets are not supported yet." << std::endl);
>> 152 std::exit(EXIT_FAILURE);
>> 153 }
149 154
150 #ifndef INCLXX_IN_GEANT4_MODE 155 #ifndef INCLXX_IN_GEANT4_MODE
151 // Fill in the global information << 156 // Echo the input parameters to the log file
152 theGlobalInfo.At = theConfig->getTargetA() << 157 INFO(theConfig->echo() << std::endl);
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 158 #endif
>> 159 }
162 160
163 fixedImpactParameter = theConfig->getImpac << 161 INCL::INCL(IPropagationModel *aPropagationModel)
>> 162 :propagationModel(aPropagationModel), theA(208), theZ(82), maxImpactParameter(0.), theConfig(NULL)
>> 163 {
>> 164 // Set the random number generator algorithm. The system can support
>> 165 // multiple different generator algorithms in a completely
>> 166 // transparent way.
>> 167 G4INCL::Random::setGenerator(new G4INCL::Ranecu());
164 } 168 }
165 169
166 INCL::~INCL() { 170 INCL::~INCL() {
167 InteractionAvatar::deleteBackupParticles() << 171 G4INCL::Pauli::deleteBlockers();
168 #ifndef INCLXX_IN_GEANT4_MODE << 172 G4INCL::CoulombDistortion::deleteCoulomb();
169 NuclearMassTable::deleteTable(); << 173 G4INCL::Random::deleteGenerator();
170 #endif << 174 G4INCL::ParticleTable::deletePDS();
171 PhaseSpaceGenerator::deletePhaseSpaceGener << 175 G4INCL::Clustering::deleteClusteringModel();
172 CrossSections::deleteCrossSections(); << 176 G4INCL::Logger::deleteLoggerSlave();
173 Pauli::deleteBlockers(); << 177 delete avatarAction;
174 CoulombDistortion::deleteCoulomb(); << 178 delete propagationAction;
175 Random::deleteGenerator(); << 179 delete eventAction;
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; 180 delete propagationModel;
185 delete theConfig; <<
186 } 181 }
187 182
188 G4bool INCL::prepareReaction(const ParticleS << 183 void INCL::setTarget(G4int A, G4int Z) {
189 if(A < 0 || A > 300 || Z < 1 || Z > 200) { << 184 if(A > 0 && A < 300 && Z > 0 && Z < 200) {
190 INCL_ERROR("Unsupported target: A = " << << 185 theA = 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; 186 theZ = Z;
247 theS = S; << 187 } else {
248 if(theConfig->isNaturalTarget()) << 188 ERROR("Unsupported target: A = " << A << " Z = " << Z << std::endl);
249 theA = ParticleTable::drawRandomNatura << 189 ERROR("Target configuration rejected." << std::endl);
250 else <<
251 theA = A; <<
252 } 190 }
253 191
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 192 // Set the maximum impact parameter
265 maxImpactParameter = CoulombDistortion::ma << 193 // TODO: for natural target abundances, make this the largest impact
266 INCL_DEBUG("Maximum impact parameter initi << 194 // parameter for all the isotopes.
267 << 195 // TODO: reduce the maximum impact parameter for Coulomb-distorted
268 // For forced CN events << 196 // trajectories. Make this dependent on the configuration choice for
269 initMaxInteractionDistance(projectileSpeci << 197 // Coulomb distortion.
270 // Set the geometric cross sectiony section << 198 NuclearDensity const * const density = NuclearDensityFactory::createDensity(A,Z);
271 if(projectileSpecies.theType == antiProton << 199 maxImpactParameter = density->getMaximumRadius();
272 G4int currentA = A; << 200 delete density;
273 if(theConfig->isNaturalTarget()){ << 201
274 currentA = ParticleTable::drawRandomNa << 202 // Set the geometric cross section
275 } << 203 theGlobalInfo.geometricCrossSection =
276 G4double kineticEnergy2=kineticEnergy; << 204 Math::tenPi*std::pow(maxImpactParameter,2);
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 205
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 } 206 }
295 207
296 G4bool INCL::initializeTarget(const G4int A, << 208 G4bool INCL::initializeTarget(G4int A, G4int Z) {
297 delete nucleus; << 209 Nucleus *previousNucleus = propagationModel->getNucleus();
>> 210 delete previousNucleus;
>> 211
>> 212 Nucleus *aNucleus = new Nucleus(A, Z, theConfig);
>> 213 aNucleus->getStore()->getBook()->reset();
>> 214 aNucleus->initializeParticles();
298 215
299 if (theAType==PType || theAType==NType) { << 216 propagationModel->setNucleus(aNucleus);
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; 217 return true;
312 } 218 }
313 219
314 const EventInfo &INCL::processEvent( << 220 const EventInfo &INCL::processEvent(Particle *projectile) {
315 ParticleSpecies const &projectileSpecies << 221 initializeTarget(theA, theZ);
316 const G4double kineticEnergy, << 222 // Usage of the projectile API:
317 const G4int targetA, << 223
318 const G4int targetZ, << 224 // Test projectile:
319 const G4int targetS << 225 // G4INCL::ThreeVector position(0.0, 0.0, 0.0);
320 ) { << 226 // G4INCL::ThreeVector momentum(0.0, 0.0, 2000.0);
321 << 227 // G4double energy = std::sqrt(momentum.mag2() + G4INCL::ProtonMass * G4INCL::ProtonMass);
322 ParticleList starlistH2; << 228 // G4INCL::Particle *projectile = new G4INCL::Particle(G4INCL::Proton, energy,
323 << 229 // momentum, position);
324 if (projectileSpecies.theType==antiProton << 230
325 << 231 // composite
326 if (targetA==1) { << 232 // G4INCL::Nucleus *projectileNucleus = new G4INCL::Nucleus(6, 12);
327 preCascade_pbarH1(projectileSpecies, k << 233 //projectileNucleus->initializeParticles();
328 } else { << 234
329 preCascade_pbarH2(projectileSpecies, k << 235 // Create a nucleus of Z = 82 and A = 208
330 theEventInfo.annihilationP = false; << 236 // G4INCL::Nucleus *theNucleus = new G4INCL::Nucleus(6, 12);
331 theEventInfo.annihilationN = false; << 237
332 << 238 // Generate the initial distribution of particles
333 G4double SpOverSn = 1.331; //from exp << 239 // theNucleus->initializeParticles();
334 << 240
335 ThreeVector dummy(0.,0.,0.); << 241 // theNucleus->shootMe(projectile);
336 G4double rndm = Random::shoot()*(SpOve << 242 // theNucleus->shootMe(projectileNucleus);
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 243
522 // ReInitialize the bias vector << 244 // Manually set the stopping time of the simulation.
523 Particle::INCLBiasVector.clear(); << 245 // propagationModel->setStoppingTime(70.0);
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 246
536 cascadeAction->beforeCascadeAction(propaga << 247 // Assign the nucleus to the propagation model
>> 248 // propagationModel->setNucleus(theNucleus);
537 249
538 const G4bool canRunCascade = preCascade(pr << 250 // Shortcut poG4inter
539 if(canRunCascade) { << 251 Nucleus *nucleus = propagationModel->getNucleus();
540 cascade(); <<
541 postCascade(projectileSpecies, kineticEn <<
542 cascadeAction->afterCascadeAction(nucleu <<
543 } <<
544 updateGlobalInfo(); <<
545 return theEventInfo; <<
546 } <<
547 252
548 G4bool INCL::preCascade(ParticleSpecies cons <<
549 // Reset theEventInfo 253 // Reset theEventInfo
550 theEventInfo.reset(); 254 theEventInfo.reset();
551 << 255
552 EventInfo::eventNumber++; 256 EventInfo::eventNumber++;
553 257
554 // Fill in the event information << 258 // Increment the global counter for the number of shots
555 theEventInfo.projectileType = projectileSp << 259 theGlobalInfo.nShots++;
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 260
591 // Randomly draw an impact parameter or us << 261 // Fill in the global information
592 // Config option << 262 // TODO: should be moved to the input processing
593 G4double impactParameter, phi; << 263 theGlobalInfo.Ap = projectile->getA();
594 if(fixedImpactParameter<0.) { << 264 theGlobalInfo.Zp = projectile->getZ();
595 impactParameter = maxImpactParameter * s << 265 theGlobalInfo.Ep = projectile->getKineticEnergy();
596 phi = Random::shoot() * Math::twoPi; << 266
597 } else { << 267 // Fill in the event information
598 impactParameter = fixedImpactParameter; << 268 theEventInfo.projectileType = projectile->getType();
599 phi = 0.; << 269 theEventInfo.Ap = projectile->getA();
600 } << 270 theEventInfo.Zp = projectile->getZ();
601 INCL_DEBUG("Selected impact parameter: " < << 271 theEventInfo.Ep = projectile->getKineticEnergy();
>> 272 theEventInfo.At = nucleus->getA();
>> 273 theEventInfo.Zt = nucleus->getZ();
602 274
>> 275 // Randomly draw an impact parameter
>> 276 G4double impactParameter = maxImpactParameter * std::sqrt(Random::shoot());
603 // Fill in the event information 277 // Fill in the event information
604 theEventInfo.impactParameter = impactParam 278 theEventInfo.impactParameter = impactParameter;
605 279
606 const G4double effectiveImpactParameter = << 280 G4bool projectileHitsTarget = propagationModel->shootProjectile(projectile, impactParameter);
607 if(effectiveImpactParameter < 0.) { << 281 if(projectileHitsTarget == false) {
>> 282 // Increment the global counter for the number of transparents
>> 283 theGlobalInfo.nTransparents++;
>> 284
608 // Fill in the event information 285 // Fill in the event information
609 theEventInfo.transparent = true; 286 theEventInfo.transparent = true;
610 return false; << 287
>> 288 // Delete the projectile!
>> 289 delete projectile;
>> 290
>> 291 return theEventInfo;
611 } 292 }
612 293
613 // Fill in the event information 294 // Fill in the event information
614 theEventInfo.transparent = false; << 295 const G4double effectiveImpactParameter =
>> 296 projectile->getTransversePosition().mag();
615 theEventInfo.effectiveImpactParameter = ef 297 theEventInfo.effectiveImpactParameter = effectiveImpactParameter;
616 298
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 { 299 do {
626 // Run book keeping actions that should 300 // Run book keeping actions that should take place before propagation:
627 cascadeAction->beforePropagationAction(p << 301 propagationAction->beforePropagationAction(propagationModel);
628 302
629 // Get the avatar with the smallest time 303 // Get the avatar with the smallest time and propagate particles
630 // to that point in time. << 304 // to that poG4int in time.
631 IAvatar *avatar = propagationModel->prop << 305 G4INCL::IAvatar *avatar = propagationModel->propagate();
632 <<
633 finalState->reset(); <<
634 306
635 // Run book keeping actions that should 307 // Run book keeping actions that should take place after propagation:
636 cascadeAction->afterPropagationAction(pr << 308 propagationAction->afterPropagationAction(propagationModel, avatar);
637 309
638 if(avatar == 0) break; // No more avatar 310 if(avatar == 0) break; // No more avatars in the avatar list.
639 311
640 // Run book keeping actions that should 312 // Run book keeping actions that should take place before avatar:
641 cascadeAction->beforeAvatarAction(avatar << 313 avatarAction->beforeAvatarAction(avatar, nucleus);
642 314
643 // Channel is responsible for calculatin 315 // Channel is responsible for calculating the outcome of the
644 // selected avatar. There are different 316 // selected avatar. There are different kinds of channels. The
645 // class IChannel is, again, an abstract << 317 // class IChannel is, again, an abstract G4interface that defines
646 // the externally observable behavior of << 318 // the externally observable behavior of all G4interaction
647 // channels. 319 // channels.
648 // The handling of the channel is transp 320 // The handling of the channel is transparent to the API.
649 // Final state tells what changed... 321 // Final state tells what changed...
650 avatar->fillFinalState(finalState); << 322 G4INCL::FinalState *finalState = avatar->getFinalState();
>> 323
651 // Run book keeping actions that should 324 // Run book keeping actions that should take place after avatar:
652 cascadeAction->afterAvatarAction(avatar, << 325 avatarAction->afterAvatarAction(avatar, nucleus, finalState);
653 326
654 // So now we must give this information 327 // So now we must give this information to the nucleus
655 nucleus->applyFinalState(finalState); 328 nucleus->applyFinalState(finalState);
656 // and now we are ready to process the n 329 // and now we are ready to process the next avatar!
657 330
658 delete avatar; 331 delete avatar;
>> 332 delete finalState;
>> 333 } while(continueCascade());
659 334
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 335 // Fill in the event information
668 theEventInfo.stoppingTime = propagationMod << 336 theEventInfo.transparent = nucleus->isEventTransparent();
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 337
691 if(theEventInfo.transparent) { 338 if(theEventInfo.transparent) {
692 ProjectileRemnant * const projectileRemn << 339 // Increment the global counter for the number of transparents
693 if(projectileRemnant) { << 340 theGlobalInfo.nTransparents++;
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 { 341 } 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 342 // Check if the nucleus contains deltas
720 theEventInfo.deltasInside = nucleus->con 343 theEventInfo.deltasInside = nucleus->containsDeltas();
721 344
722 // Take care of any remaining deltas 345 // Take care of any remaining deltas
723 theEventInfo.forcedDeltasOutside = nucle 346 theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
724 theEventInfo.forcedDeltasInside = nucleu 347 theEventInfo.forcedDeltasInside = nucleus->decayInsideDeltas();
725 348
726 // Take care of any remaining etas, omeg << 349 // Cluster decay
727 G4double timeThreshold=theConfig->getDec << 350 theEventInfo.clusterDecay = nucleus->decayOutgoingClusters();
728 theEventInfo.forcedPionResonancesOutside << 351
729 nucleus->decayOutgoingSigmaZero(timeThre <<
730 nucleus->decayOutgoingNeutralKaon(); <<
731 <<
732 // Apply Coulomb distortion, if appropri 352 // Apply Coulomb distortion, if appropriate
733 // Note that this will apply Coulomb dis 353 // Note that this will apply Coulomb distortion also on pions emitted by
734 // unphysical remnants (see decayInsideD 354 // unphysical remnants (see decayInsideDeltas). This is at variance with
735 // what INCL4.6 does, but these events a 355 // what INCL4.6 does, but these events are (should be!) so rare that
736 // whatever we do doesn't (shouldn't!) m 356 // whatever we do doesn't (shouldn't!) make any noticeable difference.
737 CoulombDistortion::distortOut(nucleus->g << 357 G4INCL::CoulombDistortion::distortOut(nucleus->getStore()->getOutgoingParticles(), nucleus);
738 358
739 // If the normal cascade predicted compl << 359 // Compute recoil momentum, energy and spin of the nucleus
740 // masses to compute the excitation ener << 360 nucleus->computeRecoilKinematics();
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 361
771 #ifndef INCLXX_IN_GEANT4_MODE << 362 // Make room for the remnant recoil by rescaling the energies of the
772 // Global checks of conservation laws << 363 // outgoing particles.
773 globalConservationChecks(false); << 364 if(nucleus->hasRemnant()) rescaleOutgoingForRecoil();
774 #endif <<
775 <<
776 // Make room for the remnant recoil by <<
777 // outgoing particles. <<
778 if(nucleus->hasRemnant()) rescaleOutgo <<
779 365
780 } <<
781 <<
782 // Cluster decay <<
783 theEventInfo.clusterDecay = nucleus->dec <<
784 <<
785 #ifndef INCLXX_IN_GEANT4_MODE <<
786 // Global checks of conservation laws 366 // Global checks of conservation laws
787 globalConservationChecks(true); << 367 globalConservationChecks();
788 #endif <<
789 368
790 // Fill the EventInfo structure 369 // Fill the EventInfo structure
791 nucleus->fillEventInfo(&theEventInfo); 370 nucleus->fillEventInfo(&theEventInfo);
792 << 371 // theEventInfo.fillFromNucleus(nucleus);
793 } << 372 theEventInfo.stoppingTime = propagationModel->getCurrentTime();
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 } 373 }
871 374
872 // assert(theCNA==nucleus->getA()); << 375 return theEventInfo;
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 } 376 }
947 377
948 void INCL::rescaleOutgoingForRecoil() { 378 void INCL::rescaleOutgoingForRecoil() {
949 RecoilCMFunctor theRecoilFunctor(nucleus, << 379 Nucleus *nucleus = propagationModel->getNucleus();
950 380
951 // Apply the root-finding algorithm << 381 G4double sumKineticEnergies = 0.0;
952 const RootFinder::Solution theSolution = R << 382 // Sum up the kinetic energies of the outgoing particles
953 if(theSolution.success) { << 383 ParticleList outgoingParticles = nucleus->getStore()->getOutgoingParticles();
954 theRecoilFunctor(theSolution.x); // Appl << 384 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i )
955 } else { << 385 sumKineticEnergies += (*i)->getKineticEnergy();
956 INCL_WARN("Couldn't accommodate remnant << 386
>> 387 // If there is too little outgoing energy, we stop here.
>> 388 if(sumKineticEnergies <= 0.001) return;
>> 389 // The rescaling factor
>> 390 G4double rescale = 1. - nucleus->getRecoilEnergy()/sumKineticEnergies;
>> 391 if(rescale < 0.0) {
>> 392 WARN("Cannot accommodate remnant recoil by scaling outgoing energies. rescale = " << rescale << std::endl);
>> 393 rescale = 0.0;
>> 394 }
>> 395
>> 396 // Rescale the energies (and the momenta) of the outgoing particles.
>> 397 ThreeVector pBalance = nucleus->getIncomingMomentum();
>> 398 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i )
>> 399 {
>> 400 const G4double mass = (*i)->getMass();
>> 401 const G4double newKineticEnergy = (*i)->getKineticEnergy() * rescale;
>> 402
>> 403 (*i)->setEnergy(mass + newKineticEnergy);
>> 404 (*i)->adjustMomentumFromEnergy();
>> 405 //nucleus->updatePotentialEnergy(*i);
>> 406
>> 407 pBalance -= (*i)->getMomentum();
>> 408 }
>> 409
>> 410 nucleus->setRecoilMomentum(pBalance);
>> 411 const G4double remnantMass = ParticleTable::getMass(nucleus->getA(),nucleus->getZ()) + nucleus->getExcitationEnergy();
>> 412 const G4double pRem2 = pBalance.mag2();
>> 413 const G4double recoilEnergy = pRem2/
>> 414 (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
>> 415 nucleus->setRecoilEnergy(recoilEnergy);
>> 416 }
>> 417
>> 418 void INCL::globalConservationChecks() {
>> 419 const Nucleus *nucleus = propagationModel->getNucleus();
>> 420
>> 421 /* FIXME: This version of the energy-conservation check only uses kinetic
>> 422 energies, to mimic what INCL4.5 does. This is unsatisfactory because it
>> 423 does not take G4into account the particle masses. At some poG4int, it would
>> 424 be nice to have real energy conservation, with real masses. When ready
>> 425 to do so, have a look at the status of the code at commit
>> 426 aad75d09b8a52d28b8eb1bd38bdf347e63b802db (or possibly simply revert the
>> 427 following commit). */
>> 428
>> 429 // Initialise balance variables with the incoming values
>> 430 G4int ZBalance = theEventInfo.Zp + theEventInfo.Zt;
>> 431 G4int ABalance = theEventInfo.Ap + theEventInfo.At;
>> 432
>> 433 G4double projectileMass = 0.0;
>> 434 // FIXME: since we are not using total energies, we must set the projectile
>> 435 // mass to zero if the projectile is a pion.
>> 436 if(theEventInfo.projectileType != PiPlus &&
>> 437 theEventInfo.projectileType != PiZero &&
>> 438 theEventInfo.projectileType != PiMinus)
>> 439 projectileMass = ParticleTable::getMass(theEventInfo.projectileType);
>> 440
>> 441 G4double EBalance = nucleus->getInitialEnergy() - ParticleTable::getMass(theEventInfo.At, theEventInfo.Zt) - projectileMass;
>> 442 ThreeVector pBalance = nucleus->getIncomingMomentum();
>> 443
>> 444 // Process outgoing particles
>> 445 ParticleList outgoingParticles = nucleus->getStore()->getOutgoingParticles();
>> 446 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i ) {
>> 447 ZBalance -= (*i)->getZ();
>> 448 ABalance -= (*i)->getA();
>> 449 if((*i)->isPion()) // Ugly: we should calculate everything using total energies! (FIXME)
>> 450 EBalance -= (*i)->getEnergy();
>> 451 else
>> 452 EBalance -= (*i)->getKineticEnergy();
>> 453 pBalance -= (*i)->getMomentum();
957 } 454 }
958 } <<
959 455
960 #ifndef INCLXX_IN_GEANT4_MODE << 456 EBalance -= nucleus->computeSeparationEnergyBalance();
961 void INCL::globalConservationChecks(G4bool a <<
962 Nucleus::ConservationBalance theBalance = <<
963 457
964 // Global conservation checks << 458 // Remnant contribution, if present
965 const G4double pLongBalance = theBalance.m << 459 if(nucleus->hasRemnant()) {
966 const G4double pTransBalance = theBalance. << 460 ZBalance -= nucleus->getZ();
967 if(theBalance.Z != 0) { << 461 ABalance -= nucleus->getA();
968 INCL_ERROR("Violation of charge conserva << 462 EBalance -= //ParticleTable::getMass(nucleus->getA(),nucleus->getZ()) +
969 } << 463 nucleus->getExcitationEnergy() + nucleus->getRecoilEnergy();
970 if(theBalance.A != 0) { << 464 pBalance -= nucleus->getRecoilMomentum();
971 INCL_ERROR("Violation of baryon-number c <<
972 } 465 }
973 if(theBalance.S != 0) { << 466
974 INCL_ERROR("Violation of strange-number << 467 // Global conservation checks
975 } << 468 const G4double pLongBalance = pBalance.getZ();
976 G4double EThreshold, pLongThreshold, pTran << 469 const G4double pTransBalance = pBalance.perp();
977 if(afterRecoil) { << 470 if(ZBalance != 0) {
978 // Less stringent checks after accommoda << 471 ERROR("Violation of charge conservation! ZBalance = " << ZBalance << std::endl);
979 EThreshold = 10.; // MeV << 472 }
980 pLongThreshold = 1.; // MeV/c << 473 if(ABalance != 0) {
981 pTransThreshold = 1.; // MeV/c << 474 ERROR("Violation of baryon-number conservation! ABalance = " << ABalance << std::endl);
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 } 475 }
988 if(std::abs(theBalance.energy)>EThreshold) << 476 if(std::abs(EBalance)>10.0) {
989 INCL_WARN("Violation of energy conservat << 477 WARN("Violation of energy conservation > 10 MeV. EBalance = " << EBalance << std::endl);
990 } 478 }
991 if(std::abs(pLongBalance)>pLongThreshold) << 479 if(std::abs(pLongBalance)>5.0) {
992 INCL_WARN("Violation of longitudinal mom << 480 WARN("Violation of longitudinal momentum conservation > 5.0 MeV. pLongBalance = " << pLongBalance << std::endl);
993 } 481 }
994 if(std::abs(pTransBalance)>pTransThreshold << 482 if(std::abs(pTransBalance)>5.0) {
995 INCL_WARN("Violation of transverse momen << 483 WARN("Violation of transverse momentum conservation > 5.0 MeV. pTransBalance = " << pTransBalance << std::endl);
996 } 484 }
997 485
998 // Feed the EventInfo variables 486 // Feed the EventInfo variables
999 theEventInfo.EBalance = theBalance.energy; << 487 theEventInfo.EBalance = EBalance;
1000 theEventInfo.pLongBalance = pLongBalance; 488 theEventInfo.pLongBalance = pLongBalance;
1001 theEventInfo.pTransBalance = pTransBalanc 489 theEventInfo.pTransBalance = pTransBalance;
1002 } 490 }
1003 #endif <<
1004 491
1005 G4bool INCL::continueCascade() { 492 G4bool INCL::continueCascade() {
>> 493 Nucleus *nucleus = propagationModel->getNucleus();
1006 // Stop if we have passed the stopping ti 494 // Stop if we have passed the stopping time
1007 if(propagationModel->getCurrentTime() > p << 495 if(propagationModel->getCurrentTime() > propagationModel->getStoppingTime()) return false;
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 496 // Stop if there are no participants and no pions inside the nucleus
1014 if(nucleus->getStore()->getBook().getCasc << 497 if(nucleus->getStore()->getBook()->getParticipants()==0 &&
1015 nucleus->getStore()->getIncomingParti << 498 nucleus->getStore()->getIncomingParticles().empty()) return false;
1016 INCL_DEBUG("No participants in the nucl << 499 // Stop if the remnant has only one nucleon
1017 return false; << 500 if(nucleus->getA() <= 1) 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; 501 return true;
1034 } 502 }
1035 503
1036 void INCL::finalizeGlobalInfo(Random::SeedV << 504 void INCL::finaliseGlobalInfo() {
1037 const G4double normalisationFactor = theG << 505 theGlobalInfo.reactionCrossSection = theGlobalInfo.geometricCrossSection *
>> 506 ((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents)) /
>> 507 ((G4double) theGlobalInfo.nShots);
>> 508 theGlobalInfo.errorReactionCrossSection = theGlobalInfo.geometricCrossSection *
>> 509 std::sqrt((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents)) /
1038 ((G4double) theGlobalInfo.nShots); 510 ((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 } 511 }
1390 512
1391 } 513 }
1392 514