ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GEMProbability.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GEMProbability.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 //
27 //---------------------------------------------------------------------
28 //
29 // Geant4 class G4GEMProbability
30 //
31 //
32 // Hadronic Process: Nuclear De-excitations
33 // by V. Lara (Sept 2001)
34 //
35 //
36 // Hadronic Process: Nuclear De-excitations
37 // by V. Lara (Sept 2001)
38 //
39 // J. M. Quesada : several fixes in total GEM width
40 // J. M. Quesada 14/07/2009 bug fixed in total emission width (hbarc)
41 // J. M. Quesada (September 2009) several fixes:
42 // -level density parameter of residual calculated at its right excitation energy.
43 // -InitialLeveldensity calculated according to the right conditions of the
44 // initial excited nucleus.
45 // J. M. Quesada 19/04/2010 fix in emission probability calculation.
46 // V.Ivanchenko 20/04/2010 added usage of G4Pow and use more safe computation
47 // V.Ivanchenko 18/05/2010 trying to speedup the most slow method
48 // by usage of G4Pow, integer Z and A; moved constructor,
49 // destructor and virtual functions to source
50 
51 #include "G4GEMProbability.hh"
52 #include "G4PairingCorrection.hh"
53 #include "G4NucleiProperties.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4Log.hh"
57 
58 G4GEMProbability:: G4GEMProbability(G4int anA, G4int aZ, G4double aSpin)
59  : G4VEmissionProbability(aZ, anA), Spin(aSpin),
60  theCoulombBarrierPtr(nullptr)
61 {
62  theEvapLDPptr = new G4EvaporationLevelDensityParameter;
63  fG4pow = G4Pow::GetInstance();
64  fPlanck= CLHEP::hbar_Planck*fG4pow->logZ(2);
65  fNucData = G4NuclearLevelData::GetInstance();
66 }
67 
68 G4GEMProbability::~G4GEMProbability()
69 {
70  delete theEvapLDPptr;
71 }
72 
73 G4double G4GEMProbability::EmissionProbability(const G4Fragment & fragment,
74  G4double MaximalKineticEnergy)
75 {
76  G4double probability = 0.0;
77 
78  if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
79  G4double CoulombBarrier = GetCoulombBarrier(fragment);
80 
81  probability =
82  CalcProbability(fragment,MaximalKineticEnergy,CoulombBarrier);
83 
84  // Next there is a loop over excited states for this channel
85  // summing probabilities
86  size_t nn = ExcitEnergies.size();
87  if (0 < nn) {
88  G4double SavedSpin = Spin;
89  for (size_t i = 0; i <nn; ++i) {
90  Spin = ExcitSpins[i];
91  // substract excitation energies
92  G4double Tmax = MaximalKineticEnergy - ExcitEnergies[i];
93  if (Tmax > 0.0) {
94  G4double width = CalcProbability(fragment,Tmax,CoulombBarrier);
95  //JMQ April 2010 added condition to prevent reported crash
96  // update probability
97  if (width > 0. && fPlanck < width*ExcitLifetimes[i]) {
98  probability += width;
99  }
100  }
101  }
102  // Restore Spin
103  Spin = SavedSpin;
104  }
105  }
106  // Normalization = probability;
107  return probability;
108 }
109 
110 G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
111  G4double MaximalKineticEnergy,
112  G4double V)
113 
114 // Calculate integrated probability (width) for evaporation channel
115 {
116  G4int A = fragment.GetA_asInt();
117  G4int Z = fragment.GetZ_asInt();
118 
119  G4int ResidualA = A - theA;
120  G4int ResidualZ = Z - theZ;
121  G4double U = fragment.GetExcitationEnergy();
122 
123  G4double NuclearMass = fragment.ComputeGroundStateMass(theZ, theA);
124 
125  G4double Alpha = CalcAlphaParam(fragment);
126  G4double Beta = CalcBetaParam(fragment);
127 
128  // ***RESIDUAL***
129  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
130 
131  G4double delta0 = fNucData->GetPairingCorrection(ResidualZ, ResidualA);
132 
133  G4double a = theEvapLDPptr->
134  LevelDensityParameter(ResidualA,ResidualZ,MaximalKineticEnergy+V-delta0);
135  G4double Ux = (2.5 + 150.0/G4double(ResidualA))*MeV;
136  G4double Ex = Ux + delta0;
137  G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
138  //JMQ fixed bug in units
139  G4double E0 = Ex - T*(G4Log(T/MeV) - G4Log(a*MeV)/4.0
140  - 1.25*G4Log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
141  // ***end RESIDUAL ***
142  // ***PARENT***
143  //JMQ (September 2009) the following quantities refer to the PARENT:
144 
145  G4double deltaCN = fNucData->GetPairingCorrection(Z, A);
146  G4double aCN = theEvapLDPptr->LevelDensityParameter(A, Z, U-deltaCN);
147  G4double UxCN = (2.5 + 150.0/G4double(A))*MeV;
148  G4double ExCN = UxCN + deltaCN;
149  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
150  // ***end PARENT***
151 
152  G4double Width;
153  G4double InitialLevelDensity;
154  G4double t = MaximalKineticEnergy/T;
155  if ( MaximalKineticEnergy < Ex ) {
156  //JMQ 190709 bug in I1 fixed (T was missing)
157  Width = (I1(t,t)*T + (Beta+V)*I0(t))/G4Exp(E0/T);
158  //JMQ 160909 fix: InitialLevelDensity has been taken away
159  //(different conditions for initial CN..)
160  } else {
161 
162  //VI minor speedup
163  G4double expE0T = G4Exp(E0/T);
164  static const G4double sqrt2 = std::sqrt(2.0);
165 
166  G4double tx = Ex/T;
167  G4double s0 = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
168  G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
169  // VI: protection against FPE exception
170  if(s0 > 350.) { s0 = 350.; }
171  Width = I1(t,tx)*T/expE0T + I3(s0,sx)*G4Exp(s0)/(sqrt2*a);
172 
173  // VI this cannot happen!
174  // For charged particles (Beta+V) = 0 beacuse Beta = -V
175  //if (theZ == 0) {
176  // Width += (Beta+V)*(I0(tx)/expE0T + 2.0*sqrt2*I2(s0,sx)*G4Exp(s0));
177  //}
178  }
179 
180  //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of
181  // hbar_planck must be used
182  // G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
183  G4double gg = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
184 
185  //JMQ 190709 fix on Rb and geometrical cross sections according to
186  // Furihata's paper (JAERI-Data/Code 2001-105, p6)
187  G4double Rb = 0.0;
188  if (theA > 4)
189  {
190  G4double Ad = fG4pow->Z13(ResidualA);
191  G4double Aj = fG4pow->Z13(theA);
192  Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
193  Rb *= fermi;
194  }
195  else if (theA>1)
196  {
197  G4double Ad = fG4pow->Z13(ResidualA);
198  G4double Aj = fG4pow->Z13(theA);
199  Rb=1.5*(Aj+Ad)*fermi;
200  }
201  else
202  {
203  G4double Ad = fG4pow->Z13(ResidualA);
204  Rb = 1.5*Ad*fermi;
205  }
206  G4double GeometricalXS = pi*Rb*Rb;
207  //end of JMQ fix on Rb by 190709
208 
209  //JMQ 160909 fix: initial level density must be calculated according to the
210  // conditions at the initial compound nucleus
211  // (it has been removed from previous "if" for the residual)
212 
213  if ( U < ExCN )
214  {
215  //JMQ fixed bug in units
216  //VI moved the computation here
217  G4double E0CN = ExCN - TCN*(G4Log(TCN/MeV) - 0.25*G4Log(aCN*MeV)
218  - 1.25*G4Log(UxCN/MeV)
219  + 2.0*std::sqrt(aCN*UxCN));
220 
221  InitialLevelDensity = (pi/12.0)*G4Exp((U-E0CN)/TCN)/TCN;
222  }
223  else
224  {
225  //VI speedup
226  G4double x = U-deltaCN;
227  G4double x1 = std::sqrt(aCN*x);
228 
229  InitialLevelDensity = (pi/12.0)*G4Exp(2*x1)/(x*std::sqrt(x1));
230  }
231 
232  //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according
233  // to Furihata's report:
234  Width *= pi*gg*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
235 
236  return Width;
237 }
238 
239 G4double G4GEMProbability::I3(G4double s0, G4double sx)
240 {
241  G4double s2 = s0*s0;
242  G4double sx2 = sx*sx;
243  G4double S = 1.0/std::sqrt(s0);
244  G4double S2 = S*S;
245  G4double Sx = 1.0/std::sqrt(sx);
246  G4double Sx2 = Sx*Sx;
247 
248  G4double p1 = S *(2.0 + S2 *( 4.0 + S2 *( 13.5 + S2 *( 60.0 + S2 * 325.125 ))));
249  G4double p2 = Sx*Sx2 *(
250  (s2-sx2) + Sx2 *(
251  (1.5*s2+0.5*sx2) + Sx2 *(
252  (3.75*s2+0.25*sx2) + Sx2 *(
253  (12.875*s2+0.625*sx2) + Sx2 *(
254  (59.0625*s2+0.9375*sx2) + Sx2 *(324.8*s2+3.28*sx2))))));
255 
256  p2 *= G4Exp(sx-s0);
257 
258  return p1-p2;
259 }
260 
261 void G4GEMProbability::Dump() const
262 {
263  G4double mass = G4NucleiProperties::GetNuclearMass(theA, theZ);
264  G4double efermi = 0.0;
265  if(theA > 1) {
266  efermi = G4NucleiProperties::GetNuclearMass(theA-1, theZ)
267  + neutron_mass_c2 - mass;
268  }
269  G4int nlev = ExcitEnergies.size();
270  G4cout << "GEM: List of Excited States for Isotope Z= "
271  << theZ << " A= " << theA << " Nlevels= " << nlev
272  << " Efermi(MeV)= " << efermi
273  << G4endl;
274  for(G4int i=0; i< nlev; ++i) {
275  G4cout << "Z= " << theZ << " A= " << theA
276  << " Mass(GeV)= " << mass/GeV
277  << " Eexc(MeV)= " << ExcitEnergies[i]
278  << " Time(ns)= " << ExcitLifetimes[i]/ns
279  << G4endl;
280  }
281  G4cout << G4endl;
282 }