| Geant4 Cross Reference |
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
43 #include "G4RDBremsstrahlungParameters.hh"
44 #include "G4RDVEMDataSet.hh"
45 #include "G4RDEMDataSet.hh"
46 #include "G4RDLogLogInterpolation.hh"
47 #include "G4Material.hh"
48 #include <fstream>
49
50
51 G4RDBremsstrahlungParameters:: G4RDBremsstrahlungParameters(const G4String& name,
52 size_t num, G4int minZ, G4int maxZ)
53 : zMin(minZ),
54 zMax(maxZ),
55 length(num)
56 {
57 LoadData(name);
58 }
59
60
61 G4RDBremsstrahlungParameters::~G4RDBremsstrahlungParameters()
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
77 G4double G4RDBremsstrahlungParameters::Parameter(G4int parameterIndex,
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
102 void G4RDBremsstrahlungParameters::LoadData(const G4String& name)
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
113 G4int nMaterials = G4Material::GetNumberOfMaterials();
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 const char* path = G4FindDataDir("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;
198 G4RDVDataSetAlgorithm* inter = new G4RDLogLogInterpolation();
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
239 G4double G4RDBremsstrahlungParameters::ParameterC(G4int id) const
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
255 void G4RDBremsstrahlungParameters::PrintData() const
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
287