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   void ComputeLPMfunctions(G4double &fXiS, G4double &fGS, G4double &fPhiS,
0136                            const G4double eps, const G4double egamma,
0137                            const G4int izet);
0138
0139 protected:
0140   static const G4int                gMaxZet;
0141   //
0142   static const G4double             gLPMconstant;
0143   //
0144   static const G4double             gXGL[8];
0145   static const G4double             gWGL[8];
0146   static const G4double             gFelLowZet[8];
0147   static const G4double             gFinelLowZet[8];
0148   //
0149   static const G4double             gXSecFactor;
0150   static const G4double             gEgLPMActivation;
0151   //
0152   static std::vector<ElementData*>  gElementData;
0153   //
0154   G4bool isFirstInstance{false};
0155   G4bool                            fIsUseLPMCorrection;
0156   G4bool                            fIsUseCompleteScreening;
0157   //
0158   G4double                          fLPMEnergy;
0159   //
0160   G4double                          fParametrizedXSectionThreshold;
0161   G4double                          fCoulombCorrectionThreshold;
0162   //
0163   G4Pow*                            fG4Calc;
0164   G4ParticleDefinition*             fTheGamma;
0165   G4ParticleDefinition*             fTheElectron;
0166   G4ParticleDefinition*             fThePositron;
0167   G4ParticleChangeForGamma*         fParticleChange;
0168 };
0169 //
0170 // Bethe screening functions for the elastic (coherent) scattering:
0171 // Bethe's phi1, phi2 coherent screening functions were computed numerically
0172 // by using (the universal) atomic form factors computed based on the Thomas-
0173 // Fermi model of the atom (using numerical solution of the Thomas-Fermi
0174 // screening function instead of Moliere's analytical approximation). The
0175 // numerical results can be well approximated (better than Butcher & Messel
0176 // especially near the delta=1 limit) by:
0177 // ## if delta <= 1.4
0178 //  phi1(delta) = 20.806 - delta*(3.190 - 0.5710*delta)
0179 //  phi2(delta) = 20.234 - delta*(2.126 - 0.0903*delta)
0180 // ## if delta  > 1.4
0181 //  phi1(delta) = phi2(delta) = 21.0190 - 4.145*ln(delta + 0.958)
0182 // with delta = 136mc^2kZ^{-1/3}/[E(Eg-E)] = 136Z^{-1/3}eps0/[eps(1-eps)] where
0183 // Eg is the initial photon energy, E is the total energy transferred to one of
0184 // the e-/e+ pair, eps0 = mc^2/Eg and eps = E/Eg.
0185
0186 inline void G4PairProductionRelModel::ComputePhi12(const G4double delta,
0187                            G4double &phi1,
0188                            G4double &phi2)
0189 {
0190     if (delta > 1.4) {
0191       phi1 = 21.0190 - 4.145*G4Log(delta + 0.958);
0192       phi2 = phi1;
0193     } else {
0194       phi1 = 20.806 - delta*(3.190 - 0.5710*delta);
0195       phi2 = 20.234 - delta*(2.126 - 0.0903*delta);
0196     }
0197 }
0198
0199 // Compute the value of the screening function 3*PHI1(delta) - PHI2(delta):
0200 inline G4double G4PairProductionRelModel::ScreenFunction1(const G4double delta)
0201 {
0202   return (delta > 1.4) ? 42.038 - 8.29*G4Log(delta + 0.958)
0203                        : 42.184 - delta*(7.444 - 1.623*delta);
0204 }
0205
0206 // Compute the value of the screening function 1.5*PHI1(delta) +0.5*PHI2(delta):
0207 inline G4double G4PairProductionRelModel::ScreenFunction2(const G4double delta)
0208 {
0209   return (delta > 1.4) ? 42.038 - 8.29*G4Log(delta + 0.958)
0210                        : 41.326 - delta*(5.848 - 0.902*delta);
0211 }
0212
0213 // Same as ScreenFunction1 and ScreenFunction2 but computes them at once
0214 inline void G4PairProductionRelModel::ScreenFunction12(const G4double delta,
0215                                                      G4double &f1, G4double &f2)
0216 {
0217   if (delta > 1.4) {
0218     f1 = 42.038 - 8.29*G4Log(delta + 0.958);
0219     f2 = f1;
0220   } else {
0221     f1 = 42.184 - delta*(7.444 - 1.623*delta);
0222     f2 = 41.326 - delta*(5.848 - 0.902*delta);
0223   }
0224 }
0225
0226 #endif