| 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 // History: 26 // History:
27 // ----------- 27 // -----------
28 // 10 Nov 2021 S. Guatelli & S. Bakr, Wollo 28 // 10 Nov 2021 S. Guatelli & S. Bakr, Wollongong University - 1st implementation
29 // 29 //
30 // Class description 30 // Class description
31 // ---------------- 31 // ----------------
32 // Computation of K, L & M shell ECPSSR ionis 32 // Computation of K, L & M shell ECPSSR ionisation cross sections for protons and alphas
33 // Based on the work of 33 // Based on the work of
34 // - S. Bakr et al. (2021) NIM B, 507:11-19. 34 // - S. Bakr et al. (2021) NIM B, 507:11-19.
35 // - S. Bakr et al (2018), NIMB B, 436: 285-2 35 // - S. Bakr et al (2018), NIMB B, 436: 285-291.
36 // ------------------------------------------- 36 // ---------------------------------------------------------------------------------------
37 37
38 #include <fstream> 38 #include <fstream>
39 #include <iomanip> 39 #include <iomanip>
40 40
41 #include "globals.hh" 41 #include "globals.hh"
42 #include "G4ios.hh" 42 #include "G4ios.hh"
43 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh"
44 44
45 #include "G4EMDataSet.hh" 45 #include "G4EMDataSet.hh"
46 #include "G4LinInterpolation.hh" 46 #include "G4LinInterpolation.hh"
47 #include "G4Proton.hh" 47 #include "G4Proton.hh"
48 #include "G4Alpha.hh" 48 #include "G4Alpha.hh"
49 49
50 #include "G4ANSTOecpssrMixsModel.hh" 50 #include "G4ANSTOecpssrMixsModel.hh"
51 51
52 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 53
54 G4ANSTOecpssrMixsModel::G4ANSTOecpssrMixsModel 54 G4ANSTOecpssrMixsModel::G4ANSTOecpssrMixsModel()
55 { 55 {
56 G4cout << "Using ANSTO M Cross Sections! "<< 56 G4cout << "Using ANSTO M Cross Sections! "<< G4endl;
57 57
58 interpolation = new G4LinInterpolation(); 58 interpolation = new G4LinInterpolation();
59 59
60 for (G4int i=67; i<93; i++) 60 for (G4int i=67; i<93; i++)
61 { 61 {
62 protonM1DataSetMap[i] = new G4EMDataSet( 62 protonM1DataSetMap[i] = new G4EMDataSet(i,interpolation);
63 protonM1DataSetMap[i]->LoadData("pixe_AN 63 protonM1DataSetMap[i]->LoadData("pixe_ANSTO/proton/m1-");
64 64
65 protonM2DataSetMap[i] = new G4EMDataSet( 65 protonM2DataSetMap[i] = new G4EMDataSet(i,interpolation);
66 protonM2DataSetMap[i]->LoadData("pixe_AN 66 protonM2DataSetMap[i]->LoadData("pixe_ANSTO/proton/m2-");
67 67
68 protonM3DataSetMap[i] = new G4EMDataSet( 68 protonM3DataSetMap[i] = new G4EMDataSet(i,interpolation);
69 protonM3DataSetMap[i]->LoadData("pixe_AN 69 protonM3DataSetMap[i]->LoadData("pixe_ANSTO/proton/m3-");
70 70
71 protonM4DataSetMap[i] = new G4EMDataSet( 71 protonM4DataSetMap[i] = new G4EMDataSet(i,interpolation);
72 protonM4DataSetMap[i]->LoadData("pixe_AN 72 protonM4DataSetMap[i]->LoadData("pixe_ANSTO/proton/m4-");
73 73
74 protonM5DataSetMap[i] = new G4EMDataSet( 74 protonM5DataSetMap[i] = new G4EMDataSet(i,interpolation);
75 protonM5DataSetMap[i]->LoadData("pixe_AN 75 protonM5DataSetMap[i]->LoadData("pixe_ANSTO/proton/m5-");
76 } 76 }
77 77
78 protonMiXsVector.push_back(protonM1DataSetMa 78 protonMiXsVector.push_back(protonM1DataSetMap);
79 protonMiXsVector.push_back(protonM2DataSetMa 79 protonMiXsVector.push_back(protonM2DataSetMap);
80 protonMiXsVector.push_back(protonM3DataSetMa 80 protonMiXsVector.push_back(protonM3DataSetMap);
81 protonMiXsVector.push_back(protonM4DataSetMa 81 protonMiXsVector.push_back(protonM4DataSetMap);
82 protonMiXsVector.push_back(protonM5DataSetMa 82 protonMiXsVector.push_back(protonM5DataSetMap);
83 83
84 84
85 for (G4int i=67; i<93; i++) 85 for (G4int i=67; i<93; i++)
86 { 86 {
87 alphaM1DataSetMap[i] = new G4EMDataSet(i 87 alphaM1DataSetMap[i] = new G4EMDataSet(i,interpolation);
88 alphaM1DataSetMap[i]->LoadData("pixe_ANS 88 alphaM1DataSetMap[i]->LoadData("pixe_ANSTO/alpha/m1-");
89 89
90 alphaM2DataSetMap[i] = new G4EMDataSet(i 90 alphaM2DataSetMap[i] = new G4EMDataSet(i,interpolation);
91 alphaM2DataSetMap[i]->LoadData("pixe_ANS 91 alphaM2DataSetMap[i]->LoadData("pixe_ANSTO/alpha/m2-");
92 92
93 alphaM3DataSetMap[i] = new G4EMDataSet(i 93 alphaM3DataSetMap[i] = new G4EMDataSet(i,interpolation);
94 alphaM3DataSetMap[i]->LoadData("pixe_ANS 94 alphaM3DataSetMap[i]->LoadData("pixe_ANSTO/alpha/m3-");
95 95
96 alphaM4DataSetMap[i] = new G4EMDataSet(i 96 alphaM4DataSetMap[i] = new G4EMDataSet(i,interpolation);
97 alphaM4DataSetMap[i]->LoadData("pixe_ANS 97 alphaM4DataSetMap[i]->LoadData("pixe_ANSTO/alpha/m4-");
98 98
99 alphaM5DataSetMap[i] = new G4EMDataSet(i 99 alphaM5DataSetMap[i] = new G4EMDataSet(i,interpolation);
100 alphaM5DataSetMap[i]->LoadData("pixe_ANS 100 alphaM5DataSetMap[i]->LoadData("pixe_ANSTO/alpha/m5-");
101 } 101 }
102 102
103 alphaMiXsVector.push_back(alphaM1DataSetMap) 103 alphaMiXsVector.push_back(alphaM1DataSetMap);
104 alphaMiXsVector.push_back(alphaM2DataSetMap) 104 alphaMiXsVector.push_back(alphaM2DataSetMap);
105 alphaMiXsVector.push_back(alphaM3DataSetMap) 105 alphaMiXsVector.push_back(alphaM3DataSetMap);
106 alphaMiXsVector.push_back(alphaM4DataSetMap) 106 alphaMiXsVector.push_back(alphaM4DataSetMap);
107 alphaMiXsVector.push_back(alphaM5DataSetMap) 107 alphaMiXsVector.push_back(alphaM5DataSetMap);
108 } 108 }
109 109
110 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 111
112 G4ANSTOecpssrMixsModel::~G4ANSTOecpssrMixsMode 112 G4ANSTOecpssrMixsModel::~G4ANSTOecpssrMixsModel()
113 { 113 {
114 protonM1DataSetMap.clear(); 114 protonM1DataSetMap.clear();
115 alphaM1DataSetMap.clear(); 115 alphaM1DataSetMap.clear();
116 116
117 protonM2DataSetMap.clear(); 117 protonM2DataSetMap.clear();
118 alphaM2DataSetMap.clear(); 118 alphaM2DataSetMap.clear();
119 119
120 protonM3DataSetMap.clear(); 120 protonM3DataSetMap.clear();
121 alphaM3DataSetMap.clear(); 121 alphaM3DataSetMap.clear();
122 122
123 protonM4DataSetMap.clear(); 123 protonM4DataSetMap.clear();
124 alphaM4DataSetMap.clear(); 124 alphaM4DataSetMap.clear();
125 125
126 protonM5DataSetMap.clear(); 126 protonM5DataSetMap.clear();
127 alphaM5DataSetMap.clear(); 127 alphaM5DataSetMap.clear();
128 128
129 delete interpolation; 129 delete interpolation;
130 } 130 }
131 131
132 //....oooOO0OOooo........oooOO0OOooo........oo 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133 133
134 G4double G4ANSTOecpssrMixsModel::CalculateMiCr 134 G4double G4ANSTOecpssrMixsModel::CalculateMiCrossSection(G4int zTarget,G4double massIncident, G4double energyIncident, G4int mShellId)
135 { 135 {
136 G4Proton* aProton = G4Proton::Proton(); 136 G4Proton* aProton = G4Proton::Proton();
137 G4Alpha* aAlpha = G4Alpha::Alpha(); 137 G4Alpha* aAlpha = G4Alpha::Alpha();
138 G4double sigma = 0; 138 G4double sigma = 0;
139 G4int mShellIndex = mShellId -1; 139 G4int mShellIndex = mShellId -1;
140 140
141 if (massIncident == aProton->GetPDGMass()) 141 if (massIncident == aProton->GetPDGMass())
142 { 142 {
143 if (energyIncident > 0.2*MeV && energyInc 143 if (energyIncident > 0.2*MeV && energyIncident < 5.*MeV && zTarget < 93 && zTarget > 66) {
144 144
145 sigma = protonMiXsVector[mShellIndex][zT 145 sigma = protonMiXsVector[mShellIndex][zTarget]->FindValue(energyIncident/MeV);
146 if (sigma !=0 && energyIncident > prot 146 if (sigma !=0 && energyIncident > protonMiXsVector[mShellIndex][zTarget]->GetEnergies(0).back()*MeV) return 0.;
147 } 147 }
148 } 148 }
149 149
150 else if (massIncident == aAlpha->GetPDGMas 150 else if (massIncident == aAlpha->GetPDGMass())
151 { 151 {
152 if (energyIncident > 0.2*MeV && en 152 if (energyIncident > 0.2*MeV && energyIncident < 10.*MeV && zTarget < 93 && zTarget > 66) {
153 153
154 sigma = alphaMiXsVector[mShellInde 154 sigma = alphaMiXsVector[mShellIndex][zTarget]->FindValue(energyIncident/MeV);
155 if (sigma !=0 && energyIncident > 155 if (sigma !=0 && energyIncident > alphaMiXsVector[mShellIndex][zTarget]->GetEnergies(0).back()*MeV) return 0.;
156 } 156 }
157 } 157 }
158 158
159 else 159 else
160 { 160 {
161 sigma = 0.; 161 sigma = 0.;
162 } 162 }
163 163
164 164
165 // sigma is in internal units: it has been c 165 // sigma is in internal units: it has been converted from
166 // the input file in barns bt the EmDataset 166 // the input file in barns bt the EmDataset
167 return sigma; 167 return sigma;
168 } 168 }
169 169
170 //....oooOO0OOooo........oooOO0OOooo........oo 170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171 171
172 G4double G4ANSTOecpssrMixsModel::CalculateM1Cr 172 G4double G4ANSTOecpssrMixsModel::CalculateM1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
173 { 173 {
174 174
175 // mShellId 175 // mShellId
176 return CalculateMiCrossSection (zTarget, ma 176 return CalculateMiCrossSection (zTarget, massIncident, energyIncident, 1);
177 177
178 } 178 }
179 179
180 //....oooOO0OOooo........oooOO0OOooo........oo 180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
181 181
182 G4double G4ANSTOecpssrMixsModel::CalculateM2Cr 182 G4double G4ANSTOecpssrMixsModel::CalculateM2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
183 { 183 {
184 184
185 // mShellId 185 // mShellId
186 return CalculateMiCrossSection (zTarget, ma 186 return CalculateMiCrossSection (zTarget, massIncident, energyIncident, 2);
187 187
188 /* 188 /*
189 189
190 G4Proton* aProton = G4Proton::Proton(); 190 G4Proton* aProton = G4Proton::Proton();
191 G4Alpha* aAlpha = G4Alpha::Alpha(); 191 G4Alpha* aAlpha = G4Alpha::Alpha();
192 G4double sigma = 0; 192 G4double sigma = 0;
193 193
194 if (energyIncident > 0.1*MeV && energyIncide 194 if (energyIncident > 0.1*MeV && energyIncident < 10*MeV && zTarget < 93 && zTarget > 61) {
195 195
196 if (massIncident == aProton->GetPDGMass()) 196 if (massIncident == aProton->GetPDGMass())
197 { 197 {
198 sigma = protonM2DataSetMap[zTarget]->FindVal 198 sigma = protonM2DataSetMap[zTarget]->FindValue(energyIncident/MeV);
199 if (sigma !=0 && energyIncident > prot 199 if (sigma !=0 && energyIncident > protonM2DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
200 } 200 }
201 else if (massIncident == aAlpha->GetPDGMas 201 else if (massIncident == aAlpha->GetPDGMass())
202 { 202 {
203 sigma = alphaM2DataSetMap[zTarget]->Fi 203 sigma = alphaM2DataSetMap[zTarget]->FindValue(energyIncident/MeV);
204 if (sigma !=0 && energyIncident > alph 204 if (sigma !=0 && energyIncident > alphaM2DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
205 } 205 }
206 else 206 else
207 { 207 {
208 sigma = 0.; 208 sigma = 0.;
209 } 209 }
210 } 210 }
211 211
212 // sigma is in internal units: it has been c 212 // sigma is in internal units: it has been converted from
213 // the input file in barns bt the EmDataset 213 // the input file in barns bt the EmDataset
214 return sigma; 214 return sigma;
215 */ 215 */
216 } 216 }
217 217
218 //....oooOO0OOooo........oooOO0OOooo........oo 218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219 219
220 G4double G4ANSTOecpssrMixsModel::CalculateM3Cr 220 G4double G4ANSTOecpssrMixsModel::CalculateM3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
221 { 221 {
222 222
223 return CalculateMiCrossSection (zTarget, ma 223 return CalculateMiCrossSection (zTarget, massIncident, energyIncident, 3);
224 /* 224 /*
225 225
226 226
227 G4Proton* aProton = G4Proton::Proton(); 227 G4Proton* aProton = G4Proton::Proton();
228 G4Alpha* aAlpha = G4Alpha::Alpha(); 228 G4Alpha* aAlpha = G4Alpha::Alpha();
229 G4double sigma = 0; 229 G4double sigma = 0;
230 230
231 if (energyIncident > 0.1*MeV && energyIncide 231 if (energyIncident > 0.1*MeV && energyIncident < 10*MeV && zTarget < 93 && zTarget > 61) {
232 232
233 if (massIncident == aProton->GetPDGMass()) 233 if (massIncident == aProton->GetPDGMass())
234 { 234 {
235 sigma = protonM3DataSetMap[zTarget]->FindVal 235 sigma = protonM3DataSetMap[zTarget]->FindValue(energyIncident/MeV);
236 if (sigma !=0 && energyIncident > prot 236 if (sigma !=0 && energyIncident > protonM3DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
237 } 237 }
238 else if (massIncident == aAlpha->GetPDGMas 238 else if (massIncident == aAlpha->GetPDGMass())
239 { 239 {
240 sigma = alphaM3DataSetMap[zTarget]->Fi 240 sigma = alphaM3DataSetMap[zTarget]->FindValue(energyIncident/MeV);
241 if (sigma !=0 && energyIncident > alph 241 if (sigma !=0 && energyIncident > alphaM3DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
242 } 242 }
243 else 243 else
244 { 244 {
245 sigma = 0.; 245 sigma = 0.;
246 } 246 }
247 } 247 }
248 248
249 // sigma is in internal units: it has been c 249 // sigma is in internal units: it has been converted from
250 // the input file in barns bt the EmDataset 250 // the input file in barns bt the EmDataset
251 return sigma; 251 return sigma;
252 */ 252 */
253 } 253 }
254 254
255 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
256 256
257 G4double G4ANSTOecpssrMixsModel::CalculateM4Cr 257 G4double G4ANSTOecpssrMixsModel::CalculateM4CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
258 { 258 {
259 259
260 return CalculateMiCrossSection (zTarget, ma 260 return CalculateMiCrossSection (zTarget, massIncident, energyIncident, 4);
261 /* 261 /*
262 G4Proton* aProton = G4Proton::Proton(); 262 G4Proton* aProton = G4Proton::Proton();
263 G4Alpha* aAlpha = G4Alpha::Alpha(); 263 G4Alpha* aAlpha = G4Alpha::Alpha();
264 G4double sigma = 0; 264 G4double sigma = 0;
265 265
266 if (energyIncident > 0.1*MeV && energyIncide 266 if (energyIncident > 0.1*MeV && energyIncident < 10*MeV && zTarget < 93 && zTarget > 61) {
267 267
268 if (massIncident == aProton->GetPDGMass()) 268 if (massIncident == aProton->GetPDGMass())
269 { 269 {
270 sigma = protonM3DataSetMap[zTarget]->FindVal 270 sigma = protonM3DataSetMap[zTarget]->FindValue(energyIncident/MeV);
271 if (sigma !=0 && energyIncident > prot 271 if (sigma !=0 && energyIncident > protonM3DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
272 } 272 }
273 else if (massIncident == aAlpha->GetPDGMas 273 else if (massIncident == aAlpha->GetPDGMass())
274 { 274 {
275 sigma = alphaM3DataSetMap[zTarget]->Fi 275 sigma = alphaM3DataSetMap[zTarget]->FindValue(energyIncident/MeV);
276 if (sigma !=0 && energyIncident > alph 276 if (sigma !=0 && energyIncident > alphaM3DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
277 } 277 }
278 else 278 else
279 { 279 {
280 sigma = 0.; 280 sigma = 0.;
281 } 281 }
282 } 282 }
283 283
284 // sigma is in internal units: it has been c 284 // sigma is in internal units: it has been converted from
285 // the input file in barns bt the EmDataset 285 // the input file in barns bt the EmDataset
286 return sigma; 286 return sigma;
287 */ 287 */
288 } 288 }
289 289
290 //....oooOO0OOooo........oooOO0OOooo........oo 290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291 291
292 G4double G4ANSTOecpssrMixsModel::CalculateM5Cr 292 G4double G4ANSTOecpssrMixsModel::CalculateM5CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
293 { 293 {
294 294
295 return CalculateMiCrossSection (zTarget, ma 295 return CalculateMiCrossSection (zTarget, massIncident, energyIncident, 5);
296 /* 296 /*
297 G4Proton* aProton = G4Proton::Proton(); 297 G4Proton* aProton = G4Proton::Proton();
298 G4Alpha* aAlpha = G4Alpha::Alpha(); 298 G4Alpha* aAlpha = G4Alpha::Alpha();
299 G4double sigma = 0; 299 G4double sigma = 0;
300 300
301 if (energyIncident > 0.1*MeV && energyIncide 301 if (energyIncident > 0.1*MeV && energyIncident < 10*MeV && zTarget < 93 && zTarget > 61) {
302 302
303 if (massIncident == aProton->GetPDGMass()) 303 if (massIncident == aProton->GetPDGMass())
304 { 304 {
305 sigma = protonM3DataSetMap[zTarget]->FindVal 305 sigma = protonM3DataSetMap[zTarget]->FindValue(energyIncident/MeV);
306 if (sigma !=0 && energyIncident > prot 306 if (sigma !=0 && energyIncident > protonM3DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
307 } 307 }
308 else if (massIncident == aAlpha->GetPDGMas 308 else if (massIncident == aAlpha->GetPDGMass())
309 { 309 {
310 sigma = alphaM3DataSetMap[zTarget]->Fi 310 sigma = alphaM3DataSetMap[zTarget]->FindValue(energyIncident/MeV);
311 if (sigma !=0 && energyIncident > alph 311 if (sigma !=0 && energyIncident > alphaM3DataSetMap[zTarget]->GetEnergies(0).back()*MeV) return 0.;
312 } 312 }
313 else 313 else
314 { 314 {
315 sigma = 0.; 315 sigma = 0.;
316 } 316 }
317 } 317 }
318 318
319 // sigma is in internal units: it has been c 319 // sigma is in internal units: it has been converted from
320 // the input file in barns bt the EmDataset 320 // the input file in barns bt the EmDataset
321 return sigma; 321 return sigma;
322 */ 322 */
323 } 323 }
324 324