Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RDBremsstrahlungParameters.cc
Go to the documentation of this file.
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 // $Id$
28 // GEANT4 tag $Name: geant4-09-01-ref-00 $
29 //
30 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
31 // V.Ivanchenko (Vladimir.Ivantchenko@cern.ch)
32 //
33 // History:
34 // -----------
35 // 31 Jul 2001 MGP Created
36 // 12.09.01 V.Ivanchenko Add activeZ and paramA
37 // 25.09.01 V.Ivanchenko Add parameter C and change interface to B
38 // 29.11.01 V.Ivanchenko Update parametrisation
39 // 18.11.02 V.Ivanchenko Fix problem of load
40 // 21.02.03 V.Ivanchenko Number of parameters is defined in the constructor
41 // 28.02.03 V.Ivanchenko Filename is defined in the constructor
42 //
43 // -------------------------------------------------------------------
44 
46 #include "G4RDVEMDataSet.hh"
47 #include "G4RDEMDataSet.hh"
49 #include "G4Material.hh"
50 #include <fstream>
51 
52 
54  size_t num, G4int minZ, G4int maxZ)
55  : zMin(minZ),
56  zMax(maxZ),
57  length(num)
58 {
59  LoadData(name);
60 }
61 
62 
64 {
65  // Reset the map of data sets: remove the data sets from the map
66  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
67 
68  for (pos = param.begin(); pos != param.end(); ++pos)
69  {
70  G4RDVEMDataSet* dataSet = (*pos).second;
71  delete dataSet;
72  }
73 
74  activeZ.clear();
75  paramC.clear();
76 }
77 
78 
80  G4int Z,
81  G4double energy) const
82 {
83  G4double value = 0.;
84  G4int id = Z*length + parameterIndex;
85  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
86 
87  pos = param.find(id);
88  if (pos!= param.end()) {
89 
90  G4RDVEMDataSet* dataSet = (*pos).second;
91  const G4DataVector ener = dataSet->GetEnergies(0);
92  G4double ee = std::max(ener.front(),std::min(ener.back(),energy));
93  value = dataSet->FindValue(ee);
94 
95  } else {
96  G4cout << "WARNING: G4RDBremsstrahlungParameters::FindValue "
97  << "did not find ID = "
98  << id << G4endl;
99  }
100 
101  return value;
102 }
103 
104 void G4RDBremsstrahlungParameters::LoadData(const G4String& name)
105 {
106  // Build the complete string identifying the file with the data set
107 
108  // define active elements
109 
110  const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
111  if (materialTable == 0)
112  G4Exception("G4RDBremsstrahlungParameters::LoadData()",
113  "DataNotFound", FatalException, "No MaterialTable found!");
114 
116 
117  G4double x = 1.e-9;
118  for (G4int mm=0; mm<100; mm++) {
119  paramC.push_back(x);
120  }
121 
122  for (G4int m=0; m<nMaterials; m++) {
123 
124  const G4Material* material= (*materialTable)[m];
125  const G4ElementVector* elementVector = material->GetElementVector();
126  const G4int nElements = material->GetNumberOfElements();
127 
128  for (G4int iEl=0; iEl<nElements; iEl++) {
129  G4Element* element = (*elementVector)[iEl];
130  G4double Z = element->GetZ();
131  G4int iz = (G4int)Z;
132  if(iz < 100)
133  paramC[iz] = 0.217635e-33*(material->GetTotNbOfElectPerVolume());
134  if (!(activeZ.contains(Z))) {
135  activeZ.push_back(Z);
136  }
137  }
138  }
139 
140  // Read parameters
141 
142  char* path = getenv("G4LEDATA");
143  if (path == 0)
144  {
145  G4String excep("G4LEDATA environment variable not set!");
146  G4Exception("G4RDBremsstrahlungParameters::LoadData()",
147  "InvalidSetup", FatalException, excep);
148  }
149 
150  G4String pathString_a(path);
151  G4String name_a = pathString_a + name;
152  std::ifstream file_a(name_a);
153  std::filebuf* lsdp_a = file_a.rdbuf();
154 
155  if (! (lsdp_a->is_open()) )
156  {
157  G4String stringConversion2("Cannot open file ");
158  G4String excep = stringConversion2 + name_a;
159  G4Exception("G4RDBremsstrahlungParameters::LoadData()",
160  "FileNotFound", FatalException, excep);
161  }
162 
163  // The file is organized into two columns:
164  // 1st column is the energy
165  // 2nd column is the corresponding value
166  // The file terminates with the pattern: -1 -1
167  // -2 -2
168 
169  G4double ener = 0.0;
170  G4double sum = 0.0;
171  G4int z = 0;
172 
173  std::vector<G4DataVector*> a;
174  for (size_t j=0; j<length; j++) {
175  G4DataVector* aa = new G4DataVector();
176  a.push_back(aa);
177  }
178  G4DataVector e;
179  e.clear();
180 
181  do {
182  file_a >> ener >> sum;
183 
184  // End of file
185  if (ener == (G4double)(-2)) {
186  break;
187 
188  // End of next element
189  } else if (ener == (G4double)(-1)) {
190 
191  z++;
192  G4double Z = (G4double)z;
193 
194  // fill map if Z is used
195  if (activeZ.contains(Z)) {
196 
197  for (size_t k=0; k<length; k++) {
198 
199  G4int id = z*length + k;
201  G4DataVector* eVector = new G4DataVector;
202  size_t eSize = e.size();
203  for (size_t s=0; s<eSize; s++) {
204  eVector->push_back(e[s]);
205  }
206  G4RDVEMDataSet* set = new G4RDEMDataSet(id,eVector,a[k],inter,1.,1.);
207  param[id] = set;
208  }
209  a.clear();
210  for (size_t j=0; j<length; j++) {
211  G4DataVector* aa = new G4DataVector();
212  a.push_back(aa);
213  }
214  } else {
215  for (size_t j=0; j<length; j++) {
216  a[j]->clear();
217  }
218  }
219  e.clear();
220 
221  } else {
222 
223  if(ener > 1000.) ener = 1000.;
224  e.push_back(ener);
225  a[length-1]->push_back(sum);
226 
227  for (size_t j=0; j<length-1; j++) {
228  G4double qRead;
229  file_a >> qRead;
230  a[j]->push_back(qRead);
231  }
232 
233  }
234  } while (ener != (G4double)(-2));
235 
236  file_a.close();
237 
238 }
239 
240 
242 {
243  G4int n = paramC.size();
244  if (id < 0 || id >= n)
245  {
246  G4String stringConversion1("Wrong id = ");
247  G4String stringConversion2(id);
248  G4String ex = stringConversion1 + stringConversion2;
249  G4Exception("G4RDBremsstrahlungParameters::ParameterC()",
250  "InvalidSetup", FatalException, ex);
251  }
252 
253  return paramC[id];
254 }
255 
256 
258 {
259 
260  G4cout << G4endl;
261  G4cout << "===== G4RDBremsstrahlungParameters =====" << G4endl;
262  G4cout << G4endl;
263  G4cout << "===== Parameters =====" << G4endl;
264  G4cout << G4endl;
265 
266  size_t nZ = activeZ.size();
267  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
268 
269  for (size_t j=0; j<nZ; j++) {
270  G4int Z = (G4int)activeZ[j];
271 
272  for (size_t i=0; i<length; i++) {
273 
274  pos = param.find(Z*length + i);
275  if (pos!= param.end()) {
276 
277  G4cout << "===== Z= " << Z
278  << " parameter[" << i << "] ====="
279  << G4endl;
280  G4RDVEMDataSet* dataSet = (*pos).second;
281  dataSet->PrintData();
282  }
283  }
284  }
285 
286  G4cout << "==========================================" << G4endl;
287 }
288 
const XML_Char * name
Definition: expat.h:151
static constexpr double mm
Definition: G4SIunits.hh:115
G4double GetTotNbOfElectPerVolume() const
Definition: G4Material.hh:212
std::vector< G4Element * > G4ElementVector
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
tuple x
Definition: test.py:50
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
string material
Definition: eplot.py:19
const XML_Char * s
Definition: expat.h:262
G4GLOB_DLL std::ostream G4cout
static constexpr double m
Definition: G4SIunits.hh:129
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4int n
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
virtual void PrintData(void) const =0
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
G4double ParameterC(G4int index) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double Parameter(G4int parameterIndex, G4int Z, G4double energy) const
virtual const G4DataVector & GetEnergies(G4int componentId) const =0
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
tuple z
Definition: test.py:28
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
G4bool contains(const G4double &) const
G4RDBremsstrahlungParameters(const G4String &name, size_t num, G4int minZ=1, G4int maxZ=99)
static const G4double pos