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