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