Geant4  10.03
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 96603 2016-04-25 13:29:51Z 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"
60 #include "G4NuclearLevelData.hh"
61 #include "G4DeexPrecoParameters.hh"
62 #include "Randomize.hh"
63 #include "G4DynamicParticle.hh"
64 #include "G4ParticleTypes.hh"
65 #include "G4ParticleTable.hh"
66 #include "G4LorentzVector.hh"
67 #include "G4Exp.hh"
68 
70 
72  : G4VPreCompoundModel(ptr,"PRECO"),theEmission(nullptr),theTransition(nullptr),
73  useSCO(false),isInitialised(false),minZ(3),minA(5)
74 {
75  //G4cout << "### NEW PrecompoundModel " << this << G4endl;
76  if(!ptr) { SetExcitationHandler(new G4ExcitationHandler()); }
77 
81 }
82 
84 
86 {
87  delete theEmission;
88  delete theTransition;
89  delete GetExcitationHandler();
90 }
91 
93 
95 {
97 }
98 
100 
102 {
103  if(isInitialised) { return; }
104  isInitialised = true;
105 
106  //G4cout << "G4PreCompoundModel::InitialiseModel() started" << G4endl;
107 
108  G4DeexPrecoParameters* param =
110 
111  // 12/pi2 factor is used in real computation
112  fLevelDensity = param->GetLevelDensity()*12.0/CLHEP::pi2;
113  fLimitEnergy = param->GetPrecoLowEnergy();
114 
115  useSCO = param->UseSoftCutoff();
116 
117  minZ = param->GetMinZForPreco();
118  minA = param->GetMinAForPreco();
119 
121  if(param->UseHETC()) { theEmission->SetHETCModel(); }
122  else { theEmission->SetDefaultModel(); }
124 
125  if(param->UseGNASH()) { theTransition = new G4GNASHTransitions; }
126  else { theTransition = new G4PreCompoundTransitions(); }
127  theTransition->UseNGB(param->NeverGoBack());
128  theTransition->UseCEMtr(param->UseCEM());
129 
131 }
132 
134 
137  G4Nucleus & theNucleus)
138 {
139  const G4ParticleDefinition* primary = thePrimary.GetDefinition();
140  if(primary != neutron && primary != proton) {
142  ed << "G4PreCompoundModel is used for ";
143  if(primary) { ed << primary->GetParticleName(); }
144  G4Exception("G4PreCompoundModel::ApplyYourself()","had0033",FatalException,
145  ed,"");
146  return 0;
147  }
148 
149  G4int Zp = 0;
150  G4int Ap = 1;
151  if(primary == proton) { Zp = 1; }
152 
153  G4int A = theNucleus.GetA_asInt();
154  G4int Z = theNucleus.GetZ_asInt();
155 
156  //G4cout << "### G4PreCompoundModel::ApplyYourself: A= " << A << " Z= " << Z
157  // << " Ap= " << Ap << " Zp= " << Zp << G4endl;
158  // 4-Momentum
159  G4LorentzVector p = thePrimary.Get4Momentum();
161  p += G4LorentzVector(0.0,0.0,0.0,mass);
162  //G4cout << "Primary 4-mom " << p << " mass= " << mass << G4endl;
163 
164  // prepare fragment
165  G4Fragment anInitialState(A + Ap, Z + Zp, p);
166  anInitialState.SetNumberOfExcitedParticle(2, 1);
167  anInitialState.SetNumberOfHoles(1,0);
168  anInitialState.SetCreationTime(thePrimary.GetGlobalTime());
169 
170  // call excitation handler
171  G4ReactionProductVector * result = DeExcite(anInitialState);
172 
173  // fill particle change
174  theResult.Clear();
176  for(G4ReactionProductVector::iterator i= result->begin();
177  i != result->end(); ++i)
178  {
179  G4DynamicParticle * aNew =
180  new G4DynamicParticle((*i)->GetDefinition(),
181  (*i)->GetTotalEnergy(),
182  (*i)->GetMomentum());
183  delete (*i);
184  theResult.AddSecondary(aNew);
185  }
186  delete result;
187 
188  //return the filled particle change
189  return &theResult;
190 }
191 
193 
195 {
196  if(!isInitialised) { InitialiseModel(); }
197 
199  G4double Eex = aFragment.GetExcitationEnergy();
200  G4int Z = aFragment.GetZ_asInt();
201  G4int A = aFragment.GetA_asInt();
202 
203  //G4cout << "### G4PreCompoundModel::DeExcite" << G4endl;
204  //G4cout << aFragment << G4endl;
205 
206  // Perform Equilibrium Emission
207  if ((Z < minZ && A < minA) || Eex < fLimitEnergy*A) {
208  PerformEquilibriumEmission(aFragment, Result);
209  return Result;
210  }
211 
212  // main loop
213  G4int count = 0;
214  static const G4int countmax = 1000;
215  for (;;) {
216  //G4cout << "### PreCompound loop over fragment" << G4endl;
217  //G4cout << aFragment << G4endl;
218  G4int EquilibriumExcitonNumber =
219  G4lrint(std::sqrt(aFragment.GetExcitationEnergy()
220  *aFragment.GetA_asInt()*fLevelDensity));
221  //
222  // G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
223  //
224  // J. M. Quesada (Jan. 08) equilibrium hole number could be used as preeq.
225  // evap. delimiter (IAEA report)
226 
227  // Loop for transitions, it is performed while there are
228  // preequilibrium transitions.
229  G4bool ThereIsTransition = false;
230 
231  // G4cout<<"----------------------------------------"<<G4endl;
232  // G4double NP=aFragment.GetNumberOfParticles();
233  // G4double NH=aFragment.GetNumberOfHoles();
234  // G4double NE=aFragment.GetNumberOfExcitons();
235  // G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
236  // G4cout<<"N. excitons="<<NE<<" N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
237  do {
238  ++count;
239  //G4cout<<"transition number .."<<count
240  // <<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
241  G4bool go_ahead = false;
242  // soft cutoff criterium as an "ad-hoc" solution to force
243  // increase in evaporation
244  G4int test = aFragment.GetNumberOfExcitons();
245  if (test <= EquilibriumExcitonNumber) { go_ahead=true; }
246 
247  //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
248  if (useSCO && go_ahead)
249  {
250  G4double x = G4double(test)/G4double(EquilibriumExcitonNumber) - 1;
251  if( G4UniformRand() < 1.0 - G4Exp(-x*x/0.32) ) { go_ahead = false; }
252  }
253 
254  // JMQ: WARNING: CalculateProbability MUST be called prior to Get!!
255  // (O values would be returned otherwise)
256  G4double TotalTransitionProbability =
261  //G4cout<<"#0 P1="<<P1<<" P2="<<P2<<" P3="<<P3<<G4endl;
262 
263  //J.M. Quesada (May 2008) Physical criterium (lamdas) PREVAILS over
264  // approximation (critical exciton number)
265  //V.Ivanchenko (May 2011) added check on number of nucleons
266  // to send a fragment to FermiBreakUp
267  if(!go_ahead || P1 <= P2+P3 ||
268  aFragment.GetZ_asInt() < minZ || aFragment.GetA_asInt() < minA ||
269  (aFragment.GetExcitationEnergy() <= fLimitEnergy*aFragment.GetA_asInt()))
270  {
271  //G4cout<<"#4 EquilibriumEmission"<<G4endl;
272  PerformEquilibriumEmission(aFragment,Result);
273  return Result;
274  }
275  else
276  {
277  //
278  // Check if number of excitons is greater than 0
279  // else perform equilibrium emission
280  if (aFragment.GetNumberOfExcitons() <= 0)
281  {
282  PerformEquilibriumEmission(aFragment,Result);
283  return Result;
284  }
285 
286  G4double TotalEmissionProbability =
287  theEmission->GetTotalProbability(aFragment);
288  //
289  //G4cout<<"#1 TotalEmissionProbability="<<TotalEmissionProbability
290  // <<" Nex= " <<aFragment.GetNumberOfExcitons()<<G4endl;
291  //J.M.Quesada (May 08) this has already been done in order to decide
292  // what to do (preeq-eq)
293  // Sum of all probabilities
294  G4double TotalProbability = TotalEmissionProbability
295  + TotalTransitionProbability;
296 
297  // Select subprocess
298  if (TotalProbability*G4UniformRand() > TotalEmissionProbability)
299  {
300  //G4cout<<"#2 Transition"<<G4endl;
301  // It will be transition to state with a new number of excitons
302  ThereIsTransition = true;
303  // Perform the transition
304  theTransition->PerformTransition(aFragment);
305  }
306  else
307  {
308  //G4cout<<"#3 Emission"<<G4endl;
309  // It will be fragment emission
310  ThereIsTransition = false;
311  Result->push_back(theEmission->PerformEmission(aFragment));
312  }
313  }
314  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
315  } while (ThereIsTransition); // end of do loop
316 
317  // stop if too many iterations
318  if(count >= countmax) {
320  ed << "G4PreCompoundModel loop over " << countmax << " iterations; "
321  << "current G4Fragment: \n" << aFragment;
322  G4Exception("G4PreCompoundModel::DeExcite()","had0034",JustWarning,
323  ed,"");
324  PerformEquilibriumEmission(aFragment, Result);
325  return Result;
326  }
327  } // end of for (;;) loop
328  return Result;
329 }
330 
332 // Initialisation
334 
336 {
337  PrintWarning("UseHETCEmission");
338 }
339 
341 {
342  PrintWarning("UseDefaultEmission");
343 }
344 
346 {
347  PrintWarning("UseGNASHTransition");
348 }
349 
351 {
352  PrintWarning("UseDefaultTransition");
353 }
354 
356 {
357  PrintWarning("UseOPTxs");
358 }
359 
361 {
362  PrintWarning("UseSICB");
363 }
364 
366 {
367  PrintWarning("UseNGB");
368 }
369 
371 {
372  PrintWarning("UseSCO");
373 }
374 
376 {
377  PrintWarning("UseCEMtr");
378 }
379 
381 {
383  ed << "Obsolete method of the preCompound model is called: "
384  << mname << "() \n Instead a corresponding method of "
385  << "G4DeexPrecoParameters class should be used";
386 
387  G4Exception("G4PreCompoundModel::ReadData()","had0803",JustWarning,ed);
388 }
389 
391 // Documentation
393 
394 void G4PreCompoundModel::ModelDescription(std::ostream& outFile) const
395 {
396  outFile
397  << "The GEANT4 precompound model is considered as an extension of the\n"
398  << "hadron kinetic model. It gives a possibility to extend the low energy range\n"
399  << "of the hadron kinetic model for nucleon-nucleus inelastic collision and it \n"
400  << "provides a ”smooth” transition from kinetic stage of reaction described by the\n"
401  << "hadron kinetic model to the equilibrium stage of reaction described by the\n"
402  << "equilibrium deexcitation models.\n"
403  << "The initial information for calculation of pre-compound nuclear stage\n"
404  << "consists of the atomic mass number A, charge Z of residual nucleus, its\n"
405  << "four momentum P0 , excitation energy U and number of excitons n, which equals\n"
406  << "the sum of the number of particles p (from them p_Z are charged) and the number of\n"
407  << "holes h.\n"
408  << "At the preequilibrium stage of reaction, we follow the exciton model approach in ref. [1],\n"
409  << "taking into account the competition among all possible nuclear transitions\n"
410  << "with ∆n = +2, −2, 0 (which are defined by their associated transition probabilities) and\n"
411  << "the emission of neutrons, protons, deutrons, thritium and helium nuclei (also defined by\n"
412  << "their associated emission probabilities according to exciton model)\n"
413  << "\n"
414  << "[1] K.K. Gudima, S.G. Mashnik, V.D. Toneev, Nucl. Phys. A401 329 (1983)\n"
415  << "\n";
416 }
417 
418 void G4PreCompoundModel::DeExciteModelDescription(std::ostream& outFile) const
419 {
420  outFile << "description of precompound model as used with DeExcite()" << "\n";
421 }
422 
virtual void PerformTransition(G4Fragment &aFragment)=0
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
void PerformEquilibriumEmission(const G4Fragment &aFragment, G4ReactionProductVector *theResult) const
static G4double GetNuclearMass(const G4double A, const G4double Z)
static const G4double * P1[nN]
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4PreCompoundEmission * theEmission
G4double GetPrecoLowEnergy() const
void PrintWarning(const G4String &mname)
const G4ParticleDefinition * neutron
virtual void ModelDescription(std::ostream &outFile) const final
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:357
virtual void DeExciteModelDescription(std::ostream &outFile) const final
G4double GetLevelDensity() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4HadFinalState theResult
virtual void BuildPhysicsTable(const G4ParticleDefinition &) final
void SetStatusChange(G4HadFinalStateStatus aS)
void SetExcitationHandler(G4ExcitationHandler *ptr)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void SetNumberOfExcitedParticle(G4int valueTot, G4int valueP)
Definition: G4Fragment.hh:338
G4ExcitationHandler * GetExcitationHandler() const
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:256
const G4ParticleDefinition * GetDefinition() const
G4VPreCompoundTransitions * theTransition
bool G4bool
Definition: G4Types.hh:79
G4double GetGlobalTime() const
G4DeexPrecoParameters * GetParameters()
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const G4LorentzVector & Get4Momentum() const
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:425
G4int GetNumberOfExcitons() const
Definition: G4Fragment.hh:322
int G4lrint(double ad)
Definition: templates.hh:163
virtual void InitialiseModel() final
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4int GetZ_asInt() const
Definition: G4Fragment.hh:261
static const G4double * P2[nN]
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
const G4ParticleDefinition * proton
double G4double
Definition: G4Types.hh:76
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment) final
static G4NuclearLevelData * GetInstance()
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:273
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus) final
void SetOPTxs(G4int opt)
G4double GetTotalProbability(const G4Fragment &aFragment)
static constexpr double pi2
Definition: G4SIunits.hh:78
CLHEP::HepLorentzVector G4LorentzVector
G4PreCompoundModel(G4ExcitationHandler *ptr=nullptr)