include/Geant4/G4PairProductionRelModel.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // -------------------------------------------------------------------
0028 //
0029 // GEANT4 Class header file
0030 //
0031 //
0032 // File name:     G4PairProductionRelModel
0033 //
0034 // Author:        Andreas Schaelicke
0035 //
0036 // Creation date: 02.04.2009
0037 //
0038 // Modifications:
0039 // 28-05-18 New version with improved screening function approximation, improved
0040 //          LPM function approximation, efficiency, documentation and cleanup.
0041 //          Corrected call to selecting target atom in the final state sampling.
0042 //          (M. Novak)
0043 //
0044 // Class Description:
0045 //
0046 // Implementation of gamma conversion to e+e- in the field of a nucleus
0047 // relativistic approximation
0048 //
0049
0050 // -------------------------------------------------------------------
0051 //
0052
0053 #ifndef G4PairProductionRelModel_h
0054 #define G4PairProductionRelModel_h 1
0055
0056 #include <CLHEP/Units/PhysicalConstants.h>
0057
0058 #include "G4VEmModel.hh"
0059 #include "G4Log.hh"
0060 #include <vector>
0061
0062 class G4ParticleChangeForGamma;
0063 class G4Pow;
0064
0065 class G4PairProductionRelModel : public G4VEmModel
0066 {
0067
0068 public:
0069
0070   explicit G4PairProductionRelModel(const G4ParticleDefinition* p = nullptr,
0071                                     const G4String& nam = "BetheHeitlerLPM");
0072
0073   ~G4PairProductionRelModel() override;
0074
0075   void Initialise(const G4ParticleDefinition*, const G4DataVector&) override;
0076
0077   void InitialiseLocal(const G4ParticleDefinition*,
0078                G4VEmModel* masterModel) override;
0079
0080   G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
0081                       G4double kinEnergy,
0082                       G4double Z,
0083                       G4double A=0.,
0084                       G4double cut=0.,
0085                       G4double emax=DBL_MAX) override;
0086
0087   void SampleSecondaries(std::vector<G4DynamicParticle*>*,
0088              const G4MaterialCutsCouple*,
0089              const G4DynamicParticle*,
0090              G4double tmin,
0091              G4double maxEnergy) override;
0092
0093   void SetupForMaterial(const G4ParticleDefinition*,
0094             const G4Material*,G4double) override;
0095
0096   inline void   SetLPMflag(G4bool val) { fIsUseLPMCorrection = val;  }
0097   inline G4bool LPMflag() const        { return fIsUseLPMCorrection; }
0098
0099   G4PairProductionRelModel & operator=
0100   (const G4PairProductionRelModel &right) = delete;
0101   G4PairProductionRelModel(const  G4PairProductionRelModel&) = delete;
0102
0103 protected:
0104
0105   // for evaluating screening related functions
0106   inline void ComputePhi12(const G4double delta,
0107                G4double &phi1, G4double &phi2);
0108   inline G4double ScreenFunction1(const G4double delta);
0109   inline G4double ScreenFunction2(const G4double delta);
0110   inline void ScreenFunction12(const G4double delta,
0111                    G4double &f1, G4double &f2);
0112   // helper methods for cross-section computation under different approximations
0113   G4double ComputeParametrizedXSectionPerAtom(G4double gammaEnergy, G4double Z);
0114   G4double ComputeXSectionPerAtom(G4double gammaEnergy, G4double Z);
0115   G4double ComputeDXSectionPerAtom(G4double eplusEnergy, G4double gammaEnergy,
0116                                    G4double Z);
0117   G4double ComputeRelDXSectionPerAtom(G4double eplusEnergy,
0118                       G4double gammaEnergy, G4double Z);
0119
0120 private:
0121
0122   // for creating some data structure per Z
0123   void InitialiseElementData();
0124   struct ElementData {
0125     G4double  fLogZ13;
0126     G4double  fCoulomb;
0127     G4double  fLradEl;
0128     G4double  fDeltaFactor;
0129     G4double  fDeltaMaxLow;
0130     G4double  fDeltaMaxHigh;
0131     G4double  fEtaValue;
0132     G4double  fLPMVarS1Cond;
0133     G4double  fLPMILVarS1Cond;
0134   };
0135   // for precomputing comp. intensive parts of LPM suppression functions and
0136   // using them at run-time
0137   void InitLPMFunctions();
0138   void ComputeLPMGsPhis(G4double &funcGS, G4double &funcPhiS,
0139                         const G4double varShat);
0140   void GetLPMFunctions(G4double &lpmGs, G4double &lpmPhis, const G4double sval);
0141   void ComputeLPMfunctions(G4double &fXiS, G4double &fGS, G4double &fPhiS,
0142                            const G4double eps, const G4double egamma,
0143                            const G4int izet);
0144   struct LPMFuncs {
0145     LPMFuncs() : fIsInitialized(false), fISDelta(100.), fSLimit(2.) {}
0146     G4bool                 fIsInitialized;
0147     G4double               fISDelta;
0148     G4double               fSLimit;
0149     std::vector<G4double>  fLPMFuncG;
0150     std::vector<G4double>  fLPMFuncPhi;
0151   };
0152
0153 protected:
0154   static const G4int                gMaxZet;
0155   //
0156   static const G4double             gLPMconstant;
0157   //
0158   static const G4double             gXGL[8];
0159   static const G4double             gWGL[8];
0160   static const G4double             gFelLowZet[8];
0161   static const G4double             gFinelLowZet[8];
0162   //
0163   static const G4double             gXSecFactor;
0164   static const G4double             gEgLPMActivation;
0165   //
0166   static std::vector<ElementData*>  gElementData;
0167   static LPMFuncs                   gLPMFuncs;
0168   //
0169   G4bool isFirstInstance{false};
0170   G4bool                            fIsUseLPMCorrection;
0171   G4bool                            fIsUseCompleteScreening;
0172   //
0173   G4double                          fLPMEnergy;
0174   //
0175   G4double                          fParametrizedXSectionThreshold;
0176   G4double                          fCoulombCorrectionThreshold;
0177   //
0178   G4Pow*                            fG4Calc;
0179   G4ParticleDefinition*             fTheGamma;
0180   G4ParticleDefinition*             fTheElectron;
0181   G4ParticleDefinition*             fThePositron;
0182   G4ParticleChangeForGamma*         fParticleChange;
0183 };
0184 //
0185 // Bethe screening functions for the elastic (coherent) scattering:
0186 // Bethe's phi1, phi2 coherent screening functions were computed numerically
0187 // by using (the universal) atomic form factors computed based on the Thomas-
0188 // Fermi model of the atom (using numerical solution of the Thomas-Fermi
0189 // screening function instead of Moliere's analytical approximation). The
0190 // numerical results can be well approximated (better than Butcher & Messel
0191 // especially near the delta=1 limit) by:
0192 // ## if delta <= 1.4
0193 //  phi1(delta) = 20.806 - delta*(3.190 - 0.5710*delta)
0194 //  phi2(delta) = 20.234 - delta*(2.126 - 0.0903*delta)
0195 // ## if delta  > 1.4
0196 //  phi1(delta) = phi2(delta) = 21.0190 - 4.145*ln(delta + 0.958)
0197 // with delta = 136mc^2kZ^{-1/3}/[E(Eg-E)] = 136Z^{-1/3}eps0/[eps(1-eps)] where
0198 // Eg is the initial photon energy, E is the total energy transferred to one of
0199 // the e-/e+ pair, eps0 = mc^2/Eg and eps = E/Eg.
0200
0201 inline void G4PairProductionRelModel::ComputePhi12(const G4double delta,
0202                            G4double &phi1,
0203                            G4double &phi2)
0204 {
0205     if (delta > 1.4) {
0206       phi1 = 21.0190 - 4.145*G4Log(delta + 0.958);
0207       phi2 = phi1;
0208     } else {
0209       phi1 = 20.806 - delta*(3.190 - 0.5710*delta);
0210       phi2 = 20.234 - delta*(2.126 - 0.0903*delta);
0211     }
0212 }
0213
0214 // Compute the value of the screening function 3*PHI1(delta) - PHI2(delta):
0215 inline G4double G4PairProductionRelModel::ScreenFunction1(const G4double delta)
0216 {
0217   return (delta > 1.4) ? 42.038 - 8.29*G4Log(delta + 0.958)
0218                        : 42.184 - delta*(7.444 - 1.623*delta);
0219 }
0220
0221 // Compute the value of the screening function 1.5*PHI1(delta) +0.5*PHI2(delta):
0222 inline G4double G4PairProductionRelModel::ScreenFunction2(const G4double delta)
0223 {
0224   return (delta > 1.4) ? 42.038 - 8.29*G4Log(delta + 0.958)
0225                        : 41.326 - delta*(5.848 - 0.902*delta);
0226 }
0227
0228 // Same as ScreenFunction1 and ScreenFunction2 but computes them at once
0229 inline void G4PairProductionRelModel::ScreenFunction12(const G4double delta,
0230                                                      G4double &f1, G4double &f2)
0231 {
0232   if (delta > 1.4) {
0233     f1 = 42.038 - 8.29*G4Log(delta + 0.958);
0234     f2 = f1;
0235   } else {
0236     f1 = 42.184 - delta*(7.444 - 1.623*delta);
0237     f2 = 41.326 - delta*(5.848 - 0.902*delta);
0238   }
0239 }
0240
0241 #endif