ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
XrayFluoHPGeDetectorType.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file XrayFluoHPGeDetectorType.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 //
27 //
28 // Author: Alfonso Mantero (Alfonso.Mantero@ge.infn.it)
29 //
30 // History:
31 // -----------
32 // 16 Jul 2003 Alfonso Mantero Created
33 //
34 // -------------------------------------------------------------------
35 
36 #include <fstream>
37 #include <sstream>
38 
40 #include "XrayFluoDataSet.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4UnitsTable.hh"
43 #include "G4DataVector.hh"
44 #include "G4LogLogInterpolation.hh"
45 #include "G4ios.hh"
46 #include "Randomize.hh"
47 
48 
50  detectorMaterial("HPGe"),efficiencySet(0)
51 {
52  LoadResponseData("response");
53 
54  LoadEfficiencyData("efficienza");
55 
56 
57 }
59 {
60  std::map<G4int,G4DataVector*,std::less<G4int> >::iterator pos;
61 
62  for (pos = energyMap.begin(); pos != energyMap.end(); pos++)
63  {
64  G4DataVector* dataSet = (*pos).second;
65  delete dataSet;
66  dataSet = 0;
67  }
68  for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
69  {
70  G4DataVector* dataSet = (*pos).second;
71  delete dataSet;
72  dataSet = 0;
73  }
74 
75  delete interpolation4;
76 
77 }
79 {
80  return detectorMaterial;
81 }
82 
84 
86 
87 {
88  if (instance == 0)
89  {
91 
92  }
93  return instance;
94 }
95 
97 {
98 
99 
100  G4double eMin = 1* keV;
101  G4double eMax = 10*keV;
102  G4double value = 0.;
103  G4double efficiency = 1.;
104 
105  const XrayFluoDataSet* dataSet = efficiencySet;
106  G4int id = 0;
107 
108  G4double random = G4UniformRand();
109 
110  if (energy>=eMin && energy <=eMax)
111  {
112  G4double infEnergy = (G4int)(energy/keV)* keV;
113 
114  G4double supEnergy = ((G4int)(energy/keV) + 1)*keV;
115 
116  G4double infData = GetInfData(energy, random, 0);// 0 is not used
117 
118  G4double supData = GetSupData(energy,random, 0); // 0 is not used
119 
120  value = (std::log10(infData)*std::log10(supEnergy/energy) +
121  std::log10(supData)*std::log10(energy/infEnergy)) /
122  std::log10(supEnergy/infEnergy);
123  value = std::pow(10,value);
124  }
125  else if (energy<eMin)
126  {
127  G4double infEnergy = eMin;
128  G4double supEnergy = eMin/keV +1*keV;
129 
130  G4double infData = GetInfData(eMin, random, 0);
131  G4double supData = GetSupData(eMin,random, 0);
132  value = (std::log10(infData)*std::log10(supEnergy/eMin) +
133  std::log10(supData)*std::log10(eMin/infEnergy)) /
134  std::log10(supEnergy/infEnergy);
135  value = std::pow(10,value);
136  value = value-eMin+ energy;
137 
138 
139  }
140  else if (energy>eMax)
141  {
142  G4double infEnergy = eMax/keV - 1. *keV;
143  G4double supEnergy = eMax;
144 
145  G4double infData = GetInfData(eMax, random, 0);
146  G4double supData = GetSupData(eMax,random, 0);
147  value = (std::log10(infData)*std::log10(supEnergy/eMax) +
148  std::log10(supData)*std::log10(eMax/infEnergy)) /
149  std::log10(supEnergy/infEnergy);
150  value = std::pow(10,value);
151  value = value+energy- eMax;
152  }
153  G4double RandomNum = G4UniformRand();
154 
155  efficiency = dataSet->FindValue(value,id);
156  if ( RandomNum>efficiency )
157  {
158  value = 0.;
159  }
160 
161  return value;
162 }
164 {
165  G4double value = 0.;
166  G4int zMin = 1;
167  G4int zMax = 10;
168 
169  G4int Z = ((G4int)(energy/keV));
170 
171  if (Z<zMin) {Z=zMin;}
172  if (Z>zMax) {Z=zMax;}
173 
174  if (Z >= zMin && Z <= zMax)
175  {
176  std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
177  pos = energyMap.find(Z);
178  std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator posData;
179  posData = dataMap.find(Z);
180  if (pos!= energyMap.end())
181  {
182  G4DataVector energySet = *((*pos).second);
183  G4DataVector dataSet = *((*posData).second);
184  G4int nData = energySet.size();
185 
186  G4double partSum = 0;
187  G4int index = 0;
188 
189  while (random> partSum)
190  {
191  partSum += dataSet[index];
192  index++;
193  }
194 
195 
196  if (index >= 0 && index < nData)
197  {
198  value = energySet[index];
199 
200  }
201 
202  }
203  }
204  return value;
205 
206 }
208 {
209  G4double value = 0.;
210  G4int zMin = 1;
211  G4int zMax = 10;
212  G4int Z = ((G4int)(energy/keV)+1);
213 
214  if (Z<zMin) {Z=zMin;}
215  if (Z>zMax) {Z=zMax;}
216  if (Z >= zMin && Z <= zMax)
217  {
218  std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
219  pos = energyMap.find(Z);
220  std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator posData;
221  posData = dataMap.find(Z);
222  if (pos!= energyMap.end())
223  {
224  G4DataVector energySet = *((*pos).second);
225  G4DataVector dataSet = *((*posData).second);
226  G4int nData = energySet.size();
227  G4double partSum = 0;
228  G4int index = 0;
229 
230  while (random> partSum)
231  {
232  partSum += dataSet[index];
233  index++;
234  }
235 
236 
237  if (index >= 0 && index < nData)
238  {
239  value = energySet[index];
240  }
241  }
242  }
243  return value;
244 
245 }
247 {
248  std::ostringstream ost;
249 
250  ost << fileName<<".dat";
251 
252  G4String name = ost.str();
253 
254  char* path = std::getenv("XRAYDATA");
255 
256  G4String pathString(path);
257  G4String dirFile = pathString + "/" + name;
258  std::ifstream file(dirFile);
259  std::filebuf* lsdp = file.rdbuf();
260 
261  if (! (lsdp->is_open()) )
262  {
264  execp << "XrayFluoHPGeDetectorType - data file: " + dirFile + " not found";
265  G4Exception("XrayFluoHPGeDetectorType::LoadResponseData()","example-xray_fluorescence02",
266  FatalException, execp);
267  }
268  G4double a = 0;
269  G4int k = 1;
270  G4int q = 0;
271 
272  G4int Z = 1;
273  G4DataVector* energies = new G4DataVector;
275 
276  do
277  {
278  file >> a;
279  G4int nColumns = 2;
280  if (a == -1)
281  {
282  if (q == 0)
283  {
284  // End of a data set
285  energyMap[Z] = energies;
286  dataMap[Z] = data;
287  // Start of new shell data set
288  energies = new G4DataVector;
289  data = new G4DataVector;
290  Z++;
291  }
292  q++;
293  if (q == nColumns)
294  {
295  q = 0;
296  }
297  }
298  else if (a == -2)
299  {
300  // End of file; delete the empty vectors
301  //created when encountering the last -1 -1 row
302  delete energies;
303  delete data;
304 
305  }
306  else
307  {
308  // 1st column is energy
309  if(k%nColumns != 0)
310  {
311  G4double e = a * keV;
312  energies->push_back(e);
313  k++;
314  }
315  else if (k%nColumns == 0)
316  {
317  // 2nd column is data
318 
319  data->push_back(a);
320  k = 1;
321  }
322  }
323  } while (a != -2); // end of file
324  file.close();
325 }
326 
328 {
329  /*
330  char* path = std::getenv("XRAYDATA");
331  G4String dirFile;
332  if (path) {
333  G4String pathString(path);
334  dirFile = pathString + "/" + fileName;
335  }
336  else{
337  path = std::getenv("PWD");
338  G4String pathString(path);
339  dirFile = pathString + "/" + fileName;
340  }
341  */
343  efficiencySet = new XrayFluoDataSet(1,fileName,interpolation4,keV,1);
344 }
345