include/Geant4/G4RiGeMuPairProductionModel.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:     G4RiGeMuPairProductionModel
0033 //
0034 // Authors:       Girardo Depaola & Ricardo Pacheco
0035 //
0036 // Creation date: 29.10.2024
0037 //
0038 //
0039 // -------------------------------------------------------------------
0040 //
0041
0042 #ifndef G4RiGeMuPairProductionModel_h
0043 #define G4RiGeMuPairProductionModel_h 1
0044
0045 #include "G4VEmModel.hh"
0046 #include "G4NistManager.hh"
0047 #include "G4ElementData.hh"
0048 #include "G4Physics2DVector.hh"
0049 #include "G4VEmAngularDistribution.hh"
0050 #include <vector>
0051
0052 class G4Element;
0053 class G4ParticleChangeForLoss;
0054 class G4RiGeAngularGenerator;
0055
0056 class G4RiGeMuPairProductionModel : public G4VEmModel
0057 {
0058 public:
0059
0060   explicit G4RiGeMuPairProductionModel(const G4ParticleDefinition* p = nullptr);
0061
0062   ~G4RiGeMuPairProductionModel() override = default;
0063
0064   void Initialise(const G4ParticleDefinition*, const G4DataVector&) override;
0065
0066   void InitialiseLocal(const G4ParticleDefinition*,
0067                        G4VEmModel* masterModel) override;
0068
0069   G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
0070               G4double kineticEnergy,
0071               G4double Z, G4double A,
0072               G4double cutEnergy,
0073               G4double maxEnergy) override;
0074
0075   G4double ComputeDEDXPerVolume(const G4Material*,
0076                                 const G4ParticleDefinition*,
0077                                 G4double kineticEnergy,
0078                                 G4double cutEnergy) override;
0079
0080   void SampleSecondaries(std::vector<G4DynamicParticle*>*,
0081                          const G4MaterialCutsCouple*,
0082                          const G4DynamicParticle*,
0083                          G4double tmin,
0084                          G4double maxEnergy) override;
0085
0086   G4double MinPrimaryEnergy(const G4Material*,
0087                             const G4ParticleDefinition*,
0088                             G4double) override;
0089
0090   G4double
0091   ComputeDMicroscopicCrossSection(G4double tkin, G4double Z,
0092           G4double pairEnergy);
0093
0094   inline void SetLowestKineticEnergy(G4double e);
0095
0096   inline void SetParticle(const G4ParticleDefinition*);
0097
0098   // hide assignment operator and copy constructor
0099   G4RiGeMuPairProductionModel& operator=
0100   (const G4RiGeMuPairProductionModel& right) = delete;
0101   G4RiGeMuPairProductionModel(const G4RiGeMuPairProductionModel&) = delete;
0102
0103 protected:
0104
0105   G4double ComputMuPairLoss(G4double Z, G4double tkin, G4double cut,
0106                             G4double tmax);
0107
0108   G4double ComputeMicroscopicCrossSection(G4double tkin,
0109                                           G4double Z,
0110                                           G4double cut);
0111
0112   G4double FindScaledEnergy(G4int Z, G4double rand, G4double logTkin,
0113           G4double yymin, G4double yymax);
0114
0115   inline G4double MaxSecondaryEnergyForElement(G4double kineticEnergy,
0116                  G4double Z);
0117
0118   void MakeSamplingTables();
0119
0120   void StoreTables() const;
0121
0122   G4bool RetrieveTables();
0123
0124   virtual void DataCorrupted(G4int Z, G4double logTkin) const;
0125
0126   G4ParticleChangeForLoss* fParticleChange = nullptr;
0127   const G4ParticleDefinition* particle = nullptr;
0128   G4NistManager* nist = nullptr;
0129
0130   G4double factorForCross;
0131   G4double sqrte;
0132   G4double particleMass = 0.0;
0133   G4double z13 = 0.0;
0134   G4double z23 = 0.0;
0135   G4double lnZ = 0.0;
0136
0137   G4double minPairEnergy;
0138   G4double lowestKinEnergy;
0139
0140   G4double emin;
0141   G4double emax;
0142   G4double ymin = -5.0;
0143   G4double dy = 0.005;
0144
0145   // Random numbers for sampling
0146   G4double randNumbs[9];
0147
0148   G4int currentZ = 0;
0149   G4int nYBinPerDecade = 4;
0150   std::size_t nbiny = 1000;
0151   std::size_t nbine = 0;
0152
0153   G4bool fTableToFile = false;
0154
0155   // static members
0156   static const G4int NZDATPAIR = 5;
0157   static const G4int NINTPAIR = 8;
0158   static const G4int ZDATPAIR[NZDATPAIR];
0159   static const G4double xgi[NINTPAIR];
0160   static const G4double wgi[NINTPAIR];
0161
0162 private:
0163
0164   G4RiGeAngularGenerator* fAngularGenerator;
0165   G4ParticleDefinition* theElectron;
0166   G4ParticleDefinition* thePositron;
0167   G4String dataName{""};
0168 };
0169
0170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0171
0172 inline void G4RiGeMuPairProductionModel::SetLowestKineticEnergy(G4double e)
0173 {
0174   lowestKinEnergy = e;
0175 }
0176
0177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0178
0179 inline
0180 void G4RiGeMuPairProductionModel::SetParticle(const G4ParticleDefinition* p)
0181 {
0182   if(nullptr == particle) {
0183     particle = p;
0184     particleMass = particle->GetPDGMass();
0185   }
0186 }
0187
0188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0189
0190 inline G4double
0191 G4RiGeMuPairProductionModel::MaxSecondaryEnergyForElement(G4double kineticEnergy,
0192                          G4double ZZ)
0193 {
0194   G4int Z = G4lrint(ZZ);
0195   if(Z != currentZ) {
0196     currentZ = Z;
0197     z13 = nist->GetZ13(Z);
0198     z23 = z13*z13;
0199     lnZ = nist->GetLOGZ(Z);
0200   }
0201   return kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
0202 }
0203
0204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0205
0206 #endif