| 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 // * 20 // * *
21 // * Parts of this code which have been devel 21 // * Parts of this code which have been developed by QinetiQ Ltd *
22 // * under contract to the European Space Agen 22 // * under contract to the European Space Agency (ESA) are the *
23 // * intellectual property of ESA. Rights to u 23 // * intellectual property of ESA. Rights to use, copy, modify and *
24 // * redistribute this software for general pu 24 // * redistribute this software for general public use are granted *
25 // * in compliance with any licensing, distrib 25 // * in compliance with any licensing, distribution and development *
26 // * policy adopted by the Geant4 Collaboratio 26 // * policy adopted by the Geant4 Collaboration. This code has been *
27 // * written by QinetiQ Ltd for the European S 27 // * written by QinetiQ Ltd for the European Space Agency, under ESA *
28 // * contract 17191/03/NL/LvH (Aurora Programm 28 // * contract 17191/03/NL/LvH (Aurora Programme). *
29 // * 29 // * *
30 // * By using, copying, modifying or distri 30 // * By using, copying, modifying or distributing the software (or *
31 // * any work based on the software) you ag 31 // * any work based on the software) you agree to acknowledge its *
32 // * use in resulting scientific publicati 32 // * use in resulting scientific publications, and indicate your *
33 // * acceptance of all terms of the Geant4 Sof 33 // * acceptance of all terms of the Geant4 Software license. *
34 // ******************************************* 34 // ********************************************************************
35 // 35 //
36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37 // 37 //
38 // MODULE: G4NuclearAbrasionGeometry.cc 38 // MODULE: G4NuclearAbrasionGeometry.cc
39 // 39 //
40 // Version: B.1 40 // Version: B.1
41 // Date: 15/04/04 41 // Date: 15/04/04
42 // Author: P R Truscott 42 // Author: P R Truscott
43 // Organisation: QinetiQ Ltd, UK 43 // Organisation: QinetiQ Ltd, UK
44 // Customer: ESA/ESTEC, NOORDWIJK 44 // Customer: ESA/ESTEC, NOORDWIJK
45 // Contract: 17191/03/NL/LvH 45 // Contract: 17191/03/NL/LvH
46 // 46 //
47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 // 48 //
49 // CHANGE HISTORY 49 // CHANGE HISTORY
50 // -------------- 50 // --------------
51 // 51 //
52 // 18 November 2003, P R Truscott, QinetiQ Ltd 52 // 18 November 2003, P R Truscott, QinetiQ Ltd, UK
53 // Created. 53 // Created.
54 // 54 //
55 // 15 March 2004, P R Truscott, QinetiQ Ltd, U 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
56 // Beta release 56 // Beta release
57 // 57 //
58 // 4 June 2004, J.P. Wellisch, CERN, Switzerla 58 // 4 June 2004, J.P. Wellisch, CERN, Switzerland
59 // resolving technical portability issues. 59 // resolving technical portability issues.
60 // 60 //
61 // 12 June 2012, A. Ribon, CERN, Switzerland 61 // 12 June 2012, A. Ribon, CERN, Switzerland
62 // Fixing trivial warning errors of shadowed v 62 // Fixing trivial warning errors of shadowed variables.
63 // 63 //
64 // 4 August 2015, A. Ribon, CERN, Switzerland 64 // 4 August 2015, A. Ribon, CERN, Switzerland
65 // Replacing std::pow with the faster G4Pow. 65 // Replacing std::pow with the faster G4Pow.
66 // 66 //
67 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 67 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
68 ////////////////////////////////////////////// 68 ////////////////////////////////////////////////////////////////////////////////
69 // 69 //
70 #include "G4NuclearAbrasionGeometry.hh" 70 #include "G4NuclearAbrasionGeometry.hh"
71 #include "G4WilsonRadius.hh" 71 #include "G4WilsonRadius.hh"
72 #include "G4PhysicalConstants.hh" 72 #include "G4PhysicalConstants.hh"
73 #include "G4SystemOfUnits.hh" 73 #include "G4SystemOfUnits.hh"
74 #include "G4Pow.hh" 74 #include "G4Pow.hh"
75 ////////////////////////////////////////////// 75 ////////////////////////////////////////////////////////////////////////////////
76 // 76 //
77 G4NuclearAbrasionGeometry::G4NuclearAbrasionGe 77 G4NuclearAbrasionGeometry::G4NuclearAbrasionGeometry (G4double AP1,
78 G4double AT1, G4double r1) 78 G4double AT1, G4double r1)
79 { 79 {
80 // 80 //
81 // 81 //
82 // Initialise variables for interaction geomet 82 // Initialise variables for interaction geometry.
83 // 83 //
84 G4WilsonRadius aR; 84 G4WilsonRadius aR;
85 AP = AP1; 85 AP = AP1;
86 AT = AT1; 86 AT = AT1;
87 rP = aR.GetWilsonRadius(AP); 87 rP = aR.GetWilsonRadius(AP);
88 rT = aR.GetWilsonRadius(AT); 88 rT = aR.GetWilsonRadius(AT);
89 r = r1; 89 r = r1;
90 n = rP / (rP + rT); 90 n = rP / (rP + rT);
91 b = r / (rP + rT); 91 b = r / (rP + rT);
92 m = rT / rP; 92 m = rT / rP;
93 Q = (1.0 - b)/n; 93 Q = (1.0 - b)/n;
94 S = Q * Q; 94 S = Q * Q;
95 T = S * Q; 95 T = S * Q;
96 R = std::sqrt(m*n); 96 R = std::sqrt(m*n);
97 U = 1.0/m - 2.0; 97 U = 1.0/m - 2.0;
98 // 98 //
99 // 99 //
100 // Initialise the threshold radius-ratio at wh 100 // Initialise the threshold radius-ratio at which interactions are considered
101 // peripheral or central. 101 // peripheral or central.
102 // 102 //
103 rth = 2.0/3.0; 103 rth = 2.0/3.0;
104 B = 10.0 * MeV; 104 B = 10.0 * MeV;
105 } 105 }
106 ////////////////////////////////////////////// 106 ////////////////////////////////////////////////////////////////////////////////
107 // 107 //
108 G4NuclearAbrasionGeometry::~G4NuclearAbrasionG 108 G4NuclearAbrasionGeometry::~G4NuclearAbrasionGeometry ()
109 {;} 109 {;}
110 ////////////////////////////////////////////// 110 ////////////////////////////////////////////////////////////////////////////////
111 // 111 //
112 void G4NuclearAbrasionGeometry::SetPeripheralT 112 void G4NuclearAbrasionGeometry::SetPeripheralThreshold (G4double rth1)
113 {if (rth1 > 0.0 && rth1 <= 1.0) rth = rth1;} 113 {if (rth1 > 0.0 && rth1 <= 1.0) rth = rth1;}
114 ////////////////////////////////////////////// 114 ////////////////////////////////////////////////////////////////////////////////
115 // 115 //
116 G4double G4NuclearAbrasionGeometry::GetPeriphe 116 G4double G4NuclearAbrasionGeometry::GetPeripheralThreshold ()
117 {return rth;} 117 {return rth;}
118 ////////////////////////////////////////////// 118 ////////////////////////////////////////////////////////////////////////////////
119 // 119 //
120 G4double G4NuclearAbrasionGeometry::P () 120 G4double G4NuclearAbrasionGeometry::P ()
121 { 121 {
122 // 122 //
123 // 123 //
124 // Initialise the value for P, then determine 124 // Initialise the value for P, then determine the actual value depending upon
125 // whether the projectile is larger or smaller 125 // whether the projectile is larger or smaller than the target and these radii
126 // in relation to the impact parameter. 126 // in relation to the impact parameter.
127 // 127 //
128 G4double valueP = 0.0; 128 G4double valueP = 0.0;
129 129
130 if (rT > rP) 130 if (rT > rP)
131 { 131 {
132 if (rT-rP<=r && r<=rT+rP) valueP = 0.125*R 132 if (rT-rP<=r && r<=rT+rP) valueP = 0.125*R*U*S - 0.125*(0.5*R*U+1.0)*T;
133 else valueP = -1.0; 133 else valueP = -1.0;
134 } 134 }
135 else 135 else
136 { 136 {
137 if (rP-rT<=r && r<=rP+rT) valueP = 0.125*R 137 if (rP-rT<=r && r<=rP+rT) valueP = 0.125*R*U*S - 0.125*(0.5*std::sqrt(n/m)*U-
138 (std::sqrt(1.0-m*m)/n - 1.0)*std::sqrt(( 138 (std::sqrt(1.0-m*m)/n - 1.0)*std::sqrt((2.0-m)/G4Pow::GetInstance()->powN(m,5)))*T;
139 else valueP = (std::s 139 else valueP = (std::sqrt(1.0-m*m)/n-1.0)*std::sqrt(1.0-b*b/n/n);
140 } 140 }
141 141
142 if (!(valueP <= 1.0 && valueP>= -1.0)) 142 if (!(valueP <= 1.0 && valueP>= -1.0))
143 { 143 {
144 if (valueP > 1.0) valueP = 1.0; 144 if (valueP > 1.0) valueP = 1.0;
145 else valueP = -1.0; 145 else valueP = -1.0;
146 } 146 }
147 return valueP; 147 return valueP;
148 } 148 }
149 ////////////////////////////////////////////// 149 ////////////////////////////////////////////////////////////////////////////////
150 // 150 //
151 G4double G4NuclearAbrasionGeometry::F () 151 G4double G4NuclearAbrasionGeometry::F ()
152 { 152 {
153 // 153 //
154 // 154 //
155 // Initialise the value for F, then determine 155 // Initialise the value for F, then determine the actual value depending upon
156 // whether the projectile is larger or smaller 156 // whether the projectile is larger or smaller than the target and these radii
157 // in relation to the impact parameter. 157 // in relation to the impact parameter.
158 // 158 //
159 G4double valueF = 0.0; 159 G4double valueF = 0.0;
160 160
161 if (rT > rP) 161 if (rT > rP)
162 { 162 {
163 if (rT-rP<=r && r<=rT+rP) valueF = 0.75*R* 163 if (rT-rP<=r && r<=rT+rP) valueF = 0.75*R*S - 0.125*(3.0*R-1.0)*T;
164 else valueF = 1.0; 164 else valueF = 1.0;
165 } 165 }
166 else 166 else
167 { 167 {
168 if (rP-rT<=r && r<=rP+rT) valueF = 0.75*R* 168 if (rP-rT<=r && r<=rP+rT) valueF = 0.75*R*S - 0.125*(3.0*std::sqrt(n/m)-
169 (1.0-G4Pow::GetInstance()->powA(1.0-m*m, 169 (1.0-G4Pow::GetInstance()->powA(1.0-m*m,3.0/2.0))*std::sqrt(1.0-G4Pow::GetInstance()->powN(1.0-m,2))/G4Pow::GetInstance()->powN(m,3))*T;
170 else valueF = (1.0-G4 170 else valueF = (1.0-G4Pow::GetInstance()->powA(1.0-m*m,3.0/2.0))*std::sqrt(1.0-b*b/n/n);
171 } 171 }
172 172
173 if (!(valueF <= 1.0 && valueF>= 0.0)) 173 if (!(valueF <= 1.0 && valueF>= 0.0))
174 { 174 {
175 if (valueF > 1.0) valueF = 1.0; 175 if (valueF > 1.0) valueF = 1.0;
176 else valueF = 0.0; 176 else valueF = 0.0;
177 } 177 }
178 return valueF; 178 return valueF;
179 } 179 }
180 ////////////////////////////////////////////// 180 ////////////////////////////////////////////////////////////////////////////////
181 // 181 //
182 G4double G4NuclearAbrasionGeometry::GetExcitat 182 G4double G4NuclearAbrasionGeometry::GetExcitationEnergyOfProjectile ()
183 { 183 {
184 G4double F1 = F(); 184 G4double F1 = F();
185 G4double P1 = P(); 185 G4double P1 = P();
186 G4double Es = 0.0; 186 G4double Es = 0.0;
187 187
188 Es = 0.95 * MeV * 4.0 * pi * rP*rP/fermi/fer 188 Es = 0.95 * MeV * 4.0 * pi * rP*rP/fermi/fermi *
189 (1.0+P1-G4Pow::GetInstance()->A23(1.0-F 189 (1.0+P1-G4Pow::GetInstance()->A23(1.0-F1));
190 // if (rT < rP && r < rP-rT) 190 // if (rT < rP && r < rP-rT)
191 if ((r-rP)/rT < rth) 191 if ((r-rP)/rT < rth)
192 { 192 {
193 G4double omega = 0.0; 193 G4double omega = 0.0;
194 if (AP < 12.0) omega = 1500.0; 194 if (AP < 12.0) omega = 1500.0;
195 else if (AP <= 16.0) omega = 1500.0 - 320. 195 else if (AP <= 16.0) omega = 1500.0 - 320.0*(AP-12.0);
196 Es *= 1.0 + F1*(5.0+omega*F1*F1); 196 Es *= 1.0 + F1*(5.0+omega*F1*F1);
197 } 197 }
198 198
199 if (Es < 0.0) 199 if (Es < 0.0)
200 Es = 0.0; 200 Es = 0.0;
201 else if (Es > B * AP) 201 else if (Es > B * AP)
202 Es = B * AP; 202 Es = B * AP;
203 return Es; 203 return Es;
204 } 204 }
205 205
206 206
207 G4double G4NuclearAbrasionGeometry::GetExcitat 207 G4double G4NuclearAbrasionGeometry::GetExcitationEnergyOfTarget ()
208 { 208 {
209 // This member function declares a new G4Nuc 209 // This member function declares a new G4NuclearAbrasionGeometry object
210 // but with the projectile and target exchan 210 // but with the projectile and target exchanged to determine the values
211 // for F and P. Determination of the excess 211 // for F and P. Determination of the excess surface area and excitation
212 // energy is as above. 212 // energy is as above.
213 213
214 G4NuclearAbrasionGeometry* revAbrasionGeomet 214 G4NuclearAbrasionGeometry* revAbrasionGeometry =
215 new G4NuclearAbrasionGeometry(AT, AP, r); 215 new G4NuclearAbrasionGeometry(AT, AP, r);
216 G4double F1 = revAbrasionGeometry->F(); 216 G4double F1 = revAbrasionGeometry->F();
217 G4double P1 = revAbrasionGeometry->P(); 217 G4double P1 = revAbrasionGeometry->P();
218 G4double Es = 0.0; 218 G4double Es = 0.0;
219 219
220 Es = 0.95 * MeV * 4.0 * pi * rT*rT/fermi/fer 220 Es = 0.95 * MeV * 4.0 * pi * rT*rT/fermi/fermi *
221 (1.0+P1-G4Pow::GetInstance()->A23(1.0-F 221 (1.0+P1-G4Pow::GetInstance()->A23(1.0-F1));
222 222
223 // if (rP < rT && r < rT-rP) 223 // if (rP < rT && r < rT-rP)
224 if ((r-rT)/rP < rth) { 224 if ((r-rT)/rP < rth) {
225 G4double omega = 0.0; 225 G4double omega = 0.0;
226 if (AT < 12.0) omega = 1500.0; 226 if (AT < 12.0) omega = 1500.0;
227 else if (AT <= 16.0) omega = 1500.0 - 320. 227 else if (AT <= 16.0) omega = 1500.0 - 320.0*(AT-12.0);
228 Es *= 1.0 + F1*(5.0+omega*F1*F1); 228 Es *= 1.0 + F1*(5.0+omega*F1*F1);
229 } 229 }
230 230
231 if (Es < 0.0) 231 if (Es < 0.0)
232 Es = 0.0; 232 Es = 0.0;
233 else if (Es > B * AT) 233 else if (Es > B * AT)
234 Es = B * AT; 234 Es = B * AT;
235 235
236 delete revAbrasionGeometry; 236 delete revAbrasionGeometry;
237 237
238 return Es; 238 return Es;
239 } 239 }
240 240