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