| 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: G4PixeCrossSectionHandler.cc 70904 2013-06-07 10:34:25Z gcosmo $
27 // 28 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 29 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // 30 //
30 // History: 31 // History:
31 // ----------- 32 // -----------
32 // 16 Jun 2008 MGP Created; Cross section ma 33 // 16 Jun 2008 MGP Created; Cross section manager for hadron impact ionization
33 // Documented in: 34 // Documented in:
34 // M.G. Pia et al., PIXE Sim 35 // M.G. Pia et al., PIXE Simulation With Geant4,
35 // IEEE Trans. Nucl. Sci., v 36 // IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 3614-3649, Dec. 2009
36 // 37 //
37 // ------------------------------------------- 38 // -------------------------------------------------------------------
38 39
39 #include "G4PixeCrossSectionHandler.hh" 40 #include "G4PixeCrossSectionHandler.hh"
40 #include "G4PhysicalConstants.hh" 41 #include "G4PhysicalConstants.hh"
41 #include "G4IInterpolator.hh" 42 #include "G4IInterpolator.hh"
42 #include "G4LogLogInterpolator.hh" 43 #include "G4LogLogInterpolator.hh"
43 #include "G4IDataSet.hh" 44 #include "G4IDataSet.hh"
44 #include "G4DataSet.hh" 45 #include "G4DataSet.hh"
45 #include "G4CompositeDataSet.hh" 46 #include "G4CompositeDataSet.hh"
46 #include "G4PixeShellDataSet.hh" 47 #include "G4PixeShellDataSet.hh"
47 #include "G4ProductionCutsTable.hh" 48 #include "G4ProductionCutsTable.hh"
48 #include "G4Material.hh" 49 #include "G4Material.hh"
49 #include "G4Element.hh" 50 #include "G4Element.hh"
50 #include "Randomize.hh" 51 #include "Randomize.hh"
51 #include "G4SystemOfUnits.hh" 52 #include "G4SystemOfUnits.hh"
52 #include "G4ParticleDefinition.hh" 53 #include "G4ParticleDefinition.hh"
53 54
54 #include <map> 55 #include <map>
55 #include <vector> 56 #include <vector>
56 #include <fstream> 57 #include <fstream>
57 #include <sstream> 58 #include <sstream>
58 #include <memory> <<
59 59
60 60
61 G4PixeCrossSectionHandler::G4PixeCrossSectionH 61 G4PixeCrossSectionHandler::G4PixeCrossSectionHandler()
62 { 62 {
63 crossSections = 0; 63 crossSections = 0;
64 interpolation = 0; 64 interpolation = 0;
65 // Initialise with default values 65 // Initialise with default values
66 Initialise(0,"","","",1.*keV,0.1*GeV,200,MeV 66 Initialise(0,"","","",1.*keV,0.1*GeV,200,MeV,barn,6,92);
67 ActiveElements(); 67 ActiveElements();
68 } 68 }
69 69
70 70
71 G4PixeCrossSectionHandler::G4PixeCrossSectionH 71 G4PixeCrossSectionHandler::G4PixeCrossSectionHandler(G4IInterpolator* algorithm,
72 const G4String& modelK, 72 const G4String& modelK,
73 const G4String& modelL, 73 const G4String& modelL,
74 const G4String& modelM, 74 const G4String& modelM,
75 G4double minE, 75 G4double minE,
76 G4double maxE, 76 G4double maxE,
77 G4int bins, 77 G4int bins,
78 G4double unitE, 78 G4double unitE,
79 G4double unitData, 79 G4double unitData,
80 G4int minZ, 80 G4int minZ,
81 G4int maxZ) 81 G4int maxZ)
82 : interpolation(algorithm), eMin(minE), eMax 82 : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
83 unit1(unitE), unit2(unitData), zMin(minZ), 83 unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
84 { 84 {
85 crossSections = 0; 85 crossSections = 0;
86 86
87 crossModel.push_back(modelK); 87 crossModel.push_back(modelK);
88 crossModel.push_back(modelL); 88 crossModel.push_back(modelL);
89 crossModel.push_back(modelM); 89 crossModel.push_back(modelM);
90 90
91 //std::cout << "PixeCrossSectionHandler cons 91 //std::cout << "PixeCrossSectionHandler constructor - crossModel[0] = "
92 // << crossModel[0] 92 // << crossModel[0]
93 // << std::endl; 93 // << std::endl;
94 94
95 ActiveElements(); 95 ActiveElements();
96 } 96 }
97 97
98 G4PixeCrossSectionHandler::~G4PixeCrossSection 98 G4PixeCrossSectionHandler::~G4PixeCrossSectionHandler()
99 { 99 {
100 delete interpolation; 100 delete interpolation;
101 interpolation = 0; 101 interpolation = 0;
102 std::map<G4int,G4IDataSet*,std::less<G4int> 102 std::map<G4int,G4IDataSet*,std::less<G4int> >::iterator pos;
103 103
104 for (pos = dataMap.begin(); pos != dataMap.e 104 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
105 { 105 {
106 // The following is a workaround for STL 106 // The following is a workaround for STL ObjectSpace implementation,
107 // which does not support the standard a 107 // which does not support the standard and does not accept
108 // the syntax pos->second 108 // the syntax pos->second
109 // G4IDataSet* dataSet = pos->second; 109 // G4IDataSet* dataSet = pos->second;
110 G4IDataSet* dataSet = (*pos).second; 110 G4IDataSet* dataSet = (*pos).second;
111 delete dataSet; 111 delete dataSet;
112 } 112 }
113 113
114 if (crossSections != 0) 114 if (crossSections != 0)
115 { 115 {
116 std::size_t n = crossSections->size(); << 116 size_t n = crossSections->size();
117 for (std::size_t i=0; i<n; ++i) << 117 for (size_t i=0; i<n; i++)
118 { 118 {
119 delete (*crossSections)[i]; 119 delete (*crossSections)[i];
120 } 120 }
121 delete crossSections; 121 delete crossSections;
122 crossSections = 0; 122 crossSections = 0;
123 } 123 }
124 } 124 }
125 125
126 void G4PixeCrossSectionHandler::Initialise(G4I 126 void G4PixeCrossSectionHandler::Initialise(G4IInterpolator* algorithm,
127 const G4String& modelK, 127 const G4String& modelK,
128 const G4String& modelL, 128 const G4String& modelL,
129 const G4String& modelM, 129 const G4String& modelM,
130 G4double minE, G4double maxE, 130 G4double minE, G4double maxE,
131 G4int numberOfBins, 131 G4int numberOfBins,
132 G4double unitE, G4double unitData 132 G4double unitE, G4double unitData,
133 G4int minZ, G4int maxZ) 133 G4int minZ, G4int maxZ)
134 { 134 {
135 if (algorithm != 0) 135 if (algorithm != 0)
136 { 136 {
137 delete interpolation; 137 delete interpolation;
138 interpolation = algorithm; 138 interpolation = algorithm;
139 } 139 }
140 else 140 else
141 { 141 {
142 interpolation = CreateInterpolation(); 142 interpolation = CreateInterpolation();
143 } 143 }
144 144
145 eMin = minE; 145 eMin = minE;
146 eMax = maxE; 146 eMax = maxE;
147 nBins = numberOfBins; 147 nBins = numberOfBins;
148 unit1 = unitE; 148 unit1 = unitE;
149 unit2 = unitData; 149 unit2 = unitData;
150 zMin = minZ; 150 zMin = minZ;
151 zMax = maxZ; 151 zMax = maxZ;
152 152
153 crossModel.push_back(modelK); 153 crossModel.push_back(modelK);
154 crossModel.push_back(modelL); 154 crossModel.push_back(modelL);
155 crossModel.push_back(modelM); 155 crossModel.push_back(modelM);
156 156
157 } 157 }
158 158
159 void G4PixeCrossSectionHandler::PrintData() co 159 void G4PixeCrossSectionHandler::PrintData() const
160 { 160 {
161 std::map<G4int,G4IDataSet*,std::less<G4int> 161 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
162 162
163 for (pos = dataMap.begin(); pos != dataMap.e 163 for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
164 { 164 {
165 // The following is a workaround for STL 165 // The following is a workaround for STL ObjectSpace implementation,
166 // which does not support the standard a 166 // which does not support the standard and does not accept
167 // the syntax pos->first or pos->second 167 // the syntax pos->first or pos->second
168 // G4int z = pos->first; 168 // G4int z = pos->first;
169 // G4IDataSet* dataSet = pos->second; 169 // G4IDataSet* dataSet = pos->second;
170 G4int z = (*pos).first; 170 G4int z = (*pos).first;
171 G4IDataSet* dataSet = (*pos).second; 171 G4IDataSet* dataSet = (*pos).second;
172 G4cout << "---- Data set for Z = " 172 G4cout << "---- Data set for Z = "
173 << z 173 << z
174 << G4endl; 174 << G4endl;
175 dataSet->PrintData(); 175 dataSet->PrintData();
176 G4cout << "----------------------------- 176 G4cout << "--------------------------------------------------" << G4endl;
177 } 177 }
178 } 178 }
179 179
180 void G4PixeCrossSectionHandler::LoadShellData( 180 void G4PixeCrossSectionHandler::LoadShellData(const G4String& fileName)
181 { 181 {
182 std::size_t nZ = activeZ.size(); << 182 size_t nZ = activeZ.size();
183 for (std::size_t i=0; i<nZ; ++i) << 183 for (size_t i=0; i<nZ; i++)
184 { 184 {
185 G4int Z = (G4int) activeZ[i]; 185 G4int Z = (G4int) activeZ[i];
186 G4IInterpolator* algo = interpolation->C 186 G4IInterpolator* algo = interpolation->Clone();
187 G4IDataSet* dataSet = new G4PixeShellDat 187 G4IDataSet* dataSet = new G4PixeShellDataSet(Z, algo,crossModel[0],crossModel[1],crossModel[2]);
188 188
189 // Degug printing 189 // Degug printing
190 //std::cout << "PixeCrossSectionHandler: 190 //std::cout << "PixeCrossSectionHandler::Load - "
191 // << Z 191 // << Z
192 // << ", modelK = " 192 // << ", modelK = "
193 // << crossModel[0] 193 // << crossModel[0]
194 // << " fileName = " 194 // << " fileName = "
195 // << fileName 195 // << fileName
196 // << std::endl; 196 // << std::endl;
197 197
198 dataSet->LoadData(fileName); 198 dataSet->LoadData(fileName);
199 dataMap[Z] = dataSet; 199 dataMap[Z] = dataSet;
200 } 200 }
201 201
202 // Build cross sections for materials if not 202 // Build cross sections for materials if not already built
203 if (! crossSections) 203 if (! crossSections)
204 { 204 {
205 BuildForMaterials(); 205 BuildForMaterials();
206 } 206 }
207 207
208 } 208 }
209 209
210 void G4PixeCrossSectionHandler::Clear() 210 void G4PixeCrossSectionHandler::Clear()
211 { 211 {
212 // Reset the map of data sets: remove the da 212 // Reset the map of data sets: remove the data sets from the map
213 std::map<G4int,G4IDataSet*,std::less<G4int> 213 std::map<G4int,G4IDataSet*,std::less<G4int> >::iterator pos;
214 214
215 if(! dataMap.empty()) 215 if(! dataMap.empty())
216 { 216 {
217 for (pos = dataMap.begin(); pos != dataM 217 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
218 { 218 {
219 // The following is a workaround for STL O 219 // The following is a workaround for STL ObjectSpace implementation,
220 // which does not support the standard and 220 // which does not support the standard and does not accept
221 // the syntax pos->first or pos->second 221 // the syntax pos->first or pos->second
222 // G4IDataSet* dataSet = pos->second; 222 // G4IDataSet* dataSet = pos->second;
223 G4IDataSet* dataSet = (*pos).second; 223 G4IDataSet* dataSet = (*pos).second;
224 delete dataSet; 224 delete dataSet;
225 dataSet = 0; 225 dataSet = 0;
226 G4int i = (*pos).first; 226 G4int i = (*pos).first;
227 dataMap[i] = 0; 227 dataMap[i] = 0;
228 } 228 }
229 dataMap.clear(); 229 dataMap.clear();
230 } 230 }
231 231
232 activeZ.clear(); 232 activeZ.clear();
233 ActiveElements(); 233 ActiveElements();
234 } 234 }
235 235
236 G4double G4PixeCrossSectionHandler::FindValue( 236 G4double G4PixeCrossSectionHandler::FindValue(G4int Z, G4double energy) const
237 { 237 {
238 G4double value = 0.; 238 G4double value = 0.;
239 239
240 std::map<G4int,G4IDataSet*,std::less<G4int> 240 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
241 pos = dataMap.find(Z); 241 pos = dataMap.find(Z);
242 if (pos!= dataMap.end()) 242 if (pos!= dataMap.end())
243 { 243 {
244 // The following is a workaround for STL 244 // The following is a workaround for STL ObjectSpace implementation,
245 // which does not support the standard a 245 // which does not support the standard and does not accept
246 // the syntax pos->first or pos->second 246 // the syntax pos->first or pos->second
247 // G4IDataSet* dataSet = pos->second; 247 // G4IDataSet* dataSet = pos->second;
248 G4IDataSet* dataSet = (*pos).second; 248 G4IDataSet* dataSet = (*pos).second;
249 value = dataSet->FindValue(energy); 249 value = dataSet->FindValue(energy);
250 } 250 }
251 else 251 else
252 { 252 {
253 G4cout << "WARNING: G4PixeCrossSectionHa 253 G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue(Z,e) did not find Z = "
254 << Z << G4endl; 254 << Z << G4endl;
255 } 255 }
256 return value; 256 return value;
257 } 257 }
258 258
259 G4double G4PixeCrossSectionHandler::FindValue( 259 G4double G4PixeCrossSectionHandler::FindValue(G4int Z, G4double energy,
260 G4int shellIndex) const 260 G4int shellIndex) const
261 { 261 {
262 G4double value = 0.; 262 G4double value = 0.;
263 263
264 std::map<G4int,G4IDataSet*,std::less<G4int> 264 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
265 pos = dataMap.find(Z); 265 pos = dataMap.find(Z);
266 if (pos!= dataMap.end()) 266 if (pos!= dataMap.end())
267 { 267 {
268 // The following is a workaround for STL 268 // The following is a workaround for STL ObjectSpace implementation,
269 // which does not support the standard a 269 // which does not support the standard and does not accept
270 // the syntax pos->first or pos->second 270 // the syntax pos->first or pos->second
271 // G4IDataSet* dataSet = pos->second; 271 // G4IDataSet* dataSet = pos->second;
272 G4IDataSet* dataSet = (*pos).second; 272 G4IDataSet* dataSet = (*pos).second;
273 if (shellIndex >= 0) 273 if (shellIndex >= 0)
274 { 274 {
275 G4int nComponents = (G4int)dataSet->Number << 275 G4int nComponents = dataSet->NumberOfComponents();
276 if(shellIndex < nComponents) 276 if(shellIndex < nComponents)
277 // The value is the cross section for sh 277 // The value is the cross section for shell component at given energy
278 value = dataSet->GetComponent(shellIndex 278 value = dataSet->GetComponent(shellIndex)->FindValue(energy);
279 else 279 else
280 { 280 {
281 G4cout << "WARNING: G4PixeCrossSection 281 G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue(Z,e,shell) did not find"
282 << " shellIndex= " << shellIndex 282 << " shellIndex= " << shellIndex
283 << " for Z= " 283 << " for Z= "
284 << Z << G4endl; 284 << Z << G4endl;
285 } 285 }
286 } else { 286 } else {
287 value = dataSet->FindValue(energy); 287 value = dataSet->FindValue(energy);
288 } 288 }
289 } 289 }
290 else 290 else
291 { 291 {
292 G4cout << "WARNING: G4PixeCrossSectionHa 292 G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue did not find Z = "
293 << Z << G4endl; 293 << Z << G4endl;
294 } 294 }
295 return value; 295 return value;
296 } 296 }
297 297
298 298
299 G4double G4PixeCrossSectionHandler::ValueForMa 299 G4double G4PixeCrossSectionHandler::ValueForMaterial(const G4Material* material,
300 G4double energy) const 300 G4double energy) const
301 { 301 {
302 G4double value = 0.; 302 G4double value = 0.;
303 303
304 const G4ElementVector* elementVector = mater 304 const G4ElementVector* elementVector = material->GetElementVector();
305 const G4double* nAtomsPerVolume = material-> 305 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
306 std::size_t nElements = material->GetNumberO << 306 G4int nElements = material->GetNumberOfElements();
307 307
308 for (std::size_t i=0 ; i<nElements ; ++i) << 308 for (G4int i=0 ; i<nElements ; i++)
309 { 309 {
310 G4int Z = (G4int) (*elementVector)[i]->G 310 G4int Z = (G4int) (*elementVector)[i]->GetZ();
311 G4double elementValue = FindValue(Z,ener 311 G4double elementValue = FindValue(Z,energy);
312 G4double nAtomsVol = nAtomsPerVolume[i]; 312 G4double nAtomsVol = nAtomsPerVolume[i];
313 value += nAtomsVol * elementValue; 313 value += nAtomsVol * elementValue;
314 } 314 }
315 315
316 return value; 316 return value;
317 } 317 }
318 318
319 /* 319 /*
320 G4IDataSet* G4PixeCrossSectionHandler::Build 320 G4IDataSet* G4PixeCrossSectionHandler::BuildMeanFreePathForMaterials(const G4DataVector* energyCuts )
321 { 321 {
322 // Builds a CompositeDataSet containing the 322 // Builds a CompositeDataSet containing the mean free path for each material
323 // in the material table 323 // in the material table
324 324
325 G4DataVector energyVector; 325 G4DataVector energyVector;
326 G4double dBin = std::log10(eMax/eMin) / nBin 326 G4double dBin = std::log10(eMax/eMin) / nBins;
327 327
328 for (G4int i=0; i<nBins+1; i++) 328 for (G4int i=0; i<nBins+1; i++)
329 { 329 {
330 energyVector.push_back(std::pow(10., std::lo 330 energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
331 } 331 }
332 332
333 // Factory method to build cross sections in 333 // Factory method to build cross sections in derived classes,
334 // related to the type of physics process 334 // related to the type of physics process
335 335
336 if (crossSections != 0) 336 if (crossSections != 0)
337 { // Reset the list of cross sections 337 { // Reset the list of cross sections
338 std::vector<G4IDataSet*>::iterator mat; 338 std::vector<G4IDataSet*>::iterator mat;
339 if (! crossSections->empty()) 339 if (! crossSections->empty())
340 { 340 {
341 for (mat = crossSections->begin(); mat!= cro 341 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
342 { 342 {
343 G4IDataSet* set = *mat; 343 G4IDataSet* set = *mat;
344 delete set; 344 delete set;
345 set = 0; 345 set = 0;
346 } 346 }
347 crossSections->clear(); 347 crossSections->clear();
348 delete crossSections; 348 delete crossSections;
349 crossSections = 0; 349 crossSections = 0;
350 } 350 }
351 } 351 }
352 352
353 crossSections = BuildCrossSectionsForMateria 353 crossSections = BuildCrossSectionsForMaterials(energyVector);
354 354
355 if (crossSections == 0) 355 if (crossSections == 0)
356 G4Exception("G4PixeCrossSectionHandler::Buil 356 G4Exception("G4PixeCrossSectionHandler::BuildMeanFreePathForMaterials",
357 "pii00000201", 357 "pii00000201",
358 FatalException, 358 FatalException,
359 "crossSections = 0"); 359 "crossSections = 0");
360 360
361 G4IInterpolator* algo = CreateInterpolation( 361 G4IInterpolator* algo = CreateInterpolation();
362 G4IDataSet* materialSet = new G4CompositeDat 362 G4IDataSet* materialSet = new G4CompositeDataSet(algo);
363 363
364 G4DataVector* energies; 364 G4DataVector* energies;
365 G4DataVector* data; 365 G4DataVector* data;
366 366
367 const G4ProductionCutsTable* theCoupleTable= 367 const G4ProductionCutsTable* theCoupleTable=
368 G4ProductionCutsTable::GetProductionCutsTabl 368 G4ProductionCutsTable::GetProductionCutsTable();
369 std::size_t numOfCouples = theCoupleTable->G << 369 size_t numOfCouples = theCoupleTable->GetTableSize();
370 370
371 371
372 for (std::size_t m=0; m<numOfCouples; m++) << 372 for (size_t m=0; m<numOfCouples; m++)
373 { 373 {
374 energies = new G4DataVector; 374 energies = new G4DataVector;
375 data = new G4DataVector; 375 data = new G4DataVector;
376 for (G4int bin=0; bin<nBins; bin++) 376 for (G4int bin=0; bin<nBins; bin++)
377 { 377 {
378 G4double energy = energyVector[bin]; 378 G4double energy = energyVector[bin];
379 energies->push_back(energy); 379 energies->push_back(energy);
380 G4IDataSet* matCrossSet = (*crossSections)[m 380 G4IDataSet* matCrossSet = (*crossSections)[m];
381 G4double materialCrossSection = 0.0; 381 G4double materialCrossSection = 0.0;
382 G4int nElm = matCrossSet->NumberOfComponents 382 G4int nElm = matCrossSet->NumberOfComponents();
383 for(G4int j=0; j<nElm; j++) { 383 for(G4int j=0; j<nElm; j++) {
384 materialCrossSection += matCrossSet->GetComp 384 materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
385 } 385 }
386 386
387 if (materialCrossSection > 0.) 387 if (materialCrossSection > 0.)
388 { 388 {
389 data->push_back(1./materialCrossSection); 389 data->push_back(1./materialCrossSection);
390 } 390 }
391 else 391 else
392 { 392 {
393 data->push_back(DBL_MAX); 393 data->push_back(DBL_MAX);
394 } 394 }
395 } 395 }
396 G4IInterpolator* algo = CreateInterpolation( 396 G4IInterpolator* algo = CreateInterpolation();
397 G4IDataSet* dataSet = new G4DataSet(m,energi 397 G4IDataSet* dataSet = new G4DataSet(m,energies,data,algo,1.,1.);
398 materialSet->AddComponent(dataSet); 398 materialSet->AddComponent(dataSet);
399 } 399 }
400 400
401 return materialSet; 401 return materialSet;
402 } 402 }
403 403
404 */ 404 */
405 405
406 void G4PixeCrossSectionHandler::BuildForMateri 406 void G4PixeCrossSectionHandler::BuildForMaterials()
407 { 407 {
408 // Builds a CompositeDataSet containing the 408 // Builds a CompositeDataSet containing the mean free path for each material
409 // in the material table 409 // in the material table
410 410
411 G4DataVector energyVector; 411 G4DataVector energyVector;
412 G4double dBin = std::log10(eMax/eMin) / nBin 412 G4double dBin = std::log10(eMax/eMin) / nBins;
413 413
414 for (G4int i=0; i<nBins+1; i++) 414 for (G4int i=0; i<nBins+1; i++)
415 { 415 {
416 energyVector.push_back(std::pow(10., std 416 energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
417 } 417 }
418 418
419 if (crossSections != 0) 419 if (crossSections != 0)
420 { // Reset the list of cross sections 420 { // Reset the list of cross sections
421 std::vector<G4IDataSet*>::iterator mat; 421 std::vector<G4IDataSet*>::iterator mat;
422 if (! crossSections->empty()) 422 if (! crossSections->empty())
423 { 423 {
424 for (mat = crossSections->begin(); mat!= c 424 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
425 { 425 {
426 G4IDataSet* set = *mat; 426 G4IDataSet* set = *mat;
427 delete set; 427 delete set;
428 set = 0; 428 set = 0;
429 } 429 }
430 crossSections->clear(); 430 crossSections->clear();
431 delete crossSections; 431 delete crossSections;
432 crossSections = 0; 432 crossSections = 0;
433 } 433 }
434 } 434 }
435 435
436 crossSections = BuildCrossSectionsForMateria 436 crossSections = BuildCrossSectionsForMaterials(energyVector);
437 437
438 if (crossSections == 0) 438 if (crossSections == 0)
439 G4Exception("G4PixeCrossSectionHandler::Bu 439 G4Exception("G4PixeCrossSectionHandler::BuildForMaterials",
440 "pii00000210", 440 "pii00000210",
441 FatalException, 441 FatalException,
442 ", crossSections = 0"); 442 ", crossSections = 0");
443 443
444 return; 444 return;
445 } 445 }
446 446
447 447
448 G4int G4PixeCrossSectionHandler::SelectRandomA 448 G4int G4PixeCrossSectionHandler::SelectRandomAtom(const G4Material* material,
449 G4double e) const 449 G4double e) const
450 { 450 {
451 // Select randomly an element within the mat 451 // Select randomly an element within the material, according to the weight
452 // determined by the cross sections in the d 452 // determined by the cross sections in the data set
453 453
454 G4int nElements = (G4int)material->GetNumber << 454 G4int nElements = material->GetNumberOfElements();
455 455
456 // Special case: the material consists of on 456 // Special case: the material consists of one element
457 if (nElements == 1) 457 if (nElements == 1)
458 { 458 {
459 G4int Z = (G4int) material->GetZ(); 459 G4int Z = (G4int) material->GetZ();
460 return Z; 460 return Z;
461 } 461 }
462 462
463 // Composite material 463 // Composite material
464 464
465 const G4ElementVector* elementVector = mater 465 const G4ElementVector* elementVector = material->GetElementVector();
466 std::size_t materialIndex = material->GetInd << 466 size_t materialIndex = material->GetIndex();
467 467
468 G4IDataSet* materialSet = (*crossSections)[m 468 G4IDataSet* materialSet = (*crossSections)[materialIndex];
469 G4double materialCrossSection0 = 0.0; 469 G4double materialCrossSection0 = 0.0;
470 G4DataVector cross; 470 G4DataVector cross;
471 cross.clear(); 471 cross.clear();
472 for ( G4int i=0; i < nElements; ++i ) << 472 for ( G4int i=0; i < nElements; i++ )
473 { 473 {
474 G4double cr = materialSet->GetComponent( 474 G4double cr = materialSet->GetComponent(i)->FindValue(e);
475 materialCrossSection0 += cr; 475 materialCrossSection0 += cr;
476 cross.push_back(materialCrossSection0); 476 cross.push_back(materialCrossSection0);
477 } 477 }
478 478
479 G4double random = G4UniformRand() * material 479 G4double random = G4UniformRand() * materialCrossSection0;
480 480
481 for (G4int k=0 ; k < nElements ; ++k ) << 481 for (G4int k=0 ; k < nElements ; k++ )
482 { 482 {
483 if (random <= cross[k]) return (G4int) ( 483 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
484 } 484 }
485 // It should never get here 485 // It should never get here
486 return 0; 486 return 0;
487 } 487 }
488 488
489 /* 489 /*
490 const G4Element* G4PixeCrossSectionHandler:: 490 const G4Element* G4PixeCrossSectionHandler::SelectRandomElement(const G4MaterialCutsCouple* couple,
491 G4double e) const 491 G4double e) const
492 { 492 {
493 // Select randomly an element within the mat 493 // Select randomly an element within the material, according to the weight determined
494 // by the cross sections in the data set 494 // by the cross sections in the data set
495 495
496 const G4Material* material = couple->GetMate 496 const G4Material* material = couple->GetMaterial();
497 G4Element* nullElement = 0; 497 G4Element* nullElement = 0;
498 G4int nElements = material->GetNumberOfEleme 498 G4int nElements = material->GetNumberOfElements();
499 const G4ElementVector* elementVector = mater 499 const G4ElementVector* elementVector = material->GetElementVector();
500 500
501 // Special case: the material consists of on 501 // Special case: the material consists of one element
502 if (nElements == 1) 502 if (nElements == 1)
503 { 503 {
504 G4Element* element = (*elementVector)[0]; 504 G4Element* element = (*elementVector)[0];
505 return element; 505 return element;
506 } 506 }
507 else 507 else
508 { 508 {
509 // Composite material 509 // Composite material
510 510
511 std::size_t materialIndex = couple->GetIndex << 511 size_t materialIndex = couple->GetIndex();
512 512
513 G4IDataSet* materialSet = (*crossSections)[m 513 G4IDataSet* materialSet = (*crossSections)[materialIndex];
514 G4double materialCrossSection0 = 0.0; 514 G4double materialCrossSection0 = 0.0;
515 G4DataVector cross; 515 G4DataVector cross;
516 cross.clear(); 516 cross.clear();
517 for (G4int i=0; i<nElements; i++) 517 for (G4int i=0; i<nElements; i++)
518 { 518 {
519 G4double cr = materialSet->GetComponent(i)-> 519 G4double cr = materialSet->GetComponent(i)->FindValue(e);
520 materialCrossSection0 += cr; 520 materialCrossSection0 += cr;
521 cross.push_back(materialCrossSection0); 521 cross.push_back(materialCrossSection0);
522 } 522 }
523 523
524 G4double random = G4UniformRand() * material 524 G4double random = G4UniformRand() * materialCrossSection0;
525 525
526 for (G4int k=0 ; k < nElements ; k++ ) 526 for (G4int k=0 ; k < nElements ; k++ )
527 { 527 {
528 if (random <= cross[k]) return (*elementVect 528 if (random <= cross[k]) return (*elementVector)[k];
529 } 529 }
530 // It should never end up here 530 // It should never end up here
531 G4cout << "G4PixeCrossSectionHandler::Select 531 G4cout << "G4PixeCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
532 return nullElement; 532 return nullElement;
533 } 533 }
534 } 534 }
535 */ 535 */
536 536
537 537
538 G4int G4PixeCrossSectionHandler::SelectRandomS 538 G4int G4PixeCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
539 { 539 {
540 // Select randomly a shell, according to the 540 // Select randomly a shell, according to the weight determined by the cross sections
541 // in the data set 541 // in the data set
542 542
543 // Note for later improvement: it would be u 543 // Note for later improvement: it would be useful to add a cache mechanism for already
544 // used shells to improve performance 544 // used shells to improve performance
545 545
546 G4int shell = 0; 546 G4int shell = 0;
547 547
548 G4double totCrossSection = FindValue(Z,e); 548 G4double totCrossSection = FindValue(Z,e);
549 G4double random = G4UniformRand() * totCross 549 G4double random = G4UniformRand() * totCrossSection;
550 G4double partialSum = 0.; 550 G4double partialSum = 0.;
551 551
552 G4IDataSet* dataSet = nullptr; << 552 G4IDataSet* dataSet = 0;
553 auto pos = dataMap.find(Z); << 553 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
>> 554 pos = dataMap.find(Z);
554 // The following is a workaround for STL Obj 555 // The following is a workaround for STL ObjectSpace implementation,
555 // which does not support the standard and d 556 // which does not support the standard and does not accept
556 // the syntax pos->first or pos->second 557 // the syntax pos->first or pos->second
557 // if (pos != dataMap.end()) dataSet = pos-> 558 // if (pos != dataMap.end()) dataSet = pos->second;
558 if (pos != dataMap.end()) dataSet = (*pos).s 559 if (pos != dataMap.end()) dataSet = (*pos).second;
559 560
560 if (dataSet != nullptr) { << 561 size_t nShells = dataSet->NumberOfComponents();
561 G4int nShells = (G4int)dataSet->NumberOfCo << 562 for (size_t i=0; i<nShells; i++)
562 for (G4int i=0; i<nShells; ++i) <<
563 { 563 {
564 const G4IDataSet* shellDataSet = dataSet 564 const G4IDataSet* shellDataSet = dataSet->GetComponent(i);
565 if (shellDataSet != 0) 565 if (shellDataSet != 0)
566 { << 566 {
567 G4double value = shellDataSet->FindVal << 567 G4double value = shellDataSet->FindValue(e);
568 partialSum += value; << 568 partialSum += value;
569 if (random <= partialSum) return i; << 569 if (random <= partialSum) return i;
570 } << 570 }
571 } 571 }
572 } <<
573 // It should never get here 572 // It should never get here
574 return shell; 573 return shell;
575 } 574 }
576 575
577 void G4PixeCrossSectionHandler::ActiveElements 576 void G4PixeCrossSectionHandler::ActiveElements()
578 { 577 {
579 const G4MaterialTable* materialTable = G4Mat 578 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
580 if (materialTable == 0) 579 if (materialTable == 0)
581 G4Exception("G4PixeCrossSectionHandler::Ac 580 G4Exception("G4PixeCrossSectionHandler::ActiveElements",
582 "pii00000220", 581 "pii00000220",
583 FatalException, 582 FatalException,
584 "no MaterialTable found"); 583 "no MaterialTable found");
585 584
586 std::size_t nMaterials = G4Material::GetNumb << 585 G4int nMaterials = G4Material::GetNumberOfMaterials();
587 586
588 for (std::size_t mat=0; mat<nMaterials; ++ma << 587 for (G4int mat=0; mat<nMaterials; mat++)
589 { 588 {
590 const G4Material* material= (*materialTa 589 const G4Material* material= (*materialTable)[mat];
591 const G4ElementVector* elementVector = m 590 const G4ElementVector* elementVector = material->GetElementVector();
592 const std::size_t nElements = material-> << 591 const G4int nElements = material->GetNumberOfElements();
593 592
594 for (std::size_t iEl=0; iEl<nElements; + << 593 for (G4int iEl=0; iEl<nElements; iEl++)
595 { 594 {
596 G4double Z = (*elementVector)[iEl]->GetZ() << 595 G4Element* element = (*elementVector)[iEl];
>> 596 G4double Z = element->GetZ();
597 if (!(activeZ.contains(Z)) && Z >= zMin && 597 if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
598 { 598 {
599 activeZ.push_back(Z); 599 activeZ.push_back(Z);
600 } 600 }
601 } 601 }
602 } 602 }
603 } 603 }
604 604
605 G4IInterpolator* G4PixeCrossSectionHandler::Cr 605 G4IInterpolator* G4PixeCrossSectionHandler::CreateInterpolation()
606 { 606 {
607 G4IInterpolator* algorithm = new G4LogLogInt 607 G4IInterpolator* algorithm = new G4LogLogInterpolator;
608 return algorithm; 608 return algorithm;
609 } 609 }
610 610
611 G4int G4PixeCrossSectionHandler::NumberOfCompo 611 G4int G4PixeCrossSectionHandler::NumberOfComponents(G4int Z) const
612 { 612 {
613 G4int n = 0; 613 G4int n = 0;
614 614
615 std::map<G4int,G4IDataSet*,std::less<G4int> 615 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
616 pos = dataMap.find(Z); 616 pos = dataMap.find(Z);
617 if (pos!= dataMap.end()) 617 if (pos!= dataMap.end())
618 { 618 {
619 G4IDataSet* dataSet = (*pos).second; 619 G4IDataSet* dataSet = (*pos).second;
620 n = (G4int)dataSet->NumberOfComponents() << 620 n = dataSet->NumberOfComponents();
621 } 621 }
622 else 622 else
623 { 623 {
624 G4cout << "WARNING: G4PixeCrossSectionHa 624 G4cout << "WARNING: G4PixeCrossSectionHandler::NumberOfComponents did not "
625 << "find Z = " 625 << "find Z = "
626 << Z << G4endl; 626 << Z << G4endl;
627 } 627 }
628 return n; 628 return n;
629 } 629 }
630 630
631 631
632 std::vector<G4IDataSet*>* 632 std::vector<G4IDataSet*>*
633 G4PixeCrossSectionHandler::BuildCrossSectionsF 633 G4PixeCrossSectionHandler::BuildCrossSectionsForMaterials(const G4DataVector& energyVector)
634 { 634 {
635 G4DataVector* energies; 635 G4DataVector* energies;
636 G4DataVector* data; 636 G4DataVector* data;
637 637
638 std::vector<G4IDataSet*>* matCrossSections = 638 std::vector<G4IDataSet*>* matCrossSections = new std::vector<G4IDataSet*>;
639 639
640 //const G4ProductionCutsTable* theCoupleTabl 640 //const G4ProductionCutsTable* theCoupleTable=G4ProductionCutsTable::GetProductionCutsTable();
641 //std::size_t numOfCouples = theCoupleTable- << 641 //size_t numOfCouples = theCoupleTable->GetTableSize();
642 642
643 std::size_t nOfBins = energyVector.size(); << 643 size_t nOfBins = energyVector.size();
644 const auto interpolationAlgo = std::unique_p << 644 const G4IInterpolator* interpolationAlgo = CreateInterpolation();
645 645
646 const G4MaterialTable* materialTable = G4Mat 646 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
647 if (materialTable == 0) 647 if (materialTable == 0)
648 G4Exception("G4PixeCrossSectionHandler::Bu 648 G4Exception("G4PixeCrossSectionHandler::BuildCrossSectionsForMaterials",
649 "pii00000230", 649 "pii00000230",
650 FatalException, 650 FatalException,
651 "no MaterialTable found"); 651 "no MaterialTable found");
652 652
653 std::size_t nMaterials = G4Material::GetNumb << 653 G4int nMaterials = G4Material::GetNumberOfMaterials();
654 654
655 for (std::size_t mat=0; mat<nMaterials; ++ma << 655 for (G4int mat=0; mat<nMaterials; mat++)
656 { 656 {
657 const G4Material* material = (*materialT 657 const G4Material* material = (*materialTable)[mat];
658 G4int nElements = (G4int)material->GetNu << 658 G4int nElements = material->GetNumberOfElements();
659 const G4ElementVector* elementVector = m 659 const G4ElementVector* elementVector = material->GetElementVector();
660 const G4double* nAtomsPerVolume = materi 660 const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
661 661
662 G4IInterpolator* algo = interpolationAlg 662 G4IInterpolator* algo = interpolationAlgo->Clone();
663 663
664 G4IDataSet* setForMat = new G4CompositeD 664 G4IDataSet* setForMat = new G4CompositeDataSet(algo,1.,1.);
665 665
666 for (G4int i=0; i<nElements; ++i) { << 666 for (G4int i=0; i<nElements; i++) {
667 667
668 G4int Z = (G4int) (*elementVector)[i]- 668 G4int Z = (G4int) (*elementVector)[i]->GetZ();
669 G4double density = nAtomsPerVolume[i]; 669 G4double density = nAtomsPerVolume[i];
670 670
671 energies = new G4DataVector; 671 energies = new G4DataVector;
672 data = new G4DataVector; 672 data = new G4DataVector;
673 673
674 674
675 for (std::size_t bin=0; bin<nOfBins; + << 675 for (size_t bin=0; bin<nOfBins; bin++)
676 { 676 {
677 G4double e = energyVector[bin]; 677 G4double e = energyVector[bin];
678 energies->push_back(e); 678 energies->push_back(e);
679 G4double cross = 0.; 679 G4double cross = 0.;
680 if (Z >= zMin && Z <= zMax) cross = dens 680 if (Z >= zMin && Z <= zMax) cross = density*FindValue(Z,e);
681 data->push_back(cross); 681 data->push_back(cross);
682 } 682 }
683 683
684 G4IInterpolator* algo1 = interpolation 684 G4IInterpolator* algo1 = interpolationAlgo->Clone();
685 G4IDataSet* elSet = new G4DataSet(i,en 685 G4IDataSet* elSet = new G4DataSet(i,energies,data,algo1,1.,1.);
686 setForMat->AddComponent(elSet); 686 setForMat->AddComponent(elSet);
687 } 687 }
688 688
689 matCrossSections->push_back(setForMat); 689 matCrossSections->push_back(setForMat);
690 } 690 }
691 return matCrossSections; 691 return matCrossSections;
692 } 692 }
693 693
694 694
695 G4double G4PixeCrossSectionHandler::Microscopi 695 G4double G4PixeCrossSectionHandler::MicroscopicCrossSection(const G4ParticleDefinition* particleDef,
696 G4double kineticEnergy, 696 G4double kineticEnergy,
697 G4double Z, 697 G4double Z,
698 G4double deltaCut) const 698 G4double deltaCut) const
699 { 699 {
700 // Cross section formula is OK for spin=0, 1 700 // Cross section formula is OK for spin=0, 1/2, 1 only !
701 // Calculates the microscopic cross section 701 // Calculates the microscopic cross section in Geant4 internal units
702 // Formula documented in Geant4 Phys. Ref. M 702 // Formula documented in Geant4 Phys. Ref. Manual
703 // ( it is called for elements, AtomicNumber 703 // ( it is called for elements, AtomicNumber = z )
704 704
705 G4double cross = 0.; 705 G4double cross = 0.;
706 706
707 // Particle mass and energy 707 // Particle mass and energy
708 G4double particleMass = particleDef->GetPD 708 G4double particleMass = particleDef->GetPDGMass();
709 G4double energy = kineticEnergy + particle 709 G4double energy = kineticEnergy + particleMass;
710 710
711 // Some kinematics 711 // Some kinematics
712 G4double gamma = energy / particleMass; 712 G4double gamma = energy / particleMass;
713 G4double beta2 = 1. - 1. / (gamma * gamma); 713 G4double beta2 = 1. - 1. / (gamma * gamma);
714 G4double var = electron_mass_c2 / particleMa 714 G4double var = electron_mass_c2 / particleMass;
715 G4double tMax = 2. * electron_mass_c2 * (gam 715 G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.*gamma*var + var*var);
716 716
717 // Calculate the total cross section 717 // Calculate the total cross section
718 718
719 if ( tMax > deltaCut ) 719 if ( tMax > deltaCut )
720 { 720 {
721 var = deltaCut / tMax; 721 var = deltaCut / tMax;
722 cross = (1. - var * (1. - beta2 * std::l 722 cross = (1. - var * (1. - beta2 * std::log(var))) / deltaCut;
723 723
724 G4double spin = particleDef->GetPDGSpin( 724 G4double spin = particleDef->GetPDGSpin() ;
725 725
726 // +term for spin=1/2 particle 726 // +term for spin=1/2 particle
727 if (spin == 0.5) 727 if (spin == 0.5)
728 { 728 {
729 cross += 0.5 * (tMax - deltaCut) / (energ 729 cross += 0.5 * (tMax - deltaCut) / (energy*energy);
730 } 730 }
731 // +term for spin=1 particle 731 // +term for spin=1 particle
732 else if (spin > 0.9 ) 732 else if (spin > 0.9 )
733 { 733 {
734 cross += -std::log(var) / (3.*deltaCut) + 734 cross += -std::log(var) / (3.*deltaCut) + (tMax-deltaCut) *
735 ((5.+1./var)*0.25 /(energy*energy) - bet 735 ((5.+1./var)*0.25 /(energy*energy) - beta2 / (tMax*deltaCut))/3.;
736 } 736 }
737 cross *= twopi_mc2_rcl2 * Z / beta2 ; 737 cross *= twopi_mc2_rcl2 * Z / beta2 ;
738 } 738 }
739 739
740 //std::cout << "Microscopic = " << cross/bar 740 //std::cout << "Microscopic = " << cross/barn
741 // << ", e = " << kineticEnergy/MeV <<std 741 // << ", e = " << kineticEnergy/MeV <<std:: endl;
742 742
743 return cross; 743 return cross;
744 } 744 }
745 745
746 746