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