| 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: G4VCrossSectionHandler.cc 66241 2012-12-13 18:34:42Z gunter $
27 // 28 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 29 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // 30 //
30 // History: 31 // History:
31 // ----------- 32 // -----------
32 // 1 Aug 2001 MGP Created 33 // 1 Aug 2001 MGP Created
33 // 09 Oct 2001 VI Add FindValue with 34 // 09 Oct 2001 VI Add FindValue with 3 parameters
34 // + NumberOfComponen 35 // + NumberOfComponents
35 // 19 Jul 2002 VI Create composite d 36 // 19 Jul 2002 VI Create composite data set for material
36 // 21 Jan 2003 VI Cut per region 37 // 21 Jan 2003 VI Cut per region
37 // 38 //
38 // 15 Jul 2009 Nicolas A. Karakatsanis 39 // 15 Jul 2009 Nicolas A. Karakatsanis
39 // 40 //
40 // - LoadNonLogData 41 // - LoadNonLogData method was created to load only the non-logarithmic data from G4EMLOW
41 // dataset. It is 42 // dataset. It is essentially performing the data loading operations as in the past.
42 // 43 //
43 // - LoadData method 44 // - LoadData method was revised in order to calculate the logarithmic values of the data
44 // It retrieves th 45 // It retrieves the data values from the G4EMLOW data files but, then, calculates the
45 // respective log 46 // respective log values and loads them to seperate data structures.
46 // The EM data set 47 // The EM data sets, initialized this way, contain both non-log and log values.
47 // These initializ 48 // These initialized data sets can enhance the computing performance of data interpolation
48 // operations 49 // operations
49 // 50 //
50 // - BuildMeanFreePa 51 // - BuildMeanFreePathForMaterials method was also revised in order to calculate the
51 // logarithmic val 52 // logarithmic values of the loaded data.
52 // It generates th 53 // It generates the data values and, then, calculates the respective log values which
53 // later load to s 54 // later load to seperate data structures.
54 // The EM data set 55 // The EM data sets, initialized this way, contain both non-log and log values.
55 // These initializ 56 // These initialized data sets can enhance the computing performance of data interpolation
56 // operations 57 // operations
57 // 58 //
58 // - LoadShellData m 59 // - LoadShellData method was revised in order to eliminate the presence of a potential
59 // memory leak ori 60 // memory leak originally identified by Riccardo Capra.
60 // Riccardo Capra 61 // Riccardo Capra Original Comment
61 // Riccardo Capra 62 // Riccardo Capra <capra@ge.infn.it>: PLEASE CHECK THE FOLLOWING PIECE OF CODE
62 // "energies" AND 63 // "energies" AND "data" G4DataVector ARE ALLOCATED, FILLED IN AND NEVER USED OR
63 // DELETED. WHATSM 64 // DELETED. WHATSMORE LOADING FILE OPERATIONS WERE DONE BY G4ShellEMDataSet
64 // EVEN BEFORE THE 65 // EVEN BEFORE THE CHANGES I DID ON THIS FILE. SO THE FOLLOWING CODE IN MY
65 // OPINION SHOULD 66 // OPINION SHOULD BE USELESS AND SHOULD PRODUCE A MEMORY LEAK.
66 // 67 //
67 // 68 //
68 // ------------------------------------------- 69 // -------------------------------------------------------------------
69 70
70 #include "G4VCrossSectionHandler.hh" 71 #include "G4VCrossSectionHandler.hh"
71 #include "G4VDataSetAlgorithm.hh" 72 #include "G4VDataSetAlgorithm.hh"
72 #include "G4LogLogInterpolation.hh" 73 #include "G4LogLogInterpolation.hh"
73 #include "G4VEMDataSet.hh" 74 #include "G4VEMDataSet.hh"
74 #include "G4EMDataSet.hh" 75 #include "G4EMDataSet.hh"
75 #include "G4CompositeEMDataSet.hh" 76 #include "G4CompositeEMDataSet.hh"
76 #include "G4ShellEMDataSet.hh" 77 #include "G4ShellEMDataSet.hh"
77 #include "G4ProductionCutsTable.hh" 78 #include "G4ProductionCutsTable.hh"
78 #include "G4Material.hh" 79 #include "G4Material.hh"
79 #include "G4Element.hh" 80 #include "G4Element.hh"
80 #include "Randomize.hh" 81 #include "Randomize.hh"
81 #include <map> 82 #include <map>
82 #include <vector> 83 #include <vector>
83 #include <fstream> 84 #include <fstream>
84 #include <sstream> 85 #include <sstream>
85 86
86 //....oooOO0OOooo........oooOO0OOooo........oo <<
87 87
88 G4VCrossSectionHandler::G4VCrossSectionHandler 88 G4VCrossSectionHandler::G4VCrossSectionHandler()
89 { 89 {
90 crossSections = 0; 90 crossSections = 0;
91 interpolation = 0; 91 interpolation = 0;
92 Initialise(); 92 Initialise();
93 ActiveElements(); 93 ActiveElements();
94 } 94 }
95 95
96 //....oooOO0OOooo........oooOO0OOooo........oo <<
97 96
98 G4VCrossSectionHandler::G4VCrossSectionHandler 97 G4VCrossSectionHandler::G4VCrossSectionHandler(G4VDataSetAlgorithm* algorithm,
99 G4double minE, 98 G4double minE,
100 G4double maxE, 99 G4double maxE,
101 G4int bins, 100 G4int bins,
102 G4double unitE, 101 G4double unitE,
103 G4double unitData, 102 G4double unitData,
104 G4int minZ, 103 G4int minZ,
105 G4int maxZ) 104 G4int maxZ)
106 : interpolation(algorithm), eMin(minE), eMax << 105 : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
107 unit1(unitE), unit2(unitData), zMin(minZ), << 106 unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
108 nBins(bins) <<
109 { 107 {
110 crossSections = nullptr; << 108 crossSections = 0;
111 ActiveElements(); 109 ActiveElements();
112 } 110 }
113 111
114 //....oooOO0OOooo........oooOO0OOooo........oo <<
115 <<
116 G4VCrossSectionHandler::~G4VCrossSectionHandle 112 G4VCrossSectionHandler::~G4VCrossSectionHandler()
117 { 113 {
118 delete interpolation; 114 delete interpolation;
119 interpolation = nullptr; << 115 interpolation = 0;
>> 116 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
120 117
121 for (auto & pos : dataMap) << 118 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
122 { 119 {
123 G4VEMDataSet* dataSet = pos.second; << 120 // The following is a workaround for STL ObjectSpace implementation,
>> 121 // which does not support the standard and does not accept
>> 122 // the syntax pos->second
>> 123 // G4VEMDataSet* dataSet = pos->second;
>> 124 G4VEMDataSet* dataSet = (*pos).second;
124 delete dataSet; 125 delete dataSet;
125 } 126 }
126 127
127 if (crossSections != nullptr) << 128 if (crossSections != 0)
128 { 129 {
129 std::size_t n = crossSections->size(); << 130 size_t n = crossSections->size();
130 for (std::size_t i=0; i<n; i++) << 131 for (size_t i=0; i<n; i++)
131 { 132 {
132 delete (*crossSections)[i]; 133 delete (*crossSections)[i];
133 } 134 }
134 delete crossSections; 135 delete crossSections;
135 crossSections = nullptr; << 136 crossSections = 0;
136 } 137 }
137 } 138 }
138 139
139 //....oooOO0OOooo........oooOO0OOooo........oo <<
140 <<
141 void G4VCrossSectionHandler::Initialise(G4VDat 140 void G4VCrossSectionHandler::Initialise(G4VDataSetAlgorithm* algorithm,
142 G4double minE, G4double maxE, 141 G4double minE, G4double maxE,
143 G4int numberOfBins, 142 G4int numberOfBins,
144 G4double unitE, G4double unitData, 143 G4double unitE, G4double unitData,
145 G4int minZ, G4int maxZ) 144 G4int minZ, G4int maxZ)
146 { 145 {
147 if (algorithm != nullptr) << 146 if (algorithm != 0)
148 { 147 {
149 delete interpolation; 148 delete interpolation;
150 interpolation = algorithm; 149 interpolation = algorithm;
151 } 150 }
152 else 151 else
153 { 152 {
154 delete interpolation; 153 delete interpolation;
155 interpolation = CreateInterpolation(); 154 interpolation = CreateInterpolation();
156 } 155 }
157 156
158 eMin = minE; 157 eMin = minE;
159 eMax = maxE; 158 eMax = maxE;
160 nBins = numberOfBins; 159 nBins = numberOfBins;
161 unit1 = unitE; 160 unit1 = unitE;
162 unit2 = unitData; 161 unit2 = unitData;
163 zMin = minZ; 162 zMin = minZ;
164 zMax = maxZ; 163 zMax = maxZ;
165 } 164 }
166 165
167 //....oooOO0OOooo........oooOO0OOooo........oo <<
168 <<
169 void G4VCrossSectionHandler::PrintData() const 166 void G4VCrossSectionHandler::PrintData() const
170 { 167 {
171 for (auto & pos : dataMap) << 168 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
>> 169
>> 170 for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
172 { 171 {
173 G4int z = pos.first; << 172 // The following is a workaround for STL ObjectSpace implementation,
174 G4VEMDataSet* dataSet = pos.second; << 173 // which does not support the standard and does not accept
>> 174 // the syntax pos->first or pos->second
>> 175 // G4int z = pos->first;
>> 176 // G4VEMDataSet* dataSet = pos->second;
>> 177 G4int z = (*pos).first;
>> 178 G4VEMDataSet* dataSet = (*pos).second;
175 G4cout << "---- Data set for Z = " 179 G4cout << "---- Data set for Z = "
176 << z 180 << z
177 << G4endl; 181 << G4endl;
178 dataSet->PrintData(); 182 dataSet->PrintData();
179 G4cout << "----------------------------- 183 G4cout << "--------------------------------------------------" << G4endl;
180 } 184 }
181 } 185 }
182 186
183 //....oooOO0OOooo........oooOO0OOooo........oo <<
184 <<
185 void G4VCrossSectionHandler::LoadData(const G4 187 void G4VCrossSectionHandler::LoadData(const G4String& fileName)
186 { 188 {
187 std::size_t nZ = activeZ.size(); << 189 size_t nZ = activeZ.size();
188 for (std::size_t i=0; i<nZ; ++i) << 190 for (size_t i=0; i<nZ; i++)
189 { 191 {
190 G4int Z = G4int(activeZ[i]); << 192 G4int Z = (G4int) activeZ[i];
191 193
192 // Build the complete string identifying << 194 // Build the complete string identifying the file with the data set
193 const char* path = G4FindDataDir("G4LEDA << 195
>> 196 char* path = getenv("G4LEDATA");
194 if (!path) 197 if (!path)
195 { 198 {
196 G4Exception("G4VCrossSectionHandler: 199 G4Exception("G4VCrossSectionHandler::LoadData",
197 "em0006",FatalException,"G4LEDATA envi 200 "em0006",FatalException,"G4LEDATA environment variable not set");
198 return; 201 return;
199 } 202 }
200 203
201 std::ostringstream ost; 204 std::ostringstream ost;
202 ost << path << '/' << fileName << Z << " 205 ost << path << '/' << fileName << Z << ".dat";
203 std::ifstream file(ost.str().c_str()); 206 std::ifstream file(ost.str().c_str());
204 std::filebuf* lsdp = file.rdbuf(); 207 std::filebuf* lsdp = file.rdbuf();
205 208
206 if (! (lsdp->is_open()) ) 209 if (! (lsdp->is_open()) )
207 { 210 {
208 G4String excep = "data file: " + ost.str() 211 G4String excep = "data file: " + ost.str() + " not found";
209 G4Exception("G4VCrossSectionHandler: 212 G4Exception("G4VCrossSectionHandler::LoadData",
210 "em0003",FatalException,excep); 213 "em0003",FatalException,excep);
211 } 214 }
212 G4double a = 0; 215 G4double a = 0;
213 G4int k = 0; 216 G4int k = 0;
214 G4int nColumns = 2; 217 G4int nColumns = 2;
215 218
216 G4DataVector* orig_reg_energies = new G4 219 G4DataVector* orig_reg_energies = new G4DataVector;
217 G4DataVector* orig_reg_data = new G4Data 220 G4DataVector* orig_reg_data = new G4DataVector;
218 G4DataVector* log_reg_energies = new G4D 221 G4DataVector* log_reg_energies = new G4DataVector;
219 G4DataVector* log_reg_data = new G4DataV 222 G4DataVector* log_reg_data = new G4DataVector;
220 223
221 do 224 do
222 { 225 {
223 file >> a; 226 file >> a;
224 227
225 if (a==0.) a=1e-300; 228 if (a==0.) a=1e-300;
226 229
227 // The file is organized into four columns 230 // The file is organized into four columns:
228 // 1st column contains the values of energ 231 // 1st column contains the values of energy
229 // 2nd column contains the corresponding d 232 // 2nd column contains the corresponding data value
230 // The file terminates with the pattern: - 233 // The file terminates with the pattern: -1 -1
231 // - 234 // -2 -2
232 // 235 //
233 if (a != -1 && a != -2) 236 if (a != -1 && a != -2)
234 { 237 {
235 if (k%nColumns == 0) 238 if (k%nColumns == 0)
236 { 239 {
237 orig_reg_energies->push_back(a*unit1); 240 orig_reg_energies->push_back(a*unit1);
238 log_reg_energies->push_back(s 241 log_reg_energies->push_back(std::log10(a)+std::log10(unit1));
239 } 242 }
240 else if (k%nColumns == 1) 243 else if (k%nColumns == 1)
241 { 244 {
242 orig_reg_data->push_back(a*unit2); 245 orig_reg_data->push_back(a*unit2);
243 log_reg_data->push_back(std:: 246 log_reg_data->push_back(std::log10(a)+std::log10(unit2));
244 } 247 }
245 k++; 248 k++;
246 } 249 }
247 } 250 }
248 while (a != -2); // End of File 251 while (a != -2); // End of File
249 252
250 file.close(); 253 file.close();
251 G4VDataSetAlgorithm* algo = interpolatio 254 G4VDataSetAlgorithm* algo = interpolation->Clone();
252 G4VEMDataSet* dataSet = new G4EMDataSet( << 255
253 log_reg_energies,log_reg_data, << 256 G4VEMDataSet* dataSet = new G4EMDataSet(Z,orig_reg_energies,orig_reg_data,log_reg_energies,log_reg_data,algo);
>> 257
254 dataMap[Z] = dataSet; 258 dataMap[Z] = dataSet;
>> 259
255 } 260 }
256 } 261 }
257 262
258 //....oooOO0OOooo........oooOO0OOooo........oo <<
259 263
260 void G4VCrossSectionHandler::LoadNonLogData(co 264 void G4VCrossSectionHandler::LoadNonLogData(const G4String& fileName)
261 { 265 {
262 std::size_t nZ = activeZ.size(); << 266 size_t nZ = activeZ.size();
263 for (std::size_t i=0; i<nZ; ++i) << 267 for (size_t i=0; i<nZ; i++)
264 { 268 {
265 G4int Z = G4int(activeZ[i]); << 269 G4int Z = (G4int) activeZ[i];
266 270
267 // Build the complete string identifying << 271 // Build the complete string identifying the file with the data set
268 const char* path = G4FindDataDir("G4LEDA << 272
>> 273 char* path = getenv("G4LEDATA");
269 if (!path) 274 if (!path)
270 { 275 {
271 G4Exception("G4VCrossSectionHandler: 276 G4Exception("G4VCrossSectionHandler::LoadNonLogData",
272 "em0006",FatalException,"G4LEDATA envi 277 "em0006",FatalException,"G4LEDATA environment variable not set");
273 return; 278 return;
274 } 279 }
275 280
276 std::ostringstream ost; 281 std::ostringstream ost;
277 ost << path << '/' << fileName << Z << " 282 ost << path << '/' << fileName << Z << ".dat";
278 std::ifstream file(ost.str().c_str()); 283 std::ifstream file(ost.str().c_str());
279 std::filebuf* lsdp = file.rdbuf(); 284 std::filebuf* lsdp = file.rdbuf();
280 285
281 if (! (lsdp->is_open()) ) 286 if (! (lsdp->is_open()) )
282 { 287 {
283 G4String excep = "data file: " + ost.str() 288 G4String excep = "data file: " + ost.str() + " not found";
284 G4Exception("G4VCrossSectionHandler: 289 G4Exception("G4VCrossSectionHandler::LoadNonLogData",
285 "em0003",FatalException,excep); 290 "em0003",FatalException,excep);
286 } 291 }
287 G4double a = 0; 292 G4double a = 0;
288 G4int k = 0; 293 G4int k = 0;
289 G4int nColumns = 2; 294 G4int nColumns = 2;
290 295
291 G4DataVector* orig_reg_energies = new G4 296 G4DataVector* orig_reg_energies = new G4DataVector;
292 G4DataVector* orig_reg_data = new G4Data 297 G4DataVector* orig_reg_data = new G4DataVector;
293 298
294 do 299 do
295 { 300 {
296 file >> a; 301 file >> a;
297 302
298 // The file is organized into four columns 303 // The file is organized into four columns:
299 // 1st column contains the values of energ 304 // 1st column contains the values of energy
300 // 2nd column contains the corresponding d 305 // 2nd column contains the corresponding data value
301 // The file terminates with the pattern: - 306 // The file terminates with the pattern: -1 -1
302 // - 307 // -2 -2
303 // 308 //
304 if (a != -1 && a != -2) 309 if (a != -1 && a != -2)
305 { 310 {
306 if (k%nColumns == 0) 311 if (k%nColumns == 0)
307 { 312 {
308 orig_reg_energies->push_back(a*unit1); 313 orig_reg_energies->push_back(a*unit1);
309 } 314 }
310 else if (k%nColumns == 1) 315 else if (k%nColumns == 1)
311 { 316 {
312 orig_reg_data->push_back(a*unit2); 317 orig_reg_data->push_back(a*unit2);
313 } 318 }
314 k++; 319 k++;
315 } 320 }
316 } 321 }
317 while (a != -2); // End of File 322 while (a != -2); // End of File
318 323
319 file.close(); 324 file.close();
320 G4VDataSetAlgorithm* algo = interpolatio 325 G4VDataSetAlgorithm* algo = interpolation->Clone();
321 326
322 G4VEMDataSet* dataSet = new G4EMDataSet( 327 G4VEMDataSet* dataSet = new G4EMDataSet(Z,orig_reg_energies,orig_reg_data,algo);
>> 328
323 dataMap[Z] = dataSet; 329 dataMap[Z] = dataSet;
>> 330
324 } 331 }
325 } 332 }
326 333
327 //....oooOO0OOooo........oooOO0OOooo........oo <<
328 <<
329 void G4VCrossSectionHandler::LoadShellData(con 334 void G4VCrossSectionHandler::LoadShellData(const G4String& fileName)
330 { 335 {
331 std::size_t nZ = activeZ.size(); << 336 size_t nZ = activeZ.size();
332 for (std::size_t i=0; i<nZ; ++i) << 337 for (size_t i=0; i<nZ; i++)
333 { 338 {
334 G4int Z = G4int(activeZ[i]); << 339 G4int Z = (G4int) activeZ[i];
335 340
336 G4VDataSetAlgorithm* algo = interpolatio 341 G4VDataSetAlgorithm* algo = interpolation->Clone();
337 G4VEMDataSet* dataSet = new G4ShellEMDat 342 G4VEMDataSet* dataSet = new G4ShellEMDataSet(Z, algo);
338 dataSet->LoadData(fileName); << 343
>> 344 dataSet->LoadData(fileName);
>> 345
339 dataMap[Z] = dataSet; 346 dataMap[Z] = dataSet;
340 } 347 }
341 } 348 }
342 349
343 //....oooOO0OOooo........oooOO0OOooo........oo << 350
>> 351
344 352
345 void G4VCrossSectionHandler::Clear() 353 void G4VCrossSectionHandler::Clear()
346 { 354 {
347 // Reset the map of data sets: remove the da 355 // Reset the map of data sets: remove the data sets from the map
>> 356 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
>> 357
348 if(! dataMap.empty()) 358 if(! dataMap.empty())
349 { 359 {
350 for (auto & pos : dataMap) << 360 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
351 { 361 {
352 G4VEMDataSet* dataSet = pos.second; << 362 // The following is a workaround for STL ObjectSpace implementation,
>> 363 // which does not support the standard and does not accept
>> 364 // the syntax pos->first or pos->second
>> 365 // G4VEMDataSet* dataSet = pos->second;
>> 366 G4VEMDataSet* dataSet = (*pos).second;
353 delete dataSet; 367 delete dataSet;
354 dataSet = nullptr; << 368 dataSet = 0;
355 G4int i = pos.first; << 369 G4int i = (*pos).first;
356 dataMap[i] = nullptr; << 370 dataMap[i] = 0;
357 } 371 }
358 dataMap.clear(); 372 dataMap.clear();
359 } 373 }
>> 374
360 activeZ.clear(); 375 activeZ.clear();
361 ActiveElements(); 376 ActiveElements();
362 } 377 }
363 378
364 //....oooOO0OOooo........oooOO0OOooo........oo <<
365 <<
366 G4double G4VCrossSectionHandler::FindValue(G4i 379 G4double G4VCrossSectionHandler::FindValue(G4int Z, G4double energy) const
367 { 380 {
368 G4double value = 0.; 381 G4double value = 0.;
369 << 382
370 auto pos = dataMap.find(Z); << 383 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
>> 384 pos = dataMap.find(Z);
371 if (pos!= dataMap.end()) 385 if (pos!= dataMap.end())
372 { 386 {
>> 387 // The following is a workaround for STL ObjectSpace implementation,
>> 388 // which does not support the standard and does not accept
>> 389 // the syntax pos->first or pos->second
>> 390 // G4VEMDataSet* dataSet = pos->second;
373 G4VEMDataSet* dataSet = (*pos).second; 391 G4VEMDataSet* dataSet = (*pos).second;
374 value = dataSet->FindValue(energy); 392 value = dataSet->FindValue(energy);
375 } 393 }
376 else 394 else
377 { 395 {
378 G4cout << "WARNING: G4VCrossSectionHandl 396 G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find Z = "
379 << Z << G4endl; 397 << Z << G4endl;
380 } 398 }
381 return value; 399 return value;
382 } 400 }
383 401
384 //....oooOO0OOooo........oooOO0OOooo........oo <<
385 <<
386 G4double G4VCrossSectionHandler::FindValue(G4i 402 G4double G4VCrossSectionHandler::FindValue(G4int Z, G4double energy,
387 G4i 403 G4int shellIndex) const
388 { 404 {
389 G4double value = 0.; 405 G4double value = 0.;
390 auto pos = dataMap.find(Z); << 406
391 if (pos!= dataMap.cend()) << 407 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
>> 408 pos = dataMap.find(Z);
>> 409 if (pos!= dataMap.end())
392 { 410 {
>> 411 // The following is a workaround for STL ObjectSpace implementation,
>> 412 // which does not support the standard and does not accept
>> 413 // the syntax pos->first or pos->second
>> 414 // G4VEMDataSet* dataSet = pos->second;
393 G4VEMDataSet* dataSet = (*pos).second; 415 G4VEMDataSet* dataSet = (*pos).second;
394 if (shellIndex >= 0) 416 if (shellIndex >= 0)
395 { 417 {
396 G4int nComponents = (G4int)dataSet->Number << 418 G4int nComponents = dataSet->NumberOfComponents();
397 if(shellIndex < nComponents) 419 if(shellIndex < nComponents)
>> 420 // - MGP - Why doesn't it use G4VEMDataSet::FindValue directly?
398 value = dataSet->GetComponent(shellIndex 421 value = dataSet->GetComponent(shellIndex)->FindValue(energy);
399 else 422 else
400 { 423 {
401 G4cout << "WARNING: G4VCrossSectionHan 424 G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find"
402 << " shellIndex= " << shellIndex 425 << " shellIndex= " << shellIndex
403 << " for Z= " 426 << " for Z= "
404 << Z << G4endl; 427 << Z << G4endl;
405 } 428 }
406 } else { 429 } else {
407 value = dataSet->FindValue(energy); 430 value = dataSet->FindValue(energy);
408 } 431 }
409 } 432 }
410 else 433 else
411 { 434 {
412 G4cout << "WARNING: G4VCrossSectionHandl 435 G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find Z = "
413 << Z << G4endl; 436 << Z << G4endl;
414 } 437 }
415 return value; 438 return value;
416 } 439 }
417 440
418 //....oooOO0OOooo........oooOO0OOooo........oo <<
419 441
420 G4double G4VCrossSectionHandler::ValueForMater 442 G4double G4VCrossSectionHandler::ValueForMaterial(const G4Material* material,
421 G4double energy) const 443 G4double energy) const
422 { 444 {
423 G4double value = 0.; 445 G4double value = 0.;
>> 446
424 const G4ElementVector* elementVector = mater 447 const G4ElementVector* elementVector = material->GetElementVector();
425 const G4double* nAtomsPerVolume = material-> 448 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
426 std::size_t nElements = material->GetNumberO << 449 G4int nElements = material->GetNumberOfElements();
427 450
428 for (std::size_t i=0 ; i<nElements ; ++i) << 451 for (G4int i=0 ; i<nElements ; i++)
429 { 452 {
430 G4int Z = (G4int) (*elementVector)[i]->G 453 G4int Z = (G4int) (*elementVector)[i]->GetZ();
431 G4double elementValue = FindValue(Z,ener 454 G4double elementValue = FindValue(Z,energy);
432 G4double nAtomsVol = nAtomsPerVolume[i]; 455 G4double nAtomsVol = nAtomsPerVolume[i];
433 value += nAtomsVol * elementValue; 456 value += nAtomsVol * elementValue;
434 } 457 }
435 458
436 return value; 459 return value;
437 } 460 }
438 461
439 //....oooOO0OOooo........oooOO0OOooo........oo <<
440 462
441 G4VEMDataSet* G4VCrossSectionHandler::BuildMea 463 G4VEMDataSet* G4VCrossSectionHandler::BuildMeanFreePathForMaterials(const G4DataVector* energyCuts)
442 { 464 {
443 // Builds a CompositeDataSet containing the 465 // Builds a CompositeDataSet containing the mean free path for each material
444 // in the material table 466 // in the material table
>> 467
445 G4DataVector energyVector; 468 G4DataVector energyVector;
446 G4double dBin = std::log10(eMax/eMin) / nBin 469 G4double dBin = std::log10(eMax/eMin) / nBins;
447 470
448 for (G4int i=0; i<nBins+1; i++) 471 for (G4int i=0; i<nBins+1; i++)
449 { 472 {
450 energyVector.push_back(std::pow(10., std 473 energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
451 } 474 }
452 475
453 // Factory method to build cross sections in 476 // Factory method to build cross sections in derived classes,
454 // related to the type of physics process 477 // related to the type of physics process
455 478
456 if (crossSections != nullptr) << 479 if (crossSections != 0)
457 { // Reset the list of cross sections 480 { // Reset the list of cross sections
>> 481 std::vector<G4VEMDataSet*>::iterator mat;
458 if (! crossSections->empty()) 482 if (! crossSections->empty())
459 { 483 {
460 for (auto mat=crossSections->begin(); mat << 484 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
461 { 485 {
462 G4VEMDataSet* set = *mat; 486 G4VEMDataSet* set = *mat;
463 delete set; 487 delete set;
464 set = nullptr; << 488 set = 0;
465 } 489 }
466 crossSections->clear(); 490 crossSections->clear();
467 delete crossSections; 491 delete crossSections;
468 crossSections = nullptr; << 492 crossSections = 0;
469 } 493 }
470 } 494 }
471 495
472 crossSections = BuildCrossSectionsForMateria 496 crossSections = BuildCrossSectionsForMaterials(energyVector,energyCuts);
473 497
474 if (crossSections == nullptr) << 498 if (crossSections == 0)
475 { 499 {
476 G4Exception("G4VCrossSectionHandler::Bui 500 G4Exception("G4VCrossSectionHandler::BuildMeanFreePathForMaterials",
477 "em1010",FatalException,"crossSections 501 "em1010",FatalException,"crossSections = 0");
478 return 0; 502 return 0;
479 } 503 }
480 504
481 G4VDataSetAlgorithm* algo = CreateInterpolat 505 G4VDataSetAlgorithm* algo = CreateInterpolation();
482 G4VEMDataSet* materialSet = new G4CompositeE 506 G4VEMDataSet* materialSet = new G4CompositeEMDataSet(algo);
>> 507 //G4cout << "G4VCrossSectionHandler new dataset " << materialSet << G4endl;
483 508
484 G4DataVector* energies; 509 G4DataVector* energies;
485 G4DataVector* data; 510 G4DataVector* data;
486 G4DataVector* log_energies; 511 G4DataVector* log_energies;
487 G4DataVector* log_data; 512 G4DataVector* log_data;
>> 513
488 514
489 const G4ProductionCutsTable* theCoupleTable= 515 const G4ProductionCutsTable* theCoupleTable=
490 G4ProductionCutsTable::GetProductionCu 516 G4ProductionCutsTable::GetProductionCutsTable();
491 G4int numOfCouples = (G4int)theCoupleTable-> << 517 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 518
492 519
493 for (G4int mLocal=0; mLocal<numOfCouples; ++ << 520 for (size_t mLocal=0; mLocal<numOfCouples; mLocal++)
494 { 521 {
495 energies = new G4DataVector; 522 energies = new G4DataVector;
496 data = new G4DataVector; 523 data = new G4DataVector;
497 log_energies = new G4DataVector; 524 log_energies = new G4DataVector;
498 log_data = new G4DataVector; 525 log_data = new G4DataVector;
499 for (G4int bin=0; bin<nBins; bin++) 526 for (G4int bin=0; bin<nBins; bin++)
500 { 527 {
501 G4double energy = energyVector[bin]; 528 G4double energy = energyVector[bin];
502 energies->push_back(energy); 529 energies->push_back(energy);
503 log_energies->push_back(std::log10(e 530 log_energies->push_back(std::log10(energy));
504 G4VEMDataSet* matCrossSet = (*crossSection 531 G4VEMDataSet* matCrossSet = (*crossSections)[mLocal];
505 G4double materialCrossSection = 0.0; 532 G4double materialCrossSection = 0.0;
506 G4int nElm = (G4int)matCrossSet->Num << 533 G4int nElm = matCrossSet->NumberOfComponents();
507 for(G4int j=0; j<nElm; ++j) { << 534 for(G4int j=0; j<nElm; j++) {
508 materialCrossSection += matCrossSe 535 materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
509 } 536 }
510 537
511 if (materialCrossSection > 0.) 538 if (materialCrossSection > 0.)
512 { 539 {
513 data->push_back(1./materialCrossSectio 540 data->push_back(1./materialCrossSection);
514 log_data->push_back(std::log10(1 541 log_data->push_back(std::log10(1./materialCrossSection));
515 } 542 }
516 else 543 else
517 { 544 {
518 data->push_back(DBL_MAX); 545 data->push_back(DBL_MAX);
519 log_data->push_back(std::log10(D 546 log_data->push_back(std::log10(DBL_MAX));
520 } 547 }
521 } 548 }
522 G4VDataSetAlgorithm* algoLocal = CreateI 549 G4VDataSetAlgorithm* algoLocal = CreateInterpolation();
>> 550
>> 551 //G4VEMDataSet* dataSet = new G4EMDataSet(m,energies,data,algo,1.,1.);
>> 552
523 G4VEMDataSet* dataSet = new G4EMDataSet( 553 G4VEMDataSet* dataSet = new G4EMDataSet(mLocal,energies,data,log_energies,log_data,algoLocal,1.,1.);
>> 554
524 materialSet->AddComponent(dataSet); 555 materialSet->AddComponent(dataSet);
525 } 556 }
>> 557
526 return materialSet; 558 return materialSet;
527 } 559 }
528 560
529 //....oooOO0OOooo........oooOO0OOooo........oo <<
530 561
531 G4int G4VCrossSectionHandler::SelectRandomAtom 562 G4int G4VCrossSectionHandler::SelectRandomAtom(const G4MaterialCutsCouple* couple,
532 563 G4double e) const
533 { 564 {
534 // Select randomly an element within the mat 565 // Select randomly an element within the material, according to the weight
535 // determined by the cross sections in the d 566 // determined by the cross sections in the data set
>> 567
536 const G4Material* material = couple->GetMate 568 const G4Material* material = couple->GetMaterial();
537 G4int nElements = (G4int)material->GetNumber << 569 G4int nElements = material->GetNumberOfElements();
538 570
539 // Special case: the material consists of on 571 // Special case: the material consists of one element
540 if (nElements == 1) 572 if (nElements == 1)
541 { 573 {
542 G4int Z = (G4int) material->GetZ(); 574 G4int Z = (G4int) material->GetZ();
543 return Z; 575 return Z;
544 } 576 }
545 577
546 // Composite material 578 // Composite material
>> 579
547 const G4ElementVector* elementVector = mater 580 const G4ElementVector* elementVector = material->GetElementVector();
548 std::size_t materialIndex = couple->GetIndex << 581 size_t materialIndex = couple->GetIndex();
549 582
550 G4VEMDataSet* materialSet = (*crossSections) 583 G4VEMDataSet* materialSet = (*crossSections)[materialIndex];
551 G4double materialCrossSection0 = 0.0; 584 G4double materialCrossSection0 = 0.0;
552 G4DataVector cross; 585 G4DataVector cross;
553 cross.clear(); 586 cross.clear();
554 for ( G4int i=0; i < nElements; i++ ) 587 for ( G4int i=0; i < nElements; i++ )
555 { 588 {
556 G4double cr = materialSet->GetComponent( 589 G4double cr = materialSet->GetComponent(i)->FindValue(e);
557 materialCrossSection0 += cr; 590 materialCrossSection0 += cr;
558 cross.push_back(materialCrossSection0); 591 cross.push_back(materialCrossSection0);
559 } 592 }
560 593
561 G4double random = G4UniformRand() * material 594 G4double random = G4UniformRand() * materialCrossSection0;
>> 595
562 for (G4int k=0 ; k < nElements ; k++ ) 596 for (G4int k=0 ; k < nElements ; k++ )
563 { 597 {
564 if (random <= cross[k]) return (G4int) ( 598 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
565 } 599 }
566 // It should never get here 600 // It should never get here
567 return 0; 601 return 0;
568 } 602 }
569 603
570 //....oooOO0OOooo........oooOO0OOooo........oo <<
571 <<
572 const G4Element* G4VCrossSectionHandler::Selec 604 const G4Element* G4VCrossSectionHandler::SelectRandomElement(const G4MaterialCutsCouple* couple,
573 G4double e) const 605 G4double e) const
574 { 606 {
575 // Select randomly an element within the mat 607 // Select randomly an element within the material, according to the weight determined
576 // by the cross sections in the data set 608 // by the cross sections in the data set
>> 609
577 const G4Material* material = couple->GetMate 610 const G4Material* material = couple->GetMaterial();
578 G4Element* nullElement = 0; 611 G4Element* nullElement = 0;
579 G4int nElements = (G4int)material->GetNumber << 612 G4int nElements = material->GetNumberOfElements();
580 const G4ElementVector* elementVector = mater 613 const G4ElementVector* elementVector = material->GetElementVector();
581 614
582 // Special case: the material consists of on 615 // Special case: the material consists of one element
583 if (nElements == 1) 616 if (nElements == 1)
584 { 617 {
585 return (*elementVector)[0]; << 618 G4Element* element = (*elementVector)[0];
>> 619 return element;
586 } 620 }
587 else 621 else
588 { 622 {
589 // Composite material 623 // Composite material
590 624
591 std::size_t materialIndex = couple->GetI << 625 size_t materialIndex = couple->GetIndex();
592 626
593 G4VEMDataSet* materialSet = (*crossSecti 627 G4VEMDataSet* materialSet = (*crossSections)[materialIndex];
594 G4double materialCrossSection0 = 0.0; 628 G4double materialCrossSection0 = 0.0;
595 G4DataVector cross; 629 G4DataVector cross;
596 cross.clear(); 630 cross.clear();
597 for (G4int i=0; i<nElements; ++i) << 631 for (G4int i=0; i<nElements; i++)
598 { 632 {
599 G4double cr = materialSet->GetCompon 633 G4double cr = materialSet->GetComponent(i)->FindValue(e);
600 materialCrossSection0 += cr; 634 materialCrossSection0 += cr;
601 cross.push_back(materialCrossSection 635 cross.push_back(materialCrossSection0);
602 } 636 }
603 637
604 G4double random = G4UniformRand() * mate 638 G4double random = G4UniformRand() * materialCrossSection0;
605 639
606 for (G4int k=0 ; k < nElements ; ++k ) << 640 for (G4int k=0 ; k < nElements ; k++ )
607 { 641 {
608 if (random <= cross[k]) return (*ele 642 if (random <= cross[k]) return (*elementVector)[k];
609 } 643 }
610 // It should never end up here 644 // It should never end up here
611 G4cout << "G4VCrossSectionHandler::Selec 645 G4cout << "G4VCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
612 return nullElement; 646 return nullElement;
613 } 647 }
614 } 648 }
615 649
616 //....oooOO0OOooo........oooOO0OOooo........oo <<
617 <<
618 G4int G4VCrossSectionHandler::SelectRandomShel 650 G4int G4VCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
619 { 651 {
620 // Select randomly a shell, according to the 652 // Select randomly a shell, according to the weight determined by the cross sections
621 // in the data set 653 // in the data set
>> 654
622 // Note for later improvement: it would be u 655 // Note for later improvement: it would be useful to add a cache mechanism for already
623 // used shells to improve performance 656 // used shells to improve performance
>> 657
624 G4int shell = 0; 658 G4int shell = 0;
625 659
626 G4double totCrossSection = FindValue(Z,e); 660 G4double totCrossSection = FindValue(Z,e);
627 G4double random = G4UniformRand() * totCross 661 G4double random = G4UniformRand() * totCrossSection;
628 G4double partialSum = 0.; 662 G4double partialSum = 0.;
629 663
630 G4VEMDataSet* dataSet = nullptr; << 664 G4VEMDataSet* dataSet = 0;
631 auto pos = dataMap.find(Z); << 665 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
>> 666 pos = dataMap.find(Z);
>> 667 // The following is a workaround for STL ObjectSpace implementation,
>> 668 // which does not support the standard and does not accept
>> 669 // the syntax pos->first or pos->second
>> 670 // if (pos != dataMap.end()) dataSet = pos->second;
632 if (pos != dataMap.end()) 671 if (pos != dataMap.end())
633 dataSet = (*pos).second; 672 dataSet = (*pos).second;
634 else 673 else
635 { 674 {
636 G4Exception("G4VCrossSectionHandler::Sel 675 G4Exception("G4VCrossSectionHandler::SelectRandomShell",
637 "em1011",FatalException,"unable to loa 676 "em1011",FatalException,"unable to load the dataSet");
638 return 0; 677 return 0;
639 } 678 }
640 679
641 G4int nShells = (G4int)dataSet->NumberOfComp << 680 size_t nShells = dataSet->NumberOfComponents();
642 for (G4int i=0; i<nShells; ++i) << 681 for (size_t i=0; i<nShells; i++)
643 { 682 {
644 const G4VEMDataSet* shellDataSet = dataS 683 const G4VEMDataSet* shellDataSet = dataSet->GetComponent(i);
645 if (shellDataSet != nullptr) << 684 if (shellDataSet != 0)
646 { 685 {
647 G4double value = shellDataSet->FindValue(e 686 G4double value = shellDataSet->FindValue(e);
648 partialSum += value; 687 partialSum += value;
649 if (random <= partialSum) return i; 688 if (random <= partialSum) return i;
650 } 689 }
651 } 690 }
652 // It should never get here 691 // It should never get here
653 return shell; 692 return shell;
654 } 693 }
655 694
656 //....oooOO0OOooo........oooOO0OOooo........oo <<
657 <<
658 void G4VCrossSectionHandler::ActiveElements() 695 void G4VCrossSectionHandler::ActiveElements()
659 { 696 {
660 const G4MaterialTable* materialTable = G4Mat 697 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
661 if (materialTable == nullptr) << 698 if (materialTable == 0)
662 G4Exception("G4VCrossSectionHandler::Act 699 G4Exception("G4VCrossSectionHandler::ActiveElements",
663 "em1001",FatalException,"no MaterialTa 700 "em1001",FatalException,"no MaterialTable found");
664 701
665 std::size_t nMaterials = G4Material::GetNumb << 702 G4int nMaterials = G4Material::GetNumberOfMaterials();
666 703
667 for (std::size_t mLocal2=0; mLocal2<nMateria << 704 for (G4int mLocal2=0; mLocal2<nMaterials; mLocal2++)
668 { 705 {
669 const G4Material* material= (*materialTa 706 const G4Material* material= (*materialTable)[mLocal2];
670 const G4ElementVector* elementVector = m 707 const G4ElementVector* elementVector = material->GetElementVector();
671 const std::size_t nElements = material-> << 708 const G4int nElements = material->GetNumberOfElements();
672 709
673 for (std::size_t iEl=0; iEl<nElements; + << 710 for (G4int iEl=0; iEl<nElements; iEl++)
674 { 711 {
675 G4double Z = (*elementVector)[iEl]->GetZ() << 712 G4Element* element = (*elementVector)[iEl];
>> 713 G4double Z = element->GetZ();
676 if (!(activeZ.contains(Z)) && Z >= zMin && 714 if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
677 { 715 {
678 activeZ.push_back(Z); 716 activeZ.push_back(Z);
679 } 717 }
680 } 718 }
681 } 719 }
682 } 720 }
683 721
684 //....oooOO0OOooo........oooOO0OOooo........oo <<
685 <<
686 G4VDataSetAlgorithm* G4VCrossSectionHandler::C 722 G4VDataSetAlgorithm* G4VCrossSectionHandler::CreateInterpolation()
687 { 723 {
688 G4VDataSetAlgorithm* algorithm = new G4LogLo 724 G4VDataSetAlgorithm* algorithm = new G4LogLogInterpolation;
689 return algorithm; 725 return algorithm;
690 } 726 }
691 727
692 //....oooOO0OOooo........oooOO0OOooo........oo <<
693 <<
694 G4int G4VCrossSectionHandler::NumberOfComponen 728 G4int G4VCrossSectionHandler::NumberOfComponents(G4int Z) const
695 { 729 {
696 G4int n = 0; 730 G4int n = 0;
697 731
698 auto pos = dataMap.find(Z); << 732 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
>> 733 pos = dataMap.find(Z);
699 if (pos!= dataMap.end()) 734 if (pos!= dataMap.end())
700 { 735 {
701 G4VEMDataSet* dataSet = (*pos).second; 736 G4VEMDataSet* dataSet = (*pos).second;
702 n = (G4int)dataSet->NumberOfComponents() << 737 n = dataSet->NumberOfComponents();
703 } 738 }
704 else 739 else
705 { 740 {
706 G4cout << "WARNING: G4VCrossSectionHandl 741 G4cout << "WARNING: G4VCrossSectionHandler::NumberOfComponents did not "
707 << "find Z = " 742 << "find Z = "
708 << Z << G4endl; 743 << Z << G4endl;
709 } 744 }
710 return n; 745 return n;
711 } 746 }
712 747
713 748
714 749