source: trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMProbability.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 15.1 KB
Line 
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: G4GEMProbability.cc,v 1.13 2010/05/19 10:21:16 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-ref-09 $
28//
29//---------------------------------------------------------------------
30//
31// Geant4 class G4GEMProbability
32//
33//
34// Hadronic Process: Nuclear De-excitations
35// by V. Lara (Sept 2001)
36//
37//
38// Hadronic Process: Nuclear De-excitations
39// by V. Lara (Sept 2001)
40//
41// J. M. Quesada : several fixes in total GEM width
42// J. M. Quesada 14/07/2009 bug fixed in total emission width (hbarc)
43// J. M. Quesada (September 2009) several fixes:
44//       -level density parameter of residual calculated at its right excitation energy.
45//       -InitialLeveldensity calculated according to the right conditions of the
46//        initial excited nucleus.
47// J. M. Quesada 19/04/2010 fix in  emission probability calculation.
48// V.Ivanchenko  20/04/2010 added usage of G4Pow and use more safe computation
49// V.Ivanchenko 18/05/2010 trying to speedup the most slow method
50//            by usage of G4Pow, integer Z and A; moved constructor,
51//            destructor and virtual functions to source
52
53#include "G4GEMProbability.hh"
54#include "G4PairingCorrection.hh"
55#include "G4Pow.hh"
56#include "G4IonTable.hh"
57
58G4GEMProbability:: G4GEMProbability() : 
59  theA(0), theZ(0), Spin(0.0), theCoulombBarrierPtr(0),
60  ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
61{
62  theEvapLDPptr = new G4EvaporationLevelDensityParameter;
63  fG4pow = G4Pow::GetInstance(); 
64  fPairCorr = G4PairingCorrection::GetInstance();
65}
66
67G4GEMProbability:: G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin) : 
68  theA(anA), theZ(aZ), Spin(aSpin), theCoulombBarrierPtr(0),
69  ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
70{
71  theEvapLDPptr = new G4EvaporationLevelDensityParameter;
72  fG4pow = G4Pow::GetInstance(); 
73  fPairCorr = G4PairingCorrection::GetInstance();
74}
75   
76G4GEMProbability::~G4GEMProbability()
77{
78  delete theEvapLDPptr;
79}
80
81G4double G4GEMProbability::CalcAlphaParam(const G4Fragment & ) const 
82{
83  return 1.0;
84}
85
86G4double G4GEMProbability::CalcBetaParam(const G4Fragment & ) const 
87{
88  return 1.0;
89}
90
91G4double G4GEMProbability::CCoeficient(const G4double ) const 
92{
93  return 0.0;
94}
95
96G4double G4GEMProbability::EmissionProbability(const G4Fragment & fragment,
97                                               const G4double MaximalKineticEnergy)
98{
99  G4double EmissionProbability = 0.0;
100   
101  if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
102    G4double CoulombBarrier = GetCoulombBarrier(fragment);
103     
104    EmissionProbability = CalcProbability(fragment,MaximalKineticEnergy,CoulombBarrier);
105    Normalization = EmissionProbability;
106    // Next there is a loop over excited states for this channel summing probabilities
107    if (ExcitationEnergies  &&  ExcitationSpins && ExcitationLifetimes) {
108      G4double SavedSpin = Spin;
109      for (size_t i = 0; i < ExcitationEnergies->size(); ++i) {
110        Spin = ExcitationSpins->operator[](i);
111        // substract excitation energies
112        G4double Tmax = MaximalKineticEnergy - ExcitationEnergies->operator[](i);
113        if (Tmax > 0.0) {
114          G4double width = CalcProbability(fragment,Tmax,CoulombBarrier);
115          //JMQ April 2010 added condition to prevent reported crash
116          // update probability
117          if (width > 0. && 
118              hbar_Planck*fG4pow->logZ(2) < width*ExcitationLifetimes->operator[](i)) {
119            EmissionProbability += width;
120          }
121        }
122      }
123      // Restore Spin
124      Spin = SavedSpin;
125    }
126  }
127  return EmissionProbability;
128}
129
130
131G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment, 
132                                           const G4double MaximalKineticEnergy,
133                                           const G4double V)
134
135// Calculate integrated probability (width) for evaporation channel
136{
137  G4int A = fragment.GetA_asInt();
138  G4int Z = fragment.GetZ_asInt();
139
140  G4int ResidualA = A - theA;
141  G4int ResidualZ = Z - theZ;
142  G4double U = fragment.GetExcitationEnergy();
143 
144  G4double NuclearMass = fragment.ComputeGroundStateMass(theZ, theA);
145 
146  G4double Alpha = CalcAlphaParam(fragment);
147  G4double Beta = CalcBetaParam(fragment);
148 
149 
150  //                       ***RESIDUAL***
151  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
152 
153  G4double delta0 = fPairCorr->GetPairingCorrection(ResidualA, ResidualZ); 
154 
155  G4double a = theEvapLDPptr->LevelDensityParameter(ResidualA,
156                                                    ResidualZ,MaximalKineticEnergy+V-delta0);
157  G4double Ux = (2.5 + 150.0/G4double(ResidualA))*MeV;
158  G4double Ex = Ux + delta0;
159  G4double T  = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
160  //JMQ fixed bug in units
161  G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
162  //                      ***end RESIDUAL ***
163 
164  //                       ***PARENT***
165  //JMQ (September 2009) the following quantities refer to the PARENT:
166     
167  G4double deltaCN=fPairCorr->GetPairingCorrection(A, Z);                                     
168  G4double aCN = theEvapLDPptr->LevelDensityParameter(A, Z, U-deltaCN);
169  G4double UxCN = (2.5 + 150.0/G4double(A))*MeV;
170  G4double ExCN = UxCN + deltaCN;
171  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
172  //JMQ fixed bug in units
173  G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
174//                       ***end PARENT***
175
176  G4double Width;
177  G4double InitialLevelDensity;
178  G4double t = MaximalKineticEnergy/T;
179  if ( MaximalKineticEnergy < Ex ) {
180    //JMQ 190709 bug in I1 fixed (T was  missing)
181    Width = (I1(t,t)*T + (Beta+V)*I0(t))/std::exp(E0/T);
182    //JMQ 160909 fix:  InitialLevelDensity has been taken away
183    //(different conditions for initial CN..)
184  } else {
185
186    //VI minor speedup
187    G4double expE0T = std::exp(E0/T);
188    const G4double sqrt2 = std::sqrt(2.0);
189
190    G4double tx = Ex/T;
191    G4double s = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
192    G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
193    Width = I1(t,tx)*T/expE0T + I3(s,sx)*std::exp(s)/(sqrt2*a);
194    // For charged particles (Beta+V) = 0 beacuse Beta = -V
195    if (theZ == 0) {
196      Width += (Beta+V)*(I0(tx)/expE0T + 2.0*sqrt2*I2(s,sx)*std::exp(s));
197    }
198  }
199 
200  //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of hbar_planck must be used
201  //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
202  G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
203 
204  //JMQ 190709 fix on Rb and  geometrical cross sections according to Furihata's paper
205  //                      (JAERI-Data/Code 2001-105, p6)
206  //    G4double RN = 0.0;
207  G4double Rb = 0.0;
208  if (theA > 4) 
209    {
210      //        G4double R1 = std::pow(ResidualA,1.0/3.0);
211      //        G4double R2 = std::pow(G4double(theA),1.0/3.0);
212      G4double Ad = fG4pow->Z13(ResidualA);
213      G4double Aj = fG4pow->Z13(theA);
214      //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
215      Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
216      Rb *= fermi;
217    }
218  else if (theA>1)
219    {
220      G4double Ad = fG4pow->Z13(ResidualA);
221      G4double Aj = fG4pow->Z13(theA);
222      Rb=1.5*(Aj+Ad)*fermi;
223    }
224  else
225    {
226      G4double Ad = fG4pow->Z13(ResidualA);
227      Rb = 1.5*Ad*fermi;
228    }
229  //   G4double GeometricalXS = pi*RN*RN*std::pow(ResidualA,2./3.);
230  G4double GeometricalXS = pi*Rb*Rb; 
231  //end of JMQ fix on Rb by 190709
232 
233
234
235//JMQ 160909 fix: initial level density must be calculated according to the
236// conditions at the initial compound nucleus
237// (it has been removed from previous "if" for the residual)
238
239   if ( U < ExCN ) 
240      {
241        InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
242      } 
243    else 
244      {
245        //        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
246        //VI speedup
247        G4double x  = U-deltaCN;
248        G4double x2 = x*x;
249        G4double x3 = x2*x;
250        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*x))/std::sqrt(std::sqrt(aCN*x2*x3));
251      }
252 //
253
254
255  //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
256  //    Width *= std::sqrt(pi)*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
257  Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity); 
258   
259
260  return Width;
261}
262
263/*
264
265G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
266                                           const G4double MaximalKineticEnergy,
267                                           const G4double V)
268// Calculate integrated probability (width) for evaporation channel
269{
270  G4double ResidualA = static_cast<G4double>(fragment.GetA() - theA);
271  G4double ResidualZ = static_cast<G4double>(fragment.GetZ() - theZ);
272  G4double U = fragment.GetExcitationEnergy();
273 
274  G4double NuclearMass = G4ParticleTable::GetParticleTable()->
275    GetIonTable()->GetNucleusMass(theZ,theA);
276 
277 
278  G4double Alpha = CalcAlphaParam(fragment);
279  G4double Beta = CalcBetaParam(fragment);
280 
281 
282  //                       ***RESIDUAL***
283  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
284 
285  G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection( static_cast<G4int>( ResidualA ) , static_cast<G4int>( ResidualZ ) ); 
286 
287  G4double a = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),
288                                                    static_cast<G4int>(ResidualZ),MaximalKineticEnergy+V-delta0);
289  G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
290  G4double Ex = Ux + delta0;
291  G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
292  //JMQ fixed bug in units
293  G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
294  //                      ***end RESIDUAL ***
295 
296  //                       ***PARENT***
297  //JMQ (September 2009) the following quantities refer to the PARENT:
298     
299  G4double deltaCN=G4PairingCorrection::GetInstance()->
300    GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));                                     
301  G4double aCN = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
302                                                      static_cast<G4int>(fragment.GetZ()),U-deltaCN);
303  G4double UxCN = (2.5 + 150.0/fragment.GetA())*MeV;
304  G4double ExCN = UxCN + deltaCN;
305  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
306  //JMQ fixed bug in units
307  G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
308//                       ***end PARENT***
309
310  G4double Width;
311  G4double InitialLevelDensity;
312  G4double t = MaximalKineticEnergy/T;
313  if ( MaximalKineticEnergy < Ex ) {
314    //JMQ 190709 bug in I1 fixed (T was  missing)
315    Width = (I1(t,t)*T + (Beta+V)*I0(t))/std::exp(E0/T);
316//JMQ 160909 fix:  InitialLevelDensity has been taken away (different conditions for initial CN..)
317  } else {
318    G4double tx = Ex/T;
319    G4double s = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
320    G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
321    Width = I1(t,tx)*T/std::exp(E0/T) + I3(s,sx)*std::exp(s)/(std::sqrt(2.0)*a);
322    // For charged particles (Beta+V) = 0 beacuse Beta = -V
323    if (theZ == 0) {
324      Width += (Beta+V)*(I0(tx)/std::exp(E0/T) + 2.0*std::sqrt(2.0)*I2(s,sx)*std::exp(s));
325    }
326  }
327 
328  //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of hbar_planck must be used
329  //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
330  G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
331 
332  //JMQ 190709 fix on Rb and  geometrical cross sections according to Furihata's paper
333  //                      (JAERI-Data/Code 2001-105, p6)
334  //    G4double RN = 0.0;
335  G4double Rb = 0.0;
336  if (theA > 4)
337    {
338      //        G4double R1 = std::pow(ResidualA,1.0/3.0);
339      //        G4double R2 = std::pow(G4double(theA),1.0/3.0);
340      G4double Ad = std::pow(ResidualA,1.0/3.0);
341      G4double Aj = std::pow(G4double(theA),1.0/3.0);
342      //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
343      Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
344      Rb *= fermi;
345    }
346  else if (theA>1)
347    {
348      G4double Ad = std::pow(ResidualA,1.0/3.0);
349      G4double Aj = std::pow(G4double(theA),1.0/3.0);
350      Rb=1.5*(Aj+Ad)*fermi;
351    }
352  else
353    {
354      G4double Ad = std::pow(ResidualA,1.0/3.0);
355      Rb = 1.5*Ad*fermi;
356    }
357  //   G4double GeometricalXS = pi*RN*RN*std::pow(ResidualA,2./3.);
358  G4double GeometricalXS = pi*Rb*Rb;
359  //end of JMQ fix on Rb by 190709
360 
361
362
363//JMQ 160909 fix: initial level density must be calculated according to the
364// conditions at the initial compound nucleus
365// (it has been removed from previous "if" for the residual)
366
367   if ( U < ExCN )
368      {
369        InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
370      }
371    else
372      {
373        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
374      }
375 //
376
377
378  //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
379  //    Width *= std::sqrt(pi)*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
380  Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
381   
382
383  return Width;
384}
385*/
386
387G4double G4GEMProbability::I3(const G4double s, const G4double sx)
388{
389  G4double s2 = s*s;
390  G4double sx2 = sx*sx;
391  G4double S = 1.0/std::sqrt(s);
392  G4double S2 = S*S;
393  G4double Sx = 1.0/std::sqrt(sx);
394  G4double Sx2 = Sx*Sx;
395 
396  G4double p1 = S *(2.0 + S2 *( 4.0 + S2 *( 13.5 + S2 *( 60.0 + S2 * 325.125 ))));
397  G4double p2 = Sx*Sx2 *(
398                         (s2-sx2) + Sx2 *(
399                                          (1.5*s2+0.5*sx2) + Sx2 *(
400                                                                   (3.75*s2+0.25*sx2) + Sx2 *(
401                                                                                              (12.875*s2+0.625*sx2) + Sx2 *(
402                                                                                                                            (59.0625*s2+0.9375*sx2) + Sx2 *(324.8*s2+3.28*sx2))))));
403 
404  p2 *= std::exp(sx-s);
405 
406  return p1-p2; 
407}
Note: See TracBrowser for help on using the repository browser.