| 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.8 2001/10/10 16:46:06 pia Exp $
>> 25 // GEANT4 tag $Name: geant4-04-00 $
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(
>> 409 const G4DataVector* energyCuts)
442 { 410 {
443 // Builds a CompositeDataSet containing the 411 // Builds a CompositeDataSet containing the mean free path for each material
444 // in the material table 412 // in the material table
>> 413
445 G4DataVector energyVector; 414 G4DataVector energyVector;
446 G4double dBin = std::log10(eMax/eMin) / nBin << 415 G4double dBin = log10(eMax/eMin) / nBins;
447 416
448 for (G4int i=0; i<nBins+1; i++) << 417 for (G4int i=0; i<nBins+1; i++)
449 { 418 {
450 energyVector.push_back(std::pow(10., std << 419 energyVector.push_back(pow(10., log10(eMin)+i*dBin));
451 } 420 }
452 421
453 // Factory method to build cross sections in << 422 // Factory method to build cross sections in derived classes,
454 // related to the type of physics process << 423 // related to the type of physics process
455 424
456 if (crossSections != nullptr) << 425 if (crossSections != 0)
457 { // Reset the list of cross sections 426 { // Reset the list of cross sections
>> 427 G4std::vector<G4VEMDataSet*>::iterator mat;
458 if (! crossSections->empty()) 428 if (! crossSections->empty())
459 { 429 {
460 for (auto mat=crossSections->begin(); mat << 430 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
461 { 431 {
462 G4VEMDataSet* set = *mat; 432 G4VEMDataSet* set = *mat;
463 delete set; 433 delete set;
464 set = nullptr; << 434 set = 0;
465 } 435 }
466 crossSections->clear(); 436 crossSections->clear();
467 delete crossSections; 437 delete crossSections;
468 crossSections = nullptr; << 438 crossSections = 0;
469 } 439 }
470 } 440 }
471 441
472 crossSections = BuildCrossSectionsForMateria 442 crossSections = BuildCrossSectionsForMaterials(energyVector,energyCuts);
473 443
474 if (crossSections == nullptr) << 444 if (crossSections == 0)
475 { << 445 G4Exception("G4VCrossSectionHandler::BuildMeanFreePathForMaterials, crossSections = 0");
476 G4Exception("G4VCrossSectionHandler::Bui <<
477 "em1010",FatalException,"crossSections <<
478 return 0; <<
479 } <<
480 446
481 G4VDataSetAlgorithm* algo = CreateInterpolat 447 G4VDataSetAlgorithm* algo = CreateInterpolation();
482 G4VEMDataSet* materialSet = new G4CompositeE 448 G4VEMDataSet* materialSet = new G4CompositeEMDataSet(algo);
483 449
484 G4DataVector* energies; 450 G4DataVector* energies;
485 G4DataVector* data; 451 G4DataVector* data;
486 G4DataVector* log_energies; <<
487 G4DataVector* log_data; <<
488 <<
489 const G4ProductionCutsTable* theCoupleTable= <<
490 G4ProductionCutsTable::GetProductionCu <<
491 G4int numOfCouples = (G4int)theCoupleTable-> <<
492 452
493 for (G4int mLocal=0; mLocal<numOfCouples; ++ << 453 size_t nMaterials = G4Material::GetNumberOfMaterials();
>> 454
>> 455 for (size_t m=0; m<nMaterials; m++)
494 { 456 {
495 energies = new G4DataVector; << 457 energies = new G4DataVector;
496 data = new G4DataVector; 458 data = new G4DataVector;
497 log_energies = new G4DataVector; <<
498 log_data = new G4DataVector; <<
499 for (G4int bin=0; bin<nBins; bin++) 459 for (G4int bin=0; bin<nBins; bin++)
500 { 460 {
501 G4double energy = energyVector[bin]; 461 G4double energy = energyVector[bin];
502 energies->push_back(energy); 462 energies->push_back(energy);
503 log_energies->push_back(std::log10(e << 463 G4VEMDataSet* matCrossSet = (*crossSections)[m];
504 G4VEMDataSet* matCrossSet = (*crossSection << 464 G4double materialCrossSection = matCrossSet->FindValue(energy);
505 G4double materialCrossSection = 0.0; << 465
506 G4int nElm = (G4int)matCrossSet->Num <<
507 for(G4int j=0; j<nElm; ++j) { <<
508 materialCrossSection += matCrossSe <<
509 } <<
510 <<
511 if (materialCrossSection > 0.) 466 if (materialCrossSection > 0.)
512 { 467 {
513 data->push_back(1./materialCrossSectio 468 data->push_back(1./materialCrossSection);
514 log_data->push_back(std::log10(1 <<
515 } 469 }
516 else 470 else
517 { 471 {
518 data->push_back(DBL_MAX); 472 data->push_back(DBL_MAX);
519 log_data->push_back(std::log10(D <<
520 } 473 }
521 } 474 }
522 G4VDataSetAlgorithm* algoLocal = CreateI << 475 G4VDataSetAlgorithm* algo = CreateInterpolation();
523 G4VEMDataSet* dataSet = new G4EMDataSet( << 476 G4VEMDataSet* dataSet = new G4EMDataSet(m,energies,data,algo,1.,1.);
524 materialSet->AddComponent(dataSet); 477 materialSet->AddComponent(dataSet);
525 } << 478 }
>> 479
526 return materialSet; 480 return materialSet;
527 } 481 }
528 482
529 //....oooOO0OOooo........oooOO0OOooo........oo << 483 G4int G4VCrossSectionHandler::SelectRandomAtom(const G4Material* material, G4double e) const
530 <<
531 G4int G4VCrossSectionHandler::SelectRandomAtom <<
532 <<
533 { 484 {
534 // Select randomly an element within the mat << 485 // Select randomly an element within the material, according to the weight
535 // determined by the cross sections in the d 486 // determined by the cross sections in the data set
536 const G4Material* material = couple->GetMate << 487
537 G4int nElements = (G4int)material->GetNumber << 488 G4int nElements = material->GetNumberOfElements();
>> 489 const G4ElementVector* elementVector = material->GetElementVector();
538 490
539 // Special case: the material consists of on 491 // Special case: the material consists of one element
540 if (nElements == 1) << 492 if (nElements == 1)
541 { 493 {
542 G4int Z = (G4int) material->GetZ(); << 494 G4int Z = (G4int) (*elementVector)[0]->GetZ();
543 return Z; 495 return Z;
544 } 496 }
545 497
546 // Composite material 498 // Composite material
547 const G4ElementVector* elementVector = mater <<
548 std::size_t materialIndex = couple->GetIndex <<
549 499
550 G4VEMDataSet* materialSet = (*crossSections) << 500 G4double materialCrossSection0 = ValueForMaterial(material,e);
551 G4double materialCrossSection0 = 0.0; << 501 // size_t materialIndex = material->GetIndex();
552 G4DataVector cross; << 502
553 cross.clear(); << 503 // G4VEMDataSet* materialSet = crossSections[materialIndex];
554 for ( G4int i=0; i < nElements; i++ ) << 504 // G4double materialCrossSection = materialSet->FindValue(e);
555 { <<
556 G4double cr = materialSet->GetComponent( <<
557 materialCrossSection0 += cr; <<
558 cross.push_back(materialCrossSection0); <<
559 } <<
560 505
561 G4double random = G4UniformRand() * material 506 G4double random = G4UniformRand() * materialCrossSection0;
562 for (G4int k=0 ; k < nElements ; k++ ) << 507 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
563 { << 508 G4double partialSumSigma = 0.;
564 if (random <= cross[k]) return (G4int) ( << 509
>> 510 for ( G4int i=0 ; i < nElements ; i++ )
>> 511 {
>> 512 G4int Z = (G4int) (*elementVector)[i]->GetZ();
>> 513 G4double crossSection = FindValue(Z,e);
>> 514 partialSumSigma += nAtomsPerVolume[i] * crossSection;
>> 515 if (random <= partialSumSigma) return Z;
565 } 516 }
566 // It should never get here 517 // It should never get here
567 return 0; 518 return 0;
568 } 519 }
569 520
570 //....oooOO0OOooo........oooOO0OOooo........oo << 521 const G4Element* G4VCrossSectionHandler::SelectRandomElement(const G4Material* material,
571 <<
572 const G4Element* G4VCrossSectionHandler::Selec <<
573 G4double e) const 522 G4double e) const
574 { 523 {
575 // Select randomly an element within the mat 524 // Select randomly an element within the material, according to the weight determined
576 // by the cross sections in the data set 525 // by the cross sections in the data set
577 const G4Material* material = couple->GetMate << 526
578 G4Element* nullElement = 0; 527 G4Element* nullElement = 0;
579 G4int nElements = (G4int)material->GetNumber << 528 G4int nElements = material->GetNumberOfElements();
580 const G4ElementVector* elementVector = mater 529 const G4ElementVector* elementVector = material->GetElementVector();
581 530
582 // Special case: the material consists of on 531 // Special case: the material consists of one element
583 if (nElements == 1) << 532 if (nElements == 1)
584 { 533 {
585 return (*elementVector)[0]; << 534 G4Element* element = (*elementVector)[0];
>> 535 return element;
586 } 536 }
587 else 537 else
588 { 538 {
589 // Composite material << 539 // Composite material
590 << 540 G4double materialCrossSection0 = ValueForMaterial(material,e);
591 std::size_t materialIndex = couple->GetI << 541 // size_t materialIndex = material->GetIndex();
592 << 542 // G4VEMDataSet* materialSet = crossSections[materialIndex];
593 G4VEMDataSet* materialSet = (*crossSecti << 543 // G4double materialCrossSection = materialSet->FindValue(e);
594 G4double materialCrossSection0 = 0.0; << 544
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 545 G4double random = G4UniformRand() * materialCrossSection0;
605 << 546 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
606 for (G4int k=0 ; k < nElements ; ++k ) << 547 G4double partialSumSigma = 0.;
607 { << 548
608 if (random <= cross[k]) return (*ele << 549 for ( G4int i=0 ; i < nElements ; i++ )
609 } << 550 {
>> 551 G4Element* element = (*elementVector)[i];
>> 552 G4int Z = (G4int) element->GetZ();
>> 553 G4double crossSection = FindValue(Z,e);
>> 554 partialSumSigma += nAtomsPerVolume[i] * crossSection;
>> 555 if (random <= partialSumSigma) return element;
>> 556 }
>> 557 }
610 // It should never end up here 558 // It should never end up here
611 G4cout << "G4VCrossSectionHandler::Selec 559 G4cout << "G4VCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
612 return nullElement; 560 return nullElement;
613 } <<
614 } 561 }
615 562
616 //....oooOO0OOooo........oooOO0OOooo........oo <<
617 <<
618 G4int G4VCrossSectionHandler::SelectRandomShel 563 G4int G4VCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
619 { 564 {
620 // Select randomly a shell, according to the 565 // Select randomly a shell, according to the weight determined by the cross sections
621 // in the data set 566 // in the data set
>> 567
622 // Note for later improvement: it would be u 568 // Note for later improvement: it would be useful to add a cache mechanism for already
623 // used shells to improve performance 569 // used shells to improve performance
>> 570
624 G4int shell = 0; 571 G4int shell = 0;
625 572
626 G4double totCrossSection = FindValue(Z,e); 573 G4double totCrossSection = FindValue(Z,e);
627 G4double random = G4UniformRand() * totCross 574 G4double random = G4UniformRand() * totCrossSection;
628 G4double partialSum = 0.; 575 G4double partialSum = 0.;
629 576
630 G4VEMDataSet* dataSet = nullptr; << 577 G4VEMDataSet* dataSet = 0;
631 auto pos = dataMap.find(Z); << 578 G4std::map<G4int,G4VEMDataSet*,G4std::less<G4int> >::const_iterator pos;
632 if (pos != dataMap.end()) << 579 pos = dataMap.find(Z);
633 dataSet = (*pos).second; << 580 // The following is a workaround for STL ObjectSpace implementation,
634 else << 581 // which does not support the standard and does not accept
635 { << 582 // the syntax pos->first or pos->second
636 G4Exception("G4VCrossSectionHandler::Sel << 583 // if (pos != dataMap.end()) dataSet = pos->second;
637 "em1011",FatalException,"unable to loa << 584 if (pos != dataMap.end()) dataSet = (*pos).second;
638 return 0; <<
639 } <<
640 585
641 G4int nShells = (G4int)dataSet->NumberOfComp << 586 size_t nShells = dataSet->NumberOfComponents();
642 for (G4int i=0; i<nShells; ++i) << 587 for (size_t i=0; i<nShells; i++)
643 { 588 {
644 const G4VEMDataSet* shellDataSet = dataS 589 const G4VEMDataSet* shellDataSet = dataSet->GetComponent(i);
645 if (shellDataSet != nullptr) << 590 if (shellDataSet != 0)
646 { 591 {
647 G4double value = shellDataSet->FindValue(e 592 G4double value = shellDataSet->FindValue(e);
648 partialSum += value; 593 partialSum += value;
649 if (random <= partialSum) return i; 594 if (random <= partialSum) return i;
650 } 595 }
651 } 596 }
652 // It should never get here 597 // It should never get here
653 return shell; 598 return shell;
654 } 599 }
655 600
656 //....oooOO0OOooo........oooOO0OOooo........oo <<
657 <<
658 void G4VCrossSectionHandler::ActiveElements() 601 void G4VCrossSectionHandler::ActiveElements()
659 { 602 {
660 const G4MaterialTable* materialTable = G4Mat 603 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
661 if (materialTable == nullptr) << 604 if (materialTable == 0)
662 G4Exception("G4VCrossSectionHandler::Act << 605 G4Exception("G4VCrossSectionHandler::ActiveElements - no MaterialTable found)");
663 "em1001",FatalException,"no MaterialTa << 606
664 << 607 G4int nMaterials = G4Material::GetNumberOfMaterials();
665 std::size_t nMaterials = G4Material::GetNumb << 608
666 << 609 for (G4int m=0; m<nMaterials; m++)
667 for (std::size_t mLocal2=0; mLocal2<nMateria <<
668 { 610 {
669 const G4Material* material= (*materialTa << 611 const G4Material* material= (*materialTable)[m];
670 const G4ElementVector* elementVector = m 612 const G4ElementVector* elementVector = material->GetElementVector();
671 const std::size_t nElements = material-> << 613 const G4int nElements = material->GetNumberOfElements();
672 << 614
673 for (std::size_t iEl=0; iEl<nElements; + << 615 for (G4int iEl=0; iEl<nElements; iEl++)
674 { 616 {
675 G4double Z = (*elementVector)[iEl]->GetZ() << 617 G4Element* element = (*elementVector)[iEl];
676 if (!(activeZ.contains(Z)) && Z >= zMin && << 618 G4double Z = element->GetZ();
>> 619 if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
677 { 620 {
678 activeZ.push_back(Z); 621 activeZ.push_back(Z);
679 } 622 }
680 } 623 }
681 } 624 }
682 } 625 }
683 626
684 //....oooOO0OOooo........oooOO0OOooo........oo <<
685 <<
686 G4VDataSetAlgorithm* G4VCrossSectionHandler::C 627 G4VDataSetAlgorithm* G4VCrossSectionHandler::CreateInterpolation()
687 { 628 {
688 G4VDataSetAlgorithm* algorithm = new G4LogLo 629 G4VDataSetAlgorithm* algorithm = new G4LogLogInterpolation;
689 return algorithm; 630 return algorithm;
690 } 631 }
691 632
692 //....oooOO0OOooo........oooOO0OOooo........oo <<
693 <<
694 G4int G4VCrossSectionHandler::NumberOfComponen 633 G4int G4VCrossSectionHandler::NumberOfComponents(G4int Z) const
695 { 634 {
696 G4int n = 0; 635 G4int n = 0;
697 << 636
698 auto pos = dataMap.find(Z); << 637 G4std::map<G4int,G4VEMDataSet*,G4std::less<G4int> >::const_iterator pos;
>> 638 pos = dataMap.find(Z);
699 if (pos!= dataMap.end()) 639 if (pos!= dataMap.end())
700 { 640 {
701 G4VEMDataSet* dataSet = (*pos).second; 641 G4VEMDataSet* dataSet = (*pos).second;
702 n = (G4int)dataSet->NumberOfComponents() << 642 n = dataSet->NumberOfComponents();
703 } 643 }
704 else 644 else
705 { 645 {
706 G4cout << "WARNING: G4VCrossSectionHandl 646 G4cout << "WARNING: G4VCrossSectionHandler::NumberOfComponents did not "
707 << "find Z = " 647 << "find Z = "
708 << Z << G4endl; 648 << Z << G4endl;
709 } 649 }
710 return n; 650 return n;
711 } 651 }
712 652
713 653
714 654