ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VCrossSectionHandler.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VCrossSectionHandler.cc
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 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 //
30 // History:
31 // -----------
32 // 1 Aug 2001 MGP Created
33 // 09 Oct 2001 VI Add FindValue with 3 parameters
34 // + NumberOfComponents
35 // 19 Jul 2002 VI Create composite data set for material
36 // 21 Jan 2003 VI Cut per region
37 //
38 // 15 Jul 2009 Nicolas A. Karakatsanis
39 //
40 // - LoadNonLogData method was created to load only the non-logarithmic data from G4EMLOW
41 // dataset. It is essentially performing the data loading operations as in the past.
42 //
43 // - LoadData method was revised in order to calculate the logarithmic values of the data
44 // It retrieves the data values from the G4EMLOW data files but, then, calculates the
45 // respective log values and loads them to seperate data structures.
46 // The EM data sets, initialized this way, contain both non-log and log values.
47 // These initialized data sets can enhance the computing performance of data interpolation
48 // operations
49 //
50 // - BuildMeanFreePathForMaterials method was also revised in order to calculate the
51 // logarithmic values of the loaded data.
52 // It generates the data values and, then, calculates the respective log values which
53 // later load to seperate data structures.
54 // The EM data sets, initialized this way, contain both non-log and log values.
55 // These initialized data sets can enhance the computing performance of data interpolation
56 // operations
57 //
58 // - LoadShellData method was revised in order to eliminate the presence of a potential
59 // memory leak originally identified by Riccardo Capra.
60 // Riccardo Capra Original Comment
61 // Riccardo Capra <capra@ge.infn.it>: PLEASE CHECK THE FOLLOWING PIECE OF CODE
62 // "energies" AND "data" G4DataVector ARE ALLOCATED, FILLED IN AND NEVER USED OR
63 // DELETED. WHATSMORE LOADING FILE OPERATIONS WERE DONE BY G4ShellEMDataSet
64 // EVEN BEFORE THE CHANGES I DID ON THIS FILE. SO THE FOLLOWING CODE IN MY
65 // OPINION SHOULD BE USELESS AND SHOULD PRODUCE A MEMORY LEAK.
66 //
67 //
68 // -------------------------------------------------------------------
69 
71 #include "G4VDataSetAlgorithm.hh"
72 #include "G4LogLogInterpolation.hh"
73 #include "G4VEMDataSet.hh"
74 #include "G4EMDataSet.hh"
75 #include "G4CompositeEMDataSet.hh"
76 #include "G4ShellEMDataSet.hh"
77 #include "G4ProductionCutsTable.hh"
78 #include "G4Material.hh"
79 #include "G4Element.hh"
80 #include "Randomize.hh"
81 #include <map>
82 #include <vector>
83 #include <fstream>
84 #include <sstream>
85 
86 
88 {
89  crossSections = 0;
90  interpolation = 0;
91  Initialise();
93 }
94 
95 
97  G4double minE,
98  G4double maxE,
99  G4int bins,
100  G4double unitE,
101  G4double unitData,
102  G4int minZ,
103  G4int maxZ)
104  : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
105  unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
106 {
107  crossSections = 0;
108  ActiveElements();
109 }
110 
112 {
113  delete interpolation;
114  interpolation = 0;
115  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
116 
117  for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
118  {
119  // The following is a workaround for STL ObjectSpace implementation,
120  // which does not support the standard and does not accept
121  // the syntax pos->second
122  // G4VEMDataSet* dataSet = pos->second;
123  G4VEMDataSet* dataSet = (*pos).second;
124  delete dataSet;
125  }
126 
127  if (crossSections != 0)
128  {
129  size_t n = crossSections->size();
130  for (size_t i=0; i<n; i++)
131  {
132  delete (*crossSections)[i];
133  }
134  delete crossSections;
135  crossSections = 0;
136  }
137 }
138 
141  G4int numberOfBins,
142  G4double unitE, G4double unitData,
143  G4int minZ, G4int maxZ)
144 {
145  if (algorithm != 0)
146  {
147  delete interpolation;
149  }
150  else
151  {
152  delete interpolation;
154  }
155 
156  eMin = minE;
157  eMax = maxE;
158  nBins = numberOfBins;
159  unit1 = unitE;
160  unit2 = unitData;
161  zMin = minZ;
162  zMax = maxZ;
163 }
164 
166 {
167  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
168 
169  for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
170  {
171  // The following is a workaround for STL ObjectSpace implementation,
172  // which does not support the standard and does not accept
173  // the syntax pos->first or pos->second
174  // G4int z = pos->first;
175  // G4VEMDataSet* dataSet = pos->second;
176  G4int z = (*pos).first;
177  G4VEMDataSet* dataSet = (*pos).second;
178  G4cout << "---- Data set for Z = "
179  << z
180  << G4endl;
181  dataSet->PrintData();
182  G4cout << "--------------------------------------------------" << G4endl;
183  }
184 }
185 
187 {
188  size_t nZ = activeZ.size();
189  for (size_t i=0; i<nZ; i++)
190  {
191  G4int Z = (G4int) activeZ[i];
192 
193  // Build the complete string identifying the file with the data set
194 
195  char* path = std::getenv("G4LEDATA");
196  if (!path)
197  {
198  G4Exception("G4VCrossSectionHandler::LoadData",
199  "em0006",FatalException,"G4LEDATA environment variable not set");
200  return;
201  }
202 
203  std::ostringstream ost;
204  ost << path << '/' << fileName << Z << ".dat";
205  std::ifstream file(ost.str().c_str());
206  std::filebuf* lsdp = file.rdbuf();
207 
208  if (! (lsdp->is_open()) )
209  {
210  G4String excep = "data file: " + ost.str() + " not found";
211  G4Exception("G4VCrossSectionHandler::LoadData",
212  "em0003",FatalException,excep);
213  }
214  G4double a = 0;
215  G4int k = 0;
216  G4int nColumns = 2;
217 
218  G4DataVector* orig_reg_energies = new G4DataVector;
219  G4DataVector* orig_reg_data = new G4DataVector;
220  G4DataVector* log_reg_energies = new G4DataVector;
221  G4DataVector* log_reg_data = new G4DataVector;
222 
223  do
224  {
225  file >> a;
226 
227  if (a==0.) a=1e-300;
228 
229  // The file is organized into four columns:
230  // 1st column contains the values of energy
231  // 2nd column contains the corresponding data value
232  // The file terminates with the pattern: -1 -1
233  // -2 -2
234  //
235  if (a != -1 && a != -2)
236  {
237  if (k%nColumns == 0)
238  {
239  orig_reg_energies->push_back(a*unit1);
240  log_reg_energies->push_back(std::log10(a)+std::log10(unit1));
241  }
242  else if (k%nColumns == 1)
243  {
244  orig_reg_data->push_back(a*unit2);
245  log_reg_data->push_back(std::log10(a)+std::log10(unit2));
246  }
247  k++;
248  }
249  }
250  while (a != -2); // End of File
251 
252  file.close();
254 
255  G4VEMDataSet* dataSet = new G4EMDataSet(Z,orig_reg_energies,orig_reg_data,log_reg_energies,log_reg_data,algo);
256 
257  dataMap[Z] = dataSet;
258 
259  }
260 }
261 
262 
264 {
265  size_t nZ = activeZ.size();
266  for (size_t i=0; i<nZ; i++)
267  {
268  G4int Z = (G4int) activeZ[i];
269 
270  // Build the complete string identifying the file with the data set
271 
272  char* path = std::getenv("G4LEDATA");
273  if (!path)
274  {
275  G4Exception("G4VCrossSectionHandler::LoadNonLogData",
276  "em0006",FatalException,"G4LEDATA environment variable not set");
277  return;
278  }
279 
280  std::ostringstream ost;
281  ost << path << '/' << fileName << Z << ".dat";
282  std::ifstream file(ost.str().c_str());
283  std::filebuf* lsdp = file.rdbuf();
284 
285  if (! (lsdp->is_open()) )
286  {
287  G4String excep = "data file: " + ost.str() + " not found";
288  G4Exception("G4VCrossSectionHandler::LoadNonLogData",
289  "em0003",FatalException,excep);
290  }
291  G4double a = 0;
292  G4int k = 0;
293  G4int nColumns = 2;
294 
295  G4DataVector* orig_reg_energies = new G4DataVector;
296  G4DataVector* orig_reg_data = new G4DataVector;
297 
298  do
299  {
300  file >> a;
301 
302  // The file is organized into four columns:
303  // 1st column contains the values of energy
304  // 2nd column contains the corresponding data value
305  // The file terminates with the pattern: -1 -1
306  // -2 -2
307  //
308  if (a != -1 && a != -2)
309  {
310  if (k%nColumns == 0)
311  {
312  orig_reg_energies->push_back(a*unit1);
313  }
314  else if (k%nColumns == 1)
315  {
316  orig_reg_data->push_back(a*unit2);
317  }
318  k++;
319  }
320  }
321  while (a != -2); // End of File
322 
323  file.close();
325 
326  G4VEMDataSet* dataSet = new G4EMDataSet(Z,orig_reg_energies,orig_reg_data,algo);
327 
328  dataMap[Z] = dataSet;
329 
330  }
331 }
332 
334 {
335  size_t nZ = activeZ.size();
336  for (size_t i=0; i<nZ; i++)
337  {
338  G4int Z = (G4int) activeZ[i];
339 
341  G4VEMDataSet* dataSet = new G4ShellEMDataSet(Z, algo);
342 
343  dataSet->LoadData(fileName);
344 
345  dataMap[Z] = dataSet;
346  }
347 }
348 
349 
350 
351 
353 {
354  // Reset the map of data sets: remove the data sets from the map
355  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
356 
357  if(! dataMap.empty())
358  {
359  for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
360  {
361  // The following is a workaround for STL ObjectSpace implementation,
362  // which does not support the standard and does not accept
363  // the syntax pos->first or pos->second
364  // G4VEMDataSet* dataSet = pos->second;
365  G4VEMDataSet* dataSet = (*pos).second;
366  delete dataSet;
367  dataSet = 0;
368  G4int i = (*pos).first;
369  dataMap[i] = 0;
370  }
371  dataMap.clear();
372  }
373 
374  activeZ.clear();
375  ActiveElements();
376 }
377 
379 {
380  G4double value = 0.;
381 
382  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
383  pos = dataMap.find(Z);
384  if (pos!= dataMap.end())
385  {
386  // The following is a workaround for STL ObjectSpace implementation,
387  // which does not support the standard and does not accept
388  // the syntax pos->first or pos->second
389  // G4VEMDataSet* dataSet = pos->second;
390  G4VEMDataSet* dataSet = (*pos).second;
391  value = dataSet->FindValue(energy);
392  }
393  else
394  {
395  G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find Z = "
396  << Z << G4endl;
397  }
398  return value;
399 }
400 
402  G4int shellIndex) const
403 {
404  G4double value = 0.;
405 
406  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
407  pos = dataMap.find(Z);
408  if (pos!= dataMap.end())
409  {
410  // The following is a workaround for STL ObjectSpace implementation,
411  // which does not support the standard and does not accept
412  // the syntax pos->first or pos->second
413  // G4VEMDataSet* dataSet = pos->second;
414  G4VEMDataSet* dataSet = (*pos).second;
415  if (shellIndex >= 0)
416  {
417  G4int nComponents = dataSet->NumberOfComponents();
418  if(shellIndex < nComponents)
419  // - MGP - Why doesn't it use G4VEMDataSet::FindValue directly?
420  value = dataSet->GetComponent(shellIndex)->FindValue(energy);
421  else
422  {
423  G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find"
424  << " shellIndex= " << shellIndex
425  << " for Z= "
426  << Z << G4endl;
427  }
428  } else {
429  value = dataSet->FindValue(energy);
430  }
431  }
432  else
433  {
434  G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find Z = "
435  << Z << G4endl;
436  }
437  return value;
438 }
439 
440 
442  G4double energy) const
443 {
444  G4double value = 0.;
445 
446  const G4ElementVector* elementVector = material->GetElementVector();
447  const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
448  G4int nElements = material->GetNumberOfElements();
449 
450  for (G4int i=0 ; i<nElements ; i++)
451  {
452  G4int Z = (G4int) (*elementVector)[i]->GetZ();
453  G4double elementValue = FindValue(Z,energy);
454  G4double nAtomsVol = nAtomsPerVolume[i];
455  value += nAtomsVol * elementValue;
456  }
457 
458  return value;
459 }
460 
461 
463 {
464  // Builds a CompositeDataSet containing the mean free path for each material
465  // in the material table
466 
467  G4DataVector energyVector;
468  G4double dBin = std::log10(eMax/eMin) / nBins;
469 
470  for (G4int i=0; i<nBins+1; i++)
471  {
472  energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
473  }
474 
475  // Factory method to build cross sections in derived classes,
476  // related to the type of physics process
477 
478  if (crossSections != 0)
479  { // Reset the list of cross sections
480  std::vector<G4VEMDataSet*>::iterator mat;
481  if (! crossSections->empty())
482  {
483  for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
484  {
485  G4VEMDataSet* set = *mat;
486  delete set;
487  set = 0;
488  }
489  crossSections->clear();
490  delete crossSections;
491  crossSections = 0;
492  }
493  }
494 
495  crossSections = BuildCrossSectionsForMaterials(energyVector,energyCuts);
496 
497  if (crossSections == 0)
498  {
499  G4Exception("G4VCrossSectionHandler::BuildMeanFreePathForMaterials",
500  "em1010",FatalException,"crossSections = 0");
501  return 0;
502  }
503 
505  G4VEMDataSet* materialSet = new G4CompositeEMDataSet(algo);
506  //G4cout << "G4VCrossSectionHandler new dataset " << materialSet << G4endl;
507 
508  G4DataVector* energies;
510  G4DataVector* log_energies;
511  G4DataVector* log_data;
512 
513 
514  const G4ProductionCutsTable* theCoupleTable=
516  size_t numOfCouples = theCoupleTable->GetTableSize();
517 
518 
519  for (size_t mLocal=0; mLocal<numOfCouples; mLocal++)
520  {
521  energies = new G4DataVector;
522  data = new G4DataVector;
523  log_energies = new G4DataVector;
524  log_data = new G4DataVector;
525  for (G4int bin=0; bin<nBins; bin++)
526  {
527  G4double energy = energyVector[bin];
528  energies->push_back(energy);
529  log_energies->push_back(std::log10(energy));
530  G4VEMDataSet* matCrossSet = (*crossSections)[mLocal];
531  G4double materialCrossSection = 0.0;
532  G4int nElm = matCrossSet->NumberOfComponents();
533  for(G4int j=0; j<nElm; j++) {
534  materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
535  }
536 
537  if (materialCrossSection > 0.)
538  {
539  data->push_back(1./materialCrossSection);
540  log_data->push_back(std::log10(1./materialCrossSection));
541  }
542  else
543  {
544  data->push_back(DBL_MAX);
545  log_data->push_back(std::log10(DBL_MAX));
546  }
547  }
549 
550  //G4VEMDataSet* dataSet = new G4EMDataSet(m,energies,data,algo,1.,1.);
551 
552  G4VEMDataSet* dataSet = new G4EMDataSet(mLocal,energies,data,log_energies,log_data,algoLocal,1.,1.);
553 
554  materialSet->AddComponent(dataSet);
555  }
556 
557  return materialSet;
558 }
559 
560 
562  G4double e) const
563 {
564  // Select randomly an element within the material, according to the weight
565  // determined by the cross sections in the data set
566 
567  const G4Material* material = couple->GetMaterial();
568  G4int nElements = material->GetNumberOfElements();
569 
570  // Special case: the material consists of one element
571  if (nElements == 1)
572  {
573  G4int Z = (G4int) material->GetZ();
574  return Z;
575  }
576 
577  // Composite material
578 
579  const G4ElementVector* elementVector = material->GetElementVector();
580  size_t materialIndex = couple->GetIndex();
581 
582  G4VEMDataSet* materialSet = (*crossSections)[materialIndex];
583  G4double materialCrossSection0 = 0.0;
585  cross.clear();
586  for ( G4int i=0; i < nElements; i++ )
587  {
588  G4double cr = materialSet->GetComponent(i)->FindValue(e);
589  materialCrossSection0 += cr;
590  cross.push_back(materialCrossSection0);
591  }
592 
593  G4double random = G4UniformRand() * materialCrossSection0;
594 
595  for (G4int k=0 ; k < nElements ; k++ )
596  {
597  if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
598  }
599  // It should never get here
600  return 0;
601 }
602 
604  G4double e) const
605 {
606  // Select randomly an element within the material, according to the weight determined
607  // by the cross sections in the data set
608 
609  const G4Material* material = couple->GetMaterial();
610  G4Element* nullElement = 0;
611  G4int nElements = material->GetNumberOfElements();
612  const G4ElementVector* elementVector = material->GetElementVector();
613 
614  // Special case: the material consists of one element
615  if (nElements == 1)
616  {
617  G4Element* element = (*elementVector)[0];
618  return element;
619  }
620  else
621  {
622  // Composite material
623 
624  size_t materialIndex = couple->GetIndex();
625 
626  G4VEMDataSet* materialSet = (*crossSections)[materialIndex];
627  G4double materialCrossSection0 = 0.0;
629  cross.clear();
630  for (G4int i=0; i<nElements; i++)
631  {
632  G4double cr = materialSet->GetComponent(i)->FindValue(e);
633  materialCrossSection0 += cr;
634  cross.push_back(materialCrossSection0);
635  }
636 
637  G4double random = G4UniformRand() * materialCrossSection0;
638 
639  for (G4int k=0 ; k < nElements ; k++ )
640  {
641  if (random <= cross[k]) return (*elementVector)[k];
642  }
643  // It should never end up here
644  G4cout << "G4VCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
645  return nullElement;
646  }
647 }
648 
650 {
651  // Select randomly a shell, according to the weight determined by the cross sections
652  // in the data set
653 
654  // Note for later improvement: it would be useful to add a cache mechanism for already
655  // used shells to improve performance
656 
657  G4int shell = 0;
658 
659  G4double totCrossSection = FindValue(Z,e);
660  G4double random = G4UniformRand() * totCrossSection;
661  G4double partialSum = 0.;
662 
663  G4VEMDataSet* dataSet = 0;
664  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
665  pos = dataMap.find(Z);
666  // The following is a workaround for STL ObjectSpace implementation,
667  // which does not support the standard and does not accept
668  // the syntax pos->first or pos->second
669  // if (pos != dataMap.end()) dataSet = pos->second;
670  if (pos != dataMap.end())
671  dataSet = (*pos).second;
672  else
673  {
674  G4Exception("G4VCrossSectionHandler::SelectRandomShell",
675  "em1011",FatalException,"unable to load the dataSet");
676  return 0;
677  }
678 
679  size_t nShells = dataSet->NumberOfComponents();
680  for (size_t i=0; i<nShells; i++)
681  {
682  const G4VEMDataSet* shellDataSet = dataSet->GetComponent(i);
683  if (shellDataSet != 0)
684  {
685  G4double value = shellDataSet->FindValue(e);
686  partialSum += value;
687  if (random <= partialSum) return i;
688  }
689  }
690  // It should never get here
691  return shell;
692 }
693 
695 {
696  const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
697  if (materialTable == 0)
698  G4Exception("G4VCrossSectionHandler::ActiveElements",
699  "em1001",FatalException,"no MaterialTable found");
700 
702 
703  for (G4int mLocal2=0; mLocal2<nMaterials; mLocal2++)
704  {
705  const G4Material* material= (*materialTable)[mLocal2];
706  const G4ElementVector* elementVector = material->GetElementVector();
707  const G4int nElements = material->GetNumberOfElements();
708 
709  for (G4int iEl=0; iEl<nElements; iEl++)
710  {
711  G4Element* element = (*elementVector)[iEl];
712  G4double Z = element->GetZ();
713  if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
714  {
715  activeZ.push_back(Z);
716  }
717  }
718  }
719 }
720 
722 {
724  return algorithm;
725 }
726 
728 {
729  G4int n = 0;
730 
731  std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
732  pos = dataMap.find(Z);
733  if (pos!= dataMap.end())
734  {
735  G4VEMDataSet* dataSet = (*pos).second;
736  n = dataSet->NumberOfComponents();
737  }
738  else
739  {
740  G4cout << "WARNING: G4VCrossSectionHandler::NumberOfComponents did not "
741  << "find Z = "
742  << Z << G4endl;
743  }
744  return n;
745 }
746 
747