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 #ifndef G4ElectroNuclearCrossSection_h
00033 #define G4ElectroNuclearCrossSection_h 1
00034
00035 #include "G4VCrossSectionDataSet.hh"
00036 #include "G4DynamicParticle.hh"
00037 #include "G4Element.hh"
00038 #include "G4ParticleTable.hh"
00039 #include "G4NucleiProperties.hh"
00040 #include <vector>
00041 #include "Randomize.hh"
00042 #include "G4Electron.hh"
00043 #include "G4Positron.hh"
00044
00045 class G4ElectroNuclearCrossSection : public G4VCrossSectionDataSet
00046 {
00047 public:
00048
00049 G4ElectroNuclearCrossSection(const G4String& name = "ElectroNuclearXS");
00050 virtual ~G4ElectroNuclearCrossSection();
00051
00052 virtual void CrossSectionDescription(std::ostream&) const;
00053
00054 virtual G4bool
00055 IsIsoApplicable(const G4DynamicParticle* aParticle, G4int ,
00056 G4int , const G4Element*, const G4Material*);
00057
00058 virtual G4double
00059 GetIsoCrossSection(const G4DynamicParticle* aParticle,
00060 G4int , G4int ,
00061 const G4Isotope*, const G4Element*, const G4Material*);
00062
00063 G4double GetEquivalentPhotonEnergy();
00064
00065 G4double GetVirtualFactor(G4double nu, G4double Q2);
00066
00067 G4double GetEquivalentPhotonQ2(G4double nu);
00068
00069 private:
00070 G4int GetFunctions(G4double a, G4double* x, G4double* y, G4double* z);
00071
00072 G4double ThresholdEnergy(G4int Z, G4int N);
00073 G4double HighEnergyJ1(G4double lE);
00074 G4double HighEnergyJ2(G4double lE);
00075 G4double HighEnergyJ3(G4double lE);
00076 G4double SolveTheEquation(G4double f);
00077 G4double Fun(G4double x);
00078 G4double DFun(G4double x);
00079
00080
00081 private:
00082 static G4int lastN;
00083 static G4int lastZ;
00084 static G4int lastF;
00085 static G4double* lastJ1;
00086 static G4double* lastJ2;
00087 static G4double* lastJ3;
00088 static G4int lastL;
00089 static G4double lastE;
00090 static G4double lastTH;
00091 static G4double lastSig;
00092 static G4double lastG;
00093 static G4double lastH;
00094
00095
00096 static std::vector <G4double*> J1;
00097
00098
00099 static std::vector <G4double*> J2;
00100
00101
00102 static std::vector <G4double*> J3;
00103 };
00104
00105
00106 inline G4double
00107 G4ElectroNuclearCrossSection::DFun(G4double x)
00108 {
00109
00110 static const G4double shd=1.0734;
00111 static const G4double poc=0.0375;
00112 static const G4double pos=16.5;
00113 static const G4double reg=.11;
00114 static const G4double mel=0.5109989;
00115 static const G4double lmel=std::log(mel);
00116 G4double y=std::exp(x-lastG-lmel);
00117 G4double flux=lastG*(2.-y*(2.-y))-1.;
00118 return (poc*(x-pos)+shd*std::exp(-reg*x))*flux;
00119 }
00120
00121
00122 inline G4double
00123 G4ElectroNuclearCrossSection::Fun(G4double x)
00124 {
00125
00126 G4double dlg1=lastG+lastG-1.;
00127 G4double lgoe=lastG/lastE;
00128 G4double HE2=HighEnergyJ2(x);
00129 return dlg1*HighEnergyJ1(x)-lgoe*(HE2+HE2-HighEnergyJ3(x)/lastE);
00130 }
00131
00132
00133 inline G4double
00134 G4ElectroNuclearCrossSection::HighEnergyJ1(G4double lEn)
00135 {
00136 static const G4double le=std::log(50000.);
00137 static const G4double le2=le*le;
00138 static const G4double a=.0375;
00139 static const G4double ha=a*.5;
00140 static const G4double ab=a*16.5;
00141 static const G4double d=0.11;
00142 static const G4double cd=1.0734/d;
00143 static const G4double ele=std::exp(-d*le);
00144 return ha*(lEn*lEn-le2)-ab*(lEn-le)-cd*(std::exp(-d*lEn)-ele);
00145 }
00146
00147
00148 inline G4double
00149 G4ElectroNuclearCrossSection::HighEnergyJ2(G4double lEn)
00150 {
00151 static const G4double e=50000.;
00152 static const G4double le=std::log(e);
00153 static const G4double le1=(le-1.)*e;
00154 static const G4double a=.0375;
00155 static const G4double ab=a*16.5;
00156 static const G4double d=1.-0.11;
00157 static const G4double cd=1.0734/d;
00158 static const G4double ele=std::exp(d*le);
00159 G4double En=std::exp(lEn);
00160 return a*((lEn-1.)*En-le1)-ab*(En-e)+cd*(std::exp(d*lEn)-ele);
00161 }
00162
00163
00164 inline G4double
00165 G4ElectroNuclearCrossSection::HighEnergyJ3(G4double lEn)
00166 {
00167 static const G4double e=50000.;
00168 static const G4double le=std::log(e);
00169 static const G4double e2=e*e;
00170 static const G4double leh=(le-.5)*e2;
00171 static const G4double ha=.0375*.5;
00172 static const G4double hab=ha*16.5;
00173 static const G4double d=2.-.11;
00174 static const G4double cd=1.0734/d;
00175 static const G4double ele=std::exp(d*le);
00176 G4double lastE2=std::exp(lEn+lEn);
00177 return ha*((lEn-.5)*lastE2-leh)-hab*(lastE2-e2)+cd*(std::exp(d*lEn)-ele);
00178 }
00179
00180 #endif