G4PixeCrossSectionHandler Class Reference

#include <G4PixeCrossSectionHandler.hh>


Public Member Functions

 G4PixeCrossSectionHandler ()
 G4PixeCrossSectionHandler (G4IInterpolator *interpolation, const G4String &modelK="ecpssr", const G4String &modelL="ecpssr", const G4String &modelM="ecpssr", G4double minE=1 *CLHEP::keV, G4double maxE=0.1 *CLHEP::GeV, G4int nBins=200, G4double unitE=CLHEP::MeV, G4double unitData=CLHEP::barn, G4int minZ=6, G4int maxZ=92)
virtual ~G4PixeCrossSectionHandler ()
void Initialise (G4IInterpolator *interpolation, const G4String &modelK="ecpssr", const G4String &modelL="ecpssr", const G4String &modelM="ecpssr", G4double minE=1 *CLHEP::keV, G4double maxE=0.1 *CLHEP::GeV, G4int nBins=200, G4double unitE=CLHEP::MeV, G4double unitData=CLHEP::barn, G4int minZ=6, G4int maxZ=92)
G4int SelectRandomAtom (const G4Material *material, G4double e) const
G4int SelectRandomShell (G4int Z, G4double e) const
G4double FindValue (G4int Z, G4double e) const
G4double FindValue (G4int Z, G4double e, G4int shellIndex) const
G4double ValueForMaterial (const G4Material *material, G4double e) const
void LoadShellData (const G4String &dataFile)
G4double MicroscopicCrossSection (const G4ParticleDefinition *particleDef, G4double kineticEnergy, G4double Z, G4double deltaCut) const
void PrintData () const
void Clear ()


Detailed Description

Definition at line 60 of file G4PixeCrossSectionHandler.hh.


Constructor & Destructor Documentation

G4PixeCrossSectionHandler::G4PixeCrossSectionHandler (  ) 

Definition at line 61 of file G4PixeCrossSectionHandler.cc.

References Initialise().

00062 {
00063   crossSections = 0;
00064   interpolation = 0;
00065   // Initialise with default values
00066   Initialise(0,"","","",1.*keV,0.1*GeV,200,MeV,barn,6,92);
00067   ActiveElements();
00068 }

G4PixeCrossSectionHandler::G4PixeCrossSectionHandler ( G4IInterpolator interpolation,
const G4String modelK = "ecpssr",
const G4String modelL = "ecpssr",
const G4String modelM = "ecpssr",
G4double  minE = 1 *CLHEP::keV,
G4double  maxE = 0.1 *CLHEP::GeV,
G4int  nBins = 200,
G4double  unitE = CLHEP::MeV,
G4double  unitData = CLHEP::barn,
G4int  minZ = 6,
G4int  maxZ = 92 
)

Definition at line 71 of file G4PixeCrossSectionHandler.cc.

00082   : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
00083     unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
00084 {
00085   crossSections = 0;
00086 
00087   crossModel.push_back(modelK);
00088   crossModel.push_back(modelL);
00089   crossModel.push_back(modelM);
00090   
00091   //std::cout << "PixeCrossSectionHandler constructor - crossModel[0] = " 
00092   //    << crossModel[0]
00093   //    << std::endl;
00094 
00095   ActiveElements();
00096 }

G4PixeCrossSectionHandler::~G4PixeCrossSectionHandler (  )  [virtual]

Definition at line 98 of file G4PixeCrossSectionHandler.cc.

References CLHEP::detail::n.

00099 {
00100   delete interpolation;
00101   interpolation = 0;
00102   std::map<G4int,G4IDataSet*,std::less<G4int> >::iterator pos;
00103 
00104   for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
00105     {
00106       // The following is a workaround for STL ObjectSpace implementation, 
00107       // which does not support the standard and does not accept 
00108       // the syntax pos->second
00109       // G4IDataSet* dataSet = pos->second;
00110       G4IDataSet* dataSet = (*pos).second;
00111       delete dataSet;
00112     }
00113 
00114   if (crossSections != 0)
00115     {
00116       size_t n = crossSections->size();
00117       for (size_t i=0; i<n; i++)
00118         {
00119           delete (*crossSections)[i];
00120         }
00121       delete crossSections;
00122       crossSections = 0;
00123     }
00124 }


Member Function Documentation

void G4PixeCrossSectionHandler::Clear (  ) 

Definition at line 210 of file G4PixeCrossSectionHandler.cc.

