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

Last change on this file since 1197 was 1196, checked in by garnier, 15 years ago

update CVS release candidate geant4.9.3.01
File size: 9.8 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//
27//
28// Hadronic Process: Nuclear De-excitations
29// by V. Lara (Sept 2001)
30//
31// J. M. Quesada : several fixes in total GEM width
32// J. M. Quesada 14/07/2009 bug fixed in total emission width (hbarc)
33// J. M. Quesada (September 2009) several fixes:
34//       -level density parameter of residual calculated at its right excitation energy.
35//       -InitialLeveldensity calculated according to the right conditions of the
36//        initial excited nucleus.
37
38#include "G4GEMProbability.hh"
39#include "G4PairingCorrection.hh"
40
41
42
43G4GEMProbability::G4GEMProbability(const G4GEMProbability &) : G4VEmissionProbability()
44{
45    throw G4HadronicException(__FILE__, __LINE__, "G4GEMProbability::copy_constructor meant to not be accessable");
46}
47
48
49
50
51const G4GEMProbability & G4GEMProbability::
52operator=(const G4GEMProbability &)
53{
54  throw G4HadronicException(__FILE__, __LINE__, "G4GEMProbability::operator= meant to not be accessable");
55  return *this;
56}
57
58
59G4bool G4GEMProbability::operator==(const G4GEMProbability &) const
60{
61  return false;
62}
63
64G4bool G4GEMProbability::operator!=(const G4GEMProbability &) const
65{
66  return true;
67}
68
69G4double G4GEMProbability::EmissionProbability(const G4Fragment & fragment,
70                                               const G4double MaximalKineticEnergy)
71{
72  G4double EmissionProbability = 0.0;
73   
74  if (MaximalKineticEnergy > 0.0 && fragment.GetExcitationEnergy() > 0.0) {
75    G4double CoulombBarrier = GetCoulombBarrier(fragment);
76     
77    EmissionProbability = CalcProbability(fragment,MaximalKineticEnergy,CoulombBarrier);
78    Normalization = EmissionProbability;
79    // Next there is a loop over excited states for this channel summing probabilities
80    if (ExcitationEnergies  &&  ExcitationSpins && ExcitationLifetimes) {
81      G4double SavedSpin = Spin;
82      for (unsigned int i = 0; i < ExcitationEnergies->size(); i++) {
83        Spin = ExcitationSpins->operator[](i);
84        // substract excitation energies
85        G4double Tmax = MaximalKineticEnergy - ExcitationEnergies->operator[](i);
86        if (Tmax > 0.0) {
87          G4double width = CalcProbability(fragment,Tmax,CoulombBarrier);
88          // update probability
89          if (hbar_Planck*std::log(2.0)/width < ExcitationLifetimes->operator[](i)) {
90            EmissionProbability += width;
91          }
92        }
93      }
94      // Restore Spin
95      Spin = SavedSpin;
96    }
97  }
98  return EmissionProbability;
99}
100
101G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment, 
102                                           const G4double MaximalKineticEnergy,
103                                           const G4double V)
104// Calculate integrated probability (width) for evaporation channel
105{
106  G4double ResidualA = static_cast<G4double>(fragment.GetA() - theA);
107  G4double ResidualZ = static_cast<G4double>(fragment.GetZ() - theZ);
108  G4double U = fragment.GetExcitationEnergy();
109 
110  G4double NuclearMass = G4ParticleTable::GetParticleTable()->
111    GetIonTable()->GetNucleusMass(theZ,theA);
112 
113 
114  G4double Alpha = CalcAlphaParam(fragment);
115  G4double Beta = CalcBetaParam(fragment);
116 
117 
118  //                       ***RESIDUAL***
119  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
120 
121  G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection( static_cast<G4int>( ResidualA ) , static_cast<G4int>( ResidualZ ) ); 
122 
123  G4double a = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),
124                                                    static_cast<G4int>(ResidualZ),MaximalKineticEnergy+V-delta0);
125  G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
126  G4double Ex = Ux + delta0;
127  G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
128  //JMQ fixed bug in units
129  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));
130  //                      ***end RESIDUAL ***
131 
132  //                       ***PARENT***
133  //JMQ (September 2009) the following quantities refer to the PARENT:
134     
135  G4double deltaCN=G4PairingCorrection::GetInstance()->
136    GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));                                     
137  G4double aCN = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
138                                                      static_cast<G4int>(fragment.GetZ()),U-deltaCN);
139  G4double UxCN = (2.5 + 150.0/fragment.GetA())*MeV;
140  G4double ExCN = UxCN + deltaCN;
141  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
142  //JMQ fixed bug in units
143  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));
144//                       ***end PARENT***
145
146  G4double Width;
147  G4double InitialLevelDensity;
148  G4double t = MaximalKineticEnergy/T;
149  if ( MaximalKineticEnergy < Ex ) {
150    //JMQ 190709 bug in I1 fixed (T was  missing)
151    Width = (I1(t,t)*T + (Beta+V)*I0(t))/std::exp(E0/T);
152//JMQ 160909 fix:  InitialLevelDensity has been taken away (different conditions for initial CN..)
153  } else {
154    G4double tx = Ex/T;
155    G4double s = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
156    G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
157    Width = I1(t,tx)*T/std::exp(E0/T) + I3(s,sx)*std::exp(s)/(std::sqrt(2.0)*a);
158    // For charged particles (Beta+V) = 0 beacuse Beta = -V
159    if (theZ == 0) {
160      Width += (Beta+V)*(I0(tx)/std::exp(E0/T) + 2.0*std::sqrt(2.0)*I2(s,sx)*std::exp(s));
161    }
162  }
163 
164  //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of hbar_planck must be used
165  //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
166  G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
167 
168  //JMQ 190709 fix on Rb and  geometrical cross sections according to Furihata's paper
169  //                      (JAERI-Data/Code 2001-105, p6)
170  //    G4double RN = 0.0;
171  G4double Rb = 0.0;
172  if (theA > 4) 
173    {
174      //        G4double R1 = std::pow(ResidualA,1.0/3.0);
175      //        G4double R2 = std::pow(G4double(theA),1.0/3.0);
176      G4double Ad = std::pow(ResidualA,1.0/3.0);
177      G4double Aj = std::pow(G4double(theA),1.0/3.0);
178      //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
179      Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
180      Rb *= fermi;
181    }
182  else if (theA>1)
183    {
184      G4double Ad = std::pow(ResidualA,1.0/3.0);
185      G4double Aj = std::pow(G4double(theA),1.0/3.0);
186      Rb=1.5*(Aj+Ad)*fermi;
187    }
188  else
189    {
190      G4double Ad = std::pow(ResidualA,1.0/3.0);
191      Rb = 1.5*Ad*fermi;
192    }
193  //   G4double GeometricalXS = pi*RN*RN*std::pow(ResidualA,2./3.);
194  G4double GeometricalXS = pi*Rb*Rb; 
195  //end of JMQ fix on Rb by 190709
196 
197
198
199//JMQ 160909 fix: initial level density must be calculated according to the
200// conditions at the initial compound nucleus
201// (it has been removed from previous "if" for the residual)
202
203   if ( U < ExCN ) 
204      {
205        InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
206      } 
207    else 
208      {
209        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);
210      }
211 //
212
213
214  //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
215  //    Width *= std::sqrt(pi)*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
216  Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity); 
217   
218
219  return Width;
220}
221
222
223
224
225G4double G4GEMProbability::I0(const G4double t)
226{
227  G4double result = (std::exp(t) - 1.0);
228  return result;
229}
230
231G4double G4GEMProbability::I1(const G4double t, const G4double tx)
232{
233  G4double result = t - tx + 1.0;
234  result *= std::exp(tx);
235  result -= (t + 1.0);
236  return result;
237}
238
239
240G4double G4GEMProbability::I2(const G4double s, const G4double sx)
241{
242  G4double S = 1.0/std::sqrt(s);
243  G4double Sx = 1.0/std::sqrt(sx);
244 
245  G4double p1 = S*S*S*( 1.0 + S*S*( 1.5 + 3.75*S*S) );
246  G4double p2 = Sx*Sx*Sx*( 1.0 + Sx*Sx*( 1.5 + 3.75*Sx*Sx) )*std::exp(sx-s);
247 
248  return p1-p2;
249}
250
251G4double G4GEMProbability::I3(const G4double s, const G4double sx)
252{
253  G4double s2 = s*s;
254  G4double sx2 = sx*sx;
255  G4double S = 1.0/std::sqrt(s);
256  G4double S2 = S*S;
257  G4double Sx = 1.0/std::sqrt(sx);
258  G4double Sx2 = Sx*Sx;
259 
260  G4double p1 = S *(2.0 + S2 *( 4.0 + S2 *( 13.5 + S2 *( 60.0 + S2 * 325.125 ))));
261  G4double p2 = Sx*Sx2 *(
262                         (s2-sx2) + Sx2 *(
263                                          (1.5*s2+0.5*sx2) + Sx2 *(
264                                                                   (3.75*s2+0.25*sx2) + Sx2 *(
265                                                                                              (12.875*s2+0.625*sx2) + Sx2 *(
266                                                                                                                            (59.0625*s2+0.9375*sx2) + Sx2 *(324.8*s2+3.28*sx2))))));
267 
268  p2 *= std::exp(sx-s);
269 
270  return p1-p2; 
271}
Note: See TracBrowser for help on using the repository browser.