ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EvaporationProbability.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EvaporationProbability.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 // J.M. Quesada (August2008). Based on:
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara (Oct 1998)
30 //
31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse
32 // cross section option
33 // JMQ (06 September 2008) Also external choices have been added for
34 // superimposed Coulomb barrier (if useSICB is set true, by default is false)
35 //
36 // JMQ (14 february 2009) bug fixed in emission width: hbarc instead of
37 // hbar_Planck in the denominator
38 //
39 // V.Ivanchenko general clean-up since 2010
40 //
41 #include "G4EvaporationProbability.hh"
42 #include "G4NuclearLevelData.hh"
43 #include "G4VCoulombBarrier.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4PairingCorrection.hh"
47 #include "G4NucleiProperties.hh"
48 #include "G4KalbachCrossSection.hh"
49 #include "G4ChatterjeeCrossSection.hh"
50 #include "Randomize.hh"
51 #include "G4Exp.hh"
52 #include "G4Log.hh"
53 #include "G4Pow.hh"
54 
55 using namespace std;
56 
57 static const G4double explim = 160.;
58 
59 G4EvaporationProbability::G4EvaporationProbability(G4int anA, G4int aZ,
60  G4double aGamma)
61  : G4VEmissionProbability(aZ, anA), fGamma(aGamma)
62 {
63  resA13 = muu = freeU = a0 = delta1 = 0.0;
64  pcoeff = fGamma*pEvapMass*CLHEP::millibarn
65  /((CLHEP::pi*CLHEP::hbarc)*(CLHEP::pi*CLHEP::hbarc));
66 
67  if(0 == theZ) { index = 0; }
68  else if(1 == theZ) { index = theA; }
69  else { index = theA + 1; }
70  if(0 == aZ) {
71  ResetIntegrator(30, 0.25*CLHEP::MeV, 0.02);
72  } else {
73  ResetIntegrator(30, 0.5*CLHEP::MeV, 0.03);
74  }
75  // G4cout << "G4EvaporationProbability: Z= " << theZ << " A= " << theA
76  // << " M(GeV)= " << pEvapMass/GeV << G4endl;
77 }
78 
79 G4EvaporationProbability::~G4EvaporationProbability()
80 {}
81 
82 G4double G4EvaporationProbability::CalcAlphaParam(const G4Fragment&)
83 {
84  return 1.0;
85 }
86 
87 G4double G4EvaporationProbability::CalcBetaParam(const G4Fragment&)
88 {
89  return 1.0;
90 }
91 
92 G4double G4EvaporationProbability::TotalProbability(
93  const G4Fragment& fragment, G4double minEnergy, G4double maxEnergy,
94  G4double CB, G4double exEnergy)
95 {
96  G4int fragA = fragment.GetA_asInt();
97  G4int fragZ = fragment.GetZ_asInt();
98  G4double U = fragment.GetExcitationEnergy();
99  a0 = pNuclearLevelData->GetLevelDensity(fragZ,fragA,U);
100  freeU = exEnergy;
101  resA13 = pG4pow->Z13(resA);
102  delta1 = pNuclearLevelData->GetPairingCorrection(resZ,resA);
103  /*
104  G4cout << "G4EvaporationProbability: Z= " << theZ << " A= " << theA
105  << " resZ= " << resZ << " resA= " << resA
106  << " fragZ= " << fragZ << " fragA= " << fragA
107  << "\n freeU= " << freeU
108  << " a0= " << a0 << " OPT= " << OPTxs << " emin= "
109  << minEnergy << " emax= " << maxEnergy
110  << " CB= " << CB << G4endl;
111  */
112  if (OPTxs==0 || (OPTxs==4 && freeU < 10.)) {
113 
114  G4double SystemEntropy = 2.0*std::sqrt(a0*freeU);
115 
116  const G4double RN2 = 2.25*CLHEP::fermi*CLHEP::fermi
117  /(CLHEP::twopi*CLHEP::hbar_Planck*hbar_Planck);
118 
119  G4double Alpha = CalcAlphaParam(fragment);
120  G4double Beta = CalcBetaParam(fragment);
121 
122  // to be checked where to use a0, where - a1
123  G4double a1 = pNuclearLevelData->GetLevelDensity(resZ,resA,freeU);
124  G4double GlobalFactor = fGamma*Alpha*pEvapMass*RN2*resA13*resA13/(a1*a1);
125 
126  G4double maxea = maxEnergy*a1;
127  G4double Term1 = Beta*a1 - 1.5 + maxea;
128  G4double Term2 = (2.0*Beta*a1-3.0)*std::sqrt(maxea) + 2*maxea;
129 
130  G4double ExpTerm1 = (SystemEntropy <= explim) ? G4Exp(-SystemEntropy) : 0.0;
131 
132  G4double ExpTerm2 = 2.*std::sqrt(maxea) - SystemEntropy;
133  ExpTerm2 = std::min(ExpTerm2, explim);
134  ExpTerm2 = G4Exp(ExpTerm2);
135 
136  pProbability = GlobalFactor*(Term1*ExpTerm1 + Term2*ExpTerm2);
137 
138  } else {
139 
140  // compute power once
141  if(0 < index) {
142  muu = G4KalbachCrossSection::ComputePowerParameter(resA, index);
143  }
144  // if Coulomb barrier cutoff is superimposed for all cross sections
145  // then the limit is the Coulomb Barrier
146  pProbability = IntegrateProbability(minEnergy, maxEnergy, CB);
147  }
148  return pProbability;
149 }
150 
151 G4double G4EvaporationProbability::ComputeProbability(G4double K, G4double CB)
152 {
153  //G4cout << "### G4EvaporationProbability::ProbabilityDistributionFunction"
154  // << G4endl;
155 
156  G4double E0 = freeU;
157  // abnormal case - should never happens
158  if(pMass < pEvapMass + pResMass) { return 0.0; }
159 
160  G4double m02 = pMass*pMass;
161  G4double m12 = pEvapMass*pEvapMass;
162  G4double mres = sqrt(m02 + m12 - 2.*pMass*(pEvapMass + K));
163 
164  G4double excRes = mres - pResMass;
165  G4double E1 = excRes - delta1;
166  if(E1 <= 0.0) { return 0.0; }
167  G4double a1 = pNuclearLevelData->GetLevelDensity(resZ,resA,excRes);
168  G4double xs = CrossSection(K, CB);
169  G4double prob = pcoeff*G4Exp(2.0*(std::sqrt(a1*E1) - std::sqrt(a0*E0)))*K*xs;
170  /*
171  G4cout << "PDF: Z= " << theZ << " A= " << theA
172  << " K= " << K << " E0= " << E0 << " E1= " << E1 << G4endl;
173  G4cout << " prob= " << prob << " pcoeff= " << pcoeff
174  << " xs= " << xs << G4endl;
175  */
176  return prob;
177 }
178 
179 G4double
180 G4EvaporationProbability::CrossSection(G4double K, G4double CB)
181 {
182  G4double res;
183  if(OPTxs <= 2) {
184  res = G4ChatterjeeCrossSection::ComputeCrossSection(K, CB, resA13, muu,
185  index, theZ, resA);
186  } else {
187  res = G4KalbachCrossSection::ComputeCrossSection(K, CB, resA13, muu,
188  index, theZ, theA, resA);
189  }
190  //G4cout << "XS: K= "<<K<<" res= "<<res<<" cb= "<<CB<<" muu= "
191  // <<muu<<" index= " << index<< G4endl;
192  return res;
193 }
194 
195 G4double
196 G4EvaporationProbability::SampleKineticEnergy(G4double minKinEnergy,
197  G4double maxKinEnergy,
198  G4double)
199 {
200  /*
201  G4cout << "### Sample probability Emin= " << minKinEnergy
202  << " Emax= " << maxKinEnergy
203  << " Z= " << theZ << " A= " << theA << G4endl;
204  */
205  G4double T = 0.0;
206  CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine();
207  if (OPTxs==0 || (OPTxs==4 && freeU < 10.)) {
208  // JMQ:
209  // It uses Dostrovsky's approximation for the inverse reaction cross
210  // in the probability for fragment emission
211  // MaximalKineticEnergy energy in the original version (V.Lara) was
212  // calculated at the Coulomb barrier.
213 
214  G4double Rb = 4.0*a0*maxKinEnergy;
215  G4double RbSqrt = std::sqrt(Rb);
216  G4double PEX1 = (RbSqrt < explim) ? G4Exp(-RbSqrt) : 0.0;
217  G4double Rk = 0.0;
218  G4double FRk = 0.0;
219  G4int nn = 0;
220  const G4int nmax = 100;
221  const G4double ssqr3 = 1.5*std::sqrt(3.0);
222  do {
223  G4double RandNumber = rndm->flat();
224  Rk = 1.0 + (1./RbSqrt)*G4Log(RandNumber + (1.0-RandNumber)*PEX1);
225  G4double Q1 = 1.0;
226  G4double Q2 = 1.0;
227  if (theZ == 0) { // for emitted neutron
228  G4double Beta = (2.12/(resA13*resA13) - 0.05)*MeV/(0.76 + 2.2/resA13);
229  Q1 = 1.0 + Beta/maxKinEnergy;
230  Q2 = Q1*std::sqrt(Q1);
231  }
232 
233  FRk = ssqr3 * Rk * (Q1 - Rk*Rk)/Q2;
234  if(nn > nmax) { break; }
235  ++nn;
236  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
237  } while (FRk < rndm->flat());
238 
239  T = std::max(maxKinEnergy * (1.0-Rk*Rk), 0.0) + minKinEnergy;
240 
241  } else {
242  T = SampleEnergy();
243  }
244  //G4cout<<"-- new Z= "<<theZ<<" A= "<< theA << " ekin= " << T << G4endl;
245  return T;
246 }