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