ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ProjectileFragmentCrossSection.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ProjectileFragmentCrossSection.hh
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 #include "G4Exp.hh"
33 #include "G4Log.hh"
34 #include "G4Pow.hh"
35 
36 // Implements Physical Review C61, 034607 (2000)
37 // Rewrite starting from EPAX Version 2
38 
39 class G4ProjectileFragmentCrossSection
40 {
41  public:
42  G4ProjectileFragmentCrossSection()
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  }
77 
78  inline G4double doit(G4double Ap, G4double Zp, G4double At, G4double Zt, G4double A, G4double Z)
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 // std::cout << "debug0 "<<S<<" "<<At13<<" "<<Ap13<<" "<<p_S[1]<<" "<<p_S[2]<<std::endl;
85  G4double p = G4Exp(p_P[2]*Ap + p_P[1]);
86  G4double yield_a = p * S * G4Exp(-p * (Ap - A));
87 // std::cout << "debug1 "<<yield_a<<std::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 // std::cout << "debug1 "<<yield_a<<std::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 // std::cout << "dq "<<A<<" "<<Ap<<" "<<p_mp[1]
125 // <<" "<<p_mp[2]<<" "<<dq<<" "<<p_mp[1] + A/Ap*p_mp[2]<<std::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 // std::cout <<"zprob "<<A<<" "<<dq<<" "<<Zp<<" "<<zbeta_p
136 // <<" "<<zbeta<<" "<<delta<<std::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 // std::cout << "1 "<<expo<<" "<<r<<" "<<zprob<<" "<<Z<<" "<<u_p<<std::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 // std::cout << "debug "<<fract<<" "<<yield_a<<std::endl;
185  G4double epaxv2=fract*yield_a;
186  return epaxv2;
187  }
188 
189  void testMe()
190  {
191  G4ProjectileFragmentCrossSection i;
192 // std::cout << i.doit(58, 28, 9, 4, 49, 28) << std::endl;
193 // Sigma = 9.800163E-13 b
194  }
195  private:
196  G4double p_S[3];
197  G4double p_P[3];
198  G4double p_Delta[5];
199  G4double p_R[3];
200  G4double p_Un[2];
201  G4double p_Up[4];
202  G4double p_mn[3];
203  G4double p_mp[3];
204  G4double corr_d[3];
205  G4double corr_r[3];
206  G4double corr_y[3];
207 };
208 #endif