ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLInteractionAvatar.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLInteractionAvatar.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 /* \file G4INCLInteractionAvatar.cc
39  * \brief Virtual class for interaction avatars.
40  *
41  * This class is inherited by decay and collision avatars. The goal is to
42  * provide a uniform treatment of common physics, such as Pauli blocking,
43  * enforcement of energy conservation, etc.
44  *
45  * \date Mar 1st, 2011
46  * \author Davide Mancusi
47  */
48 
50 #include "G4INCLKinematicsUtils.hh"
51 #include "G4INCLCrossSections.hh"
52 #include "G4INCLPauliBlocking.hh"
53 #include "G4INCLRootFinder.hh"
54 #include "G4INCLLogger.hh"
55 #include "G4INCLConfigEnums.hh"
56 // #include <cassert>
57 
58 namespace G4INCL {
59 
64 
66  : IAvatar(time), theNucleus(n),
67  particle1(p1), particle2(NULL),
68  isPiN(false),
69  weight(1.),
70  violationEFunctor(NULL)
71  {
72  }
73 
75  G4INCL::Particle *p2)
76  : IAvatar(time), theNucleus(n),
77  particle1(p1), particle2(p2),
78  isPiN((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon())),
79  weight(1.),
80  violationEFunctor(NULL)
81  {
82  }
83 
85  }
86 
88  delete backupParticle1;
89  if(backupParticle2)
90  delete backupParticle2;
91  backupParticle1 = NULL;
92  backupParticle2 = NULL;
93  }
94 
96  if(backupParticle1)
97  (*backupParticle1) = (*particle1);
98  else
100 
101  if(particle2) {
102  if(backupParticle2)
103  (*backupParticle2) = (*particle2);
104  else
106 
110  } else {
112  }
113  }
114 
116  if(!theNucleus || p->isMeson()) return; // Local energy does not make any sense without a nucleus
117 
120  }
121 
124 
126 
127  if(particle2) {
131  } else {
133  }
135  }
136 
138  if(!theNucleus)
139  return false;
140 
141  ThreeVector pos = p->getPosition();
142  p->rpCorrelate();
143  G4double pos2 = pos.mag2();
145  short iterations=0;
146  const short maxIterations=50;
147 
148  if(pos2 < r*r) return true;
149 
150  while( pos2 >= r*r && iterations<maxIterations ) /* Loop checking, 10.07.2015, D.Mancusi */
151  {
152  pos *= std::sqrt(r*r*0.9801/pos2); // 0.9801 == 0.99*0.99
153  pos2 = pos.mag2();
154  iterations++;
155  }
156  if( iterations < maxIterations)
157  {
158  INCL_DEBUG("Particle position vector length was : " << p->getPosition().mag() << ", rescaled to: " << pos.mag() << '\n');
159  p->setPosition(pos);
160  return true;
161  }
162  else
163  return false;
164  }
165 
167  INCL_DEBUG("postInteraction: final state: " << '\n' << fs->print() << '\n');
172  modifiedAndCreated.insert(modifiedAndCreated.end(), created.begin(), created.end());
174  ModifiedAndDestroyed.insert(ModifiedAndDestroyed.end(), Destroyed.begin(), Destroyed.end());
175 
176  // Boost back to lab
178 
179  // If there is no Nucleus, just return
180  if(!theNucleus) return;
181 
182  // Mark pions and kaons that have been created outside their well (we will force them
183  // to be emitted later).
184  for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
185  if(((*i)->isPion() || (*i)->isKaon() || (*i)->isAntiKaon()) && (*i)->getPosition().mag() > theNucleus->getSurfaceRadius(*i)) {
186  (*i)->makeParticipant();
187  (*i)->setOutOfWell();
188  fs->addOutgoingParticle(*i);
189  INCL_DEBUG("Pion was created outside its potential well." << '\n'
190  << (*i)->print());
191  }
192 
193  // Try to enforce energy conservation
195  G4bool success = enforceEnergyConservation(fs);
196  if(!success) {
197  INCL_DEBUG("Enforcing energy conservation: failed!" << '\n');
198 
199  // Restore the state of the initial particles
201 
202  // Delete newly created particles
203  for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
204  delete *i;
205 
206  fs->reset();
209 
210  return; // Interaction is blocked. Return an empty final state.
211  }
212  INCL_DEBUG("Enforcing energy conservation: success!" << '\n');
213 
214  INCL_DEBUG("postInteraction after energy conservation: final state: " << '\n' << fs->print() << '\n');
215 
216  // Check that outgoing delta resonances can decay to pi-N
217  for(ParticleIter i=modified.begin(), e=modified.end(); i!=e; ++i )
218  if((*i)->isDelta() &&
219  (*i)->getMass() < ParticleTable::minDeltaMass) {
220  INCL_DEBUG("Mass of the produced delta below decay threshold; forbidding collision. deltaMass=" <<
221  (*i)->getMass() << '\n');
222 
223  // Restore the state of the initial particles
225 
226  // Delete newly created particles
227  for(ParticleIter j=created.begin(), end=created.end(); j!=end; ++j )
228  delete *j;
229 
230  fs->reset();
233 
234  return; // Interaction is blocked. Return an empty final state.
235  }
236 
237  INCL_DEBUG("Random seeds before Pauli blocking: " << Random::getSeeds() << '\n');
238  // Test Pauli blocking
240 
241  if(isBlocked) {
242  INCL_DEBUG("Pauli: Blocked!" << '\n');
243 
244  // Restore the state of the initial particles
246 
247  // Delete newly created particles
248  for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
249  delete *i;
250 
251  fs->reset();
252  fs->makePauliBlocked();
254 
255  return; // Interaction is blocked. Return an empty final state.
256  }
257  INCL_DEBUG("Pauli: Allowed!" << '\n');
258 
259  // Test CDPP blocking
261 
262  if(isCDPPBlocked) {
263  INCL_DEBUG("CDPP: Blocked!" << '\n');
264 
265  // Restore the state of the initial particles
267 
268  // Delete newly created particles
269  for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
270  delete *i;
271 
272  fs->reset();
273  fs->makePauliBlocked();
275 
276  return; // Interaction is blocked. Return an empty final state.
277  }
278  INCL_DEBUG("CDPP: Allowed!" << '\n');
279 
280  // If all went well, try to bring particles inside the nucleus...
281  for(ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i )
282  {
283  // ...except for pions beyond their surface radius.
284  if((*i)->isOutOfWell()) continue;
285 
286  const G4bool successBringParticlesInside = bringParticleInside(*i);
287  if( !successBringParticlesInside ) {
288  INCL_ERROR("Failed to bring particle inside the nucleus!" << '\n');
289  }
290  }
291 
292  // Collision accepted!
293  // Biasing of the final state
294  std::vector<G4int> newBiasCollisionVector;
295  newBiasCollisionVector = ModifiedAndDestroyed.getParticleListBiasVector();
296  if(std::fabs(weight-1.) > 1E-6){
297  newBiasCollisionVector.push_back(Particle::nextBiasedCollisionID);
299  weight = 1.; // useless?
300  }
301  for(ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i ) {
302  (*i)->setBiasCollisionVector(newBiasCollisionVector);
303  if(!(*i)->isOutOfWell()) {
304  // Decide if the particle should be made into a spectator
305  // (Back to spectator)
306  G4bool goesBackToSpectator = false;
307  if((*i)->isNucleon() && theNucleus->getStore()->getConfig()->getBackToSpectator()) {
308  G4double threshold = (*i)->getPotentialEnergy();
309  if((*i)->getType()==Proton)
311  if((*i)->getKineticEnergy() < threshold)
312  goesBackToSpectator = true;
313  }
314  // Thaw the particle propagation
315  (*i)->thawPropagation();
316 
317  // Increment or decrement the participant counters
318  if(goesBackToSpectator) {
319  INCL_DEBUG("The following particle goes back to spectator:" << '\n'
320  << (*i)->print() << '\n');
321  if(!(*i)->isTargetSpectator()) {
323  }
324  (*i)->makeTargetSpectator();
325  } else {
326  if((*i)->isTargetSpectator()) {
328  }
329  (*i)->makeParticipant();
330  }
331  }
332  }
333  ParticleList destroyed = fs->getDestroyedParticles();
334  for(ParticleIter i=destroyed.begin(), e=destroyed.end(); i!=e; ++i )
335  if(!(*i)->isTargetSpectator())
337  return;
338  }
339 
341  (*particle1) = (*backupParticle1);
342  if(particle2)
343  (*particle2) = (*backupParticle2);
344  }
345 
347  if(!theNucleus) return false;
348  LocalEnergyType theLocalEnergyType;
349  if(getType()==DecayAvatarType || isPiN)
350  theLocalEnergyType = theNucleus->getStore()->getConfig()->getLocalEnergyPiType();
351  else
352  theLocalEnergyType = theNucleus->getStore()->getConfig()->getLocalEnergyBBType();
353 
354  const G4bool firstAvatar = (theNucleus->getStore()->getBook().getAcceptedCollisions() == 0);
355  return ((theLocalEnergyType == FirstCollisionLocalEnergy && firstAvatar) ||
356  theLocalEnergyType == AlwaysLocalEnergy);
357  }
358 
360  // Set up the violationE calculation
361  const G4bool manyBodyFinalState = (modifiedAndCreated.size() > 1);
362 
363  if(manyBodyFinalState)
365  else {
366  Particle * const p = modified.front();
367  // The following condition is necessary for the functor to work
368  // correctly. A similar condition exists in INCL4.6.
370  return false;
372  }
373 
374  // Apply the root-finding algorithm
375  const RootFinder::Solution theSolution = RootFinder::solve(violationEFunctor, 1.0);
376  if(theSolution.success) { // Apply the solution
377  (*violationEFunctor)(theSolution.x);
378  } else if(theNucleus){
379  INCL_DEBUG("Couldn't enforce energy conservation after an interaction, root-finding algorithm failed." << '\n');
381  }
382  delete violationEFunctor;
383  violationEFunctor = NULL;
384  return theSolution.success;
385  }
386 
387  /* *** ***
388  * *** InteractionAvatar::ViolationEMomentumFunctor methods ***
389  * *** ***/
390 
391  InteractionAvatar::ViolationEMomentumFunctor::ViolationEMomentumFunctor(Nucleus * const nucleus, ParticleList const &modAndCre, const G4double totalEnergyBeforeInteraction, ThreeVector const &boost, const G4bool localE) :
392  RootFunctor(0., 1E6),
393  finalParticles(modAndCre),
394  initialEnergy(totalEnergyBeforeInteraction),
395  theNucleus(nucleus),
396  boostVector(boost),
397  shouldUseLocalEnergy(localE)
398  {
399  // Store the particle momenta (necessary for the calls to
400  // scaleParticleMomenta() to work)
401  for(ParticleIter i=finalParticles.begin(), e=finalParticles.end(); i!=e; ++i) {
402  (*i)->boost(boostVector);
403  particleMomenta.push_back((*i)->getMomentum());
404  }
405  }
406 
408  particleMomenta.clear();
409  }
410 
412  scaleParticleMomenta(alpha);
413 
414  G4double deltaE = 0.0;
415  for(ParticleIter i=finalParticles.begin(), e=finalParticles.end(); i!=e; ++i)
416  deltaE += (*i)->getEnergy() - (*i)->getPotentialEnergy();
417  deltaE -= initialEnergy;
418  return deltaE;
419  }
420 
422 
423  std::vector<ThreeVector>::const_iterator iP = particleMomenta.begin();
424  for(ParticleIter i=finalParticles.begin(), e=finalParticles.end(); i!=e; ++i, ++iP) {
425  (*i)->setMomentum((*iP)*alpha);
426  (*i)->adjustEnergyFromMomentum();
427  (*i)->rpCorrelate();
428  (*i)->boost(-boostVector);
429  if(theNucleus)
431  else
432  (*i)->setPotentialEnergy(0.);
433 
434 //jcd if(shouldUseLocalEnergy && !(*i)->isPion()) { // This translates AECSVT's loops 1, 3 and 4
435  if(shouldUseLocalEnergy && !(*i)->isPion() && !(*i)->isEta() && !(*i)->isOmega() &&
436  !(*i)->isKaon() && !(*i)->isAntiKaon() && !(*i)->isSigma() && !(*i)->isLambda()) { // This translates AECSVT's loops 1, 3 and 4
437 // assert(theNucleus); // Local energy without a nucleus doesn't make sense
438  const G4double energy = (*i)->getEnergy(); // Store the energy of the particle
439  G4double locE = KinematicsUtils::getLocalEnergy(theNucleus, *i); // Initial value of local energy
440  G4double locEOld;
441  G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
442  for(G4int iterLocE=0;
444  ++iterLocE) {
445  locEOld = locE;
446  (*i)->setEnergy(energy + locE); // Update the energy of the particle...
447  (*i)->adjustMomentumFromEnergy();
448  theNucleus->updatePotentialEnergy(*i); // ...update its potential energy...
449  locE = KinematicsUtils::getLocalEnergy(theNucleus, *i); // ...and recompute locE.
450  deltaLocE = std::abs(locE-locEOld);
451  }
452  }
453 
454 //jlrs For lambdas and nuclei with masses higher than 19 also local energy
455  if(shouldUseLocalEnergy && (*i)->isLambda() && theNucleus->getA()>19) {
456 // assert(theNucleus); // Local energy without a nucleus doesn't make sense
457  const G4double energy = (*i)->getEnergy(); // Store the energy of the particle
458  G4double locE = KinematicsUtils::getLocalEnergy(theNucleus, *i); // Initial value of local energy
459  G4double locEOld;
460  G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
461  for(G4int iterLocE=0;
463  ++iterLocE) {
464  locEOld = locE;
465  (*i)->setEnergy(energy + locE); // Update the energy of the particle...
466  (*i)->adjustMomentumFromEnergy();
467  theNucleus->updatePotentialEnergy(*i); // ...update its potential energy...
468  locE = KinematicsUtils::getLocalEnergy(theNucleus, *i); // ...and recompute locE.
469  deltaLocE = std::abs(locE-locEOld);
470  }
471  }
472  }
473  }
474 
476  if(!success)
477  scaleParticleMomenta(1.);
478  }
479 
480  /* *** ***
481  * *** InteractionAvatar::ViolationEEnergyFunctor methods ***
482  * *** ***/
483 
484  InteractionAvatar::ViolationEEnergyFunctor::ViolationEEnergyFunctor(Nucleus * const nucleus, Particle * const aParticle, const G4double totalEnergyBeforeInteraction, const G4bool localE) :
485  RootFunctor(0., 1E6),
486  initialEnergy(totalEnergyBeforeInteraction),
487  theNucleus(nucleus),
488  theParticle(aParticle),
489  theEnergy(theParticle->getEnergy()),
490  theMomentum(theParticle->getMomentum()),
491  energyThreshold(KinematicsUtils::energy(theMomentum,ParticleTable::minDeltaMass)),
492  shouldUseLocalEnergy(localE)
493  {
494 // assert(theParticle->isDelta());
495  }
496 
498  setParticleEnergy(alpha);
499  return theParticle->getEnergy() - theParticle->getPotentialEnergy() - initialEnergy;
500  }
501 
503 
504  G4double locE;
506 // assert(theNucleus); // Local energy without a nucleus doesn't make sense
507  locE = KinematicsUtils::getLocalEnergy(theNucleus, theParticle); // Initial value of local energy
508  } else
509  locE = 0.;
510  G4double locEOld;
511  G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
512  for(G4int iterLocE=0;
514  ++iterLocE) {
515  locEOld = locE;
516  G4double particleEnergy = energyThreshold + locE + alpha*(theEnergy-energyThreshold);
517  const G4double theMass2 = std::pow(particleEnergy,2.)-theMomentum.mag2();
518  G4double theMass;
519  if(theMass2>ParticleTable::minDeltaMass2)
520  theMass = std::sqrt(theMass2);
521  else {
522  theMass = ParticleTable::minDeltaMass;
523  particleEnergy = energyThreshold;
524  }
525  theParticle->setMass(theMass);
526  theParticle->setEnergy(particleEnergy); // Update the energy of the particle...
527  if(theNucleus) {
528  theNucleus->updatePotentialEnergy(theParticle); // ...update its potential energy...
530  locE = KinematicsUtils::getLocalEnergy(theNucleus, theParticle); // ...and recompute locE.
531  else
532  locE = 0.;
533  } else
534  locE = 0.;
535  deltaLocE = std::abs(locE-locEOld);
536  }
537 
538  }
539 
541  if(!success)
542  setParticleEnergy(1.);
543  }
544 
545 }