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