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