00211 {
00212   // Reset the map of data sets: remove the data sets from the map 
00213   std::map<G4int,G4IDataSet*,std::less<G4int> >::iterator pos;
00214 
00215   if(! dataMap.empty())
00216     {
00217       for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
00218         {
00219           // The following is a workaround for STL ObjectSpace implementation, 
00220           // which does not support the standard and does not accept
00221           // the syntax pos->first or pos->second
00222           // G4IDataSet* dataSet = pos->second;
00223           G4IDataSet* dataSet = (*pos).second;
00224           delete dataSet;
00225           dataSet = 0;
00226           G4int i = (*pos).first;
00227           dataMap[i] = 0;
00228         }
00229       dataMap.clear();
00230     }
00231 
00232   activeZ.clear();
00233   ActiveElements();
00234 }

G4double G4PixeCrossSectionHandler::FindValue ( G4int  Z,
G4double  e,
G4int  shellIndex 
) const

Definition at line 259 of file G4PixeCrossSectionHandler.cc.

References G4cout, and G4endl.

00261 {
00262   G4double value = 0.;
00263 
00264   std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00265   pos = dataMap.find(Z);
00266   if (pos!= dataMap.end())
00267     {
00268       // The following is a workaround for STL ObjectSpace implementation, 
00269       // which does not support the standard and does not accept 
00270       // the syntax pos->first or pos->second
00271       // G4IDataSet* dataSet = pos->second;
00272       G4IDataSet* dataSet = (*pos).second;
00273       if (shellIndex >= 0) 
00274         {
00275           G4int nComponents = dataSet->NumberOfComponents();
00276           if(shellIndex < nComponents)    
00277             // The value is the cross section for shell component at given energy
00278             value = dataSet->GetComponent(shellIndex)->FindValue(energy);
00279           else 
00280             {
00281               G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue(Z,e,shell) did not find"
00282                      << " shellIndex= " << shellIndex
00283                      << " for  Z= "
00284                      << Z << G4endl;
00285             }
00286         } else {
00287         value = dataSet->FindValue(energy);
00288       }
00289     }
00290   else
00291     {
00292       G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue did not find Z = "
00293              << Z << G4endl;
00294     }
00295   return value;
00296 }

G4double G4PixeCrossSectionHandler::FindValue ( G4int  Z,
G4double  e 
) const

Definition at line 236 of file G4PixeCrossSectionHandler.cc.

References G4cout, and G4endl.

Referenced by SelectRandomShell(), and ValueForMaterial().

00237 {
00238   G4double value = 0.;
00239   
00240   std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00241   pos = dataMap.find(Z);
00242   if (pos!= dataMap.end())
00243     {
00244       // The following is a workaround for STL ObjectSpace implementation, 
00245       // which does not support the standard and does not accept 
00246       // the syntax pos->first or pos->second
00247       // G4IDataSet* dataSet = pos->second;
00248       G4IDataSet* dataSet = (*pos).second;
00249       value = dataSet->FindValue(energy);
00250     }
00251   else
00252     {
00253       G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue(Z,e) did not find Z = "
00254              << Z << G4endl;
00255     }
00256   return value;
00257 }

void G4PixeCrossSectionHandler::Initialise ( G4IInterpolator interpolation,
const G4String modelK = "ecpssr",
const G4String modelL = "ecpssr",
const G4String modelM = "ecpssr",
G4double  minE = 1 *CLHEP::keV,
G4double  maxE = 0.1 *CLHEP::GeV,
G4int  nBins = 200,
G4double  unitE = CLHEP::MeV,
G4double  unitData = CLHEP::barn,
G4int  minZ = 6,
G4int  maxZ = 92 
)

Definition at line 126 of file G4PixeCrossSectionHandler.cc.

Referenced by G4PixeCrossSectionHandler().

00134 {
00135   if (algorithm != 0) 
00136     {
00137       delete interpolation;
00138       interpolation = algorithm;
00139     }
00140   else
00141     {
00142       interpolation = CreateInterpolation();
00143     }
00144 
00145   eMin = minE;
00146   eMax = maxE;
00147   nBins = numberOfBins;
00148   unit1 = unitE;
00149   unit2 = unitData;
00150   zMin = minZ;
00151   zMax = maxZ;
00152 
00153   crossModel.push_back(modelK);
00154   crossModel.push_back(modelL);
00155   crossModel.push_back(modelM);
00156 
00157 }

