Geant4  10.01.p03
G4PreCompoundModel.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 // $Id: G4PreCompoundModel.cc 90591 2015-06-04 13:45:29Z gcosmo $
27 //
28 // by V. Lara
29 //
30 // Modified:
31 // 01.04.2008 J.M.Quesada Several changes. Soft cut-off switched off.
32 // 01.05.2008 J.M.Quesada Protection against non-physical preeq.
33 // transitional regime has been set
34 // 03.09.2008 J.M.Quesada for external choice of inverse cross section option
35 // 06.09.2008 J.M.Quesada Also external choices have been added for:
36 // - superimposed Coulomb barrier (useSICB=true)
37 // - "never go back" hipothesis (useNGB=true)
38 // - soft cutoff from preeq. to equlibrium (useSCO=true)
39 // - CEM transition probabilities (useCEMtr=true)
40 // 20.08.2010 V.Ivanchenko Cleanup of the code:
41 // - integer Z and A;
42 // - emission and transition classes created at
43 // initialisation
44 // - options are set at initialisation
45 // - do not use copy-constructors for G4Fragment
46 // 03.01.2012 V.Ivanchenko Added pointer to G4ExcitationHandler to the
47 // constructor
48 
49 #include "G4PreCompoundModel.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4PreCompoundEmission.hh"
54 #include "G4GNASHTransitions.hh"
55 #include "G4ParticleDefinition.hh"
56 #include "G4Proton.hh"
57 #include "G4Neutron.hh"
58 
59 #include "G4NucleiProperties.hh"
61 #include "Randomize.hh"
62 #include "G4DynamicParticle.hh"
63 #include "G4ParticleTypes.hh"
64 #include "G4ParticleTable.hh"
65 #include "G4LorentzVector.hh"
66 #include "G4Exp.hh"
67 
69 
71  : G4VPreCompoundModel(ptr,"PRECO"), useHETCEmission(false),
72  useGNASHTransition(false), OPTxs(3), useSICB(false),
73  useNGB(false), useSCO(false), useCEMtr(true), maxZ(3), maxA(5)
74 {
75  if(!ptr) { SetExcitationHandler(new G4ExcitationHandler()); }
77 
78  // 12/pi2 factor is used in real computation
79  fLevelDensity = param.GetLevelDensity()*12.0/CLHEP::pi2;
80 
83  else { theEmission->SetDefaultModel(); }
86 
88  else { theTransition = new G4PreCompoundTransitions(); }
91 
94 }
95 
97 
99 {
100  delete theEmission;
101  delete theTransition;
102  delete GetExcitationHandler();
103 }
104 
106 
109  G4Nucleus & theNucleus)
110 {
111  const G4ParticleDefinition* primary = thePrimary.GetDefinition();
112  if(primary != neutron && primary != proton) {
114  ed << "G4PreCompoundModel is used for ";
115  if(primary) { ed << primary->GetParticleName(); }
116  G4Exception("G4PreCompoundModel::ApplyYourself()","had0033",FatalException,
117  ed,"");
118  return 0;
119  }
120 
121  G4int Zp = 0;
122  G4int Ap = 1;
123  if(primary == proton) { Zp = 1; }
124 
125  G4int A = theNucleus.GetA_asInt();
126  G4int Z = theNucleus.GetZ_asInt();
127 
128  //G4cout << "### G4PreCompoundModel::ApplyYourself: A= " << A << " Z= " << Z
129  // << " Ap= " << Ap << " Zp= " << Zp << G4endl;
130  // 4-Momentum
131  G4LorentzVector p = thePrimary.Get4Momentum();
133  p += G4LorentzVector(0.0,0.0,0.0,mass);
134  //G4cout << "Primary 4-mom " << p << " mass= " << mass << G4endl;
135 
136  // prepare fragment
137  G4Fragment anInitialState(A + Ap, Z + Zp, p);
138  anInitialState.SetNumberOfExcitedParticle(2, 1);
139  anInitialState.SetNumberOfHoles(1,0);
140  anInitialState.SetCreationTime(thePrimary.GetGlobalTime());
141 
142  // call excitation handler
143  G4ReactionProductVector * result = DeExcite(anInitialState);
144 
145  // fill particle change
146  theResult.Clear();
148  for(G4ReactionProductVector::iterator i= result->begin();
149  i != result->end(); ++i)
150  {
151  G4DynamicParticle * aNew =
152  new G4DynamicParticle((*i)->GetDefinition(),
153  (*i)->GetTotalEnergy(),
154  (*i)->GetMomentum());
155  delete (*i);
156  theResult.AddSecondary(aNew);
157  }
158  delete result;
159 
160  //return the filled particle change
161  return &theResult;
162 }
163 
165 
167 {
169  G4double Eex = aFragment.GetExcitationEnergy();
170  G4int Z = aFragment.GetZ_asInt();
171  G4int A = aFragment.GetA_asInt();
172 
173  //G4cout << "### G4PreCompoundModel::DeExcite" << G4endl;
174  //G4cout << aFragment << G4endl;
175 
176  // Perform Equilibrium Emission
177  if ((Z < maxZ && A < maxA) || Eex < MeV /*|| Eex > 3.*MeV*A*/) {
178  PerformEquilibriumEmission(aFragment, Result);
179  return Result;
180  }
181 
182  // main loop
183  G4int count = 0;
184  const G4int countmax = 10000;
185  for (;;) {
186  //G4cout << "### PreCompound loop over fragment" << G4endl;
187  //G4cout << aFragment << G4endl;
188 
189  theEmission->Initialize(aFragment);
190 
191  G4int EquilibriumExcitonNumber =
192  G4lrint(std::sqrt(aFragment.GetExcitationEnergy()
193  *aFragment.GetA_asInt()*fLevelDensity));
194  //
195  // G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
196  //
197  // J. M. Quesada (Jan. 08) equilibrium hole number could be used as preeq.
198  // evap. delimiter (IAEA report)
199 
200  // Loop for transitions, it is performed while there are
201  // preequilibrium transitions.
202  G4bool ThereIsTransition = false;
203 
204  // G4cout<<"----------------------------------------"<<G4endl;
205  // G4double NP=aFragment.GetNumberOfParticles();
206  // G4double NH=aFragment.GetNumberOfHoles();
207  // G4double NE=aFragment.GetNumberOfExcitons();
208  // G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
209  // G4cout<<"N. excitons="<<NE<<" N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
210  do {
211  ++count;
212  //G4cout<<"transition number .."<<count
213  // <<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
214  G4bool go_ahead = false;
215  // soft cutoff criterium as an "ad-hoc" solution to force
216  // increase in evaporation
217  G4int test = aFragment.GetNumberOfExcitons();
218  if (test <= EquilibriumExcitonNumber) { go_ahead=true; }
219 
220  //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
221  if (useSCO && go_ahead)
222  {
223  G4double x = G4double(test)/G4double(EquilibriumExcitonNumber) - 1;
224  if( G4UniformRand() < 1.0 - G4Exp(-x*x/0.32) ) { go_ahead = false; }
225  }
226 
227  // JMQ: WARNING: CalculateProbability MUST be called prior to Get!!
228  // (O values would be returned otherwise)
229  G4double TotalTransitionProbability =
234  //G4cout<<"#0 P1="<<P1<<" P2="<<P2<<" P3="<<P3<<G4endl;
235 
236  //J.M. Quesada (May 2008) Physical criterium (lamdas) PREVAILS over
237  // approximation (critical exciton number)
238  //V.Ivanchenko (May 2011) added check on number of nucleons
239  // to send a fragment to FermiBreakUp
240  if(!go_ahead || P1 <= P2+P3 ||
241  (aFragment.GetZ_asInt() < maxZ && aFragment.GetA_asInt() < maxA) )
242  {
243  //G4cout<<"#4 EquilibriumEmission"<<G4endl;
244  PerformEquilibriumEmission(aFragment,Result);
245  return Result;
246  }
247  else
248  {
249  G4double TotalEmissionProbability =
250  theEmission->GetTotalProbability(aFragment);
251  //
252  //G4cout<<"#1 TotalEmissionProbability="<<TotalEmissionProbability
253  // <<" Nex= " <<aFragment.GetNumberOfExcitons()<<G4endl;
254  //
255  // Check if number of excitons is greater than 0
256  // else perform equilibrium emission
257  if (aFragment.GetNumberOfExcitons() <= 0)
258  {
259  PerformEquilibriumEmission(aFragment,Result);
260  return Result;
261  }
262 
263  //J.M.Quesada (May 08) this has already been done in order to decide
264  // what to do (preeq-eq)
265  // Sum of all probabilities
266  G4double TotalProbability = TotalEmissionProbability
267  + TotalTransitionProbability;
268 
269  // Select subprocess
270  if (TotalProbability*G4UniformRand() > TotalEmissionProbability)
271  {
272  //G4cout<<"#2 Transition"<<G4endl;
273  // It will be transition to state with a new number of excitons
274  ThereIsTransition = true;
275  // Perform the transition
276  theTransition->PerformTransition(aFragment);
277  }
278  else
279  {
280  //G4cout<<"#3 Emission"<<G4endl;
281  // It will be fragment emission
282  ThereIsTransition = false;
283  Result->push_back(theEmission->PerformEmission(aFragment));
284  }
285  }
286  } while (ThereIsTransition); // end of do loop
287  if(count >= countmax) {
289  ed << "G4PreCompoundModel loop over " << countmax << " iterations; "
290  << "current G4Fragment: \n" << aFragment;
291  G4Exception("G4PreCompoundModel::DeExcite()","had0034",JustWarning,
292  ed,"");
293  count = 0;
294  }
295  } // end of for (;;) loop
296  return Result;
297 }
298 
300 // Initialisation
302 
304 {
305  useHETCEmission = true;
307 }
308 
310 {
311  useHETCEmission = false;
313 }
314 
316  useGNASHTransition = true;
317  delete theTransition;
321 }
322 
324  useGNASHTransition = false;
325  delete theTransition;
329 }
330 
332 {
333  OPTxs = opt;
335 }
336 
338 {
339  useSICB = true;
341 }
342 
344 {
345  useNGB = true;
346 }
347 
349 {
350  useSCO = true;
351 }
352 
354 {
355  useCEMtr = true;
356 }
357 
359 // Documentation
361 
362 void G4PreCompoundModel::ModelDescription(std::ostream& outFile) const
363 {
364  outFile << "The GEANT4 precompound model is considered as an extension of the\n"
365  << "hadron kinetic model. It gives a possibility to extend the low energy range\n"
366  << "of the hadron kinetic model for nucleon-nucleus inelastic collision and it \n"
367  << "provides a ”smooth” transition from kinetic stage of reaction described by the\n"
368  << "hadron kinetic model to the equilibrium stage of reaction described by the\n"
369  << "equilibrium deexcitation models.\n"
370  << "The initial information for calculation of pre-compound nuclear stage\n"
371  << "consists of the atomic mass number A, charge Z of residual nucleus, its\n"
372  << "four momentum P0 , excitation energy U and number of excitons n, which equals\n"
373  << "the sum of the number of particles p (from them p_Z are charged) and the number of\n"
374  << "holes h.\n"
375  << "At the preequilibrium stage of reaction, we follow the exciton model approach in ref. [1],\n"
376  << "taking into account the competition among all possible nuclear transitions\n"
377  << "with ∆n = +2, −2, 0 (which are defined by their associated transition probabilities) and\n"
378  << "the emission of neutrons, protons, deutrons, thritium and helium nuclei (also defined by\n"
379  << "their associated emission probabilities according to exciton model)\n"
380  << "\n"
381  << "[1] K.K. Gudima, S.G. Mashnik, V.D. Toneev, Nucl. Phys. A401 329 (1983)\n"
382  << "\n";
383 }
384 
385 void G4PreCompoundModel::DeExciteModelDescription(std::ostream& outFile) const
386 {
387  outFile << "description of precompound model as used with DeExcite()"
388  << "\n";
389 }
390 
virtual void PerformTransition(G4Fragment &aFragment)=0
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
void PerformEquilibriumEmission(const G4Fragment &aFragment, G4ReactionProductVector *theResult) const
static const double MeV
Definition: G4SIunits.hh:193
static G4double GetNuclearMass(const G4double A, const G4double Z)
static const G4double * P1[nN]
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4PreCompoundEmission * theEmission
G4PreCompoundModel(G4ExcitationHandler *ptr=0)
const G4ParticleDefinition * neutron
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:360
int G4int
Definition: G4Types.hh:78
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)
const G4String & GetParticleName() const
G4HadFinalState theResult
void SetStatusChange(G4HadFinalStateStatus aS)
void SetExcitationHandler(G4ExcitationHandler *ptr)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void SetNumberOfExcitedParticle(G4int valueTot, G4int valueP)
Definition: G4Fragment.hh:341
G4ExcitationHandler * GetExcitationHandler() const
#define G4UniformRand()
Definition: Randomize.hh:93
G4int GetA_asInt() const
Definition: G4Fragment.hh:243
const G4ParticleDefinition * GetDefinition() const
G4VPreCompoundTransitions * theTransition
bool G4bool
Definition: G4Types.hh:79
G4double GetGlobalTime() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static const G4double A[nN]
const G4LorentzVector & Get4Momentum() const
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual G4double CalculateProbability(const G4Fragment &aFragment)=0
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4ReactionProduct * PerformEmission(G4Fragment &aFragment)
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:418
G4int GetNumberOfExcitons() const
Definition: G4Fragment.hh:325
int G4lrint(double ad)
Definition: templates.hh:163
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
virtual void ModelDescription(std::ostream &outFile) const
G4int GetZ_asInt() const
Definition: G4Fragment.hh:248
static const G4double * P2[nN]
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
const G4ParticleDefinition * proton
double G4double
Definition: G4Types.hh:76
virtual void DeExciteModelDescription(std::ostream &outFile) const
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:260
void SetOPTxs(G4int opt)
void Initialize(const G4Fragment &aFragment)
G4double GetTotalProbability(const G4Fragment &aFragment)
CLHEP::HepLorentzVector G4LorentzVector