ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RDBremsstrahlungParameters.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RDBremsstrahlungParameters.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 //
41 // -------------------------------------------------------------------
42 
44 #include "G4RDVEMDataSet.hh"
45 #include "G4RDEMDataSet.hh"
47 #include "G4Material.hh"
48 #include <fstream>
49 
50 
52  size_t num, G4int minZ, G4int maxZ)
53  : zMin(minZ),
54  zMax(maxZ),
55  length(num)
56 {
57  LoadData(name);
58 }
59 
60 
62 {
63  // Reset the map of data sets: remove the data sets from the map
64  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
65 
66  for (pos = param.begin(); pos != param.end(); ++pos)
67  {
68  G4RDVEMDataSet* dataSet = (*pos).second;
69  delete dataSet;
70  }
71 
72  activeZ.clear();
73  paramC.clear();
74 }
75 
76 
78  G4int Z,
79  G4double energy) const
80 {
81  G4double value = 0.;
82  G4int id = Z*length + parameterIndex;
83  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
84 
85  pos = param.find(id);
86  if (pos!= param.end()) {
87 
88  G4RDVEMDataSet* dataSet = (*pos).second;
89  const G4DataVector ener = dataSet->GetEnergies(0);
90  G4double ee = std::max(ener.front(),std::min(ener.back(),energy));
91  value = dataSet->FindValue(ee);
92 
93  } else {
94  G4cout << "WARNING: G4RDBremsstrahlungParameters::FindValue "
95  << "did not find ID = "
96  << id << G4endl;
97  }
98 
99  return value;
100 }
101 
103 {
104  // Build the complete string identifying the file with the data set
105 
106  // define active elements
107 
108  const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
109  if (materialTable == 0)
110  G4Exception("G4RDBremsstrahlungParameters::LoadData()",
111  "DataNotFound", FatalException, "No MaterialTable found!");
112 
114 
115  G4double x = 1.e-9;
116  for (G4int mm=0; mm<100; mm++) {
117  paramC.push_back(x);
118  }
119 
120  for (G4int m=0; m<nMaterials; m++) {
121 
122  const G4Material* material= (*materialTable)[m];
123  const G4ElementVector* elementVector = material->GetElementVector();
124  const G4int nElements = material->GetNumberOfElements();
125 
126  for (G4int iEl=0; iEl<nElements; iEl++) {
127  G4Element* element = (*elementVector)[iEl];
128  G4double Z = element->GetZ();
129  G4int iz = (G4int)Z;
130  if(iz < 100)
131  paramC[iz] = 0.217635e-33*(material->GetTotNbOfElectPerVolume());
132  if (!(activeZ.contains(Z))) {
133  activeZ.push_back(Z);
134  }
135  }
136  }
137 
138  // Read parameters
139 
140  char* path = std::getenv("G4LEDATA");
141  if (path == 0)
142  {
143  G4String excep("G4LEDATA environment variable not set!");
144  G4Exception("G4RDBremsstrahlungParameters::LoadData()",
145  "InvalidSetup", FatalException, excep);
146  }
147 
148  G4String pathString_a(path);
149  G4String name_a = pathString_a + name;
150  std::ifstream file_a(name_a);
151  std::filebuf* lsdp_a = file_a.rdbuf();
152 
153  if (! (lsdp_a->is_open()) )
154  {
155  G4String stringConversion2("Cannot open file ");
156  G4String excep = stringConversion2 + name_a;
157  G4Exception("G4RDBremsstrahlungParameters::LoadData()",
158  "FileNotFound", FatalException, excep);
159  }
160 
161  // The file is organized into two columns:
162  // 1st column is the energy
163  // 2nd column is the corresponding value
164  // The file terminates with the pattern: -1 -1
165  // -2 -2
166 
167  G4double ener = 0.0;
168  G4double sum = 0.0;
169  G4int z = 0;
170 
171  std::vector<G4DataVector*> a;
172  for (size_t j=0; j<length; j++) {
173  G4DataVector* aa = new G4DataVector();
174  a.push_back(aa);
175  }
176  G4DataVector e;
177  e.clear();
178 
179  do {
180  file_a >> ener >> sum;
181 
182  // End of file
183  if (ener == (G4double)(-2)) {
184  break;
185 
186  // End of next element
187  } else if (ener == (G4double)(-1)) {
188 
189  z++;
190  G4double Z = (G4double)z;
191 
192  // fill map if Z is used
193  if (activeZ.contains(Z)) {
194 
195  for (size_t k=0; k<length; k++) {
196 
197  G4int id = z*length + k;
199  G4DataVector* eVector = new G4DataVector;
200  size_t eSize = e.size();
201  for (size_t s=0; s<eSize; s++) {
202  eVector->push_back(e[s]);
203  }
204  G4RDVEMDataSet* set = new G4RDEMDataSet(id,eVector,a[k],inter,1.,1.);
205  param[id] = set;
206  }
207  a.clear();
208  for (size_t j=0; j<length; j++) {
209  G4DataVector* aa = new G4DataVector();
210  a.push_back(aa);
211  }
212  } else {
213  for (size_t j=0; j<length; j++) {
214  a[j]->clear();
215  }
216  }
217  e.clear();
218 
219  } else {
220 
221  if(ener > 1000.) ener = 1000.;
222  e.push_back(ener);
223  a[length-1]->push_back(sum);
224 
225  for (size_t j=0; j<length-1; j++) {
226  G4double qRead;
227  file_a >> qRead;
228  a[j]->push_back(qRead);
229  }
230 
231  }
232  } while (ener != (G4double)(-2));
233 
234  file_a.close();
235 
236 }
237 
238 
240 {
241  G4int n = paramC.size();
242  if (id < 0 || id >= n)
243  {
244  G4String stringConversion1("Wrong id = ");
245  G4String stringConversion2(id);
246  G4String ex = stringConversion1 + stringConversion2;
247  G4Exception("G4RDBremsstrahlungParameters::ParameterC()",
248  "InvalidSetup", FatalException, ex);
249  }
250 
251  return paramC[id];
252 }
253 
254 
256 {
257 
258  G4cout << G4endl;
259  G4cout << "===== G4RDBremsstrahlungParameters =====" << G4endl;
260  G4cout << G4endl;
261  G4cout << "===== Parameters =====" << G4endl;
262  G4cout << G4endl;
263 
264  size_t nZ = activeZ.size();
265  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
266 
267  for (size_t j=0; j<nZ; j++) {
268  G4int Z = (G4int)activeZ[j];
269 
270  for (size_t i=0; i<length; i++) {
271 
272  pos = param.find(Z*length + i);
273  if (pos!= param.end()) {
274 
275  G4cout << "===== Z= " << Z
276  << " parameter[" << i << "] ====="
277  << G4endl;
278  G4RDVEMDataSet* dataSet = (*pos).second;
279  dataSet->PrintData();
280  }
281  }
282  }
283 
284  G4cout << "==========================================" << G4endl;
285 }
286