| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 // V.Ivanchenko 26.02.2018 28 // V.Ivanchenko 26.02.2018
29 // 29 //
30 30
31 #include "G4KalbachCrossSection.hh" 31 #include "G4KalbachCrossSection.hh"
32 #include "G4Exp.hh" 32 #include "G4Exp.hh"
33 #include "G4Pow.hh" 33 #include "G4Pow.hh"
34 34
35 //from subroutine sigpar of PRECO-2000 by Con 35 //from subroutine sigpar of PRECO-2000 by Constance Kalbach Walker
36 // Calculate optical model reaction cross 36 // Calculate optical model reaction cross sections
37 // using the empirical parameterization 37 // using the empirical parameterization
38 // of Narasimha Murthy, Chaterjee, and Gup 38 // of Narasimha Murthy, Chaterjee, and Gupta
39 // going over to the geometrical limit at 39 // going over to the geometrical limit at high energy.
40 // 40 //
41 // Proton cross sections scaled down 41 // Proton cross sections scaled down with signor for a<100
42 // (appropriate for becchetti-greenl 42 // (appropriate for becchetti-greenlees potential).
43 // p2 reduced and global red'n facto 43 // p2 reduced and global red'n factor introduced below Bc
44 // Neutron cross sections scaled dow 44 // Neutron cross sections scaled down with signor for a<40
45 // Scaled up for A>210 (added June ' 45 // Scaled up for A>210 (added June '98 to conform with
46 // my published papers) 46 // my published papers)
47 // (appropriate for Mani et al poten 47 // (appropriate for Mani et al potential)
48 // 48 //
49 49
50 // index: 0-neutron, 1-proton, 2-deuteron, 3-t 50 // index: 0-neutron, 1-proton, 2-deuteron, 3-triton, 4-He3, 5-He4
51 // parameters: p0, p1, p2, lambda0, lambda1, m 51 // parameters: p0, p1, p2, lambda0, lambda1, mu0, mu1, nu0, nu1, nu2, ra
52 52
53 static const G4double paramK[6][11] = { 53 static const G4double paramK[6][11] = {
54 // n from mani, melkanoff and iori 54 // n from mani, melkanoff and iori
55 {-312., 0., 0., 12.10, -11.27, 234.1, 55 {-312., 0., 0., 12.10, -11.27, 234.1, 38.26, 1.55, -106.1, 1280.8, 0.0},
56 // p from becchetti and greenlees (but modi 56 // p from becchetti and greenlees (but modified with sub-barrier
57 // correction function and p2 changed from - 57 // correction function and p2 changed from -449)
58 {15.72, 9.65, -300., 0.00437,-16.58, 244.7, 58 {15.72, 9.65, -300., 0.00437,-16.58, 244.7, 0.503, 273.1, -182.4, -1.872, 0.0},
59 // d from o.m. of perey and perey 59 // d from o.m. of perey and perey
60 {0.798, 420.3,-1651., 0.00619, -7.54, 583.5, 60 {0.798, 420.3,-1651., 0.00619, -7.54, 583.5, 0.337, 421.8, -474.5, -3.592, 0.8},
61 // t from o.m. of hafele, flynn et al 61 // t from o.m. of hafele, flynn et al
62 {-21.45,484.7,-1608., 0.0186, -8.9, 686.3, 62 {-21.45,484.7,-1608., 0.0186, -8.9, 686.3, 0.325, 368.9, -522.2, -4.998, 0.8},
63 // 3he from o.m. of gibson et al 63 // 3he from o.m. of gibson et al
64 {-2.88,205.6, -1487.,0.00459,-8.93, 611.2, 64 {-2.88,205.6, -1487.,0.00459,-8.93, 611.2, 0.35 , 473.8, -468.2, -2.225, 0.8},
65 // alpha from huizenga and igo 65 // alpha from huizenga and igo
66 { 10.95,-85.2, 1146., 0.0643,-13.96, 781.2, 66 { 10.95,-85.2, 1146., 0.0643,-13.96, 781.2, 0.29, -304.7,-470.0, -8.580, 1.2}
67 }; 67 };
68 68
69 G4double G4KalbachCrossSection::ComputePowerPa 69 G4double G4KalbachCrossSection::ComputePowerParameter(G4int resA, G4int idx)
70 { 70 {
71 return G4Pow::GetInstance()->powZ(resA, para 71 return G4Pow::GetInstance()->powZ(resA, paramK[idx][6]);
72 } 72 }
73 73
74 G4double 74 G4double
75 G4KalbachCrossSection::ComputeCrossSection(G4d 75 G4KalbachCrossSection::ComputeCrossSection(G4double K, G4double cb,
76 G4double resA13, G4double a 76 G4double resA13, G4double amu1,
77 G4int idx, G4int Z, G4int A 77 G4int idx, G4int Z, G4int A,
78 G4int resA) 78 G4int resA)
79 { 79 {
80 G4double sig = 0.0; 80 G4double sig = 0.0;
81 G4double signor = 1.0; 81 G4double signor = 1.0;
82 G4double lambda, mu, nu; 82 G4double lambda, mu, nu;
83 G4double ec = std::min(4.0, 100./(G4double)r << 83 G4double ec = 0.5;
84 if(0 < Z) { ec = cb; } 84 if(0 < Z) { ec = cb; }
85 << 85 //JMQ 13.02.2009 tuning for improving cluster emission ddxs
>> 86 // (spallation benchmark)
>> 87 /*
>> 88 G4double xx = 1.7;
>> 89 if(1 == A) { xx = 1.5; }
>> 90 ec = 1.44 * Z * resZ / (xx*resA13 + paramK[idx][10]);
>> 91 }
>> 92 */
86 G4double ecsq = ec*ec; 93 G4double ecsq = ec*ec;
87 G4double elab = K * (A + resA) / G4double(re 94 G4double elab = K * (A + resA) / G4double(resA);
88 95
89 if(idx == 0) { // parameterization for neutr 96 if(idx == 0) { // parameterization for neutron
90 97
91 if(resA < 40) { signor =0.7 + resA*0 98 if(resA < 40) { signor =0.7 + resA*0.0075; }
92 else if(resA > 210) { signor = 1. + (resA- 99 else if(resA > 210) { signor = 1. + (resA-210)*0.004; }
93 lambda = paramK[idx][3]/resA13 + paramK[id 100 lambda = paramK[idx][3]/resA13 + paramK[idx][4];
94 mu = (paramK[idx][5] + paramK[idx][6]*resA 101 mu = (paramK[idx][5] + paramK[idx][6]*resA13)*resA13;
95 // JMQ 20.11.2008 very low energy behaviou 102 // JMQ 20.11.2008 very low energy behaviour corrected
96 // (problem for A (apprx.)> 103 // (problem for A (apprx.)>60) fix for avoiding
97 // neutron xs going to zero 104 // neutron xs going to zero at very low energies
98 nu = std::abs((paramK[idx][7]*resA + param 105 nu = std::abs((paramK[idx][7]*resA + paramK[idx][8]*resA13)*resA13
99 + paramK[idx][9]); 106 + paramK[idx][9]);
100 107
101 } else { // parameterization for charged 108 } else { // parameterization for charged
102 // proton correction 109 // proton correction
103 if(idx == 1) { 110 if(idx == 1) {
104 if (resA <= 60) { signor = 0.92; } 111 if (resA <= 60) { signor = 0.92; }
105 else if (resA < 100) { signor = 0.8 + re 112 else if (resA < 100) { signor = 0.8 + resA*0.002; }
106 } 113 }
107 lambda = paramK[idx][3]*resA + paramK[idx] 114 lambda = paramK[idx][3]*resA + paramK[idx][4];
108 mu = paramK[idx][5]*amu1; 115 mu = paramK[idx][5]*amu1;
109 nu = amu1* (paramK[idx][7] + paramK[idx][8 116 nu = amu1* (paramK[idx][7] + paramK[idx][8]*ec + paramK[idx][9]*ecsq);
110 } 117 }
111 /* 118 /*
112 G4cout << "## idx=" << idx << " K=" << K < << 119 G4cout << "## idx= " << idx << " K= " << K << " elab= " << elab << " ec= " << ec
113 << " ec=" << ec << " lambda=" << la << 120 << " lambda= " << lambda << " mu= " << mu << " nu= " << nu << G4endl;
114 << " nu=" << nu << G4endl; <<
115 */ 121 */
116 // threashold cross section 122 // threashold cross section
117 if(elab < ec) { 123 if(elab < ec) {
118 G4double p = paramK[idx][0]; 124 G4double p = paramK[idx][0];
119 if(0 < Z) { p += paramK[idx][1]/ec + param 125 if(0 < Z) { p += paramK[idx][1]/ec + paramK[idx][2]/ecsq; }
120 G4double a = -2*p*ec + lambda - nu/ecsq; 126 G4double a = -2*p*ec + lambda - nu/ecsq;
121 G4double b = p*ecsq + mu + 2*nu/ec; 127 G4double b = p*ecsq + mu + 2*nu/ec;
>> 128 G4double ecut;
122 G4double det = a*a - 4*p*b; 129 G4double det = a*a - 4*p*b;
123 G4double ecut = (det > 0.0) ? (std::sqrt(d << 130 if (det > 0.0) { ecut = (std::sqrt(det) - a)/(2*p); }
>> 131 else { ecut = -a/(2*p); }
124 132
125 //G4cout << " elab= " << elab << " ecut= 133 //G4cout << " elab= " << elab << " ecut= " << ecut << " sig= " << sig
126 // << " sig1= " << (p*elab*elab + a 134 // << " sig1= " << (p*elab*elab + a*elab + b)*signor << G4endl;
127 // If ecut>0, sig=0 at elab=ecut 135 // If ecut>0, sig=0 at elab=ecut
128 if(0 == idx) { 136 if(0 == idx) {
129 sig = (lambda*ec + mu + nu/ec)*signor*st 137 sig = (lambda*ec + mu + nu/ec)*signor*std::sqrt(elab/ec);
130 } else if(elab >= ecut) { 138 } else if(elab >= ecut) {
131 sig = (p*elab*elab + a*elab + b)*signor; 139 sig = (p*elab*elab + a*elab + b)*signor;
132 140
133 // extra proton correction 141 // extra proton correction
134 if(1 == idx) { 142 if(1 == idx) {
135 // c and w are for global correction factor 143 // c and w are for global correction factor for
136 // they are scaled down for light targets wh 144 // they are scaled down for light targets where ec is low.
137 G4double cc = std::min(3.15, ec*0.5); 145 G4double cc = std::min(3.15, ec*0.5);
138 G4double signor2 = (ec - elab - cc) *3.15/ ( 146 G4double signor2 = (ec - elab - cc) *3.15/ (0.7*cc);
139 sig /= (1. + G4Exp(signor2)); 147 sig /= (1. + G4Exp(signor2));
140 } 148 }
141 } 149 }
142 //G4cout << " ecut= " << ecut << " a 150 //G4cout << " ecut= " << ecut << " a= " << a << " b= " << b
143 // << " signor= " << signor << " sig= 151 // << " signor= " << signor << " sig= " << sig << G4endl;
144 152
145 // high energy cross section 153 // high energy cross section
146 } else { 154 } else {
147 // etest is the energy above which the rxn 155 // etest is the energy above which the rxn cross section is
148 // compared with the geometrical limit and 156 // compared with the geometrical limit and the max taken.
149 157
150 // neutron parameters 158 // neutron parameters
151 G4double etest = 32.; 159 G4double etest = 32.;
152 G4double xnulam = 1.0; 160 G4double xnulam = 1.0;
153 161
154 // parameters for charged 162 // parameters for charged
155 static const G4double flow = 1.e-18; 163 static const G4double flow = 1.e-18;
156 static const G4double spill= 1.e+18; 164 static const G4double spill= 1.e+18;
157 if(0 < Z) { 165 if(0 < Z) {
158 etest = 0.0; 166 etest = 0.0;
159 xnulam = nu / lambda; 167 xnulam = nu / lambda;
160 xnulam = std::min(xnulam, spill); 168 xnulam = std::min(xnulam, spill);
161 if (xnulam >= flow) { 169 if (xnulam >= flow) {
162 if(1 == idx) { etest = std::sqrt(xnulam) + 7 170 if(1 == idx) { etest = std::sqrt(xnulam) + 7.; }
163 else { etest = 1.2 *std::sqrt(xnulam 171 else { etest = 1.2 *std::sqrt(xnulam); }
164 } 172 }
165 } 173 }
166 // ** For xnulam.gt.0, sig reaches a maxim 174 // ** For xnulam.gt.0, sig reaches a maximum at sqrt(xnulam).
167 sig = (lambda*elab + mu + nu/elab)*signor; 175 sig = (lambda*elab + mu + nu/elab)*signor;
168 if (xnulam >= flow && elab >= etest) { 176 if (xnulam >= flow && elab >= etest) {
169 G4double geom = std::sqrt(A*K); 177 G4double geom = std::sqrt(A*K);
170 geom = 1.23*resA13 + paramK[idx][10] + 4 178 geom = 1.23*resA13 + paramK[idx][10] + 4.573/geom;
171 geom = 31.416 * geom * geom; 179 geom = 31.416 * geom * geom;
172 sig = std::max(sig, geom); 180 sig = std::max(sig, geom);
173 } 181 }
174 } 182 }
175 sig = std::max(sig, 0.0); 183 sig = std::max(sig, 0.0);
176 //G4cout << " ---- sig= " << sig << G4endl; 184 //G4cout << " ---- sig= " << sig << G4endl;
177 return sig; 185 return sig;
178 } 186 }
179 187