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