void G4PixeCrossSectionHandler::LoadShellData ( const G4String dataFile  ) 

Definition at line 180 of file G4PixeCrossSectionHandler.cc.

References G4IInterpolator::Clone(), and G4IDataSet::LoadData().

Referenced by G4hImpactIonisation::PostStepDoIt().

00181 {
00182   size_t nZ = activeZ.size();
00183   for (size_t i=0; i<nZ; i++)
00184     {
00185       G4int Z = (G4int) activeZ[i];     
00186       G4IInterpolator* algo = interpolation->Clone();
00187       G4IDataSet* dataSet = new G4PixeShellDataSet(Z, algo,crossModel[0],crossModel[1],crossModel[2]);
00188 
00189       // Degug printing
00190       //std::cout << "PixeCrossSectionHandler::Load - "
00191       //        << Z
00192       //        << ", modelK = "
00193       //        << crossModel[0]
00194       //        << " fileName = "
00195       //        << fileName
00196       //        << std::endl;
00197 
00198       dataSet->LoadData(fileName);
00199       dataMap[Z] = dataSet;
00200     }
00201 
00202   // Build cross sections for materials if not already built
00203   if (! crossSections)
00204     {
00205       BuildForMaterials();
00206     }
00207 
00208 }

G4double G4PixeCrossSectionHandler::MicroscopicCrossSection ( const G4ParticleDefinition particleDef,
G4double  kineticEnergy,
G4double  Z,
G4double  deltaCut 
) const

Definition at line 695 of file G4PixeCrossSectionHandler.cc.

References G4ParticleDefinition::GetPDGMass(), and G4ParticleDefinition::GetPDGSpin().

00699 {
00700   // Cross section formula is OK for spin=0, 1/2, 1 only !
00701   // Calculates the microscopic cross section in Geant4 internal units
00702   // Formula documented in Geant4 Phys. Ref. Manual
00703   // ( it is called for elements, AtomicNumber = z )
00704 
00705     G4double cross = 0.;
00706 
00707   // Particle mass and energy
00708     G4double particleMass = particleDef->GetPDGMass();
00709     G4double energy = kineticEnergy + particleMass;
00710 
00711   // Some kinematics
00712   G4double gamma = energy / particleMass;
00713   G4double beta2 = 1. - 1. / (gamma * gamma);
00714   G4double var = electron_mass_c2 / particleMass;
00715   G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.*gamma*var + var*var);
00716 
00717   // Calculate the total cross section
00718 
00719   if ( tMax > deltaCut ) 
00720     {
00721       var = deltaCut / tMax;
00722       cross = (1. - var * (1. - beta2 * std::log(var))) / deltaCut;
00723       
00724       G4double spin = particleDef->GetPDGSpin() ;
00725       
00726       // +term for spin=1/2 particle
00727       if (spin == 0.5) 
00728         {
00729           cross +=  0.5 * (tMax - deltaCut) / (energy*energy);
00730         }
00731       // +term for spin=1 particle
00732       else if (spin > 0.9 )
00733         {
00734           cross += -std::log(var) / (3.*deltaCut) + (tMax-deltaCut) * 
00735             ((5.+1./var)*0.25 /(energy*energy) - beta2 / (tMax*deltaCut))/3.;
00736         }
00737       cross *= twopi_mc2_rcl2 * Z / beta2 ;
00738     }
00739 
00740   //std::cout << "Microscopic = " << cross/barn 
00741   //    << ", e = " << kineticEnergy/MeV <<std:: endl; 
00742 
00743   return cross;
00744 }

void G4PixeCrossSectionHandler::PrintData (  )  const

Definition at line 159 of file G4PixeCrossSectionHandler.cc.

References G4cout, and G4endl.

00160 {
00161   std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00162 
00163   for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
00164     {
00165       // The following is a workaround for STL ObjectSpace implementation, 
00166       // which does not support the standard and does not accept 
00167       // the syntax pos->first or pos->second
00168       // G4int z = pos->first;
00169       // G4IDataSet* dataSet = pos->second;
00170       G4int z = (*pos).first;
00171       G4IDataSet* dataSet = (*pos).second;     
00172       G4cout << "---- Data set for Z = "
00173              << z
00174              << G4endl;
00175       dataSet->PrintData();
00176       G4cout << "--------------------------------------------------" << G4endl;
00177     }
00178 }

