| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /*
39 * StandardPropagationModel.cpp
40 *
41 * \date 4 juin 2009
42 * \author Pekka Kaitaniemi
43 */
44
45 #include "G4INCLStandardPropagationModel.hh"
46 #include "G4INCLPbarAtrestEntryChannel.hh"
47 #include "G4INCLSurfaceAvatar.hh"
48 #include "G4INCLBinaryCollisionAvatar.hh"
49 #include "G4INCLDecayAvatar.hh"
50 #include "G4INCLCrossSections.hh"
51 #include "G4INCLRandom.hh"
52 #include <iostream>
53 #include <algorithm>
54 #include "G4INCLLogger.hh"
55 #include "G4INCLGlobals.hh"
56 #include "G4INCLKinematicsUtils.hh"
57 #include "G4INCLCoulombDistortion.hh"
58 #include "G4INCLDeltaDecayChannel.hh"
59 #include "G4INCLSigmaZeroDecayChannel.hh"
60 #include "G4INCLPionResonanceDecayChannel.hh"
61 #include "G4INCLParticleEntryAvatar.hh"
62 #include "G4INCLIntersection.hh"
63 #include <vector>
64
65 namespace G4INCL {
66
67 StandardPropagationModel::StandardPropagat
68 :theNucleus(0), maximumTime(70.0), curre
69 hadronizationTime(hTime),
70 firstAvatar(true),
71 theLocalEnergyType(localEnergyType),
72 theLocalEnergyDeltaType(localEnergyDelta
73 {
74 }
75
76 StandardPropagationModel::~StandardPropaga
77 {
78 delete theNucleus;
79 }
80
81 G4INCL::Nucleus* StandardPropagationModel:
82 {
83 return theNucleus;
84 }
85
86 //D
87
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;
184
185 // Create the projectile particle
186 const G4double projectileMass = Particle
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
199 TLab = p->getKineticEnergy()/p->getA()
200 }
201
202 // energy-dependent stopping time above
203 if(TLab>2000.)
204 temfin *= (5.8E4-TLab)/5.6E4;
205
206 maximumTime = temfin;
207
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
218 // parameter larger than the maximum imp
219 // account Coulomb distortion.
220 if(impactParameter>CoulombDistortion::ma
221 INCL_DEBUG("impactParameter>CoulombDis
222 delete p;
223 return -1.;
224 }
225
226 ThreeVector position(impactParameter * s
227 impactParameter * std::sin(phi),
228 0.);
229 p->setPosition(position);
230
231 // Fill in the relevant kinematic variab
232 theNucleus->setIncomingAngularMomentum(p
233 theNucleus->setIncomingMomentum(p->getMo
234 theNucleus->setInitialEnergy(p->getEnerg
235 + ParticleTable::getTableMass(theNuc
236
237 // Reset the particle kinematics to the
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
251 return p->getTransversePosition().mag(
252 } else {
253 delete p;
254 return -1.;
255 }
256 }
257
258 G4double StandardPropagationModel::shootCo
259 theNucleus->setNucleusNucleusCollision()
260 currentTime = 0.0;
261
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 }
331
332 G4double StandardPropagationModel::getStop
333 return maximumTime;
334 }
335
336 void StandardPropagationModel::setStopping
337 // assert(time>0.0);
338 maximumTime = time;
339 }
340
341 G4double StandardPropagationModel::getCurr
342 return currentTime;
343 }
344
345 void StandardPropagationModel::setNucleus(
346 {
347 theNucleus = nucleus;
348 }
349
350 void StandardPropagationModel::registerAva
351 {
352 if(anAvatar) theNucleus->getStore()->add
353 }
354
355 IAvatar *StandardPropagationModel::generat
356 // Is either particle a participant?
357 if(!p1->isParticipant() && !p2->isPartic
358
359 // Is it a pi-resonance collision (we do
360 if((p1->isResonance() && p2->isPion()) |
361 return NULL;
362
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
368 G4double minDistOfApproachSquared = 0.0;
369 G4double t = getTime(p1, p2, &minDistOfA
370 if(t>maximumTime || t<currentTime+hadron
371
372 // Local energy. Jump through some hoops
373 // at the collision point, and clean up
374 G4bool hasLocalEnergy;
375 if(p1->isPion() || p2->isPion())
376 hasLocalEnergy = ((theLocalEnergyDelta
377 theNucleus->getStore()->getBook(
378 theLocalEnergyDeltaType == AlwaysL
379 else
380 hasLocalEnergy = ((theLocalEnergyType
381 theNucleus->getStore()->getBook(
382 theLocalEnergyType == AlwaysLocalE
383 const G4bool p1HasLocalEnergy = (hasLoca
384 const G4bool p2HasLocalEnergy = (hasLoca
385
386 if(p1HasLocalEnergy) {
387 backupParticle1 = *p1;
388 p1->propagate(t - currentTime);
389 if(p1->getPosition().mag() > theNucleu
390 *p1 = backupParticle1;
391 return NULL;
392 }
393 KinematicsUtils::transformToLocalEnerg
394 }
395 if(p2HasLocalEnergy) {
396 backupParticle2 = *p2;
397 p2->propagate(t - currentTime);
398 if(p2->getPosition().mag() > theNucleu
399 *p2 = backupParticle2;
400 if(p1HasLocalEnergy) {
401 *p1 = backupParticle1;
402 }
403 return NULL;
404 }
405 KinematicsUtils::transformToLocalEnerg
406 }
407
408 // Compute the total cross section
409 const G4double totalCrossSection = Cross
410 const G4double squareTotalEnergyInCM = K
411
412 // Restore particles to their state befo
413 if(p1HasLocalEnergy) {
414 *p1 = backupParticle1;
415 }
416 if(p2HasLocalEnergy) {
417 *p2 = backupParticle2;
418 }
419
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
426 if(Math::tenPi*minDistOfApproachSquared
427
428 // Bomb out if the two collision partner
429 // assert(p1->getID() != p2->getID());
430
431 // Return a new avatar, then!
432 return new G4INCL::BinaryCollisionAvatar
433 }
434
435 G4double StandardPropagationModel::getRefl
436 Intersection theIntersection(
437 IntersectionFactory::getLaterTraject
438 aParticle->getPosition(),
439 aParticle->getPropagationVelocity(
440 theNucleus->getSurfaceRadius(aPart
441 G4double time;
442 if(theIntersection.exists) {
443 time = currentTime + theIntersection.t
444 } else {
445 INCL_ERROR("Imaginary reflection time
446 << aParticle->print());
447 time = 10000.0;
448 }
449 return time;
450 }
451
452 G4double StandardPropagationModel::getTime
453 G4INCL::Particle const * const
454 {
455 G4double time;
456 G4INCL::ThreeVector t13 = particleA->get
457 t13 -= particleB->getPropagationVelocity
458 G4INCL::ThreeVector distance = particleA
459 distance -= particleB->getPosition();
460 const G4double t7 = t13.dot(distance);
461 const G4double dt = t13.mag2();
462 if(dt <= 1.0e-10) {
463 (*minDistOfApproach) = 100000.0;
464 return currentTime + 100000.0;
465 } else {
466 time = -t7/dt;
467 }
468 (*minDistOfApproach) = distance.mag2() +
469 return currentTime + time;
470 }
471
472 void StandardPropagationModel::generateUpd
473
474 // Loop over all the updated particles
475 for(ParticleIter updated=updatedParticle
476 {
477 // Loop over all the particles
478 for(ParticleIter particle=particles.be
479 {
480 /* Consider the generation of a coll
481 * is not one of the updated particl
482 * The criterion makes sure that you
483 * updated particles. */
484 if(updatedParticles.contains(*partic
485
486 registerAvatar(generateBinaryCollisi
487 }
488 }
489 }
490
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
502
503 const G4bool haveExcept = !except.empty(
504
505 // Loop over all the particles
506 for(ParticleIter p1=particles.begin(), e
507 {
508 // Loop over the rest of the particles
509 ParticleIter p2 = p1;
510 for(++p2; p2 != particles.end(); ++p2)
511 {
512 // Skip the collision if both partic
513 if(haveExcept && except.contains(*p1
514
515 registerAvatar(generateBinaryCollisi
516 }
517 }
518
519 }
520
521 void StandardPropagationModel::updateAvata
522
523 for(ParticleIter iter=particles.begin(),
524 G4double time = this->getReflectionTim
525 if(time <= maximumTime) registerAvatar
526 }
527 ParticleList const &p = theNucleus->getS
528 generateUpdatedCollisions(particles, p);
529 }
530
531 void StandardPropagationModel::generateAll
532 ParticleList const &particles = theNucle
533 // assert(!particles.empty());
534 G4double time;
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
549 if(time <= maximumTime) registerAvatar
550 }
551 ParticleList except = fs->getModifiedPar
552 ParticleList const &entering = fs->getEn
553 except.insert(except.end(), entering.beg
554 generateCollisions(particles,except);
555 generateDecays(particles);
556 }
557 #endif
558
559 void StandardPropagationModel::generateDec
560 for(ParticleIter i=particles.begin(), e=
561 if((*i)->isDelta()) {
562 G4double decayTime = DeltaDecayCh
563 G4double time = currentTime + decay
564 if(time <= maximumTime) {
565 registerAvatar(new DecayAvatar((*
566 }
567 }
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 }
583 }
584
585 G4INCL::IAvatar* StandardPropagationModel:
586 {
587 if(fs) {
588 // We update only the information rela
589 // by the previous avatar.
590 #ifdef INCL_REGENERATE_AVATARS
591 #warning "The INCL_REGENERATE_AVATARS code has
592 if(!fs->getModifiedParticles().empty()
593 // Regenerates the entire avatar lis
594 // updated particles
595 theNucleus->getStore()->clearAvatars
596 theNucleus->getStore()->initialisePa
597 generateAllAvatarsExceptUpdated(fs);
598 }
599 #else
600 ParticleList const &updatedParticles =
601 if(fs->getValidity()==PauliBlockedFS)
602 // This final state might represents
603 // decay
604 // assert(updatedParticles.empty() || (updated
605 // assert(fs->getEnteringParticles().empty() &
606 generateDecays(updatedParticles);
607 } else {
608 ParticleList const &entering = fs->g
609 generateDecays(updatedParticles);
610 generateDecays(entering);
611
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
623 }
624
625 G4INCL::IAvatar *theAvatar = theNucleus-
626 if(theAvatar == 0) return 0; // Avatar l
627 // theAvatar->dispose();
628
629 if(theAvatar->getTime() < currentTime) {
630 INCL_ERROR("Avatar time = " << theAvat
631 return 0;
632 } else if(theAvatar->getTime() > current
633 theNucleus->getStore()->timeStep(theAv
634
635 currentTime = theAvatar->getTime();
636 theNucleus->getStore()->getBook().setC
637 }
638
639 return theAvatar;
640 }
641 }
642