Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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;
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  // Forced CN?
380  INCL_DEBUG("Trying compound nucleus" << '\n');
383  // Global checks of conservation laws
384 #ifndef INCLXX_IN_GEANT4_MODE
385  if(!theEventInfo.transparent) globalConservationChecks(true);
386 #endif
387  return;
388  }
389 
391 
393  ProjectileRemnant * const projectileRemnant = nucleus->getProjectileRemnant();
394  if(projectileRemnant) {
395  // Clear the incoming list (particles will be deleted by the ProjectileRemnant)
397  } else {
398  // Delete particles in the incoming list
400  }
401  } else {
402 
403  // Check if the nucleus contains strange particles
408 
409  // Capture antiKaons and Sigmas and produce Lambda instead
411 
412  // Emit antiKaons and Sigmas still inside the nucleus
415  // Should be activated only for geant4
416 #ifdef INCLXX_IN_GEANT4_MODE
418 #endif // INCLXX_IN_GEANT4_MODE
419 
420  // The event bias
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;
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::random_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  break;
561  case PauliBlockedFS:
563  default:
564  success = false;
565  break;
566  }
567  }
568 
569  if(!success || !atLeastOneNucleonEntering) {
570  INCL_DEBUG("No nucleon entering in forced CN, forcing a transparent" << '\n');
571  forceTransparent = true;
572  return;
573  }
574 
575 // assert(theCNA==nucleus->getA());
576 // assert(theCNA<=theEventInfo.At+theEventInfo.Ap);
577 // assert(theCNZ<=theEventInfo.Zt+theEventInfo.Zp);
578 
579  // Update the kinematics of the CN
580  theCNEnergy -= theProjectileRemnant->getEnergy();
581  theCNMomentum -= theProjectileRemnant->getMomentum();
582 
583  // Deal with the projectile remnant
585 
586  // Subtract the angular momentum of the projectile remnant
587 // assert(nucleus->getStore()->getOutgoingParticles().empty());
588  theCNSpin -= theProjectileRemnant->getAngularMomentum();
589 
590  // Compute the excitation energy of the CN
591  const G4double theCNMass = ParticleTable::getTableMass(theCNA,theCNZ);
592  const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.mag2();
593  if(theCNInvariantMassSquared<0.) {
594  // Negative invariant mass squared, return a transparent
595  forceTransparent = true;
596  return;
597  }
598  const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
599  if(theCNExcitationEnergy<0.) {
600  // Negative excitation energy, return a transparent
601  INCL_DEBUG("CN excitation energy is negative, forcing a transparent" << '\n'
602  << " theCNA = " << theCNA << '\n'
603  << " theCNZ = " << theCNZ << '\n'
604  << " theCNEnergy = " << theCNEnergy << '\n'
605  << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
606  << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
607  << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
608  );
609  forceTransparent = true;
610  return;
611  } else {
612  // Positive excitation energy, can make a CN
613  INCL_DEBUG("CN excitation energy is positive, forcing a CN" << '\n'
614  << " theCNA = " << theCNA << '\n'
615  << " theCNZ = " << theCNZ << '\n'
616  << " theCNEnergy = " << theCNEnergy << '\n'
617  << " theCNMomentum = (" << theCNMomentum.getX() << ", "<< theCNMomentum.getY() << ", " << theCNMomentum.getZ() << ")" << '\n'
618  << " theCNExcitationEnergy = " << theCNExcitationEnergy << '\n'
619  << " theCNSpin = (" << theCNSpin.getX() << ", "<< theCNSpin.getY() << ", " << theCNSpin.getZ() << ")" << '\n'
620  );
621  nucleus->setA(theCNA);
622  nucleus->setZ(theCNZ);
623  nucleus->setMomentum(theCNMomentum);
624  nucleus->setEnergy(theCNEnergy);
625  nucleus->setExcitationEnergy(theCNExcitationEnergy);
626  nucleus->setMass(theCNMass+theCNExcitationEnergy);
627  nucleus->setSpin(theCNSpin); // neglects any orbital angular momentum of the CN
628 
629  // Take care of any remaining deltas
631 
632  // Take care of any remaining etas and/or omegas
633  G4double timeThreshold=theConfig->getDecayTimeThreshold();
635 
636  // Cluster decay
638 
639  // Fill the EventInfo structure
641  }
642  }
643 
645  RecoilCMFunctor theRecoilFunctor(nucleus, theEventInfo);
646 
647  // Apply the root-finding algorithm
648  const RootFinder::Solution theSolution = RootFinder::solve(&theRecoilFunctor, 1.0);
649  if(theSolution.success) {
650  theRecoilFunctor(theSolution.x); // Apply the solution
651  } else {
652  INCL_WARN("Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." << '\n');
653  }
654 
655  }
656 
657 #ifndef INCLXX_IN_GEANT4_MODE
658  void INCL::globalConservationChecks(G4bool afterRecoil) {
660 
661  // Global conservation checks
662  const G4double pLongBalance = theBalance.momentum.getZ();
663  const G4double pTransBalance = theBalance.momentum.perp();
664  if(theBalance.Z != 0) {
665  INCL_ERROR("Violation of charge conservation! ZBalance = " << theBalance.Z << '\n');
666  }
667  if(theBalance.A != 0) {
668  INCL_ERROR("Violation of baryon-number conservation! ABalance = " << theBalance.A << " Emit Lambda=" << theEventInfo.emitLambda << " eventNumber=" << theEventInfo.eventNumber << '\n');
669  }
670  G4double EThreshold, pLongThreshold, pTransThreshold;
671  if(afterRecoil) {
672  // Less stringent checks after accommodating recoil
673  EThreshold = 10.; // MeV
674  pLongThreshold = 1.; // MeV/c
675  pTransThreshold = 1.; // MeV/c
676  } else {
677  // More stringent checks before accommodating recoil
678  EThreshold = 0.1; // MeV
679  pLongThreshold = 0.1; // MeV/c
680  pTransThreshold = 0.1; // MeV/c
681  }
682  if(std::abs(theBalance.energy)>EThreshold) {
683  INCL_WARN("Violation of energy conservation > " << EThreshold << " MeV. EBalance = " << theBalance.energy << " Emit Lambda=" << theEventInfo.emitLambda << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
684  }
685  if(std::abs(pLongBalance)>pLongThreshold) {
686  INCL_WARN("Violation of longitudinal momentum conservation > " << pLongThreshold << " MeV/c. pLongBalance = " << pLongBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
687  }
688  if(std::abs(pTransBalance)>pTransThreshold) {
689  INCL_WARN("Violation of transverse momentum conservation > " << pTransThreshold << " MeV/c. pTransBalance = " << pTransBalance << " afterRecoil = " << afterRecoil << " eventNumber=" << theEventInfo.eventNumber << '\n');
690  }
691 
692  // Feed the EventInfo variables
693  theEventInfo.EBalance = theBalance.energy;
694  theEventInfo.pLongBalance = pLongBalance;
695  theEventInfo.pTransBalance = pTransBalance;
696  }
697 #endif
698 
700  // Stop if we have passed the stopping time
702  INCL_DEBUG("Cascade time (" << propagationModel->getCurrentTime()
703  << ") exceeded stopping time (" << propagationModel->getStoppingTime()
704  << "), stopping cascade" << '\n');
705  return false;
706  }
707  // Stop if there are no participants and no pions inside the nucleus
708  if(nucleus->getStore()->getBook().getCascading()==0 &&
709  nucleus->getStore()->getIncomingParticles().empty()) {
710  INCL_DEBUG("No participants in the nucleus and no incoming particles left, stopping cascade" << '\n');
711  return false;
712  }
713  // Stop if the remnant is smaller than minRemnantSize
714  if(nucleus->getA() <= minRemnantSize) {
715  INCL_DEBUG("Remnant size (" << nucleus->getA()
716  << ") smaller than or equal to minimum (" << minRemnantSize
717  << "), stopping cascade" << '\n');
718  return false;
719  }
720  // Stop if we have to try and make a compound nucleus or if we have to
721  // force a transparent
723  INCL_DEBUG("Trying to make a compound nucleus, stopping cascade" << '\n');
724  return false;
725  }
726 
727  return true;
728  }
729 
730  void INCL::finalizeGlobalInfo(Random::SeedVector const &initialSeeds) {
731  const G4double normalisationFactor = theGlobalInfo.geometricCrossSection /
733  theGlobalInfo.nucleonAbsorptionCrossSection = normalisationFactor *
735  theGlobalInfo.pionAbsorptionCrossSection = normalisationFactor *
737  theGlobalInfo.reactionCrossSection = normalisationFactor *
739  theGlobalInfo.errorReactionCrossSection = normalisationFactor *
741  theGlobalInfo.forcedCNCrossSection = normalisationFactor *
743  theGlobalInfo.errorForcedCNCrossSection = normalisationFactor *
745  theGlobalInfo.completeFusionCrossSection = normalisationFactor *
747  theGlobalInfo.errorCompleteFusionCrossSection = normalisationFactor *
748  std::sqrt((G4double) (theGlobalInfo.nCompleteFusion));
751 
752  theGlobalInfo.initialRandomSeeds.assign(initialSeeds.begin(), initialSeeds.end());
753 
755  theGlobalInfo.finalRandomSeeds.assign(theSeeds.begin(), theSeeds.end());
756  }
757 
759  // Do nothing if this is not a nucleus-nucleus reaction
761  return 0;
762 
763  // Get the spectators (geometrical+dynamical) from the Store
766 
767  G4int nUnmergedSpectators = 0;
768 
769  // If there are no spectators, do nothing
770  if(dynSpectators.empty() && geomSpectators.empty()) {
771  return 0;
772  } else if(dynSpectators.size()==1 && geomSpectators.empty()) {
773  // No geometrical spectators, one dynamical spectator
774  // Just put it back in the outgoing list
775  nucleus->getStore()->addToOutgoing(dynSpectators.front());
776  } else {
777  // Make a cluster out of the geometrical spectators
778  ProjectileRemnant *theProjectileRemnant = nucleus->getProjectileRemnant();
779 
780  // Add the dynamical spectators to the bunch
781  ParticleList rejected = theProjectileRemnant->addAllDynamicalSpectators(dynSpectators);
782  // Put back the rejected spectators into the outgoing list
783  nUnmergedSpectators = rejected.size();
784  nucleus->getStore()->addToOutgoing(rejected);
785 
786  // Deal with the projectile remnant
788 
789  }
790 
791  return nUnmergedSpectators;
792  }
793 
795  if(projectileSpecies.theType != Composite) {
797  return;
798  }
799 
802 
803  const G4double theNNDistance = CrossSections::interactionDistanceNN(projectileSpecies, kineticEnergy);
804  maxInteractionDistance = r0 + theNNDistance;
805  INCL_DEBUG("Initialised interaction distance: r0 = " << r0 << '\n'
806  << " theNNDistance = " << theNNDistance << '\n'
807  << " maxInteractionDistance = " << maxInteractionDistance << '\n');
808  }
809 
811  G4double rMax = 0.0;
812  if(A==0) {
813  IsotopicDistribution const &anIsotopicDistribution =
815  IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
816  for(IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
817  const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, i->theA, Z);
818  const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, i->theA, Z);
819  const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
820  rMax = std::max(maximumRadius, rMax);
821  }
822  } else {
823  const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
824  const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
825  const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
826  rMax = std::max(maximumRadius, rMax);
827  }
828  if(p.theType==Composite || p.theType==Proton || p.theType==Neutron) {
831  } else if(p.theType==PiPlus
832  || p.theType==PiZero
833  || p.theType==PiMinus) {
836  } else if(p.theType==KPlus
837  || p.theType==KZero) {
840  } else if(p.theType==KZeroBar
841  || p.theType==KMinus) {
844  } else if(p.theType==Lambda
845  ||p.theType==SigmaPlus
846  || p.theType==SigmaZero
847  || p.theType==SigmaMinus) {
850  }
851  INCL_DEBUG("Initialised universe radius: " << maxUniverseRadius << '\n');
852  }
853 
855  // Increment the global counter for the number of shots
857 
859  // Increment the global counter for the number of transparents
861  // Increment the global counter for the number of forced transparents
862  if(forceTransparent)
864  return;
865  }
866 
867  // Check if we have an absorption:
870 
871  // Count complete-fusion events
873 
876 
877  // Counters for the number of violations of energy conservation in
878  // collisions
880  }
881 
882 }
void afterCascadeAction(Nucleus *)
void finalizeGlobalInfo(Random::SeedVector const &initialSeeds)
virtual G4INCL::IAvatar * propagate(FinalState const *const fs)=0
Float_t effectiveImpactParameter
Effective (Coulomb-distorted) impact parameter [fm].
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual void setNucleus(G4INCL::Nucleus *nucleus)=0
void setZ(const G4int Z)
Set the charge number of the cluster.
Intersection-point structure.
G4bool decayInsideDeltas()
Force the decay of deltas inside the nucleus.
G4double Double_t
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
void beforeAvatarAction(IAvatar *a, Nucleus *n)
G4double interactionDistanceKN(const G4double projectileKineticEnergy)
Compute the &quot;interaction distance&quot;.
G4double getProjectileKineticEnergy() const
Get the projectile kinetic energy.
G4double computeExcitationEnergy() const
Compute the current excitation energy.
G4int emitInsideLambda()
Force emission of all Lambda (desexitation code with strangeness not implanted yet) ...
G4int minRemnantSize
Remnant size below which cascade stops.
void finalizeProjectileRemnant(const G4double emissionTime)
Finalise the projectile remnant.
G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy)
Initialise the cascade.
G4double getEnergy() const
G4double shoot0()
void deleteClusteringModel()
Delete the clustering model.
CascadeActionType getCascadeActionType() const
Get the cascade-action type.
static std::vector< G4double > INCLBiasVector
Time ordered vector of all bias applied.
Nucleus * nucleus
static void setCutNN(const G4double c)
Int_t nTransparents
Number of transparent shots.
G4bool containsKaon()
Returns true if the nucleus contains any Kaons.
Bool_t kaonsInside
Event involved kaons in the nucleus at the end of the cascade.
virtual G4double shoot(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)=0
Short_t St
Strangeness number of the target nucleus.
void initialize(Config const *const)
Initialize generator according to a Config object.
CascadeAction * cascadeAction
Float_t errorCompleteFusionCrossSection
Error on the calculated complete-fusion cross section (nParticles==0)
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
void setMass(G4double mass)
INCL(Config const *const config)
void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z)
Initialize the universe radius.
static std::string const getVersionString()
Get the INCL version string.
const G4INCL::ThreeVector & getMomentum() const
LocalEnergyType getLocalEnergyBBType() const
Get the type of local energy for N-N avatars.
G4int getS() const
Returns the strangeness number.
Short_t At
Target mass number given as input.
G4bool decayInsideStrangeParticles()
Force the transformation of strange particles into a Lambda;.
void makeCompoundNucleus()
Make a compound nucleus.
ParticleList addAllDynamicalSpectators(ParticleList const &pL)
Add back all dynamical spectators to the projectile remnant.
G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S)
IsotopeVector const & getIsotopes() const
Get the isotope vector.
const char * p
Definition: xmltok.h:285
FinalStateValidity getValidity() const
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
Float_t errorReactionCrossSection
Error on the calculated reaction cross section.
G4double getDecayTimeThreshold() const
Get the decay time threshold time.
Static class for cluster formation.
void reset()
Reset the EventInfo members.
Book & getBook()
Definition: G4INCLStore.hh:259
static void deleteBackupParticles()
Release the memory allocated for the backup particles.
Short_t Zt
Target charge number given as input.
G4double fixedImpactParameter
G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n)
Return the maximum impact parameter for Coulomb-distorted trajectories.
Float_t reactionCrossSection
Calculated reaction cross section.
G4double interactionDistanceYN(const G4double projectileKineticEnergy)
Compute the &quot;interaction distance&quot;.
void initialize(Config const *const aConfig)
Initialise blockers according to a Config object.
Bool_t sigmasInside
Event involved sigmas in the nucleus at the end of the cascade.
G4bool forceTransparent
Functions that encapsulate a mass table.
void applyFinalState(FinalState *)
G4int getTargetZ() const
Get the target charge number.
Definition: G4INCLConfig.hh:97
Float_t Ep
Projectile kinetic energy given as input.
Float_t completeFusionCrossSection
Calculated complete-fusion cross section (nParticles==0)
void rescaleOutgoingForRecoil()
Rescale the energies of the outgoing particles.
G4double getImpactParameter() const
G4double maxInteractionDistance
G4int drawRandomNaturalIsotope(const G4int Z)
FinalState * getFinalState()
std::string getDeExcitationString() const
Get the de-excitation string.
void deleteGenerator()
G4bool hasRemnant() const
Does the nucleus give a cascade remnant?
G4double getHadronizationTime() const
Get the hadronization time.
double S(double temp)
Int_t nForcedTransparents
Number of forced transparents.
ParticleList const & getIncomingParticles() const
Definition: G4INCLStore.hh:217
void initVerbosityLevelFromEnvvar()
Alternative CascadeAction for dumping avatars.
void fillFinalState(FinalState *fs)
#define INCL_DEBUG(x)
void addParticleEntryAvatar(IAvatar *a)
Add one ParticleEntry avatar.
Definition: G4INCLStore.cc:66
void initializeParticles()
G4double interactionDistanceNN(const ParticleSpecies &aSpecies, const G4double kineticEnergy)
Compute the &quot;interaction distance&quot;.
G4double getZ() const
Float_t pLongBalance
Longitudinal momentum-conservation balance [MeV/c].
Float_t Z
G4int getTargetA() const
Get the target mass number.
Definition: G4INCLConfig.hh:94
static G4double getTotalBias()
General bias vector function.
void beforePropagationAction(IPropagationModel *pm)
Int_t nCompleteFusion
Number of complete-fusion events (nParticles==0)
void fillEventInfo(EventInfo *eventInfo)
Config const *const theConfig
Int_t nUnmergedSpectators
Number of dynamical spectators that were merged back into the projectile remnant. ...
Adapter const & getAdapter()
ParticleSpecies getProjectileSpecies() const
Get the projectile species.
ParticleList const & getParticles() const
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
void deleteCoulomb()
Delete the Coulomb-distortion object.
Int_t nShots
Number of shots.
Bool_t pionAbsorption
True if the event is a pion absorption.
Class that stores isotopic abundances for a given element.
const G4double tenPi
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
G4bool targetInitSuccess
void cascade()
The actual cascade loop.
G4double shoot()
Definition: G4INCLRandom.cc:93
Float_t impactParameter
Impact parameter [fm].
Solution solve(RootFunctor const *const f, const G4double x0)
Numerically solve a one-dimensional equation.
void reset()
Definition: G4INCLBook.hh:52
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
Float_t EBalance
Energy-conservation balance [MeV].
Short_t Sp
Strangeness number of the projectile nucleus.
SeedVector getSeeds()
Definition: G4INCLRandom.cc:89
Int_t nPionAbsorptions
Number of nucleon absorptions (no outcoming pions)
#define INCL_ERROR(x)
void useFusionKinematics()
Adjust the kinematics for complete-fusion events.
G4int makeProjectileRemnant()
Make a projectile pre-fragment out of geometrical spectators.
void initialize(Config const *const theConfig)
Initialize the clustering model based on the Config object.
Short_t Zt
Charge number of the target nucleus.
Class to adjust remnant recoil in the reaction CM system.
Bool_t forcedDeltasOutside
Event involved forced delta decays outside the nucleus.
void afterPropagationAction(IPropagationModel *pm, IAvatar *avatar)
double A(double temperature)
void deleteBlockers()
Delete blockers.
G4int getZ() const
Returns the charge number.
void setSpin(const ThreeVector &j)
Set the spin of the nucleus.
Bool_t antikaonsInside
Event involved antikaons in the nucleus at the end of the cascade.
G4double getX() const
void clearIncoming()
Clear the incoming list.
Definition: G4INCLStore.hh:145
void clearOutgoing()
Definition: G4INCLStore.cc:223
void beforeRunAction(Config const *config)
G4int getA() const
Returns the baryon number.
G4double perp() const
LocalEnergyType getLocalEnergyPiType() const
Get the type of local energy for pi-N and decay avatars.
G4bool decayMe()
Force the phase-space decay of the Nucleus.
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
void emitInsideStrangeParticles()
Force emission of all strange particles inside the nucleus.
Float_t eventBias
Event bias.
G4bool isEventTransparent() const
Is the event transparent?
G4bool containsSigma()
Returns true if the nucleus contains any Sigma.
virtual G4double getStoppingTime()=0
void reset()
Reset the projectile remnant to the state at the beginning of the cascade.
IsotopeVector::iterator IsotopeIter
std::vector< Int_t > initialRandomSeeds
Initial seeds for the pseudo-random-number generator.
G4bool emitInsideKaon()
Force emission of all Kaon inside the nucleus.
Simple container for output of calculation-wide results.
G4bool isNaturalTarget() const
Natural targets.
Definition: G4INCLConfig.hh:87
Float_t pTransBalance
Transverse momentum-conservation balance [MeV/c].
const G4double twoPi
Struct for conservation laws.
G4double getCutNN() const
static G4ThreadLocal Int_t eventNumber
Number of the event.
const EventInfo & processEvent()
Float_t Ep
Projectile kinetic energy given as input.
void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy)
Initialise the maximum interaction distance.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
Bool_t emitKaon
Event involved forced Kaon emission.
IPropagationModel * propagationModel
void updateGlobalInfo()
Update global counters and other members of theGlobalInfo object.
#define INCL_WARN(x)
int G4int
Definition: G4Types.hh:78
void deleteIncoming()
Clear the incoming list and delete the particles.
Definition: G4INCLStore.hh:150
G4bool decayOutgoingSigmaZero(G4double timeThreshold)
Force the decay of outgoing Neutral Sigma.
Bool_t forcedPionResonancesOutside
Event involved forced eta/omega decays outside the nucleus.
Simple class for computing intersections between a straight line and a sphere.
Int_t nForcedCompoundNucleus
Number of forced compound-nucleus events.
G4double interactionDistancePiN(const G4double projectileKineticEnergy)
Compute the &quot;interaction distance&quot;.
Short_t nCascadeParticles
Number of cascade particles.
void setExcitationEnergy(const G4double e)
Set the excitation energy of the cluster.
Bool_t nucleonAbsorption
True if the event is a nucleon absorption.
G4bool decayOutgoingDeltas()
Force the decay of outgoing deltas.
void initialize(Config const *const theConfig)
Initialize the Coulomb-distortion algorithm.
Int_t emitLambda
Number of forced Lambda emit out of the nucleus.
Float_t errorForcedCNCrossSection
Error on the calculated forced-compound-nucleus cross section.
GlobalInfo theGlobalInfo
G4bool containsDeltas()
Returns true if the nucleus contains any deltas.
Bool_t deltasInside
Event involved deltas in the nucleus at the end of the cascade.
G4double mag2() const
Short_t Ap
Mass number of the projectile nucleus.
void postCascade()
Finalise the cascade and clean up.
Float_t nucleonAbsorptionCrossSection
Nucleon absorption cross section.
G4double interactionDistanceKbarN(const G4double projectileKineticEnergy)
Compute the &quot;interaction distance&quot;.
void setEnergy(G4double energy)
void addToOutgoing(Particle *p)
add the particle to the outgoing particle list.
Definition: G4INCLStore.hh:190
std::string cascadeModel
Name of the cascade model.
Intersection getEarlierTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r)
Compute the first intersection of a straight particle trajectory with a sphere.
std::string deexcitationModel
Name of the de-excitation model.
Bool_t forcedDeltasInside
Event involved forced delta decays inside the nucleus.
void clearCache()
Clear the INuclearPotential cache.
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
Class containing default actions to be performed at intermediate cascade steps.
virtual G4double getCurrentTime()=0
void distortOut(ParticleList const &pL, Nucleus const *const n)
Modify the momentum of an outgoing particle.
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
static G4ThreadLocal G4int nextBiasedCollisionID
Float_t stoppingTime
Cascade stopping time [fm/c].
void initialize(Config const *const theConfig)
Short_t Zp
Charge number of the projectile nucleus.
std::vector< Int_t > finalRandomSeeds
Final seeds for the pseudo-random-number generator.
G4bool decayOutgoingNeutralKaon()
Force the transformation of outgoing Neutral Kaon into propation eigenstate.
Int_t projectileType
Projectile particle type.
G4bool containsAntiKaon()
Returns true if the nucleus contains any anti Kaons.
Float_t energyViolationInteractionCrossSection
Cross section for attempted collisions/decays for which the energy-conservation algorithm failed to f...
Static class for selecting Coulomb distortion.
Short_t Ap
Projectile mass number given as input.
Short_t Zp
Projectile charge number given as input.
G4bool containsLambda()
Returns true if the nucleus contains any Lambda.
Bool_t lambdasInside
Event involved lambdas in the nucleus at the end of the cascade.
Bool_t absorbedStrangeParticle
Event involved forced strange absorbtion inside the nucleus.
Float_t pionAbsorptionCrossSection
Pion absorption cross section.
void initialize(Config const *const theConfig)
ParticleList const & getOutgoingParticles() const
Definition: G4INCLStore.hh:223
G4double getBias() const
Get the bias.
void setA(const G4int A)
Set the mass number of the cluster.
ParticleList extractDynamicalSpectators()
Returns a list of dynamical spectators.
Definition: G4INCLStore.hh:232
Abstract interface to the nuclear potential.
G4double getExcitationEnergy() const
Get the excitation energy of the nucleus.
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
void beforeCascadeAction(IPropagationModel *)
std::vector< Isotope > IsotopeVector
G4double maxImpactParameter
Short_t At
Mass number of the target nucleus.
G4double getY() const
Store * getStore() const
void initialize(Config const *const theConfig=0)
Initialize the particle table.
const ThreeVector & getIncomingAngularMomentum() const
Get the incoming angular-momentum vector.
Float_t geometricCrossSection
Geometric cross section.
G4double maxUniverseRadius
G4int getCascading() const
Definition: G4INCLBook.hh:105
G4bool getTryCompoundNucleus()
G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
G4bool continueCascade()
Stopping criterion for the cascade.
Int_t nNucleonAbsorptions
Number of nucleon absorptions (no outcoming particles)
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
EventInfo theEventInfo
ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const
Compute charge, mass, energy and momentum balance.
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static void setBias(const G4double b)
Set the bias.
G4bool decayOutgoingPionResonances(G4double timeThreshold)
Force the decay of outgoing PionResonances (eta/omega).
Bool_t clusterDecay
Event involved cluster decay.
G4bool initializeTarget(const G4int A, const G4int Z, const G4int S)
void afterAvatarAction(IAvatar *a, Nucleus *n, FinalState *fs)
Bool_t transparent
True if the event is transparent.
Float_t forcedCNCrossSection
Calculated forced-compound-nucleus cross section.