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 #include "G4eIonisationCrossSectionHandler.hh"
00060 #include "G4SystemOfUnits.hh"
00061 #include "G4VEnergySpectrum.hh"
00062 #include "G4DataVector.hh"
00063 #include "G4CompositeEMDataSet.hh"
00064 #include "G4VDataSetAlgorithm.hh"
00065 #include "G4LinLogLogInterpolation.hh"
00066 #include "G4SemiLogInterpolation.hh"
00067 #include "G4VEMDataSet.hh"
00068 #include "G4EMDataSet.hh"
00069 #include "G4Material.hh"
00070 #include "G4ProductionCutsTable.hh"
00071
00072
00073 G4eIonisationCrossSectionHandler::G4eIonisationCrossSectionHandler(
00074 const G4VEnergySpectrum* spec, G4VDataSetAlgorithm* alg,
00075 G4double emin, G4double emax, G4int nbin)
00076 : G4VCrossSectionHandler(),
00077 theParam(spec),verbose(0)
00078 {
00079 G4VCrossSectionHandler::Initialise(alg, emin, emax, nbin);
00080 interp = new G4LinLogLogInterpolation();
00081 }
00082
00083
00084 G4eIonisationCrossSectionHandler::~G4eIonisationCrossSectionHandler()
00085 {
00086 delete interp;
00087 }
00088
00089
00090 std::vector<G4VEMDataSet*>* G4eIonisationCrossSectionHandler::BuildCrossSectionsForMaterials(
00091 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->GetAtomicNumDensityVector();
00114 G4int nElements = material->GetNumberOfElements();
00115
00116 if(verbose > 0)
00117 {
00118 G4cout << "eIonisation CS for " << mLocal << "th material "
00119 << material->GetName()
00120 << " eEl= " << nElements << G4endl;
00121 }
00122
00123 G4double tcut = (*energyCuts)[mLocal];
00124
00125 G4VDataSetAlgorithm* algo = interp->Clone();
00126 G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
00127
00128 for (G4int i=0; i<nElements; i++) {
00129
00130 G4int Z = (G4int) (*elementVector)[i]->GetZ();
00131 G4int nShells = NumberOfComponents(Z);
00132
00133 energies = new G4DataVector;
00134 cs = new G4DataVector;
00135
00136 log_energies = new G4DataVector;
00137 log_cs = new G4DataVector;
00138
00139 G4double density = nAtomsPerVolume[i];
00140
00141 for (G4int bin=0; bin<nOfBins; bin++) {
00142
00143 G4double e = energyVector[bin];
00144 energies->push_back(e);
00145 log_energies->push_back(std::log10(e));
00146 G4double value = 0.0;
00147 G4double log_value = -300;
00148
00149 if(e > tcut) {
00150 for (G4int n=0; n<nShells; n++) {
00151 G4double cross = FindValue(Z, e, n);
00152 G4double p = theParam->Probability(Z, tcut, e, e, n);
00153 value += cross * p * density;
00154
00155 if(verbose>0 && mLocal == 0 && e>=1. && e<=0.)
00156 {
00157 G4cout << "G4eIonCrossSH: e(MeV)= " << e/MeV
00158 << " n= " << n
00159 << " cross= " << cross
00160 << " p= " << p
00161 << " value= " << value
00162 << " tcut(MeV)= " << tcut/MeV
00163 << " rho= " << density
00164 << " Z= " << Z
00165 << G4endl;
00166 }
00167
00168 }
00169 if (value == 0.) value = 1e-300;
00170 log_value = std::log10(value);
00171 }
00172 cs->push_back(value);
00173 log_cs->push_back(log_value);
00174 }
00175 G4VDataSetAlgorithm* algoLocal = interp->Clone();
00176
00177
00178
00179 G4VEMDataSet* elSet = new G4EMDataSet(i,energies,cs,log_energies,log_cs,algoLocal,1.,1.);
00180
00181 setForMat->AddComponent(elSet);
00182 }
00183 set->push_back(setForMat);
00184 }
00185
00186 return set;
00187 }
00188
00189 G4double G4eIonisationCrossSectionHandler::GetCrossSectionAboveThresholdForElement(G4double energy,
00190 G4double cutEnergy,
00191 G4int Z)
00192 {
00193 G4int nShells = NumberOfComponents(Z);
00194 G4double value = 0.;
00195 if(energy > cutEnergy)
00196 {
00197 for (G4int n=0; n<nShells; n++) {
00198 G4double cross = FindValue(Z, energy, n);
00199 G4double p = theParam->Probability(Z, cutEnergy, energy, energy, n);
00200 value += cross * p;
00201 }
00202 }
00203 return value;
00204 }