include/Geant4/G4CoherentPairProduction.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 // Author:      Alexei Sytov
0027 // Co-author:   Gianfranco Paterno (testing)
0028 // Using the key points of G4BaierKatkov and developments of V.V. Tikhomirov,
0029 // partially described in L. Bandiera et al. Eur. Phys. J. C 82, 699 (2022)
0030
0031 #ifndef G4CoherentPairProduction_h
0032 #define G4CoherentPairProduction_h 1
0033
0034 #include "G4VDiscreteProcess.hh"
0035
0036 #include <vector>
0037 #include <CLHEP/Units/SystemOfUnits.h>
0038 #include <CLHEP/Units/PhysicalConstants.h>
0039 #include <CLHEP/Vector/TwoVector.h>
0040
0041 #include "G4ChannelingFastSimCrystalData.hh"
0042 #include "G4LogicalVolume.hh"
0043 #include "G4ParticleTable.hh"
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 class G4CoherentPairProduction : public G4VDiscreteProcess
0048 {
0049 public:
0050     G4CoherentPairProduction(const G4String& processName = "cpp",
0051                              G4ProcessType aType = fElectromagnetic);
0052
0053     ~G4CoherentPairProduction() = default;
0054
0055     G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&) override;
0056
0057     G4bool IsApplicable(const G4ParticleDefinition& aPD) override
0058     {
0059         return(aPD.GetParticleName() == "gamma");
0060     }
0061
0062     // print documentation in html format
0063     void ProcessDescription(std::ostream&) const override;
0064
0065     ///special functions
0066     void Input(const G4Material* crystal,
0067                const G4String &lattice)
0068     {Input(crystal,lattice,"");}
0069
0070     void Input(const G4Material* crystal,
0071                const G4String &lattice,
0072                const G4String &filePath);
0073
0074     // an option to use crystal data already created outside this class
0075     void Input(const G4ChannelingFastSimCrystalData* crystalData);
0076
0077     ///activate incoherent scattering
0078     ///(standard gamma conversion should be switched off in physics list)
0079     void ActivateIncoherentScattering(){fIncoherentScattering = true;}
0080
0081     G4ChannelingFastSimCrystalData* GetCrystalData() {return fCrystalData;}
0082
0083     ///get cuts
0084     // minimal energy for non-zero cross section
0085     G4double ModelMinPrimaryEnergy() { return fLowEnergyLimit;}
0086     G4double GetHighAngleLimit() {return fHighAngleLimit;}
0087     G4double GetPPKineticEnergyCut() {return fPPKineticEnergyCut;}
0088
0089     /// get the number of pairs in sampling of Baier-Katkov Integral
0090     /// (MC integration by e+- energy and angles <=> e+- momentum)
0091     G4int GetSamplingPairsNumber(){return fNMCPairs;}
0092
0093     /// get the number of particle angles 1/gamma in pair production
0094     /// defining the width of the angular distribution of pair sampling
0095     /// in the Baier-Katkov Integral
0096     G4double GetChargeParticleAngleFactor(){return fChargeParticleAngleFactor;}
0097
0098     /// get number of trajectory steps of a single particle (e- or e+)
0099     G4double GetNTrajectorySteps(){return fNTrajectorySteps;}
0100
0101     /// get effective radiation length
0102     /// (due to coherent process of pair production)
0103     /// simulated for the current photon
0104     G4double GetEffectiveLrad(){return fEffectiveLrad;}
0105
0106     ///get the name of G4Region in which the model is applicable
0107     G4String GetG4RegionName() {return fG4RegionName;}
0108
0109     ///set cuts
0110     void SetLowEnergyLimit(G4double energy){fLowEnergyLimit=energy;}
0111     void SetHighAngleLimit(G4double angle) {fHighAngleLimit=angle;}
0112     void SetPPKineticEnergyCut(G4double kineticEnergyCut) {fPPKineticEnergyCut=kineticEnergyCut;}
0113
0114     /// set the number of pairs in sampling of Baier-Katkov Integral
0115     /// (MC integration by e+- energy and angles <=> e+- momentum)
0116     void SetSamplingPairsNumber(G4int nPairs){fNMCPairs = nPairs;}
0117
0118     /// set the number of particle angles 1/gamma in pair production
0119     /// defining the width of the angular distribution of pair sampling
0120     /// in the Baier-Katkov Integral
0121     void SetChargeParticleAngleFactor(G4double chargeParticleAngleFactor)
0122     {fChargeParticleAngleFactor = chargeParticleAngleFactor;}
0123
0124     /// set number of trajectory steps of a single particle (e- or e+)
0125     void SetNTrajectorySteps(G4int nTrajectorySteps)
0126     {fNTrajectorySteps = nTrajectorySteps;}
0127
0128     ///set the name of G4Region in which the model is applicable
0129     void SetG4RegionName(const G4String& nameG4Region){fG4RegionName=nameG4Region;}
0130
0131     G4double GetMeanFreePath(const G4Track& aTrack,
0132                              G4double,
0133                              G4ForceCondition* condition) override;
0134
0135 private:
0136
0137     G4int FindVectorIndex(std::vector<G4double> &myvector, G4double value);
0138
0139     G4ChannelingFastSimCrystalData* fCrystalData{nullptr};
0140
0141     //collection of etotal
0142     std::vector <CLHEP::Hep2Vector> fullVectorEtotal;
0143
0144     //collection of x
0145     std::vector <CLHEP::Hep2Vector> fullVectorX;
0146
0147     //collection of y
0148     std::vector <CLHEP::Hep2Vector> fullVectorY;
0149
0150     //collection of tx
0151     std::vector <CLHEP::Hep2Vector> fullVectorTX;
0152
0153     //collection of tx
0154     std::vector <CLHEP::Hep2Vector> fullVectorTY;
0155
0156     //the vector of the discrete CDF of the production of sampling e+e- pairs
0157     //(in reality per distance along the photon direction)
0158     std::vector <G4double> fPairProductionCDFdz;
0159
0160     G4double fLowEnergyLimit = 1*CLHEP::GeV;
0161     G4double fHighAngleLimit = 50*CLHEP::mrad;
0162
0163     ///minimal kinetic energy of a charged particle produced
0164     G4double fPPKineticEnergyCut = 1*CLHEP::MeV;
0165
0166     ///Monte Carlo statistics of e+- pair sampling in Baier-Katkov for 1 photon
0167     G4int fNMCPairs = 150;
0168
0169     G4double fChargeParticleAngleFactor = 4; // number of particle angles 1/gamma:
0170         // more fChargeParticleAngleFactor => higher paramParticleAngle
0171
0172     ///number of trajectory steps of a single particle (e- or e+)
0173     G4int fNTrajectorySteps=250;
0174
0175     ///effective radiation length (due to coherent process of pair production)
0176     G4double fEffectiveLrad = 0.;
0177
0178     ///the name of G4Region in which the model is applicable
0179     G4String fG4RegionName = "Crystal";
0180
0181     ///charged particle mass
0182     const G4double fMass = CLHEP::electron_mass_c2;
0183
0184     ///flag of simulation of incoherent scattering
0185     G4bool fIncoherentScattering = false;
0186
0187 };
0188
0189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0190
0191 #endif
0192