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