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