Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ProjectileFragmentCrossSection.hh
Go to the documentation of this file.
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 #ifndef G4ProjectileFragmentCrossSection_h
28 #define G4ProjectileFragmentCrossSection_h 1
29 
30 #include <cmath>
31 #include <iostream>
32 
33 // Implements Physical Review C61, 034607 (2000)
34 // Rewrite starting from EPAX Version 2
35 
36 class G4ProjectileFragmentCrossSection
37 {
38  public:
39  G4ProjectileFragmentCrossSection()
40  {
41  p_S[1] = -2.38; // scale factor for xsect in barn
42  p_S[2] = 0.27;
43 
44  p_P[1] = -2.5840E+00; // slope of mass yield curve
45  p_P[2] = -7.5700E-03;
46 
47  p_Delta[1] = -1.0870E+00; // centroid rel. to beta-stability
48  p_Delta[2] = +3.0470E-02;
49  p_Delta[3] = +2.1353E-04;
50  p_Delta[4] = +7.1350E+01;
51 
52  p_R[1] = +0.885E+00; // width parameter R
53  p_R[2] = -9.8160E-03;
54 
55  p_Un[1] = 1.65; // slope par. n-rich ride of Z distr.
56 
57  p_Up[1] = 1.7880; // slope par. p-rich ride of Z distr.
58  p_Up[2] = +4.7210E-03;
59  p_Up[3] = -1.3030E-05;
60 
61  p_mn[1] = 0.400; // memory effect n-rich projectiles
62  p_mn[2] = 0.600;
63 
64  p_mp[1] = -10.25; // memory effect p-rich projectiles
65  p_mp[2] = +10.1;
66 
67  corr_d[1] = -25.0; // correction close to proj.: centroid dzp
68  corr_d[2] = 0.800;
69  corr_r[1] = +20.0; // correction close to proj.: width R
70  corr_r[2] = 0.820;
71  corr_y[1] = 200.0; // correction close to proj.: Yield_a
72  corr_y[2] = 0.90;
73  }
74 
75  inline G4double doit(G4double Ap, G4double Zp, G4double At, G4double Zt, G4double A, G4double Z)
76  {
77 // calculate mass yield
78  G4double Ap13 = std::pow(Ap, 1./3.);
79  G4double At13 = std::pow(At, 1./3.);
80  G4double S = p_S[2] * (At13 + Ap13 + p_S[1]);
81 // cout << "debug0 "<<S<<" "<<At13<<" "<<Ap13<<" "<<p_S[1]<<" "<<p_S[2]<<endl;
82  G4double p = std::exp(p_P[2]*Ap + p_P[1]);
83  G4double yield_a = p * S * std::exp(-p * (Ap - A));
84  cout << "debug1 "<<yield_a<<endl;
85 // modification close to projectile
86  G4double f_mod_y=1.0;
87  if (A/Ap > corr_y[2])
88  {
89  f_mod_y=corr_y[1]*std::pow(A/Ap-corr_y[2], 2) + 1.0;
90  }
91  yield_a= yield_a * f_mod_y;
92  cout << "debug1 "<<yield_a<<endl;
93 
94 // calculate maximum of charge dispersion zprob
95  G4double zbeta = A/(1.98+0.0155*std::pow(A, (2./3.)));
96  G4double zbeta_p = Ap/(1.98+0.0155*std::pow(Ap, (2./3.)));
97  G4double delta;
98  if(A > p_Delta[4])
99  {
100  delta = p_Delta[1] + p_Delta[2]*A;
101  }
102  else
103  {
104  delta = p_Delta[3]*A*A;
105  }
106 
107 // modification close to projectile
108  G4double f_mod=1.0;
109  if(A/Ap > corr_d[2])
110  {
111  f_mod = corr_d[1]*std::pow(A/Ap-corr_d[2], 2) + 1.0;
112  }
113  delta = delta*f_mod;
114  G4double zprob = zbeta+delta;
115 
116 // correction for proton- and neutron-rich projectiles
117  G4double dq;
118  if((Zp-zbeta_p)>0)
119  {
120  dq = std::exp(p_mp[1] + G4double(A)/G4double(Ap)*p_mp[2]);
121  cout << "dq "<<A<<" "<<Ap<<" "<<p_mp[1]
122  <<" "<<p_mp[2]<<" "<<dq<<" "<<p_mp[1] + A/Ap*p_mp[2]<<endl;
123  }
124  else
125  {
126  dq = p_mn[1]*std::pow(A/Ap, 2.0) + p_mn[2]*std::pow(A/Ap, 4.0);
127  }
128  zprob = zprob + dq * (Zp-zbeta_p);
129 
130 // small corr. since Xe-129 and Pb-208 are not on Z_beta line
131  zprob = zprob + 0.0020*A;
132  cout <<"zprob "<<A<<" "<<dq<<" "<<Zp<<" "<<zbeta_p
133  <<" "<<zbeta<<" "<<delta<<endl;
134 
135 // calculate width parameter R
136  G4double r = std::exp(p_R[1] + p_R[2]*A);
137 
138 // modification close to projectile
139  f_mod=1.0;
140  if (A/Ap > corr_r[2])
141  {
142  f_mod = corr_r[1]*Ap*std::pow(A/Ap-corr_r[2], 4.0)+1.0;
143  }
144  r = r*f_mod;
145 
146 // change width according to dev. from beta-stability
147  if ((Zp-zbeta_p) < 0.0)
148  {
149  r=r*(1.0-0.0833*std::abs(Zp-zbeta_p));
150  }
151 
152 // calculate slope parameters u_n, u_p
153  G4double u_n = p_Un[1];
154  G4double u_p = p_Up[1] + p_Up[2]*A + p_Up[3]*A*A;
155 
156 // calculate charge dispersion
157  G4double expo, fract;
158  if((zprob-Z) > 0)
159  {
160 // neutron-rich
161  expo = -r*std::pow(std::abs(zprob-Z), u_n);
162  fract = std::exp(expo)*std::sqrt(r/3.14159);
163  }
164  else
165  {
166 // proton-rich
167  expo = -r*std::pow(std::abs(zprob-Z), u_p);
168  fract = std::exp(expo)*std::sqrt(r/3.14159);
169  cout << "1 "<<expo<<" "<<r<<" "<<zprob<<" "<<Z<<" "<<u_p<<endl;
170 // go to exponential slope
171  G4double dfdz = 1.2 + 0.647*std::pow(A/2.,0.3);
172  G4double z_exp = zprob + dfdz * std::log(10.) / (2.*r);
173  if( Z>z_exp )
174  {
175  expo = -r*std::pow(std::abs(zprob-z_exp), u_p);
176  fract = std::exp(expo)*std::sqrt(r/3.14159)
177  / std::pow(std::pow(10, dfdz), Z-z_exp);
178  }
179  }
180 
181  cout << "debug "<<fract<<" "<<yield_a<<endl;
182  G4double epaxv2=fract*yield_a;
183  return epaxv2;
184  }
185 
186  void testMe()
187  {
188  G4ProjectileFragmentCrossSection i;
189  cout << i.doit(58, 28, 9, 4, 49, 28) << endl;
190  // Sigma = 9.800163E-13 b
191  }
192  private:
193  G4double p_S[3];
194  G4double p_P[3];
195  G4double p_Delta[5];
196  G4double p_R[3];
197  G4double p_Un[2];
198  G4double p_Up[4];
199  G4double p_mn[3];
200  G4double p_mp[3];
201  G4double corr_d[3];
202  G4double corr_r[3];
203  G4double corr_y[3];
204 };
205 #endif
Generated on Sat May 25 2013 14:33:40 for Geant4 by   doxygen 1.8.4