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