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