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