G4int G4PixeCrossSectionHandler::SelectRandomAtom ( const G4Material material,
G4double  e 
) const

Definition at line 448 of file G4PixeCrossSectionHandler.cc.

References G4IDataSet::FindValue(), G4UniformRand, G4IDataSet::GetComponent(), G4Material::GetElementVector(), G4Material::GetIndex(), G4Material::GetNumberOfElements(), and G4Material::GetZ().

Referenced by G4hImpactIonisation::PostStepDoIt().

00450 {
00451   // Select randomly an element within the material, according to the weight
00452   // determined by the cross sections in the data set
00453 
00454   G4int nElements = material->GetNumberOfElements();
00455 
00456   // Special case: the material consists of one element
00457   if (nElements == 1)
00458     {
00459       G4int Z = (G4int) material->GetZ();
00460       return Z;
00461     }
00462 
00463   // Composite material
00464 
00465   const G4ElementVector* elementVector = material->GetElementVector();
00466   size_t materialIndex = material->GetIndex();
00467 
00468   G4IDataSet* materialSet = (*crossSections)[materialIndex];
00469   G4double materialCrossSection0 = 0.0;
00470   G4DataVector cross;
00471   cross.clear();
00472   for ( G4int i=0; i < nElements; i++ )
00473     {
00474       G4double cr = materialSet->GetComponent(i)->FindValue(e);
00475       materialCrossSection0 += cr;
00476       cross.push_back(materialCrossSection0);
00477     }
00478 
00479   G4double random = G4UniformRand() * materialCrossSection0;
00480 
00481   for (G4int k=0 ; k < nElements ; k++ )
00482     {
00483       if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
00484     }
00485   // It should never get here
00486   return 0;
00487 }

G4int G4PixeCrossSectionHandler::SelectRandomShell ( G4int  Z,
G4double  e 
) const

Definition at line 538 of file G4PixeCrossSectionHandler.cc.

References FindValue(), G4UniformRand, G4IDataSet::GetComponent(), and G4IDataSet::NumberOfComponents().

Referenced by G4hImpactIonisation::PostStepDoIt().

00539 {
00540   // Select randomly a shell, according to the weight determined by the cross sections
00541   // in the data set
00542 
00543   // Note for later improvement: it would be useful to add a cache mechanism for already
00544   // used shells to improve performance
00545 
00546   G4int shell = 0;
00547 
00548   G4double totCrossSection = FindValue(Z,e);
00549   G4double random = G4UniformRand() * totCrossSection;
00550   G4double partialSum = 0.;
00551 
00552   G4IDataSet* dataSet = 0;
00553   std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00554   pos = dataMap.find(Z);
00555   // The following is a workaround for STL ObjectSpace implementation,
00556   // which does not support the standard and does not accept
00557   // the syntax pos->first or pos->second
00558   // if (pos != dataMap.end()) dataSet = pos->second;
00559   if (pos != dataMap.end()) dataSet = (*pos).second;
00560 
00561   size_t nShells = dataSet->NumberOfComponents();
00562   for (size_t i=0; i<nShells; i++)
00563     {
00564       const G4IDataSet* shellDataSet = dataSet->GetComponent(i);
00565       if (shellDataSet != 0)
00566         {
00567           G4double value = shellDataSet->FindValue(e);
00568           partialSum += value;
00569           if (random <= partialSum) return i;
00570         }
00571     }
00572   // It should never get here
00573   return shell;
00574 }

G4double G4PixeCrossSectionHandler::ValueForMaterial ( const G4Material material,
G4double  e 
) const

Definition at line 299 of file G4PixeCrossSectionHandler.cc.

References FindValue(), G4Material::GetElementVector(), G4Material::GetNumberOfElements(), and G4Material::GetVecNbOfAtomsPerVolume().

00301 {
00302   G4double value = 0.;
00303 
00304   const G4ElementVector* elementVector = material->GetElementVector();
00305   const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
00306   G4int nElements = material->GetNumberOfElements();
00307 
00308   for (G4int i=0 ; i<nElements ; i++)
00309     {
00310       G4int Z = (G4int) (*elementVector)[i]->GetZ();
00311       G4double elementValue = FindValue(Z,energy);
00312       G4double nAtomsVol = nAtomsPerVolume[i];
00313       value += nAtomsVol * elementValue;
00314     }
00315 
00316   return value;
00317 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:56 2013 for Geant4 by  doxygen 1.4.7