ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLStandardPropagationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLStandardPropagationModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
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 Helsinki Institute of Physics, Finland
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 
46 #include "G4INCLSurfaceAvatar.hh"
48 #include "G4INCLDecayAvatar.hh"
49 #include "G4INCLCrossSections.hh"
50 #include "G4INCLRandom.hh"
51 #include <iostream>
52 #include <algorithm>
53 #include "G4INCLLogger.hh"
54 #include "G4INCLGlobals.hh"
55 #include "G4INCLKinematicsUtils.hh"
61 #include "G4INCLIntersection.hh"
62 
63 namespace G4INCL {
64 
66  :theNucleus(0), maximumTime(70.0), currentTime(0.0),
67  hadronizationTime(hTime),
68  firstAvatar(true),
69  theLocalEnergyType(localEnergyType),
70  theLocalEnergyDeltaType(localEnergyDeltaType)
71  {
72  }
73 
75  {
76  delete theNucleus;
77  }
78 
80  {
81  return theNucleus;
82  }
83 
84  G4double StandardPropagationModel::shoot(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
85  if(projectileSpecies.theType==Composite)
86  return shootComposite(projectileSpecies, kineticEnergy, impactParameter, phi);
87  else
88  return shootParticle(projectileSpecies.theType, kineticEnergy, impactParameter, phi);
89  }
90 
91  G4double StandardPropagationModel::shootParticle(ParticleType const type, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
93  currentTime = 0.0;
94 
95  // Create the projectile particle
96  const G4double projectileMass = ParticleTable::getTableParticleMass(type);
97  G4double energy = kineticEnergy + projectileMass;
98  G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
99  ThreeVector momentum(0.0, 0.0, momentumZ);
100  Particle *p= new G4INCL::Particle(type, energy, momentum, ThreeVector());
101 
102  G4double temfin;
103  G4double TLab;
104  if( p->isMeson()) {
105  temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
106  TLab = p->getKineticEnergy();
107  } else {
108  temfin = 29.8 * std::pow(theNucleus->getA(), 0.16);
109  TLab = p->getKineticEnergy()/p->getA();
110  }
111 
112  // energy-dependent stopping time above 2 AGeV
113  if(TLab>2000.)
114  temfin *= (5.8E4-TLab)/5.6E4;
115 
116  maximumTime = temfin;
117 
118  // If the incoming particle is slow, use a larger stopping time
119  const G4double rMax = theNucleus->getUniverseRadius();
120  const G4double distance = 2.*rMax;
121  const G4double projectileVelocity = p->boostVector().mag();
122  const G4double traversalTime = distance / projectileVelocity;
123  if(maximumTime < traversalTime)
124  maximumTime = traversalTime;
125  INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
126 
127  // If Coulomb is activated, do not process events with impact
128  // parameter larger than the maximum impact parameter, taking into
129  // account Coulomb distortion.
130  if(impactParameter>CoulombDistortion::maxImpactParameter(p->getSpecies(), kineticEnergy, theNucleus)) {
131  INCL_DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << '\n');
132  delete p;
133  return -1.;
134  }
135 
136  ThreeVector position(impactParameter * std::cos(phi),
137  impactParameter * std::sin(phi),
138  0.);
139  p->setPosition(position);
140 
141  // Fill in the relevant kinematic variables
146 
147  // Reset the particle kinematics to the INCL values
148  p->setINCLMass();
149  p->setEnergy(p->getMass() + kineticEnergy);
151 
154  firstAvatar = false;
155 
156  // Get the entry avatars from Coulomb and put them in the Store
158  if(theEntryAvatar) {
159  theNucleus->getStore()->addParticleEntryAvatar(theEntryAvatar);
160 
161  return p->getTransversePosition().mag();
162  } else {
163  delete p;
164  return -1.;
165  }
166  }
167 
168  G4double StandardPropagationModel::shootComposite(ParticleSpecies const &species, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
170  currentTime = 0.0;
171 
172  // Create the ProjectileRemnant object
173  ProjectileRemnant *pr = new ProjectileRemnant(species, kineticEnergy);
174 
175  // Same stopping time as for nucleon-nucleus
176  maximumTime = 29.8 * std::pow(theNucleus->getA(), 0.16);
177 
178  // If the incoming cluster is slow, use a larger stopping time
179  const G4double rms = ParticleTable::getLargestNuclearRadius(pr->getA(), pr->getZ());
180  const G4double rMax = theNucleus->getUniverseRadius();
181  const G4double distance = 2.*rMax + 2.725*rms;
182  const G4double projectileVelocity = pr->boostVector().mag();
183  const G4double traversalTime = distance / projectileVelocity;
184  if(maximumTime < traversalTime)
185  maximumTime = traversalTime;
186  INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
187 
188  // If Coulomb is activated, do not process events with impact
189  // parameter larger than the maximum impact parameter, taking into
190  // account Coulomb distortion.
191  if(impactParameter>CoulombDistortion::maxImpactParameter(pr,theNucleus)) {
192  INCL_DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << '\n');
193  delete pr;
194  return -1.;
195  }
196 
197  // Position the cluster at the right impact parameter
198  ThreeVector position(impactParameter * std::cos(phi),
199  impactParameter * std::sin(phi),
200  0.);
201  pr->setPosition(position);
202 
203  // Fill in the relevant kinematic variables
208 
210  firstAvatar = false;
211 
212  // Get the entry avatars from Coulomb
213  IAvatarList theAvatarList
215 
216  if(theAvatarList.empty()) {
217  INCL_DEBUG("No ParticleEntryAvatar found, transparent event" << '\n');
218  delete pr;
219  return -1.;
220  }
221 
222  /* Store the internal kinematics of the projectile remnant.
223  *
224  * Note that this is at variance with the Fortran version, which stores
225  * the initial kinematics of the particles *after* putting the spectators
226  * on mass shell, but *before* removing the same energy from the
227  * participants. Due to the different code flow, doing so in the C++
228  * version leads to wrong excitation energies for the forced compound
229  * nucleus.
230  */
231  pr->storeComponents();
232 
233  // Tell the Nucleus about the ProjectileRemnant
235 
236  // Register the ParticleEntryAvatars
237  theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
238 
239  return pr->getTransversePosition().mag();
240  }
241 
243  return maximumTime;
244  }
245 
247 // assert(time>0.0);
248  maximumTime = time;
249  }
250 
252  return currentTime;
253  }
254 
256  {
257  theNucleus = nucleus;
258  }
259 
261  {
262  if(anAvatar) theNucleus->getStore()->add(anAvatar);
263  }
264 
266  // Is either particle a participant?
267  if(!p1->isParticipant() && !p2->isParticipant() && p1->getParticipantType()==p2->getParticipantType()) return NULL;
268 
269  // Is it a pi-resonance collision (we don't treat them)?
270  if((p1->isResonance() && p2->isPion()) || (p1->isPion() && p2->isResonance()))
271  return NULL;
272 
273  // Will the avatar take place between now and the end of the cascade?
274  G4double minDistOfApproachSquared = 0.0;
275  G4double t = getTime(p1, p2, &minDistOfApproachSquared);
276  if(t>maximumTime || t<currentTime+hadronizationTime) return NULL;
277 
278  // Local energy. Jump through some hoops to calculate the cross section
279  // at the collision point, and clean up after yourself afterwards.
280  G4bool hasLocalEnergy;
281  if(p1->isPion() || p2->isPion())
282  hasLocalEnergy = ((theLocalEnergyDeltaType == FirstCollisionLocalEnergy &&
285  else
286  hasLocalEnergy = ((theLocalEnergyType == FirstCollisionLocalEnergy &&
289  const G4bool p1HasLocalEnergy = (hasLocalEnergy && !p1->isMeson());
290  const G4bool p2HasLocalEnergy = (hasLocalEnergy && !p2->isMeson());
291 
292  if(p1HasLocalEnergy) {
293  backupParticle1 = *p1;
294  p1->propagate(t - currentTime);
295  if(p1->getPosition().mag() > theNucleus->getSurfaceRadius(p1)) {
296  *p1 = backupParticle1;
297  return NULL;
298  }
300  }
301  if(p2HasLocalEnergy) {
302  backupParticle2 = *p2;
303  p2->propagate(t - currentTime);
304  if(p2->getPosition().mag() > theNucleus->getSurfaceRadius(p2)) {
305  *p2 = backupParticle2;
306  if(p1HasLocalEnergy) {
307  *p1 = backupParticle1;
308  }
309  return NULL;
310  }
312  }
313 
314  // Compute the total cross section
315  const G4double totalCrossSection = CrossSections::total(p1, p2);
317 
318  // Restore particles to their state before the local-energy tweak
319  if(p1HasLocalEnergy) {
320  *p1 = backupParticle1;
321  }
322  if(p2HasLocalEnergy) {
323  *p2 = backupParticle2;
324  }
325 
326  // Is the CM energy > cutNN? (no cutNN on the first collision)
328  && p1->isNucleon() && p2->isNucleon()
329  && squareTotalEnergyInCM < BinaryCollisionAvatar::getCutNNSquared()) return NULL;
330 
331  // Do the particles come close enough to each other?
332  if(Math::tenPi*minDistOfApproachSquared > totalCrossSection) return NULL;
333 
334  // Bomb out if the two collision partners are the same particle
335 // assert(p1->getID() != p2->getID());
336 
337  // Return a new avatar, then!
338  return new G4INCL::BinaryCollisionAvatar(t, totalCrossSection, theNucleus, p1, p2);
339  }
340 
342  Intersection theIntersection(
344  aParticle->getPosition(),
345  aParticle->getPropagationVelocity(),
346  theNucleus->getSurfaceRadius(aParticle)));
347  G4double time;
348  if(theIntersection.exists) {
349  time = currentTime + theIntersection.time;
350  } else {
351  INCL_ERROR("Imaginary reflection time for particle: " << '\n'
352  << aParticle->print());
353  time = 10000.0;
354  }
355  return time;
356  }
357 
359  G4INCL::Particle const * const particleB, G4double *minDistOfApproach) const
360  {
361  G4double time;
362  G4INCL::ThreeVector t13 = particleA->getPropagationVelocity();
363  t13 -= particleB->getPropagationVelocity();
364  G4INCL::ThreeVector distance = particleA->getPosition();
365  distance -= particleB->getPosition();
366  const G4double t7 = t13.dot(distance);
367  const G4double dt = t13.mag2();
368  if(dt <= 1.0e-10) {
369  (*minDistOfApproach) = 100000.0;
370  return currentTime + 100000.0;
371  } else {
372  time = -t7/dt;
373  }
374  (*minDistOfApproach) = distance.mag2() + time * t7;
375  return currentTime + time;
376  }
377 
378  void StandardPropagationModel::generateUpdatedCollisions(const ParticleList &updatedParticles, const ParticleList &particles) {
379 
380  // Loop over all the updated particles
381  for(ParticleIter updated=updatedParticles.begin(), e=updatedParticles.end(); updated!=e; ++updated)
382  {
383  // Loop over all the particles
384  for(ParticleIter particle=particles.begin(), end=particles.end(); particle!=end; ++particle)
385  {
386  /* Consider the generation of a collision avatar only if (*particle)
387  * is not one of the updated particles.
388  * The criterion makes sure that you don't generate avatars between
389  * updated particles. */
390  if(updatedParticles.contains(*particle)) continue;
391 
393  }
394  }
395  }
396 
398  // Loop over all the particles
399  for(ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1) {
400  // Loop over the rest of the particles
401  for(ParticleIter p2 = p1 + 1; p2 != particles.end(); ++p2) {
403  }
404  }
405  }
406 
408 
409  const G4bool haveExcept = !except.empty();
410 
411  // Loop over all the particles
412  for(ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1)
413  {
414  // Loop over the rest of the particles
415  ParticleIter p2 = p1;
416  for(++p2; p2 != particles.end(); ++p2)
417  {
418  // Skip the collision if both particles must be excluded
419  if(haveExcept && except.contains(*p1) && except.contains(*p2)) continue;
420 
422  }
423  }
424 
425  }
426 
428 
429  for(ParticleIter iter=particles.begin(), e=particles.end(); iter!=e; ++iter) {
430  G4double time = this->getReflectionTime(*iter);
431  if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*iter, time, theNucleus));
432  }
434  generateUpdatedCollisions(particles, p); // Predict collisions with spectators and participants
435  }
436 
438  ParticleList const &particles = theNucleus->getStore()->getParticles();
439 // assert(!particles.empty());
440  G4double time;
441  for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
442  time = this->getReflectionTime(*i);
443  if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
444  }
445  generateCollisions(particles);
446  generateDecays(particles);
447  }
448 
449 #ifdef INCL_REGENERATE_AVATARS
450  void StandardPropagationModel::generateAllAvatarsExceptUpdated(FinalState const * const fs) {
451  ParticleList const &particles = theNucleus->getStore()->getParticles();
452 // assert(!particles.empty());
453  for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
454  G4double time = this->getReflectionTime(*i);
455  if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
456  }
457  ParticleList except = fs->getModifiedParticles();
458  ParticleList const &entering = fs->getEnteringParticles();
459  except.insert(except.end(), entering.begin(), entering.end());
460  generateCollisions(particles,except);
461  generateDecays(particles);
462  }
463 #endif
464 
466  for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
467  if((*i)->isDelta()) {
468  G4double decayTime = DeltaDecayChannel::computeDecayTime((*i)); // time in fm/c
469  G4double time = currentTime + decayTime;
470  if(time <= maximumTime) {
471  registerAvatar(new DecayAvatar((*i), time, theNucleus));
472  }
473  }
474  else if((*i)->getType() == SigmaZero) {
475  G4double decayTime = SigmaZeroDecayChannel::computeDecayTime((*i)); // time in fm/c
476  G4double time = currentTime + decayTime;
477  if(time <= maximumTime) {
478  registerAvatar(new DecayAvatar((*i), time, theNucleus));
479  }
480  }
481  if((*i)->isOmega()) {
483  G4double timeOmega = currentTime + decayTimeOmega;
484  if(timeOmega <= maximumTime) {
485  registerAvatar(new DecayAvatar((*i), timeOmega, theNucleus));
486  }
487  }
488  }
489  }
490 
492  {
493  if(fs) {
494  // We update only the information related to particles that were updated
495  // by the previous avatar.
496 #ifdef INCL_REGENERATE_AVATARS
497 #warning "The INCL_REGENERATE_AVATARS code has not been tested in a while. Use it at your peril."
498  if(!fs->getModifiedParticles().empty() || !fs->getEnteringParticles().empty() || !fs->getCreatedParticles().empty()) {
499  // Regenerates the entire avatar list, skipping collisions between
500  // updated particles
502  theNucleus->getStore()->initialiseParticleAvatarConnections();
503  generateAllAvatarsExceptUpdated(fs);
504  }
505 #else
506  ParticleList const &updatedParticles = fs->getModifiedParticles();
507  if(fs->getValidity()==PauliBlockedFS) {
508  // This final state might represents the outcome of a Pauli-blocked delta
509  // decay
510 // assert(updatedParticles.empty() || (updatedParticles.size()==1 && updatedParticles.front()->isResonance()));
511 // assert(fs->getEnteringParticles().empty() && fs->getCreatedParticles().empty() && fs->getOutgoingParticles().empty() && fs->getDestroyedParticles().empty());
512  generateDecays(updatedParticles);
513  } else {
514  ParticleList const &entering = fs->getEnteringParticles();
515  generateDecays(updatedParticles);
516  generateDecays(entering);
517 
518  ParticleList const &created = fs->getCreatedParticles();
519  if(created.empty() && entering.empty())
520  updateAvatars(updatedParticles);
521  else {
522  ParticleList updatedParticlesCopy = updatedParticles;
523  updatedParticlesCopy.insert(updatedParticlesCopy.end(), entering.begin(), entering.end());
524  updatedParticlesCopy.insert(updatedParticlesCopy.end(), created.begin(), created.end());
525  updateAvatars(updatedParticlesCopy);
526  }
527  }
528 #endif
529  }
530 
532  if(theAvatar == 0) return 0; // Avatar list is empty
533  // theAvatar->dispose();
534 
535  if(theAvatar->getTime() < currentTime) {
536  INCL_ERROR("Avatar time = " << theAvatar->getTime() << ", currentTime = " << currentTime << '\n');
537  return 0;
538  } else if(theAvatar->getTime() > currentTime) {
539  theNucleus->getStore()->timeStep(theAvatar->getTime() - currentTime);
540 
541  currentTime = theAvatar->getTime();
543  }
544 
545  return theAvatar;
546  }
547 }