Geant4_10
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 74903 2013-10-23 16:47:40Z 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 initialisation
43 // - options are set at initialisation
44 // - do not use copy-constructors for G4Fragment
45 // 03.01.2012 V.Ivanchenko Added pointer to G4ExcitationHandler to the
46 // constructor
47 
48 #include "G4PreCompoundModel.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4PreCompoundEmission.hh"
53 #include "G4GNASHTransitions.hh"
54 #include "G4ParticleDefinition.hh"
55 #include "G4Proton.hh"
56 #include "G4Neutron.hh"
57 
58 #include "G4NucleiProperties.hh"
60 #include "Randomize.hh"
61 #include "G4DynamicParticle.hh"
62 #include "G4ParticleTypes.hh"
63 #include "G4ParticleTable.hh"
64 #include "G4LorentzVector.hh"
65 #include "G4Exp.hh"
66 
68  : G4VPreCompoundModel(ptr,"PRECO"), useHETCEmission(false),
69  useGNASHTransition(false), OPTxs(3), useSICB(false),
70  useNGB(false), useSCO(false), useCEMtr(true), maxZ(3), maxA(5)
71  //maxZ(9), maxA(17)
72 {
73  if(!ptr) { SetExcitationHandler(new G4ExcitationHandler()); }
75  fLevelDensity = param.GetLevelDensity();
76 
77  theEmission = new G4PreCompoundEmission();
78  if(useHETCEmission) { theEmission->SetHETCModel(); }
79  else { theEmission->SetDefaultModel(); }
80  theEmission->SetOPTxs(OPTxs);
81  theEmission->UseSICB(useSICB);
82 
83  if(useGNASHTransition) { theTransition = new G4GNASHTransitions; }
84  else { theTransition = new G4PreCompoundTransitions(); }
85  theTransition->UseNGB(useNGB);
86  theTransition->UseCEMtr(useCEMtr);
87 
88  proton = G4Proton::Proton();
89  neutron = G4Neutron::Neutron();
90 }
91 
93 {
94  delete theEmission;
95  delete theTransition;
96  delete GetExcitationHandler();
97 }
98 
100 {
101  outFile << "The GEANT4 precompound model is considered as an extension of the\n"
102  << "hadron kinetic model. It gives a possibility to extend the low energy range\n"
103  << "of the hadron kinetic model for nucleon-nucleus inelastic collision and it \n"
104  << "provides a ”smooth” transition from kinetic stage of reaction described by the\n"
105  << "hadron kinetic model to the equilibrium stage of reaction described by the\n"
106  << "equilibrium deexcitation models.\n"
107  << "The initial information for calculation of pre-compound nuclear stage\n"
108  << "consists of the atomic mass number A, charge Z of residual nucleus, its\n"
109  << "four momentum P0 , excitation energy U and number of excitons n, which equals\n"
110  << "the sum of the number of particles p (from them p_Z are charged) and the number of\n"
111  << "holes h.\n"
112  << "At the preequilibrium stage of reaction, we follow the exciton model approach in ref. [1],\n"
113  << "taking into account the competition among all possible nuclear transitions\n"
114  << "with ∆n = +2, −2, 0 (which are defined by their associated transition probabilities) and\n"
115  << "the emission of neutrons, protons, deutrons, thritium and helium nuclei (also defined by\n"
116  << "their associated emission probabilities according to exciton model)\n"
117  << "\n"
118  << "[1] K.K. Gudima, S.G. Mashnik, V.D. Toneev, Nucl. Phys. A401 329 (1983)\n"
119  << std::endl;
120 }
122 
124  G4Nucleus & theNucleus)
125 {
126  const G4ParticleDefinition* primary = thePrimary.GetDefinition();
127  if(primary != neutron && primary != proton) {
128  std::ostringstream errOs;
129  errOs << "BAD primary type in G4PreCompoundModel: "
130  << primary->GetParticleName() <<G4endl;
131  throw G4HadronicException(__FILE__, __LINE__, errOs.str());
132  }
133 
134  G4int Zp = 0;
135  G4int Ap = 1;
136  if(primary == proton) { Zp = 1; }
137 
138  G4int A = theNucleus.GetA_asInt();
139  G4int Z = theNucleus.GetZ_asInt();
140 
141  //G4cout << "### G4PreCompoundModel::ApplyYourself: A= " << A << " Z= " << Z
142  // << " Ap= " << Ap << " Zp= " << Zp << G4endl;
143  // 4-Momentum
144  G4LorentzVector p = thePrimary.Get4Momentum();
146  p += G4LorentzVector(0.0,0.0,0.0,mass);
147  //G4cout << "Primary 4-mom " << p << " mass= " << mass << G4endl;
148 
149  // prepare fragment
150  G4Fragment anInitialState(A + Ap, Z + Zp, p);
151 
152  // projectile and target nucleons
153  // Add nucleon on which interaction happens
154  //++Ap;
155  //if(A*G4UniformRand() <= G4double(Z)) { Zp += 1; }
156  anInitialState.SetNumberOfExcitedParticle(2, 1);
157  anInitialState.SetNumberOfHoles(1,0);
158  // anInitialState.SetNumberOfExcitedParticle(Ap, Zp);
159  // anInitialState.SetNumberOfHoles(Ap,Zp);
160 
161  anInitialState.SetCreationTime(thePrimary.GetGlobalTime());
162 
163  // call excitation handler
164  G4ReactionProductVector * result = DeExcite(anInitialState);
165 
166  // fill particle change
167  theResult.Clear();
168  theResult.SetStatusChange(stopAndKill);
169  for(G4ReactionProductVector::iterator i= result->begin(); i != result->end(); ++i)
170  {
171  G4DynamicParticle * aNew =
172  new G4DynamicParticle((*i)->GetDefinition(),
173  (*i)->GetTotalEnergy(),
174  (*i)->GetMomentum());
175  delete (*i);
176  theResult.AddSecondary(aNew);
177  }
178  delete result;
179 
180  //return the filled particle change
181  return &theResult;
182 }
183 
185 
187 {
189  G4double Eex = aFragment.GetExcitationEnergy();
190  G4int Z = aFragment.GetZ_asInt();
191  G4int A = aFragment.GetA_asInt();
192 
193  //G4cout << "### G4PreCompoundModel::DeExcite" << G4endl;
194  //G4cout << aFragment << G4endl;
195 
196  // Perform Equilibrium Emission
197  if ((Z < maxZ && A < maxA) || Eex < MeV /*|| Eex > 3.*MeV*A*/) {
198  PerformEquilibriumEmission(aFragment, Result);
199  return Result;
200  }
201 
202  // main loop
203  for (;;) {
204 
205  //fragment++;
206  //G4cout<<"-------------------------------------------------------------------"<<G4endl;
207  //G4cout<<"Fragment number .. "<<fragment<<G4endl;
208 
209  // Initialize fragment according with the nucleus parameters
210  //G4cout << "### Loop over fragment" << G4endl;
211  //G4cout << aFragment << G4endl;
212 
213  theEmission->Initialize(aFragment);
214 
215  G4double gg = (6.0/pi2)*aFragment.GetA_asInt()*fLevelDensity;
216 
217  G4int EquilibriumExcitonNumber =
218  G4lrint(std::sqrt(2*gg*aFragment.GetExcitationEnergy()));
219  //
220  // G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
221  //
222  // J. M. Quesada (Jan. 08) equilibrium hole number could be used as preeq.
223  // evap. delimiter (IAEA report)
224 
225  // Loop for transitions, it is performed while there are
226  // preequilibrium transitions.
227  G4bool ThereIsTransition = false;
228 
229  // G4cout<<"----------------------------------------"<<G4endl;
230  // G4double NP=aFragment.GetNumberOfParticles();
231  // G4double NH=aFragment.GetNumberOfHoles();
232  // G4double NE=aFragment.GetNumberOfExcitons();
233  // G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
234  // G4cout<<"N. excitons="<<NE<<" N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
235  //G4int transition=0;
236  do {
237  //transition++;
238  //G4cout<<"transition number .."<<transition<<G4endl;
239  //G4cout<<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
240  G4bool go_ahead = false;
241  // soft cutoff criterium as an "ad-hoc" solution to force
242  // increase in evaporation
243  G4int test = aFragment.GetNumberOfExcitons();
244  if (test <= EquilibriumExcitonNumber) { go_ahead=true; }
245 
246  //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
247  if (useSCO && go_ahead)
248  {
249  G4double x = G4double(test)/G4double(EquilibriumExcitonNumber) - 1;
250  if( G4UniformRand() < 1.0 - G4Exp(-x*x/0.32) ) { go_ahead = false; }
251  }
252 
253  // JMQ: WARNING: CalculateProbability MUST be called prior to Get methods !!
254  // (O values would be returned otherwise)
255  G4double TotalTransitionProbability =
256  theTransition->CalculateProbability(aFragment);
257  G4double P1 = theTransition->GetTransitionProb1();
258  G4double P2 = theTransition->GetTransitionProb2();
259  G4double P3 = theTransition->GetTransitionProb3();
260  //G4cout<<"#0 P1="<<P1<<" P2="<<P2<<" P3="<<P3<<G4endl;
261 
262  //J.M. Quesada (May 2008) Physical criterium (lamdas) PREVAILS over
263  // approximation (critical exciton number)
264  //V.Ivanchenko (May 2011) added check on number of nucleons
265  // to send a fragment to FermiBreakUp
266  if(!go_ahead || P1 <= P2+P3 ||
267  (aFragment.GetZ_asInt() < maxZ && aFragment.GetA_asInt() < maxA) )
268  {
269  //G4cout<<"#4 EquilibriumEmission"<<G4endl;
270  PerformEquilibriumEmission(aFragment,Result);
271  return Result;
272  }
273  else
274  {
275  G4double TotalEmissionProbability =
276  theEmission->GetTotalProbability(aFragment);
277  //
278  //G4cout<<"#1 TotalEmissionProbability="<<TotalEmissionProbability
279  // <<" Nex= " <<aFragment.GetNumberOfExcitons()<<G4endl;
280  //
281  // Check if number of excitons is greater than 0
282  // else perform equilibrium emission
283  if (aFragment.GetNumberOfExcitons() <= 0)
284  {
285  PerformEquilibriumEmission(aFragment,Result);
286  return Result;
287  }
288 
289  //J.M.Quesada (May 08) this has already been done in order to decide
290  // what to do (preeq-eq)
291  // Sum of all probabilities
292  G4double TotalProbability = TotalEmissionProbability
293  + TotalTransitionProbability;
294 
295  // Select subprocess
296  if (TotalProbability*G4UniformRand() > TotalEmissionProbability)
297  {
298  //G4cout<<"#2 Transition"<<G4endl;
299  // It will be transition to state with a new number of excitons
300  ThereIsTransition = true;
301  // Perform the transition
302  theTransition->PerformTransition(aFragment);
303  }
304  else
305  {
306  //G4cout<<"#3 Emission"<<G4endl;
307  // It will be fragment emission
308  ThereIsTransition = false;
309  Result->push_back(theEmission->PerformEmission(aFragment));
310  }
311  }
312  } while (ThereIsTransition); // end of do loop
313  } // end of for (;;) loop
314  return Result;
315 }
316 
318 // Initialisation
320 
322 {
323  useHETCEmission = true;
324  theEmission->SetHETCModel();
325 }
326 
328 {
329  useHETCEmission = false;
330  theEmission->SetDefaultModel();
331 }
332 
334  useGNASHTransition = true;
335  delete theTransition;
336  theTransition = new G4GNASHTransitions;
337  theTransition->UseNGB(useNGB);
338  theTransition->UseCEMtr(useCEMtr);
339 }
340 
342  useGNASHTransition = false;
343  delete theTransition;
344  theTransition = new G4PreCompoundTransitions();
345  theTransition->UseNGB(useNGB);
346  theTransition->UseCEMtr(useCEMtr);
347 }
348 
350 {
351  OPTxs = opt;
352  theEmission->SetOPTxs(OPTxs);
353 }
354 
356 {
357  useSICB = true;
358  theEmission->UseSICB(useSICB);
359 }
360 
362 {
363  useNGB = true;
364 }
365 
367 {
368  useSCO = true;
369 }
370 
372 {
373  useCEMtr = true;
374 }
375 
virtual void PerformTransition(G4Fragment &aFragment)=0
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static G4double GetNuclearMass(const G4double A, const G4double Z)
std::ofstream outFile
Definition: GammaRayTel.cc:68
G4PreCompoundModel(G4ExcitationHandler *ptr=0)
const char * p
Definition: xmltok.h:285
G4double G4NeutronHPJENDLHEData::G4double result
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:335
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)
const G4String & GetParticleName() const
void SetStatusChange(G4HadFinalStateStatus aS)
void SetExcitationHandler(G4ExcitationHandler *ptr)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void SetNumberOfExcitedParticle(G4int valueTot, G4int valueP)
Definition: G4Fragment.hh:316
G4ExcitationHandler * GetExcitationHandler() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4int GetA_asInt() const
Definition: G4Fragment.hh:218
const G4ParticleDefinition * GetDefinition() const
Float_t Z
Definition: plot.C:39
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)
virtual G4double CalculateProbability(const G4Fragment &aFragment)=0
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4ReactionProduct * PerformEmission(G4Fragment &aFragment)
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:383
G4int GetNumberOfExcitons() const
Definition: G4Fragment.hh:300
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:223
#define G4endl
Definition: G4ios.hh:61
def test
Definition: mcscore.py:117
double G4double
Definition: G4Types.hh:76
void AddSecondary(G4DynamicParticle *aP)
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:235
void SetOPTxs(G4int opt)
void Initialize(const G4Fragment &aFragment)
G4double GetTotalProbability(const G4Fragment &aFragment)
CLHEP::HepLorentzVector G4LorentzVector