ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PreCompoundModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PreCompoundModel.cc
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 // by V. Lara
27 //
28 // Modified:
29 // 01.04.2008 J.M.Quesada Several changes. Soft cut-off switched off.
30 // 01.05.2008 J.M.Quesada Protection against non-physical preeq.
31 // transitional regime has been set
32 // 03.09.2008 J.M.Quesada for external choice of inverse cross section option
33 // 06.09.2008 J.M.Quesada Also external choices have been added for:
34 // - superimposed Coulomb barrier (useSICB=true)
35 // - "never go back" hipothesis (useNGB=true)
36 // - soft cutoff from preeq. to equlibrium (useSCO=true)
37 // - CEM transition probabilities (useCEMtr=true)
38 // 20.08.2010 V.Ivanchenko Cleanup of the code:
39 // - integer Z and A;
40 // - emission and transition classes created at
41 // initialisation
42 // - options are set at initialisation
43 // - do not use copy-constructors for G4Fragment
44 // 03.01.2012 V.Ivanchenko Added pointer to G4ExcitationHandler to the
45 // constructor
46 
47 #include "G4PreCompoundModel.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4PreCompoundEmission.hh"
52 #include "G4GNASHTransitions.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "G4Proton.hh"
55 #include "G4Neutron.hh"
56 
57 #include "G4NucleiProperties.hh"
58 #include "G4NuclearLevelData.hh"
59 #include "G4DeexPrecoParameters.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 
70  : G4VPreCompoundModel(ptr,"PRECO"),theEmission(nullptr),theTransition(nullptr),
71  useSCO(false),isInitialised(false),isActive(true),minZ(3),minA(5)
72 {
73  //G4cout << "### NEW PrecompoundModel " << this << G4endl;
74  if(!ptr) { SetExcitationHandler(new G4ExcitationHandler()); }
75 
79  fLowLimitExc = 0.0;
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 
109 
110  fLowLimitExc = param->GetPrecoLowEnergy();
112 
113  useSCO = param->UseSoftCutoff();
114 
115  minZ = param->GetMinZForPreco();
116  minA = param->GetMinAForPreco();
117 
119  if(param->UseHETC()) { theEmission->SetHETCModel(); }
120  //else { theEmission->SetDefaultModel(); }
122 
123  if(param->UseGNASH()) { theTransition = new G4GNASHTransitions; }
124  else { theTransition = new G4PreCompoundTransitions(); }
125  theTransition->UseNGB(param->NeverGoBack());
126  theTransition->UseCEMtr(param->UseCEM());
127 
128  if(param->PrecoDummy()) { isActive = false; }
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 nullptr;
147  }
148 
149  G4int Zp = 0;
150  G4int Ap = 1;
151  if(primary == proton) { Zp = 1; }
152 
153  G4double timePrimary=thePrimary.GetGlobalTime();
154 
155  G4int A = theNucleus.GetA_asInt();
156  G4int Z = theNucleus.GetZ_asInt();
157 
158  //G4cout << "### G4PreCompoundModel::ApplyYourself: A= " << A << " Z= " << Z
159  // << " Ap= " << Ap << " Zp= " << Zp << G4endl;
160  // 4-Momentum
161  G4LorentzVector p = thePrimary.Get4Momentum();
163  p += G4LorentzVector(0.0,0.0,0.0,mass);
164  //G4cout << "Primary 4-mom " << p << " mass= " << mass << G4endl;
165 
166  // prepare fragment
167  G4Fragment anInitialState(A + Ap, Z + Zp, p);
168  anInitialState.SetNumberOfExcitedParticle(2, 1);
169  anInitialState.SetNumberOfHoles(1,0);
170  anInitialState.SetCreationTime(thePrimary.GetGlobalTime());
171 
172  // call excitation handler
173  G4ReactionProductVector * result = DeExcite(anInitialState);
174 
175  // fill particle change
176  theResult.Clear();
178  for(G4ReactionProductVector::iterator i= result->begin();
179  i != result->end(); ++i)
180  {
181  G4DynamicParticle * aNewDP =
182  new G4DynamicParticle((*i)->GetDefinition(),
183  (*i)->GetTotalEnergy(),
184  (*i)->GetMomentum());
185  G4HadSecondary aNew = G4HadSecondary(aNewDP);
186  G4double time=(*i)->GetFormationTime();
187  if(time < 0.0) { time = 0.0; }
188  aNew.SetTime(timePrimary + time);
189  aNew.SetCreatorModelType((*i)->GetCreatorModel());
190  delete (*i);
191  theResult.AddSecondary(aNew);
192  }
193  delete result;
194 
195  //return the filled particle change
196  return &theResult;
197 }
198 
200 
202 {
203  if(!isInitialised) { InitialiseModel(); }
204 
206  G4double U = aFragment.GetExcitationEnergy();
207  G4int Z = aFragment.GetZ_asInt();
208  G4int A = aFragment.GetA_asInt();
209 
210  //G4cout << "### G4PreCompoundModel::DeExcite" << G4endl;
211  //G4cout << aFragment << G4endl;
212 
213  // Perform Equilibrium Emission
214  if (!isActive || (Z < minZ && A < minA) ||
215  U < fLowLimitExc*A || U > A*fHighLimitExc) {
216  PerformEquilibriumEmission(aFragment, Result);
217  return Result;
218  }
219 
220  // main loop
221  G4int count = 0;
222  const G4double ldfact = 12.0/CLHEP::pi2;
223  const G4int countmax = 1000;
224  for (;;) {
225  //G4cout << "### PreCompound loop over fragment" << G4endl;
226  //G4cout << aFragment << G4endl;
227  U = aFragment.GetExcitationEnergy();
228  Z = aFragment.GetZ_asInt();
229  A = aFragment.GetA_asInt();
230  G4int eqExcitonNumber =
231  G4lrint(std::sqrt(ldfact*U*fNuclData->GetLevelDensity(Z, A, U)));
232  //
233  // G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
234  //
235  // J. M. Quesada (Jan. 08) equilibrium hole number could be used as preeq.
236  // evap. delimiter (IAEA report)
237 
238  // Loop for transitions, it is performed while there are
239  // preequilibrium transitions.
240  G4bool isTransition = false;
241 
242  // G4cout<<"----------------------------------------"<<G4endl;
243  // G4double NP=aFragment.GetNumberOfParticles();
244  // G4double NH=aFragment.GetNumberOfHoles();
245  // G4double NE=aFragment.GetNumberOfExcitons();
246  // G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
247  // G4cout<<"N. excitons="<<NE<<" N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
248  do {
249  ++count;
250  //G4cout<<"transition number .."<<count
251  // <<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
252  // soft cutoff criterium as an "ad-hoc" solution to force
253  // increase in evaporation
254  G4int ne = aFragment.GetNumberOfExcitons();
255  G4bool go_ahead = (ne <= eqExcitonNumber);
256 
257  //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
258  if (useSCO && go_ahead) {
259  G4double x = (G4double)(ne - eqExcitonNumber)/(G4double)eqExcitonNumber;
260  if( G4UniformRand() < 1.0 - G4Exp(-x*x/0.32) ) { go_ahead = false; }
261  }
262 
263  // JMQ: WARNING: CalculateProbability MUST be called prior to Get!!
264  // (O values would be returned otherwise)
265  G4double transProbability = theTransition->CalculateProbability(aFragment);
269  //G4cout<<"#0 P1="<<P1<<" P2="<<P2<<" P3="<<P3<<G4endl;
270 
271  //J.M. Quesada (May 2008) Physical criterium (lamdas) PREVAILS over
272  // approximation (critical exciton number)
273  //V.Ivanchenko (May 2011) added check on number of nucleons
274  // to send a fragment to FermiBreakUp
275  // or check on limits of excitation
276  if(!go_ahead || P1 <= P2+P3 || Z < minZ || A < minA ||
277  U <= fLowLimitExc*A || U > A*fHighLimitExc ||
278  aFragment.GetNumberOfExcitons() <= 0) {
279  //G4cout<<"#4 EquilibriumEmission"<<G4endl;
280  PerformEquilibriumEmission(aFragment,Result);
281  return Result;
282  }
283  G4double emissionProbability = theEmission->GetTotalProbability(aFragment);
284 
285  //G4cout<<"#1 TotalEmissionProbability="<<TotalEmissionProbability
286  // <<" Nex= " <<aFragment.GetNumberOfExcitons()<<G4endl;
287  //J.M.Quesada (May 08) this has already been done in order to decide
288  // what to do (preeq-eq)
289  // Sum of all probabilities
290  G4double TotalProbability = emissionProbability + transProbability;
291 
292  // Select subprocess
293  if (TotalProbability*G4UniformRand() > emissionProbability) {
294  //G4cout<<"#2 Transition"<<G4endl;
295  // It will be transition to state with a new number of excitons
296  isTransition = true;
297  // Perform the transition
298  theTransition->PerformTransition(aFragment);
299  } else {
300  //G4cout<<"#3 Emission"<<G4endl;
301  // It will be fragment emission
302  isTransition = false;
303  Result->push_back(theEmission->PerformEmission(aFragment));
304  }
305  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
306  } while (isTransition); // end of do loop
307 
308  // stop if too many iterations
309  if(count >= countmax) {
311  ed << "G4PreCompoundModel loop over " << countmax << " iterations; "
312  << "current G4Fragment: \n" << aFragment;
313  G4Exception("G4PreCompoundModel::DeExcite()","had0034",JustWarning,
314  ed,"");
315  PerformEquilibriumEmission(aFragment, Result);
316  return Result;
317  }
318  } // end of for (;;) loop
319  return Result;
320 }
321 
323 // Initialisation
325 
327 {
328  PrintWarning("UseHETCEmission");
329 }
330 
332 {
333  PrintWarning("UseDefaultEmission");
334 }
335 
337 {
338  PrintWarning("UseGNASHTransition");
339 }
340 
342 {
343  PrintWarning("UseDefaultTransition");
344 }
345 
347 {
348  PrintWarning("UseOPTxs");
349 }
350 
352 {
353  PrintWarning("UseSICB");
354 }
355 
357 {
358  PrintWarning("UseNGB");
359 }
360 
362 {
363  PrintWarning("UseSCO");
364 }
365 
367 {
368  PrintWarning("UseCEMtr");
369 }
370 
372 {
374  ed << "Obsolete method of the preCompound model is called: "
375  << mname << "() \n Instead a corresponding method of "
376  << "G4DeexPrecoParameters class should be used";
377 
378  G4Exception("G4PreCompoundModel::ReadData()","had0803",JustWarning,ed);
379 }
380 
382 // Documentation
384 
385 void G4PreCompoundModel::ModelDescription(std::ostream& outFile) const
386 {
387  outFile
388  << "The GEANT4 precompound model is considered as an extension of the\n"
389  << "hadron kinetic model. It gives a possibility to extend the low energy range\n"
390  << "of the hadron kinetic model for nucleon-nucleus inelastic collision and it \n"
391  << "provides a ”smooth” transition from kinetic stage of reaction described by the\n"
392  << "hadron kinetic model to the equilibrium stage of reaction described by the\n"
393  << "equilibrium deexcitation models.\n"
394  << "The initial information for calculation of pre-compound nuclear stage\n"
395  << "consists of the atomic mass number A, charge Z of residual nucleus, its\n"
396  << "four momentum P0 , excitation energy U and number of excitons n, which equals\n"
397  << "the sum of the number of particles p (from them p_Z are charged) and the number of\n"
398  << "holes h.\n"
399  << "At the preequilibrium stage of reaction, we follow the exciton model approach in ref. [1],\n"
400  << "taking into account the competition among all possible nuclear transitions\n"
401  << "with ∆n = +2, −2, 0 (which are defined by their associated transition probabilities) and\n"
402  << "the emission of neutrons, protons, deuterons, thritium and helium nuclei (also defined by\n"
403  << "their associated emission probabilities according to exciton model)\n"
404  << "\n"
405  << "[1] K.K. Gudima, S.G. Mashnik, V.D. Toneev, Nucl. Phys. A401 329 (1983)\n"
406  << "\n";
407 }
408 
409 void G4PreCompoundModel::DeExciteModelDescription(std::ostream& outFile) const
410 {
411  outFile << "description of precompound model as used with DeExcite()" << "\n";
412 }
413