Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCascade.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 
42 #include "G4INCLCascade.hh"
43 #include "G4INCLRandom.hh"
45 #include "G4INCLParticleTable.hh"
47 #include "G4INCLGlobalInfo.hh"
48 
49 #include "G4INCLPauliBlocking.hh"
50 
51 #include "G4INCLCrossSections.hh"
52 
54 
55 #include "G4INCLLogger.hh"
56 #include "G4INCLGlobals.hh"
58 
60 
62 
63 #include "G4INCLClustering.hh"
64 
65 #include "G4INCLIntersection.hh"
66 
68 
69 #include "G4INCLCascadeAction.hh"
71 
72 #include <cstring>
73 #include <cstdlib>
74 #include <numeric>
75 
76 namespace G4INCL {
77 
78  INCL::INCL(Config const * const config)
79  :propagationModel(0), theA(208), theZ(82),
80  targetInitSuccess(false),
82  maxUniverseRadius(0.),
83  maxInteractionDistance(0.),
84  fixedImpactParameter(0.),
85  theConfig(config),
86  nucleus(NULL),
87  forceTransparent(false),
88  minRemnantSize(4)
89  {
90  // Set the logger object.
91 #ifdef INCLXX_IN_GEANT4_MODE
93 #else // INCLXX_IN_GEANT4_MODE
94  Logger::initialize(theConfig);
95 #endif // INCLXX_IN_GEANT4_MODE
96 
97  // Set the random number generator algorithm. The system can support
98  // multiple different generator algorithms in a completely
99  // transparent way.
100  Random::initialize(theConfig);
101 
102  // Select the Pauli and CDPP blocking algorithms
103  Pauli::initialize(theConfig);
104 
105  // Set the cross-section set
106  CrossSections::initialize(theConfig);
107 
108  // Set the phase-space generator
110 
111  // Select the Coulomb-distortion algorithm:
113 
114  // Select the clustering algorithm:
115  Clustering::initialize(theConfig);
116 
117  // Initialize the INCL particle table:
118  ParticleTable::initialize(theConfig);
119 
120  // Initialize the value of cutNN in BinaryCollisionAvatar
122 
123  // Propagation model is responsible for finding avatars and
124  // transporting the particles. In principle this step is "hidden"
125  // behind an abstract interface and the rest of the system does not
126  // care how the transportation and avatar finding is done. This
127  // should allow us to "easily" experiment with different avatar
128  // finding schemes and even to support things like curved
129  // trajectories in the future.
130  propagationModel = new StandardPropagationModel(theConfig->getLocalEnergyBBType(),theConfig->getLocalEnergyPiType(),theConfig->getHadronizationTime());
131  if(theConfig->getCascadeActionType() == AvatarDumpActionType)
132  cascadeAction = new AvatarDumpAction();
133  else
134  cascadeAction = new CascadeAction();
135  cascadeAction->beforeRunAction(theConfig);
136 
137  theGlobalInfo.cascadeModel = theConfig->getVersionString();
138  theGlobalInfo.deexcitationModel = theConfig->getDeExcitationString();
139 #ifdef INCL_ROOT_USE
140  theGlobalInfo.rootSelection = theConfig->getROOTSelectionString();
141 #endif
142 
143 #ifndef INCLXX_IN_GEANT4_MODE
144  // Fill in the global information
145  theGlobalInfo.At = theConfig->getTargetA();
146  theGlobalInfo.Zt = theConfig->getTargetZ();
147  const ParticleSpecies theSpecies = theConfig->getProjectileSpecies();
148  theGlobalInfo.Ap = theSpecies.theA;
149  theGlobalInfo.Zp = theSpecies.theZ;
150  theGlobalInfo.Ep = theConfig->getProjectileKineticEnergy();
151 #endif
152 
153  fixedImpactParameter = theConfig->getImpactParameter();
154  }
155 
158 #ifndef INCLXX_IN_GEANT4_MODE
159  NuclearMassTable::deleteTable();
160 #endif
167 #ifndef INCLXX_IN_GEANT4_MODE
168  Logger::deleteLoggerSlave();
169 #endif
172  cascadeAction->afterRunAction();
173  delete cascadeAction;
174  delete propagationModel;
175  delete theConfig;
176  }
177 
178  G4bool INCL::prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z) {
179  if(A < 0 || A > 300 || Z < 1 || Z > 200) {
180  INCL_ERROR("Unsupported target: A = " << A << " Z = " << Z << '\n'
181  << "Target configuration rejected." << '\n');
182  return false;
183  }
184  if(projectileSpecies.theType==Composite &&
185  (projectileSpecies.theZ==projectileSpecies.theA || projectileSpecies.theZ==0)) {
186  INCL_ERROR("Unsupported projectile: A = " << projectileSpecies.theA << " Z = " << projectileSpecies.theZ << '\n'
187  << "Projectile configuration rejected." << '\n');
188  return false;
189  }
190 
191  // Reset the forced-transparent flag
192  forceTransparent = false;
193 
194  // Initialise the maximum universe radius
195  initUniverseRadius(projectileSpecies, kineticEnergy, A, Z);
196 
197  // Initialise the nucleus
198  theZ = Z;
199  if(theConfig->isNaturalTarget())
201  else
202  theA = A;
203  initializeTarget(theA, theZ);
204 
205  // Set the maximum impact parameter
206  maxImpactParameter = CoulombDistortion::maxImpactParameter(projectileSpecies, kineticEnergy, nucleus);
207  INCL_DEBUG("Maximum impact parameter initialised: " << maxImpactParameter << '\n');
208 
209  // For forced CN events
210  initMaxInteractionDistance(projectileSpecies, kineticEnergy);
211 
212  // Set the geometric cross section
213  theGlobalInfo.geometricCrossSection =
214  Math::tenPi*std::pow(maxImpactParameter,2);
215 
216  // Set the minimum remnant size
217  if(projectileSpecies.theA > 0)
218  minRemnantSize = std::min(theA, 4);
219  else
220  minRemnantSize = std::min(theA-1, 4);
221 
222  return true;
223  }
224 
226  delete nucleus;
227 
228  nucleus = new Nucleus(A, Z, theConfig, maxUniverseRadius);
229  nucleus->getStore()->getBook().reset();
230  nucleus->initializeParticles();
231 
232  propagationModel->setNucleus(nucleus);
233  return true;
234  }
235 
237  ParticleSpecies const &projectileSpecies,
238  const G4double kineticEnergy,
239  const G4int targetA,
240  const G4int targetZ
241  ) {
242  // Set the target and the projectile
243  targetInitSuccess = prepareReaction(projectileSpecies, kineticEnergy, targetA, targetZ);
244 
245  if(!targetInitSuccess) {
246  INCL_WARN("Target initialisation failed for A=" << targetA << ", Z=" << targetZ << '\n');
247  theEventInfo.transparent=true;
248  return theEventInfo;
249  }
250 
251  cascadeAction->beforeCascadeAction(propagationModel);
252 
253  const G4bool canRunCascade = preCascade(projectileSpecies, kineticEnergy);
254  if(canRunCascade) {
255  cascade();
256  postCascade();
257  cascadeAction->afterCascadeAction(nucleus);
258  }
259  updateGlobalInfo();
260  return theEventInfo;
261  }
262 
263  G4bool INCL::preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
264  // Reset theEventInfo
265  theEventInfo.reset();
266 
268 
269  // Fill in the event information
270  theEventInfo.projectileType = projectileSpecies.theType;
271  theEventInfo.Ap = projectileSpecies.theA;
272  theEventInfo.Zp = projectileSpecies.theZ;
273  theEventInfo.Ep = kineticEnergy;
274  theEventInfo.At = nucleus->getA();
275  theEventInfo.Zt = nucleus->getZ();
276 
277  // Do nothing below the Coulomb barrier
278  if(maxImpactParameter<=0.) {
279  // Fill in the event information
280  theEventInfo.transparent = true;
281 
282  return false;
283  }
284 
285  // Randomly draw an impact parameter or use a fixed value, depending on the
286  // Config option
287  G4double impactParameter, phi;
288  if(fixedImpactParameter<0.) {
289  impactParameter = maxImpactParameter * std::sqrt(Random::shoot0());
290  phi = Random::shoot() * Math::twoPi;
291  } else {
292  impactParameter = fixedImpactParameter;
293  phi = 0.;
294  }
295  INCL_DEBUG("Selected impact parameter: " << impactParameter << '\n');
296 
297  // Fill in the event information
298  theEventInfo.impactParameter = impactParameter;
299 
300  const G4double effectiveImpactParameter = propagationModel->shoot(projectileSpecies, kineticEnergy, impactParameter, phi);
301  if(effectiveImpactParameter < 0.) {
302  // Fill in the event information
303  theEventInfo.transparent = true;
304 
305  return false;
306  }
307 
308  // Fill in the event information
309  theEventInfo.transparent = false;
310  theEventInfo.effectiveImpactParameter = effectiveImpactParameter;
311 
312  return true;
313  }
314 
315  void INCL::cascade() {
316  FinalState *finalState = new FinalState;
317 
318  unsigned long loopCounter = 0;
319  const unsigned long maxLoopCounter = 10000000;
320  do {
321  // Run book keeping actions that should take place before propagation:
322  cascadeAction->beforePropagationAction(propagationModel);
323 
324  // Get the avatar with the smallest time and propagate particles
325  // to that point in time.
326  IAvatar *avatar = propagationModel->propagate(finalState);
327 
328  finalState->reset();
329 
330  // Run book keeping actions that should take place after propagation:
331  cascadeAction->afterPropagationAction(propagationModel, avatar);
332 
333  if(avatar == 0) break; // No more avatars in the avatar list.
334 
335  // Run book keeping actions that should take place before avatar:
336  cascadeAction->beforeAvatarAction(avatar, nucleus);
337 
338  // Channel is responsible for calculating the outcome of the
339  // selected avatar. There are different kinds of channels. The
340  // class IChannel is, again, an abstract interface that defines
341  // the externally observable behavior of all interaction
342  // channels.
343  // The handling of the channel is transparent to the API.
344  // Final state tells what changed...
345  avatar->fillFinalState(finalState);
346 
347  // Run book keeping actions that should take place after avatar:
348  cascadeAction->afterAvatarAction(avatar, nucleus, finalState);
349 
350  // So now we must give this information to the nucleus
351  nucleus->applyFinalState(finalState);
352  // and now we are ready to process the next avatar!
353 
354  delete avatar;
355 
356  ++loopCounter;
357  } while(continueCascade() && loopCounter<maxLoopCounter); /* Loop checking, 10.07.2015, D.Mancusi */
358 
359  delete finalState;
360  }
361 
362  void INCL::postCascade() {
363  // Fill in the event information
364  theEventInfo.stoppingTime = propagationModel->getCurrentTime();
365 
366  // Forced CN?
367  if(nucleus->getTryCompoundNucleus()) {
368  INCL_DEBUG("Trying compound nucleus" << '\n');
369  makeCompoundNucleus();
370  theEventInfo.transparent = forceTransparent;
371  // Global checks of conservation laws
372 #ifndef INCLXX_IN_GEANT4_MODE
373  if(!theEventInfo.transparent) globalConservationChecks(true);
374 #endif
375  return;
376  }
377 
378  theEventInfo.transparent = forceTransparent || nucleus->isEventTransparent();
379 
380  if(theEventInfo.transparent) {
381  ProjectileRemnant * const projectileRemnant = nucleus->getProjectileRemnant();
382  if(projectileRemnant) {
383  // Clear the incoming list (particles will be deleted by the ProjectileRemnant)
384  nucleus->getStore()->clearIncoming();
385  } else {
386  // Delete particles in the incoming list
387  nucleus->getStore()->deleteIncoming();
388  }
389  } else {
390  // Check if the nucleus contains deltas
391  theEventInfo.deltasInside = nucleus->containsDeltas();
392 
393  // Take care of any remaining deltas
394  theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
395  theEventInfo.forcedDeltasInside = nucleus->decayInsideDeltas();
396 
397  // Take care of any remaining etas and/or omegas
398 /* if (theConfig->getDecayTimeThreshold() > 7.e-19) {
399  theEventInfo.forcedPionResonancesOutside = nucleus->decayOutgoingPionResonances();
400  }*/
401  G4double timeThreshold=theConfig->getDecayTimeThreshold();
402  theEventInfo.forcedPionResonancesOutside = nucleus->decayOutgoingPionResonances(timeThreshold);
403 
404  // Apply Coulomb distortion, if appropriate
405  // Note that this will apply Coulomb distortion also on pions emitted by
406  // unphysical remnants (see decayInsideDeltas). This is at variance with
407  // what INCL4.6 does, but these events are (should be!) so rare that
408  // whatever we do doesn't (shouldn't!) make any noticeable difference.
410 
411  // If the normal cascade predicted complete fusion, use the tabulated
412  // masses to compute the excitation energy, the recoil, etc.
413  if(nucleus->getStore()->getOutgoingParticles().size()==0
414  && (!nucleus->getProjectileRemnant()
415  || nucleus->getProjectileRemnant()->getParticles().size()==0)) {
416 
417  INCL_DEBUG("Cascade resulted in complete fusion, using realistic fusion kinematics" << '\n');
418 
419  nucleus->useFusionKinematics();
420 
421  if(nucleus->getExcitationEnergy()<0.) {
422  // Complete fusion is energetically impossible, return a transparent
423  INCL_WARN("Complete-fusion kinematics yields negative excitation energy, returning a transparent!" << '\n');
424  theEventInfo.transparent = true;
425  return;
426  }
427 
428  } else { // Normal cascade here
429 
430  // Set the excitation energy
431  nucleus->setExcitationEnergy(nucleus->computeExcitationEnergy());
432 
433  // Make a projectile pre-fragment out of the geometrical and dynamical
434  // spectators
435  theEventInfo.nUnmergedSpectators = makeProjectileRemnant();
436 
437  // Compute recoil momentum, energy and spin of the nucleus
438  if(nucleus->getA()==1 && minRemnantSize>1) {
439  INCL_ERROR("Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << '\n');
440  }
441  nucleus->computeRecoilKinematics();
442 
443 #ifndef INCLXX_IN_GEANT4_MODE
444  // Global checks of conservation laws
445  globalConservationChecks(false);
446 #endif
447 
448  // Make room for the remnant recoil by rescaling the energies of the
449  // outgoing particles.
450  if(nucleus->hasRemnant()) rescaleOutgoingForRecoil();
451 
452  }
453 
454  // Cluster decay
455  theEventInfo.clusterDecay = nucleus->decayOutgoingClusters() | nucleus->decayMe();
456 
457 #ifndef INCLXX_IN_GEANT4_MODE
458  // Global checks of conservation laws
459  globalConservationChecks(true);
460 #endif
461 
462  // Fill the EventInfo structure
463  nucleus->fillEventInfo(&theEventInfo);
464 
465  }
466  }
467 
468  void INCL::makeCompoundNucleus() {
469  // If this is not a nucleus-nucleus collision, don't attempt to make a
470  // compound nucleus.
471  //
472  // Yes, even nucleon-nucleus collisions can lead to particles entering
473  // below the Fermi level. Take e.g. 1-MeV p + He4.
474  if(!nucleus->isNucleusNucleusCollision()) {
475  forceTransparent = true;
476  return;
477  }
478 
479  // Reset the internal Nucleus variables
480  nucleus->getStore()->clearIncoming();
481  nucleus->getStore()->clearOutgoing();
482  nucleus->getProjectileRemnant()->reset();
483  nucleus->setA(theEventInfo.At);
484  nucleus->setZ(theEventInfo.Zt);
485 
486  // CN kinematical variables
487  // Note: the CN orbital angular momentum is neglected in what follows. We
488  // should actually take it into account!
489  ThreeVector theCNMomentum = nucleus->getIncomingMomentum();
490  ThreeVector theCNSpin = nucleus->getIncomingAngularMomentum();
491  const G4double theTargetMass = ParticleTable::getTableMass(theEventInfo.At, theEventInfo.Zt);
492  G4int theCNA=theEventInfo.At, theCNZ=theEventInfo.Zt;
493  Cluster * const theProjectileRemnant = nucleus->getProjectileRemnant();
494  G4double theCNEnergy = theTargetMass + theProjectileRemnant->getEnergy();
495 
496  // Loop over the potential participants
497  ParticleList const &initialProjectileComponents = theProjectileRemnant->getParticles();
498  std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
499  // Shuffle the list of potential participants
500  std::random_shuffle(shuffledComponents.begin(), shuffledComponents.end(), Random::getAdapter());
501 
502  G4bool success = true;
503  G4bool atLeastOneNucleonEntering = false;
504  for(std::vector<Particle*>::const_iterator p=shuffledComponents.begin(), e=shuffledComponents.end(); p!=e; ++p) {
505  // Skip particles that miss the interaction distance
507  (*p)->getPosition(),
508  (*p)->getPropagationVelocity(),
509  maxInteractionDistance));
510  if(!intersectionInteractionDistance.exists)
511  continue;
512 
513  // Build an entry avatar for this nucleon
514  atLeastOneNucleonEntering = true;
515  ParticleEntryAvatar *theAvatar = new ParticleEntryAvatar(0.0, nucleus, *p);
516  nucleus->getStore()->addParticleEntryAvatar(theAvatar);
517  FinalState *fs = theAvatar->getFinalState();
518  nucleus->applyFinalState(fs);
519  FinalStateValidity validity = fs->getValidity();
520  delete fs;
521  switch(validity) {
522  case ValidFS:
524  case ParticleBelowZeroFS:
525  // Add the particle to the CN
526  theCNA++;
527  theCNZ += (*p)->getZ();
528  break;
529  case PauliBlockedFS:
531  default:
532  success = false;
533  break;
534  }
535  }
536 
537  if(!success || !atLeastOneNucleonEntering) {
538  INCL_DEBUG("No nucleon entering in forced CN, forcing a transparent" << '\n');
539  forceTransparent = true;
540  return;
541  }
542 
543 // assert(theCNA==nucleus->getA());
544 // assert(theCNA<=theEventInfo.At+theEventInfo.Ap);
545 // assert(theCNZ<=theEventInfo.Zt+theEventInfo.Zp);
546 
547  // Update the kinematics of the CN
548  theCNEnergy -= theProjectileRemnant->getEnergy();
549  theCNMomentum -= theProjectileRemnant->getMomentum();
550 
551  // Deal with the projectile remnant
552  nucleus->finalizeProjectileRemnant(propagationModel->getCurrentTime());
553 
554  // Subtract the angular momentum of the projectile remnant
555 // assert(nucleus->getStore()->getOutgoingParticles().empty());
556  theCNSpin -= theProjectileRemnant->getAngularMomentum();
557 
558  // Compute the excitation energy of the CN
559  const G4double theCNMass = ParticleTable::getTableMass(theCNA,theCNZ);
560  const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.mag2();
561  if(theCNInvariantMassSquared<0.) {
562  // Negative invariant mass squared, return a transparent
563  forceTransparent = true;
564  return;
565  }
566  const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
567  if(theCNExcitationEnergy<0.) {
568  // Negative excitation energy, return a transparent
569  INCL_DEBUG("CN excitation energy is negative, forcing a transparent" << '\n'
570  << " theCNA = " << theCNA << '\n'
571  << " theCNZ = " << theCNZ << '\n'
572  << " theCNEnergy = " << theCNEnergy << '\n'
573  << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
574  << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
575  << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
576  );
577  forceTransparent = true;
578  return;
579  } else {
580  // Positive excitation energy, can make a CN
581  INCL_DEBUG("CN excitation energy is positive, forcing a CN" << '\n'
582  << " theCNA = " << theCNA << '\n'
583  << " theCNZ = " << theCNZ << '\n'
584  << " theCNEnergy = " << theCNEnergy << '\n'
585  << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
586  << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
587  << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
588  );
589  nucleus->setA(theCNA);
590  nucleus->setZ(theCNZ);
591  nucleus->setMomentum(theCNMomentum);
592  nucleus->setEnergy(theCNEnergy);
593  nucleus->setExcitationEnergy(theCNExcitationEnergy);
594  nucleus->setMass(theCNMass+theCNExcitationEnergy);
595  nucleus->setSpin(theCNSpin); // neglects any orbital angular momentum of the CN
596 
597  // Take care of any remaining deltas
598  theEventInfo.forcedDeltasOutside = nucleus->decayOutgoingDeltas();
599 
600  // Take care of any remaining etas and/or omegas
601  G4double timeThreshold=theConfig->getDecayTimeThreshold();
602  theEventInfo.forcedPionResonancesOutside = nucleus->decayOutgoingPionResonances(timeThreshold);
603 
604  // Cluster decay
605  theEventInfo.clusterDecay = nucleus->decayOutgoingClusters() | nucleus->decayMe();
606 
607  // Fill the EventInfo structure
608  nucleus->fillEventInfo(&theEventInfo);
609  }
610  }
611 
612  void INCL::rescaleOutgoingForRecoil() {
613  RecoilCMFunctor theRecoilFunctor(nucleus, theEventInfo);
614 
615  // Apply the root-finding algorithm
616  const RootFinder::Solution theSolution = RootFinder::solve(&theRecoilFunctor, 1.0);
617  if(theSolution.success) {
618  theRecoilFunctor(theSolution.x); // Apply the solution
619  } else {
620  INCL_WARN("Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." << '\n');
621  }
622 
623  }
624 
625 #ifndef INCLXX_IN_GEANT4_MODE
626  void INCL::globalConservationChecks(G4bool afterRecoil) {
627  Nucleus::ConservationBalance theBalance = nucleus->getConservationBalance(theEventInfo,afterRecoil);
628 
629  // Global conservation checks
630  const G4double pLongBalance = theBalance.momentum.getZ();
631  const G4double pTransBalance = theBalance.momentum.perp();
632  if(theBalance.Z != 0) {
633  INCL_ERROR("Violation of charge conservation! ZBalance = " << theBalance.Z << '\n');
634  }
635  if(theBalance.A != 0) {
636  INCL_ERROR("Violation of baryon-number conservation! ABalance = " << theBalance.A << '\n');
637  }
638  G4double EThreshold, pLongThreshold, pTransThreshold;
639  if(afterRecoil) {
640  // Less stringent checks after accommodating recoil
641  EThreshold = 10.; // MeV
642  pLongThreshold = 1.; // MeV/c
643  pTransThreshold = 1.; // MeV/c
644  } else {
645  // More stringent checks before accommodating recoil
646  EThreshold = 0.1; // MeV
647  pLongThreshold = 0.1; // MeV/c
648  pTransThreshold = 0.1; // MeV/c
649  }
650  if(std::abs(theBalance.energy)>EThreshold) {
651  INCL_WARN("Violation of energy conservation > " << EThreshold << " MeV. EBalance = " << theBalance.energy << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
652  }
653  if(std::abs(pLongBalance)>pLongThreshold) {
654  INCL_WARN("Violation of longitudinal momentum conservation > " << pLongThreshold << " MeV/c. pLongBalance = " << pLongBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
655  }
656  if(std::abs(pTransBalance)>pTransThreshold) {
657  INCL_WARN("Violation of transverse momentum conservation > " << pTransThreshold << " MeV/c. pTransBalance = " << pTransBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
658  }
659 
660  // Feed the EventInfo variables
661  theEventInfo.EBalance = theBalance.energy;
662  theEventInfo.pLongBalance = pLongBalance;
663  theEventInfo.pTransBalance = pTransBalance;
664  }
665 #endif
666 
667  G4bool INCL::continueCascade() {
668  // Stop if we have passed the stopping time
669  if(propagationModel->getCurrentTime() > propagationModel->getStoppingTime()) {
670  INCL_DEBUG("Cascade time (" << propagationModel->getCurrentTime()
671  << ") exceeded stopping time (" << propagationModel->getStoppingTime()
672  << "), stopping cascade" << '\n');
673  return false;
674  }
675  // Stop if there are no participants and no pions inside the nucleus
676  if(nucleus->getStore()->getBook().getCascading()==0 &&
677  nucleus->getStore()->getIncomingParticles().empty()) {
678  INCL_DEBUG("No participants in the nucleus and no incoming particles left, stopping cascade" << '\n');
679  return false;
680  }
681  // Stop if the remnant is smaller than minRemnantSize
682  if(nucleus->getA() <= minRemnantSize) {
683  INCL_DEBUG("Remnant size (" << nucleus->getA()
684  << ") smaller than or equal to minimum (" << minRemnantSize
685  << "), stopping cascade" << '\n');
686  return false;
687  }
688  // Stop if we have to try and make a compound nucleus or if we have to
689  // force a transparent
690  if(nucleus->getTryCompoundNucleus()) {
691  INCL_DEBUG("Trying to make a compound nucleus, stopping cascade" << '\n');
692  return false;
693  }
694 
695  return true;
696  }
697 
698  void INCL::finalizeGlobalInfo(Random::SeedVector const &initialSeeds) {
699  const G4double normalisationFactor = theGlobalInfo.geometricCrossSection /
700  ((G4double) theGlobalInfo.nShots);
701  theGlobalInfo.nucleonAbsorptionCrossSection = normalisationFactor *
702  ((G4double) theGlobalInfo.nNucleonAbsorptions);
703  theGlobalInfo.pionAbsorptionCrossSection = normalisationFactor *
704  ((G4double) theGlobalInfo.nPionAbsorptions);
705  theGlobalInfo.reactionCrossSection = normalisationFactor *
706  ((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents));
707  theGlobalInfo.errorReactionCrossSection = normalisationFactor *
708  std::sqrt((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents));
709  theGlobalInfo.forcedCNCrossSection = normalisationFactor *
710  ((G4double) theGlobalInfo.nForcedCompoundNucleus);
711  theGlobalInfo.errorForcedCNCrossSection = normalisationFactor *
712  std::sqrt((G4double) (theGlobalInfo.nForcedCompoundNucleus));
713  theGlobalInfo.completeFusionCrossSection = normalisationFactor *
714  ((G4double) theGlobalInfo.nCompleteFusion);
715  theGlobalInfo.errorCompleteFusionCrossSection = normalisationFactor *
716  std::sqrt((G4double) (theGlobalInfo.nCompleteFusion));
717  theGlobalInfo.energyViolationInteractionCrossSection = normalisationFactor *
718  ((G4double) theGlobalInfo.nEnergyViolationInteraction);
719 
720  theGlobalInfo.initialRandomSeeds.assign(initialSeeds.begin(), initialSeeds.end());
721 
723  theGlobalInfo.finalRandomSeeds.assign(theSeeds.begin(), theSeeds.end());
724  }
725 
726  G4int INCL::makeProjectileRemnant() {
727  // Do nothing if this is not a nucleus-nucleus reaction
728  if(!nucleus->getProjectileRemnant())
729  return 0;
730 
731  // Get the spectators (geometrical+dynamical) from the Store
732  ParticleList geomSpectators(nucleus->getProjectileRemnant()->getParticles());
733  ParticleList dynSpectators(nucleus->getStore()->extractDynamicalSpectators());
734 
735  G4int nUnmergedSpectators = 0;
736 
737  // If there are no spectators, do nothing
738  if(dynSpectators.empty() && geomSpectators.empty()) {
739  return 0;
740  } else if(dynSpectators.size()==1 && geomSpectators.empty()) {
741  // No geometrical spectators, one dynamical spectator
742  // Just put it back in the outgoing list
743  nucleus->getStore()->addToOutgoing(dynSpectators.front());
744  } else {
745  // Make a cluster out of the geometrical spectators
746  ProjectileRemnant *theProjectileRemnant = nucleus->getProjectileRemnant();
747 
748  // Add the dynamical spectators to the bunch
749  ParticleList rejected = theProjectileRemnant->addAllDynamicalSpectators(dynSpectators);
750  // Put back the rejected spectators into the outgoing list
751  nUnmergedSpectators = rejected.size();
752  nucleus->getStore()->addToOutgoing(rejected);
753 
754  // Deal with the projectile remnant
755  nucleus->finalizeProjectileRemnant(propagationModel->getCurrentTime());
756 
757  }
758 
759  return nUnmergedSpectators;
760  }
761 
762  void INCL::initMaxInteractionDistance(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
763  if(projectileSpecies.theType != Composite) {
764  maxInteractionDistance = 0.;
765  return;
766  }
767 
770 
771  const G4double theNNDistance = CrossSections::interactionDistanceNN(projectileSpecies, kineticEnergy);
772  maxInteractionDistance = r0 + theNNDistance;
773  INCL_DEBUG("Initialised interaction distance: r0 = " << r0 << '\n'
774  << " theNNDistance = " << theNNDistance << '\n'
775  << " maxInteractionDistance = " << maxInteractionDistance << '\n');
776  }
777 
778  void INCL::initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z) {
779  G4double rMax = 0.0;
780  if(A==0) {
781  IsotopicDistribution const &anIsotopicDistribution =
783  IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
784  for(IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
785  const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, i->theA, Z);
786  const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, i->theA, Z);
787  const G4double maximumRadius = std::min(pMaximumRadius, nMaximumRadius);
788  rMax = std::max(maximumRadius, rMax);
789  }
790  } else {
791  const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
792  const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
793  rMax = std::min(pMaximumRadius, nMaximumRadius);
794  }
795  if(p.theType==Composite || p.theType==Proton || p.theType==Neutron) {
797  maxUniverseRadius = rMax + interactionDistanceNN;
798  } else if(p.theType==PiPlus
799  || p.theType==PiZero
800  || p.theType==PiMinus) {
802  maxUniverseRadius = rMax + interactionDistancePiN;
803  }
804  INCL_DEBUG("Initialised universe radius: " << maxUniverseRadius << '\n');
805  }
806 
807  void INCL::updateGlobalInfo() {
808  // Increment the global counter for the number of shots
809  theGlobalInfo.nShots++;
810 
811  if(theEventInfo.transparent) {
812  // Increment the global counter for the number of transparents
813  theGlobalInfo.nTransparents++;
814  // Increment the global counter for the number of forced transparents
815  if(forceTransparent)
816  theGlobalInfo.nForcedTransparents++;
817  return;
818  }
819 
820  // Check if we have an absorption:
821  if(theEventInfo.nucleonAbsorption) theGlobalInfo.nNucleonAbsorptions++;
822  if(theEventInfo.pionAbsorption) theGlobalInfo.nPionAbsorptions++;
823 
824  // Count complete-fusion events
825  if(theEventInfo.nCascadeParticles==0) theGlobalInfo.nCompleteFusion++;
826 
827  if(nucleus->getTryCompoundNucleus())
828  theGlobalInfo.nForcedCompoundNucleus++;
829 
830  // Counters for the number of violations of energy conservation in
831  // collisions
832  theGlobalInfo.nEnergyViolationInteraction += theEventInfo.nEnergyViolationInteraction;
833  }
834 
835 }
Short_t Ap
Projectile mass number given as input.
Short_t Zp
Charge number of the projectile nucleus.
G4int getA() const
Returns the baryon number.
Float_t effectiveImpactParameter
Effective (Coulomb-distorted) impact parameter [fm].
CascadeActionType getCascadeActionType() const
Get the cascade-action type.
ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const
Compute charge, mass, energy and momentum balance.
void initializeParticles()
Float_t nucleonAbsorptionCrossSection
Nucleon absorption cross section.
Bool_t nucleonAbsorption
True if the event is a nucleon absorption.
G4bool containsDeltas()
Returns true if the nucleus contains any deltas.
IsotopeVector::iterator IsotopeIter
void setMass(G4double mass)
static void setCutNN(const G4double c)
void deleteGenerator()
G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z)
Float_t pLongBalance
Longitudinal momentum-conservation balance [MeV/c].
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
G4int getTargetZ() const
Get the target charge number.
Definition: G4INCLConfig.hh:97
Int_t nForcedCompoundNucleus
Number of forced compound-nucleus events.
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
std::string cascadeModel
Name of the cascade model.
ParticleList const & getParticles() const
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
Abstract interface to the nuclear potential.
virtual G4INCL::IAvatar * propagate(FinalState const *const fs)=0
Float_t errorReactionCrossSection
Error on the calculated reaction cross section.
const char * p
Definition: xmltok.h:285
const G4double tenPi
void applyFinalState(FinalState *)
#define INCL_ERROR(x)
void setZ(const G4int Z)
Set the charge number of the cluster.
void initialize(Config const *const)
Initialize generator according to a Config object.
const ThreeVector & getIncomingAngularMomentum() const
Get the incoming angular-momentum vector.
Store * getStore() const
G4double interactionDistanceNN(const ParticleSpecies &aSpecies, const G4double kineticEnergy)
Compute the &quot;interaction distance&quot;.
Float_t forcedCNCrossSection
Calculated forced-compound-nucleus cross section.
Float_t pionAbsorptionCrossSection
Pion absorption cross section.
std::vector< Int_t > finalRandomSeeds
Final seeds for the pseudo-random-number generator.
Bool_t forcedDeltasOutside
Event involved forced delta decays outside the nucleus.
G4double getEnergy() const
#define INCL_WARN(x)
Bool_t pionAbsorption
True if the event is a pion absorption.
G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n)
Return the maximum impact parameter for Coulomb-distorted trajectories.
G4bool isNaturalTarget() const
Natural targets.
Definition: G4INCLConfig.hh:87
Bool_t clusterDecay
Event involved cluster decay.
G4bool hasRemnant() const
Does the nucleus give a cascade remnant?
void deleteIncoming()
Clear the incoming list and delete the particles.
Definition: G4INCLStore.hh:150
ParticleList const & getIncomingParticles() const
Definition: G4INCLStore.hh:217
Float_t Ep
Projectile kinetic energy given as input.
Short_t At
Target mass number given as input.
Int_t nTransparents
Number of transparent shots.
G4double shoot0()
void reset()
Reset the projectile remnant to the state at the beginning of the cascade.
std::vector< Int_t > initialRandomSeeds
Initial seeds for the pseudo-random-number generator.
int G4int
Definition: G4Types.hh:78
Float_t stoppingTime
Cascade stopping time [fm/c].
Class containing default actions to be performed at intermediate cascade steps.
Float_t Ep
Projectile kinetic energy given as input.
Short_t At
Mass number of the target nucleus.
G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
void setEnergy(G4double energy)
static G4ThreadLocal Int_t eventNumber
Number of the event.
SeedVector getSeeds()
Definition: G4INCLRandom.cc:89
void clearCache()
Clear the INuclearPotential cache.
G4double getImpactParameter() const
void finalizeGlobalInfo(Random::SeedVector const &initialSeeds)
Short_t Zt
Charge number of the target nucleus.
Float_t impactParameter
Impact parameter [fm].
void useFusionKinematics()
Adjust the kinematics for complete-fusion events.
void deleteCoulomb()
Delete the Coulomb-distortion object.
Static class for selecting Coulomb distortion.
INCL(Config const *const config)
static std::string const getVersionString()
Get the INCL version string.
Simple container for output of calculation-wide results.
double A(double temperature)
Float_t pTransBalance
Transverse momentum-conservation balance [MeV/c].
void deleteBlockers()
Delete blockers.
Float_t EBalance
Energy-conservation balance [MeV].
G4bool decayOutgoingDeltas()
Force the decay of outgoing deltas.
ParticleList const & getOutgoingParticles() const
Definition: G4INCLStore.hh:223
static void deleteBackupParticles()
Release the memory allocated for the backup particles.
G4int getCascading() const
Definition: G4INCLBook.hh:105
virtual G4double shoot(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)=0
Book & getBook()
Definition: G4INCLStore.hh:259
Bool_t deltasInside
Event involved deltas in the nucleus at the end of the cascade.
ParticleList extractDynamicalSpectators()
Returns a list of dynamical spectators.
Definition: G4INCLStore.hh:232
Int_t nForcedTransparents
Number of forced transparents.
std::string deexcitationModel
Name of the de-excitation model.
bool G4bool
Definition: G4Types.hh:79
void initialize(Config const *const theConfig)
G4double getDecayTimeThreshold() const
Get the hadronization time.
G4int getTargetA() const
Get the target mass number.
Definition: G4INCLConfig.hh:94
Intersection getEarlierTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r)
Compute the first intersection of a straight particle trajectory with a sphere.
const EventInfo & processEvent()
G4int getZ() const
Returns the charge number.
Alternative CascadeAction for dumping avatars.
G4double getExcitationEnergy() const
Get the excitation energy of the nucleus.
G4double computeExcitationEnergy() const
Compute the current excitation energy.
Float_t completeFusionCrossSection
Calculated complete-fusion cross section (nParticles==0)
Int_t nShots
Number of shots.
Int_t nUnmergedSpectators
Number of dynamical spectators that were merged back into the projectile remnant. ...
G4bool decayMe()
Force the phase-space decay of the Nucleus.
void initialize(Config const *const theConfig)
Initialize the clustering model based on the Config object.
void clearOutgoing()
Definition: G4INCLStore.cc:223
std::vector< Isotope > IsotopeVector
G4double getHadronizationTime() const
Get the hadronization time.
Functions that encapsulate a mass table.
Float_t reactionCrossSection
Calculated reaction cross section.
G4bool getTryCompoundNucleus()
Float_t geometricCrossSection
Geometric cross section.
G4bool initializeTarget(const G4int A, const G4int Z)
Short_t Zp
Projectile charge number given as input.
Adapter const & getAdapter()
void initVerbosityLevelFromEnvvar()
Bool_t forcedPionResonancesOutside
Event involved forced eta/omega decays outside the nucleus.
void addParticleEntryAvatar(IAvatar *a)
Add one ParticleEntry avatar.
Definition: G4INCLStore.cc:66
Float_t errorCompleteFusionCrossSection
Error on the calculated complete-fusion cross section (nParticles==0)
void reset()
Definition: G4INCLBook.hh:52
void setA(const G4int A)
Set the mass number of the cluster.
ParticleList addAllDynamicalSpectators(ParticleList const &pL)
Add back all dynamical spectators to the projectile remnant.
G4bool decayInsideDeltas()
Force the decay of deltas inside the nucleus.
void finalizeProjectileRemnant(const G4double emissionTime)
Finalise the projectile remnant.
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
Static class for cluster formation.
void deleteClusteringModel()
Delete the clustering model.
std::string getDeExcitationString() const
Get the de-excitation string.
ParticleSpecies getProjectileSpecies() const
Get the projectile species.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual void setNucleus(G4INCL::Nucleus *nucleus)=0
Int_t nPionAbsorptions
Number of nucleon absorptions (no outcoming pions)
G4double interactionDistancePiN(const G4double projectileKineticEnergy)
Compute the &quot;interaction distance&quot;.
LocalEnergyType getLocalEnergyPiType() const
Get the type of local energy for pi-N and decay avatars.
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
LocalEnergyType getLocalEnergyBBType() const
Get the type of local energy for N-N avatars.
void initialize(Config const *const theConfig)
void reset()
Reset the EventInfo members.
virtual G4double getStoppingTime()=0
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4bool isEventTransparent() const
Is the event transparent?
const G4double twoPi
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
Solution solve(RootFunctor const *const f, const G4double x0)
Numerically solve a one-dimensional equation.
void fillEventInfo(EventInfo *eventInfo)
G4int drawRandomNaturalIsotope(const G4int Z)
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
G4double shoot()
Definition: G4INCLRandom.cc:93
void initialize(Config const *const aConfig)
Initialise blockers according to a Config object.
Short_t Zt
Target charge number given as input.
Bool_t transparent
True if the event is transparent.
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
double G4double
Definition: G4Types.hh:76
#define INCL_DEBUG(x)
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
virtual G4double getCurrentTime()=0
G4double getProjectileKineticEnergy() const
Get the projectile kinetic energy.
Int_t nNucleonAbsorptions
Number of nucleon absorptions (no outcoming particles)
Float_t energyViolationInteractionCrossSection
Cross section for attempted collisions/decays for which the energy-conservation algorithm failed to f...
void addToOutgoing(Particle *p)
add the particle to the outgoing particle list.
Definition: G4INCLStore.hh:190
Bool_t forcedDeltasInside
Event involved forced delta decays inside the nucleus.
Short_t Ap
Mass number of the projectile nucleus.
Short_t nCascadeParticles
Number of cascade particles.
G4bool decayOutgoingPionResonances(G4double timeThreshold)
Force the decay of outgoing PionResonances (eta/omega).
void setExcitationEnergy(const G4double e)
Set the excitation energy of the cluster.
void initialize(Config const *const theConfig=0)
Initialize the particle table.
void setSpin(const ThreeVector &j)
Set the spin of the nucleus.
Float_t errorForcedCNCrossSection
Error on the calculated forced-compound-nucleus cross section.
G4double getCutNN() const
G4double getZ() const
void distortOut(ParticleList const &pL, Nucleus const *const n)
Modify the momentum of an outgoing particle.
Simple class for computing intersections between a straight line and a sphere.
void initialize(Config const *const theConfig)
Initialize the Coulomb-distortion algorithm.
void clearIncoming()
Clear the incoming list.
Definition: G4INCLStore.hh:145
Int_t nCompleteFusion
Number of complete-fusion events (nParticles==0)
Int_t projectileType
Projectile particle type.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)