ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BremsstrahlungParameters.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BremsstrahlungParameters.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: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // V.Ivanchenko (Vladimir.Ivantchenko@cern.ch)
30 //
31 // History:
32 // -----------
33 // 31 Jul 2001 MGP Created
34 // 12.09.01 V.Ivanchenko Add activeZ and paramA
35 // 25.09.01 V.Ivanchenko Add parameter C and change interface to B
36 // 29.11.01 V.Ivanchenko Update parametrisation
37 // 18.11.02 V.Ivanchenko Fix problem of load
38 // 21.02.03 V.Ivanchenko Number of parameters is defined in the constructor
39 // 28.02.03 V.Ivanchenko Filename is defined in the constructor
40 // 03.12.10 V.Ivanchenko Fixed memory leak in LoadData
41 //
42 // -------------------------------------------------------------------
43 
44 #include <fstream>
45 
47 #include "G4PhysicalConstants.hh"
48 #include "G4VEMDataSet.hh"
49 #include "G4EMDataSet.hh"
50 #include "G4LogLogInterpolation.hh"
51 #include "G4Material.hh"
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,G4VEMDataSet*,std::less<G4int> >::iterator pos;
67 
68  for (pos = param.begin(); pos != param.end(); ++pos)
69  {
70  G4VEMDataSet* 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,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
86 
87  pos = param.find(id);
88  if (pos!= param.end()) {
89 
90  G4VEMDataSet* 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: G4BremsstrahlungParameters::FindValue "
97  << "did not find ID = "
98  << id << G4endl;
99  }
100 
101  return value;
102 }
103 
105 {
106  const G4double mConst =
108 
109  // Build the complete string identifying the file with the data set
110  // define active elements
111 
112  const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
113  if (materialTable == 0)
114  G4Exception("G4BremsstrahlungParameters::LoadData",
115  "em1001",FatalException,"Unable to find MaterialTable");
116 
118 
119  G4double x = 1.e-9;
120  for (G4int mmLocal=0; mmLocal<100; mmLocal++) {
121  paramC.push_back(x);
122  }
123 
124  for (G4int mLocal=0; mLocal<nMaterials; mLocal++) {
125 
126  const G4Material* material= (*materialTable)[mLocal];
127  const G4ElementVector* elementVector = material->GetElementVector();
128  const G4int nElements = material->GetNumberOfElements();
129 
130  for (G4int iEl=0; iEl<nElements; iEl++) {
131  G4Element* element = (*elementVector)[iEl];
132  G4double Z = element->GetZ();
133  G4int iz = (G4int)Z;
134  if(iz < 100) {
135  paramC[iz] = mConst*material->GetTotNbOfElectPerVolume();
136  //paramC[iz] = 0.217635e-33*(material->GetTotNbOfElectPerVolume());
137  }
138  if (!(activeZ.contains(Z))) {
139  activeZ.push_back(Z);
140  }
141  }
142  }
143 
144  // Read parameters
145 
146  char* path = std::getenv("G4LEDATA");
147  if (path == 0)
148  {
149  G4Exception("G4BremsstrahlungParameters::LoadData",
150  "em0006",FatalException,"G4LEDATA environment variable not set");
151  return;
152  }
153 
154  G4String pathString_a(path);
155  G4String name_a = pathString_a + name;
156  std::ifstream file_a(name_a);
157  std::filebuf* lsdp_a = file_a.rdbuf();
158 
159  if (! (lsdp_a->is_open()) )
160  {
161  G4String stringConversion2("G4BremsstrahlungParameters::LoadData");
162  G4String excep = stringConversion2 + name_a;
163  G4Exception("G4BremsstrahlungParameters::LoadData",
164  "em0003",FatalException,excep);
165  }
166 
167  // The file is organized into two columns:
168  // 1st column is the energy
169  // 2nd column is the corresponding value
170  // The file terminates with the pattern: -1 -1
171  // -2 -2
172 
173  G4double ener = 0.0;
174  G4double sum = 0.0;
175  G4int z = 0;
176 
177  std::vector<G4DataVector*> a;
178  a.resize(length);
179 
180  G4DataVector e;
181  e.clear();
182 
183  G4bool isReady = false;
184 
185  do {
186  file_a >> ener >> sum;
187 
188  // End of file
189  if (ener == (G4double)(-2)) {
190  break;
191 
192  // End of next element
193  } else if (ener == (G4double)(-1)) {
194 
195  ++z;
196  G4double Z = (G4double)z;
197 
198  // fill map if Z is used
199  if (activeZ.contains(Z)) {
200 
201  for (size_t k=0; k<length; ++k) {
202 
203  G4int id = z*length + k;
205  G4DataVector* eVector = new G4DataVector;
206  size_t eSize = e.size();
207  for (size_t sLocal=0; sLocal<eSize; sLocal++) {
208  eVector->push_back(e[sLocal]);
209  }
210  G4VEMDataSet* set = new G4EMDataSet(id,eVector,a[k],inter,1.,1.);
211  param[id] = set;
212  }
213  } else {
214  for (size_t j=0; j<length; j++) {
215  delete a[j];
216  }
217  }
218  isReady = false;
219 
220  } else {
221 
222  if(!isReady) {
223  isReady = true;
224  e.clear();
225  for (size_t j=0; j<length; ++j) {
226  a[j] = new G4DataVector();
227  }
228  }
229 
230  if(ener > 1000.) ener = 1000.;
231  e.push_back(ener);
232  a[length-1]->push_back(sum);
233 
234  for (size_t j=0; j<length-1; j++) {
235  G4double qRead;
236  file_a >> qRead;
237  a[j]->push_back(qRead);
238  }
239 
240  }
241  } while (ener != (G4double)(-2));
242 
243  file_a.close();
244 
245 }
246 
247 
249 {
250  G4int n = paramC.size();
251  if (id < 0 || id >= n)
252  {
253  G4String stringConversion2(id);
254  G4String ex = "Wrong id " + stringConversion2;
255  G4Exception("G4BremsstrahlungParameters::ParameterC",
256  "em1002",FatalException,ex);
257 
258  }
259 
260  return paramC[id];
261 }
262 
263 
265 {
266 
267  G4cout << G4endl;
268  G4cout << "===== G4BremsstrahlungParameters =====" << G4endl;
269  G4cout << G4endl;
270  G4cout << "===== Parameters =====" << G4endl;
271  G4cout << G4endl;
272 
273  size_t nZ = activeZ.size();
274  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
275 
276  for (size_t j=0; j<nZ; j++) {
277  G4int Z = (G4int)activeZ[j];
278 
279  for (size_t i=0; i<length; i++) {
280 
281  pos = param.find(Z*length + i);
282  if (pos!= param.end()) {
283 
284  G4cout << "===== Z= " << Z
285  << " parameter[" << i << "] ====="
286  << G4endl;
287  G4VEMDataSet* dataSet = (*pos).second;
288  dataSet->PrintData();
289  }
290  }
291  }
292 
293  G4cout << "==========================================" << G4endl;
294 }
295