| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA
29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA
31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege
32 // Davide Mancusi, CEA-Saclay, France <<
33 // 32 //
34 #define INCLXX_IN_GEANT4_MODE 1 33 #define INCLXX_IN_GEANT4_MODE 1
35 34
36 #include "globals.hh" 35 #include "globals.hh"
37 36
38 /* 37 /*
39 * StandardPropagationModel.cpp 38 * StandardPropagationModel.cpp
40 * 39 *
41 * \date 4 juin 2009 40 * \date 4 juin 2009
42 * \author Pekka Kaitaniemi 41 * \author Pekka Kaitaniemi
43 */ 42 */
44 43
45 #include "G4INCLStandardPropagationModel.hh" 44 #include "G4INCLStandardPropagationModel.hh"
46 #include "G4INCLPbarAtrestEntryChannel.hh" <<
47 #include "G4INCLSurfaceAvatar.hh" 45 #include "G4INCLSurfaceAvatar.hh"
48 #include "G4INCLBinaryCollisionAvatar.hh" 46 #include "G4INCLBinaryCollisionAvatar.hh"
49 #include "G4INCLDecayAvatar.hh" 47 #include "G4INCLDecayAvatar.hh"
50 #include "G4INCLCrossSections.hh" 48 #include "G4INCLCrossSections.hh"
51 #include "G4INCLRandom.hh" 49 #include "G4INCLRandom.hh"
52 #include <iostream> 50 #include <iostream>
53 #include <algorithm> 51 #include <algorithm>
54 #include "G4INCLLogger.hh" 52 #include "G4INCLLogger.hh"
55 #include "G4INCLGlobals.hh" 53 #include "G4INCLGlobals.hh"
56 #include "G4INCLKinematicsUtils.hh" 54 #include "G4INCLKinematicsUtils.hh"
57 #include "G4INCLCoulombDistortion.hh" 55 #include "G4INCLCoulombDistortion.hh"
58 #include "G4INCLDeltaDecayChannel.hh" 56 #include "G4INCLDeltaDecayChannel.hh"
59 #include "G4INCLSigmaZeroDecayChannel.hh" <<
60 #include "G4INCLPionResonanceDecayChannel.hh" <<
61 #include "G4INCLParticleEntryAvatar.hh" 57 #include "G4INCLParticleEntryAvatar.hh"
62 #include "G4INCLIntersection.hh" 58 #include "G4INCLIntersection.hh"
63 #include <vector> <<
64 59
65 namespace G4INCL { 60 namespace G4INCL {
66 61
67 StandardPropagationModel::StandardPropagat << 62 StandardPropagationModel::StandardPropagationModel(LocalEnergyType localEnergyType, LocalEnergyType localEnergyDeltaType)
68 :theNucleus(0), maximumTime(70.0), curre << 63 :theNucleus(0), maximumTime(70.0), currentTime(0.0), firstAvatar(true),
69 hadronizationTime(hTime), <<
70 firstAvatar(true), <<
71 theLocalEnergyType(localEnergyType), 64 theLocalEnergyType(localEnergyType),
72 theLocalEnergyDeltaType(localEnergyDelta 65 theLocalEnergyDeltaType(localEnergyDeltaType)
73 { 66 {
74 } 67 }
75 68
76 StandardPropagationModel::~StandardPropaga 69 StandardPropagationModel::~StandardPropagationModel()
77 { 70 {
78 delete theNucleus; 71 delete theNucleus;
79 } 72 }
80 73
81 G4INCL::Nucleus* StandardPropagationModel: 74 G4INCL::Nucleus* StandardPropagationModel::getNucleus()
82 { 75 {
83 return theNucleus; 76 return theNucleus;
84 } 77 }
85 78
86 //D << 79 G4double StandardPropagationModel::shoot(ParticleSpecies const projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
87 << 80 if(projectileSpecies.theType==Composite)
88 G4double StandardPropagationModel::shoot(P <<
89 if(projectileSpecies.theType==Composite) <<
90 return shootComposite(projectileSpecie 81 return shootComposite(projectileSpecies, kineticEnergy, impactParameter, phi);
91 } << 82 else
92 else if(projectileSpecies.theType==antiP <<
93 return shootAtrest(projectileSpecies.t <<
94 } <<
95 else{ <<
96 return shootParticle(projectileSpecies 83 return shootParticle(projectileSpecies.theType, kineticEnergy, impactParameter, phi);
97 } <<
98 } 84 }
99 85
100 //D <<
101 <<
102 G4double StandardPropagationModel::shootAt <<
103 theNucleus->setParticleNucleusCollision( <<
104 currentTime = 0.0; <<
105 <<
106 // Create final state particles <<
107 const G4double projectileMass = Particle <<
108 G4double energy = kineticEnergy + projec <<
109 G4double momentumZ = std::sqrt(energy*en <<
110 ThreeVector momentum(0.0, 0.0, momentumZ <<
111 Particle *pb = new G4INCL::Particle(t, e <<
112 PbarAtrestEntryChannel *obj = new PbarAt <<
113 ParticleList fslist = obj->makeMesonStar <<
114 const G4bool isProton = obj->ProtonIsThe <<
115 delete pb; <<
116 <<
117 //set Stopping time according to highest <<
118 G4double temfin; <<
119 G4double TLab; <<
120 std::vector<G4double> energies; <<
121 std::vector<G4double> projections; <<
122 ThreeVector ab, cd; <<
123 <<
124 for(ParticleIter pit = fslist.begin(), e <<
125 energies.push_back((*pit)->getKineticE <<
126 ab = (*pit)->boostVector(); <<
127 cd = (*pit)->getPosition(); <<
128 projections.push_back(ab.dot(cd)); //p <<
129 }// make vector of energies <<
130 <<
131 temfin = 30.18 * std::pow(theNucleus->ge <<
132 TLab = *max_element(energies.begin(), en <<
133 <<
134 // energy-dependent stopping time above <<
135 if(TLab>2000.) <<
136 temfin *= (5.8E4-TLab)/5.6E4; <<
137 <<
138 maximumTime = temfin; <<
139 <<
140 // If the incoming particle is slow, use <<
141 const G4double rMax = theNucleus->getUni <<
142 const G4double distance = 2.*rMax; <<
143 const G4double maxMesonVelocityProjectio <<
144 const G4double traversalTime = distance <<
145 if(maximumTime < traversalTime) <<
146 maximumTime = traversalTime; <<
147 INCL_DEBUG("Cascade stopping time is " < <<
148 <<
149 <<
150 // Fill in the relevant kinematic variab <<
151 theNucleus->setIncomingAngularMomentum(G <<
152 theNucleus->setIncomingMomentum(G4INCL:: <<
153 if(isProton){ <<
154 theNucleus->setInitialEnergy(pb->getMa <<
155 + ParticleTable::getTableMass(theNuc <<
156 } <<
157 else{ <<
158 theNucleus->setInitialEnergy(pb->getMa <<
159 + ParticleTable::getTableMass(theNuc <<
160 } <<
161 //kinetic energy excluded from the balan <<
162 <<
163 for(ParticleIter p = fslist.begin(), e = <<
164 (*p)->makeProjectileSpectator(); <<
165 } <<
166 <<
167 generateAllAvatars(); <<
168 firstAvatar = false; <<
169 <<
170 // Get the entry avatars for mesons <<
171 IAvatarList theAvatarList = obj->bringMe <<
172 delete obj; <<
173 theNucleus->getStore()->addParticleEntry <<
174 INCL_DEBUG("Avatars added" << '\n'); <<
175 <<
176 return 99.; <<
177 } <<
178 <<
179 //D <<
180 <<
181 G4double StandardPropagationModel::shootPa 86 G4double StandardPropagationModel::shootParticle(ParticleType const type, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
182 theNucleus->setParticleNucleusCollision( 87 theNucleus->setParticleNucleusCollision();
183 currentTime = 0.0; 88 currentTime = 0.0;
184 89
185 // Create the projectile particle 90 // Create the projectile particle
186 const G4double projectileMass = Particle 91 const G4double projectileMass = ParticleTable::getTableParticleMass(type);
187 G4double energy = kineticEnergy + projec 92 G4double energy = kineticEnergy + projectileMass;
188 G4double momentumZ = std::sqrt(energy*en 93 G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
189 ThreeVector momentum(0.0, 0.0, momentumZ 94 ThreeVector momentum(0.0, 0.0, momentumZ);
190 Particle *p= new G4INCL::Particle(type, 95 Particle *p= new G4INCL::Particle(type, energy, momentum, ThreeVector());
191 96
192 G4double temfin; << 97 G4double temfin = 0.0;
193 G4double TLab; << 98 if( p->isNucleon() )
194 if( p->isMeson()) { <<
195 temfin = 30.18 * std::pow(theNucleus-> <<
196 TLab = p->getKineticEnergy(); <<
197 } else { <<
198 temfin = 29.8 * std::pow(theNucleus->g 99 temfin = 29.8 * std::pow(theNucleus->getA(), 0.16);
199 TLab = p->getKineticEnergy()/p->getA() << 100 else {
>> 101 const G4double tlab = p->getEnergy() - p->getMass();
>> 102 temfin = 30.18 * std::pow(theNucleus->getA(), 0.17*(1.0 - 5.7E-5*tlab));
200 } 103 }
201 104
202 // energy-dependent stopping time above <<
203 if(TLab>2000.) <<
204 temfin *= (5.8E4-TLab)/5.6E4; <<
205 <<
206 maximumTime = temfin; 105 maximumTime = temfin;
207 106
208 // If the incoming particle is slow, use 107 // If the incoming particle is slow, use a larger stopping time
209 const G4double rMax = theNucleus->getUni 108 const G4double rMax = theNucleus->getUniverseRadius();
210 const G4double distance = 2.*rMax; 109 const G4double distance = 2.*rMax;
211 const G4double projectileVelocity = p->b 110 const G4double projectileVelocity = p->boostVector().mag();
212 const G4double traversalTime = distance 111 const G4double traversalTime = distance / projectileVelocity;
213 if(maximumTime < traversalTime) 112 if(maximumTime < traversalTime)
214 maximumTime = traversalTime; 113 maximumTime = traversalTime;
215 INCL_DEBUG("Cascade stopping time is " < << 114 INCL_DEBUG("Cascade stopping time is " << maximumTime << std::endl);
216 115
217 // If Coulomb is activated, do not proce 116 // If Coulomb is activated, do not process events with impact
218 // parameter larger than the maximum imp 117 // parameter larger than the maximum impact parameter, taking into
219 // account Coulomb distortion. 118 // account Coulomb distortion.
220 if(impactParameter>CoulombDistortion::ma 119 if(impactParameter>CoulombDistortion::maxImpactParameter(p->getSpecies(), kineticEnergy, theNucleus)) {
221 INCL_DEBUG("impactParameter>CoulombDis << 120 INCL_DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << std::endl);
222 delete p; 121 delete p;
223 return -1.; 122 return -1.;
224 } 123 }
225 124
226 ThreeVector position(impactParameter * s 125 ThreeVector position(impactParameter * std::cos(phi),
227 impactParameter * std::sin(phi), 126 impactParameter * std::sin(phi),
228 0.); 127 0.);
229 p->setPosition(position); 128 p->setPosition(position);
230 129
231 // Fill in the relevant kinematic variab 130 // Fill in the relevant kinematic variables
232 theNucleus->setIncomingAngularMomentum(p 131 theNucleus->setIncomingAngularMomentum(p->getAngularMomentum());
233 theNucleus->setIncomingMomentum(p->getMo 132 theNucleus->setIncomingMomentum(p->getMomentum());
234 theNucleus->setInitialEnergy(p->getEnerg 133 theNucleus->setInitialEnergy(p->getEnergy()
235 + ParticleTable::getTableMass(theNuc << 134 + ParticleTable::getTableMass(theNucleus->getA(),theNucleus->getZ()));
236 135
237 // Reset the particle kinematics to the 136 // Reset the particle kinematics to the INCL values
238 p->setINCLMass(); 137 p->setINCLMass();
239 p->setEnergy(p->getMass() + kineticEnerg 138 p->setEnergy(p->getMass() + kineticEnergy);
240 p->adjustMomentumFromEnergy(); 139 p->adjustMomentumFromEnergy();
241 140
242 p->makeProjectileSpectator(); 141 p->makeProjectileSpectator();
243 generateAllAvatars(); 142 generateAllAvatars();
244 firstAvatar = false; 143 firstAvatar = false;
245 144
246 // Get the entry avatars from Coulomb an 145 // Get the entry avatars from Coulomb and put them in the Store
247 ParticleEntryAvatar *theEntryAvatar = Co 146 ParticleEntryAvatar *theEntryAvatar = CoulombDistortion::bringToSurface(p, theNucleus);
248 if(theEntryAvatar) { 147 if(theEntryAvatar) {
249 theNucleus->getStore()->addParticleEnt 148 theNucleus->getStore()->addParticleEntryAvatar(theEntryAvatar);
250 149
>> 150 theNucleus->setProjectileChargeNumber(p->getZ());
>> 151 theNucleus->setProjectileMassNumber(p->getA());
>> 152
251 return p->getTransversePosition().mag( 153 return p->getTransversePosition().mag();
252 } else { 154 } else {
253 delete p; 155 delete p;
254 return -1.; 156 return -1.;
255 } 157 }
256 } 158 }
257 159
258 G4double StandardPropagationModel::shootCo << 160 G4double StandardPropagationModel::shootComposite(ParticleSpecies const species, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
259 theNucleus->setNucleusNucleusCollision() 161 theNucleus->setNucleusNucleusCollision();
260 currentTime = 0.0; 162 currentTime = 0.0;
261 163
262 // Create the ProjectileRemnant object 164 // Create the ProjectileRemnant object
263 ProjectileRemnant *pr = new ProjectileRe 165 ProjectileRemnant *pr = new ProjectileRemnant(species, kineticEnergy);
264 166
265 // Same stopping time as for nucleon-nuc 167 // Same stopping time as for nucleon-nucleus
266 maximumTime = 29.8 * std::pow(theNucleus 168 maximumTime = 29.8 * std::pow(theNucleus->getA(), 0.16);
267 169
268 // If the incoming cluster is slow, use 170 // If the incoming cluster is slow, use a larger stopping time
269 const G4double rms = ParticleTable::getL 171 const G4double rms = ParticleTable::getLargestNuclearRadius(pr->getA(), pr->getZ());
270 const G4double rMax = theNucleus->getUni 172 const G4double rMax = theNucleus->getUniverseRadius();
271 const G4double distance = 2.*rMax + 2.72 173 const G4double distance = 2.*rMax + 2.725*rms;
272 const G4double projectileVelocity = pr-> 174 const G4double projectileVelocity = pr->boostVector().mag();
273 const G4double traversalTime = distance 175 const G4double traversalTime = distance / projectileVelocity;
274 if(maximumTime < traversalTime) 176 if(maximumTime < traversalTime)
275 maximumTime = traversalTime; 177 maximumTime = traversalTime;
276 INCL_DEBUG("Cascade stopping time is " < << 178 INCL_DEBUG("Cascade stopping time is " << maximumTime << std::endl);
277 179
278 // If Coulomb is activated, do not proce 180 // If Coulomb is activated, do not process events with impact
279 // parameter larger than the maximum imp 181 // parameter larger than the maximum impact parameter, taking into
280 // account Coulomb distortion. 182 // account Coulomb distortion.
281 if(impactParameter>CoulombDistortion::ma 183 if(impactParameter>CoulombDistortion::maxImpactParameter(pr,theNucleus)) {
282 INCL_DEBUG("impactParameter>CoulombDis << 184 INCL_DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << std::endl);
283 delete pr; 185 delete pr;
284 return -1.; 186 return -1.;
285 } 187 }
286 188
287 // Position the cluster at the right imp 189 // Position the cluster at the right impact parameter
288 ThreeVector position(impactParameter * s 190 ThreeVector position(impactParameter * std::cos(phi),
289 impactParameter * std::sin(phi), 191 impactParameter * std::sin(phi),
290 0.); 192 0.);
291 pr->setPosition(position); 193 pr->setPosition(position);
292 194
293 // Fill in the relevant kinematic variab 195 // Fill in the relevant kinematic variables
294 theNucleus->setIncomingAngularMomentum(p 196 theNucleus->setIncomingAngularMomentum(pr->getAngularMomentum());
295 theNucleus->setIncomingMomentum(pr->getM 197 theNucleus->setIncomingMomentum(pr->getMomentum());
296 theNucleus->setInitialEnergy(pr->getEner 198 theNucleus->setInitialEnergy(pr->getEnergy()
297 + ParticleTable::getTableMass(theNuc << 199 + ParticleTable::getTableMass(theNucleus->getA(),theNucleus->getZ()));
298 200
299 generateAllAvatars(); 201 generateAllAvatars();
300 firstAvatar = false; 202 firstAvatar = false;
301 203
302 // Get the entry avatars from Coulomb 204 // Get the entry avatars from Coulomb
303 IAvatarList theAvatarList 205 IAvatarList theAvatarList
304 = CoulombDistortion::bringToSurface(pr 206 = CoulombDistortion::bringToSurface(pr, theNucleus);
305 207
306 if(theAvatarList.empty()) { 208 if(theAvatarList.empty()) {
307 INCL_DEBUG("No ParticleEntryAvatar fou << 209 INCL_DEBUG("No ParticleEntryAvatar found, transparent event" << std::endl);
308 delete pr; 210 delete pr;
309 return -1.; 211 return -1.;
310 } 212 }
311 213
312 /* Store the internal kinematics of the 214 /* Store the internal kinematics of the projectile remnant.
313 * 215 *
314 * Note that this is at variance with th 216 * Note that this is at variance with the Fortran version, which stores
315 * the initial kinematics of the particl 217 * the initial kinematics of the particles *after* putting the spectators
316 * on mass shell, but *before* removing 218 * on mass shell, but *before* removing the same energy from the
317 * participants. Due to the different co 219 * participants. Due to the different code flow, doing so in the C++
318 * version leads to wrong excitation ene 220 * version leads to wrong excitation energies for the forced compound
319 * nucleus. 221 * nucleus.
320 */ 222 */
321 pr->storeComponents(); 223 pr->storeComponents();
322 224
323 // Tell the Nucleus about the Projectile 225 // Tell the Nucleus about the ProjectileRemnant
324 theNucleus->setProjectileRemnant(pr); 226 theNucleus->setProjectileRemnant(pr);
325 227
>> 228 // Set the number of projectile particles
>> 229 theNucleus->setProjectileChargeNumber(pr->getZ());
>> 230 theNucleus->setProjectileMassNumber(pr->getA());
>> 231
326 // Register the ParticleEntryAvatars 232 // Register the ParticleEntryAvatars
327 theNucleus->getStore()->addParticleEntry 233 theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
328 234
329 return pr->getTransversePosition().mag() 235 return pr->getTransversePosition().mag();
330 } 236 }
331 237
332 G4double StandardPropagationModel::getStop 238 G4double StandardPropagationModel::getStoppingTime() {
333 return maximumTime; 239 return maximumTime;
334 } 240 }
335 241
336 void StandardPropagationModel::setStopping 242 void StandardPropagationModel::setStoppingTime(G4double time) {
337 // assert(time>0.0); 243 // assert(time>0.0);
338 maximumTime = time; 244 maximumTime = time;
339 } 245 }
340 246
341 G4double StandardPropagationModel::getCurr 247 G4double StandardPropagationModel::getCurrentTime() {
342 return currentTime; 248 return currentTime;
343 } 249 }
344 250
345 void StandardPropagationModel::setNucleus( 251 void StandardPropagationModel::setNucleus(G4INCL::Nucleus *nucleus)
346 { 252 {
347 theNucleus = nucleus; 253 theNucleus = nucleus;
348 } 254 }
349 255
350 void StandardPropagationModel::registerAva 256 void StandardPropagationModel::registerAvatar(G4INCL::IAvatar *anAvatar)
351 { 257 {
352 if(anAvatar) theNucleus->getStore()->add 258 if(anAvatar) theNucleus->getStore()->add(anAvatar);
353 } 259 }
354 260
355 IAvatar *StandardPropagationModel::generat 261 IAvatar *StandardPropagationModel::generateBinaryCollisionAvatar(Particle * const p1, Particle * const p2) {
356 // Is either particle a participant? 262 // Is either particle a participant?
357 if(!p1->isParticipant() && !p2->isPartic 263 if(!p1->isParticipant() && !p2->isParticipant() && p1->getParticipantType()==p2->getParticipantType()) return NULL;
358 264
359 // Is it a pi-resonance collision (we do 265 // Is it a pi-resonance collision (we don't treat them)?
360 if((p1->isResonance() && p2->isPion()) | 266 if((p1->isResonance() && p2->isPion()) || (p1->isPion() && p2->isResonance()))
361 return NULL; 267 return NULL;
362 268
363 // Is it a photon collision (we don't tr <<
364 if(p1->isPhoton() || p2->isPhoton()) <<
365 return NULL; <<
366 <<
367 // Will the avatar take place between no 269 // Will the avatar take place between now and the end of the cascade?
368 G4double minDistOfApproachSquared = 0.0; 270 G4double minDistOfApproachSquared = 0.0;
369 G4double t = getTime(p1, p2, &minDistOfA 271 G4double t = getTime(p1, p2, &minDistOfApproachSquared);
370 if(t>maximumTime || t<currentTime+hadron << 272 if(t>maximumTime || t<currentTime) return NULL;
371 273
372 // Local energy. Jump through some hoops 274 // Local energy. Jump through some hoops to calculate the cross section
373 // at the collision point, and clean up 275 // at the collision point, and clean up after yourself afterwards.
374 G4bool hasLocalEnergy; 276 G4bool hasLocalEnergy;
375 if(p1->isPion() || p2->isPion()) 277 if(p1->isPion() || p2->isPion())
376 hasLocalEnergy = ((theLocalEnergyDelta 278 hasLocalEnergy = ((theLocalEnergyDeltaType == FirstCollisionLocalEnergy &&
377 theNucleus->getStore()->getBook( 279 theNucleus->getStore()->getBook().getAcceptedCollisions()==0) ||
378 theLocalEnergyDeltaType == AlwaysL 280 theLocalEnergyDeltaType == AlwaysLocalEnergy);
379 else 281 else
380 hasLocalEnergy = ((theLocalEnergyType 282 hasLocalEnergy = ((theLocalEnergyType == FirstCollisionLocalEnergy &&
381 theNucleus->getStore()->getBook( 283 theNucleus->getStore()->getBook().getAcceptedCollisions()==0) ||
382 theLocalEnergyType == AlwaysLocalE 284 theLocalEnergyType == AlwaysLocalEnergy);
383 const G4bool p1HasLocalEnergy = (hasLoca << 285 const G4bool p1HasLocalEnergy = (hasLocalEnergy && !p1->isPion());
384 const G4bool p2HasLocalEnergy = (hasLoca << 286 const G4bool p2HasLocalEnergy = (hasLocalEnergy && !p2->isPion());
385 287
386 if(p1HasLocalEnergy) { 288 if(p1HasLocalEnergy) {
387 backupParticle1 = *p1; 289 backupParticle1 = *p1;
388 p1->propagate(t - currentTime); 290 p1->propagate(t - currentTime);
389 if(p1->getPosition().mag() > theNucleu 291 if(p1->getPosition().mag() > theNucleus->getSurfaceRadius(p1)) {
390 *p1 = backupParticle1; 292 *p1 = backupParticle1;
391 return NULL; 293 return NULL;
392 } 294 }
393 KinematicsUtils::transformToLocalEnerg 295 KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p1);
394 } << 296 }
395 if(p2HasLocalEnergy) { 297 if(p2HasLocalEnergy) {
396 backupParticle2 = *p2; 298 backupParticle2 = *p2;
397 p2->propagate(t - currentTime); 299 p2->propagate(t - currentTime);
398 if(p2->getPosition().mag() > theNucleu 300 if(p2->getPosition().mag() > theNucleus->getSurfaceRadius(p2)) {
399 *p2 = backupParticle2; 301 *p2 = backupParticle2;
400 if(p1HasLocalEnergy) { 302 if(p1HasLocalEnergy) {
401 *p1 = backupParticle1; 303 *p1 = backupParticle1;
402 } 304 }
403 return NULL; 305 return NULL;
404 } 306 }
405 KinematicsUtils::transformToLocalEnerg 307 KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p2);
406 } 308 }
407 309
408 // Compute the total cross section 310 // Compute the total cross section
409 const G4double totalCrossSection = Cross 311 const G4double totalCrossSection = CrossSections::total(p1, p2);
410 const G4double squareTotalEnergyInCM = K 312 const G4double squareTotalEnergyInCM = KinematicsUtils::squareTotalEnergyInCM(p1,p2);
411 313
412 // Restore particles to their state befo 314 // Restore particles to their state before the local-energy tweak
413 if(p1HasLocalEnergy) { 315 if(p1HasLocalEnergy) {
414 *p1 = backupParticle1; 316 *p1 = backupParticle1;
415 } 317 }
416 if(p2HasLocalEnergy) { 318 if(p2HasLocalEnergy) {
417 *p2 = backupParticle2; 319 *p2 = backupParticle2;
418 } 320 }
419 321
420 // Is the CM energy > cutNN? (no cutNN o 322 // Is the CM energy > cutNN? (no cutNN on the first collision)
421 if(theNucleus->getStore()->getBook().get 323 if(theNucleus->getStore()->getBook().getAcceptedCollisions()>0
422 && p1->isNucleon() && p2->isNucleon( 324 && p1->isNucleon() && p2->isNucleon()
423 && squareTotalEnergyInCM < BinaryCol 325 && squareTotalEnergyInCM < BinaryCollisionAvatar::getCutNNSquared()) return NULL;
424 326
425 // Do the particles come close enough to 327 // Do the particles come close enough to each other?
426 if(Math::tenPi*minDistOfApproachSquared 328 if(Math::tenPi*minDistOfApproachSquared > totalCrossSection) return NULL;
427 329
428 // Bomb out if the two collision partner 330 // Bomb out if the two collision partners are the same particle
429 // assert(p1->getID() != p2->getID()); 331 // assert(p1->getID() != p2->getID());
430 332
431 // Return a new avatar, then! 333 // Return a new avatar, then!
432 return new G4INCL::BinaryCollisionAvatar 334 return new G4INCL::BinaryCollisionAvatar(t, totalCrossSection, theNucleus, p1, p2);
433 } 335 }
434 336
435 G4double StandardPropagationModel::getRefl 337 G4double StandardPropagationModel::getReflectionTime(G4INCL::Particle const * const aParticle) {
436 Intersection theIntersection( 338 Intersection theIntersection(
437 IntersectionFactory::getLaterTraject 339 IntersectionFactory::getLaterTrajectoryIntersection(
438 aParticle->getPosition(), 340 aParticle->getPosition(),
439 aParticle->getPropagationVelocity( 341 aParticle->getPropagationVelocity(),
440 theNucleus->getSurfaceRadius(aPart 342 theNucleus->getSurfaceRadius(aParticle)));
441 G4double time; 343 G4double time;
442 if(theIntersection.exists) { 344 if(theIntersection.exists) {
443 time = currentTime + theIntersection.t 345 time = currentTime + theIntersection.time;
444 } else { 346 } else {
445 INCL_ERROR("Imaginary reflection time << 347 INCL_ERROR("Imaginary reflection time for particle: " << std::endl
446 << aParticle->print()); 348 << aParticle->print());
447 time = 10000.0; 349 time = 10000.0;
448 } 350 }
449 return time; 351 return time;
450 } 352 }
451 353
452 G4double StandardPropagationModel::getTime 354 G4double StandardPropagationModel::getTime(G4INCL::Particle const * const particleA,
453 G4INCL::Particle const * const 355 G4INCL::Particle const * const particleB, G4double *minDistOfApproach) const
454 { 356 {
455 G4double time; 357 G4double time;
456 G4INCL::ThreeVector t13 = particleA->get 358 G4INCL::ThreeVector t13 = particleA->getPropagationVelocity();
457 t13 -= particleB->getPropagationVelocity 359 t13 -= particleB->getPropagationVelocity();
458 G4INCL::ThreeVector distance = particleA 360 G4INCL::ThreeVector distance = particleA->getPosition();
459 distance -= particleB->getPosition(); 361 distance -= particleB->getPosition();
460 const G4double t7 = t13.dot(distance); 362 const G4double t7 = t13.dot(distance);
461 const G4double dt = t13.mag2(); 363 const G4double dt = t13.mag2();
462 if(dt <= 1.0e-10) { 364 if(dt <= 1.0e-10) {
463 (*minDistOfApproach) = 100000.0; 365 (*minDistOfApproach) = 100000.0;
464 return currentTime + 100000.0; 366 return currentTime + 100000.0;
465 } else { 367 } else {
466 time = -t7/dt; 368 time = -t7/dt;
467 } 369 }
468 (*minDistOfApproach) = distance.mag2() + 370 (*minDistOfApproach) = distance.mag2() + time * t7;
469 return currentTime + time; 371 return currentTime + time;
470 } 372 }
471 373
472 void StandardPropagationModel::generateUpd 374 void StandardPropagationModel::generateUpdatedCollisions(const ParticleList &updatedParticles, const ParticleList &particles) {
473 375
474 // Loop over all the updated particles 376 // Loop over all the updated particles
475 for(ParticleIter updated=updatedParticle 377 for(ParticleIter updated=updatedParticles.begin(), e=updatedParticles.end(); updated!=e; ++updated)
476 { 378 {
477 // Loop over all the particles 379 // Loop over all the particles
478 for(ParticleIter particle=particles.be 380 for(ParticleIter particle=particles.begin(), end=particles.end(); particle!=end; ++particle)
479 { 381 {
480 /* Consider the generation of a coll 382 /* Consider the generation of a collision avatar only if (*particle)
481 * is not one of the updated particl 383 * is not one of the updated particles.
482 * The criterion makes sure that you 384 * The criterion makes sure that you don't generate avatars between
483 * updated particles. */ 385 * updated particles. */
484 if(updatedParticles.contains(*partic << 386 if((*particle)->isInList(updatedParticles)) continue;
485 387
486 registerAvatar(generateBinaryCollisi 388 registerAvatar(generateBinaryCollisionAvatar(*particle,*updated));
487 } 389 }
488 } 390 }
489 } 391 }
490 392
491 void StandardPropagationModel::generateCol 393 void StandardPropagationModel::generateCollisions(const ParticleList &particles) {
492 // Loop over all the particles 394 // Loop over all the particles
493 for(ParticleIter p1=particles.begin(), e 395 for(ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1) {
494 // Loop over the rest of the particles 396 // Loop over the rest of the particles
495 for(ParticleIter p2 = p1 + 1; p2 != pa 397 for(ParticleIter p2 = p1 + 1; p2 != particles.end(); ++p2) {
496 registerAvatar(generateBinaryCollisi 398 registerAvatar(generateBinaryCollisionAvatar(*p1,*p2));
497 } 399 }
498 } 400 }
499 } 401 }
500 402
501 void StandardPropagationModel::generateCol 403 void StandardPropagationModel::generateCollisions(const ParticleList &particles, const ParticleList &except) {
502 404
503 const G4bool haveExcept = !except.empty( 405 const G4bool haveExcept = !except.empty();
504 406
505 // Loop over all the particles 407 // Loop over all the particles
506 for(ParticleIter p1=particles.begin(), e 408 for(ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1)
507 { 409 {
508 // Loop over the rest of the particles 410 // Loop over the rest of the particles
509 ParticleIter p2 = p1; 411 ParticleIter p2 = p1;
510 for(++p2; p2 != particles.end(); ++p2) 412 for(++p2; p2 != particles.end(); ++p2)
511 { 413 {
512 // Skip the collision if both partic 414 // Skip the collision if both particles must be excluded
513 if(haveExcept && except.contains(*p1 << 415 if(haveExcept && (*p1)->isInList(except) && (*p2)->isInList(except)) continue;
514 416
515 registerAvatar(generateBinaryCollisi 417 registerAvatar(generateBinaryCollisionAvatar(*p1,*p2));
516 } 418 }
517 } 419 }
518 420
519 } 421 }
520 422
521 void StandardPropagationModel::updateAvata 423 void StandardPropagationModel::updateAvatars(const ParticleList &particles) {
522 424
523 for(ParticleIter iter=particles.begin(), 425 for(ParticleIter iter=particles.begin(), e=particles.end(); iter!=e; ++iter) {
524 G4double time = this->getReflectionTim 426 G4double time = this->getReflectionTime(*iter);
525 if(time <= maximumTime) registerAvatar 427 if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*iter, time, theNucleus));
526 } 428 }
527 ParticleList const &p = theNucleus->getS 429 ParticleList const &p = theNucleus->getStore()->getParticles();
528 generateUpdatedCollisions(particles, p); 430 generateUpdatedCollisions(particles, p); // Predict collisions with spectators and participants
529 } 431 }
530 432
531 void StandardPropagationModel::generateAll 433 void StandardPropagationModel::generateAllAvatars() {
532 ParticleList const &particles = theNucle 434 ParticleList const &particles = theNucleus->getStore()->getParticles();
533 // assert(!particles.empty()); 435 // assert(!particles.empty());
534 G4double time; <<
535 for(ParticleIter i=particles.begin(), e= 436 for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
536 time = this->getReflectionTime(*i); << 437 G4double time = this->getReflectionTime(*i);
537 if(time <= maximumTime) registerAvatar 438 if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
538 } 439 }
539 generateCollisions(particles); 440 generateCollisions(particles);
540 generateDecays(particles); 441 generateDecays(particles);
541 } 442 }
542 443
543 #ifdef INCL_REGENERATE_AVATARS 444 #ifdef INCL_REGENERATE_AVATARS
544 void StandardPropagationModel::generateAll << 445 void StandardPropagationModel::generateAllAvatarsExceptUpdated() {
545 ParticleList const &particles = theNucle 446 ParticleList const &particles = theNucleus->getStore()->getParticles();
546 // assert(!particles.empty()); 447 // assert(!particles.empty());
547 for(ParticleIter i=particles.begin(), e= 448 for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
548 G4double time = this->getReflectionTim 449 G4double time = this->getReflectionTime(*i);
549 if(time <= maximumTime) registerAvatar 450 if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
550 } 451 }
551 ParticleList except = fs->getModifiedPar << 452 ParticleList except = theNucleus->getUpdatedParticles();
552 ParticleList const &entering = fs->getEn <<
553 except.insert(except.end(), entering.beg <<
554 generateCollisions(particles,except); 453 generateCollisions(particles,except);
555 generateDecays(particles); 454 generateDecays(particles);
556 } 455 }
557 #endif 456 #endif
558 457
559 void StandardPropagationModel::generateDec 458 void StandardPropagationModel::generateDecays(const ParticleList &particles) {
560 for(ParticleIter i=particles.begin(), e= 459 for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
561 if((*i)->isDelta()) { << 460 if((*i)->isDelta()) {
562 G4double decayTime = DeltaDecayCh << 461 G4double decayTime = DeltaDecayChannel::computeDecayTime((*i));
563 G4double time = currentTime + decay << 462 G4double time = currentTime + decayTime;
564 if(time <= maximumTime) { << 463 if(time <= maximumTime) {
565 registerAvatar(new DecayAvatar((* << 464 registerAvatar(new DecayAvatar((*i), time, theNucleus));
566 } << 465 }
567 } << 466 }
568 else if((*i)->getType() == SigmaZero) <<
569 G4double decayTime = SigmaZeroDecay <<
570 G4double time = currentTime + decay <<
571 if(time <= maximumTime) { <<
572 registerAvatar(new DecayAvatar((* <<
573 } <<
574 } <<
575 if((*i)->isOmega()) { <<
576 G4double decayTimeOmega = PionResona <<
577 G4double timeOmega = currentTime + d <<
578 if(timeOmega <= maximumTime) { <<
579 registerAvatar(new DecayAvatar((*i <<
580 } <<
581 } <<
582 } 467 }
583 } 468 }
584 469
585 G4INCL::IAvatar* StandardPropagationModel: << 470 G4INCL::IAvatar* StandardPropagationModel::propagate()
586 { 471 {
587 if(fs) { << 472 // We update only the information related to particles that were updated
588 // We update only the information rela << 473 // by the previous avatar.
589 // by the previous avatar. <<
590 #ifdef INCL_REGENERATE_AVATARS 474 #ifdef INCL_REGENERATE_AVATARS
591 #warning "The INCL_REGENERATE_AVATARS code has 475 #warning "The INCL_REGENERATE_AVATARS code has not been tested in a while. Use it at your peril."
592 if(!fs->getModifiedParticles().empty() << 476 if(theNucleus->getUpdatedParticles().size()!=0 || theNucleus->getCreatedParticles().size()!=0) {
593 // Regenerates the entire avatar lis << 477 // Regenerates the entire avatar list, skipping collisions between
594 // updated particles << 478 // updated particles
595 theNucleus->getStore()->clearAvatars << 479 theNucleus->getStore()->clearAvatars();
596 theNucleus->getStore()->initialisePa << 480 theNucleus->getStore()->initialiseParticleAvatarConnections();
597 generateAllAvatarsExceptUpdated(fs); << 481 generateAllAvatarsExceptUpdated();
598 } << 482 }
599 #else 483 #else
600 ParticleList const &updatedParticles = << 484 // Deltas are created by transforming nucleon into a delta for
601 if(fs->getValidity()==PauliBlockedFS) << 485 // efficiency reasons
602 // This final state might represents << 486 Particle * const blockedDelta = theNucleus->getBlockedDelta();
603 // decay << 487 ParticleList updatedParticles = theNucleus->getUpdatedParticles();
604 // assert(updatedParticles.empty() || (updated << 488 if(blockedDelta)
605 // assert(fs->getEnteringParticles().empty() & << 489 updatedParticles.push_back(blockedDelta);
606 generateDecays(updatedParticles); << 490 generateDecays(updatedParticles);
607 } else { << 491
608 ParticleList const &entering = fs->g << 492 ParticleList needNewAvatars = theNucleus->getUpdatedParticles();
609 generateDecays(updatedParticles); << 493 ParticleList const &created = theNucleus->getCreatedParticles();
610 generateDecays(entering); << 494 needNewAvatars.insert(needNewAvatars.end(), created.begin(), created.end());
611 << 495 updateAvatars(needNewAvatars);
612 ParticleList const &created = fs->ge <<
613 if(created.empty() && entering.empty <<
614 updateAvatars(updatedParticles); <<
615 else { <<
616 ParticleList updatedParticlesCopy <<
617 updatedParticlesCopy.insert(update <<
618 updatedParticlesCopy.insert(update <<
619 updateAvatars(updatedParticlesCopy <<
620 } <<
621 } <<
622 #endif 496 #endif
623 } <<
624 497
625 G4INCL::IAvatar *theAvatar = theNucleus- 498 G4INCL::IAvatar *theAvatar = theNucleus->getStore()->findSmallestTime();
626 if(theAvatar == 0) return 0; // Avatar l 499 if(theAvatar == 0) return 0; // Avatar list is empty
627 // theAvatar->dispose(); 500 // theAvatar->dispose();
628 501
629 if(theAvatar->getTime() < currentTime) { 502 if(theAvatar->getTime() < currentTime) {
630 INCL_ERROR("Avatar time = " << theAvat << 503 INCL_ERROR("Avatar time = " << theAvatar->getTime() << ", currentTime = " << currentTime << std::endl);
631 return 0; 504 return 0;
632 } else if(theAvatar->getTime() > current 505 } else if(theAvatar->getTime() > currentTime) {
633 theNucleus->getStore()->timeStep(theAv 506 theNucleus->getStore()->timeStep(theAvatar->getTime() - currentTime);
634 507
635 currentTime = theAvatar->getTime(); 508 currentTime = theAvatar->getTime();
636 theNucleus->getStore()->getBook().setC 509 theNucleus->getStore()->getBook().setCurrentTime(currentTime);
637 } 510 }
638 511
639 return theAvatar; 512 return theAvatar;
>> 513 }
>> 514
>> 515 void StandardPropagationModel::putSpectatorsOnShell(IAvatarList const &entryAvatars, ParticleList const &spectators) {
>> 516 G4double deltaE = 0.0;
>> 517 for(ParticleIter p=spectators.begin(), e=spectators.end(); p!=e; ++p) {
>> 518 // put the spectators on shell (conserving their momentum)
>> 519 const G4double oldEnergy = (*p)->getEnergy();
>> 520 (*p)->setTableMass();
>> 521 (*p)->adjustEnergyFromMomentum();
>> 522 deltaE += (*p)->getEnergy() - oldEnergy;
>> 523 }
>> 524
>> 525 deltaE /= entryAvatars.size(); // energy to remove from each participant
>> 526
>> 527 for(IAvatarIter a=entryAvatars.begin(), e=entryAvatars.end(); a!=e; ++a) {
>> 528 // remove the energy from the participant
>> 529 Particle *p = (*a)->getParticles().front();
>> 530 ParticleType const t = p->getType();
>> 531 // we also need to slightly correct the participant mass
>> 532 const G4double energy = p->getEnergy() - deltaE
>> 533 - ParticleTable::getTableParticleMass(t) + ParticleTable::getINCLMass(t);
>> 534 p->setEnergy(energy);
>> 535 const G4double newMass = std::sqrt(energy*energy - p->getMomentum().mag2());
>> 536 p->setMass(newMass);
>> 537 }
640 } 538 }
641 } 539 }
642 540