Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLInteractionAvatar.cc
Go to the documentation of this file.
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  violationEFunctor(NULL)
70  {
71  }
72 
74  G4INCL::Particle *p2)
75  : IAvatar(time), theNucleus(n),
76  particle1(p1), particle2(p2),
77  isPiN((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon())),
78  violationEFunctor(NULL)
79  {
80  }
81 
83  }
84 
86  delete backupParticle1;
87  if(backupParticle2)
88  delete backupParticle2;
89  backupParticle1 = NULL;
90  backupParticle2 = NULL;
91  }
92 
94  if(backupParticle1)
95  (*backupParticle1) = (*particle1);
96  else
98 
99  if(particle2) {
100  if(backupParticle2)
101  (*backupParticle2) = (*particle2);
102  else
104 
108  } else {
110  }
111  }
112 
114  if(!theNucleus || p->isPion()) return; // Local energy does not make any sense without a nucleus
115 
118  }
119 
122 
124 
125  if(particle2) {
129  } else {
131  }
133  }
134 
136  if(!theNucleus)
137  return false;
138 
139  ThreeVector pos = p->getPosition();
140  p->rpCorrelate();
141  G4double pos2 = pos.mag2();
142  const G4double r = theNucleus->getSurfaceRadius(p);
143  short iterations=0;
144  const short maxIterations=50;
145 
146  if(pos2 < r*r) return true;
147 
148  while( pos2 >= r*r && iterations<maxIterations ) /* Loop checking, 10.07.2015, D.Mancusi */
149  {
150  pos *= std::sqrt(r*r*0.9801/pos2); // 0.9801 == 0.99*0.99
151  pos2 = pos.mag2();
152  iterations++;
153  }
154  if( iterations < maxIterations)
155  {
156  INCL_DEBUG("Particle position vector length was : " << p->getPosition().mag() << ", rescaled to: " << pos.mag() << '\n');
157  p->setPosition(pos);
158  return true;
159  }
160  else
161  return false;
162  }
163 
165  INCL_DEBUG("postInteraction: final state: " << '\n' << fs->print() << '\n');
169  modifiedAndCreated.insert(modifiedAndCreated.end(), created.begin(), created.end());
170 
171  // Boost back to lab
173 
174  // If there is no Nucleus, just return
175  if(!theNucleus) return;
176 
177  // Mark pions that have been created outside their well (we will force them
178  // to be emitted later).
179  for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
180  if((*i)->isPion() && (*i)->getPosition().mag() > theNucleus->getSurfaceRadius(*i)) {
181  (*i)->makeParticipant();
182  (*i)->setOutOfWell();
183  fs->addOutgoingParticle(*i);
184  INCL_DEBUG("Pion was created outside its potential well." << '\n'
185  << (*i)->print());
186  }
187 
188  // Try to enforce energy conservation
190  G4bool success = enforceEnergyConservation(fs);
191  if(!success) {
192  INCL_DEBUG("Enforcing energy conservation: failed!" << '\n');
193 
194  // Restore the state of the initial particles
196 
197  // Delete newly created particles
198  for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
199  delete *i;
200 
201  fs->reset();
204 
205  return; // Interaction is blocked. Return an empty final state.
206  }
207  INCL_DEBUG("Enforcing energy conservation: success!" << '\n');
208 
209  INCL_DEBUG("postInteraction after energy conservation: final state: " << '\n' << fs->print() << '\n');
210 
211  // Check that outgoing delta resonances can decay to pi-N
212  for(ParticleIter i=modified.begin(), e=modified.end(); i!=e; ++i )
213  if((*i)->isDelta() &&
214  (*i)->getMass() < ParticleTable::minDeltaMass) {
215  INCL_DEBUG("Mass of the produced delta below decay threshold; forbidding collision. deltaMass=" <<
216  (*i)->getMass() << '\n');
217 
218  // Restore the state of the initial particles
220 
221  // Delete newly created particles
222  for(ParticleIter j=created.begin(), end=created.end(); j!=end; ++j )
223  delete *j;
224 
225  fs->reset();
228 
229  return; // Interaction is blocked. Return an empty final state.
230  }
231 
232  INCL_DEBUG("Random seeds before Pauli blocking: " << Random::getSeeds() << '\n');
233  // Test Pauli blocking
235 
236  if(isBlocked) {
237  INCL_DEBUG("Pauli: Blocked!" << '\n');
238 
239  // Restore the state of the initial particles
241 
242  // Delete newly created particles
243  for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
244  delete *i;
245 
246  fs->reset();
247  fs->makePauliBlocked();
249 
250  return; // Interaction is blocked. Return an empty final state.
251  }
252  INCL_DEBUG("Pauli: Allowed!" << '\n');
253 
254  // Test CDPP blocking
256 
257  if(isCDPPBlocked) {
258  INCL_DEBUG("CDPP: Blocked!" << '\n');
259 
260  // Restore the state of the initial particles
262 
263  // Delete newly created particles
264  for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i )
265  delete *i;
266 
267  fs->reset();
268  fs->makePauliBlocked();
270 
271  return; // Interaction is blocked. Return an empty final state.
272  }
273  INCL_DEBUG("CDPP: Allowed!" << '\n');
274 
275  // If all went well, try to bring particles inside the nucleus...
276  for(ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i )
277  {
278  // ...except for pions beyond their surface radius.
279  if((*i)->isOutOfWell()) continue;
280 
281  const G4bool successBringParticlesInside = bringParticleInside(*i);
282  if( !successBringParticlesInside ) {
283  INCL_ERROR("Failed to bring particle inside the nucleus!" << '\n');
284  }
285  }
286 
287  // Collision accepted!
288  for(ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i ) {
289  if(!(*i)->isOutOfWell()) {
290  // Decide if the particle should be made into a spectator
291  // (Back to spectator)
292  G4bool goesBackToSpectator = false;
293  if((*i)->isNucleon() && theNucleus->getStore()->getConfig()->getBackToSpectator()) {
294  G4double threshold = (*i)->getPotentialEnergy();
295  if((*i)->getType()==Proton)
297  if((*i)->getKineticEnergy() < threshold)
298  goesBackToSpectator = true;
299  }
300 
301  // Thaw the particle propagation
302  (*i)->thawPropagation();
303 
304  // Increment or decrement the participant counters
305  if(goesBackToSpectator) {
306  INCL_DEBUG("The following particle goes back to spectator:" << '\n'
307  << (*i)->print() << '\n');
308  if(!(*i)->isTargetSpectator()) {
310  }
311  (*i)->makeTargetSpectator();
312  } else {
313  if((*i)->isTargetSpectator()) {
315  }
316  (*i)->makeParticipant();
317  }
318  }
319  }
320  ParticleList destroyed = fs->getDestroyedParticles();
321  for(ParticleIter i=destroyed.begin(), e=destroyed.end(); i!=e; ++i )
322  if(!(*i)->isTargetSpectator())
324 
325  return;
326  }
327 
329  (*particle1) = (*backupParticle1);
330  if(particle2)
331  (*particle2) = (*backupParticle2);
332  }
333 
335  if(!theNucleus) return false;
336  LocalEnergyType theLocalEnergyType;
337  if(getType()==DecayAvatarType || isPiN)
338  theLocalEnergyType = theNucleus->getStore()->getConfig()->getLocalEnergyPiType();
339  else
340  theLocalEnergyType = theNucleus->getStore()->getConfig()->getLocalEnergyBBType();
341 
342  const G4bool firstAvatar = (theNucleus->getStore()->getBook().getAcceptedCollisions() == 0);
343  return ((theLocalEnergyType == FirstCollisionLocalEnergy && firstAvatar) ||
344  theLocalEnergyType == AlwaysLocalEnergy);
345  }
346 
348  // Set up the violationE calculation
349  const G4bool manyBodyFinalState = (modifiedAndCreated.size() > 1);
350  if(manyBodyFinalState)
351  violationEFunctor = new ViolationEMomentumFunctor(theNucleus, modifiedAndCreated, fs->getTotalEnergyBeforeInteraction(), boostVector, shouldUseLocalEnergy());
352  else {
353  Particle * const p = modified.front();
354  // The following condition is necessary for the functor to work
355  // correctly. A similar condition exists in INCL4.6.
357  return false;
358  violationEFunctor = new ViolationEEnergyFunctor(theNucleus, p, fs->getTotalEnergyBeforeInteraction(), shouldUseLocalEnergy());
359  }
360 
361  // Apply the root-finding algorithm
362  const RootFinder::Solution theSolution = RootFinder::solve(violationEFunctor, 1.0);
363  if(theSolution.success) { // Apply the solution
364  (*violationEFunctor)(theSolution.x);
365  } else if(theNucleus){
366  INCL_DEBUG("Couldn't enforce energy conservation after an interaction, root-finding algorithm failed." << '\n');
368  }
369  delete violationEFunctor;
370  violationEFunctor = NULL;
371  return theSolution.success;
372  }
373 
374  /* *** ***
375  * *** InteractionAvatar::ViolationEMomentumFunctor methods ***
376  * *** ***/
377 
378  InteractionAvatar::ViolationEMomentumFunctor::ViolationEMomentumFunctor(Nucleus * const nucleus, ParticleList const &modAndCre, const G4double totalEnergyBeforeInteraction, ThreeVector const &boost, const G4bool localE) :
379  RootFunctor(0., 1E6),
380  finalParticles(modAndCre),
381  initialEnergy(totalEnergyBeforeInteraction),
382  theNucleus(nucleus),
383  boostVector(boost),
384  shouldUseLocalEnergy(localE)
385  {
386  // Store the particle momenta (necessary for the calls to
387  // scaleParticleMomenta() to work)
388  for(ParticleIter i=finalParticles.begin(), e=finalParticles.end(); i!=e; ++i) {
389  (*i)->boost(boostVector);
390  particleMomenta.push_back((*i)->getMomentum());
391  }
392  }
393 
394  InteractionAvatar::ViolationEMomentumFunctor::~ViolationEMomentumFunctor() {
395  particleMomenta.clear();
396  }
397 
398  G4double InteractionAvatar::ViolationEMomentumFunctor::operator()(const G4double alpha) const {
399  scaleParticleMomenta(alpha);
400 
401  G4double deltaE = 0.0;
402  for(ParticleIter i=finalParticles.begin(), e=finalParticles.end(); i!=e; ++i)
403  deltaE += (*i)->getEnergy() - (*i)->getPotentialEnergy();
404  deltaE -= initialEnergy;
405  return deltaE;
406  }
407 
408  void InteractionAvatar::ViolationEMomentumFunctor::scaleParticleMomenta(const G4double alpha) const {
409 
410  std::vector<ThreeVector>::const_iterator iP = particleMomenta.begin();
411  for(ParticleIter i=finalParticles.begin(), e=finalParticles.end(); i!=e; ++i, ++iP) {
412  (*i)->setMomentum((*iP)*alpha);
413  (*i)->adjustEnergyFromMomentum();
414  (*i)->rpCorrelate();
415  (*i)->boost(-boostVector);
416  if(theNucleus)
417  theNucleus->updatePotentialEnergy(*i);
418  else
419  (*i)->setPotentialEnergy(0.);
420 
421 //jcd if(shouldUseLocalEnergy && !(*i)->isPion()) { // This translates AECSVT's loops 1, 3 and 4
422  if(shouldUseLocalEnergy && !(*i)->isPion() && !(*i)->isEta() && !(*i)->isOmega()) { // This translates AECSVT's loops 1, 3 and 4
423 // assert(theNucleus); // Local energy without a nucleus doesn't make sense
424  const G4double energy = (*i)->getEnergy(); // Store the energy of the particle
425  G4double locE = KinematicsUtils::getLocalEnergy(theNucleus, *i); // Initial value of local energy
426  G4double locEOld;
427  G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
428  for(G4int iterLocE=0;
430  ++iterLocE) {
431  locEOld = locE;
432  (*i)->setEnergy(energy + locE); // Update the energy of the particle...
433  (*i)->adjustMomentumFromEnergy();
434  theNucleus->updatePotentialEnergy(*i); // ...update its potential energy...
435  locE = KinematicsUtils::getLocalEnergy(theNucleus, *i); // ...and recompute locE.
436  deltaLocE = std::abs(locE-locEOld);
437  }
438  }
439  }
440  }
441 
442  void InteractionAvatar::ViolationEMomentumFunctor::cleanUp(const G4bool success) const {
443  if(!success)
444  scaleParticleMomenta(1.);
445  }
446 
447  /* *** ***
448  * *** InteractionAvatar::ViolationEEnergyFunctor methods ***
449  * *** ***/
450 
451  InteractionAvatar::ViolationEEnergyFunctor::ViolationEEnergyFunctor(Nucleus * const nucleus, Particle * const aParticle, const G4double totalEnergyBeforeInteraction, const G4bool localE) :
452  RootFunctor(0., 1E6),
453  initialEnergy(totalEnergyBeforeInteraction),
454  theNucleus(nucleus),
455  theParticle(aParticle),
456  theEnergy(theParticle->getEnergy()),
457  theMomentum(theParticle->getMomentum()),
458  energyThreshold(KinematicsUtils::energy(theMomentum,ParticleTable::minDeltaMass)),
459  shouldUseLocalEnergy(localE)
460  {
461 // assert(theParticle->isDelta());
462  }
463 
464  G4double InteractionAvatar::ViolationEEnergyFunctor::operator()(const G4double alpha) const {
465  setParticleEnergy(alpha);
466  return theParticle->getEnergy() - theParticle->getPotentialEnergy() - initialEnergy;
467  }
468 
469  void InteractionAvatar::ViolationEEnergyFunctor::setParticleEnergy(const G4double alpha) const {
470 
471  G4double locE;
472  if(shouldUseLocalEnergy) {
473 // assert(theNucleus); // Local energy without a nucleus doesn't make sense
474  locE = KinematicsUtils::getLocalEnergy(theNucleus, theParticle); // Initial value of local energy
475  } else
476  locE = 0.;
477  G4double locEOld;
478  G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3;
479  for(G4int iterLocE=0;
481  ++iterLocE) {
482  locEOld = locE;
483  G4double particleEnergy = energyThreshold + locE + alpha*(theEnergy-energyThreshold);
484  const G4double theMass2 = std::pow(particleEnergy,2.)-theMomentum.mag2();
485  G4double theMass;
486  if(theMass2>ParticleTable::minDeltaMass2)
487  theMass = std::sqrt(theMass2);
488  else {
489  theMass = ParticleTable::minDeltaMass;
490  particleEnergy = energyThreshold;
491  }
492  theParticle->setMass(theMass);
493  theParticle->setEnergy(particleEnergy); // Update the energy of the particle...
494  if(theNucleus) {
495  theNucleus->updatePotentialEnergy(theParticle); // ...update its potential energy...
496  if(shouldUseLocalEnergy)
497  locE = KinematicsUtils::getLocalEnergy(theNucleus, theParticle); // ...and recompute locE.
498  else
499  locE = 0.;
500  } else
501  locE = 0.;
502  deltaLocE = std::abs(locE-locEOld);
503  }
504 
505  }
506 
507  void InteractionAvatar::ViolationEEnergyFunctor::cleanUp(const G4bool success) const {
508  if(!success)
509  setParticleEnergy(1.);
510  }
511 
512 }
G4bool isBlocked(ParticleList const &p, Nucleus const *const n)
Check Pauli blocking.
static const G4double locEAccuracy
Target accuracy in the determination of the local-energy Q-value.
G4bool isCDPPBlocked(ParticleList const &p, Nucleus const *const n)
Check CDPP blocking.
G4double getMass() const
Get the cached particle mass.
void transformToLocalEnergyFrame(Nucleus const *const n, Particle *const p)
G4int getAcceptedCollisions() const
Definition: G4INCLBook.hh:100
void preInteractionLocalEnergy(Particle *const p)
Apply local-energy transformation, if appropriate.
G4ThreadLocal G4double minDeltaMass2
G4bool shouldUseLocalEnergy() const
true if the given avatar should use local energy
ParticleList const & getModifiedParticles() const
static const G4int maxIterLocE
Max number of iterations for the determination of the local-energy Q-value.
const char * p
Definition: xmltok.h:285
#define INCL_ERROR(x)
void boost(const ThreeVector &aBoostVector)
const G4INCL::ThreeVector & getMomentum() const
Config const * getConfig()
Definition: G4INCLStore.hh:273
Store * getStore() const
G4double getEnergy() const
G4double getSurfaceRadius(Particle const *const particle) const
Get the maximum allowed radius for a given particle.
#define G4ThreadLocal
Definition: tls.hh:89
ThreeVector makeBoostVector(Particle const *const p1, Particle const *const p2)
int G4int
Definition: G4Types.hh:78
G4double mag2() const
G4double getPotentialEnergy() const
Get the particle potential energy.
static G4ThreadLocal Particle * backupParticle2
void rpCorrelate()
Make the particle follow a strict r-p correlation.
InteractionAvatar(G4double, G4INCL::Nucleus *, G4INCL::Particle *)
SeedVector getSeeds()
Definition: G4INCLRandom.cc:89
ParticleList const & getCreatedParticles() const
void boost(const ThreeVector &b) const
static void deleteBackupParticles()
Release the memory allocated for the backup particles.
Book & getBook()
Definition: G4INCLStore.hh:259
G4bool getBackToSpectator() const
Get back-to-spectator.
bool G4bool
Definition: G4Types.hh:79
virtual void setPosition(const G4INCL::ThreeVector &position)
void addOutgoingParticle(Particle *p)
void incrementEnergyViolationInteraction()
Definition: G4INCLBook.hh:80
void preInteractionBlocking()
Store the state of the particles before the interaction.
static G4ThreadLocal Particle * backupParticle1
const G4int n
void setTotalEnergyBeforeInteraction(G4double E)
G4double getTotalEnergyBeforeInteraction() const
void incrementCascading()
Definition: G4INCLBook.hh:77
void decrementCascading()
Definition: G4INCLBook.hh:78
G4double total(Particle const *const p1, Particle const *const p2)
std::string print() const
G4double getTransmissionBarrier(Particle const *const p)
Get the transmission barrier.
G4ThreadLocal G4double minDeltaMass
LocalEnergyType getLocalEnergyPiType() const
Get the type of local energy for pi-N and decay avatars.
G4double energy(const ThreeVector &p, const G4double m)
LocalEnergyType getLocalEnergyBBType() const
Get the type of local energy for N-N avatars.
ParticleList const & getDestroyedParticles() const
const G4INCL::ThreeVector & getPosition() const
Solution solve(RootFunctor const *const f, const G4double x0)
Numerically solve a one-dimensional equation.
G4bool enforceEnergyConservation(FinalState *const fs)
Enforce energy conservation.
void restoreParticles() const
Restore the state of both particles.
const G4double twoThirds
G4double mag() const
double G4double
Definition: G4Types.hh:76
G4bool bringParticleInside(Particle *const p)
#define INCL_DEBUG(x)
static const G4double alpha
G4bool isPion() const
Is this a pion?
static const G4double pos
AvatarType getType() const
ParticleList::const_iterator ParticleIter
Static root-finder algorithm.
G4double getLocalEnergy(Nucleus const *const n, Particle *const p)