| 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 // 29 //
30 // History: 30 // History:
31 // ----------- 31 // -----------
32 // 1 Aug 2001 MGP Created 32 // 1 Aug 2001 MGP Created
33 // 09 Oct 2001 VI Add FindValue with 33 // 09 Oct 2001 VI Add FindValue with 3 parameters
34 // + NumberOfComponen 34 // + NumberOfComponents
35 // 19 Jul 2002 VI Create composite d 35 // 19 Jul 2002 VI Create composite data set for material
36 // 21 Jan 2003 VI Cut per region 36 // 21 Jan 2003 VI Cut per region
37 // 37 //
38 // ------------------------------------------- 38 // -------------------------------------------------------------------
39 39
40 #include "G4RDVCrossSectionHandler.hh" 40 #include "G4RDVCrossSectionHandler.hh"
41 #include "G4RDVDataSetAlgorithm.hh" 41 #include "G4RDVDataSetAlgorithm.hh"
42 #include "G4RDLogLogInterpolation.hh" 42 #include "G4RDLogLogInterpolation.hh"
43 #include "G4RDVEMDataSet.hh" 43 #include "G4RDVEMDataSet.hh"
44 #include "G4RDEMDataSet.hh" 44 #include "G4RDEMDataSet.hh"
45 #include "G4RDCompositeEMDataSet.hh" 45 #include "G4RDCompositeEMDataSet.hh"
46 #include "G4RDShellEMDataSet.hh" 46 #include "G4RDShellEMDataSet.hh"
47 #include "G4ProductionCutsTable.hh" 47 #include "G4ProductionCutsTable.hh"
48 #include "G4Material.hh" 48 #include "G4Material.hh"
49 #include "G4Element.hh" 49 #include "G4Element.hh"
50 #include "Randomize.hh" 50 #include "Randomize.hh"
51 #include <map> 51 #include <map>
52 #include <vector> 52 #include <vector>
53 #include <fstream> 53 #include <fstream>
54 #include <sstream> 54 #include <sstream>
55 55
56 56
57 G4RDVCrossSectionHandler::G4RDVCrossSectionHan 57 G4RDVCrossSectionHandler::G4RDVCrossSectionHandler()
58 { 58 {
59 crossSections = 0; 59 crossSections = 0;
60 interpolation = 0; 60 interpolation = 0;
61 Initialise(); 61 Initialise();
62 ActiveElements(); 62 ActiveElements();
63 } 63 }
64 64
65 65
66 G4RDVCrossSectionHandler::G4RDVCrossSectionHan 66 G4RDVCrossSectionHandler::G4RDVCrossSectionHandler(G4RDVDataSetAlgorithm* algorithm,
67 G4double minE, 67 G4double minE,
68 G4double maxE, 68 G4double maxE,
69 G4int bins, 69 G4int bins,
70 G4double unitE, 70 G4double unitE,
71 G4double unitData, 71 G4double unitData,
72 G4int minZ, 72 G4int minZ,
73 G4int maxZ) 73 G4int maxZ)
74 : interpolation(algorithm), eMin(minE), eMax 74 : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
75 unit1(unitE), unit2(unitData), zMin(minZ), 75 unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
76 { 76 {
77 crossSections = 0; 77 crossSections = 0;
78 ActiveElements(); 78 ActiveElements();
79 } 79 }
80 80
81 G4RDVCrossSectionHandler::~G4RDVCrossSectionHa 81 G4RDVCrossSectionHandler::~G4RDVCrossSectionHandler()
82 { 82 {
83 delete interpolation; 83 delete interpolation;
84 interpolation = 0; 84 interpolation = 0;
85 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 85 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
86 86
87 for (pos = dataMap.begin(); pos != dataMap.e 87 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
88 { 88 {
89 // The following is a workaround for STL 89 // The following is a workaround for STL ObjectSpace implementation,
90 // which does not support the standard a 90 // which does not support the standard and does not accept
91 // the syntax pos->second 91 // the syntax pos->second
92 // G4RDVEMDataSet* dataSet = pos->second 92 // G4RDVEMDataSet* dataSet = pos->second;
93 G4RDVEMDataSet* dataSet = (*pos).second; 93 G4RDVEMDataSet* dataSet = (*pos).second;
94 delete dataSet; 94 delete dataSet;
95 } 95 }
96 96
97 if (crossSections != 0) 97 if (crossSections != 0)
98 { 98 {
99 size_t n = crossSections->size(); 99 size_t n = crossSections->size();
100 for (size_t i=0; i<n; i++) 100 for (size_t i=0; i<n; i++)
101 { 101 {
102 delete (*crossSections)[i]; 102 delete (*crossSections)[i];
103 } 103 }
104 delete crossSections; 104 delete crossSections;
105 crossSections = 0; 105 crossSections = 0;
106 } 106 }
107 } 107 }
108 108
109 void G4RDVCrossSectionHandler::Initialise(G4RD 109 void G4RDVCrossSectionHandler::Initialise(G4RDVDataSetAlgorithm* algorithm,
110 G4double minE, G4double maxE, 110 G4double minE, G4double maxE,
111 G4int numberOfBins, 111 G4int numberOfBins,
112 G4double unitE, G4double unitData, 112 G4double unitE, G4double unitData,
113 G4int minZ, G4int maxZ) 113 G4int minZ, G4int maxZ)
114 { 114 {
115 if (algorithm != 0) 115 if (algorithm != 0)
116 { 116 {
117 delete interpolation; 117 delete interpolation;
118 interpolation = algorithm; 118 interpolation = algorithm;
119 } 119 }
120 else 120 else
121 { 121 {
122 interpolation = CreateInterpolation(); 122 interpolation = CreateInterpolation();
123 } 123 }
124 124
125 eMin = minE; 125 eMin = minE;
126 eMax = maxE; 126 eMax = maxE;
127 nBins = numberOfBins; 127 nBins = numberOfBins;
128 unit1 = unitE; 128 unit1 = unitE;
129 unit2 = unitData; 129 unit2 = unitData;
130 zMin = minZ; 130 zMin = minZ;
131 zMax = maxZ; 131 zMax = maxZ;
132 } 132 }
133 133
134 void G4RDVCrossSectionHandler::PrintData() con 134 void G4RDVCrossSectionHandler::PrintData() const
135 { 135 {
136 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 136 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
137 137
138 for (pos = dataMap.begin(); pos != dataMap.e 138 for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
139 { 139 {
140 // The following is a workaround for STL 140 // The following is a workaround for STL ObjectSpace implementation,
141 // which does not support the standard a 141 // which does not support the standard and does not accept
142 // the syntax pos->first or pos->second 142 // the syntax pos->first or pos->second
143 // G4int z = pos->first; 143 // G4int z = pos->first;
144 // G4RDVEMDataSet* dataSet = pos->second 144 // G4RDVEMDataSet* dataSet = pos->second;
145 G4int z = (*pos).first; 145 G4int z = (*pos).first;
146 G4RDVEMDataSet* dataSet = (*pos).second; 146 G4RDVEMDataSet* dataSet = (*pos).second;
147 G4cout << "---- Data set for Z = " 147 G4cout << "---- Data set for Z = "
148 << z 148 << z
149 << G4endl; 149 << G4endl;
150 dataSet->PrintData(); 150 dataSet->PrintData();
151 G4cout << "----------------------------- 151 G4cout << "--------------------------------------------------" << G4endl;
152 } 152 }
153 } 153 }
154 154
155 void G4RDVCrossSectionHandler::LoadData(const 155 void G4RDVCrossSectionHandler::LoadData(const G4String& fileName)
156 { 156 {
157 size_t nZ = activeZ.size(); 157 size_t nZ = activeZ.size();
158 for (size_t i=0; i<nZ; i++) 158 for (size_t i=0; i<nZ; i++)
159 { 159 {
160 G4int Z = (G4int) activeZ[i]; 160 G4int Z = (G4int) activeZ[i];
161 161
162 // Build the complete string identifying 162 // Build the complete string identifying the file with the data set
163 163
164 const char* path = G4FindDataDir("G4LEDA 164 const char* path = G4FindDataDir("G4LEDATA");
165 if (!path) 165 if (!path)
166 { 166 {
167 G4String excep = "G4LEDATA environment var 167 G4String excep = "G4LEDATA environment variable not set!";
168 G4Exception("G4RDVCrossSectionHandler::Loa 168 G4Exception("G4RDVCrossSectionHandler::LoadData()",
169 "InvalidSetup", FatalExc 169 "InvalidSetup", FatalException, excep);
170 } 170 }
171 171
172 std::ostringstream ost; 172 std::ostringstream ost;
173 ost << path << '/' << fileName << Z << " 173 ost << path << '/' << fileName << Z << ".dat";
174 std::ifstream file(ost.str().c_str()); 174 std::ifstream file(ost.str().c_str());
175 std::filebuf* lsdp = file.rdbuf(); 175 std::filebuf* lsdp = file.rdbuf();
176 176
177 if (! (lsdp->is_open()) ) 177 if (! (lsdp->is_open()) )
178 { 178 {
179 G4String excep = "Data file: " + ost.str() 179 G4String excep = "Data file: " + ost.str() + " not found!";
180 G4Exception("G4RDVCrossSectionHandler::Loa 180 G4Exception("G4RDVCrossSectionHandler::LoadData()",
181 "DataNotFound", FatalExc 181 "DataNotFound", FatalException, excep);
182 } 182 }
183 G4double a = 0; 183 G4double a = 0;
184 G4int k = 1; 184 G4int k = 1;
185 G4DataVector* energies = new G4DataVecto 185 G4DataVector* energies = new G4DataVector;
186 G4DataVector* data = new G4DataVector; 186 G4DataVector* data = new G4DataVector;
187 do 187 do
188 { 188 {
189 file >> a; 189 file >> a;
190 G4int nColumns = 2; 190 G4int nColumns = 2;
191 // The file is organized into two columns: 191 // The file is organized into two columns:
192 // 1st column is the energy 192 // 1st column is the energy
193 // 2nd column is the corresponding value 193 // 2nd column is the corresponding value
194 // The file terminates with the pattern: - 194 // The file terminates with the pattern: -1 -1
195 // - 195 // -2 -2
196 if (a == -1 || a == -2) 196 if (a == -1 || a == -2)
197 { 197 {
198 } 198 }
199 else 199 else
200 { 200 {
201 if (k%nColumns != 0) 201 if (k%nColumns != 0)
202 { 202 {
203 G4double e = a * unit1; 203 G4double e = a * unit1;
204 energies->push_back(e); 204 energies->push_back(e);
205 k++; 205 k++;
206 } 206 }
207 else if (k%nColumns == 0) 207 else if (k%nColumns == 0)
208 { 208 {
209 G4double value = a * unit2; 209 G4double value = a * unit2;
210 data->push_back(value); 210 data->push_back(value);
211 k = 1; 211 k = 1;
212 } 212 }
213 } 213 }
214 } while (a != -2); // end of file 214 } while (a != -2); // end of file
215 215
216 file.close(); 216 file.close();
217 G4RDVDataSetAlgorithm* algo = interpolat 217 G4RDVDataSetAlgorithm* algo = interpolation->Clone();
218 G4RDVEMDataSet* dataSet = new G4RDEMData 218 G4RDVEMDataSet* dataSet = new G4RDEMDataSet(Z,energies,data,algo);
219 dataMap[Z] = dataSet; 219 dataMap[Z] = dataSet;
220 } 220 }
221 } 221 }
222 222
223 void G4RDVCrossSectionHandler::LoadShellData(c 223 void G4RDVCrossSectionHandler::LoadShellData(const G4String& fileName)
224 { 224 {
225 size_t nZ = activeZ.size(); 225 size_t nZ = activeZ.size();
226 for (size_t i=0; i<nZ; i++) 226 for (size_t i=0; i<nZ; i++)
227 { 227 {
228 G4int Z = (G4int) activeZ[i]; 228 G4int Z = (G4int) activeZ[i];
229 229
230 // Riccardo Capra <capra@ge.infn.it>: PL 230 // Riccardo Capra <capra@ge.infn.it>: PLEASE CHECK THE FOLLOWING PIECE OF CODE
231 // "energies" AND "data" G4DataVector AR 231 // "energies" AND "data" G4DataVector ARE ALLOCATED, FILLED IN AND NEVER USED OR
232 // DELETED. WHATSMORE LOADING FILE OPERA 232 // DELETED. WHATSMORE LOADING FILE OPERATIONS WERE DONE BY G4RDShellEMDataSet
233 // EVEN BEFORE THE CHANGES I DID ON THIS 233 // EVEN BEFORE THE CHANGES I DID ON THIS FILE. SO THE FOLLOWING CODE IN MY
234 // OPINION SHOULD BE USELESS AND SHOULD 234 // OPINION SHOULD BE USELESS AND SHOULD PRODUCE A MEMORY LEAK.
235 235
236 // Build the complete string identifying 236 // Build the complete string identifying the file with the data set
237 237
238 const char* path = G4FindDataDir("G4LEDA 238 const char* path = G4FindDataDir("G4LEDATA");
239 if (!path) 239 if (!path)
240 { 240 {
241 G4String excep = "G4LEDATA environment var 241 G4String excep = "G4LEDATA environment variable not set!";
242 G4Exception("G4RDVCrossSectionHandler::Loa 242 G4Exception("G4RDVCrossSectionHandler::LoadShellData()",
243 "InvalidSetup", FatalExc 243 "InvalidSetup", FatalException, excep);
244 } 244 }
245 245
246 std::ostringstream ost; 246 std::ostringstream ost;
247 247
248 ost << path << '/' << fileName << Z << " 248 ost << path << '/' << fileName << Z << ".dat";
249 249
250 std::ifstream file(ost.str().c_str()); 250 std::ifstream file(ost.str().c_str());
251 std::filebuf* lsdp = file.rdbuf(); 251 std::filebuf* lsdp = file.rdbuf();
252 252
253 if (! (lsdp->is_open()) ) 253 if (! (lsdp->is_open()) )
254 { 254 {
255 G4String excep = "Data file: " + ost.str() 255 G4String excep = "Data file: " + ost.str() + " not found!";
256 G4Exception("G4RDVCrossSectionHandler::Loa 256 G4Exception("G4RDVCrossSectionHandler::LoadShellData()",
257 "DataNotFound", FatalExc 257 "DataNotFound", FatalException, excep);
258 } 258 }
259 G4double a = 0; 259 G4double a = 0;
260 G4int k = 1; 260 G4int k = 1;
261 G4DataVector* energies = new G4DataVecto 261 G4DataVector* energies = new G4DataVector;
262 G4DataVector* data = new G4DataVector; 262 G4DataVector* data = new G4DataVector;
263 do 263 do
264 { 264 {
265 file >> a; 265 file >> a;
266 G4int nColumns = 2; 266 G4int nColumns = 2;
267 // The file is organized into two columns: 267 // The file is organized into two columns:
268 // 1st column is the energy 268 // 1st column is the energy
269 // 2nd column is the corresponding value 269 // 2nd column is the corresponding value
270 // The file terminates with the pattern: - 270 // The file terminates with the pattern: -1 -1
271 // - 271 // -2 -2
272 if (a == -1 || a == -2) 272 if (a == -1 || a == -2)
273 { 273 {
274 } 274 }
275 else 275 else
276 { 276 {
277 if (k%nColumns != 0) 277 if (k%nColumns != 0)
278 { 278 {
279 G4double e = a * unit1; 279 G4double e = a * unit1;
280 energies->push_back(e); 280 energies->push_back(e);
281 k++; 281 k++;
282 } 282 }
283 else if (k%nColumns == 0) 283 else if (k%nColumns == 0)
284 { 284 {
285 G4double value = a * unit2; 285 G4double value = a * unit2;
286 data->push_back(value); 286 data->push_back(value);
287 k = 1; 287 k = 1;
288 } 288 }
289 } 289 }
290 } while (a != -2); // end of file 290 } while (a != -2); // end of file
291 291
292 file.close(); 292 file.close();
293 293
294 // Riccardo Capra <capra@ge.infn.it>: EN 294 // Riccardo Capra <capra@ge.infn.it>: END OF CODE THAT IN MY OPINION SHOULD BE
295 // REMOVED. 295 // REMOVED.
296 296
297 G4RDVDataSetAlgorithm* algo = interpolat 297 G4RDVDataSetAlgorithm* algo = interpolation->Clone();
298 G4RDVEMDataSet* dataSet = new G4RDShellE 298 G4RDVEMDataSet* dataSet = new G4RDShellEMDataSet(Z, algo);
299 dataSet->LoadData(fileName); 299 dataSet->LoadData(fileName);
300 dataMap[Z] = dataSet; 300 dataMap[Z] = dataSet;
301 } 301 }
302 } 302 }
303 303
304 void G4RDVCrossSectionHandler::Clear() 304 void G4RDVCrossSectionHandler::Clear()
305 { 305 {
306 // Reset the map of data sets: remove the da 306 // Reset the map of data sets: remove the data sets from the map
307 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 307 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
308 308
309 if(! dataMap.empty()) 309 if(! dataMap.empty())
310 { 310 {
311 for (pos = dataMap.begin(); pos != dat 311 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
312 { 312 {
313 // The following is a workaround for STL O 313 // The following is a workaround for STL ObjectSpace implementation,
314 // which does not support the standard and 314 // which does not support the standard and does not accept
315 // the syntax pos->first or pos->second 315 // the syntax pos->first or pos->second
316 // G4RDVEMDataSet* dataSet = pos->second; 316 // G4RDVEMDataSet* dataSet = pos->second;
317 G4RDVEMDataSet* dataSet = (*pos).second; 317 G4RDVEMDataSet* dataSet = (*pos).second;
318 delete dataSet; 318 delete dataSet;
319 dataSet = 0; 319 dataSet = 0;
320 G4int i = (*pos).first; 320 G4int i = (*pos).first;
321 dataMap[i] = 0; 321 dataMap[i] = 0;
322 } 322 }
323 dataMap.clear(); 323 dataMap.clear();
324 } 324 }
325 325
326 activeZ.clear(); 326 activeZ.clear();
327 ActiveElements(); 327 ActiveElements();
328 } 328 }
329 329
330 G4double G4RDVCrossSectionHandler::FindValue(G 330 G4double G4RDVCrossSectionHandler::FindValue(G4int Z, G4double energy) const
331 { 331 {
332 G4double value = 0.; 332 G4double value = 0.;
333 333
334 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 334 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
335 pos = dataMap.find(Z); 335 pos = dataMap.find(Z);
336 if (pos!= dataMap.end()) 336 if (pos!= dataMap.end())
337 { 337 {
338 // The following is a workaround for STL 338 // The following is a workaround for STL ObjectSpace implementation,
339 // which does not support the standard a 339 // which does not support the standard and does not accept
340 // the syntax pos->first or pos->second 340 // the syntax pos->first or pos->second
341 // G4RDVEMDataSet* dataSet = pos->second 341 // G4RDVEMDataSet* dataSet = pos->second;
342 G4RDVEMDataSet* dataSet = (*pos).second; 342 G4RDVEMDataSet* dataSet = (*pos).second;
343 value = dataSet->FindValue(energy); 343 value = dataSet->FindValue(energy);
344 } 344 }
345 else 345 else
346 { 346 {
347 G4cout << "WARNING: G4RDVCrossSectionHan 347 G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find Z = "
348 << Z << G4endl; 348 << Z << G4endl;
349 } 349 }
350 return value; 350 return value;
351 } 351 }
352 352
353 G4double G4RDVCrossSectionHandler::FindValue(G 353 G4double G4RDVCrossSectionHandler::FindValue(G4int Z, G4double energy,
354 G4i 354 G4int shellIndex) const
355 { 355 {
356 G4double value = 0.; 356 G4double value = 0.;
357 357
358 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 358 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
359 pos = dataMap.find(Z); 359 pos = dataMap.find(Z);
360 if (pos!= dataMap.end()) 360 if (pos!= dataMap.end())
361 { 361 {
362 // The following is a workaround for STL 362 // The following is a workaround for STL ObjectSpace implementation,
363 // which does not support the standard a 363 // which does not support the standard and does not accept
364 // the syntax pos->first or pos->second 364 // the syntax pos->first or pos->second
365 // G4RDVEMDataSet* dataSet = pos->second 365 // G4RDVEMDataSet* dataSet = pos->second;
366 G4RDVEMDataSet* dataSet = (*pos).second; 366 G4RDVEMDataSet* dataSet = (*pos).second;
367 if (shellIndex >= 0) 367 if (shellIndex >= 0)
368 { 368 {
369 G4int nComponents = dataSet->NumberOfCompo 369 G4int nComponents = dataSet->NumberOfComponents();
370 if(shellIndex < nComponents) 370 if(shellIndex < nComponents)
371 // - MGP - Why doesn't it use G4RDVEMDat 371 // - MGP - Why doesn't it use G4RDVEMDataSet::FindValue directly?
372 value = dataSet->GetComponent(shellIndex 372 value = dataSet->GetComponent(shellIndex)->FindValue(energy);
373 else 373 else
374 { 374 {
375 G4cout << "WARNING: G4RDVCrossSectionH 375 G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find"
376 << " shellIndex= " << shellIndex 376 << " shellIndex= " << shellIndex
377 << " for Z= " 377 << " for Z= "
378 << Z << G4endl; 378 << Z << G4endl;
379 } 379 }
380 } else { 380 } else {
381 value = dataSet->FindValue(energy); 381 value = dataSet->FindValue(energy);
382 } 382 }
383 } 383 }
384 else 384 else
385 { 385 {
386 G4cout << "WARNING: G4RDVCrossSectionHan 386 G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find Z = "
387 << Z << G4endl; 387 << Z << G4endl;
388 } 388 }
389 return value; 389 return value;
390 } 390 }
391 391
392 392
393 G4double G4RDVCrossSectionHandler::ValueForMat 393 G4double G4RDVCrossSectionHandler::ValueForMaterial(const G4Material* material,
394 G4double energy) const 394 G4double energy) const
395 { 395 {
396 G4double value = 0.; 396 G4double value = 0.;
397 397
398 const G4ElementVector* elementVector = mater 398 const G4ElementVector* elementVector = material->GetElementVector();
399 const G4double* nAtomsPerVolume = material-> 399 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
400 G4int nElements = material->GetNumberOfEleme 400 G4int nElements = material->GetNumberOfElements();
401 401
402 for (G4int i=0 ; i<nElements ; i++) 402 for (G4int i=0 ; i<nElements ; i++)
403 { 403 {
404 G4int Z = (G4int) (*elementVector)[i]->G 404 G4int Z = (G4int) (*elementVector)[i]->GetZ();
405 G4double elementValue = FindValue(Z,ener 405 G4double elementValue = FindValue(Z,energy);
406 G4double nAtomsVol = nAtomsPerVolume[i]; 406 G4double nAtomsVol = nAtomsPerVolume[i];
407 value += nAtomsVol * elementValue; 407 value += nAtomsVol * elementValue;
408 } 408 }
409 409
410 return value; 410 return value;
411 } 411 }
412 412
413 413
414 G4RDVEMDataSet* G4RDVCrossSectionHandler::Buil 414 G4RDVEMDataSet* G4RDVCrossSectionHandler::BuildMeanFreePathForMaterials(const G4DataVector* energyCuts)
415 { 415 {
416 // Builds a CompositeDataSet containing the 416 // Builds a CompositeDataSet containing the mean free path for each material
417 // in the material table 417 // in the material table
418 418
419 G4DataVector energyVector; 419 G4DataVector energyVector;
420 G4double dBin = std::log10(eMax/eMin) / nBin 420 G4double dBin = std::log10(eMax/eMin) / nBins;
421 421
422 for (G4int i=0; i<nBins+1; i++) 422 for (G4int i=0; i<nBins+1; i++)
423 { 423 {
424 energyVector.push_back(std::pow(10., std 424 energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
425 } 425 }
426 426
427 // Factory method to build cross sections in 427 // Factory method to build cross sections in derived classes,
428 // related to the type of physics process 428 // related to the type of physics process
429 429
430 if (crossSections != 0) 430 if (crossSections != 0)
431 { // Reset the list of cross sections 431 { // Reset the list of cross sections
432 std::vector<G4RDVEMDataSet*>::iterator m 432 std::vector<G4RDVEMDataSet*>::iterator mat;
433 if (! crossSections->empty()) 433 if (! crossSections->empty())
434 { 434 {
435 for (mat = crossSections->begin(); mat!= c 435 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
436 { 436 {
437 G4RDVEMDataSet* set = *mat; 437 G4RDVEMDataSet* set = *mat;
438 delete set; 438 delete set;
439 set = 0; 439 set = 0;
440 } 440 }
441 crossSections->clear(); 441 crossSections->clear();
442 delete crossSections; 442 delete crossSections;
443 crossSections = 0; 443 crossSections = 0;
444 } 444 }
445 } 445 }
446 446
447 crossSections = BuildCrossSectionsForMateria 447 crossSections = BuildCrossSectionsForMaterials(energyVector,energyCuts);
448 448
449 if (crossSections == 0) 449 if (crossSections == 0)
450 G4Exception("G4RDVCrossSectionHandler::Bui 450 G4Exception("G4RDVCrossSectionHandler::BuildMeanFreePathForMaterials()",
451 "InvalidCondition", FatalExcep 451 "InvalidCondition", FatalException, "CrossSections = 0!");
452 452
453 G4RDVDataSetAlgorithm* algo = CreateInterpol 453 G4RDVDataSetAlgorithm* algo = CreateInterpolation();
454 G4RDVEMDataSet* materialSet = new G4RDCompos 454 G4RDVEMDataSet* materialSet = new G4RDCompositeEMDataSet(algo);
455 455
456 G4DataVector* energies; 456 G4DataVector* energies;
457 G4DataVector* data; 457 G4DataVector* data;
458 458
459 const G4ProductionCutsTable* theCoupleTable= 459 const G4ProductionCutsTable* theCoupleTable=
460 G4ProductionCutsTable::GetProductionCu 460 G4ProductionCutsTable::GetProductionCutsTable();
461 size_t numOfCouples = theCoupleTable->GetTab 461 size_t numOfCouples = theCoupleTable->GetTableSize();
462 462
463 463
464 for (size_t m=0; m<numOfCouples; m++) 464 for (size_t m=0; m<numOfCouples; m++)
465 { 465 {
466 energies = new G4DataVector; 466 energies = new G4DataVector;
467 data = new G4DataVector; 467 data = new G4DataVector;
468 for (G4int bin=0; bin<nBins; bin++) 468 for (G4int bin=0; bin<nBins; bin++)
469 { 469 {
470 G4double energy = energyVector[bin]; 470 G4double energy = energyVector[bin];
471 energies->push_back(energy); 471 energies->push_back(energy);
472 G4RDVEMDataSet* matCrossSet = (*crossSecti 472 G4RDVEMDataSet* matCrossSet = (*crossSections)[m];
473 G4double materialCrossSection = 0.0; 473 G4double materialCrossSection = 0.0;
474 G4int nElm = matCrossSet->NumberOfCo 474 G4int nElm = matCrossSet->NumberOfComponents();
475 for(G4int j=0; j<nElm; j++) { 475 for(G4int j=0; j<nElm; j++) {
476 materialCrossSection += matCrossSe 476 materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
477 } 477 }
478 478
479 if (materialCrossSection > 0.) 479 if (materialCrossSection > 0.)
480 { 480 {
481 data->push_back(1./materialCrossSectio 481 data->push_back(1./materialCrossSection);
482 } 482 }
483 else 483 else
484 { 484 {
485 data->push_back(DBL_MAX); 485 data->push_back(DBL_MAX);
486 } 486 }
487 } 487 }
488 G4RDVDataSetAlgorithm* algo = CreateInte 488 G4RDVDataSetAlgorithm* algo = CreateInterpolation();
489 G4RDVEMDataSet* dataSet = new G4RDEMData 489 G4RDVEMDataSet* dataSet = new G4RDEMDataSet(m,energies,data,algo,1.,1.);
490 materialSet->AddComponent(dataSet); 490 materialSet->AddComponent(dataSet);
491 } 491 }
492 492
493 return materialSet; 493 return materialSet;
494 } 494 }
495 495
496 G4int G4RDVCrossSectionHandler::SelectRandomAt 496 G4int G4RDVCrossSectionHandler::SelectRandomAtom(const G4MaterialCutsCouple* couple,
497 497 G4double e) const
498 { 498 {
499 // Select randomly an element within the mat 499 // Select randomly an element within the material, according to the weight
500 // determined by the cross sections in the d 500 // determined by the cross sections in the data set
501 501
502 const G4Material* material = couple->GetMate 502 const G4Material* material = couple->GetMaterial();
503 G4int nElements = material->GetNumberOfEleme 503 G4int nElements = material->GetNumberOfElements();
504 504
505 // Special case: the material consists of on 505 // Special case: the material consists of one element
506 if (nElements == 1) 506 if (nElements == 1)
507 { 507 {
508 G4int Z = (G4int) material->GetZ(); 508 G4int Z = (G4int) material->GetZ();
509 return Z; 509 return Z;
510 } 510 }
511 511
512 // Composite material 512 // Composite material
513 513
514 const G4ElementVector* elementVector = mater 514 const G4ElementVector* elementVector = material->GetElementVector();
515 size_t materialIndex = couple->GetIndex(); 515 size_t materialIndex = couple->GetIndex();
516 516
517 G4RDVEMDataSet* materialSet = (*crossSection 517 G4RDVEMDataSet* materialSet = (*crossSections)[materialIndex];
518 G4double materialCrossSection0 = 0.0; 518 G4double materialCrossSection0 = 0.0;
519 G4DataVector cross; 519 G4DataVector cross;
520 cross.clear(); 520 cross.clear();
521 for ( G4int i=0; i < nElements; i++ ) 521 for ( G4int i=0; i < nElements; i++ )
522 { 522 {
523 G4double cr = materialSet->GetComponent( 523 G4double cr = materialSet->GetComponent(i)->FindValue(e);
524 materialCrossSection0 += cr; 524 materialCrossSection0 += cr;
525 cross.push_back(materialCrossSection0); 525 cross.push_back(materialCrossSection0);
526 } 526 }
527 527
528 G4double random = G4UniformRand() * material 528 G4double random = G4UniformRand() * materialCrossSection0;
529 529
530 for (G4int k=0 ; k < nElements ; k++ ) 530 for (G4int k=0 ; k < nElements ; k++ )
531 { 531 {
532 if (random <= cross[k]) return (G4int) ( 532 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
533 } 533 }
534 // It should never get here 534 // It should never get here
535 return 0; 535 return 0;
536 } 536 }
537 537
538 const G4Element* G4RDVCrossSectionHandler::Sel 538 const G4Element* G4RDVCrossSectionHandler::SelectRandomElement(const G4MaterialCutsCouple* couple,
539 G4double e) const 539 G4double e) const
540 { 540 {
541 // Select randomly an element within the mat 541 // Select randomly an element within the material, according to the weight determined
542 // by the cross sections in the data set 542 // by the cross sections in the data set
543 543
544 const G4Material* material = couple->GetMate 544 const G4Material* material = couple->GetMaterial();
545 G4Element* nullElement = 0; 545 G4Element* nullElement = 0;
546 G4int nElements = material->GetNumberOfEleme 546 G4int nElements = material->GetNumberOfElements();
547 const G4ElementVector* elementVector = mater 547 const G4ElementVector* elementVector = material->GetElementVector();
548 548
549 // Special case: the material consists of on 549 // Special case: the material consists of one element
550 if (nElements == 1) 550 if (nElements == 1)
551 { 551 {
552 G4Element* element = (*elementVector)[0] 552 G4Element* element = (*elementVector)[0];
553 return element; 553 return element;
554 } 554 }
555 else 555 else
556 { 556 {
557 // Composite material 557 // Composite material
558 558
559 size_t materialIndex = couple->GetIndex( 559 size_t materialIndex = couple->GetIndex();
560 560
561 G4RDVEMDataSet* materialSet = (*crossSec 561 G4RDVEMDataSet* materialSet = (*crossSections)[materialIndex];
562 G4double materialCrossSection0 = 0.0; 562 G4double materialCrossSection0 = 0.0;
563 G4DataVector cross; 563 G4DataVector cross;
564 cross.clear(); 564 cross.clear();
565 for (G4int i=0; i<nElements; i++) 565 for (G4int i=0; i<nElements; i++)
566 { 566 {
567 G4double cr = materialSet->GetCompon 567 G4double cr = materialSet->GetComponent(i)->FindValue(e);
568 materialCrossSection0 += cr; 568 materialCrossSection0 += cr;
569 cross.push_back(materialCrossSection 569 cross.push_back(materialCrossSection0);
570 } 570 }
571 571
572 G4double random = G4UniformRand() * mate 572 G4double random = G4UniformRand() * materialCrossSection0;
573 573
574 for (G4int k=0 ; k < nElements ; k++ ) 574 for (G4int k=0 ; k < nElements ; k++ )
575 { 575 {
576 if (random <= cross[k]) return (*ele 576 if (random <= cross[k]) return (*elementVector)[k];
577 } 577 }
578 // It should never end up here 578 // It should never end up here
579 G4cout << "G4RDVCrossSectionHandler::Sel 579 G4cout << "G4RDVCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
580 return nullElement; 580 return nullElement;
581 } 581 }
582 } 582 }
583 583
584 G4int G4RDVCrossSectionHandler::SelectRandomSh 584 G4int G4RDVCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
585 { 585 {
586 // Select randomly a shell, according to the 586 // Select randomly a shell, according to the weight determined by the cross sections
587 // in the data set 587 // in the data set
588 588
589 // Note for later improvement: it would be u 589 // Note for later improvement: it would be useful to add a cache mechanism for already
590 // used shells to improve performance 590 // used shells to improve performance
591 591
592 G4int shell = 0; 592 G4int shell = 0;
593 593
594 G4double totCrossSection = FindValue(Z,e); 594 G4double totCrossSection = FindValue(Z,e);
595 G4double random = G4UniformRand() * totCross 595 G4double random = G4UniformRand() * totCrossSection;
596 G4double partialSum = 0.; 596 G4double partialSum = 0.;
597 597
598 G4RDVEMDataSet* dataSet = 0; 598 G4RDVEMDataSet* dataSet = 0;
599 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 599 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
600 pos = dataMap.find(Z); 600 pos = dataMap.find(Z);
601 // The following is a workaround for STL Obj 601 // The following is a workaround for STL ObjectSpace implementation,
602 // which does not support the standard and d 602 // which does not support the standard and does not accept
603 // the syntax pos->first or pos->second 603 // the syntax pos->first or pos->second
604 // if (pos != dataMap.end()) dataSet = pos-> 604 // if (pos != dataMap.end()) dataSet = pos->second;
605 if (pos != dataMap.end()) dataSet = (*pos).s 605 if (pos != dataMap.end()) dataSet = (*pos).second;
606 606
607 size_t nShells = dataSet->NumberOfComponents 607 size_t nShells = dataSet->NumberOfComponents();
608 for (size_t i=0; i<nShells; i++) 608 for (size_t i=0; i<nShells; i++)
609 { 609 {
610 const G4RDVEMDataSet* shellDataSet = dat 610 const G4RDVEMDataSet* shellDataSet = dataSet->GetComponent(i);
611 if (shellDataSet != 0) 611 if (shellDataSet != 0)
612 { 612 {
613 G4double value = shellDataSet->FindValue(e 613 G4double value = shellDataSet->FindValue(e);
614 partialSum += value; 614 partialSum += value;
615 if (random <= partialSum) return i; 615 if (random <= partialSum) return i;
616 } 616 }
617 } 617 }
618 // It should never get here 618 // It should never get here
619 return shell; 619 return shell;
620 } 620 }
621 621
622 void G4RDVCrossSectionHandler::ActiveElements( 622 void G4RDVCrossSectionHandler::ActiveElements()
623 { 623 {
624 const G4MaterialTable* materialTable = G4Mat 624 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
625 if (materialTable == 0) 625 if (materialTable == 0)
626 G4Exception("G4RDVCrossSectionHandler::Act 626 G4Exception("G4RDVCrossSectionHandler::ActiveElements",
627 "InvalidSetup", FatalException 627 "InvalidSetup", FatalException, "No MaterialTable found!");
628 628
629 G4int nMaterials = G4Material::GetNumberOfMa 629 G4int nMaterials = G4Material::GetNumberOfMaterials();
630 630
631 for (G4int m=0; m<nMaterials; m++) 631 for (G4int m=0; m<nMaterials; m++)
632 { 632 {
633 const G4Material* material= (*materialTa 633 const G4Material* material= (*materialTable)[m];
634 const G4ElementVector* elementVector = m 634 const G4ElementVector* elementVector = material->GetElementVector();
635 const G4int nElements = material->GetNum 635 const G4int nElements = material->GetNumberOfElements();
636 636
637 for (G4int iEl=0; iEl<nElements; iEl++) 637 for (G4int iEl=0; iEl<nElements; iEl++)
638 { 638 {
639 G4Element* element = (*elementVector)[iEl] 639 G4Element* element = (*elementVector)[iEl];
640 G4double Z = element->GetZ(); 640 G4double Z = element->GetZ();
641 if (!(activeZ.contains(Z)) && Z >= zMin && 641 if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
642 { 642 {
643 activeZ.push_back(Z); 643 activeZ.push_back(Z);
644 } 644 }
645 } 645 }
646 } 646 }
647 } 647 }
648 648
649 G4RDVDataSetAlgorithm* G4RDVCrossSectionHandle 649 G4RDVDataSetAlgorithm* G4RDVCrossSectionHandler::CreateInterpolation()
650 { 650 {
651 G4RDVDataSetAlgorithm* algorithm = new G4RDL 651 G4RDVDataSetAlgorithm* algorithm = new G4RDLogLogInterpolation;
652 return algorithm; 652 return algorithm;
653 } 653 }
654 654
655 G4int G4RDVCrossSectionHandler::NumberOfCompon 655 G4int G4RDVCrossSectionHandler::NumberOfComponents(G4int Z) const
656 { 656 {
657 G4int n = 0; 657 G4int n = 0;
658 658
659 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 659 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
660 pos = dataMap.find(Z); 660 pos = dataMap.find(Z);
661 if (pos!= dataMap.end()) 661 if (pos!= dataMap.end())
662 { 662 {
663 G4RDVEMDataSet* dataSet = (*pos).second; 663 G4RDVEMDataSet* dataSet = (*pos).second;
664 n = dataSet->NumberOfComponents(); 664 n = dataSet->NumberOfComponents();
665 } 665 }
666 else 666 else
667 { 667 {
668 G4cout << "WARNING: G4RDVCrossSectionHan 668 G4cout << "WARNING: G4RDVCrossSectionHandler::NumberOfComponents did not "
669 << "find Z = " 669 << "find Z = "
670 << Z << G4endl; 670 << Z << G4endl;
671 } 671 }
672 return n; 672 return n;
673 } 673 }
674 674
675 675
676 676