Geant4  10.02
G4KalbachCrossSection.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 // $Id: G4KalbachCrossSection.hh 66241 2012-12-13 18:34:42Z gunter $
28 //
29 // V.Ivanchenko 13.04.2015
30 //
31 // J.M. Quesada 22.04.2015 several fixes
32 
33 #ifndef G4KalbachCrossSection_h
34 #define G4KalbachCrossSection_h 1
35 
36 #include "globals.hh"
37 #include "G4Exp.hh"
38 #include "G4Pow.hh"
39 
40 //from subroutine sigpar of PRECO-2000 by Constance Kalbach Walker
41 // Calculate optical model reaction cross sections
42 // using the empirical parameterization
43 // of Narasimha Murthy, Chaterjee, and Gupta
44 // going over to the geometrical limit at high energy.
45 //
46 // Proton cross sections scaled down with signor for a<100
47 // (appropriate for becchetti-greenlees potential).
48 // p2 reduced and global red'n factor introduced below Bc
49 // Neutron cross sections scaled down with signor for a<40
50 // Scaled up for A>210 (added June '98 to conform with
51 // my published papers)
52 // (appropriate for Mani et al potential)
53 //
54 
55 static const G4double flow = 1.e-18;
56 static const G4double spill= 1.e+18;
57 
58 // index: 0-neutron, 1-proton, 2-deuteron, 3-triton, 4-He3, 5-He4
59 // parameters: p0, p1, p2, lambda0, lambda1, mu0, mu1, nu0, nu1, nu2, ra
60 
61 static const G4double paramK[6][11] = {
62 // n from mani, melkanoff and iori
63  {-312., 0., 0., 12.10, -11.27, 234.1, 38.26, 1.55, -106.1, 1280.8, 0.0},
64 // p from becchetti and greenlees (but modified with sub-barrier
65 // correction function and p2 changed from -449)
66  {15.72, 9.65, -300., 0.00437,-16.58, 244.7, 0.503, 273.1, -182.4, -1.872, 0.0},
67 // d from o.m. of perey and perey
68  {0.798, 420.3,-1651., 0.00619, -7.54, 583.5, 0.337, 421.8, -474.5, -3.592, 0.8},
69 // t from o.m. of hafele, flynn et al
70  {-21.45,484.7,-1608., 0.0186, -8.9, 686.3, 0.325, 368.9, -522.2, -4.998, 0.8},
71 // 3he from o.m. of gibson et al
72  {-2.88,205.6, -1487.,0.00459,-8.93, 611.2, 0.35 , 473.8, -468.2, -2.225, 0.8},
73 // alpha from huizenga and igo
74  { 10.95,-85.2, 1146., 0.0643,-13.96, 781.2, 0.29, -304.7,-470.0, -8.580, 1.2}
75 };
76 
78 {
79 public:
80 
82  {
83  return G4Pow::GetInstance()->powZ(resA, paramK[idx][6]);
84  }
85 
87  G4int idx, G4int Z, G4int A,
88  G4int resZ, G4int resA)
89  {
90  G4double sig = 0.0;
91  G4double signor = 1.0;
92  G4double lambda, mu, nu;
93  G4double ec = 0.5;
94  if(0 < Z) {
95  //JMQ 13.02.2009 tuning for improving cluster emission ddxs
96  // (spallation benchmark)
97  G4double xx = 1.7;
98  if(1 == A) { xx = 1.5; }
99  ec = 1.44 * Z * resZ / (xx*resA13 + paramK[idx][10]);
100  }
101  G4double ecsq = ec*ec;
102  G4double elab = K * (A + resA) / G4double(resA);
103 
104  if(idx == 0) { // parameterization for neutron
105 
106  if(resA < 40) { signor =0.7 + resA*0.0075; }
107  else if(resA > 210) { signor = 1. + (resA-210)*0.004; }
108  lambda = paramK[idx][3]/resA13 + paramK[idx][4];
109  mu = paramK[idx][5]*resA13 + paramK[idx][6]*resA13*resA13;
110  // JMQ 20.11.2008 very low energy behaviour corrected
111  // (problem for A (apprx.)>60) fix for avoiding
112  // neutron xs going to zero at very low energies
113  nu = std::abs(paramK[idx][7]*resA13*resA + paramK[idx][8]*resA13*resA13
114  + paramK[idx][9]);
115 
116  } else { // parameterization for charged
117  // proton correction
118  if(idx == 1) {
119  if (resA <= 60) { signor = 0.92; }
120  else if (resA < 100) { signor = 0.8 + resA*0.002; }
121  }
122  lambda = paramK[idx][3]*resA + paramK[idx][4];
123  mu = paramK[idx][5]*amu1;
124  nu = amu1* (paramK[idx][7] + paramK[idx][8]*ec + paramK[idx][9]*ecsq);
125  }
126 
127  // threashold cross section
128  if(elab <= ec) {
129  G4double p = paramK[idx][0];
130  if(0 < Z) { p += paramK[idx][1]/ec + paramK[idx][2]/ecsq; }
131  G4double a = -2*p*ec + lambda - nu/ecsq;
132  G4double b = p*ecsq + mu + 2*nu/ec;
133  G4double ecut;
134  G4double det = a*a - 4*p*b;
135  if (det > 0.0) { ecut = (std::sqrt(det) - a)/(p + p); }
136  else { ecut = a/(p + p); }
137 
138  // If ecut>0, sig=0 at elab=ecut
139  if(elab > ecut) {
140  sig = (p*elab*elab + a*elab + b)*signor;
141 
142  // extra proton correction
143  if(1 == idx) {
144  // c and w are for global correction factor for
145  // they are scaled down for light targets where ec is low.
146  G4double cc = std::min(3.15, ec*0.5);
147  G4double signor2 = (ec - elab - cc) *3.15/ (0.7*cc);
148  sig /= (1. + G4Exp(signor2));
149  }
150  }
151 
152  // high energy cross section
153  } else {
154  // etest is the energy above which the rxn cross section is
155  // compared with the geometrical limit and the max taken.
156 
157  // neutron parameters
158  G4double etest = 32.;
159  G4double xnulam = 1.0;
160  // parameters for charged
161  if(0 < Z) {
162  etest = 0.0;
163  xnulam = nu / lambda;
164  if(xnulam > spill) { xnulam= 0.0; }
165  else if (xnulam >= flow) {
166  if(1 == idx) { etest = std::sqrt(xnulam) + 7.; }
167  else { etest = 1.2 *std::sqrt(xnulam); }
168  }
169  }
170  // ** For xnulam.gt.0, sig reaches a maximum at sqrt(xnulam).
171  sig = (lambda*elab + mu + nu/elab)*signor;
172  if (xnulam >= flow && elab >= etest) {
173  G4double geom = std::sqrt(A*K);
174  geom = 1.23*resA13 + paramK[idx][10] + 4.573/geom;
175  geom = 31.416 * geom * geom;
176  sig = std::max(sig, geom);
177  }
178  }
179  sig = std::max(sig, 0.0);
180  return sig;
181  }
182 };
183 
184 #endif
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static G4double ComputePowerParameter(G4int resA, G4int idx)
G4double a
Definition: TRTMaterials.hh:39
int G4int
Definition: G4Types.hh:78
double A(double temperature)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:254
double G4double
Definition: G4Types.hh:76
static const G4double paramK[6][11]
static const G4double spill
static G4double ComputeCrossSection(G4double K, G4double resA13, G4double amu1, G4int idx, G4int Z, G4int A, G4int resZ, G4int resA)
static const G4double flow