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 #ifndef G4QMuonNuclearCrossSection_h
00046 #define G4QMuonNuclearCrossSection_h 1
00047
00048 #include "G4ParticleTable.hh"
00049 #include "G4NucleiProperties.hh"
00050 #include <vector>
00051 #include "Randomize.hh"
00052 #include "G4MuonPlus.hh"
00053 #include "G4MuonMinus.hh"
00054 #include "G4VQCrossSection.hh"
00055
00056 class G4QMuonNuclearCrossSection : public G4VQCrossSection
00057 {
00058 protected:
00059
00060 G4QMuonNuclearCrossSection() {};
00061
00062 public:
00063
00064 ~G4QMuonNuclearCrossSection();
00065
00066 static G4VQCrossSection* GetPointer();
00067
00068 G4double ThresholdEnergy(G4int Z, G4int N, G4int PDG=13);
00069
00070
00071 virtual G4double GetCrossSection(G4bool fCS, G4double pMom, G4int tgZ, G4int tgN,
00072 G4int pPDG=0);
00073
00074 G4double CalculateCrossSection(G4bool CS, G4int F, G4int I, G4int PDG, G4int Z, G4int N,
00075 G4double Momentum);
00076
00077 G4int GetExchangePDGCode();
00078
00079 G4double GetExchangeEnergy();
00080
00081 G4double GetVirtualFactor(G4double nu, G4double Q2);
00082
00083 G4double GetExchangeQ2(G4double nu);
00084
00085 private:
00086 G4int GetFunctions(G4double a, G4double* x, G4double* y, G4double* z);
00087 G4double HighEnergyJ1(G4double lE);
00088 G4double HighEnergyJ2(G4double lE);
00089 G4double HighEnergyJ3(G4double lE);
00090 G4double SolveTheEquation(G4double f);
00091 G4double Fun(G4double x);
00092 G4double DFun(G4double x);
00093
00094
00095 private:
00096 static G4bool onlyCS;
00097 static G4double lastSig;
00098 static G4int lastL;
00099 static G4double lastE;
00100 static G4int lastF;
00101 static G4double lastG;
00102 static G4double lastH;
00103 static G4double* lastJ1;
00104 static G4double* lastJ2;
00105 static G4double* lastJ3;
00106 static G4int lastPDG;
00107 static G4int lastN;
00108 static G4int lastZ;
00109 static G4double lastP;
00110 static G4double lastTH;
00111 static G4double lastCS;
00112 static G4int lastI;
00113 static std::vector <G4double*>* J1;
00114 static std::vector <G4double*>* J2;
00115 static std::vector <G4double*>* J3;
00116 };
00117
00118 inline G4double G4QMuonNuclearCrossSection::DFun(G4double x)
00119 {
00120 static const G4double shd=1.0734;
00121 static const G4double poc=0.0375;
00122 static const G4double pos=16.5;
00123 static const G4double reg=.11;
00124 static const G4double mmu=105.65839;
00125 static const G4double lmmu=std::log(mmu);
00126 G4double y=std::exp(x-lastG-lmmu);
00127 G4double flux=lastG*(2.-y*(2.-y))-1.;
00128 return (poc*(x-pos)+shd*std::exp(-reg*x))*flux;
00129 }
00130
00131 inline G4double G4QMuonNuclearCrossSection::Fun(G4double x)
00132 {
00133 G4double dlg1=lastG+lastG-1.;
00134 G4double lgoe=lastG/lastE;
00135 G4double HE2=HighEnergyJ2(x);
00136 return dlg1*HighEnergyJ1(x)-lgoe*(HE2+HE2-HighEnergyJ3(x)/lastE);
00137 }
00138
00139 inline G4double G4QMuonNuclearCrossSection::HighEnergyJ1(G4double lEn)
00140 {
00141 static const G4double le=std::log(50000.);
00142 static const G4double le2=le*le;
00143 static const G4double a=.0375;
00144 static const G4double ha=a*.5;
00145 static const G4double ab=a*16.5;
00146 static const G4double d=0.11;
00147 static const G4double cd=1.0734/d;
00148 static const G4double ele=std::exp(-d*le);
00149 return ha*(lEn*lEn-le2)-ab*(lEn-le)-cd*(std::exp(-d*lEn)-ele);
00150 }
00151
00152 inline G4double G4QMuonNuclearCrossSection::HighEnergyJ2(G4double lEn)
00153 {
00154 static const G4double e=50000.;
00155 static const G4double le=std::log(e);
00156 static const G4double le1=(le-1.)*e;
00157 static const G4double a=.0375;
00158 static const G4double ab=a*16.5;
00159 static const G4double d=1.-0.11;
00160 static const G4double cd=1.0734/d;
00161 static const G4double ele=std::exp(d*le);
00162 G4double En=std::exp(lEn);
00163 return a*((lEn-1.)*En-le1)-ab*(En-e)+cd*(std::exp(d*lEn)-ele);
00164 }
00165
00166 inline G4double G4QMuonNuclearCrossSection::HighEnergyJ3(G4double lEn)
00167 {
00168 static const G4double e=50000.;
00169 static const G4double le=std::log(e);
00170 static const G4double e2=e*e;
00171 static const G4double leh=(le-.5)*e2;
00172 static const G4double ha=.0375*.5;
00173 static const G4double hab=ha*16.5;
00174 static const G4double d=2.-.11;
00175 static const G4double cd=1.0734/d;
00176 static const G4double ele=std::exp(d*le);
00177 G4double lastE2=std::exp(lEn+lEn);
00178 return ha*((lEn-.5)*lastE2-leh)-hab*(lastE2-e2)+cd*(std::exp(d*lEn)-ele);
00179 }
00180
00181 #endif