ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RDVCrossSectionHandler.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RDVCrossSectionHandler.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 // -------------------------------------------------------------------
39 
41 #include "G4RDVDataSetAlgorithm.hh"
43 #include "G4RDVEMDataSet.hh"
44 #include "G4RDEMDataSet.hh"
46 #include "G4RDShellEMDataSet.hh"
47 #include "G4ProductionCutsTable.hh"
48 #include "G4Material.hh"
49 #include "G4Element.hh"
50 #include "Randomize.hh"
51 #include <map>
52 #include <vector>
53 #include <fstream>
54 #include <sstream>
55 
56 
58 {
59  crossSections = 0;
60  interpolation = 0;
61  Initialise();
63 }
64 
65 
67  G4double minE,
68  G4double maxE,
69  G4int bins,
70  G4double unitE,
71  G4double unitData,
72  G4int minZ,
73  G4int maxZ)
74  : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
75  unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
76 {
77  crossSections = 0;
79 }
80 
82 {
83  delete interpolation;
84  interpolation = 0;
85  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
86 
87  for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
88  {
89  // The following is a workaround for STL ObjectSpace implementation,
90  // which does not support the standard and does not accept
91  // the syntax pos->second
92  // G4RDVEMDataSet* dataSet = pos->second;
93  G4RDVEMDataSet* dataSet = (*pos).second;
94  delete dataSet;
95  }
96 
97  if (crossSections != 0)
98  {
99  size_t n = crossSections->size();
100  for (size_t i=0; i<n; i++)
101  {
102  delete (*crossSections)[i];
103  }
104  delete crossSections;
105  crossSections = 0;
106  }
107 }
108 
111  G4int numberOfBins,
112  G4double unitE, G4double unitData,
113  G4int minZ, G4int maxZ)
114 {
115  if (algorithm != 0)
116  {
117  delete interpolation;
119  }
120  else
121  {
123  }
124 
125  eMin = minE;
126  eMax = maxE;
127  nBins = numberOfBins;
128  unit1 = unitE;
129  unit2 = unitData;
130  zMin = minZ;
131  zMax = maxZ;
132 }
133 
135 {
136  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
137 
138  for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
139  {
140  // The following is a workaround for STL ObjectSpace implementation,
141  // which does not support the standard and does not accept
142  // the syntax pos->first or pos->second
143  // G4int z = pos->first;
144  // G4RDVEMDataSet* dataSet = pos->second;
145  G4int z = (*pos).first;
146  G4RDVEMDataSet* dataSet = (*pos).second;
147  G4cout << "---- Data set for Z = "
148  << z
149  << G4endl;
150  dataSet->PrintData();
151  G4cout << "--------------------------------------------------" << G4endl;
152  }
153 }
154 
156 {
157  size_t nZ = activeZ.size();
158  for (size_t i=0; i<nZ; i++)
159  {
160  G4int Z = (G4int) activeZ[i];
161 
162  // Build the complete string identifying the file with the data set
163 
164  char* path = std::getenv("G4LEDATA");
165  if (!path)
166  {
167  G4String excep = "G4LEDATA environment variable not set!";
168  G4Exception("G4RDVCrossSectionHandler::LoadData()",
169  "InvalidSetup", FatalException, excep);
170  }
171 
172  std::ostringstream ost;
173  ost << path << '/' << fileName << Z << ".dat";
174  std::ifstream file(ost.str().c_str());
175  std::filebuf* lsdp = file.rdbuf();
176 
177  if (! (lsdp->is_open()) )
178  {
179  G4String excep = "Data file: " + ost.str() + " not found!";
180  G4Exception("G4RDVCrossSectionHandler::LoadData()",
181  "DataNotFound", FatalException, excep);
182  }
183  G4double a = 0;
184  G4int k = 1;
185  G4DataVector* energies = new G4DataVector;
187  do
188  {
189  file >> a;
190  G4int nColumns = 2;
191  // The file is organized into two columns:
192  // 1st column is the energy
193  // 2nd column is the corresponding value
194  // The file terminates with the pattern: -1 -1
195  // -2 -2
196  if (a == -1 || a == -2)
197  {
198  }
199  else
200  {
201  if (k%nColumns != 0)
202  {
203  G4double e = a * unit1;
204  energies->push_back(e);
205  k++;
206  }
207  else if (k%nColumns == 0)
208  {
209  G4double value = a * unit2;
210  data->push_back(value);
211  k = 1;
212  }
213  }
214  } while (a != -2); // end of file
215 
216  file.close();
218  G4RDVEMDataSet* dataSet = new G4RDEMDataSet(Z,energies,data,algo);
219  dataMap[Z] = dataSet;
220  }
221 }
222 
224 {
225  size_t nZ = activeZ.size();
226  for (size_t i=0; i<nZ; i++)
227  {
228  G4int Z = (G4int) activeZ[i];
229 
230  // Riccardo Capra <capra@ge.infn.it>: PLEASE CHECK THE FOLLOWING PIECE OF CODE
231  // "energies" AND "data" G4DataVector ARE ALLOCATED, FILLED IN AND NEVER USED OR
232  // DELETED. WHATSMORE LOADING FILE OPERATIONS WERE DONE BY G4RDShellEMDataSet
233  // EVEN BEFORE THE CHANGES I DID ON THIS FILE. SO THE FOLLOWING CODE IN MY
234  // OPINION SHOULD BE USELESS AND SHOULD PRODUCE A MEMORY LEAK.
235 
236  // Build the complete string identifying the file with the data set
237 
238  char* path = std::getenv("G4LEDATA");
239  if (!path)
240  {
241  G4String excep = "G4LEDATA environment variable not set!";
242  G4Exception("G4RDVCrossSectionHandler::LoadShellData()",
243  "InvalidSetup", FatalException, excep);
244  }
245 
246  std::ostringstream ost;
247 
248  ost << path << '/' << fileName << Z << ".dat";
249 
250  std::ifstream file(ost.str().c_str());
251  std::filebuf* lsdp = file.rdbuf();
252 
253  if (! (lsdp->is_open()) )
254  {
255  G4String excep = "Data file: " + ost.str() + " not found!";
256  G4Exception("G4RDVCrossSectionHandler::LoadShellData()",
257  "DataNotFound", FatalException, excep);
258  }
259  G4double a = 0;
260  G4int k = 1;
261  G4DataVector* energies = new G4DataVector;
263  do
264  {
265  file >> a;
266  G4int nColumns = 2;
267  // The file is organized into two columns:
268  // 1st column is the energy
269  // 2nd column is the corresponding value
270  // The file terminates with the pattern: -1 -1
271  // -2 -2
272  if (a == -1 || a == -2)
273  {
274  }
275  else
276  {
277  if (k%nColumns != 0)
278  {
279  G4double e = a * unit1;
280  energies->push_back(e);
281  k++;
282  }
283  else if (k%nColumns == 0)
284  {
285  G4double value = a * unit2;
286  data->push_back(value);
287  k = 1;
288  }
289  }
290  } while (a != -2); // end of file
291 
292  file.close();
293 
294  // Riccardo Capra <capra@ge.infn.it>: END OF CODE THAT IN MY OPINION SHOULD BE
295  // REMOVED.
296 
298  G4RDVEMDataSet* dataSet = new G4RDShellEMDataSet(Z, algo);
299  dataSet->LoadData(fileName);
300  dataMap[Z] = dataSet;
301  }
302 }
303 
305 {
306  // Reset the map of data sets: remove the data sets from the map
307  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
308 
309  if(! dataMap.empty())
310  {
311  for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
312  {
313  // The following is a workaround for STL ObjectSpace implementation,
314  // which does not support the standard and does not accept
315  // the syntax pos->first or pos->second
316  // G4RDVEMDataSet* dataSet = pos->second;
317  G4RDVEMDataSet* dataSet = (*pos).second;
318  delete dataSet;
319  dataSet = 0;
320  G4int i = (*pos).first;
321  dataMap[i] = 0;
322  }
323  dataMap.clear();
324  }
325 
326  activeZ.clear();
327  ActiveElements();
328 }
329 
331 {
332  G4double value = 0.;
333 
334  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
335  pos = dataMap.find(Z);
336  if (pos!= dataMap.end())
337  {
338  // The following is a workaround for STL ObjectSpace implementation,
339  // which does not support the standard and does not accept
340  // the syntax pos->first or pos->second
341  // G4RDVEMDataSet* dataSet = pos->second;
342  G4RDVEMDataSet* dataSet = (*pos).second;
343  value = dataSet->FindValue(energy);
344  }
345  else
346  {
347  G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find Z = "
348  << Z << G4endl;
349  }
350  return value;
351 }
352 
354  G4int shellIndex) const
355 {
356  G4double value = 0.;
357 
358  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
359  pos = dataMap.find(Z);
360  if (pos!= dataMap.end())
361  {
362  // The following is a workaround for STL ObjectSpace implementation,
363  // which does not support the standard and does not accept
364  // the syntax pos->first or pos->second
365  // G4RDVEMDataSet* dataSet = pos->second;
366  G4RDVEMDataSet* dataSet = (*pos).second;
367  if (shellIndex >= 0)
368  {
369  G4int nComponents = dataSet->NumberOfComponents();
370  if(shellIndex < nComponents)
371  // - MGP - Why doesn't it use G4RDVEMDataSet::FindValue directly?
372  value = dataSet->GetComponent(shellIndex)->FindValue(energy);
373  else
374  {
375  G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find"
376  << " shellIndex= " << shellIndex
377  << " for Z= "
378  << Z << G4endl;
379  }
380  } else {
381  value = dataSet->FindValue(energy);
382  }
383  }
384  else
385  {
386  G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find Z = "
387  << Z << G4endl;
388  }
389  return value;
390 }
391 
392 
394  G4double energy) const
395 {
396  G4double value = 0.;
397 
398  const G4ElementVector* elementVector = material->GetElementVector();
399  const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
400  G4int nElements = material->GetNumberOfElements();
401 
402  for (G4int i=0 ; i<nElements ; i++)
403  {
404  G4int Z = (G4int) (*elementVector)[i]->GetZ();
405  G4double elementValue = FindValue(Z,energy);
406  G4double nAtomsVol = nAtomsPerVolume[i];
407  value += nAtomsVol * elementValue;
408  }
409 
410  return value;
411 }
412 
413 
415 {
416  // Builds a CompositeDataSet containing the mean free path for each material
417  // in the material table
418 
419  G4DataVector energyVector;
420  G4double dBin = std::log10(eMax/eMin) / nBins;
421 
422  for (G4int i=0; i<nBins+1; i++)
423  {
424  energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
425  }
426 
427  // Factory method to build cross sections in derived classes,
428  // related to the type of physics process
429 
430  if (crossSections != 0)
431  { // Reset the list of cross sections
432  std::vector<G4RDVEMDataSet*>::iterator mat;
433  if (! crossSections->empty())
434  {
435  for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
436  {
437  G4RDVEMDataSet* set = *mat;
438  delete set;
439  set = 0;
440  }
441  crossSections->clear();
442  delete crossSections;
443  crossSections = 0;
444  }
445  }
446 
447  crossSections = BuildCrossSectionsForMaterials(energyVector,energyCuts);
448 
449  if (crossSections == 0)
450  G4Exception("G4RDVCrossSectionHandler::BuildMeanFreePathForMaterials()",
451  "InvalidCondition", FatalException, "CrossSections = 0!");
452 
454  G4RDVEMDataSet* materialSet = new G4RDCompositeEMDataSet(algo);
455 
456  G4DataVector* energies;
458 
459  const G4ProductionCutsTable* theCoupleTable=
461  size_t numOfCouples = theCoupleTable->GetTableSize();
462 
463 
464  for (size_t m=0; m<numOfCouples; m++)
465  {
466  energies = new G4DataVector;
467  data = new G4DataVector;
468  for (G4int bin=0; bin<nBins; bin++)
469  {
470  G4double energy = energyVector[bin];
471  energies->push_back(energy);
472  G4RDVEMDataSet* matCrossSet = (*crossSections)[m];
473  G4double materialCrossSection = 0.0;
474  G4int nElm = matCrossSet->NumberOfComponents();
475  for(G4int j=0; j<nElm; j++) {
476  materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
477  }
478 
479  if (materialCrossSection > 0.)
480  {
481  data->push_back(1./materialCrossSection);
482  }
483  else
484  {
485  data->push_back(DBL_MAX);
486  }
487  }
489  G4RDVEMDataSet* dataSet = new G4RDEMDataSet(m,energies,data,algo,1.,1.);
490  materialSet->AddComponent(dataSet);
491  }
492 
493  return materialSet;
494 }
495 
497  G4double e) const
498 {
499  // Select randomly an element within the material, according to the weight
500  // determined by the cross sections in the data set
501 
502  const G4Material* material = couple->GetMaterial();
503  G4int nElements = material->GetNumberOfElements();
504 
505  // Special case: the material consists of one element
506  if (nElements == 1)
507  {
508  G4int Z = (G4int) material->GetZ();
509  return Z;
510  }
511 
512  // Composite material
513 
514  const G4ElementVector* elementVector = material->GetElementVector();
515  size_t materialIndex = couple->GetIndex();
516 
517  G4RDVEMDataSet* materialSet = (*crossSections)[materialIndex];
518  G4double materialCrossSection0 = 0.0;
520  cross.clear();
521  for ( G4int i=0; i < nElements; i++ )
522  {
523  G4double cr = materialSet->GetComponent(i)->FindValue(e);
524  materialCrossSection0 += cr;
525  cross.push_back(materialCrossSection0);
526  }
527 
528  G4double random = G4UniformRand() * materialCrossSection0;
529 
530  for (G4int k=0 ; k < nElements ; k++ )
531  {
532  if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
533  }
534  // It should never get here
535  return 0;
536 }
537 
539  G4double e) const
540 {
541  // Select randomly an element within the material, according to the weight determined
542  // by the cross sections in the data set
543 
544  const G4Material* material = couple->GetMaterial();
545  G4Element* nullElement = 0;
546  G4int nElements = material->GetNumberOfElements();
547  const G4ElementVector* elementVector = material->GetElementVector();
548 
549  // Special case: the material consists of one element
550  if (nElements == 1)
551  {
552  G4Element* element = (*elementVector)[0];
553  return element;
554  }
555  else
556  {
557  // Composite material
558 
559  size_t materialIndex = couple->GetIndex();
560 
561  G4RDVEMDataSet* materialSet = (*crossSections)[materialIndex];
562  G4double materialCrossSection0 = 0.0;
564  cross.clear();
565  for (G4int i=0; i<nElements; i++)
566  {
567  G4double cr = materialSet->GetComponent(i)->FindValue(e);
568  materialCrossSection0 += cr;
569  cross.push_back(materialCrossSection0);
570  }
571 
572  G4double random = G4UniformRand() * materialCrossSection0;
573 
574  for (G4int k=0 ; k < nElements ; k++ )
575  {
576  if (random <= cross[k]) return (*elementVector)[k];
577  }
578  // It should never end up here
579  G4cout << "G4RDVCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
580  return nullElement;
581  }
582 }
583 
585 {
586  // Select randomly a shell, according to the weight determined by the cross sections
587  // in the data set
588 
589  // Note for later improvement: it would be useful to add a cache mechanism for already
590  // used shells to improve performance
591 
592  G4int shell = 0;
593 
594  G4double totCrossSection = FindValue(Z,e);
595  G4double random = G4UniformRand() * totCrossSection;
596  G4double partialSum = 0.;
597 
598  G4RDVEMDataSet* dataSet = 0;
599  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
600  pos = dataMap.find(Z);
601  // The following is a workaround for STL ObjectSpace implementation,
602  // which does not support the standard and does not accept
603  // the syntax pos->first or pos->second
604  // if (pos != dataMap.end()) dataSet = pos->second;
605  if (pos != dataMap.end()) dataSet = (*pos).second;
606 
607  size_t nShells = dataSet->NumberOfComponents();
608  for (size_t i=0; i<nShells; i++)
609  {
610  const G4RDVEMDataSet* shellDataSet = dataSet->GetComponent(i);
611  if (shellDataSet != 0)
612  {
613  G4double value = shellDataSet->FindValue(e);
614  partialSum += value;
615  if (random <= partialSum) return i;
616  }
617  }
618  // It should never get here
619  return shell;
620 }
621 
623 {
624  const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
625  if (materialTable == 0)
626  G4Exception("G4RDVCrossSectionHandler::ActiveElements",
627  "InvalidSetup", FatalException, "No MaterialTable found!");
628 
630 
631  for (G4int m=0; m<nMaterials; m++)
632  {
633  const G4Material* material= (*materialTable)[m];
634  const G4ElementVector* elementVector = material->GetElementVector();
635  const G4int nElements = material->GetNumberOfElements();
636 
637  for (G4int iEl=0; iEl<nElements; iEl++)
638  {
639  G4Element* element = (*elementVector)[iEl];
640  G4double Z = element->GetZ();
641  if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
642  {
643  activeZ.push_back(Z);
644  }
645  }
646  }
647 }
648 
650 {
652  return algorithm;
653 }
654 
656 {
657  G4int n = 0;
658 
659  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
660  pos = dataMap.find(Z);
661  if (pos!= dataMap.end())
662  {
663  G4RDVEMDataSet* dataSet = (*pos).second;
664  n = dataSet->NumberOfComponents();
665  }
666  else
667  {
668  G4cout << "WARNING: G4RDVCrossSectionHandler::NumberOfComponents did not "
669  << "find Z = "
670  << Z << G4endl;
671  }
672  return n;
673 }
674 
675