| 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 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 /** \file G4INCLCascade.cc 38 /** \file G4INCLCascade.cc
39 * 39 *
40 * INCL Cascade 40 * INCL Cascade
41 */ 41 */
42 #include "G4INCLCascade.hh" 42 #include "G4INCLCascade.hh"
43 #include "G4INCLRandom.hh" 43 #include "G4INCLRandom.hh"
44 #include "G4INCLStandardPropagationModel.hh" 44 #include "G4INCLStandardPropagationModel.hh"
45 #include "G4INCLParticleTable.hh" 45 #include "G4INCLParticleTable.hh"
46 #include "G4INCLParticle.hh" <<
47 #include "G4INCLNuclearMassTable.hh" 46 #include "G4INCLNuclearMassTable.hh"
48 #include "G4INCLGlobalInfo.hh" 47 #include "G4INCLGlobalInfo.hh"
49 #include "G4INCLNucleus.hh" <<
50 48
51 #include "G4INCLPauliBlocking.hh" 49 #include "G4INCLPauliBlocking.hh"
52 50
53 #include "G4INCLCrossSections.hh" 51 #include "G4INCLCrossSections.hh"
54 52
55 #include "G4INCLPhaseSpaceGenerator.hh" 53 #include "G4INCLPhaseSpaceGenerator.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" 59 #include "G4INCLINuclearPotential.hh"
62 60
63 #include "G4INCLCoulombDistortion.hh" 61 #include "G4INCLCoulombDistortion.hh"
64 62
65 #include "G4INCLClustering.hh" 63 #include "G4INCLClustering.hh"
66 64
67 #include "G4INCLIntersection.hh" 65 #include "G4INCLIntersection.hh"
68 66
69 #include "G4INCLBinaryCollisionAvatar.hh" 67 #include "G4INCLBinaryCollisionAvatar.hh"
70 68
71 #include "G4INCLCascadeAction.hh" 69 #include "G4INCLCascadeAction.hh"
72 #include "G4INCLAvatarDumpAction.hh" 70 #include "G4INCLAvatarDumpAction.hh"
73 71
74 #include <cstring> << 72 #include <cstring>
75 #include <cstdlib> 73 #include <cstdlib>
76 #include <numeric> 74 #include <numeric>
77 75
78 #include "G4INCLPbarAtrestEntryChannel.hh" <<
79 <<
80 namespace G4INCL { 76 namespace G4INCL {
81 << 77
82 INCL::INCL(Config const * const config) 78 INCL::INCL(Config const * const config)
83 :propagationModel(0), theA(208), theZ(82), << 79 :propagationModel(0), theA(208), theZ(82),
84 targetInitSuccess(false), 80 targetInitSuccess(false),
85 maxImpactParameter(0.), 81 maxImpactParameter(0.),
86 maxUniverseRadius(0.), 82 maxUniverseRadius(0.),
87 maxInteractionDistance(0.), 83 maxInteractionDistance(0.),
88 fixedImpactParameter(0.), 84 fixedImpactParameter(0.),
89 theConfig(config), 85 theConfig(config),
90 nucleus(NULL), 86 nucleus(NULL),
91 forceTransparent(false), 87 forceTransparent(false),
92 minRemnantSize(4) 88 minRemnantSize(4)
93 { 89 {
94 // Set the logger object. 90 // Set the logger object.
95 #ifdef INCLXX_IN_GEANT4_MODE 91 #ifdef INCLXX_IN_GEANT4_MODE
96 Logger::initVerbosityLevelFromEnvvar(); 92 Logger::initVerbosityLevelFromEnvvar();
97 #else // INCLXX_IN_GEANT4_MODE 93 #else // INCLXX_IN_GEANT4_MODE
98 Logger::initialize(theConfig); 94 Logger::initialize(theConfig);
99 #endif // INCLXX_IN_GEANT4_MODE 95 #endif // INCLXX_IN_GEANT4_MODE
100 96
101 // Set the random number generator algorit 97 // Set the random number generator algorithm. The system can support
102 // multiple different generator algorithms 98 // multiple different generator algorithms in a completely
103 // transparent way. 99 // transparent way.
104 Random::initialize(theConfig); 100 Random::initialize(theConfig);
105 101
106 // Select the Pauli and CDPP blocking algo 102 // Select the Pauli and CDPP blocking algorithms
107 Pauli::initialize(theConfig); 103 Pauli::initialize(theConfig);
108 104
109 // Set the cross-section set 105 // Set the cross-section set
110 CrossSections::initialize(theConfig); 106 CrossSections::initialize(theConfig);
111 107
112 // Set the phase-space generator 108 // Set the phase-space generator
113 PhaseSpaceGenerator::initialize(theConfig) 109 PhaseSpaceGenerator::initialize(theConfig);
114 110
115 // Select the Coulomb-distortion algorithm 111 // Select the Coulomb-distortion algorithm:
116 CoulombDistortion::initialize(theConfig); 112 CoulombDistortion::initialize(theConfig);
117 113
118 // Select the clustering algorithm: 114 // Select the clustering algorithm:
119 Clustering::initialize(theConfig); 115 Clustering::initialize(theConfig);
120 116
121 // Initialize the INCL particle table: 117 // Initialize the INCL particle table:
122 ParticleTable::initialize(theConfig); 118 ParticleTable::initialize(theConfig);
123 119
124 // Initialize the value of cutNN in Binary 120 // Initialize the value of cutNN in BinaryCollisionAvatar
125 BinaryCollisionAvatar::setCutNN(theConfig- 121 BinaryCollisionAvatar::setCutNN(theConfig->getCutNN());
126 122
127 // Initialize the value of strange cross s <<
128 BinaryCollisionAvatar::setBias(theConfig-> <<
129 <<
130 // Propagation model is responsible for fi 123 // Propagation model is responsible for finding avatars and
131 // transporting the particles. In principl 124 // transporting the particles. In principle this step is "hidden"
132 // behind an abstract interface and the re 125 // behind an abstract interface and the rest of the system does not
133 // care how the transportation and avatar 126 // care how the transportation and avatar finding is done. This
134 // should allow us to "easily" experiment 127 // should allow us to "easily" experiment with different avatar
135 // finding schemes and even to support thi 128 // finding schemes and even to support things like curved
136 // trajectories in the future. 129 // trajectories in the future.
137 propagationModel = new StandardPropagation 130 propagationModel = new StandardPropagationModel(theConfig->getLocalEnergyBBType(),theConfig->getLocalEnergyPiType(),theConfig->getHadronizationTime());
138 if(theConfig->getCascadeActionType() == Av 131 if(theConfig->getCascadeActionType() == AvatarDumpActionType)
139 cascadeAction = new AvatarDumpAction(); 132 cascadeAction = new AvatarDumpAction();
140 else 133 else
141 cascadeAction = new CascadeAction(); 134 cascadeAction = new CascadeAction();
142 cascadeAction->beforeRunAction(theConfig); 135 cascadeAction->beforeRunAction(theConfig);
143 136
144 theGlobalInfo.cascadeModel = theConfig->ge 137 theGlobalInfo.cascadeModel = theConfig->getVersionString();
145 theGlobalInfo.deexcitationModel = theConfi 138 theGlobalInfo.deexcitationModel = theConfig->getDeExcitationString();
146 #ifdef INCL_ROOT_USE 139 #ifdef INCL_ROOT_USE
147 theGlobalInfo.rootSelection = theConfig->g 140 theGlobalInfo.rootSelection = theConfig->getROOTSelectionString();
148 #endif 141 #endif
149 142
150 #ifndef INCLXX_IN_GEANT4_MODE 143 #ifndef INCLXX_IN_GEANT4_MODE
151 // Fill in the global information 144 // Fill in the global information
152 theGlobalInfo.At = theConfig->getTargetA() 145 theGlobalInfo.At = theConfig->getTargetA();
153 theGlobalInfo.Zt = theConfig->getTargetZ() 146 theGlobalInfo.Zt = theConfig->getTargetZ();
154 theGlobalInfo.St = theConfig->getTargetS() <<
155 const ParticleSpecies theSpecies = theConf 147 const ParticleSpecies theSpecies = theConfig->getProjectileSpecies();
156 theGlobalInfo.Ap = theSpecies.theA; 148 theGlobalInfo.Ap = theSpecies.theA;
157 theGlobalInfo.Zp = theSpecies.theZ; 149 theGlobalInfo.Zp = theSpecies.theZ;
158 theGlobalInfo.Sp = theSpecies.theS; <<
159 theGlobalInfo.Ep = theConfig->getProjectil 150 theGlobalInfo.Ep = theConfig->getProjectileKineticEnergy();
160 theGlobalInfo.biasFactor = theConfig->getB <<
161 #endif 151 #endif
162 152
163 fixedImpactParameter = theConfig->getImpac 153 fixedImpactParameter = theConfig->getImpactParameter();
164 } 154 }
165 155
166 INCL::~INCL() { 156 INCL::~INCL() {
167 InteractionAvatar::deleteBackupParticles() 157 InteractionAvatar::deleteBackupParticles();
168 #ifndef INCLXX_IN_GEANT4_MODE 158 #ifndef INCLXX_IN_GEANT4_MODE
169 NuclearMassTable::deleteTable(); 159 NuclearMassTable::deleteTable();
170 #endif 160 #endif
171 PhaseSpaceGenerator::deletePhaseSpaceGener 161 PhaseSpaceGenerator::deletePhaseSpaceGenerator();
172 CrossSections::deleteCrossSections(); 162 CrossSections::deleteCrossSections();
173 Pauli::deleteBlockers(); 163 Pauli::deleteBlockers();
174 CoulombDistortion::deleteCoulomb(); 164 CoulombDistortion::deleteCoulomb();
175 Random::deleteGenerator(); 165 Random::deleteGenerator();
176 Clustering::deleteClusteringModel(); 166 Clustering::deleteClusteringModel();
177 #ifndef INCLXX_IN_GEANT4_MODE 167 #ifndef INCLXX_IN_GEANT4_MODE
178 Logger::deleteLoggerSlave(); 168 Logger::deleteLoggerSlave();
179 #endif 169 #endif
180 NuclearDensityFactory::clearCache(); 170 NuclearDensityFactory::clearCache();
181 NuclearPotential::clearCache(); 171 NuclearPotential::clearCache();
182 cascadeAction->afterRunAction(); 172 cascadeAction->afterRunAction();
183 delete cascadeAction; 173 delete cascadeAction;
184 delete propagationModel; 174 delete propagationModel;
185 delete theConfig; 175 delete theConfig;
186 } 176 }
187 177
188 G4bool INCL::prepareReaction(const ParticleS << 178 G4bool INCL::prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z) {
189 if(A < 0 || A > 300 || Z < 1 || Z > 200) { 179 if(A < 0 || A > 300 || Z < 1 || Z > 200) {
190 INCL_ERROR("Unsupported target: A = " << << 180 INCL_ERROR("Unsupported target: A = " << A << " Z = " << Z << '\n'
191 << "Target configuration reje 181 << "Target configuration rejected." << '\n');
192 return false; 182 return false;
193 } 183 }
194 if(projectileSpecies.theType==Composite && 184 if(projectileSpecies.theType==Composite &&
195 (projectileSpecies.theZ==projectileSpec 185 (projectileSpecies.theZ==projectileSpecies.theA || projectileSpecies.theZ==0)) {
196 INCL_ERROR("Unsupported projectile: A = << 186 INCL_ERROR("Unsupported projectile: A = " << projectileSpecies.theA << " Z = " << projectileSpecies.theZ << '\n'
197 << "Projectile configuration 187 << "Projectile configuration rejected." << '\n');
198 return false; 188 return false;
199 } 189 }
200 190
201 // Reset the forced-transparent flag 191 // Reset the forced-transparent flag
202 forceTransparent = false; 192 forceTransparent = false;
203 193
204 // Initialise the maximum universe radius 194 // Initialise the maximum universe radius
205 initUniverseRadius(projectileSpecies, kine 195 initUniverseRadius(projectileSpecies, kineticEnergy, A, Z);
206 // Initialise the nucleus <<
207 196
208 //D << 197 // Initialise the nucleus
209 //reset << 198 theZ = Z;
210 G4bool ProtonIsTheVictim = false; << 199 if(theConfig->isNaturalTarget())
211 G4bool NeutronIsTheVictim = false; << 200 theA = ParticleTable::drawRandomNaturalIsotope(Z);
212 theEventInfo.annihilationP = false; << 201 else
213 theEventInfo.annihilationN = false; << 202 theA = A;
214 << 203 initializeTarget(theA, theZ);
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 204
262 initializeTarget(theA, theZ, theS, theATyp <<
263 <<
264 // Set the maximum impact parameter 205 // Set the maximum impact parameter
265 maxImpactParameter = CoulombDistortion::ma 206 maxImpactParameter = CoulombDistortion::maxImpactParameter(projectileSpecies, kineticEnergy, nucleus);
266 INCL_DEBUG("Maximum impact parameter initi 207 INCL_DEBUG("Maximum impact parameter initialised: " << maxImpactParameter << '\n');
267 208
268 // For forced CN events 209 // For forced CN events
269 initMaxInteractionDistance(projectileSpeci 210 initMaxInteractionDistance(projectileSpecies, kineticEnergy);
270 // Set the geometric cross sectiony section << 211
271 if(projectileSpecies.theType == antiProton << 212 // Set the geometric cross section
272 G4int currentA = A; << 213 theGlobalInfo.geometricCrossSection =
273 if(theConfig->isNaturalTarget()){ << 214 Math::tenPi*std::pow(maxImpactParameter,2);
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 215
288 // Set the minimum remnant size 216 // Set the minimum remnant size
289 if(projectileSpecies.theA > 0) 217 if(projectileSpecies.theA > 0)
290 minRemnantSize = std::min(theA, 4); 218 minRemnantSize = std::min(theA, 4);
291 else 219 else
292 minRemnantSize = std::min(theA-1, 4); 220 minRemnantSize = std::min(theA-1, 4);
>> 221
293 return true; 222 return true;
294 } 223 }
295 224
296 G4bool INCL::initializeTarget(const G4int A, << 225 G4bool INCL::initializeTarget(const G4int A, const G4int Z) {
297 delete nucleus; 226 delete nucleus;
298 227
299 if (theAType==PType || theAType==NType) { << 228 nucleus = new Nucleus(A, Z, theConfig, maxUniverseRadius);
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(); 229 nucleus->getStore()->getBook().reset();
309 nucleus->initializeParticles(); 230 nucleus->initializeParticles();
>> 231
310 propagationModel->setNucleus(nucleus); 232 propagationModel->setNucleus(nucleus);
311 return true; 233 return true;
312 } 234 }
313 235
314 const EventInfo &INCL::processEvent( 236 const EventInfo &INCL::processEvent(
315 ParticleSpecies const &projectileSpecies 237 ParticleSpecies const &projectileSpecies,
316 const G4double kineticEnergy, 238 const G4double kineticEnergy,
317 const G4int targetA, 239 const G4int targetA,
318 const G4int targetZ, << 240 const G4int targetZ
319 const G4int targetS <<
320 ) { 241 ) {
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 242 // Set the target and the projectile
528 targetInitSuccess = prepareReaction(projec << 243 targetInitSuccess = prepareReaction(projectileSpecies, kineticEnergy, targetA, targetZ);
529 244
530 if(!targetInitSuccess) { 245 if(!targetInitSuccess) {
531 INCL_WARN("Target initialisation failed << 246 INCL_WARN("Target initialisation failed for A=" << targetA << ", Z=" << targetZ << '\n');
532 theEventInfo.transparent=true; 247 theEventInfo.transparent=true;
533 return theEventInfo; 248 return theEventInfo;
534 } 249 }
535 250
536 cascadeAction->beforeCascadeAction(propaga 251 cascadeAction->beforeCascadeAction(propagationModel);
537 252
538 const G4bool canRunCascade = preCascade(pr 253 const G4bool canRunCascade = preCascade(projectileSpecies, kineticEnergy);
539 if(canRunCascade) { 254 if(canRunCascade) {
540 cascade(); 255 cascade();
541 postCascade(projectileSpecies, kineticEn << 256 postCascade();
542 cascadeAction->afterCascadeAction(nucleu 257 cascadeAction->afterCascadeAction(nucleus);
543 } 258 }
544 updateGlobalInfo(); 259 updateGlobalInfo();
545 return theEventInfo; 260 return theEventInfo;
546 } 261 }
547 262
548 G4bool INCL::preCascade(ParticleSpecies cons 263 G4bool INCL::preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
549 // Reset theEventInfo 264 // Reset theEventInfo
550 theEventInfo.reset(); 265 theEventInfo.reset();
551 << 266
552 EventInfo::eventNumber++; 267 EventInfo::eventNumber++;
553 268
554 // Fill in the event information 269 // Fill in the event information
555 theEventInfo.projectileType = projectileSp 270 theEventInfo.projectileType = projectileSpecies.theType;
556 theEventInfo.Ap = (Short_t)projectileSpeci << 271 theEventInfo.Ap = projectileSpecies.theA;
557 theEventInfo.Zp = (Short_t)projectileSpeci << 272 theEventInfo.Zp = projectileSpecies.theZ;
558 theEventInfo.Sp = (Short_t)projectileSpeci <<
559 theEventInfo.Ep = kineticEnergy; 273 theEventInfo.Ep = kineticEnergy;
560 theEventInfo.St = (Short_t)nucleus->getS() << 274 theEventInfo.At = nucleus->getA();
>> 275 theEventInfo.Zt = nucleus->getZ();
561 276
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 277 // Do nothing below the Coulomb barrier
577 if(maxImpactParameter<=0.) { 278 if(maxImpactParameter<=0.) {
578 // Fill in the event information 279 // Fill in the event information
579 //Particle *pbar = new Particle; << 280 theEventInfo.transparent = true;
580 //PbarAtrestEntryChannel *obj = new PbarAt << 281
581 if(projectileSpecies.theType == antiProt << 282 return false;
582 INCL_DEBUG("at rest annihilation" << ' <<
583 //theEventInfo.transparent = false; <<
584 } else { <<
585 theEventInfo.transparent = true; <<
586 return false; <<
587 } <<
588 } 283 }
589 <<
590 284
591 // Randomly draw an impact parameter or us 285 // Randomly draw an impact parameter or use a fixed value, depending on the
592 // Config option 286 // Config option
593 G4double impactParameter, phi; 287 G4double impactParameter, phi;
594 if(fixedImpactParameter<0.) { 288 if(fixedImpactParameter<0.) {
595 impactParameter = maxImpactParameter * s 289 impactParameter = maxImpactParameter * std::sqrt(Random::shoot0());
596 phi = Random::shoot() * Math::twoPi; 290 phi = Random::shoot() * Math::twoPi;
597 } else { 291 } else {
598 impactParameter = fixedImpactParameter; 292 impactParameter = fixedImpactParameter;
599 phi = 0.; 293 phi = 0.;
600 } 294 }
601 INCL_DEBUG("Selected impact parameter: " < 295 INCL_DEBUG("Selected impact parameter: " << impactParameter << '\n');
602 296
603 // Fill in the event information 297 // Fill in the event information
604 theEventInfo.impactParameter = impactParam 298 theEventInfo.impactParameter = impactParameter;
605 299
606 const G4double effectiveImpactParameter = 300 const G4double effectiveImpactParameter = propagationModel->shoot(projectileSpecies, kineticEnergy, impactParameter, phi);
607 if(effectiveImpactParameter < 0.) { 301 if(effectiveImpactParameter < 0.) {
608 // Fill in the event information 302 // Fill in the event information
609 theEventInfo.transparent = true; 303 theEventInfo.transparent = true;
>> 304
610 return false; 305 return false;
611 } 306 }
612 307
613 // Fill in the event information 308 // Fill in the event information
614 theEventInfo.transparent = false; 309 theEventInfo.transparent = false;
615 theEventInfo.effectiveImpactParameter = ef 310 theEventInfo.effectiveImpactParameter = effectiveImpactParameter;
616 311
617 return true; 312 return true;
618 } 313 }
619 314
620 void INCL::cascade() { 315 void INCL::cascade() {
621 FinalState *finalState = new FinalState; 316 FinalState *finalState = new FinalState;
622 317
623 unsigned long loopCounter = 0; 318 unsigned long loopCounter = 0;
624 const unsigned long maxLoopCounter = 10000 319 const unsigned long maxLoopCounter = 10000000;
625 do { 320 do {
626 // Run book keeping actions that should 321 // Run book keeping actions that should take place before propagation:
627 cascadeAction->beforePropagationAction(p 322 cascadeAction->beforePropagationAction(propagationModel);
628 323
629 // Get the avatar with the smallest time 324 // Get the avatar with the smallest time and propagate particles
630 // to that point in time. 325 // to that point in time.
631 IAvatar *avatar = propagationModel->prop 326 IAvatar *avatar = propagationModel->propagate(finalState);
632 327
633 finalState->reset(); 328 finalState->reset();
634 329
635 // Run book keeping actions that should 330 // Run book keeping actions that should take place after propagation:
636 cascadeAction->afterPropagationAction(pr 331 cascadeAction->afterPropagationAction(propagationModel, avatar);
637 332
638 if(avatar == 0) break; // No more avatar 333 if(avatar == 0) break; // No more avatars in the avatar list.
639 334
640 // Run book keeping actions that should 335 // Run book keeping actions that should take place before avatar:
641 cascadeAction->beforeAvatarAction(avatar 336 cascadeAction->beforeAvatarAction(avatar, nucleus);
642 337
643 // Channel is responsible for calculatin 338 // Channel is responsible for calculating the outcome of the
644 // selected avatar. There are different 339 // selected avatar. There are different kinds of channels. The
645 // class IChannel is, again, an abstract 340 // class IChannel is, again, an abstract interface that defines
646 // the externally observable behavior of 341 // the externally observable behavior of all interaction
647 // channels. 342 // channels.
648 // The handling of the channel is transp 343 // The handling of the channel is transparent to the API.
649 // Final state tells what changed... 344 // Final state tells what changed...
650 avatar->fillFinalState(finalState); 345 avatar->fillFinalState(finalState);
>> 346
651 // Run book keeping actions that should 347 // Run book keeping actions that should take place after avatar:
652 cascadeAction->afterAvatarAction(avatar, 348 cascadeAction->afterAvatarAction(avatar, nucleus, finalState);
653 349
654 // So now we must give this information 350 // So now we must give this information to the nucleus
655 nucleus->applyFinalState(finalState); 351 nucleus->applyFinalState(finalState);
656 // and now we are ready to process the n 352 // and now we are ready to process the next avatar!
657 353
658 delete avatar; 354 delete avatar;
659 355
660 ++loopCounter; 356 ++loopCounter;
661 } while(continueCascade() && loopCounter<m 357 } while(continueCascade() && loopCounter<maxLoopCounter); /* Loop checking, 10.07.2015, D.Mancusi */
662 << 358
663 delete finalState; 359 delete finalState;
664 } 360 }
665 361
666 void INCL::postCascade(const ParticleSpecies << 362 void INCL::postCascade() {
667 // Fill in the event information 363 // Fill in the event information
668 theEventInfo.stoppingTime = propagationMod 364 theEventInfo.stoppingTime = propagationModel->getCurrentTime();
669 365
670 // The event bias <<
671 theEventInfo.eventBias = (Double_t) Partic <<
672 <<
673 // Forced CN? 366 // Forced CN?
674 if(!(projectileSpecies.theType==antiProton << 367 if(nucleus->getTryCompoundNucleus()) {
675 if(nucleus->getTryCompoundNucleus()) { << 368 INCL_DEBUG("Trying compound nucleus" << '\n');
676 INCL_DEBUG("Trying compound nucleus" < << 369 makeCompoundNucleus();
677 makeCompoundNucleus(); << 370 theEventInfo.transparent = forceTransparent;
678 theEventInfo.transparent = forceTransp <<
679 // Global checks of conservation laws 371 // Global checks of conservation laws
680 #ifndef INCLXX_IN_GEANT4_MODE 372 #ifndef INCLXX_IN_GEANT4_MODE
681 if(!theEventInfo.transparent) globalCons 373 if(!theEventInfo.transparent) globalConservationChecks(true);
682 #endif 374 #endif
683 return; 375 return;
684 } <<
685 } 376 }
686 377
687 if(!(projectileSpecies.theType==antiProton << 378 theEventInfo.transparent = forceTransparent || nucleus->isEventTransparent();
688 theEventInfo.transparent = forceTranspar <<
689 } <<
690 379
691 if(theEventInfo.transparent) { 380 if(theEventInfo.transparent) {
692 ProjectileRemnant * const projectileRemn 381 ProjectileRemnant * const projectileRemnant = nucleus->getProjectileRemnant();
693 if(projectileRemnant) { 382 if(projectileRemnant) {
694 // Clear the incoming list (particles 383 // Clear the incoming list (particles will be deleted by the ProjectileRemnant)
695 nucleus->getStore()->clearIncoming(); 384 nucleus->getStore()->clearIncoming();
696 } else { 385 } else {
697 // Delete particles in the incoming li 386 // Delete particles in the incoming list
698 nucleus->getStore()->deleteIncoming(); 387 nucleus->getStore()->deleteIncoming();
699 } 388 }
700 } else { 389 } 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 390 // Check if the nucleus contains deltas
720 theEventInfo.deltasInside = nucleus->con 391 theEventInfo.deltasInside = nucleus->containsDeltas();
721 392
722 // Take care of any remaining deltas 393 // Take care of any remaining deltas
723 theEventInfo.forcedDeltasOutside = nucle 394 theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
724 theEventInfo.forcedDeltasInside = nucleu 395 theEventInfo.forcedDeltasInside = nucleus->decayInsideDeltas();
725 396
726 // Take care of any remaining etas, omeg << 397 // Take care of any remaining etas and/or omegas
727 G4double timeThreshold=theConfig->getDec << 398 /* if (theConfig->getDecayTimeThreshold() > 7.e-19) {
728 theEventInfo.forcedPionResonancesOutside << 399 theEventInfo.forcedPionResonancesOutside = nucleus->decayOutgoingPionResonances();
729 nucleus->decayOutgoingSigmaZero(timeThre << 400 }*/
730 nucleus->decayOutgoingNeutralKaon(); << 401 G4double timeThreshold=theConfig->getDecayTimeThreshold();
731 << 402 theEventInfo.forcedPionResonancesOutside = nucleus->decayOutgoingPionResonances(timeThreshold);
>> 403
732 // Apply Coulomb distortion, if appropri 404 // Apply Coulomb distortion, if appropriate
733 // Note that this will apply Coulomb dis 405 // Note that this will apply Coulomb distortion also on pions emitted by
734 // unphysical remnants (see decayInsideD 406 // unphysical remnants (see decayInsideDeltas). This is at variance with
735 // what INCL4.6 does, but these events a 407 // what INCL4.6 does, but these events are (should be!) so rare that
736 // whatever we do doesn't (shouldn't!) m 408 // whatever we do doesn't (shouldn't!) make any noticeable difference.
737 CoulombDistortion::distortOut(nucleus->g 409 CoulombDistortion::distortOut(nucleus->getStore()->getOutgoingParticles(), nucleus);
738 410
739 // If the normal cascade predicted compl 411 // If the normal cascade predicted complete fusion, use the tabulated
740 // masses to compute the excitation ener 412 // masses to compute the excitation energy, the recoil, etc.
741 if(nucleus->getStore()->getOutgoingParti 413 if(nucleus->getStore()->getOutgoingParticles().size()==0
742 && (!nucleus->getProjectileRemnant() 414 && (!nucleus->getProjectileRemnant()
743 || nucleus->getProjectileRemnant( 415 || nucleus->getProjectileRemnant()->getParticles().size()==0)) {
744 416
745 INCL_DEBUG("Cascade resulted in comple 417 INCL_DEBUG("Cascade resulted in complete fusion, using realistic fusion kinematics" << '\n');
746 418
747 nucleus->useFusionKinematics(); 419 nucleus->useFusionKinematics();
748 420
749 if(nucleus->getExcitationEnergy()<0.) 421 if(nucleus->getExcitationEnergy()<0.) {
750 // Complete fusion is energetically 422 // Complete fusion is energetically impossible, return a transparent
751 INCL_WARN("Complete-fusion kinematic 423 INCL_WARN("Complete-fusion kinematics yields negative excitation energy, returning a transparent!" << '\n');
752 theEventInfo.transparent = true; 424 theEventInfo.transparent = true;
753 return; 425 return;
754 } 426 }
755 427
756 } else { // Normal cascade here 428 } else { // Normal cascade here
757 429
758 // Set the excitation energy 430 // Set the excitation energy
759 nucleus->setExcitationEnergy(nucleus-> 431 nucleus->setExcitationEnergy(nucleus->computeExcitationEnergy());
760 432
761 // Make a projectile pre-fragment out 433 // Make a projectile pre-fragment out of the geometrical and dynamical
762 // spectators 434 // spectators
763 theEventInfo.nUnmergedSpectators = mak 435 theEventInfo.nUnmergedSpectators = makeProjectileRemnant();
764 436
765 // Compute recoil momentum, energy and 437 // Compute recoil momentum, energy and spin of the nucleus
766 if(nucleus->getA()==1 && minRemnantSiz 438 if(nucleus->getA()==1 && minRemnantSize>1) {
767 INCL_ERROR("Computing one-nucleon re 439 INCL_ERROR("Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << '\n');
768 } 440 }
769 nucleus->computeRecoilKinematics(); 441 nucleus->computeRecoilKinematics();
770 442
771 #ifndef INCLXX_IN_GEANT4_MODE 443 #ifndef INCLXX_IN_GEANT4_MODE
772 // Global checks of conservation laws 444 // Global checks of conservation laws
773 globalConservationChecks(false); 445 globalConservationChecks(false);
774 #endif 446 #endif
775 447
776 // Make room for the remnant recoil by 448 // Make room for the remnant recoil by rescaling the energies of the
777 // outgoing particles. 449 // outgoing particles.
778 if(nucleus->hasRemnant()) rescaleOutgo 450 if(nucleus->hasRemnant()) rescaleOutgoingForRecoil();
779 451
780 } 452 }
781 453
782 // Cluster decay 454 // Cluster decay
783 theEventInfo.clusterDecay = nucleus->dec << 455 theEventInfo.clusterDecay = nucleus->decayOutgoingClusters() | nucleus->decayMe();
784 456
785 #ifndef INCLXX_IN_GEANT4_MODE 457 #ifndef INCLXX_IN_GEANT4_MODE
786 // Global checks of conservation laws 458 // Global checks of conservation laws
787 globalConservationChecks(true); 459 globalConservationChecks(true);
788 #endif 460 #endif
789 461
790 // Fill the EventInfo structure 462 // Fill the EventInfo structure
791 nucleus->fillEventInfo(&theEventInfo); 463 nucleus->fillEventInfo(&theEventInfo);
792 464
793 } 465 }
794 } 466 }
795 467
796 void INCL::makeCompoundNucleus() { 468 void INCL::makeCompoundNucleus() {
797 // If this is not a nucleus-nucleus collis 469 // If this is not a nucleus-nucleus collision, don't attempt to make a
798 // compound nucleus. 470 // compound nucleus.
799 // 471 //
800 // Yes, even nucleon-nucleus collisions ca 472 // Yes, even nucleon-nucleus collisions can lead to particles entering
801 // below the Fermi level. Take e.g. 1-MeV 473 // below the Fermi level. Take e.g. 1-MeV p + He4.
802 if(!nucleus->isNucleusNucleusCollision()) 474 if(!nucleus->isNucleusNucleusCollision()) {
803 forceTransparent = true; 475 forceTransparent = true;
804 return; 476 return;
805 } 477 }
806 478
807 // Reset the internal Nucleus variables 479 // Reset the internal Nucleus variables
808 nucleus->getStore()->clearIncoming(); 480 nucleus->getStore()->clearIncoming();
809 nucleus->getStore()->clearOutgoing(); 481 nucleus->getStore()->clearOutgoing();
810 nucleus->getProjectileRemnant()->reset(); 482 nucleus->getProjectileRemnant()->reset();
811 nucleus->setA(theEventInfo.At); 483 nucleus->setA(theEventInfo.At);
812 nucleus->setZ(theEventInfo.Zt); 484 nucleus->setZ(theEventInfo.Zt);
813 485
814 // CN kinematical variables 486 // CN kinematical variables
815 // Note: the CN orbital angular momentum i 487 // Note: the CN orbital angular momentum is neglected in what follows. We
816 // should actually take it into account! 488 // should actually take it into account!
817 ThreeVector theCNMomentum = nucleus->getIn 489 ThreeVector theCNMomentum = nucleus->getIncomingMomentum();
818 ThreeVector theCNSpin = nucleus->getIncomi 490 ThreeVector theCNSpin = nucleus->getIncomingAngularMomentum();
819 const G4double theTargetMass = ParticleTab << 491 const G4double theTargetMass = ParticleTable::getTableMass(theEventInfo.At, theEventInfo.Zt);
820 G4int theCNA=theEventInfo.At, theCNZ=theEv << 492 G4int theCNA=theEventInfo.At, theCNZ=theEventInfo.Zt;
821 Cluster * const theProjectileRemnant = nuc 493 Cluster * const theProjectileRemnant = nucleus->getProjectileRemnant();
822 G4double theCNEnergy = theTargetMass + the 494 G4double theCNEnergy = theTargetMass + theProjectileRemnant->getEnergy();
823 495
824 // Loop over the potential participants 496 // Loop over the potential participants
825 ParticleList const &initialProjectileCompo 497 ParticleList const &initialProjectileComponents = theProjectileRemnant->getParticles();
826 std::vector<Particle *> shuffledComponents 498 std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
827 // Shuffle the list of potential participa 499 // Shuffle the list of potential participants
828 std::shuffle(shuffledComponents.begin(), s << 500 std::random_shuffle(shuffledComponents.begin(), shuffledComponents.end(), Random::getAdapter());
829 501
830 G4bool success = true; 502 G4bool success = true;
831 G4bool atLeastOneNucleonEntering = false; 503 G4bool atLeastOneNucleonEntering = false;
832 for(std::vector<Particle*>::const_iterator 504 for(std::vector<Particle*>::const_iterator p=shuffledComponents.begin(), e=shuffledComponents.end(); p!=e; ++p) {
833 // Skip particles that miss the interact 505 // Skip particles that miss the interaction distance
834 Intersection intersectionInteractionDist 506 Intersection intersectionInteractionDistance(IntersectionFactory::getEarlierTrajectoryIntersection(
835 (*p)->getPosition(), 507 (*p)->getPosition(),
836 (*p)->getPropagationVelocity(), 508 (*p)->getPropagationVelocity(),
837 maxInteractionDistance)); 509 maxInteractionDistance));
838 if(!intersectionInteractionDistance.exis 510 if(!intersectionInteractionDistance.exists)
839 continue; 511 continue;
840 512
841 // Build an entry avatar for this nucleo 513 // Build an entry avatar for this nucleon
842 atLeastOneNucleonEntering = true; 514 atLeastOneNucleonEntering = true;
843 ParticleEntryAvatar *theAvatar = new Par 515 ParticleEntryAvatar *theAvatar = new ParticleEntryAvatar(0.0, nucleus, *p);
844 nucleus->getStore()->addParticleEntryAva 516 nucleus->getStore()->addParticleEntryAvatar(theAvatar);
845 FinalState *fs = theAvatar->getFinalStat 517 FinalState *fs = theAvatar->getFinalState();
846 nucleus->applyFinalState(fs); 518 nucleus->applyFinalState(fs);
847 FinalStateValidity validity = fs->getVal 519 FinalStateValidity validity = fs->getValidity();
848 delete fs; 520 delete fs;
849 switch(validity) { 521 switch(validity) {
850 case ValidFS: 522 case ValidFS:
851 case ParticleBelowFermiFS: 523 case ParticleBelowFermiFS:
852 case ParticleBelowZeroFS: 524 case ParticleBelowZeroFS:
853 // Add the particle to the CN 525 // Add the particle to the CN
854 theCNA++; 526 theCNA++;
855 theCNZ += (*p)->getZ(); 527 theCNZ += (*p)->getZ();
856 theCNS += (*p)->getS(); <<
857 break; 528 break;
858 case PauliBlockedFS: 529 case PauliBlockedFS:
859 case NoEnergyConservationFS: 530 case NoEnergyConservationFS:
860 default: 531 default:
861 success = false; 532 success = false;
862 break; 533 break;
863 } 534 }
864 } 535 }
865 536
866 if(!success || !atLeastOneNucleonEntering) 537 if(!success || !atLeastOneNucleonEntering) {
867 INCL_DEBUG("No nucleon entering in force 538 INCL_DEBUG("No nucleon entering in forced CN, forcing a transparent" << '\n');
868 forceTransparent = true; 539 forceTransparent = true;
869 return; 540 return;
870 } 541 }
871 542
872 // assert(theCNA==nucleus->getA()); 543 // assert(theCNA==nucleus->getA());
873 // assert(theCNA<=theEventInfo.At+theEventInfo 544 // assert(theCNA<=theEventInfo.At+theEventInfo.Ap);
874 // assert(theCNZ<=theEventInfo.Zt+theEventInfo 545 // assert(theCNZ<=theEventInfo.Zt+theEventInfo.Zp);
875 // assert(theCNS>=theEventInfo.St+theEventInfo <<
876 546
877 // Update the kinematics of the CN 547 // Update the kinematics of the CN
878 theCNEnergy -= theProjectileRemnant->getEn 548 theCNEnergy -= theProjectileRemnant->getEnergy();
879 theCNMomentum -= theProjectileRemnant->get 549 theCNMomentum -= theProjectileRemnant->getMomentum();
880 550
881 // Deal with the projectile remnant 551 // Deal with the projectile remnant
882 nucleus->finalizeProjectileRemnant(propaga 552 nucleus->finalizeProjectileRemnant(propagationModel->getCurrentTime());
883 553
884 // Subtract the angular momentum of the pr 554 // Subtract the angular momentum of the projectile remnant
885 // assert(nucleus->getStore()->getOutgoingPart 555 // assert(nucleus->getStore()->getOutgoingParticles().empty());
886 theCNSpin -= theProjectileRemnant->getAngu 556 theCNSpin -= theProjectileRemnant->getAngularMomentum();
887 557
888 // Compute the excitation energy of the CN 558 // Compute the excitation energy of the CN
889 const G4double theCNMass = ParticleTable:: << 559 const G4double theCNMass = ParticleTable::getTableMass(theCNA,theCNZ);
890 const G4double theCNInvariantMassSquared = 560 const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.mag2();
891 if(theCNInvariantMassSquared<0.) { 561 if(theCNInvariantMassSquared<0.) {
892 // Negative invariant mass squared, retu 562 // Negative invariant mass squared, return a transparent
893 forceTransparent = true; 563 forceTransparent = true;
894 return; 564 return;
895 } 565 }
896 const G4double theCNExcitationEnergy = std 566 const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
897 if(theCNExcitationEnergy<0.) { 567 if(theCNExcitationEnergy<0.) {
898 // Negative excitation energy, return a 568 // Negative excitation energy, return a transparent
899 INCL_DEBUG("CN excitation energy is nega 569 INCL_DEBUG("CN excitation energy is negative, forcing a transparent" << '\n'
900 << " theCNA = " << theCNA << '\n' 570 << " theCNA = " << theCNA << '\n'
901 << " theCNZ = " << theCNZ << '\n' 571 << " theCNZ = " << theCNZ << '\n'
902 << " theCNS = " << theCNS << '\n' <<
903 << " theCNEnergy = " << theCNEner 572 << " theCNEnergy = " << theCNEnergy << '\n'
904 << " theCNMomentum = (" << theCNM 573 << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
905 << " theCNExcitationEnergy = " << 574 << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
906 << " theCNSpin = (" << theCNSpin. 575 << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
907 ); 576 );
908 forceTransparent = true; 577 forceTransparent = true;
909 return; 578 return;
910 } else { 579 } else {
911 // Positive excitation energy, can make 580 // Positive excitation energy, can make a CN
912 INCL_DEBUG("CN excitation energy is posi 581 INCL_DEBUG("CN excitation energy is positive, forcing a CN" << '\n'
913 << " theCNA = " << theCNA << '\n' 582 << " theCNA = " << theCNA << '\n'
914 << " theCNZ = " << theCNZ << '\n' 583 << " theCNZ = " << theCNZ << '\n'
915 << " theCNS = " << theCNS << '\n' <<
916 << " theCNEnergy = " << theCNEner 584 << " theCNEnergy = " << theCNEnergy << '\n'
917 << " theCNMomentum = (" << theCNM 585 << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
918 << " theCNExcitationEnergy = " << 586 << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
919 << " theCNSpin = (" << theCNSpin. 587 << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
920 ); 588 );
921 nucleus->setA(theCNA); 589 nucleus->setA(theCNA);
922 nucleus->setZ(theCNZ); 590 nucleus->setZ(theCNZ);
923 nucleus->setS(theCNS); <<
924 nucleus->setMomentum(theCNMomentum); 591 nucleus->setMomentum(theCNMomentum);
925 nucleus->setEnergy(theCNEnergy); 592 nucleus->setEnergy(theCNEnergy);
926 nucleus->setExcitationEnergy(theCNExcita 593 nucleus->setExcitationEnergy(theCNExcitationEnergy);
927 nucleus->setMass(theCNMass+theCNExcitati 594 nucleus->setMass(theCNMass+theCNExcitationEnergy);
928 nucleus->setSpin(theCNSpin); // neglects 595 nucleus->setSpin(theCNSpin); // neglects any orbital angular momentum of the CN
929 596
930 // Take care of any remaining deltas 597 // Take care of any remaining deltas
931 theEventInfo.forcedDeltasOutside = nucle 598 theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
932 599
933 // Take care of any remaining etas and/o << 600 // Take care of any remaining etas and/or omegas
934 G4double timeThreshold=theConfig->getDec << 601 G4double timeThreshold=theConfig->getDecayTimeThreshold();
935 theEventInfo.forcedPionResonancesOutside 602 theEventInfo.forcedPionResonancesOutside = nucleus->decayOutgoingPionResonances(timeThreshold);
936 << 603
937 // Take care of any remaining Kaons <<
938 theEventInfo.emitKaon = nucleus->emitIns <<
939 <<
940 // Cluster decay 604 // Cluster decay
941 theEventInfo.clusterDecay = nucleus->dec << 605 theEventInfo.clusterDecay = nucleus->decayOutgoingClusters() | nucleus->decayMe();
942 606
943 // Fill the EventInfo structure 607 // Fill the EventInfo structure
944 nucleus->fillEventInfo(&theEventInfo); 608 nucleus->fillEventInfo(&theEventInfo);
945 } 609 }
946 } 610 }
947 611
948 void INCL::rescaleOutgoingForRecoil() { 612 void INCL::rescaleOutgoingForRecoil() {
949 RecoilCMFunctor theRecoilFunctor(nucleus, 613 RecoilCMFunctor theRecoilFunctor(nucleus, theEventInfo);
950 614
951 // Apply the root-finding algorithm 615 // Apply the root-finding algorithm
952 const RootFinder::Solution theSolution = R 616 const RootFinder::Solution theSolution = RootFinder::solve(&theRecoilFunctor, 1.0);
953 if(theSolution.success) { 617 if(theSolution.success) {
954 theRecoilFunctor(theSolution.x); // Appl 618 theRecoilFunctor(theSolution.x); // Apply the solution
955 } else { 619 } else {
956 INCL_WARN("Couldn't accommodate remnant 620 INCL_WARN("Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." << '\n');
957 } 621 }
>> 622
958 } 623 }
959 624
960 #ifndef INCLXX_IN_GEANT4_MODE 625 #ifndef INCLXX_IN_GEANT4_MODE
961 void INCL::globalConservationChecks(G4bool a 626 void INCL::globalConservationChecks(G4bool afterRecoil) {
962 Nucleus::ConservationBalance theBalance = 627 Nucleus::ConservationBalance theBalance = nucleus->getConservationBalance(theEventInfo,afterRecoil);
963 628
964 // Global conservation checks 629 // Global conservation checks
965 const G4double pLongBalance = theBalance.m 630 const G4double pLongBalance = theBalance.momentum.getZ();
966 const G4double pTransBalance = theBalance. 631 const G4double pTransBalance = theBalance.momentum.perp();
967 if(theBalance.Z != 0) { 632 if(theBalance.Z != 0) {
968 INCL_ERROR("Violation of charge conserva << 633 INCL_ERROR("Violation of charge conservation! ZBalance = " << theBalance.Z << '\n');
969 } 634 }
970 if(theBalance.A != 0) { 635 if(theBalance.A != 0) {
971 INCL_ERROR("Violation of baryon-number c << 636 INCL_ERROR("Violation of baryon-number conservation! ABalance = " << theBalance.A << '\n');
972 } <<
973 if(theBalance.S != 0) { <<
974 INCL_ERROR("Violation of strange-number <<
975 } 637 }
976 G4double EThreshold, pLongThreshold, pTran 638 G4double EThreshold, pLongThreshold, pTransThreshold;
977 if(afterRecoil) { 639 if(afterRecoil) {
978 // Less stringent checks after accommoda 640 // Less stringent checks after accommodating recoil
979 EThreshold = 10.; // MeV 641 EThreshold = 10.; // MeV
980 pLongThreshold = 1.; // MeV/c 642 pLongThreshold = 1.; // MeV/c
981 pTransThreshold = 1.; // MeV/c 643 pTransThreshold = 1.; // MeV/c
982 } else { 644 } else {
983 // More stringent checks before accommod 645 // More stringent checks before accommodating recoil
984 EThreshold = 0.1; // MeV 646 EThreshold = 0.1; // MeV
985 pLongThreshold = 0.1; // MeV/c 647 pLongThreshold = 0.1; // MeV/c
986 pTransThreshold = 0.1; // MeV/c 648 pTransThreshold = 0.1; // MeV/c
987 } 649 }
988 if(std::abs(theBalance.energy)>EThreshold) 650 if(std::abs(theBalance.energy)>EThreshold) {
989 INCL_WARN("Violation of energy conservat << 651 INCL_WARN("Violation of energy conservation > " << EThreshold << " MeV. EBalance = " << theBalance.energy << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
990 } 652 }
991 if(std::abs(pLongBalance)>pLongThreshold) 653 if(std::abs(pLongBalance)>pLongThreshold) {
992 INCL_WARN("Violation of longitudinal mom 654 INCL_WARN("Violation of longitudinal momentum conservation > " << pLongThreshold << " MeV/c. pLongBalance = " << pLongBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
993 } 655 }
994 if(std::abs(pTransBalance)>pTransThreshold 656 if(std::abs(pTransBalance)>pTransThreshold) {
995 INCL_WARN("Violation of transverse momen 657 INCL_WARN("Violation of transverse momentum conservation > " << pTransThreshold << " MeV/c. pTransBalance = " << pTransBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
996 } 658 }
997 659
998 // Feed the EventInfo variables 660 // Feed the EventInfo variables
999 theEventInfo.EBalance = theBalance.energy; 661 theEventInfo.EBalance = theBalance.energy;
1000 theEventInfo.pLongBalance = pLongBalance; 662 theEventInfo.pLongBalance = pLongBalance;
1001 theEventInfo.pTransBalance = pTransBalanc 663 theEventInfo.pTransBalance = pTransBalance;
1002 } 664 }
1003 #endif 665 #endif
1004 666
1005 G4bool INCL::continueCascade() { 667 G4bool INCL::continueCascade() {
1006 // Stop if we have passed the stopping ti 668 // Stop if we have passed the stopping time
1007 if(propagationModel->getCurrentTime() > p 669 if(propagationModel->getCurrentTime() > propagationModel->getStoppingTime()) {
1008 INCL_DEBUG("Cascade time (" << propagat 670 INCL_DEBUG("Cascade time (" << propagationModel->getCurrentTime()
1009 << ") exceeded stopping time (" << 671 << ") exceeded stopping time (" << propagationModel->getStoppingTime()
1010 << "), stopping cascade" << '\n'); 672 << "), stopping cascade" << '\n');
1011 return false; 673 return false;
1012 } 674 }
1013 // Stop if there are no participants and 675 // Stop if there are no participants and no pions inside the nucleus
1014 if(nucleus->getStore()->getBook().getCasc 676 if(nucleus->getStore()->getBook().getCascading()==0 &&
1015 nucleus->getStore()->getIncomingParti 677 nucleus->getStore()->getIncomingParticles().empty()) {
1016 INCL_DEBUG("No participants in the nucl 678 INCL_DEBUG("No participants in the nucleus and no incoming particles left, stopping cascade" << '\n');
1017 return false; 679 return false;
1018 } 680 }
1019 // Stop if the remnant is smaller than mi 681 // Stop if the remnant is smaller than minRemnantSize
1020 if(nucleus->getA() <= minRemnantSize) { 682 if(nucleus->getA() <= minRemnantSize) {
1021 INCL_DEBUG("Remnant size (" << nucleus- 683 INCL_DEBUG("Remnant size (" << nucleus->getA()
1022 << ") smaller than or equal to mini 684 << ") smaller than or equal to minimum (" << minRemnantSize
1023 << "), stopping cascade" << '\n'); 685 << "), stopping cascade" << '\n');
1024 return false; 686 return false;
1025 } 687 }
1026 // Stop if we have to try and make a comp 688 // Stop if we have to try and make a compound nucleus or if we have to
1027 // force a transparent 689 // force a transparent
1028 if(nucleus->getTryCompoundNucleus()) { 690 if(nucleus->getTryCompoundNucleus()) {
1029 INCL_DEBUG("Trying to make a compound n 691 INCL_DEBUG("Trying to make a compound nucleus, stopping cascade" << '\n');
1030 return false; 692 return false;
1031 } 693 }
1032 694
1033 return true; 695 return true;
1034 } 696 }
1035 697
1036 void INCL::finalizeGlobalInfo(Random::SeedV 698 void INCL::finalizeGlobalInfo(Random::SeedVector const &initialSeeds) {
1037 const G4double normalisationFactor = theG 699 const G4double normalisationFactor = theGlobalInfo.geometricCrossSection /
1038 ((G4double) theGlobalInfo.nShots); 700 ((G4double) theGlobalInfo.nShots);
1039 theGlobalInfo.nucleonAbsorptionCrossSecti 701 theGlobalInfo.nucleonAbsorptionCrossSection = normalisationFactor *
1040 ((G4double) theGlobalInfo.nNucleonAbsor 702 ((G4double) theGlobalInfo.nNucleonAbsorptions);
1041 theGlobalInfo.pionAbsorptionCrossSection 703 theGlobalInfo.pionAbsorptionCrossSection = normalisationFactor *
1042 ((G4double) theGlobalInfo.nPionAbsorpti 704 ((G4double) theGlobalInfo.nPionAbsorptions);
1043 theGlobalInfo.reactionCrossSection = norm 705 theGlobalInfo.reactionCrossSection = normalisationFactor *
1044 ((G4double) (theGlobalInfo.nShots - the 706 ((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents));
1045 theGlobalInfo.errorReactionCrossSection = 707 theGlobalInfo.errorReactionCrossSection = normalisationFactor *
1046 std::sqrt((G4double) (theGlobalInfo.nSh 708 std::sqrt((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents));
1047 theGlobalInfo.forcedCNCrossSection = norm 709 theGlobalInfo.forcedCNCrossSection = normalisationFactor *
1048 ((G4double) theGlobalInfo.nForcedCompou 710 ((G4double) theGlobalInfo.nForcedCompoundNucleus);
1049 theGlobalInfo.errorForcedCNCrossSection = 711 theGlobalInfo.errorForcedCNCrossSection = normalisationFactor *
1050 std::sqrt((G4double) (theGlobalInfo.nFo 712 std::sqrt((G4double) (theGlobalInfo.nForcedCompoundNucleus));
1051 theGlobalInfo.completeFusionCrossSection 713 theGlobalInfo.completeFusionCrossSection = normalisationFactor *
1052 ((G4double) theGlobalInfo.nCompleteFusi 714 ((G4double) theGlobalInfo.nCompleteFusion);
1053 theGlobalInfo.errorCompleteFusionCrossSec 715 theGlobalInfo.errorCompleteFusionCrossSection = normalisationFactor *
1054 std::sqrt((G4double) (theGlobalInfo.nCo 716 std::sqrt((G4double) (theGlobalInfo.nCompleteFusion));
1055 theGlobalInfo.energyViolationInteractionC 717 theGlobalInfo.energyViolationInteractionCrossSection = normalisationFactor *
1056 ((G4double) theGlobalInfo.nEnergyViolat 718 ((G4double) theGlobalInfo.nEnergyViolationInteraction);
1057 719
1058 theGlobalInfo.initialRandomSeeds.assign(i 720 theGlobalInfo.initialRandomSeeds.assign(initialSeeds.begin(), initialSeeds.end());
1059 721
1060 Random::SeedVector theSeeds = Random::get 722 Random::SeedVector theSeeds = Random::getSeeds();
1061 theGlobalInfo.finalRandomSeeds.assign(the 723 theGlobalInfo.finalRandomSeeds.assign(theSeeds.begin(), theSeeds.end());
1062 } 724 }
1063 725
1064 G4int INCL::makeProjectileRemnant() { 726 G4int INCL::makeProjectileRemnant() {
1065 // Do nothing if this is not a nucleus-nu 727 // Do nothing if this is not a nucleus-nucleus reaction
1066 if(!nucleus->getProjectileRemnant()) 728 if(!nucleus->getProjectileRemnant())
1067 return 0; 729 return 0;
1068 730
1069 // Get the spectators (geometrical+dynami 731 // Get the spectators (geometrical+dynamical) from the Store
1070 ParticleList geomSpectators(nucleus->getP 732 ParticleList geomSpectators(nucleus->getProjectileRemnant()->getParticles());
1071 ParticleList dynSpectators(nucleus->getSt 733 ParticleList dynSpectators(nucleus->getStore()->extractDynamicalSpectators());
1072 734
1073 G4int nUnmergedSpectators = 0; 735 G4int nUnmergedSpectators = 0;
1074 736
1075 // If there are no spectators, do nothing 737 // If there are no spectators, do nothing
1076 if(dynSpectators.empty() && geomSpectator 738 if(dynSpectators.empty() && geomSpectators.empty()) {
1077 return 0; 739 return 0;
1078 } else if(dynSpectators.size()==1 && geom 740 } else if(dynSpectators.size()==1 && geomSpectators.empty()) {
1079 // No geometrical spectators, one dynam 741 // No geometrical spectators, one dynamical spectator
1080 // Just put it back in the outgoing lis 742 // Just put it back in the outgoing list
1081 nucleus->getStore()->addToOutgoing(dynS 743 nucleus->getStore()->addToOutgoing(dynSpectators.front());
1082 } else { 744 } else {
1083 // Make a cluster out of the geometrica 745 // Make a cluster out of the geometrical spectators
1084 ProjectileRemnant *theProjectileRemnant 746 ProjectileRemnant *theProjectileRemnant = nucleus->getProjectileRemnant();
1085 747
1086 // Add the dynamical spectators to the 748 // Add the dynamical spectators to the bunch
1087 ParticleList rejected = theProjectileRe 749 ParticleList rejected = theProjectileRemnant->addAllDynamicalSpectators(dynSpectators);
1088 // Put back the rejected spectators int 750 // Put back the rejected spectators into the outgoing list
1089 nUnmergedSpectators = (G4int)rejected.s << 751 nUnmergedSpectators = rejected.size();
1090 nucleus->getStore()->addToOutgoing(reje 752 nucleus->getStore()->addToOutgoing(rejected);
1091 753
1092 // Deal with the projectile remnant 754 // Deal with the projectile remnant
1093 nucleus->finalizeProjectileRemnant(prop 755 nucleus->finalizeProjectileRemnant(propagationModel->getCurrentTime());
1094 756
1095 } 757 }
1096 758
1097 return nUnmergedSpectators; 759 return nUnmergedSpectators;
1098 } 760 }
1099 761
1100 void INCL::initMaxInteractionDistance(Parti 762 void INCL::initMaxInteractionDistance(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
1101 if(projectileSpecies.theType != Composite 763 if(projectileSpecies.theType != Composite) {
1102 maxInteractionDistance = 0.; 764 maxInteractionDistance = 0.;
1103 return; 765 return;
1104 } 766 }
1105 767
1106 const G4double r0 = std::max(ParticleTabl 768 const G4double r0 = std::max(ParticleTable::getNuclearRadius(Proton, theA, theZ),
1107 ParticleTable: 769 ParticleTable::getNuclearRadius(Neutron, theA, theZ));
1108 770
1109 const G4double theNNDistance = CrossSecti 771 const G4double theNNDistance = CrossSections::interactionDistanceNN(projectileSpecies, kineticEnergy);
1110 maxInteractionDistance = r0 + theNNDistan 772 maxInteractionDistance = r0 + theNNDistance;
1111 INCL_DEBUG("Initialised interaction dista 773 INCL_DEBUG("Initialised interaction distance: r0 = " << r0 << '\n'
1112 << " theNNDistance = " << theNND 774 << " theNNDistance = " << theNNDistance << '\n'
1113 << " maxInteractionDistance = " 775 << " maxInteractionDistance = " << maxInteractionDistance << '\n');
1114 } 776 }
1115 777
1116 void INCL::initUniverseRadius(ParticleSpeci 778 void INCL::initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z) {
1117 G4double rMax = 0.0; 779 G4double rMax = 0.0;
1118 if(A==0) { 780 if(A==0) {
1119 IsotopicDistribution const &anIsotopicD 781 IsotopicDistribution const &anIsotopicDistribution =
1120 ParticleTable::getNaturalIsotopicDist 782 ParticleTable::getNaturalIsotopicDistribution(Z);
1121 IsotopeVector theIsotopes = anIsotopicD 783 IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
1122 for(IsotopeIter i=theIsotopes.begin(), 784 for(IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
1123 const G4double pMaximumRadius = Parti 785 const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, i->theA, Z);
1124 const G4double nMaximumRadius = Parti 786 const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, i->theA, Z);
1125 const G4double maximumRadius = std::m << 787 const G4double maximumRadius = std::min(pMaximumRadius, nMaximumRadius);
1126 rMax = std::max(maximumRadius, rMax); 788 rMax = std::max(maximumRadius, rMax);
1127 } 789 }
1128 } else { 790 } else {
1129 const G4double pMaximumRadius = Particl 791 const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
1130 const G4double nMaximumRadius = Particl 792 const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
1131 const G4double maximumRadius = std::max << 793 rMax = std::min(pMaximumRadius, nMaximumRadius);
1132 rMax = std::max(maximumRadius, rMax); <<
1133 } 794 }
1134 if(p.theType==Composite || p.theType==Pro 795 if(p.theType==Composite || p.theType==Proton || p.theType==Neutron) {
1135 const G4double interactionDistanceNN = 796 const G4double interactionDistanceNN = CrossSections::interactionDistanceNN(p, kineticEnergy);
1136 maxUniverseRadius = rMax + interactionD 797 maxUniverseRadius = rMax + interactionDistanceNN;
1137 } else if(p.theType==PiPlus 798 } else if(p.theType==PiPlus
1138 || p.theType==PiZero 799 || p.theType==PiZero
1139 || p.theType==PiMinus) { 800 || p.theType==PiMinus) {
1140 const G4double interactionDistancePiN = 801 const G4double interactionDistancePiN = CrossSections::interactionDistancePiN(kineticEnergy);
1141 maxUniverseRadius = rMax + interactionD 802 maxUniverseRadius = rMax + interactionDistancePiN;
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 } 803 }
1160 INCL_DEBUG("Initialised universe radius: 804 INCL_DEBUG("Initialised universe radius: " << maxUniverseRadius << '\n');
1161 } 805 }
1162 806
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() { 807 void INCL::updateGlobalInfo() {
1187 // Increment the global counter for the n 808 // Increment the global counter for the number of shots
1188 theGlobalInfo.nShots++; 809 theGlobalInfo.nShots++;
1189 810
1190 if(theEventInfo.transparent) { 811 if(theEventInfo.transparent) {
1191 // Increment the global counter for the 812 // Increment the global counter for the number of transparents
1192 theGlobalInfo.nTransparents++; 813 theGlobalInfo.nTransparents++;
1193 // Increment the global counter for the 814 // Increment the global counter for the number of forced transparents
1194 if(forceTransparent) 815 if(forceTransparent)
1195 theGlobalInfo.nForcedTransparents++; 816 theGlobalInfo.nForcedTransparents++;
1196 return; 817 return;
1197 } 818 }
1198 819
1199 // Check if we have an absorption: 820 // Check if we have an absorption:
1200 if(theEventInfo.nucleonAbsorption) theGlo 821 if(theEventInfo.nucleonAbsorption) theGlobalInfo.nNucleonAbsorptions++;
1201 if(theEventInfo.pionAbsorption) theGlobal 822 if(theEventInfo.pionAbsorption) theGlobalInfo.nPionAbsorptions++;
1202 823
1203 // Count complete-fusion events 824 // Count complete-fusion events
1204 if(theEventInfo.nCascadeParticles==0) the 825 if(theEventInfo.nCascadeParticles==0) theGlobalInfo.nCompleteFusion++;
1205 826
1206 if(nucleus->getTryCompoundNucleus()) 827 if(nucleus->getTryCompoundNucleus())
1207 theGlobalInfo.nForcedCompoundNucleus++; 828 theGlobalInfo.nForcedCompoundNucleus++;
1208 829
1209 // Counters for the number of violations 830 // Counters for the number of violations of energy conservation in
1210 // collisions 831 // collisions
1211 theGlobalInfo.nEnergyViolationInteraction 832 theGlobalInfo.nEnergyViolationInteraction += theEventInfo.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 } 833 }
1390 834
1391 } 835 }
1392 836