include/Geant4/G4MuPairProductionModel.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:     G4MuPairProductionModel
0033 //
0034 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
0035 //
0036 // Creation date: 18.05.2002
0037 //
0038 // Modifications:
0039 //
0040 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
0041 // 27-01-03 Make models region aware (V.Ivanchenko)
0042 // 13-02-03 Add name (V.Ivanchenko)
0043 // 10-02-04 Update parameterisation using R.Kokoulin model (V.Ivanchenko)
0044 // 10-02-04 Add lowestKinEnergy (V.Ivanchenko)
0045 // 13-02-06 Add ComputeCrossSectionPerAtom (mma)
0046 // 12-05-06 Add parameter to SelectRandomAtom (A.Bogdanov)
0047 // 11-10-07 Add ignoreCut flag (V.Ivanchenko)
0048 // 28-02-08 Reorganized protected methods and members (V.Ivanchenko)
0049 //
0050 // Class Description:
0051 //
0052 // Implementation of e+e- pair production by muons
0053 // A.G. Bogdanov et al., IEEE Trans. Nuc. Sci., Vol.53, No.2, 2006
0054 // Base class for all pair production models for muons and hadrons.
0055 //
0056 // -------------------------------------------------------------------
0057 //
0058
0059 #ifndef G4MuPairProductionModel_h
0060 #define G4MuPairProductionModel_h 1
0061
0062 #include "G4VEmModel.hh"
0063 #include "G4NistManager.hh"
0064 #include "G4ElementData.hh"
0065 #include "G4Physics2DVector.hh"
0066 #include <vector>
0067
0068 class G4Element;
0069 class G4ParticleChangeForLoss;
0070 class G4ParticleChangeForGamma;
0071
0072 class G4MuPairProductionModel : public G4VEmModel
0073 {
0074 public:
0075
0076   explicit G4MuPairProductionModel(const G4ParticleDefinition* p = nullptr,
0077                                    const G4String& nam = "muPairProd");
0078
0079   ~G4MuPairProductionModel() = default;
0080
0081   void Initialise(const G4ParticleDefinition*, const G4DataVector&) override;
0082
0083   void InitialiseLocal(const G4ParticleDefinition*,
0084                        G4VEmModel* masterModel) override;
0085
0086   G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
0087                       G4double kineticEnergy,
0088                       G4double Z, G4double A,
0089                       G4double cutEnergy,
0090                       G4double maxEnergy) override;
0091
0092   G4double ComputeDEDXPerVolume(const G4Material*,
0093                                 const G4ParticleDefinition*,
0094                                 G4double kineticEnergy,
0095                                 G4double cutEnergy) override;
0096
0097   void SampleSecondaries(std::vector<G4DynamicParticle*>*,
0098                 const G4MaterialCutsCouple*,
0099                 const G4DynamicParticle*,
0100                 G4double tmin,
0101                 G4double maxEnergy) override;
0102
0103   G4double MinPrimaryEnergy(const G4Material*,
0104                             const G4ParticleDefinition*,
0105                             G4double) override;
0106
0107   virtual G4double
0108   ComputeDMicroscopicCrossSection(G4double tkin, G4double Z,
0109                   G4double pairEnergy);
0110
0111   inline void SetLowestKineticEnergy(G4double e);
0112
0113   inline void SetParticle(const G4ParticleDefinition*);
0114
0115   // hide assignment operator and copy constructor
0116   G4MuPairProductionModel & operator=
0117   (const G4MuPairProductionModel &right) = delete;
0118   G4MuPairProductionModel(const  G4MuPairProductionModel&) = delete;
0119
0120 protected:
0121
0122   G4double ComputMuPairLoss(G4double Z, G4double tkin, G4double cut,
0123                             G4double tmax);
0124
0125   G4double ComputeMicroscopicCrossSection(G4double tkin,
0126                                           G4double Z,
0127                                           G4double cut);
0128
0129   G4double FindScaledEnergy(G4int Z, G4double rand, G4double logTkin,
0130                 G4double yymin, G4double yymax);
0131
0132   inline G4double MaxSecondaryEnergyForElement(G4double kineticEnergy,
0133                            G4double Z);
0134
0135   void MakeSamplingTables();
0136
0137   void StoreTables() const;
0138
0139   G4bool RetrieveTables();
0140
0141   virtual void DataCorrupted(G4int Z, G4double logTkin) const;
0142
0143   G4ParticleChangeForLoss* fParticleChange = nullptr;
0144   const G4ParticleDefinition* particle = nullptr;
0145   G4NistManager* nist = nullptr;
0146
0147   G4double factorForCross;
0148   G4double sqrte;
0149   G4double particleMass = 0.0;
0150   G4double z13 = 0.0;
0151   G4double z23 = 0.0;
0152   G4double lnZ = 0.0;
0153
0154   G4double minPairEnergy;
0155   G4double lowestKinEnergy;
0156
0157   G4double emin;
0158   G4double emax;
0159   G4double ymin = -5.0;
0160   G4double dy = 0.005;
0161
0162   G4int currentZ = 0;
0163   G4int nYBinPerDecade = 4;
0164   std::size_t nbiny = 1000;
0165   std::size_t nbine = 0;
0166
0167   G4bool fTableToFile = false;
0168
0169   // static members
0170   static const G4int NZDATPAIR = 5;
0171   static const G4int NINTPAIR = 8;
0172   static const G4int ZDATPAIR[NZDATPAIR];
0173   static const G4double xgi[NINTPAIR];
0174   static const G4double wgi[NINTPAIR];
0175
0176 private:
0177
0178   G4ParticleDefinition* theElectron;
0179   G4ParticleDefinition* thePositron;
0180   G4String dataName{""};
0181 };
0182
0183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0184
0185 inline void G4MuPairProductionModel::SetLowestKineticEnergy(G4double e)
0186 {
0187   lowestKinEnergy = e;
0188 }
0189
0190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0191
0192 inline
0193 void G4MuPairProductionModel::SetParticle(const G4ParticleDefinition* p)
0194 {
0195   if(nullptr == particle) {
0196     particle = p;
0197     particleMass = particle->GetPDGMass();
0198   }
0199 }
0200
0201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0202
0203 inline G4double
0204 G4MuPairProductionModel::MaxSecondaryEnergyForElement(G4double kineticEnergy,
0205                               G4double ZZ)
0206 {
0207   G4int Z = G4lrint(ZZ);
0208   if(Z != currentZ) {
0209     currentZ = Z;
0210     z13 = nist->GetZ13(Z);
0211     z23 = z13*z13;
0212     lnZ = nist->GetLOGZ(Z);
0213   }
0214   return kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
0215 }
0216
0217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0218
0219 #endif