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
00047
00048
00049
00050 #ifndef G4PairProductionRelModel_h
00051 #define G4PairProductionRelModel_h 1
00052
00053 #include <CLHEP/Units/PhysicalConstants.h>
00054
00055 #include "G4VEmModel.hh"
00056 #include "G4PhysicsTable.hh"
00057 #include "G4NistManager.hh"
00058
00059 class G4ParticleChangeForGamma;
00060
00061 class G4PairProductionRelModel : public G4VEmModel
00062 {
00063
00064 public:
00065
00066 G4PairProductionRelModel(const G4ParticleDefinition* p = 0,
00067 const G4String& nam = "BetheHeitlerLPM");
00068
00069 virtual ~G4PairProductionRelModel();
00070
00071 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
00072
00073 virtual G4double ComputeCrossSectionPerAtom(
00074 const G4ParticleDefinition*,
00075 G4double kinEnergy,
00076 G4double Z,
00077 G4double A=0.,
00078 G4double cut=0.,
00079 G4double emax=DBL_MAX);
00080
00081 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
00082 const G4MaterialCutsCouple*,
00083 const G4DynamicParticle*,
00084 G4double tmin,
00085 G4double maxEnergy);
00086
00087 virtual void SetupForMaterial(const G4ParticleDefinition*,
00088 const G4Material*,G4double);
00089
00090
00091 inline void SetCurrentElement(G4double );
00092
00093
00094 inline void SetLPMconstant(G4double val);
00095 inline G4double LPMconstant() const;
00096
00097 inline void SetLPMflag(G4bool);
00098 inline G4bool LPMflag() const;
00099
00100 protected:
00101
00102 inline G4double Phi1(G4double delta) const;
00103 inline G4double Phi2(G4double delta) const;
00104 inline G4double DeltaMax() const;
00105 inline G4double DeltaMin(G4double) const;
00106
00107 void CalcLPMFunctions(G4double k, G4double eplus);
00108
00109
00110 G4double ScreenFunction1(G4double ScreenVariable);
00111 G4double ScreenFunction2(G4double ScreenVariable);
00112
00113 G4double ComputeXSectionPerAtom(G4double totalEnergy, G4double Z);
00114
00115 G4double ComputeDXSectionPerAtom(G4double eplusEnergy, G4double totalEnergy, G4double Z);
00116 G4double ComputeRelDXSectionPerAtom(G4double eplusEnergy, G4double totalEnergy, G4double Z);
00117
00118
00119 G4PairProductionRelModel & operator=(const G4PairProductionRelModel &right);
00120 G4PairProductionRelModel(const G4PairProductionRelModel&);
00121
00122 G4NistManager* nist;
00123
00124 G4ParticleDefinition* theGamma;
00125 G4ParticleDefinition* theElectron;
00126 G4ParticleDefinition* thePositron;
00127 G4ParticleChangeForGamma* fParticleChange;
00128
00129 G4double fLPMconstant;
00130 G4bool fLPMflag;
00131
00132
00133 G4double z13, z23, lnZ;
00134 G4double Fel, Finel, fCoulomb;
00135 G4double currentZ;
00136
00137
00138 G4double lpmEnergy;
00139 G4double xiLPM, phiLPM, gLPM;
00140
00141
00142 G4bool use_completescreening;
00143
00144 static const G4double xgi[8], wgi[8];
00145 static const G4double Fel_light[5];
00146 static const G4double Finel_light[5];
00147 static const G4double facFel;
00148 static const G4double facFinel;
00149
00150 static const G4double preS1, logTwo;
00151
00152 };
00153
00154
00155
00156
00157 inline
00158 void G4PairProductionRelModel::SetLPMconstant(G4double val)
00159 {
00160 fLPMconstant = val;
00161 }
00162
00163
00164
00165 inline
00166 G4double G4PairProductionRelModel::LPMconstant() const
00167 {
00168 return fLPMconstant;
00169 }
00170
00171
00172
00173 inline
00174 void G4PairProductionRelModel::SetLPMflag(G4bool val)
00175 {
00176 fLPMflag = val;
00177 }
00178
00179
00180
00181 inline
00182 G4bool G4PairProductionRelModel::LPMflag() const
00183 {
00184 return fLPMflag;
00185 }
00186
00187
00188
00189 inline void G4PairProductionRelModel::SetCurrentElement(G4double Z)
00190 {
00191 if(Z != currentZ) {
00192 currentZ = Z;
00193
00194 G4int iz = G4int(Z);
00195 z13 = nist->GetZ13(iz);
00196 z23 = z13*z13;
00197 lnZ = nist->GetLOGZ(iz);
00198
00199 if (iz <= 4) {
00200 Fel = Fel_light[iz];
00201 Finel = Finel_light[iz] ;
00202 }
00203 else {
00204 Fel = facFel - lnZ/3. ;
00205 Finel = facFinel - 2.*lnZ/3. ;
00206 }
00207 fCoulomb=GetCurrentElement()->GetfCoulomb();
00208 }
00209 }
00210
00211
00212
00213 inline G4double G4PairProductionRelModel::Phi1(G4double delta) const
00214 {
00215 G4double screenVal;
00216
00217 if (delta > 1.)
00218 screenVal = 21.12 - 4.184*std::log(delta+0.952);
00219 else
00220 screenVal = 20.868 - delta*(3.242 - 0.625*delta);
00221
00222 return screenVal;
00223 }
00224
00225
00226
00227 inline G4double G4PairProductionRelModel::Phi2(G4double delta) const
00228 {
00229 G4double screenVal;
00230
00231 if (delta > 1.)
00232 screenVal = 21.12 - 4.184*std::log(delta+0.952);
00233 else
00234 screenVal = 20.209 - delta*(1.930 + 0.086*delta);
00235
00236 return screenVal;
00237 }
00238
00239 inline G4double G4PairProductionRelModel::ScreenFunction1(G4double ScreenVariable)
00240
00241
00242
00243 {
00244 G4double screenVal;
00245
00246 if (ScreenVariable > 1.)
00247 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
00248 else
00249 screenVal = 42.392 - ScreenVariable*(7.796 - 1.961*ScreenVariable);
00250
00251 return screenVal;
00252 }
00253
00254
00255
00256 inline G4double G4PairProductionRelModel::ScreenFunction2(G4double ScreenVariable)
00257
00258
00259
00260 {
00261 G4double screenVal;
00262
00263 if (ScreenVariable > 1.)
00264 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
00265 else
00266 screenVal = 41.405 - ScreenVariable*(5.828 - 0.8945*ScreenVariable);
00267
00268 return screenVal;
00269 }
00270
00271
00272
00273
00274 inline G4double G4PairProductionRelModel::DeltaMax() const
00275 {
00276
00277 G4double FZ = 8.*(lnZ/3. + fCoulomb);
00278 return std::exp( (42.24-FZ)/8.368 ) + 0.952;
00279 }
00280
00281 inline G4double G4PairProductionRelModel::DeltaMin(G4double k) const
00282 {
00283 return 4.*136./z13*(CLHEP::electron_mass_c2/k);
00284 }
00285
00286
00287
00288 #endif