| 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 << 32 //
>> 33 // INCL++ revision: v5.0_rc3
33 // 34 //
34 #define INCLXX_IN_GEANT4_MODE 1 35 #define INCLXX_IN_GEANT4_MODE 1
35 36
36 #include "globals.hh" 37 #include "globals.hh"
37 38
38 /* 39 /*
39 * StandardPropagationModel.cpp 40 * StandardPropagationModel.cpp
40 * 41 *
41 * \date 4 juin 2009 << 42 * Created on: 4 juin 2009
42 * \author Pekka Kaitaniemi << 43 * Author: Pekka Kaitaniemi
43 */ 44 */
44 45
45 #include "G4INCLStandardPropagationModel.hh" 46 #include "G4INCLStandardPropagationModel.hh"
46 #include "G4INCLPbarAtrestEntryChannel.hh" <<
47 #include "G4INCLSurfaceAvatar.hh" 47 #include "G4INCLSurfaceAvatar.hh"
48 #include "G4INCLBinaryCollisionAvatar.hh" 48 #include "G4INCLBinaryCollisionAvatar.hh"
49 #include "G4INCLDecayAvatar.hh" 49 #include "G4INCLDecayAvatar.hh"
50 #include "G4INCLCrossSections.hh" 50 #include "G4INCLCrossSections.hh"
51 #include "G4INCLRandom.hh" 51 #include "G4INCLRandom.hh"
52 #include <iostream> 52 #include <iostream>
53 #include <algorithm> <<
54 #include "G4INCLLogger.hh" 53 #include "G4INCLLogger.hh"
55 #include "G4INCLGlobals.hh" 54 #include "G4INCLGlobals.hh"
56 #include "G4INCLKinematicsUtils.hh" 55 #include "G4INCLKinematicsUtils.hh"
57 #include "G4INCLCoulombDistortion.hh" 56 #include "G4INCLCoulombDistortion.hh"
58 #include "G4INCLDeltaDecayChannel.hh" 57 #include "G4INCLDeltaDecayChannel.hh"
59 #include "G4INCLSigmaZeroDecayChannel.hh" <<
60 #include "G4INCLPionResonanceDecayChannel.hh" <<
61 #include "G4INCLParticleEntryAvatar.hh" <<
62 #include "G4INCLIntersection.hh" <<
63 #include <vector> <<
64 58
65 namespace G4INCL { 59 namespace G4INCL {
66 60
67 StandardPropagationModel::StandardPropagat << 61 StandardPropagationModel::StandardPropagationModel(LocalEnergyType localEnergyType, LocalEnergyType localEnergyDeltaType)
68 :theNucleus(0), maximumTime(70.0), curre << 62 :theNucleus(0), maximumTime(70.0), currentTime(0.0), firstAvatar(true),
69 hadronizationTime(hTime), <<
70 firstAvatar(true), <<
71 theLocalEnergyType(localEnergyType), 63 theLocalEnergyType(localEnergyType),
72 theLocalEnergyDeltaType(localEnergyDelta 64 theLocalEnergyDeltaType(localEnergyDeltaType)
73 { 65 {
74 } 66 }
75 67
76 StandardPropagationModel::~StandardPropaga 68 StandardPropagationModel::~StandardPropagationModel()
77 { 69 {
78 delete theNucleus; 70 delete theNucleus;
79 } 71 }
80 72
81 G4INCL::Nucleus* StandardPropagationModel: 73 G4INCL::Nucleus* StandardPropagationModel::getNucleus()
82 { 74 {
83 return theNucleus; 75 return theNucleus;
84 } 76 }
85 77
86 //D << 78 G4bool StandardPropagationModel::shootProjectile(G4INCL::Particle *p, G4double impactParameter) {
87 << 79 firstAvatar = true;
88 G4double StandardPropagationModel::shoot(P <<
89 if(projectileSpecies.theType==Composite) <<
90 return shootComposite(projectileSpecie <<
91 } <<
92 else if(projectileSpecies.theType==antiP <<
93 return shootAtrest(projectileSpecies.t <<
94 } <<
95 else{ <<
96 return shootParticle(projectileSpecies <<
97 } <<
98 } <<
99 <<
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 <<
182 theNucleus->setParticleNucleusCollision( <<
183 currentTime = 0.0; 80 currentTime = 0.0;
184 81
185 // Create the projectile particle << 82 G4double temfin = 0.0;
186 const G4double projectileMass = Particle << 83 if( p->isNucleon() )
187 G4double energy = kineticEnergy + projec <<
188 G4double momentumZ = std::sqrt(energy*en <<
189 ThreeVector momentum(0.0, 0.0, momentumZ <<
190 Particle *p= new G4INCL::Particle(type, <<
191 <<
192 G4double temfin; <<
193 G4double TLab; <<
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 84 temfin = 29.8 * std::pow(theNucleus->getA(), 0.16);
199 TLab = p->getKineticEnergy()/p->getA() << 85 else {
>> 86 const G4double tlab = p->getEnergy() - p->getMass();
>> 87 temfin = 30.18 * std::pow(theNucleus->getA(), 0.17*(1.0 - 5.7E-5*tlab));
200 } 88 }
201 89
202 // energy-dependent stopping time above <<
203 if(TLab>2000.) <<
204 temfin *= (5.8E4-TLab)/5.6E4; <<
205 <<
206 maximumTime = temfin; 90 maximumTime = temfin;
207 91
208 // If the incoming particle is slow, use <<
209 const G4double rMax = theNucleus->getUni <<
210 const G4double distance = 2.*rMax; <<
211 const G4double projectileVelocity = p->b <<
212 const G4double traversalTime = distance <<
213 if(maximumTime < traversalTime) <<
214 maximumTime = traversalTime; <<
215 INCL_DEBUG("Cascade stopping time is " < <<
216 <<
217 // If Coulomb is activated, do not proce 92 // If Coulomb is activated, do not process events with impact
218 // parameter larger than the maximum imp << 93 // parameter larger than the maximum impact parameter, taking G4into
219 // account Coulomb distortion. 94 // account Coulomb distortion.
220 if(impactParameter>CoulombDistortion::ma << 95 if(impactParameter>CoulombDistortion::maxImpactParameter(p,theNucleus))
221 INCL_DEBUG("impactParameter>CoulombDis << 96 return false;
222 delete p; <<
223 return -1.; <<
224 } <<
225 97
226 ThreeVector position(impactParameter * s << 98 const G4double tbid = Random::shoot() * Math::twoPi;
227 impactParameter * std::sin(phi), << 99 ThreeVector position(impactParameter * std::cos(tbid),
228 0.); << 100 impactParameter * std::sin(tbid),
>> 101 -1.E3);
229 p->setPosition(position); 102 p->setPosition(position);
230 103
231 // Fill in the relevant kinematic variab <<
232 theNucleus->setIncomingAngularMomentum(p 104 theNucleus->setIncomingAngularMomentum(p->getAngularMomentum());
233 theNucleus->setIncomingMomentum(p->getMo 105 theNucleus->setIncomingMomentum(p->getMomentum());
234 theNucleus->setInitialEnergy(p->getEnerg << 106 theNucleus->setInitialEnergy(p->getEnergy() + ParticleTable::getMass(theNucleus->getA(),theNucleus->getZ()));
235 + ParticleTable::getTableMass(theNuc <<
236 107
237 // Reset the particle kinematics to the << 108 CoulombDistortion::bringToSurface(p, theNucleus);
238 p->setINCLMass(); <<
239 p->setEnergy(p->getMass() + kineticEnerg <<
240 p->adjustMomentumFromEnergy(); <<
241 <<
242 p->makeProjectileSpectator(); <<
243 generateAllAvatars(); <<
244 firstAvatar = false; <<
245 <<
246 // Get the entry avatars from Coulomb an <<
247 ParticleEntryAvatar *theEntryAvatar = Co <<
248 if(theEntryAvatar) { <<
249 theNucleus->getStore()->addParticleEnt <<
250 109
251 return p->getTransversePosition().mag( << 110 theNucleus->getStore()->addIncomingParticle(p); // puts the particle in the waiting list
252 } else { << 111 theNucleus->particleEnters(p); // removes the particle from the waiting list and changes its kinetic energy
253 delete p; << 112 theNucleus->insertParticipant(p);
254 return -1.; << 113 return true;
255 } <<
256 } 114 }
257 115
258 G4double StandardPropagationModel::shootCo << 116 G4bool StandardPropagationModel::shootProjectile(G4INCL::Nucleus * /* p */, G4double /* impactParameter */) {
259 theNucleus->setNucleusNucleusCollision() << 117 firstAvatar = true;
260 currentTime = 0.0; 118 currentTime = 0.0;
261 << 119 return true;
262 // Create the ProjectileRemnant object <<
263 ProjectileRemnant *pr = new ProjectileRe <<
264 <<
265 // Same stopping time as for nucleon-nuc <<
266 maximumTime = 29.8 * std::pow(theNucleus <<
267 <<
268 // If the incoming cluster is slow, use <<
269 const G4double rms = ParticleTable::getL <<
270 const G4double rMax = theNucleus->getUni <<
271 const G4double distance = 2.*rMax + 2.72 <<
272 const G4double projectileVelocity = pr-> <<
273 const G4double traversalTime = distance <<
274 if(maximumTime < traversalTime) <<
275 maximumTime = traversalTime; <<
276 INCL_DEBUG("Cascade stopping time is " < <<
277 <<
278 // If Coulomb is activated, do not proce <<
279 // parameter larger than the maximum imp <<
280 // account Coulomb distortion. <<
281 if(impactParameter>CoulombDistortion::ma <<
282 INCL_DEBUG("impactParameter>CoulombDis <<
283 delete pr; <<
284 return -1.; <<
285 } <<
286 <<
287 // Position the cluster at the right imp <<
288 ThreeVector position(impactParameter * s <<
289 impactParameter * std::sin(phi), <<
290 0.); <<
291 pr->setPosition(position); <<
292 <<
293 // Fill in the relevant kinematic variab <<
294 theNucleus->setIncomingAngularMomentum(p <<
295 theNucleus->setIncomingMomentum(pr->getM <<
296 theNucleus->setInitialEnergy(pr->getEner <<
297 + ParticleTable::getTableMass(theNuc <<
298 <<
299 generateAllAvatars(); <<
300 firstAvatar = false; <<
301 <<
302 // Get the entry avatars from Coulomb <<
303 IAvatarList theAvatarList <<
304 = CoulombDistortion::bringToSurface(pr <<
305 <<
306 if(theAvatarList.empty()) { <<
307 INCL_DEBUG("No ParticleEntryAvatar fou <<
308 delete pr; <<
309 return -1.; <<
310 } <<
311 <<
312 /* Store the internal kinematics of the <<
313 * <<
314 * Note that this is at variance with th <<
315 * the initial kinematics of the particl <<
316 * on mass shell, but *before* removing <<
317 * participants. Due to the different co <<
318 * version leads to wrong excitation ene <<
319 * nucleus. <<
320 */ <<
321 pr->storeComponents(); <<
322 <<
323 // Tell the Nucleus about the Projectile <<
324 theNucleus->setProjectileRemnant(pr); <<
325 <<
326 // Register the ParticleEntryAvatars <<
327 theNucleus->getStore()->addParticleEntry <<
328 <<
329 return pr->getTransversePosition().mag() <<
330 } 120 }
331 121
332 G4double StandardPropagationModel::getStop 122 G4double StandardPropagationModel::getStoppingTime() {
333 return maximumTime; 123 return maximumTime;
334 } 124 }
335 125
336 void StandardPropagationModel::setStopping 126 void StandardPropagationModel::setStoppingTime(G4double time) {
337 // assert(time>0.0); << 127 if(time > 0.0) {
338 maximumTime = time; << 128 maximumTime = time;
>> 129 } else {
>> 130 ERROR("new stopping time is smaller than 0!" << std::endl);
>> 131 }
339 } 132 }
340 133
341 G4double StandardPropagationModel::getCurr 134 G4double StandardPropagationModel::getCurrentTime() {
342 return currentTime; 135 return currentTime;
343 } 136 }
344 137
345 void StandardPropagationModel::setNucleus( 138 void StandardPropagationModel::setNucleus(G4INCL::Nucleus *nucleus)
346 { 139 {
347 theNucleus = nucleus; 140 theNucleus = nucleus;
348 } 141 }
349 142
350 void StandardPropagationModel::registerAva 143 void StandardPropagationModel::registerAvatar(G4INCL::IAvatar *anAvatar)
351 { 144 {
352 if(anAvatar) theNucleus->getStore()->add 145 if(anAvatar) theNucleus->getStore()->add(anAvatar);
353 } 146 }
354 147
355 IAvatar *StandardPropagationModel::generat << 148 IAvatar *StandardPropagationModel::generateBinaryCollisionAvatar(Particle * const p1, Particle * const p2) const {
356 // Is either particle a participant? 149 // Is either particle a participant?
357 if(!p1->isParticipant() && !p2->isPartic << 150 if(!p1->isParticipant() && !p2->isParticipant()) return NULL;
358 151
359 // Is it a pi-resonance collision (we do 152 // Is it a pi-resonance collision (we don't treat them)?
360 if((p1->isResonance() && p2->isPion()) | 153 if((p1->isResonance() && p2->isPion()) || (p1->isPion() && p2->isResonance()))
361 return NULL; 154 return NULL;
362 155
363 // Is it a photon collision (we don't tr << 156 // 2N < Tf
364 if(p1->isPhoton() || p2->isPhoton()) << 157 if(
>> 158 (p1->isNucleon() && p1->getKineticEnergy()<theNucleus->getPotential()->getFermiEnergy(p1->getType())) &&
>> 159 (p2->isNucleon() && p2->getKineticEnergy()<theNucleus->getPotential()->getFermiEnergy(p2->getType()))
>> 160 )
365 return NULL; 161 return NULL;
366 162
>> 163 // Is the CM energy > cutNN? (no cutNN on the first collision)
>> 164 if(theNucleus->getStore()->getBook()->getAcceptedCollisions()>0
>> 165 && p1->isNucleon() && p2->isNucleon()
>> 166 && KinematicsUtils::squareTotalEnergyInCM(p1,p2) < BinaryCollisionAvatar::cutNNSquared) return NULL;
>> 167
367 // Will the avatar take place between no 168 // Will the avatar take place between now and the end of the cascade?
368 G4double minDistOfApproachSquared = 0.0; 169 G4double minDistOfApproachSquared = 0.0;
369 G4double t = getTime(p1, p2, &minDistOfA 170 G4double t = getTime(p1, p2, &minDistOfApproachSquared);
370 if(t>maximumTime || t<currentTime+hadron << 171 if(t>maximumTime || t<currentTime) return NULL;
371 172
372 // Local energy. Jump through some hoops 173 // Local energy. Jump through some hoops to calculate the cross section
373 // at the collision point, and clean up << 174 // at the collision poG4int, and clean up after yourself afterwards.
>> 175 ThreeVector mom1, mom2, pos1, pos2;
>> 176 G4double energy1 = 0.0, energy2 = 0.0;
374 G4bool hasLocalEnergy; 177 G4bool hasLocalEnergy;
375 if(p1->isPion() || p2->isPion()) 178 if(p1->isPion() || p2->isPion())
376 hasLocalEnergy = ((theLocalEnergyDelta 179 hasLocalEnergy = ((theLocalEnergyDeltaType == FirstCollisionLocalEnergy &&
377 theNucleus->getStore()->getBook( << 180 theNucleus->getStore()->getBook()->getAcceptedCollisions()==0) ||
378 theLocalEnergyDeltaType == AlwaysL 181 theLocalEnergyDeltaType == AlwaysLocalEnergy);
379 else 182 else
380 hasLocalEnergy = ((theLocalEnergyType 183 hasLocalEnergy = ((theLocalEnergyType == FirstCollisionLocalEnergy &&
381 theNucleus->getStore()->getBook( << 184 theNucleus->getStore()->getBook()->getAcceptedCollisions()==0) ||
382 theLocalEnergyType == AlwaysLocalE 185 theLocalEnergyType == AlwaysLocalEnergy);
383 const G4bool p1HasLocalEnergy = (hasLoca << 186 const G4bool p1HasLocalEnergy = (hasLocalEnergy && !p1->isPion());
384 const G4bool p2HasLocalEnergy = (hasLoca << 187 const G4bool p2HasLocalEnergy = (hasLocalEnergy && !p2->isPion());
385 188
386 if(p1HasLocalEnergy) { 189 if(p1HasLocalEnergy) {
387 backupParticle1 = *p1; << 190 mom1 = p1->getMomentum();
>> 191 pos1 = p1->getPosition();
>> 192 energy1 = p1->getEnergy();
388 p1->propagate(t - currentTime); 193 p1->propagate(t - currentTime);
389 if(p1->getPosition().mag() > theNucleu 194 if(p1->getPosition().mag() > theNucleus->getSurfaceRadius(p1)) {
390 *p1 = backupParticle1; << 195 p1->setPosition(pos1);
>> 196 p1->setMomentum(mom1);
>> 197 p1->setEnergy(energy1);
391 return NULL; 198 return NULL;
392 } 199 }
393 KinematicsUtils::transformToLocalEnerg 200 KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p1);
394 } << 201 }
395 if(p2HasLocalEnergy) { 202 if(p2HasLocalEnergy) {
396 backupParticle2 = *p2; << 203 energy2 = p2->getEnergy();
>> 204 mom2 = p2->getMomentum();
>> 205 pos2 = p2->getPosition();
397 p2->propagate(t - currentTime); 206 p2->propagate(t - currentTime);
398 if(p2->getPosition().mag() > theNucleu 207 if(p2->getPosition().mag() > theNucleus->getSurfaceRadius(p2)) {
399 *p2 = backupParticle2; << 208 p2->setPosition(pos2);
>> 209 p2->setMomentum(mom2);
>> 210 p2->setEnergy(energy2);
400 if(p1HasLocalEnergy) { 211 if(p1HasLocalEnergy) {
401 *p1 = backupParticle1; << 212 p1->setPosition(pos1);
>> 213 p1->setMomentum(mom1);
>> 214 p1->setEnergy(energy1);
402 } 215 }
403 return NULL; 216 return NULL;
404 } 217 }
405 KinematicsUtils::transformToLocalEnerg 218 KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p2);
406 } 219 }
407 220
408 // Compute the total cross section 221 // Compute the total cross section
409 const G4double totalCrossSection = Cross 222 const G4double totalCrossSection = CrossSections::total(p1, p2);
410 const G4double squareTotalEnergyInCM = K <<
411 223
412 // Restore particles to their state befo 224 // Restore particles to their state before the local-energy tweak
413 if(p1HasLocalEnergy) { 225 if(p1HasLocalEnergy) {
414 *p1 = backupParticle1; << 226 p1->setPosition(pos1);
>> 227 p1->setMomentum(mom1);
>> 228 p1->setEnergy(energy1);
415 } 229 }
416 if(p2HasLocalEnergy) { 230 if(p2HasLocalEnergy) {
417 *p2 = backupParticle2; << 231 p2->setPosition(pos2);
>> 232 p2->setMomentum(mom2);
>> 233 p2->setEnergy(energy2);
418 } 234 }
419 235
420 // Is the CM energy > cutNN? (no cutNN o <<
421 if(theNucleus->getStore()->getBook().get <<
422 && p1->isNucleon() && p2->isNucleon( <<
423 && squareTotalEnergyInCM < BinaryCol <<
424 <<
425 // Do the particles come close enough to 236 // Do the particles come close enough to each other?
426 if(Math::tenPi*minDistOfApproachSquared 237 if(Math::tenPi*minDistOfApproachSquared > totalCrossSection) return NULL;
427 238
428 // Bomb out if the two collision partner << 239 // Warn if the two collision partners are the same particle
429 // assert(p1->getID() != p2->getID()); << 240 if(p1->getID() == p2->getID()) {
>> 241 ERROR("At BinaryCollisonAvatar generation, ID1 (" << p1->getID()
>> 242 << ") == ID2 (" << p2->getID() <<")." << std::endl);
>> 243 }
430 244
431 // Return a new avatar, then! 245 // Return a new avatar, then!
432 return new G4INCL::BinaryCollisionAvatar 246 return new G4INCL::BinaryCollisionAvatar(t, totalCrossSection, theNucleus, p1, p2);
433 } 247 }
434 248
435 G4double StandardPropagationModel::getRefl 249 G4double StandardPropagationModel::getReflectionTime(G4INCL::Particle const * const aParticle) {
436 Intersection theIntersection( << 250 G4double time = 0.0;
437 IntersectionFactory::getLaterTraject << 251 const G4double T2 = aParticle->getMomentum().mag2();
438 aParticle->getPosition(), << 252 const G4double T4 = aParticle->getPosition().mag2();
439 aParticle->getPropagationVelocity( << 253
440 theNucleus->getSurfaceRadius(aPart << 254 const G4double r = theNucleus->getSurfaceRadius(aParticle);
441 G4double time; << 255 const G4double T1 = aParticle->getMomentum().dot(aParticle->getPosition());
442 if(theIntersection.exists) { << 256 const G4double T3 = T1/T2;
443 time = currentTime + theIntersection.t << 257 const G4double T5 = T3*T3 + (r*r-T4)/T2;
444 } else { << 258 if(T5 < 0.0) {
445 INCL_ERROR("Imaginary reflection time << 259 ERROR("Imaginary reflection time! Delta = " << T5 << " for particle: " << std::endl
446 << aParticle->print()); << 260 << aParticle->prG4int());
447 time = 10000.0; 261 time = 10000.0;
>> 262 } else {
>> 263 time = currentTime + (-T3 + std::sqrt(T5)) * aParticle->getEnergy();
448 } 264 }
449 return time; 265 return time;
450 } 266 }
451 267
452 G4double StandardPropagationModel::getTime 268 G4double StandardPropagationModel::getTime(G4INCL::Particle const * const particleA,
453 G4INCL::Particle const * const 269 G4INCL::Particle const * const particleB, G4double *minDistOfApproach) const
454 { 270 {
455 G4double time; 271 G4double time;
456 G4INCL::ThreeVector t13 = particleA->get << 272 G4INCL::ThreeVector t13 = particleA->getMomentum()/particleA->getEnergy();
457 t13 -= particleB->getPropagationVelocity << 273 t13 -= particleB->getMomentum()/particleB->getEnergy();
458 G4INCL::ThreeVector distance = particleA 274 G4INCL::ThreeVector distance = particleA->getPosition();
459 distance -= particleB->getPosition(); 275 distance -= particleB->getPosition();
460 const G4double t7 = t13.dot(distance); 276 const G4double t7 = t13.dot(distance);
461 const G4double dt = t13.mag2(); 277 const G4double dt = t13.mag2();
462 if(dt <= 1.0e-10) { 278 if(dt <= 1.0e-10) {
463 (*minDistOfApproach) = 100000.0; 279 (*minDistOfApproach) = 100000.0;
464 return currentTime + 100000.0; 280 return currentTime + 100000.0;
465 } else { 281 } else {
466 time = -t7/dt; 282 time = -t7/dt;
467 } 283 }
468 (*minDistOfApproach) = distance.mag2() + 284 (*minDistOfApproach) = distance.mag2() + time * t7;
469 return currentTime + time; 285 return currentTime + time;
470 } 286 }
471 287
>> 288 void StandardPropagationModel::checkCollisions(const ParticleList &participants,
>> 289 const ParticleList &particles)
>> 290 {
>> 291 G4int iind = 1, jind = 1;
>> 292 for(ParticleIter i = participants.begin(); i != participants.end(); ++i) {
>> 293 jind = 0;
>> 294 for(ParticleIter j = particles.begin(); j != particles.end(); ++j) {
>> 295 if( (*j)->isParticipant() )
>> 296 {
>> 297 ++jind;
>> 298 if(jind >= iind) continue;
>> 299 }
>> 300 if((*i)->getID() == (*j)->getID()) continue; // Do not process the collision of a particle with itself
>> 301
>> 302 registerAvatar(generateBinaryCollisionAvatar(*i,*j));
>> 303 }
>> 304 ++iind;
>> 305 }
>> 306 }
>> 307
472 void StandardPropagationModel::generateUpd 308 void StandardPropagationModel::generateUpdatedCollisions(const ParticleList &updatedParticles, const ParticleList &particles) {
473 309
474 // Loop over all the updated particles 310 // Loop over all the updated particles
475 for(ParticleIter updated=updatedParticle << 311 for(ParticleIter updated = updatedParticles.begin(); updated != updatedParticles.end(); ++updated)
476 { 312 {
477 // Loop over all the particles 313 // Loop over all the particles
478 for(ParticleIter particle=particles.be << 314 for(ParticleIter particle = particles.begin(); particle != particles.end(); ++particle)
479 { 315 {
480 /* Consider the generation of a coll 316 /* Consider the generation of a collision avatar only if (*particle)
481 * is not one of the updated particl 317 * is not one of the updated particles.
482 * The criterion makes sure that you 318 * The criterion makes sure that you don't generate avatars between
483 * updated particles. */ 319 * updated particles. */
484 if(updatedParticles.contains(*partic << 320 if((*particle)->isInList(updatedParticles)) continue;
485 321
486 registerAvatar(generateBinaryCollisi 322 registerAvatar(generateBinaryCollisionAvatar(*particle,*updated));
487 } 323 }
488 } 324 }
489 } 325 }
490 326
491 void StandardPropagationModel::generateCol <<
492 // Loop over all the particles <<
493 for(ParticleIter p1=particles.begin(), e <<
494 // Loop over the rest of the particles <<
495 for(ParticleIter p2 = p1 + 1; p2 != pa <<
496 registerAvatar(generateBinaryCollisi <<
497 } <<
498 } <<
499 } <<
500 <<
501 void StandardPropagationModel::generateCol 327 void StandardPropagationModel::generateCollisions(const ParticleList &particles, const ParticleList &except) {
502 328
503 const G4bool haveExcept = !except.empty( << 329 G4bool haveExcept;
>> 330 haveExcept=(except.size()!=0);
504 331
505 // Loop over all the particles 332 // Loop over all the particles
506 for(ParticleIter p1=particles.begin(), e << 333 for(ParticleIter p1 = particles.begin(); p1 != particles.end(); ++p1)
507 { 334 {
508 // Loop over the rest of the particles 335 // Loop over the rest of the particles
509 ParticleIter p2 = p1; 336 ParticleIter p2 = p1;
510 for(++p2; p2 != particles.end(); ++p2) 337 for(++p2; p2 != particles.end(); ++p2)
511 { 338 {
512 // Skip the collision if both partic 339 // Skip the collision if both particles must be excluded
513 if(haveExcept && except.contains(*p1 << 340 if(haveExcept && (*p1)->isInList(except) && (*p2)->isInList(except)) continue;
514 341
515 registerAvatar(generateBinaryCollisi 342 registerAvatar(generateBinaryCollisionAvatar(*p1,*p2));
516 } 343 }
517 } 344 }
518 345
519 } 346 }
520 347
521 void StandardPropagationModel::updateAvata 348 void StandardPropagationModel::updateAvatars(const ParticleList &particles) {
522 349
523 for(ParticleIter iter=particles.begin(), << 350 for(ParticleIter iter = particles.begin(); iter != particles.end(); ++iter) {
524 G4double time = this->getReflectionTim 351 G4double time = this->getReflectionTime(*iter);
525 if(time <= maximumTime) registerAvatar 352 if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*iter, time, theNucleus));
526 } 353 }
527 ParticleList const &p = theNucleus->getS << 354 ParticleList p = theNucleus->getStore()->getParticles();
528 generateUpdatedCollisions(particles, p); 355 generateUpdatedCollisions(particles, p); // Predict collisions with spectators and participants
529 } 356 }
530 357
531 void StandardPropagationModel::generateAll << 358 void StandardPropagationModel::generateAllAvatars(G4bool excludeUpdated) {
532 ParticleList const &particles = theNucle << 359 ParticleList particles = theNucleus->getStore()->getParticles();
533 // assert(!particles.empty()); << 360 if(particles.empty()) { ERROR("No particles inside the nucleus!" << std::endl); }
534 G4double time; << 361 for(ParticleIter i = particles.begin(); i != particles.end(); ++i) {
535 for(ParticleIter i=particles.begin(), e= <<
536 time = this->getReflectionTime(*i); <<
537 if(time <= maximumTime) registerAvatar <<
538 } <<
539 generateCollisions(particles); <<
540 generateDecays(particles); <<
541 } <<
542 <<
543 #ifdef INCL_REGENERATE_AVATARS <<
544 void StandardPropagationModel::generateAll <<
545 ParticleList const &particles = theNucle <<
546 // assert(!particles.empty()); <<
547 for(ParticleIter i=particles.begin(), e= <<
548 G4double time = this->getReflectionTim 362 G4double time = this->getReflectionTime(*i);
549 if(time <= maximumTime) registerAvatar 363 if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
550 } 364 }
551 ParticleList except = fs->getModifiedPar << 365 ParticleList except;
552 ParticleList const &entering = fs->getEn << 366 if(excludeUpdated)
553 except.insert(except.end(), entering.beg << 367 except = theNucleus->getUpdatedParticles();
554 generateCollisions(particles,except); 368 generateCollisions(particles,except);
555 generateDecays(particles); 369 generateDecays(particles);
556 } 370 }
557 #endif <<
558 371
559 void StandardPropagationModel::generateDec 372 void StandardPropagationModel::generateDecays(const ParticleList &particles) {
560 for(ParticleIter i=particles.begin(), e= << 373 for(ParticleIter i = particles.begin(); i != particles.end(); ++i) {
561 if((*i)->isDelta()) { << 374 if((*i)->isDelta()) {
562 G4double decayTime = DeltaDecayCh << 375 G4double decayTime = DeltaDecayChannel::computeDecayTime((*i));
563 G4double time = currentTime + decay << 376 G4double time = currentTime + decayTime;
564 if(time <= maximumTime) { << 377 if(time <= maximumTime) {
565 registerAvatar(new DecayAvatar((* << 378 registerAvatar(new DecayAvatar((*i), time, theNucleus));
566 } << 379 }
567 } << 380 }
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 } 381 }
583 } 382 }
584 383
585 G4INCL::IAvatar* StandardPropagationModel: << 384 G4INCL::IAvatar* StandardPropagationModel::propagate()
586 { 385 {
587 if(fs) { << 386 if(firstAvatar) { // When we propagate particles for the first time we create the full list of avatars.
588 // We update only the information rela << 387 generateAllAvatars();
>> 388 /* if(!theNucleus->getStore()->containsCollisions()) {
>> 389 theNucleus->forceTransparent();
>> 390 return NULL;
>> 391 }*/
>> 392 firstAvatar = false;
>> 393 } else { // For subsequent avatars we update only the
>> 394 // information related to particles that were updated
589 // by the previous avatar. 395 // by the previous avatar.
590 #ifdef INCL_REGENERATE_AVATARS 396 #ifdef INCL_REGENERATE_AVATARS
591 #warning "The INCL_REGENERATE_AVATARS code has 397 #warning "The INCL_REGENERATE_AVATARS code has not been tested in a while. Use it at your peril."
592 if(!fs->getModifiedParticles().empty() << 398 // Regenerates the entire avatar list, skipping collisions between
593 // Regenerates the entire avatar lis << 399 // updated particles
594 // updated particles << 400 if(theNucleus->getUpdatedParticles().size()!=0 || theNucleus->getCreatedParticles().size()!=0) {
595 theNucleus->getStore()->clearAvatars 401 theNucleus->getStore()->clearAvatars();
596 theNucleus->getStore()->initialisePa 402 theNucleus->getStore()->initialiseParticleAvatarConnections();
597 generateAllAvatarsExceptUpdated(fs); << 403 generateAllAvatars(true);
598 } 404 }
599 #else 405 #else
600 ParticleList const &updatedParticles = << 406 // Deltas are created by transforming nucleon G4into a delta for
601 if(fs->getValidity()==PauliBlockedFS) << 407 // efficiency reasons
602 // This final state might represents << 408 Particle * const blockedDelta = theNucleus->getBlockedDelta();
603 // decay << 409 ParticleList updatedParticles = theNucleus->getUpdatedParticles();
604 // assert(updatedParticles.empty() || (updated << 410 if(blockedDelta)
605 // assert(fs->getEnteringParticles().empty() & << 411 updatedParticles.push_back(blockedDelta);
606 generateDecays(updatedParticles); << 412 generateDecays(updatedParticles);
607 } else { << 413
608 ParticleList const &entering = fs->g << 414 ParticleList needNewAvatars = theNucleus->getUpdatedParticles();
609 generateDecays(updatedParticles); << 415 ParticleList created = theNucleus->getCreatedParticles();
610 generateDecays(entering); << 416 needNewAvatars.splice(needNewAvatars.end(), created);
611 << 417 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 418 #endif
623 } 419 }
624 420
625 G4INCL::IAvatar *theAvatar = theNucleus- 421 G4INCL::IAvatar *theAvatar = theNucleus->getStore()->findSmallestTime();
626 if(theAvatar == 0) return 0; // Avatar l 422 if(theAvatar == 0) return 0; // Avatar list is empty
627 // theAvatar->dispose(); 423 // theAvatar->dispose();
628 424
629 if(theAvatar->getTime() < currentTime) { << 425 theNucleus->getStore()->timeStep(theAvatar->getTime() - currentTime);
630 INCL_ERROR("Avatar time = " << theAvat <<
631 return 0; <<
632 } else if(theAvatar->getTime() > current <<
633 theNucleus->getStore()->timeStep(theAv <<
634 426
>> 427 if(theAvatar->getTime() <= currentTime) {
>> 428 ERROR("Avatar time = " << theAvatar->getTime() << ", currentTime = " << currentTime << std::endl);
>> 429 return 0;
>> 430 } else {
635 currentTime = theAvatar->getTime(); 431 currentTime = theAvatar->getTime();
636 theNucleus->getStore()->getBook().setC << 432 theNucleus->getStore()->getBook()->setCurrentTime(currentTime);
637 } 433 }
638 434
639 return theAvatar; 435 return theAvatar;
640 } 436 }
641 } 437 }
642 438