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