Geant4  10.02
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
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.
101 
102  // Select the Pauli and CDPP blocking algorithms
104 
105  // Set the cross-section set
107 
108  // Set the phase-space generator
110 
111  // Select the Coulomb-distortion algorithm:
113 
114  // Select the clustering algorithm:
116 
117  // Initialize the INCL particle table:
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.
133  else
136 
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
147  const ParticleSpecies theSpecies = theConfig->getProjectileSpecies();
148  theGlobalInfo.Ap = theSpecies.theA;
149  theGlobalInfo.Zp = theSpecies.theZ;
151 #endif
152 
154  }
155 
158 #ifndef INCLXX_IN_GEANT4_MODE
159  NuclearMassTable::deleteTable();
160 #endif
167 #ifndef INCLXX_IN_GEANT4_MODE
168  Logger::deleteLoggerSlave();
169 #endif
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;
201  else
202  theA = A;
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
214  Math::tenPi*std::pow(maxImpactParameter,2);
215 
216  // Set the minimum remnant size
217  if(projectileSpecies.theA > 0)
219  else
220  minRemnantSize = std::min(theA-1, 4);
221 
222  return true;
223  }
224 
226  delete nucleus;
227 
229  nucleus->getStore()->getBook().reset();
231 
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');
248  return theEventInfo;
249  }
250 
252 
253  const G4bool canRunCascade = preCascade(projectileSpecies, kineticEnergy);
254  if(canRunCascade) {
255  cascade();
256  postCascade();
258  }
260  return theEventInfo;
261  }
262 
263  G4bool INCL::preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
264  // Reset theEventInfo
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;
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:
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:
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:
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 
363  // Fill in the event information
365 
366  // Forced CN?
368  INCL_DEBUG("Trying compound nucleus" << '\n');
371  // Global checks of conservation laws
372 #ifndef INCLXX_IN_GEANT4_MODE
373  if(!theEventInfo.transparent) globalConservationChecks(true);
374 #endif
375  return;
376  }
377 
379 
381  ProjectileRemnant * const projectileRemnant = nucleus->getProjectileRemnant();
382  if(projectileRemnant) {
383  // Clear the incoming list (particles will be deleted by the ProjectileRemnant)
385  } else {
386  // Delete particles in the incoming list
388  }
389  } else {
390  // Check if the nucleus contains deltas
392 
393  // Take care of any remaining deltas
396 
397  // Apply Coulomb distortion, if appropriate
398  // Note that this will apply Coulomb distortion also on pions emitted by
399  // unphysical remnants (see decayInsideDeltas). This is at variance with
400  // what INCL4.6 does, but these events are (should be!) so rare that
401  // whatever we do doesn't (shouldn't!) make any noticeable difference.
403 
404  // If the normal cascade predicted complete fusion, use the tabulated
405  // masses to compute the excitation energy, the recoil, etc.
406  if(nucleus->getStore()->getOutgoingParticles().size()==0
408  || nucleus->getProjectileRemnant()->getParticles().size()==0)) {
409 
410  INCL_DEBUG("Cascade resulted in complete fusion, using realistic fusion kinematics" << '\n');
411 
413 
414  if(nucleus->getExcitationEnergy()<0.) {
415  // Complete fusion is energetically impossible, return a transparent
416  INCL_WARN("Complete-fusion kinematics yields negative excitation energy, returning a transparent!" << '\n');
417  theEventInfo.transparent = true;
418  return;
419  }
420 
421  } else { // Normal cascade here
422 
423  // Set the excitation energy
425 
426  // Make a projectile pre-fragment out of the geometrical and dynamical
427  // spectators
429 
430  // Compute recoil momentum, energy and spin of the nucleus
431  if(nucleus->getA()==1 && minRemnantSize>1) {
432  INCL_ERROR("Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << '\n');
433  }
435 
436 #ifndef INCLXX_IN_GEANT4_MODE
437  // Global checks of conservation laws
438  globalConservationChecks(false);
439 #endif
440 
441  // Make room for the remnant recoil by rescaling the energies of the
442  // outgoing particles.
444 
445  }
446 
447  // Cluster decay
449 
450 #ifndef INCLXX_IN_GEANT4_MODE
451  // Global checks of conservation laws
452  globalConservationChecks(true);
453 #endif
454 
455  // Fill the EventInfo structure
457 
458  }
459  }
460 
462  // If this is not a nucleus-nucleus collision, don't attempt to make a
463  // compound nucleus.
464  //
465  // Yes, even nucleon-nucleus collisions can lead to particles entering
466  // below the Fermi level. Take e.g. 1-MeV p + He4.
468  forceTransparent = true;
469  return;
470  }
471 
472  // Reset the internal Nucleus variables
478 
479  // CN kinematical variables
480  // Note: the CN orbital angular momentum is neglected in what follows. We
481  // should actually take it into account!
482  ThreeVector theCNMomentum = nucleus->getIncomingMomentum();
485  G4int theCNA=theEventInfo.At, theCNZ=theEventInfo.Zt;
486  Cluster * const theProjectileRemnant = nucleus->getProjectileRemnant();
487  G4double theCNEnergy = theTargetMass + theProjectileRemnant->getEnergy();
488 
489  // Loop over the potential participants
490  ParticleList const &initialProjectileComponents = theProjectileRemnant->getParticles();
491  std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
492  // Shuffle the list of potential participants
493  std::random_shuffle(shuffledComponents.begin(), shuffledComponents.end(), Random::getAdapter());
494 
495  G4bool success = true;
496  G4bool atLeastOneNucleonEntering = false;
497  for(std::vector<Particle*>::const_iterator p=shuffledComponents.begin(), e=shuffledComponents.end(); p!=e; ++p) {
498  // Skip particles that miss the interaction distance
500  (*p)->getPosition(),
501  (*p)->getPropagationVelocity(),
503  if(!intersectionInteractionDistance.exists)
504  continue;
505 
506  // Build an entry avatar for this nucleon
507  atLeastOneNucleonEntering = true;
508  ParticleEntryAvatar *theAvatar = new ParticleEntryAvatar(0.0, nucleus, *p);
509  nucleus->getStore()->addParticleEntryAvatar(theAvatar);
510  FinalState *fs = theAvatar->getFinalState();
512  FinalStateValidity validity = fs->getValidity();
513  delete fs;
514  switch(validity) {
515  case ValidFS:
517  case ParticleBelowZeroFS:
518  // Add the particle to the CN
519  theCNA++;
520  theCNZ += (*p)->getZ();
521  break;
522  case PauliBlockedFS:
524  default:
525  success = false;
526  break;
527  }
528  }
529 
530  if(!success || !atLeastOneNucleonEntering) {
531  INCL_DEBUG("No nucleon entering in forced CN, forcing a transparent" << '\n');
532  forceTransparent = true;
533  return;
534  }
535 
536 // assert(theCNA==nucleus->getA());
537 // assert(theCNA<=theEventInfo.At+theEventInfo.Ap);
538 // assert(theCNZ<=theEventInfo.Zt+theEventInfo.Zp);
539 
540  // Update the kinematics of the CN
541  theCNEnergy -= theProjectileRemnant->getEnergy();
542  theCNMomentum -= theProjectileRemnant->getMomentum();
543 
544  // Deal with the projectile remnant
546 
547  // Subtract the angular momentum of the projectile remnant
548 // assert(nucleus->getStore()->getOutgoingParticles().empty());
549  theCNSpin -= theProjectileRemnant->getAngularMomentum();
550 
551  // Compute the excitation energy of the CN
552  const G4double theCNMass = ParticleTable::getTableMass(theCNA,theCNZ);
553  const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.mag2();
554  if(theCNInvariantMassSquared<0.) {
555  // Negative invariant mass squared, return a transparent
556  forceTransparent = true;
557  return;
558  }
559  const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
560  if(theCNExcitationEnergy<0.) {
561  // Negative excitation energy, return a transparent
562  INCL_DEBUG("CN excitation energy is negative, forcing a transparent" << '\n'
563  << " theCNA = " << theCNA << '\n'
564  << " theCNZ = " << theCNZ << '\n'
565  << " theCNEnergy = " << theCNEnergy << '\n'
566  << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
567  << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
568  << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
569  );
570  forceTransparent = true;
571  return;
572  } else {
573  // Positive excitation energy, can make a CN
574  INCL_DEBUG("CN excitation energy is positive, forcing a CN" << '\n'
575  << " theCNA = " << theCNA << '\n'
576  << " theCNZ = " << theCNZ << '\n'
577  << " theCNEnergy = " << theCNEnergy << '\n'
578  << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
579  << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
580  << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
581  );
582  nucleus->setA(theCNA);
583  nucleus->setZ(theCNZ);
584  nucleus->setMomentum(theCNMomentum);
585  nucleus->setEnergy(theCNEnergy);
586  nucleus->setExcitationEnergy(theCNExcitationEnergy);
587  nucleus->setMass(theCNMass+theCNExcitationEnergy);
588  nucleus->setSpin(theCNSpin); // neglects any orbital angular momentum of the CN
589 
590  // Take care of any remaining deltas
592 
593  // Cluster decay
595 
596  // Fill the EventInfo structure
598  }
599  }
600 
602  RecoilCMFunctor theRecoilFunctor(nucleus, theEventInfo);
603 
604  // Apply the root-finding algorithm
605  const RootFinder::Solution theSolution = RootFinder::solve(&theRecoilFunctor, 1.0);
606  if(theSolution.success) {
607  theRecoilFunctor(theSolution.x); // Apply the solution
608  } else {
609  INCL_WARN("Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." << '\n');
610  }
611 
612  }
613 
614 #ifndef INCLXX_IN_GEANT4_MODE
615  void INCL::globalConservationChecks(G4bool afterRecoil) {
617 
618  // Global conservation checks
619  const G4double pLongBalance = theBalance.momentum.getZ();
620  const G4double pTransBalance = theBalance.momentum.perp();
621  if(theBalance.Z != 0) {
622  INCL_ERROR("Violation of charge conservation! ZBalance = " << theBalance.Z << '\n');
623  }
624  if(theBalance.A != 0) {
625  INCL_ERROR("Violation of baryon-number conservation! ABalance = " << theBalance.A << '\n');
626  }
627  G4double EThreshold, pLongThreshold, pTransThreshold;
628  if(afterRecoil) {
629  // Less stringent checks after accommodating recoil
630  EThreshold = 10.; // MeV
631  pLongThreshold = 1.; // MeV/c
632  pTransThreshold = 1.; // MeV/c
633  } else {
634  // More stringent checks before accommodating recoil
635  EThreshold = 0.1; // MeV
636  pLongThreshold = 0.1; // MeV/c
637  pTransThreshold = 0.1; // MeV/c
638  }
639  if(std::abs(theBalance.energy)>EThreshold) {
640  INCL_WARN("Violation of energy conservation > " << EThreshold << " MeV. EBalance = " << theBalance.energy << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
641  }
642  if(std::abs(pLongBalance)>pLongThreshold) {
643  INCL_WARN("Violation of longitudinal momentum conservation > " << pLongThreshold << " MeV/c. pLongBalance = " << pLongBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
644  }
645  if(std::abs(pTransBalance)>pTransThreshold) {
646  INCL_WARN("Violation of transverse momentum conservation > " << pTransThreshold << " MeV/c. pTransBalance = " << pTransBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
647  }
648 
649  // Feed the EventInfo variables
650  theEventInfo.EBalance = theBalance.energy;
651  theEventInfo.pLongBalance = pLongBalance;
652  theEventInfo.pTransBalance = pTransBalance;
653  }
654 #endif
655 
657  // Stop if we have passed the stopping time
659  INCL_DEBUG("Cascade time (" << propagationModel->getCurrentTime()
660  << ") exceeded stopping time (" << propagationModel->getStoppingTime()
661  << "), stopping cascade" << '\n');
662  return false;
663  }
664  // Stop if there are no participants and no pions inside the nucleus
665  if(nucleus->getStore()->getBook().getCascading()==0 &&
666  nucleus->getStore()->getIncomingParticles().empty()) {
667  INCL_DEBUG("No participants in the nucleus and no incoming particles left, stopping cascade" << '\n');
668  return false;
669  }
670  // Stop if the remnant is smaller than minRemnantSize
671  if(nucleus->getA() <= minRemnantSize) {
672  INCL_DEBUG("Remnant size (" << nucleus->getA()
673  << ") smaller than or equal to minimum (" << minRemnantSize
674  << "), stopping cascade" << '\n');
675  return false;
676  }
677  // Stop if we have to try and make a compound nucleus or if we have to
678  // force a transparent
680  INCL_DEBUG("Trying to make a compound nucleus, stopping cascade" << '\n');
681  return false;
682  }
683 
684  return true;
685  }
686 
687  void INCL::finalizeGlobalInfo(Random::SeedVector const &initialSeeds) {
688  const G4double normalisationFactor = theGlobalInfo.geometricCrossSection /
690  theGlobalInfo.nucleonAbsorptionCrossSection = normalisationFactor *
692  theGlobalInfo.pionAbsorptionCrossSection = normalisationFactor *
694  theGlobalInfo.reactionCrossSection = normalisationFactor *
696  theGlobalInfo.errorReactionCrossSection = normalisationFactor *
698  theGlobalInfo.forcedCNCrossSection = normalisationFactor *
700  theGlobalInfo.errorForcedCNCrossSection = normalisationFactor *
702  theGlobalInfo.completeFusionCrossSection = normalisationFactor *
704  theGlobalInfo.errorCompleteFusionCrossSection = normalisationFactor *
705  std::sqrt((G4double) (theGlobalInfo.nCompleteFusion));
708 
709  theGlobalInfo.initialRandomSeeds.assign(initialSeeds.begin(), initialSeeds.end());
710 
712  theGlobalInfo.finalRandomSeeds.assign(theSeeds.begin(), theSeeds.end());
713  }
714 
716  // Do nothing if this is not a nucleus-nucleus reaction
718  return 0;
719 
720  // Get the spectators (geometrical+dynamical) from the Store
723 
724  G4int nUnmergedSpectators = 0;
725 
726  // If there are no spectators, do nothing
727  if(dynSpectators.empty() && geomSpectators.empty()) {
728  return 0;
729  } else if(dynSpectators.size()==1 && geomSpectators.empty()) {
730  // No geometrical spectators, one dynamical spectator
731  // Just put it back in the outgoing list
732  nucleus->getStore()->addToOutgoing(dynSpectators.front());
733  } else {
734  // Make a cluster out of the geometrical spectators
735  ProjectileRemnant *theProjectileRemnant = nucleus->getProjectileRemnant();
736 
737  // Add the dynamical spectators to the bunch
738  ParticleList rejected = theProjectileRemnant->addAllDynamicalSpectators(dynSpectators);
739  // Put back the rejected spectators into the outgoing list
740  nUnmergedSpectators = rejected.size();
741  nucleus->getStore()->addToOutgoing(rejected);
742 
743  // Deal with the projectile remnant
745 
746  }
747 
748  return nUnmergedSpectators;
749  }
750 
751  void INCL::initMaxInteractionDistance(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
752  if(projectileSpecies.theType != Composite) {
754  return;
755  }
756 
759 
760  const G4double theNNDistance = CrossSections::interactionDistanceNN(projectileSpecies, kineticEnergy);
761  maxInteractionDistance = r0 + theNNDistance;
762  INCL_DEBUG("Initialised interaction distance: r0 = " << r0 << '\n'
763  << " theNNDistance = " << theNNDistance << '\n'
764  << " maxInteractionDistance = " << maxInteractionDistance << '\n');
765  }
766 
767  void INCL::initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z) {
768  G4double rMax = 0.0;
769  if(A==0) {
770  IsotopicDistribution const &anIsotopicDistribution =
772  IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
773  for(IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
774  const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, i->theA, Z);
775  const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, i->theA, Z);
776  const G4double maximumRadius = std::min(pMaximumRadius, nMaximumRadius);
777  rMax = std::max(maximumRadius, rMax);
778  }
779  } else {
780  const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
781  const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
782  rMax = std::min(pMaximumRadius, nMaximumRadius);
783  }
784  if(p.theType==Composite || p.theType==Proton || p.theType==Neutron) {
787  } else if(p.theType==PiPlus
788  || p.theType==PiZero
789  || p.theType==PiMinus) {
792  }
793  INCL_DEBUG("Initialised universe radius: " << maxUniverseRadius << '\n');
794  }
795 
797  // Increment the global counter for the number of shots
799 
801  // Increment the global counter for the number of transparents
803  // Increment the global counter for the number of forced transparents
804  if(forceTransparent)
806  return;
807  }
808 
809  // Check if we have an absorption:
812 
813  // Count complete-fusion events
815 
818 
819  // Counters for the number of violations of energy conservation in
820  // collisions
822  }
823 
824 }
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].
G4double maxUniverseRadius
void fillFinalState(FinalState *fs)
G4bool targetInitSuccess
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()
Call the Cluster method to generate the initial distribution of particles.
Float_t nucleonAbsorptionCrossSection
Nucleon absorption cross section.
G4int minRemnantSize
Remnant size below which cascade stops.
Bool_t nucleonAbsorption
True if the event is a nucleon absorption.
G4bool containsDeltas()
Returns true if the nucleus contains any deltas.
IsotopeVector::iterator IsotopeIter
The INCL configuration object.
Definition: G4INCLConfig.hh:60
void setMass(G4double mass)
Set the mass of the particle in MeV/c^2.
static void setCutNN(const G4double c)
FinalStateValidity getValidity() const
void deleteGenerator()
Delete the generator.
G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z)
G4bool continueCascade()
Stopping criterion for the cascade.
void cascade()
The actual cascade loop.
Float_t pLongBalance
Longitudinal momentum-conservation balance [MeV/c].
G4double maxInteractionDistance
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
void postCascade()
Finalise the cascade and clean up.
G4int getTargetZ() const
Get the target charge number.
Definition: G4INCLConfig.hh:97
G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy)
Initialise the cascade.
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
Get the list of particles in the cluster.
void beforeRunAction(Config const *config)
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
Propagate the particles and get the next avatar.
Float_t errorReactionCrossSection
Error on the calculated reaction cross section.
void beforeCascadeAction(IPropagationModel *)
const G4double tenPi
Config const *const theConfig
void applyFinalState(FinalState *)
Apply reaction final state information to the nucleus.
#define INCL_ERROR(x)
void updateGlobalInfo()
Update global counters and other members of theGlobalInfo object.
void setZ(const G4int Z)
Set the charge number of the cluster.
const G4INCL::ThreeVector & getMomentum() const
Get the momentum vector.
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 "interaction distance".
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
Get the energy of the particle in MeV.
Struct for conservation laws.
#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
Return the list of incoming particles (i.e.
Definition: G4INCLStore.hh:217
Float_t Ep
Projectile kinetic energy given as input.
Short_t At
Target mass number given as input.
Cluster is a particle (inherits from the Particle class) that is actually a collection of elementary ...
Int_t nTransparents
Number of transparent shots.
G4double shoot0()
Return a random number in the ]0,1] interval.
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.
CascadeAction * cascadeAction
int G4int
Definition: G4Types.hh:78
Float_t stoppingTime
Cascade stopping time [fm/c].
G4double mag2() const
Get the square of the length.
Class containing default actions to be performed at intermediate cascade steps.
Float_t Ep
Projectile kinetic energy given as input.
IPropagationModel * propagationModel
Short_t At
Mass number of the target nucleus.
G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
void setEnergy(G4double energy)
Set the energy of the particle in MeV.
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
static G4ThreadLocal Int_t eventNumber
Number of the event.
SeedVector getSeeds()
Get the seeds of the current generator.
Definition: G4INCLRandom.cc:89
void clearCache()
Clear the INuclearPotential cache.
void beforeAvatarAction(IAvatar *a, Nucleus *n)
G4double getImpactParameter() const
void finalizeGlobalInfo(Random::SeedVector const &initialSeeds)
Short_t Zt
Charge number of the target nucleus.
Float_t impactParameter
Impact parameter [fm].
G4double fixedImpactParameter
void rescaleOutgoingForRecoil()
Rescale the energies of the outgoing particles.
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)
void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy)
Initialise the maximum interaction distance.
static std::string const getVersionString()
Get the INCL version string.
Final state of an interaction.
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
Return the list of outgoing particles (i.e.
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()
Return the pointer to the Book object which keeps track of various counters.
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
GlobalInfo theGlobalInfo
Int_t nForcedTransparents
Number of forced transparents.
Class that stores isotopic abundances for a given element.
std::string deexcitationModel
Name of the de-excitation model.
void beforePropagationAction(IPropagationModel *pm)
bool G4bool
Definition: G4Types.hh:79
void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z)
Initialize the universe radius.
void initialize(Config const *const theConfig)
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.
Class to adjust remnant recoil in the reaction CM system.
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. ...
void afterCascadeAction(Nucleus *)
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()
Clear all outgoing particles from the store.
Definition: G4INCLStore.cc:223
std::vector< Isotope > IsotopeVector
IsotopeVector const & getIsotopes() const
Get the isotope vector.
G4double getHadronizationTime() const
Get the hadronization time.
Functions that encapsulate a mass table.
Float_t reactionCrossSection
Calculated reaction cross section.
Standard INCL4 particle propagation and avatar prediction.
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()
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.
G4double maxImpactParameter
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
Set the nucleus for the propagation model.
Int_t nPionAbsorptions
Number of nucleon absorptions (no outcoming pions)
G4double interactionDistancePiN(const G4double projectileKineticEnergy)
Compute the "interaction distance".
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.
G4double perp() const
void makeCompoundNucleus()
Make a compound nucleus.
virtual G4double getStoppingTime()=0
Get the current stopping time.
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4bool isEventTransparent() const
Is the event transparent?
FinalState * getFinalState()
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)
Fill the event info which contains INCL output data.
G4int drawRandomNaturalIsotope(const G4int Z)
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
G4double getX() const
G4bool forceTransparent
G4double shoot()
Generate flat distribution of random numbers.
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.
Intersection-point structure.
Bool_t transparent
True if the event is transparent.
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
double G4double
Definition: G4Types.hh:76
void afterPropagationAction(IPropagationModel *pm, IAvatar *avatar)
#define INCL_DEBUG(x)
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
virtual G4double getCurrentTime()=0
Returns the current global time of the system.
Nucleus * nucleus
void afterAvatarAction(IAvatar *a, Nucleus *n, FinalState *fs)
G4double getProjectileKineticEnergy() const
Get the projectile kinetic energy.
G4int makeProjectileRemnant()
Make a projectile pre-fragment out of geometrical spectators.
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.
EventInfo theEventInfo
Short_t Ap
Mass number of the projectile nucleus.
Short_t nCascadeParticles
Number of cascade particles.
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
const G4double r0
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)
Set the momentum vector.
G4double getY() const