ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ecpssrFormFactorMixsModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ecpssrFormFactorMixsModel.cc
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 // 01 Oct 2011 A.M., S.I. - 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 A. Taborda et al.
34 // EXRS2012 proceedings
35 // ---------------------------------------------------------------------------------------
36 
37 #include <fstream>
38 #include <iomanip>
39 
40 #include "globals.hh"
41 #include "G4ios.hh"
42 #include "G4SystemOfUnits.hh"
43 
44 #include "G4EMDataSet.hh"
45 #include "G4LinInterpolation.hh"
46 #include "G4Proton.hh"
47 #include "G4Alpha.hh"
48 
49 #include "G4ecpssrFormFactorMixsModel.hh"
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 
53 G4ecpssrFormFactorMixsModel::G4ecpssrFormFactorMixsModel()
54 {
55  interpolation = new G4LinInterpolation();
56 
57  for (G4int i=29; i<93; i++)
58  {
59  protonM1DataSetMap[i] = new G4EMDataSet(i,interpolation);
60  protonM1DataSetMap[i]->LoadData("pixe/ecpssr/proton/m1-i01m001c01-");
61 
62  protonM2DataSetMap[i] = new G4EMDataSet(i,interpolation);
63  protonM2DataSetMap[i]->LoadData("pixe/ecpssr/proton/m2-i01m001c01-");
64 
65  protonM3DataSetMap[i] = new G4EMDataSet(i,interpolation);
66  protonM3DataSetMap[i]->LoadData("pixe/ecpssr/proton/m3-i01m001c01-");
67 
68  protonM4DataSetMap[i] = new G4EMDataSet(i,interpolation);
69  protonM4DataSetMap[i]->LoadData("pixe/ecpssr/proton/m4-i01m001c01-");
70 
71  protonM5DataSetMap[i] = new G4EMDataSet(i,interpolation);
72  protonM5DataSetMap[i]->LoadData("pixe/ecpssr/proton/m5-i01m001c01-");
73  }
74 
75  protonMiXsVector.push_back(protonM1DataSetMap);
76  protonMiXsVector.push_back(protonM2DataSetMap);
77  protonMiXsVector.push_back(protonM3DataSetMap);
78  protonMiXsVector.push_back(protonM4DataSetMap);
79  protonMiXsVector.push_back(protonM5DataSetMap);
80 
81 
82  for (G4int i=29; i<93; i++)
83  {
84  alphaM1DataSetMap[i] = new G4EMDataSet(i,interpolation);
85  alphaM1DataSetMap[i]->LoadData("pixe/ecpssr/alpha/m1-i02m004c02-");
86 
87  alphaM2DataSetMap[i] = new G4EMDataSet(i,interpolation);
88  alphaM2DataSetMap[i]->LoadData("pixe/ecpssr/alpha/m2-i02m004c02-");
89 
90  alphaM3DataSetMap[i] = new G4EMDataSet(i,interpolation);
91  alphaM3DataSetMap[i]->LoadData("pixe/ecpssr/alpha/m3-i02m004c02-");
92 
93  alphaM4DataSetMap[i] = new G4EMDataSet(i,interpolation);
94  alphaM4DataSetMap[i]->LoadData("pixe/ecpssr/alpha/m4-i02m004c02-");
95 
96  alphaM5DataSetMap[i] = new G4EMDataSet(i,interpolation);
97  alphaM5DataSetMap[i]->LoadData("pixe/ecpssr/alpha/m5-i02m004c02-");
98 
99  }
100 
101  alphaMiXsVector.push_back(alphaM1DataSetMap);
102  alphaMiXsVector.push_back(alphaM2DataSetMap);
103  alphaMiXsVector.push_back(alphaM3DataSetMap);
104  alphaMiXsVector.push_back(alphaM4DataSetMap);
105  alphaMiXsVector.push_back(alphaM5DataSetMap);
106 
107 
108 
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112 
113 G4ecpssrFormFactorMixsModel::~G4ecpssrFormFactorMixsModel()
114 {
115  protonM1DataSetMap.clear();
116  alphaM1DataSetMap.clear();
117 
118  protonM2DataSetMap.clear();
119  alphaM2DataSetMap.clear();
120 
121  protonM3DataSetMap.clear();
122  alphaM3DataSetMap.clear();
123 
124  protonM4DataSetMap.clear();
125  alphaM4DataSetMap.clear();
126 
127  protonM5DataSetMap.clear();
128  alphaM5DataSetMap.clear();
129 
130  delete interpolation;
131 }
132 
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 
135 G4double G4ecpssrFormFactorMixsModel::CalculateMiCrossSection(G4int zTarget,G4double massIncident, G4double energyIncident, G4int mShellId)
136 {
137  G4Proton* aProton = G4Proton::Proton();
138  G4Alpha* aAlpha = G4Alpha::Alpha();
139  G4double sigma = 0;
140  G4int mShellIndex = mShellId -1;
141 
142  if (energyIncident > 0.1*MeV && energyIncident < 100*MeV && zTarget < 93 && zTarget > 28) {
143 
144  if (massIncident == aProton->GetPDGMass())
145  {
146  sigma = protonMiXsVector[mShellIndex][zTarget]->FindValue(energyIncident/MeV);
147  if (sigma !=0 && energyIncident > protonMiXsVector[mShellIndex][zTarget]->GetEnergies(0).back()*MeV) return 0.;
148  }
149  else if (massIncident == aAlpha->GetPDGMass())
150  {
151  sigma = alphaMiXsVector[mShellIndex][zTarget]->FindValue(energyIncident/MeV);
152  if (sigma !=0 && energyIncident > alphaMiXsVector[mShellIndex][zTarget]->GetEnergies(0).back()*MeV) return 0.;
153  }
154  else
155  {
156  sigma = 0.;
157  }
158  }
159 
160  // sigma is in internal units: it has been converted from
161  // the input file in barns bt the EmDataset
162  return sigma;
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166 
167 G4double G4ecpssrFormFactorMixsModel::CalculateM1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
168 {
169 
170  // mShellId
171  return CalculateMiCrossSection (zTarget, massIncident, energyIncident, 1);
172 
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176 
177 G4double G4ecpssrFormFactorMixsModel::CalculateM2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
178 {
179 
180  // mShellId
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  */
211 }
212 
213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
214 
215 G4double G4ecpssrFormFactorMixsModel::CalculateM3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
216 {
217 
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  */
248 }
249 
250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
251 
252 G4double G4ecpssrFormFactorMixsModel::CalculateM4CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
253 {
254 
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  */
283 }
284 
285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
286 
287 G4double G4ecpssrFormFactorMixsModel::CalculateM5CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
288 {
289 
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  */
318 }