00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 #include "G4BremsstrahlungCrossSectionHandler.hh"
00062 #include "G4eBremsstrahlungSpectrum.hh"
00063 #include "G4DataVector.hh"
00064 #include "G4CompositeEMDataSet.hh"
00065 #include "G4VDataSetAlgorithm.hh"
00066 #include "G4SemiLogInterpolation.hh"
00067 #include "G4VEMDataSet.hh"
00068 #include "G4EMDataSet.hh"
00069 #include "G4Material.hh"
00070 #include "G4ProductionCutsTable.hh"
00071
00072
00073
00074 G4BremsstrahlungCrossSectionHandler::G4BremsstrahlungCrossSectionHandler(const G4VEnergySpectrum* spec,
00075 G4VDataSetAlgorithm* )
00076 : theBR(spec)
00077 {
00078 interp = new G4SemiLogInterpolation();
00079 }
00080
00081
00082
00083 G4BremsstrahlungCrossSectionHandler::~G4BremsstrahlungCrossSectionHandler()
00084 {
00085 delete interp;
00086 }
00087
00088
00089
00090 std::vector<G4VEMDataSet*>*
00091 G4BremsstrahlungCrossSectionHandler::BuildCrossSectionsForMaterials(const G4DataVector& energyVector,
00092 const G4DataVector* energyCuts)
00093 {
00094 std::vector<G4VEMDataSet*>* set = new std::vector<G4VEMDataSet*>;
00095
00096 G4DataVector* energies;
00097 G4DataVector* cs;
00098
00099 G4DataVector* log_energies;
00100 G4DataVector* log_cs;
00101
00102 G4int nOfBins = energyVector.size();
00103
00104 const G4ProductionCutsTable* theCoupleTable=
00105 G4ProductionCutsTable::GetProductionCutsTable();
00106 size_t numOfCouples = theCoupleTable->GetTableSize();
00107
00108 for (size_t mLocal=0; mLocal<numOfCouples; mLocal++) {
00109
00110 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(mLocal);
00111 const G4Material* material= couple->GetMaterial();
00112 const G4ElementVector* elementVector = material->GetElementVector();
00113 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
00114 G4int nElements = material->GetNumberOfElements();
00115
00116 G4double tcut = (*energyCuts)[mLocal];
00117
00118 G4VDataSetAlgorithm* algo = interp->Clone();
00119 G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
00120
00121 for (G4int i=0; i<nElements; i++) {
00122
00123 G4int Z = (G4int) ((*elementVector)[i]->GetZ());
00124
00125 energies = new G4DataVector;
00126 cs = new G4DataVector;
00127
00128 log_energies = new G4DataVector;
00129 log_cs = new G4DataVector;
00130
00131 G4double density = nAtomsPerVolume[i];
00132
00133 for (G4int bin=0; bin<nOfBins; bin++) {
00134
00135 G4double e = energyVector[bin];
00136 energies->push_back(e);
00137 if (e==0.) e=1e-300;
00138 log_energies->push_back(std::log10(e));
00139 G4double value = 0.0;
00140
00141 if(e > tcut) {
00142 G4double elemCs = FindValue(Z, e);
00143
00144 value = theBR->Probability(Z, tcut, e, e);
00145
00146 value *= elemCs*density;
00147 }
00148 cs->push_back(value);
00149
00150 if (value==0.) value=1e-300;
00151 log_cs->push_back(std::log10(value));
00152 }
00153 G4VDataSetAlgorithm* algol = interp->Clone();
00154
00155
00156
00157 G4VEMDataSet* elSet = new G4EMDataSet(i,energies,cs,log_energies,log_cs,algol,1.,1.);
00158
00159 setForMat->AddComponent(elSet);
00160 }
00161 set->push_back(setForMat);
00162 }
00163
00164 return set;
00165 }
00166
00167
00168
00169 G4double G4BremsstrahlungCrossSectionHandler::GetCrossSectionAboveThresholdForElement(G4double energy,
00170 G4double cutEnergy,
00171 G4int Z)
00172 {
00173 G4double value = 0.;
00174 if(energy > cutEnergy)
00175 {
00176 G4double elemCs = FindValue(Z, energy);
00177 value = theBR->Probability(Z,cutEnergy, energy, energy);
00178 value *= elemCs;
00179 }
00180 return value;
00181 }