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