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