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