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