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 #include "G4KokoulinMuonNuclearXS.hh"
00040
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4PhysicsTable.hh"
00044 #include "G4PhysicsLogVector.hh"
00045
00046
00047 G4KokoulinMuonNuclearXS::G4KokoulinMuonNuclearXS()
00048 :G4VCrossSectionDataSet("KokoulinMuonNuclearXS"), theCrossSectionTable(0),
00049 LowestKineticEnergy(1*GeV), HighestKineticEnergy(1*PeV),
00050 TotBin(60), CutFixed(0.2*GeV)
00051 {
00052 BuildCrossSectionTable();
00053 }
00054
00055 G4KokoulinMuonNuclearXS::~G4KokoulinMuonNuclearXS()
00056 {
00057 if (theCrossSectionTable) {
00058 theCrossSectionTable->clearAndDestroy();
00059 delete theCrossSectionTable;
00060 }
00061 }
00062
00063 G4bool
00064 G4KokoulinMuonNuclearXS::IsElementApplicable(const G4DynamicParticle*,
00065 G4int, const G4Material*)
00066 {
00067 return true;
00068 }
00069
00070
00071 void G4KokoulinMuonNuclearXS::BuildCrossSectionTable()
00072 {
00073 G4double lowEdgeEnergy;
00074 G4PhysicsLogVector* ptrVector;
00075
00076 G4int nEl = G4Element::GetNumberOfElements();
00077 const G4ElementTable* theElementTable = G4Element::GetElementTable();
00078
00079 theCrossSectionTable = new G4PhysicsTable(nEl);
00080
00081 G4double AtomicNumber;
00082 G4double AtomicWeight;
00083 G4double Value;
00084
00085 for (G4int j = 0; j < nEl; j++) {
00086 ptrVector = new G4PhysicsLogVector(LowestKineticEnergy,
00087 HighestKineticEnergy, TotBin);
00088 AtomicNumber = (*theElementTable)[j]->GetZ();
00089 AtomicWeight = (*theElementTable)[j]->GetA();
00090
00091 for (G4int i = 0; i <= TotBin; ++i) {
00092 lowEdgeEnergy = ptrVector->Energy(i);
00093 Value = ComputeMicroscopicCrossSection(lowEdgeEnergy,
00094 AtomicNumber, AtomicWeight);
00095 ptrVector->PutValue(i,Value);
00096 }
00097
00098 theCrossSectionTable->insertAt(j, ptrVector);
00099 }
00100
00101
00102 for (G4int j = 0; j < nEl; j++) {
00103 G4int ZZ = G4int((*theElementTable)[j]->GetZ());
00104 zelMap[ZZ] = j;
00105 }
00106 }
00107
00108 G4double G4KokoulinMuonNuclearXS::
00109 ComputeMicroscopicCrossSection(G4double KineticEnergy,
00110 G4double AtomicNumber, G4double AtomicWeight)
00111 {
00112
00113
00114
00115 const G4double xgi[] = {0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801};
00116 const G4double wgi[] = {0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506};
00117 const G4double ak1 = 6.9;
00118 const G4double ak2 = 1.0;
00119
00120 G4double Mass = G4MuonMinus::MuonMinus()->GetPDGMass();
00121
00122 G4double CrossSection = 0.0;
00123 if (AtomicNumber < 1.) return CrossSection;
00124 if (KineticEnergy <= CutFixed) return CrossSection;
00125
00126 G4double epmin = CutFixed;
00127 G4double epmax = KineticEnergy + Mass - 0.5*proton_mass_c2;
00128 if (epmax <= epmin) return CrossSection;
00129
00130 G4double aaa = std::log(epmin) ;
00131 G4double bbb = std::log(epmax) ;
00132 G4int kkk = G4int((bbb-aaa)/ak1 +ak2) ;
00133 G4double hhh = (bbb-aaa)/kkk ;
00134 G4double epln;
00135 G4double ep;
00136 G4double x;
00137
00138 for (G4int l = 0; l < kkk; l++) {
00139 x = aaa + hhh*l;
00140 for (G4int ll = 0; ll < 8; ll++) {
00141 epln=x+xgi[ll]*hhh;
00142 ep = std::exp(epln);
00143 CrossSection += ep*wgi[ll]*ComputeDDMicroscopicCrossSection(KineticEnergy,
00144 AtomicNumber, AtomicWeight, ep);
00145 }
00146 }
00147
00148 CrossSection *= hhh ;
00149 if (CrossSection < 0.) CrossSection = 0.;
00150 return CrossSection;
00151 }
00152
00153 G4double G4KokoulinMuonNuclearXS::
00154 ComputeDDMicroscopicCrossSection(G4double KineticEnergy,
00155 G4double, G4double AtomicWeight,
00156 G4double epsilon)
00157 {
00158
00159
00160
00161
00162 const G4double alam2 = 0.400*GeV*GeV;
00163 const G4double alam = 0.632456*GeV;
00164 const G4double coeffn = fine_structure_const/pi;
00165
00166 G4double ParticleMass = G4MuonMinus::MuonMinus()->GetPDGMass();
00167 G4double TotalEnergy = KineticEnergy + ParticleMass;
00168
00169 G4double DCrossSection = 0.;
00170
00171 if ((epsilon >= TotalEnergy - 0.5*proton_mass_c2) ||
00172 (epsilon <= CutFixed) ) return DCrossSection;
00173
00174 G4double ep = epsilon/GeV;
00175 G4double a = AtomicWeight/(g/mole);
00176 G4double aeff = 0.22*a+0.78*std::exp(0.89*std::log(a));
00177 G4double sigph = (49.2+11.1*std::log(ep)+151.8/std::sqrt(ep))*microbarn;
00178
00179 G4double v = epsilon/TotalEnergy;
00180 G4double v1 = 1.-v;
00181 G4double v2 = v*v;
00182 G4double mass2 = ParticleMass*ParticleMass;
00183
00184 G4double up = TotalEnergy*TotalEnergy*v1/mass2*(1.+mass2*v2/(alam2*v1));
00185 G4double down = 1.+epsilon/alam*(1.+alam/(2.*proton_mass_c2)+epsilon/alam);
00186
00187 DCrossSection = coeffn*aeff*sigph/epsilon*
00188 (-v1+(v1+0.5*v2*(1.+2.*mass2/alam2))*std::log(up/down));
00189
00190 if (DCrossSection < 0.) DCrossSection = 0.;
00191 return DCrossSection;
00192 }
00193
00194 G4double G4KokoulinMuonNuclearXS::
00195 GetElementCrossSection(const G4DynamicParticle* aPart,
00196 G4int ZZ, const G4Material*)
00197 {
00198 G4int j = zelMap[ZZ];
00199 return (*theCrossSectionTable)[j]->Value(aPart->GetKineticEnergy());
00200 }
00201