#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 () |
Definition at line 60 of file G4PixeCrossSectionHandler.hh.
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 }
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 }
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 }
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 }
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 }