| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // V.Ivanchenko (Vladimir.Ivantchenko@ 29 // V.Ivanchenko (Vladimir.Ivantchenko@cern.ch)
30 // 30 //
31 // History: 31 // History:
32 // ----------- 32 // -----------
33 // 31 Jul 2001 MGP Created 33 // 31 Jul 2001 MGP Created
34 // 12.09.01 V.Ivanchenko Add activeZ and pa 34 // 12.09.01 V.Ivanchenko Add activeZ and paramA
35 // 25.09.01 V.Ivanchenko Add parameter C an 35 // 25.09.01 V.Ivanchenko Add parameter C and change interface to B
36 // 29.11.01 V.Ivanchenko Update parametrisa 36 // 29.11.01 V.Ivanchenko Update parametrisation
37 // 18.11.02 V.Ivanchenko Fix problem of loa 37 // 18.11.02 V.Ivanchenko Fix problem of load
38 // 21.02.03 V.Ivanchenko Number of paramete 38 // 21.02.03 V.Ivanchenko Number of parameters is defined in the constructor
39 // 28.02.03 V.Ivanchenko Filename is define 39 // 28.02.03 V.Ivanchenko Filename is defined in the constructor
40 // 40 //
41 // ------------------------------------------- 41 // -------------------------------------------------------------------
42 42
43 #include "G4RDBremsstrahlungParameters.hh" 43 #include "G4RDBremsstrahlungParameters.hh"
44 #include "G4RDVEMDataSet.hh" 44 #include "G4RDVEMDataSet.hh"
45 #include "G4RDEMDataSet.hh" 45 #include "G4RDEMDataSet.hh"
46 #include "G4RDLogLogInterpolation.hh" 46 #include "G4RDLogLogInterpolation.hh"
47 #include "G4Material.hh" 47 #include "G4Material.hh"
48 #include <fstream> 48 #include <fstream>
49 49
50 50
51 G4RDBremsstrahlungParameters:: G4RDBremsstrahl 51 G4RDBremsstrahlungParameters:: G4RDBremsstrahlungParameters(const G4String& name,
52 size_t num, G4int minZ, G4int maxZ) 52 size_t num, G4int minZ, G4int maxZ)
53 : zMin(minZ), 53 : zMin(minZ),
54 zMax(maxZ), 54 zMax(maxZ),
55 length(num) 55 length(num)
56 { 56 {
57 LoadData(name); 57 LoadData(name);
58 } 58 }
59 59
60 60
61 G4RDBremsstrahlungParameters::~G4RDBremsstrahl 61 G4RDBremsstrahlungParameters::~G4RDBremsstrahlungParameters()
62 { 62 {
63 // Reset the map of data sets: remove the da 63 // Reset the map of data sets: remove the data sets from the map
64 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 64 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
65 65
66 for (pos = param.begin(); pos != param.end() 66 for (pos = param.begin(); pos != param.end(); ++pos)
67 { 67 {
68 G4RDVEMDataSet* dataSet = (*pos).second; 68 G4RDVEMDataSet* dataSet = (*pos).second;
69 delete dataSet; 69 delete dataSet;
70 } 70 }
71 71
72 activeZ.clear(); 72 activeZ.clear();
73 paramC.clear(); 73 paramC.clear();
74 } 74 }
75 75
76 76
77 G4double G4RDBremsstrahlungParameters::Paramet 77 G4double G4RDBremsstrahlungParameters::Parameter(G4int parameterIndex,
78 78 G4int Z,
79 79 G4double energy) const
80 { 80 {
81 G4double value = 0.; 81 G4double value = 0.;
82 G4int id = Z*length + parameterIndex; 82 G4int id = Z*length + parameterIndex;
83 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 83 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
84 84
85 pos = param.find(id); 85 pos = param.find(id);
86 if (pos!= param.end()) { 86 if (pos!= param.end()) {
87 87
88 G4RDVEMDataSet* dataSet = (*pos).second; 88 G4RDVEMDataSet* dataSet = (*pos).second;
89 const G4DataVector ener = dataSet->GetEner 89 const G4DataVector ener = dataSet->GetEnergies(0);
90 G4double ee = std::max(ener.front(),std::m 90 G4double ee = std::max(ener.front(),std::min(ener.back(),energy));
91 value = dataSet->FindValue(ee); 91 value = dataSet->FindValue(ee);
92 92
93 } else { 93 } else {
94 G4cout << "WARNING: G4RDBremsstrahlungPara 94 G4cout << "WARNING: G4RDBremsstrahlungParameters::FindValue "
95 << "did not find ID = " 95 << "did not find ID = "
96 << id << G4endl; 96 << id << G4endl;
97 } 97 }
98 98
99 return value; 99 return value;
100 } 100 }
101 101
102 void G4RDBremsstrahlungParameters::LoadData(co 102 void G4RDBremsstrahlungParameters::LoadData(const G4String& name)
103 { 103 {
104 // Build the complete string identifying the 104 // Build the complete string identifying the file with the data set
105 105
106 // define active elements 106 // define active elements
107 107
108 const G4MaterialTable* materialTable = G4Mat 108 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
109 if (materialTable == 0) 109 if (materialTable == 0)
110 G4Exception("G4RDBremsstrahlungParameters 110 G4Exception("G4RDBremsstrahlungParameters::LoadData()",
111 "DataNotFound", FatalExceptio 111 "DataNotFound", FatalException, "No MaterialTable found!");
112 112
113 G4int nMaterials = G4Material::GetNumberOfMa 113 G4int nMaterials = G4Material::GetNumberOfMaterials();
114 114
115 G4double x = 1.e-9; 115 G4double x = 1.e-9;
116 for (G4int mm=0; mm<100; mm++) { 116 for (G4int mm=0; mm<100; mm++) {
117 paramC.push_back(x); 117 paramC.push_back(x);
118 } 118 }
119 119
120 for (G4int m=0; m<nMaterials; m++) { 120 for (G4int m=0; m<nMaterials; m++) {
121 121
122 const G4Material* material= (*materialTabl 122 const G4Material* material= (*materialTable)[m];
123 const G4ElementVector* elementVector = mat 123 const G4ElementVector* elementVector = material->GetElementVector();
124 const G4int nElements = material->GetNumbe 124 const G4int nElements = material->GetNumberOfElements();
125 125
126 for (G4int iEl=0; iEl<nElements; iEl++) { 126 for (G4int iEl=0; iEl<nElements; iEl++) {
127 G4Element* element = (*elementVector)[iE 127 G4Element* element = (*elementVector)[iEl];
128 G4double Z = element->GetZ(); 128 G4double Z = element->GetZ();
129 G4int iz = (G4int)Z; 129 G4int iz = (G4int)Z;
130 if(iz < 100) 130 if(iz < 100)
131 paramC[iz] = 0.217635e-33*(materia 131 paramC[iz] = 0.217635e-33*(material->GetTotNbOfElectPerVolume());
132 if (!(activeZ.contains(Z))) { 132 if (!(activeZ.contains(Z))) {
133 activeZ.push_back(Z); 133 activeZ.push_back(Z);
134 } 134 }
135 } 135 }
136 } 136 }
137 137
138 // Read parameters 138 // Read parameters
139 139
140 const char* path = G4FindDataDir("G4LEDATA") << 140 char* path = std::getenv("G4LEDATA");
141 if (path == 0) 141 if (path == 0)
142 { 142 {
143 G4String excep("G4LEDATA environment var 143 G4String excep("G4LEDATA environment variable not set!");
144 G4Exception("G4RDBremsstrahlungParameter 144 G4Exception("G4RDBremsstrahlungParameters::LoadData()",
145 "InvalidSetup", FatalExcepti 145 "InvalidSetup", FatalException, excep);
146 } 146 }
147 147
148 G4String pathString_a(path); 148 G4String pathString_a(path);
149 G4String name_a = pathString_a + name; 149 G4String name_a = pathString_a + name;
150 std::ifstream file_a(name_a); 150 std::ifstream file_a(name_a);
151 std::filebuf* lsdp_a = file_a.rdbuf(); 151 std::filebuf* lsdp_a = file_a.rdbuf();
152 152
153 if (! (lsdp_a->is_open()) ) 153 if (! (lsdp_a->is_open()) )
154 { 154 {
155 G4String stringConversion2("Cannot open 155 G4String stringConversion2("Cannot open file ");
156 G4String excep = stringConversion2 + nam 156 G4String excep = stringConversion2 + name_a;
157 G4Exception("G4RDBremsstrahlungParameter 157 G4Exception("G4RDBremsstrahlungParameters::LoadData()",
158 "FileNotFound", FatalExcepti 158 "FileNotFound", FatalException, excep);
159 } 159 }
160 160
161 // The file is organized into two columns: 161 // The file is organized into two columns:
162 // 1st column is the energy 162 // 1st column is the energy
163 // 2nd column is the corresponding value 163 // 2nd column is the corresponding value
164 // The file terminates with the pattern: -1 164 // The file terminates with the pattern: -1 -1
165 // -2 165 // -2 -2
166 166
167 G4double ener = 0.0; 167 G4double ener = 0.0;
168 G4double sum = 0.0; 168 G4double sum = 0.0;
169 G4int z = 0; 169 G4int z = 0;
170 170
171 std::vector<G4DataVector*> a; 171 std::vector<G4DataVector*> a;
172 for (size_t j=0; j<length; j++) { 172 for (size_t j=0; j<length; j++) {
173 G4DataVector* aa = new G4DataVector(); 173 G4DataVector* aa = new G4DataVector();
174 a.push_back(aa); 174 a.push_back(aa);
175 } 175 }
176 G4DataVector e; 176 G4DataVector e;
177 e.clear(); 177 e.clear();
178 178
179 do { 179 do {
180 file_a >> ener >> sum; 180 file_a >> ener >> sum;
181 181
182 // End of file 182 // End of file
183 if (ener == (G4double)(-2)) { 183 if (ener == (G4double)(-2)) {
184 break; 184 break;
185 185
186 // End of next element 186 // End of next element
187 } else if (ener == (G4double)(-1)) { 187 } else if (ener == (G4double)(-1)) {
188 188
189 z++; 189 z++;
190 G4double Z = (G4double)z; 190 G4double Z = (G4double)z;
191 191
192 // fill map if Z is used 192 // fill map if Z is used
193 if (activeZ.contains(Z)) { 193 if (activeZ.contains(Z)) {
194 194
195 for (size_t k=0; k<length; k++) { 195 for (size_t k=0; k<length; k++) {
196 196
197 G4int id = z*length + k; 197 G4int id = z*length + k;
198 G4RDVDataSetAlgorithm* inter = new G4RD 198 G4RDVDataSetAlgorithm* inter = new G4RDLogLogInterpolation();
199 G4DataVector* eVector = new G4DataVector 199 G4DataVector* eVector = new G4DataVector;
200 size_t eSize = e.size(); 200 size_t eSize = e.size();
201 for (size_t s=0; s<eSize; s++) { 201 for (size_t s=0; s<eSize; s++) {
202 eVector->push_back(e[s]); 202 eVector->push_back(e[s]);
203 } 203 }
204 G4RDVEMDataSet* set = new G4RDEMDataSet( 204 G4RDVEMDataSet* set = new G4RDEMDataSet(id,eVector,a[k],inter,1.,1.);
205 param[id] = set; 205 param[id] = set;
206 } 206 }
207 a.clear(); 207 a.clear();
208 for (size_t j=0; j<length; j++) { 208 for (size_t j=0; j<length; j++) {
209 G4DataVector* aa = new G4DataVecto 209 G4DataVector* aa = new G4DataVector();
210 a.push_back(aa); 210 a.push_back(aa);
211 } 211 }
212 } else { 212 } else {
213 for (size_t j=0; j<length; j++) { 213 for (size_t j=0; j<length; j++) {
214 a[j]->clear(); 214 a[j]->clear();
215 } 215 }
216 } 216 }
217 e.clear(); 217 e.clear();
218 218
219 } else { 219 } else {
220 220
221 if(ener > 1000.) ener = 1000.; 221 if(ener > 1000.) ener = 1000.;
222 e.push_back(ener); 222 e.push_back(ener);
223 a[length-1]->push_back(sum); 223 a[length-1]->push_back(sum);
224 224
225 for (size_t j=0; j<length-1; j++) { 225 for (size_t j=0; j<length-1; j++) {
226 G4double qRead; 226 G4double qRead;
227 file_a >> qRead; 227 file_a >> qRead;
228 a[j]->push_back(qRead); 228 a[j]->push_back(qRead);
229 } 229 }
230 230
231 } 231 }
232 } while (ener != (G4double)(-2)); 232 } while (ener != (G4double)(-2));
233 233
234 file_a.close(); 234 file_a.close();
235 235
236 } 236 }
237 237
238 238
239 G4double G4RDBremsstrahlungParameters::Paramet 239 G4double G4RDBremsstrahlungParameters::ParameterC(G4int id) const
240 { 240 {
241 G4int n = paramC.size(); 241 G4int n = paramC.size();
242 if (id < 0 || id >= n) 242 if (id < 0 || id >= n)
243 { 243 {
244 G4String stringConversion1("Wrong id = " 244 G4String stringConversion1("Wrong id = ");
245 G4String stringConversion2(id); 245 G4String stringConversion2(id);
246 G4String ex = stringConversion1 + string 246 G4String ex = stringConversion1 + stringConversion2;
247 G4Exception("G4RDBremsstrahlungParameter 247 G4Exception("G4RDBremsstrahlungParameters::ParameterC()",
248 "InvalidSetup", FatalExcepti 248 "InvalidSetup", FatalException, ex);
249 } 249 }
250 250
251 return paramC[id]; 251 return paramC[id];
252 } 252 }
253 253
254 254
255 void G4RDBremsstrahlungParameters::PrintData() 255 void G4RDBremsstrahlungParameters::PrintData() const
256 { 256 {
257 257
258 G4cout << G4endl; 258 G4cout << G4endl;
259 G4cout << "===== G4RDBremsstrahlungParameter 259 G4cout << "===== G4RDBremsstrahlungParameters =====" << G4endl;
260 G4cout << G4endl; 260 G4cout << G4endl;
261 G4cout << "===== Parameters =====" << G4endl 261 G4cout << "===== Parameters =====" << G4endl;
262 G4cout << G4endl; 262 G4cout << G4endl;
263 263
264 size_t nZ = activeZ.size(); 264 size_t nZ = activeZ.size();
265 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 265 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
266 266
267 for (size_t j=0; j<nZ; j++) { 267 for (size_t j=0; j<nZ; j++) {
268 G4int Z = (G4int)activeZ[j]; 268 G4int Z = (G4int)activeZ[j];
269 269
270 for (size_t i=0; i<length; i++) { 270 for (size_t i=0; i<length; i++) {
271 271
272 pos = param.find(Z*length + i); 272 pos = param.find(Z*length + i);
273 if (pos!= param.end()) { 273 if (pos!= param.end()) {
274 274
275 G4cout << "===== Z= " << Z 275 G4cout << "===== Z= " << Z
276 << " parameter[" << i << "] 276 << " parameter[" << i << "] ====="
277 << G4endl; 277 << G4endl;
278 G4RDVEMDataSet* dataSet = (*pos).secon 278 G4RDVEMDataSet* dataSet = (*pos).second;
279 dataSet->PrintData(); 279 dataSet->PrintData();
280 } 280 }
281 } 281 }
282 } 282 }
283 283
284 G4cout << "================================= 284 G4cout << "==========================================" << G4endl;
285 } 285 }
286 286
287 287