Geant4  10.02.p01
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 91837 2015-08-07 07:27:08Z 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
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 = 1000;
185  for (;;) {
186  //G4cout << "### PreCompound loop over fragment" << G4endl;
187  //G4cout << aFragment << G4endl;
188  G4int EquilibriumExcitonNumber =
189  G4lrint(std::sqrt(aFragment.GetExcitationEnergy()
190  *aFragment.GetA_asInt()*fLevelDensity));
191  //
192  // G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
193  //
194  // J. M. Quesada (Jan. 08) equilibrium hole number could be used as preeq.
195  // evap. delimiter (IAEA report)
196 
197  // Loop for transitions, it is performed while there are
198  // preequilibrium transitions.
199  G4bool ThereIsTransition = false;
200 
201  // G4cout<<"----------------------------------------"<<G4endl;
202  // G4double NP=aFragment.GetNumberOfParticles();
203  // G4double NH=aFragment.GetNumberOfHoles();
204  // G4double NE=aFragment.GetNumberOfExcitons();
205  // G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
206  // G4cout<<"N. excitons="<<NE<<" N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
207  do {
208  ++count;
209  //G4cout<<"transition number .."<<count
210  // <<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
211  G4bool go_ahead = false;
212  // soft cutoff criterium as an "ad-hoc" solution to force
213  // increase in evaporation
214  G4int test = aFragment.GetNumberOfExcitons();
215  if (test <= EquilibriumExcitonNumber) { go_ahead=true; }
216 
217  //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
218  if (useSCO && go_ahead)
219  {
220  G4double x = G4double(test)/G4double(EquilibriumExcitonNumber) - 1;
221  if( G4UniformRand() < 1.0 - G4Exp(-x*x/0.32) ) { go_ahead = false; }
222  }
223 
224  // JMQ: WARNING: CalculateProbability MUST be called prior to Get!!
225  // (O values would be returned otherwise)
226  G4double TotalTransitionProbability =
231  //G4cout<<"#0 P1="<<P1<<" P2="<<P2<<" P3="<<P3<<G4endl;
232 
233  //J.M. Quesada (May 2008) Physical criterium (lamdas) PREVAILS over
234  // approximation (critical exciton number)
235  //V.Ivanchenko (May 2011) added check on number of nucleons
236  // to send a fragment to FermiBreakUp
237  if(!go_ahead || P1 <= P2+P3 ||
238  (aFragment.GetZ_asInt() < maxZ && aFragment.GetA_asInt() < maxA) )
239  {
240  //G4cout<<"#4 EquilibriumEmission"<<G4endl;
241  PerformEquilibriumEmission(aFragment,Result);
242  return Result;
243  }
244  else
245  {
246  //
247  // Check if number of excitons is greater than 0
248  // else perform equilibrium emission
249  if (aFragment.GetNumberOfExcitons() <= 0)
250  {
251  PerformEquilibriumEmission(aFragment,Result);
252  return Result;
253  }
254 
255  G4double TotalEmissionProbability =
256  theEmission->GetTotalProbability(aFragment);
257  //
258  //G4cout<<"#1 TotalEmissionProbability="<<TotalEmissionProbability
259  // <<" Nex= " <<aFragment.GetNumberOfExcitons()<<G4endl;
260  //J.M.Quesada (May 08) this has already been done in order to decide
261  // what to do (preeq-eq)
262  // Sum of all probabilities
263  G4double TotalProbability = TotalEmissionProbability
264  + TotalTransitionProbability;
265 
266  // Select subprocess
267  if (TotalProbability*G4UniformRand() > TotalEmissionProbability)
268  {
269  //G4cout<<"#2 Transition"<<G4endl;
270  // It will be transition to state with a new number of excitons
271  ThereIsTransition = true;
272  // Perform the transition
273  theTransition->PerformTransition(aFragment);
274  }
275  else
276  {
277  //G4cout<<"#3 Emission"<<G4endl;
278  // It will be fragment emission
279  ThereIsTransition = false;
280  Result->push_back(theEmission->PerformEmission(aFragment));
281  }
282  }
283  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
284  } while (ThereIsTransition); // end of do loop
285 
286  // stop if too many iterations
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  PerformEquilibriumEmission(aFragment, Result);
294  return Result;
295  }
296  } // end of for (;;) loop
297  return Result;
298 }
299 
301 // Initialisation
303 
305 {
306  useHETCEmission = true;
308 }
309 
311 {
312  useHETCEmission = false;
314 }
315 
317  useGNASHTransition = true;
318  delete theTransition;
322 }
323 
325  useGNASHTransition = false;
326  delete theTransition;
330 }
331 
333 {
334  OPTxs = opt;
336 }
337 
339 {
340  useSICB = true;
342 }
343 
345 {
346  useNGB = true;
347 }
348 
350 {
351  useSCO = true;
352 }
353 
355 {
356  useCEMtr = true;
357 }
358 
360 // Documentation
362 
363 void G4PreCompoundModel::ModelDescription(std::ostream& outFile) const
364 {
365  outFile << "The GEANT4 precompound model is considered as an extension of the\n"
366  << "hadron kinetic model. It gives a possibility to extend the low energy range\n"
367  << "of the hadron kinetic model for nucleon-nucleus inelastic collision and it \n"
368  << "provides a ”smooth” transition from kinetic stage of reaction described by the\n"
369  << "hadron kinetic model to the equilibrium stage of reaction described by the\n"
370  << "equilibrium deexcitation models.\n"
371  << "The initial information for calculation of pre-compound nuclear stage\n"
372  << "consists of the atomic mass number A, charge Z of residual nucleus, its\n"
373  << "four momentum P0 , excitation energy U and number of excitons n, which equals\n"
374  << "the sum of the number of particles p (from them p_Z are charged) and the number of\n"
375  << "holes h.\n"
376  << "At the preequilibrium stage of reaction, we follow the exciton model approach in ref. [1],\n"
377  << "taking into account the competition among all possible nuclear transitions\n"
378  << "with ∆n = +2, −2, 0 (which are defined by their associated transition probabilities) and\n"
379  << "the emission of neutrons, protons, deutrons, thritium and helium nuclei (also defined by\n"
380  << "their associated emission probabilities according to exciton model)\n"
381  << "\n"
382  << "[1] K.K. Gudima, S.G. Mashnik, V.D. Toneev, Nucl. Phys. A401 329 (1983)\n"
383  << "\n";
384 }
385 
386 void G4PreCompoundModel::DeExciteModelDescription(std::ostream& outFile) const
387 {
388  outFile << "description of precompound model as used with DeExcite()"
389  << "\n";
390 }
391 
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:211
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:352
static const double pi2
Definition: G4SIunits.hh:77
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:333
G4ExcitationHandler * GetExcitationHandler() const
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:251
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
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:410
G4int GetNumberOfExcitons() const
Definition: G4Fragment.hh:317
int G4lrint(double ad)
Definition: templates.hh:163
const G4double x[NPOINTSGL]
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
virtual void ModelDescription(std::ostream &outFile) const
G4int GetZ_asInt() const
Definition: G4Fragment.hh:256
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:268
void SetOPTxs(G4int opt)
G4double GetTotalProbability(const G4Fragment &aFragment)
CLHEP::HepLorentzVector G4LorentzVector