ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCascade.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLCascade.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
42 #include "G4INCLCascade.hh"
43 #include "G4INCLRandom.hh"
45 #include "G4INCLParticleTable.hh"
46 #include "G4INCLParticle.hh"
48 #include "G4INCLGlobalInfo.hh"
49 
50 #include "G4INCLPauliBlocking.hh"
51 
52 #include "G4INCLCrossSections.hh"
53 
55 
56 #include "G4INCLLogger.hh"
57 #include "G4INCLGlobals.hh"
59 
61 
63 
64 #include "G4INCLClustering.hh"
65 
66 #include "G4INCLIntersection.hh"
67 
69 
70 #include "G4INCLCascadeAction.hh"
72 
73 #include <cstring>
74 #include <cstdlib>
75 #include <numeric>
76 
77 namespace G4INCL {
78 
79  INCL::INCL(Config const * const config)
80  :propagationModel(0), theA(208), theZ(82), theS(0),
81  targetInitSuccess(false),
83  maxUniverseRadius(0.),
84  maxInteractionDistance(0.),
85  fixedImpactParameter(0.),
86  theConfig(config),
87  nucleus(NULL),
88  forceTransparent(false),
89  minRemnantSize(4)
90  {
91  // Set the logger object.
92 #ifdef INCLXX_IN_GEANT4_MODE
94 #else // INCLXX_IN_GEANT4_MODE
96 #endif // INCLXX_IN_GEANT4_MODE
97 
98  // Set the random number generator algorithm. The system can support
99  // multiple different generator algorithms in a completely
100  // transparent way.
102 
103  // Select the Pauli and CDPP blocking algorithms
105 
106  // Set the cross-section set
108 
109  // Set the phase-space generator
111 
112  // Select the Coulomb-distortion algorithm:
114 
115  // Select the clustering algorithm:
117 
118  // Initialize the INCL particle table:
120 
121  // Initialize the value of cutNN in BinaryCollisionAvatar
123 
124  // Initialize the value of strange cross section bias
126 
127  // Propagation model is responsible for finding avatars and
128  // transporting the particles. In principle this step is "hidden"
129  // behind an abstract interface and the rest of the system does not
130  // care how the transportation and avatar finding is done. This
131  // should allow us to "easily" experiment with different avatar
132  // finding schemes and even to support things like curved
133  // trajectories in the future.
137  else
140 
143 #ifdef INCL_ROOT_USE
144  theGlobalInfo.rootSelection = theConfig->getROOTSelectionString();
145 #endif
146 
147 #ifndef INCLXX_IN_GEANT4_MODE
148  // Fill in the global information
151  const ParticleSpecies theSpecies = theConfig->getProjectileSpecies();
152  theGlobalInfo.Ap = theSpecies.theA;
153  theGlobalInfo.Zp = theSpecies.theZ;
155 #endif
156 
158  }
159 
162 #ifndef INCLXX_IN_GEANT4_MODE
163  NuclearMassTable::deleteTable();
164 #endif
171 #ifndef INCLXX_IN_GEANT4_MODE
172  Logger::deleteLoggerSlave();
173 #endif
177  delete cascadeAction;
178  delete propagationModel;
179  delete theConfig;
180  }
181 
182  G4bool INCL::prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S) {
183  if(A < 0 || A > 300 || Z < 1 || Z > 200) {
184  INCL_ERROR("Unsupported target: A = " << A << " Z = " << Z << " S = " << S << '\n'
185  << "Target configuration rejected." << '\n');
186  return false;
187  }
188  if(projectileSpecies.theType==Composite &&
189  (projectileSpecies.theZ==projectileSpecies.theA || projectileSpecies.theZ==0)) {
190  INCL_ERROR("Unsupported projectile: A = " << projectileSpecies.theA << " Z = " << projectileSpecies.theZ << " S = " << projectileSpecies.theS << '\n'
191  << "Projectile configuration rejected." << '\n');
192  return false;
193  }
194 
195  // Reset the forced-transparent flag
196  forceTransparent = false;
197 
198  // Initialise the maximum universe radius
199  initUniverseRadius(projectileSpecies, kineticEnergy, A, Z);
200 
201  // Initialise the nucleus
202  theZ = Z;
203  theS = S;
206  else
207  theA = A;
209 
210  // Set the maximum impact parameter
211  maxImpactParameter = CoulombDistortion::maxImpactParameter(projectileSpecies, kineticEnergy, nucleus);
212  INCL_DEBUG("Maximum impact parameter initialised: " << maxImpactParameter << '\n');
213 
214  // For forced CN events
215  initMaxInteractionDistance(projectileSpecies, kineticEnergy);
216 
217  // Set the geometric cross section
219  Math::tenPi*std::pow(maxImpactParameter,2);
220 
221  // Set the minimum remnant size
222  if(projectileSpecies.theA > 0)
224  else
225  minRemnantSize = std::min(theA-1, 4);
226 
227  return true;
228  }
229 
231  delete nucleus;
232 
233  nucleus = new Nucleus(A, Z, S, theConfig, maxUniverseRadius);
234  nucleus->getStore()->getBook().reset();
236 
238  return true;
239  }
240 
242  ParticleSpecies const &projectileSpecies,
243  const G4double kineticEnergy,
244  const G4int targetA,
245  const G4int targetZ,
246  const G4int targetS
247  ) {
248  // ReInitialize the bias vector
249  Particle::INCLBiasVector.clear();
250  //Particle::INCLBiasVector.Clear();
252 
253  // Set the target and the projectile
254  targetInitSuccess = prepareReaction(projectileSpecies, kineticEnergy, targetA, targetZ, targetS);
255 
256  if(!targetInitSuccess) {
257  INCL_WARN("Target initialisation failed for A=" << targetA << ", Z=" << targetZ << ", S=" << targetS << '\n');
259  return theEventInfo;
260  }
261 
263 
264  const G4bool canRunCascade = preCascade(projectileSpecies, kineticEnergy);
265  if(canRunCascade) {
266  cascade();
267  postCascade();
269  }
271  return theEventInfo;
272  }
273 
274  G4bool INCL::preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
275  // Reset theEventInfo
277 
279 
280  // Fill in the event information
281  theEventInfo.projectileType = projectileSpecies.theType;
282  theEventInfo.Ap = projectileSpecies.theA;
283  theEventInfo.Zp = projectileSpecies.theZ;
284  theEventInfo.Sp = projectileSpecies.theS;
285  theEventInfo.Ep = kineticEnergy;
289 
290  // Do nothing below the Coulomb barrier
291  if(maxImpactParameter<=0.) {
292  // Fill in the event information
293  theEventInfo.transparent = true;
294 
295  return false;
296  }
297 
298  // Randomly draw an impact parameter or use a fixed value, depending on the
299  // Config option
300  G4double impactParameter, phi;
301  if(fixedImpactParameter<0.) {
302  impactParameter = maxImpactParameter * std::sqrt(Random::shoot0());
303  phi = Random::shoot() * Math::twoPi;
304  } else {
305  impactParameter = fixedImpactParameter;
306  phi = 0.;
307  }
308  INCL_DEBUG("Selected impact parameter: " << impactParameter << '\n');
309 
310  // Fill in the event information
311  theEventInfo.impactParameter = impactParameter;
312 
313  const G4double effectiveImpactParameter = propagationModel->shoot(projectileSpecies, kineticEnergy, impactParameter, phi);
314  if(effectiveImpactParameter < 0.) {
315  // Fill in the event information
316  theEventInfo.transparent = true;
317 
318  return false;
319  }
320 
321  // Fill in the event information
322  theEventInfo.transparent = false;
323  theEventInfo.effectiveImpactParameter = effectiveImpactParameter;
324 
325  return true;
326  }
327 
328  void INCL::cascade() {
329  FinalState *finalState = new FinalState;
330 
331  unsigned long loopCounter = 0;
332  const unsigned long maxLoopCounter = 10000000;
333  do {
334  // Run book keeping actions that should take place before propagation:
336 
337  // Get the avatar with the smallest time and propagate particles
338  // to that point in time.
339  IAvatar *avatar = propagationModel->propagate(finalState);
340 
341  finalState->reset();
342 
343  // Run book keeping actions that should take place after propagation:
345 
346  if(avatar == 0) break; // No more avatars in the avatar list.
347 
348  // Run book keeping actions that should take place before avatar:
350 
351  // Channel is responsible for calculating the outcome of the
352  // selected avatar. There are different kinds of channels. The
353  // class IChannel is, again, an abstract interface that defines
354  // the externally observable behavior of all interaction
355  // channels.
356  // The handling of the channel is transparent to the API.
357  // Final state tells what changed...
358  avatar->fillFinalState(finalState);
359  // Run book keeping actions that should take place after avatar:
360  cascadeAction->afterAvatarAction(avatar, nucleus, finalState);
361 
362  // So now we must give this information to the nucleus
363  nucleus->applyFinalState(finalState);
364  // and now we are ready to process the next avatar!
365 
366  delete avatar;
367 
368  ++loopCounter;
369  } while(continueCascade() && loopCounter<maxLoopCounter); /* Loop checking, 10.07.2015, D.Mancusi */
370 
371  delete finalState;
372  }
373 
375  // Fill in the event information
377 
378  // The event bias
380 
381  // Forced CN?
383  INCL_DEBUG("Trying compound nucleus" << '\n');
386  // Global checks of conservation laws
387 #ifndef INCLXX_IN_GEANT4_MODE
388  if(!theEventInfo.transparent) globalConservationChecks(true);
389 #endif
390  return;
391  }
392 
394 
396  ProjectileRemnant * const projectileRemnant = nucleus->getProjectileRemnant();
397  if(projectileRemnant) {
398  // Clear the incoming list (particles will be deleted by the ProjectileRemnant)
400  } else {
401  // Delete particles in the incoming list
403  }
404  } else {
405 
406  // Check if the nucleus contains strange particles
411 
412  // Capture antiKaons and Sigmas and produce Lambda instead
414 
415  // Emit strange particles still inside the nucleus
418 
419 #ifdef INCLXX_IN_GEANT4_MODE
421 #endif // INCLXX_IN_GEANT4_MODE
422 
423  // Check if the nucleus contains deltas
425 
426  // Take care of any remaining deltas
429 
430  // Take care of any remaining etas, omegas, neutral Sigmas and/or neutral kaons
431  G4double timeThreshold=theConfig->getDecayTimeThreshold();
433  nucleus->decayOutgoingSigmaZero(timeThreshold);
435 
436  // Apply Coulomb distortion, if appropriate
437  // Note that this will apply Coulomb distortion also on pions emitted by
438  // unphysical remnants (see decayInsideDeltas). This is at variance with
439  // what INCL4.6 does, but these events are (should be!) so rare that
440  // whatever we do doesn't (shouldn't!) make any noticeable difference.
442 
443  // If the normal cascade predicted complete fusion, use the tabulated
444  // masses to compute the excitation energy, the recoil, etc.
445  if(nucleus->getStore()->getOutgoingParticles().size()==0
447  || nucleus->getProjectileRemnant()->getParticles().size()==0)) {
448 
449  INCL_DEBUG("Cascade resulted in complete fusion, using realistic fusion kinematics" << '\n');
450 
452 
453  if(nucleus->getExcitationEnergy()<0.) {
454  // Complete fusion is energetically impossible, return a transparent
455  INCL_WARN("Complete-fusion kinematics yields negative excitation energy, returning a transparent!" << '\n');
456  theEventInfo.transparent = true;
457  return;
458  }
459 
460  } else { // Normal cascade here
461 
462  // Set the excitation energy
464 
465  // Make a projectile pre-fragment out of the geometrical and dynamical
466  // spectators
468 
469  // Compute recoil momentum, energy and spin of the nucleus
470  if(nucleus->getA()==1 && minRemnantSize>1) {
471  INCL_ERROR("Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << '\n');
472  }
474 
475 #ifndef INCLXX_IN_GEANT4_MODE
476  // Global checks of conservation laws
477  globalConservationChecks(false);
478 #endif
479 
480  // Make room for the remnant recoil by rescaling the energies of the
481  // outgoing particles.
483 
484  }
485 
486  // Cluster decay
488 
489 #ifndef INCLXX_IN_GEANT4_MODE
490  // Global checks of conservation laws
491  globalConservationChecks(true);
492 #endif
493 
494  // Fill the EventInfo structure
496 
497  }
498  }
499 
501  // If this is not a nucleus-nucleus collision, don't attempt to make a
502  // compound nucleus.
503  //
504  // Yes, even nucleon-nucleus collisions can lead to particles entering
505  // below the Fermi level. Take e.g. 1-MeV p + He4.
507  forceTransparent = true;
508  return;
509  }
510 
511  // Reset the internal Nucleus variables
517 
518  // CN kinematical variables
519  // Note: the CN orbital angular momentum is neglected in what follows. We
520  // should actually take it into account!
521  ThreeVector theCNMomentum = nucleus->getIncomingMomentum();
524  G4int theCNA=theEventInfo.At, theCNZ=theEventInfo.Zt, theCNS=theEventInfo.St;
525  Cluster * const theProjectileRemnant = nucleus->getProjectileRemnant();
526  G4double theCNEnergy = theTargetMass + theProjectileRemnant->getEnergy();
527 
528  // Loop over the potential participants
529  ParticleList const &initialProjectileComponents = theProjectileRemnant->getParticles();
530  std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
531  // Shuffle the list of potential participants
532  std::shuffle(shuffledComponents.begin(), shuffledComponents.end(), Random::getAdapter());
533 
534  G4bool success = true;
535  G4bool atLeastOneNucleonEntering = false;
536  for(std::vector<Particle*>::const_iterator p=shuffledComponents.begin(), e=shuffledComponents.end(); p!=e; ++p) {
537  // Skip particles that miss the interaction distance
539  (*p)->getPosition(),
540  (*p)->getPropagationVelocity(),
542  if(!intersectionInteractionDistance.exists)
543  continue;
544 
545  // Build an entry avatar for this nucleon
546  atLeastOneNucleonEntering = true;
547  ParticleEntryAvatar *theAvatar = new ParticleEntryAvatar(0.0, nucleus, *p);
548  nucleus->getStore()->addParticleEntryAvatar(theAvatar);
549  FinalState *fs = theAvatar->getFinalState();
551  FinalStateValidity validity = fs->getValidity();
552  delete fs;
553  switch(validity) {
554  case ValidFS:
556  case ParticleBelowZeroFS:
557  // Add the particle to the CN
558  theCNA++;
559  theCNZ += (*p)->getZ();
560  theCNS += (*p)->getS();
561  break;
562  case PauliBlockedFS:
564  default:
565  success = false;
566  break;
567  }
568  }
569 
570  if(!success || !atLeastOneNucleonEntering) {
571  INCL_DEBUG("No nucleon entering in forced CN, forcing a transparent" << '\n');
572  forceTransparent = true;
573  return;
574  }
575 
576 // assert(theCNA==nucleus->getA());
577 // assert(theCNA<=theEventInfo.At+theEventInfo.Ap);
578 // assert(theCNZ<=theEventInfo.Zt+theEventInfo.Zp);
579 // assert(theCNS>=theEventInfo.St+theEventInfo.Sp);
580 
581  // Update the kinematics of the CN
582  theCNEnergy -= theProjectileRemnant->getEnergy();
583  theCNMomentum -= theProjectileRemnant->getMomentum();
584 
585  // Deal with the projectile remnant
587 
588  // Subtract the angular momentum of the projectile remnant
589 // assert(nucleus->getStore()->getOutgoingParticles().empty());
590  theCNSpin -= theProjectileRemnant->getAngularMomentum();
591 
592  // Compute the excitation energy of the CN
593  const G4double theCNMass = ParticleTable::getTableMass(theCNA,theCNZ,theCNS);
594  const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.mag2();
595  if(theCNInvariantMassSquared<0.) {
596  // Negative invariant mass squared, return a transparent
597  forceTransparent = true;
598  return;
599  }
600  const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
601  if(theCNExcitationEnergy<0.) {
602  // Negative excitation energy, return a transparent
603  INCL_DEBUG("CN excitation energy is negative, forcing a transparent" << '\n'
604  << " theCNA = " << theCNA << '\n'
605  << " theCNZ = " << theCNZ << '\n'
606  << " theCNS = " << theCNS << '\n'
607  << " theCNEnergy = " << theCNEnergy << '\n'
608  << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
609  << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
610  << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
611  );
612  forceTransparent = true;
613  return;
614  } else {
615  // Positive excitation energy, can make a CN
616  INCL_DEBUG("CN excitation energy is positive, forcing a CN" << '\n'
617  << " theCNA = " << theCNA << '\n'
618  << " theCNZ = " << theCNZ << '\n'
619  << " theCNS = " << theCNS << '\n'
620  << " theCNEnergy = " << theCNEnergy << '\n'
621  << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
622  << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
623  << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
624  );
625  nucleus->setA(theCNA);
626  nucleus->setZ(theCNZ);
627  nucleus->setS(theCNS);
628  nucleus->setMomentum(theCNMomentum);
629  nucleus->setEnergy(theCNEnergy);
630  nucleus->setExcitationEnergy(theCNExcitationEnergy);
631  nucleus->setMass(theCNMass+theCNExcitationEnergy);
632  nucleus->setSpin(theCNSpin); // neglects any orbital angular momentum of the CN
633 
634  // Take care of any remaining deltas
636 
637  // Take care of any remaining etas and/or omegas
638  G4double timeThreshold=theConfig->getDecayTimeThreshold();
640 
641  // Take care of any remaining Kaons
643 
644  // Cluster decay
646 
647  // Fill the EventInfo structure
649  }
650  }
651 
653  RecoilCMFunctor theRecoilFunctor(nucleus, theEventInfo);
654 
655  // Apply the root-finding algorithm
656  const RootFinder::Solution theSolution = RootFinder::solve(&theRecoilFunctor, 1.0);
657  if(theSolution.success) {
658  theRecoilFunctor(theSolution.x); // Apply the solution
659  } else {
660  INCL_WARN("Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." << '\n');
661  }
662 
663  }
664 
665 #ifndef INCLXX_IN_GEANT4_MODE
666  void INCL::globalConservationChecks(G4bool afterRecoil) {
668 
669  // Global conservation checks
670  const G4double pLongBalance = theBalance.momentum.getZ();
671  const G4double pTransBalance = theBalance.momentum.perp();
672  if(theBalance.Z != 0) {
673  INCL_ERROR("Violation of charge conservation! ZBalance = " << theBalance.Z << " eventNumber=" << theEventInfo.eventNumber << '\n');
674  }
675  if(theBalance.A != 0) {
676  INCL_ERROR("Violation of baryon-number conservation! ABalance = " << theBalance.A << " Emit Lambda=" << theEventInfo.emitLambda << " eventNumber=" << theEventInfo.eventNumber << '\n');
677  }
678  if(theBalance.S != 0) {
679  INCL_ERROR("Violation of strange-number conservation! SBalance = " << theBalance.S << " eventNumber=" << theEventInfo.eventNumber << '\n');
680  }
681  G4double EThreshold, pLongThreshold, pTransThreshold;
682  if(afterRecoil) {
683  // Less stringent checks after accommodating recoil
684  EThreshold = 10.; // MeV
685  pLongThreshold = 1.; // MeV/c
686  pTransThreshold = 1.; // MeV/c
687  } else {
688  // More stringent checks before accommodating recoil
689  EThreshold = 0.1; // MeV
690  pLongThreshold = 0.1; // MeV/c
691  pTransThreshold = 0.1; // MeV/c
692  }
693  if(std::abs(theBalance.energy)>EThreshold) {
694  INCL_WARN("Violation of energy conservation > " << EThreshold << " MeV. EBalance = " << theBalance.energy << " Emit Lambda=" << theEventInfo.emitLambda << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
695  }
696  if(std::abs(pLongBalance)>pLongThreshold) {
697  INCL_WARN("Violation of longitudinal momentum conservation > " << pLongThreshold << " MeV/c. pLongBalance = " << pLongBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
698  }
699  if(std::abs(pTransBalance)>pTransThreshold) {
700  INCL_WARN("Violation of transverse momentum conservation > " << pTransThreshold << " MeV/c. pTransBalance = " << pTransBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
701  }
702 
703  // Feed the EventInfo variables
704  theEventInfo.EBalance = theBalance.energy;
705  theEventInfo.pLongBalance = pLongBalance;
706  theEventInfo.pTransBalance = pTransBalance;
707  }
708 #endif
709 
711  // Stop if we have passed the stopping time
713  INCL_DEBUG("Cascade time (" << propagationModel->getCurrentTime()
714  << ") exceeded stopping time (" << propagationModel->getStoppingTime()
715  << "), stopping cascade" << '\n');
716  return false;
717  }
718  // Stop if there are no participants and no pions inside the nucleus
719  if(nucleus->getStore()->getBook().getCascading()==0 &&
720  nucleus->getStore()->getIncomingParticles().empty()) {
721  INCL_DEBUG("No participants in the nucleus and no incoming particles left, stopping cascade" << '\n');
722  return false;
723  }
724  // Stop if the remnant is smaller than minRemnantSize
725  if(nucleus->getA() <= minRemnantSize) {
726  INCL_DEBUG("Remnant size (" << nucleus->getA()
727  << ") smaller than or equal to minimum (" << minRemnantSize
728  << "), stopping cascade" << '\n');
729  return false;
730  }
731  // Stop if we have to try and make a compound nucleus or if we have to
732  // force a transparent
734  INCL_DEBUG("Trying to make a compound nucleus, stopping cascade" << '\n');
735  return false;
736  }
737 
738  return true;
739  }
740 
741  void INCL::finalizeGlobalInfo(Random::SeedVector const &initialSeeds) {
742  const G4double normalisationFactor = theGlobalInfo.geometricCrossSection /
744  theGlobalInfo.nucleonAbsorptionCrossSection = normalisationFactor *
746  theGlobalInfo.pionAbsorptionCrossSection = normalisationFactor *
748  theGlobalInfo.reactionCrossSection = normalisationFactor *
750  theGlobalInfo.errorReactionCrossSection = normalisationFactor *
752  theGlobalInfo.forcedCNCrossSection = normalisationFactor *
754  theGlobalInfo.errorForcedCNCrossSection = normalisationFactor *
756  theGlobalInfo.completeFusionCrossSection = normalisationFactor *
758  theGlobalInfo.errorCompleteFusionCrossSection = normalisationFactor *
759  std::sqrt((G4double) (theGlobalInfo.nCompleteFusion));
762 
763  theGlobalInfo.initialRandomSeeds.assign(initialSeeds.begin(), initialSeeds.end());
764 
766  theGlobalInfo.finalRandomSeeds.assign(theSeeds.begin(), theSeeds.end());
767  }
768 
770  // Do nothing if this is not a nucleus-nucleus reaction
772  return 0;
773 
774  // Get the spectators (geometrical+dynamical) from the Store
777 
778  G4int nUnmergedSpectators = 0;
779 
780  // If there are no spectators, do nothing
781  if(dynSpectators.empty() && geomSpectators.empty()) {
782  return 0;
783  } else if(dynSpectators.size()==1 && geomSpectators.empty()) {
784  // No geometrical spectators, one dynamical spectator
785  // Just put it back in the outgoing list
786  nucleus->getStore()->addToOutgoing(dynSpectators.front());
787  } else {
788  // Make a cluster out of the geometrical spectators
789  ProjectileRemnant *theProjectileRemnant = nucleus->getProjectileRemnant();
790 
791  // Add the dynamical spectators to the bunch
792  ParticleList rejected = theProjectileRemnant->addAllDynamicalSpectators(dynSpectators);
793  // Put back the rejected spectators into the outgoing list
794  nUnmergedSpectators = rejected.size();
795  nucleus->getStore()->addToOutgoing(rejected);
796 
797  // Deal with the projectile remnant
799 
800  }
801 
802  return nUnmergedSpectators;
803  }
804 
805  void INCL::initMaxInteractionDistance(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy) {
806  if(projectileSpecies.theType != Composite) {
808  return;
809  }
810 
813 
814  const G4double theNNDistance = CrossSections::interactionDistanceNN(projectileSpecies, kineticEnergy);
815  maxInteractionDistance = r0 + theNNDistance;
816  INCL_DEBUG("Initialised interaction distance: r0 = " << r0 << '\n'
817  << " theNNDistance = " << theNNDistance << '\n'
818  << " maxInteractionDistance = " << maxInteractionDistance << '\n');
819  }
820 
821  void INCL::initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z) {
822  G4double rMax = 0.0;
823  if(A==0) {
824  IsotopicDistribution const &anIsotopicDistribution =
826  IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
827  for(IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
828  const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, i->theA, Z);
829  const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, i->theA, Z);
830  const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
831  rMax = std::max(maximumRadius, rMax);
832  }
833  } else {
834  const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
835  const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
836  const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
837  rMax = std::max(maximumRadius, rMax);
838  }
839  if(p.theType==Composite || p.theType==Proton || p.theType==Neutron) {
842  } else if(p.theType==PiPlus
843  || p.theType==PiZero
844  || p.theType==PiMinus) {
847  } else if(p.theType==KPlus
848  || p.theType==KZero) {
851  } else if(p.theType==KZeroBar
852  || p.theType==KMinus) {
855  } else if(p.theType==Lambda
856  ||p.theType==SigmaPlus
857  || p.theType==SigmaZero
858  || p.theType==SigmaMinus) {
861  }
862  INCL_DEBUG("Initialised universe radius: " << maxUniverseRadius << '\n');
863  }
864 
866  // Increment the global counter for the number of shots
868 
870  // Increment the global counter for the number of transparents
872  // Increment the global counter for the number of forced transparents
873  if(forceTransparent)
875  return;
876  }
877 
878  // Check if we have an absorption:
881 
882  // Count complete-fusion events
884 
887 
888  // Counters for the number of violations of energy conservation in
889  // collisions
891  }
892 
893 }