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 #include "G4RToEConvForGamma.hh"
00036 #include "G4ParticleDefinition.hh"
00037 #include "G4ParticleTable.hh"
00038 #include "G4Material.hh"
00039 #include "G4PhysicsLogVector.hh"
00040
00041 #include "G4ios.hh"
00042 #include "G4SystemOfUnits.hh"
00043
00044 G4RToEConvForGamma::G4RToEConvForGamma() : G4VRangeToEnergyConverter()
00045 {
00046 theParticle = G4ParticleTable::GetParticleTable()->FindParticle("gamma");
00047 if (theParticle ==0) {
00048 #ifdef G4VERBOSE
00049 if (GetVerboseLevel()>0) {
00050 G4cout << " G4RToEConvForGamma::G4RToEConvForGamma() ";
00051 G4cout << " Gamma is not defined !!" << G4endl;
00052 }
00053 #endif
00054 }
00055 }
00056
00057 G4RToEConvForGamma::~G4RToEConvForGamma()
00058 {
00059 }
00060
00061
00062
00063
00064
00065 void G4RToEConvForGamma::BuildAbsorptionLengthVector(
00066 const G4Material* aMaterial,
00067 G4RangeVector* absorptionLengthVector )
00068 {
00069
00070
00071
00072
00073
00074 const G4CrossSectionTable* aCrossSectionTable = (G4CrossSectionTable*)(theLossTable);
00075 const G4ElementVector* elementVector = aMaterial->GetElementVector();
00076 const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
00077
00078
00079 G4int NumEl = aMaterial->GetNumberOfElements();
00080 G4double absorptionLengthMax = 0.0;
00081 for (size_t ibin=0; ibin<size_t(TotBin); ibin++) {
00082 G4double SIGMA = 0. ;
00083 for (size_t iel=0; iel<size_t(NumEl); iel++) {
00084 G4int IndEl = (*elementVector)[iel]->GetIndex();
00085 SIGMA += atomicNumDensityVector[iel]*
00086 (*((*aCrossSectionTable)[IndEl]))[ibin];
00087 }
00088
00089 absorptionLengthVector->PutValue(ibin, 5./SIGMA);
00090 if (absorptionLengthMax < 5./SIGMA ) absorptionLengthMax = 5./SIGMA;
00091 }
00092 }
00093
00094
00095
00096
00097
00098
00099 G4double G4RToEConvForGamma::ComputeCrossSection(G4double AtomicNumber,
00100 G4double KineticEnergy) const
00101 {
00102
00103
00104
00105 static G4double Z;
00106 const G4double t1keV = 1.*keV;
00107 const G4double t200keV = 200.*keV;
00108 const G4double t100MeV = 100.*MeV;
00109
00110 static G4double s200keV, s1keV;
00111 static G4double tmin, tlow;
00112 static G4double smin, slow;
00113 static G4double cmin, clow, chigh;
00114
00115 if(std::abs(AtomicNumber-Z)>0.1) {
00116 Z = AtomicNumber;
00117 G4double Zsquare = Z*Z;
00118 G4double Zlog = std::log(Z);
00119 G4double Zlogsquare = Zlog*Zlog;
00120
00121 s200keV = (0.2651-0.1501*Zlog+0.02283*Zlogsquare)*Zsquare;
00122 tmin = (0.552+218.5/Z+557.17/Zsquare)*MeV;
00123 smin = (0.01239+0.005585*Zlog-0.000923*Zlogsquare)*std::exp(1.5*Zlog);
00124 cmin=std::log(s200keV/smin)/(std::log(tmin/t200keV)*std::log(tmin/t200keV));
00125 tlow = 0.2*std::exp(-7.355/std::sqrt(Z))*MeV;
00126 slow = s200keV*std::exp(0.042*Z*std::log(t200keV/tlow)*std::log(t200keV/tlow));
00127 s1keV = 300.*Zsquare;
00128 clow =std::log(s1keV/slow)/std::log(tlow/t1keV);
00129
00130 chigh=(7.55e-5-0.0542e-5*Z)*Zsquare*Z/std::log(t100MeV/tmin);
00131 }
00132
00133
00134 G4double xs;
00135 if ( KineticEnergy<tlow ) {
00136 if(KineticEnergy<t1keV) xs = slow*std::exp(clow*std::log(tlow/t1keV));
00137 else xs = slow*std::exp(clow*std::log(tlow/KineticEnergy));
00138
00139 } else if ( KineticEnergy<t200keV ) {
00140 xs = s200keV
00141 * std::exp(0.042*Z*std::log(t200keV/KineticEnergy)*std::log(t200keV/KineticEnergy));
00142
00143 } else if( KineticEnergy<tmin ){
00144 xs = smin
00145 * std::exp(cmin*std::log(tmin/KineticEnergy)*std::log(tmin/KineticEnergy));
00146
00147 } else {
00148 xs = smin + chigh*std::log(KineticEnergy/tmin);
00149
00150 }
00151 return xs * barn;
00152 }
00153