source: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc @ 1315

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 10.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//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
34#include "G4ProtonEvaporationProbability.hh"
35
36G4ProtonEvaporationProbability::G4ProtonEvaporationProbability() :
37    G4EvaporationProbability(1,1,2,&theCoulombBarrier) // A,Z,Gamma,&theCoulombBarrier
38{
39   
40}
41
42G4ProtonEvaporationProbability::G4ProtonEvaporationProbability(const G4ProtonEvaporationProbability &) : G4EvaporationProbability()
43{
44    throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationProbability::copy_constructor meant to not be accessable");
45}
46
47const G4ProtonEvaporationProbability & G4ProtonEvaporationProbability::
48operator=(const G4ProtonEvaporationProbability &)
49{
50    throw G4HadronicException(__FILE__, __LINE__, "G4ProtonEvaporationProbability::operator= meant to not be accessable");
51    return *this;
52}
53
54
55G4bool G4ProtonEvaporationProbability::operator==(const G4ProtonEvaporationProbability &) const
56{
57    return false;
58}
59
60G4bool G4ProtonEvaporationProbability::operator!=(const G4ProtonEvaporationProbability &) const
61{
62    return true;
63}
64
65  G4double G4ProtonEvaporationProbability::CalcAlphaParam(const G4Fragment & fragment) 
66  { return 1.0 + CCoeficient(static_cast<G4double>(fragment.GetZ()-GetZ()));}
67       
68  G4double G4ProtonEvaporationProbability::CalcBetaParam(const G4Fragment & ) 
69  { return 0.0; }
70
71  G4double G4ProtonEvaporationProbability::CCoeficient(const G4double aZ) 
72{
73    // Data comes from
74    // Dostrovsky, Fraenkel and Friedlander
75    // Physical Review, vol 116, num. 3 1959
76    //
77    // const G4int size = 5;
78    // G4double Zlist[5] = { 10.0, 20.0, 30.0, 50.0, 70.0};
79    // G4double Cp[5] = { 0.50, 0.28, 0.20, 0.15, 0.10};
80    G4double C = 0.0;
81       
82    if (aZ >= 70) {
83        C = 0.10;
84    } else {
85        C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
86    }
87       
88    return C;
89       
90}
91
92///////////////////////////////////////////////////////////////////////////////////
93//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections for protons
94//OPT=0 Dostrovski's parameterization
95//OPT=1 Chatterjee's parameterization
96//OPT=2,4 Wellisch's parameterization
97//OPT=3 Kalbach's parameterization
98//
99G4double G4ProtonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
100{
101//  G4cout<<" In G4ProtonEVaporationProbability OPTxs="<<OPTxs<<G4endl;
102//  G4cout<<" In G4ProtonEVaporationProbability useSICB="<<useSICB<<G4endl;
103
104  theA=GetA();
105  theZ=GetZ();
106  ResidualA=fragment.GetA()-theA;
107  ResidualZ=fragment.GetZ()-theZ; 
108 
109  ResidualAthrd=std::pow(ResidualA,0.33333);
110  FragmentA=fragment.GetA();
111  FragmentAthrd=std::pow(FragmentA,0.33333);
112  U=fragment.GetExcitationEnergy();
113
114  if (OPTxs==0) {std::ostringstream errOs;
115    errOs << "We should'n be here (OPT =0) at evaporation cross section calculation (protons)!!"  <<G4endl;
116    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
117    return 0.;}
118  else if( OPTxs==1 ) return GetOpt1( K);
119  else if( OPTxs==2 ||OPTxs==4) return GetOpt2( K);
120  else if (OPTxs==3 )  return GetOpt3( K);
121  else{
122    std::ostringstream errOs;
123    errOs << "BAD PROTON CROSS SECTION OPTION AT EVAPORATION!!"  <<G4endl;
124    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
125    return 0.;
126  }
127}
128//********************* OPT=1 : Chatterjee's cross section ************************
129//(fitting to cross section from Bechetti & Greenles OM potential)
130
131G4double G4ProtonEvaporationProbability::GetOpt1(const  G4double K)
132{
133  G4double Kc=K; 
134
135// JMQ  xsec is set constat above limit of validity
136  if (K>50)  Kc=50;
137
138  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2,xs;
139  G4double p, p0, p1, p2,Ec,delta,q,r,ji;
140 
141  p0 = 15.72;
142  p1 = 9.65;
143  p2 = -449.0;
144  landa0 = 0.00437;
145  landa1 = -16.58;
146  mu0 = 244.7;
147  mu1 = 0.503;
148  nu0 = 273.1;
149  nu1 = -182.4;
150  nu2 = -1.872; 
151  delta=0.; 
152
153  Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
154  p = p0 + p1/Ec + p2/(Ec*Ec);
155  landa = landa0*ResidualA + landa1;
156  mu = mu0*std::pow(ResidualA,mu1);
157  nu = std::pow(ResidualA,mu1)*(nu0 + nu1*Ec + nu2*(Ec*Ec));
158  q = landa - nu/(Ec*Ec) - 2*p*Ec;
159  r = mu + 2*nu/Ec + p*(Ec*Ec);
160
161  ji=std::max(Kc,Ec);
162  if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
163  else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
164   if (xs <0.0) {xs=0.0;}
165
166   return xs; 
167
168}
169
170
171
172//************* OPT=2 : Wellisch's proton reaction cross section ************************
173
174G4double G4ProtonEvaporationProbability::GetOpt2(const  G4double K)
175{
176
177  G4double rnpro,rnneu,eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
178 
179//This is redundant when the Coulomb  barrier is overimposed to all cross sections
180//It should be kept when Coulomb barrier only imposed at OPTxs=2, this is why ..
181
182         if(!useSICB && K <= theCoulombBarrier.GetCoulombBarrier(G4lrint(ResidualA),G4lrint(ResidualZ),U)) return xine_th=0.0;
183
184  eekin=K;
185  rnpro=ResidualZ;
186  rnneu=ResidualA-ResidualZ;
187  ekin=eekin/1000;
188
189  r0=1.36*1.e-15;
190  fac=pi*r0*r0;
191  b0=2.247-0.915*(1.-1./ResidualAthrd);
192  fac1=b0*(1.-1./ResidualAthrd);
193  fac2=1.;
194  if(rnneu > 1.5) fac2=std::log(rnneu);
195  xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1);
196  xine_th=(1.-0.15*std::exp(-ekin))*xine_th/(1.00-0.0007*ResidualA);   
197  ff1=0.70-0.0020*ResidualA ;
198  ff2=1.00+1/ResidualA;
199  ff3=0.8+18/ResidualA-0.002*ResidualA;
200  fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2))));
201  xine_th=xine_th*(1.+ff3*fac);
202  ff1=1.-1/ResidualA-0.001*ResidualA;
203  ff2=1.17-2.7/ResidualA-0.0014*ResidualA;
204  fac=-8.*ff1*(std::log10(ekin)+2.0*ff2);
205  fac=1./(1.+std::exp(fac));
206  xine_th=xine_th*fac;               
207  if (xine_th < 0.0){
208    std::ostringstream errOs;
209    G4cout<<"WARNING:  negative Wellisch cross section "<<G4endl; 
210    errOs << "RESIDUAL: A=" << ResidualA << " Z=" << ResidualZ <<G4endl;
211    errOs <<"  xsec("<<ekin<<" MeV) ="<<xine_th <<G4endl;
212    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
213  }
214
215  return xine_th;
216           
217}
218
219
220// *********** OPT=3 : Kalbach's cross sections (from PRECO code)*************
221G4double G4ProtonEvaporationProbability::GetOpt3(const  G4double K)
222{
223//     ** p from  becchetti and greenlees (but modified with sub-barrier
224//     ** correction function and xp2 changed from -449)
225
226  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
227  G4double p, p0, p1, p2;
228  p0 = 15.72;
229  p1 = 9.65;
230  p2 = -300.;
231  landa0 = 0.00437;
232  landa1 = -16.58;
233  mu0 = 244.7;
234  mu1 = 0.503;
235  nu0 = 273.1;
236  nu1 = -182.4;
237  nu2 = -1.872;
238
239// parameters for  proton cross section refinement
240  G4double afit,bfit,a2,b2;
241  afit=-0.0785656;
242  bfit=5.10789;
243  a2= -0.00089076;
244  b2= 0.0231597; 
245
246  G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig; 
247  G4double b,ecut,cut,ecut2,geom,elab;
248
249
250  G4double      flow = 1.e-18;
251  G4double       spill= 1.e+18; 
252
253
254  if (ResidualA <= 60.)  signor = 0.92;
255  else if (ResidualA < 100.) signor = 0.8 + ResidualA*0.002;
256
257
258  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
259  ecsq = ec * ec;
260  p = p0 + p1/ec + p2/ecsq;
261  landa = landa0*ResidualA + landa1;
262  a = std::pow(ResidualA,mu1);
263  mu = mu0 * a;
264  nu = a* (nu0+nu1*ec+nu2*ecsq);
265 
266  c =std::min(3.15,ec*0.5);
267  w = 0.7 * c / 3.15; 
268
269  xnulam = nu / landa;
270  if (xnulam > spill) xnulam=0.;
271  if (xnulam >= flow) etest =std::sqrt(xnulam) + 7.;
272
273  a = -2.*p*ec + landa - nu/ecsq;
274  b = p*ecsq + mu + 2.*nu/ec;
275  ecut = 0.;
276  cut = a*a - 4.*p*b;
277  if (cut > 0.) ecut = std::sqrt(cut);
278  ecut = (ecut-a) / (p+p);
279  ecut2 = ecut;
280//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
281//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
282//  if (cut < 0.) ecut2 = ecut - 2.;
283 if (cut < 0.) ecut2 = ecut;
284  elab = K * FragmentA / ResidualA;
285  sig = 0.;
286  if (elab <= ec) { //start for E<Ec
287    if (elab > ecut2)  sig = (p*elab*elab+a*elab+b) * signor;
288    signor2 = (ec-elab-c) / w;
289    signor2 = 1. + std::exp(signor2);
290    sig = sig / signor2;     
291                       }       //end for E<=Ec
292  else            {           //start for  E>Ec
293    sig = (landa*elab+mu+nu/elab) * signor;
294    geom = 0.;
295
296    if (xnulam < flow || elab < etest)
297     {
298        if (sig <0.0) {sig=0.0;}
299        return sig;
300      }
301    geom = std::sqrt(theA*K);
302    geom = 1.23*ResidualAthrd + ra + 4.573/geom;
303    geom = 31.416 * geom * geom;
304    sig = std::max(geom,sig);
305
306  }   //end for E>Ec
307 return sig;}
308
309
310
311//   ************************** end of cross sections *******************************
312
Note: See TracBrowser for help on using the repository browser.