Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ProjectileFragmentCrossSection Class Reference

#include <G4ProjectileFragmentCrossSection.hh>

Public Member Functions

 G4ProjectileFragmentCrossSection ()
 
G4double doit (G4double Ap, G4double Zp, G4double At, G4double Zt, G4double A, G4double Z)
 
void testMe ()
 

Detailed Description

Definition at line 39 of file G4ProjectileFragmentCrossSection.hh.

Constructor & Destructor Documentation

G4ProjectileFragmentCrossSection::G4ProjectileFragmentCrossSection ( )
inline

Definition at line 42 of file G4ProjectileFragmentCrossSection.hh.

43  {
44  p_S[1] = -2.38; // scale factor for xsect in barn
45  p_S[2] = 0.27;
46 
47  p_P[1] = -2.5840E+00; // slope of mass yield curve
48  p_P[2] = -7.5700E-03;
49 
50  p_Delta[1] = -1.0870E+00; // centroid rel. to beta-stability
51  p_Delta[2] = +3.0470E-02;
52  p_Delta[3] = +2.1353E-04;
53  p_Delta[4] = +7.1350E+01;
54 
55  p_R[1] = +0.885E+00; // width parameter R
56  p_R[2] = -9.8160E-03;
57 
58  p_Un[1] = 1.65; // slope par. n-rich ride of Z distr.
59 
60  p_Up[1] = 1.7880; // slope par. p-rich ride of Z distr.
61  p_Up[2] = +4.7210E-03;
62  p_Up[3] = -1.3030E-05;
63 
64  p_mn[1] = 0.400; // memory effect n-rich projectiles
65  p_mn[2] = 0.600;
66 
67  p_mp[1] = -10.25; // memory effect p-rich projectiles
68  p_mp[2] = +10.1;
69 
70  corr_d[1] = -25.0; // correction close to proj.: centroid dzp
71  corr_d[2] = 0.800;
72  corr_r[1] = +20.0; // correction close to proj.: width R
73  corr_r[2] = 0.820;
74  corr_y[1] = 200.0; // correction close to proj.: Yield_a
75  corr_y[2] = 0.90;
76  }

Member Function Documentation

G4double G4ProjectileFragmentCrossSection::doit ( G4double  Ap,
G4double  Zp,
G4double  At,
G4double  Zt,
G4double  A,
G4double  Z 
)
inline

Definition at line 78 of file G4ProjectileFragmentCrossSection.hh.

