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