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