79  {
80 // calculate mass yield
81  G4double Ap13 = G4Pow::GetInstance()->powA(Ap, 1./3.);
82  G4double At13 = G4Pow::GetInstance()->powA(At, 1./3.);
83  G4double S = p_S[2] * (At13 + Ap13 + p_S[1]);
84 // cout << "debug0 "<<S<<" "<<At13<<" "<<Ap13<<" "<<p_S[1]<<" "<<p_S[2]<<endl;
85  G4double p = G4Exp(p_P[2]*Ap + p_P[1]);
86  G4double yield_a = p * S * G4Exp(-p * (Ap - A));
87  cout << "debug1 "<<yield_a<<endl;
88 // modification close to projectile
89  G4double f_mod_y=1.0;
90  if (A/Ap > corr_y[2])
91  {
92  f_mod_y=corr_y[1]*G4Pow::GetInstance()->powN(A/Ap-corr_y[2], 2) + 1.0;
93  }
94  yield_a= yield_a * f_mod_y;
95  cout << "debug1 "<<yield_a<<endl;
96 
97 // calculate maximum of charge dispersion zprob
98  G4double zbeta = A/(1.98+0.0155*G4Pow::GetInstance()->powA(A, (2./3.)));
99  G4double zbeta_p = Ap/(1.98+0.0155*G4Pow::GetInstance()->powA(Ap, (2./3.)));
100  G4double delta;
101  if(A > p_Delta[4])
102  {
103  delta = p_Delta[1] + p_Delta[2]*A;
104  }
105  else
106  {
107  delta = p_Delta[3]*A*A;
108  }
109 
110 // modification close to projectile
111  G4double f_mod=1.0;
112  if(A/Ap > corr_d[2])
113  {
114  f_mod = corr_d[1]*G4Pow::GetInstance()->powN(A/Ap-corr_d[2], 2) + 1.0;
115  }
116  delta = delta*f_mod;
117  G4double zprob = zbeta+delta;
118 
119 // correction for proton- and neutron-rich projectiles
120  G4double dq;
121  if((Zp-zbeta_p)>0)
122  {
123  dq = G4Exp(p_mp[1] + G4double(A)/G4double(Ap)*p_mp[2]);
124  cout << "dq "<<A<<" "<<Ap<<" "<<p_mp[1]
125  <<" "<<p_mp[2]<<" "<<dq<<" "<<p_mp[1] + A/Ap*p_mp[2]<<endl;
126  }
127  else
128  {
129  dq = p_mn[1]*G4Pow::GetInstance()->powN(A/Ap, 2) + p_mn[2]*G4Pow::GetInstance()->powN(A/Ap, 4);
130  }
131  zprob = zprob + dq * (Zp-zbeta_p);
132 
133 // small corr. since Xe-129 and Pb-208 are not on Z_beta line
134  zprob = zprob + 0.0020*A;
135  cout <<"zprob "<<A<<" "<<dq<<" "<<Zp<<" "<<zbeta_p
136  <<" "<<zbeta<<" "<<delta<<endl;
137 
138 // calculate width parameter R
139  G4double r = G4Exp(p_R[1] + p_R[2]*A);
140 
141 // modification close to projectile
142  f_mod=1.0;
143  if (A/Ap > corr_r[2])
144  {
145  f_mod = corr_r[1]*Ap*G4Pow::GetInstance()->powN(A/Ap-corr_r[2], 4)+1.0;
146  }
147  r = r*f_mod;
148 
149 // change width according to dev. from beta-stability
150  if ((Zp-zbeta_p) < 0.0)
151  {
152  r=r*(1.0-0.0833*std::abs(Zp-zbeta_p));
153  }
154 
155 // calculate slope parameters u_n, u_p
156  G4double u_n = p_Un[1];
157  G4double u_p = p_Up[1] + p_Up[2]*A + p_Up[3]*A*A;
158 
159 // calculate charge dispersion
160  G4double expo, fract;
161  if((zprob-Z) > 0)
162  {
163 // neutron-rich
164  expo = -r*G4Pow::GetInstance()->powA(std::abs(zprob-Z), u_n);
165  fract = G4Exp(expo)*std::sqrt(r/3.14159);
166  }
167  else
168  {
169 // proton-rich
170  expo = -r*G4Pow::GetInstance()->powA(std::abs(zprob-Z), u_p);
171  fract = G4Exp(expo)*std::sqrt(r/3.14159);
172  cout << "1 "<<expo<<" "<<r<<" "<<zprob<<" "<<Z<<" "<<u_p<<endl;
173 // go to exponential slope
174  G4double dfdz = 1.2 + 0.647*G4Pow::GetInstance()->powA(A/2.,0.3);
175  G4double z_exp = zprob + dfdz * G4Log(10.) / (2.*r);
176  if( Z>z_exp )
177  {
178  expo = -r*G4Pow::GetInstance()->powA(std::abs(zprob-z_exp), u_p);
179  fract = G4Exp(expo)*std::sqrt(r/3.14159)
180  / G4Pow::GetInstance()->powA(G4Pow::GetInstance()->powA(10, dfdz), Z-z_exp);
181  }
182  }
183 
184  cout << "debug "<<fract<<" "<<yield_a<<endl;
185  G4double epaxv2=fract*yield_a;
186  return epaxv2;
187  }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:128
double S(double temp)
const char * p
Definition: xmltok.h:285
double A(double temperature)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ProjectileFragmentCrossSection::testMe ( )
inline

Definition at line 189 of file G4ProjectileFragmentCrossSection.hh.

190  {
192  cout << i.doit(58, 28, 9, 4, 49, 28) << endl;
193  // Sigma = 9.800163E-13 b
194  }
G4double doit(G4double Ap, G4double Zp, G4double At, G4double Zt, G4double A, G4double Z)

Here is the call graph for this function:


The documentation for this class was generated from the following file: