| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 //
32 // File name: G4RiGeMuPairProductionModel
33 //
34 // Authors: Girardo Depaola & Ricardo Pacheco
35 //
36 // Creation date: 29.10.2024
37 //
38 //
39 // -------------------------------------------------------------------
40 //
41
42 #ifndef G4RiGeMuPairProductionModel_h
43 #define G4RiGeMuPairProductionModel_h 1
44
45 #include "G4VEmModel.hh"
46 #include "G4NistManager.hh"
47 #include "G4ElementData.hh"
48 #include "G4Physics2DVector.hh"
49 #include "G4VEmAngularDistribution.hh"
50 #include <vector>
51
52 class G4Element;
53 class G4ParticleChangeForLoss;
54 class G4RiGeAngularGenerator;
55
56 class G4RiGeMuPairProductionModel : public G4VEmModel
57 {
58 public:
59
60 explicit G4RiGeMuPairProductionModel(const G4ParticleDefinition* p = nullptr);
61
62 ~G4RiGeMuPairProductionModel() override = default;
63
64 void Initialise(const G4ParticleDefinition*, const G4DataVector&) override;
65
66 void InitialiseLocal(const G4ParticleDefinition*,
67 G4VEmModel* masterModel) override;
68
69 G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
70 G4double kineticEnergy,
71 G4double Z, G4double A,
72 G4double cutEnergy,
73 G4double maxEnergy) override;
74
75 G4double ComputeDEDXPerVolume(const G4Material*,
76 const G4ParticleDefinition*,
77 G4double kineticEnergy,
78 G4double cutEnergy) override;
79
80 void SampleSecondaries(std::vector<G4DynamicParticle*>*,
81 const G4MaterialCutsCouple*,
82 const G4DynamicParticle*,
83 G4double tmin,
84 G4double maxEnergy) override;
85
86 G4double MinPrimaryEnergy(const G4Material*,
87 const G4ParticleDefinition*,
88 G4double) override;
89
90 G4double
91 ComputeDMicroscopicCrossSection(G4double tkin, G4double Z,
92 G4double pairEnergy);
93
94 inline void SetLowestKineticEnergy(G4double e);
95
96 inline void SetParticle(const G4ParticleDefinition*);
97
98 // hide assignment operator and copy constructor
99 G4RiGeMuPairProductionModel& operator=
100 (const G4RiGeMuPairProductionModel& right) = delete;
101 G4RiGeMuPairProductionModel(const G4RiGeMuPairProductionModel&) = delete;
102
103 protected:
104
105 G4double ComputMuPairLoss(G4double Z, G4double tkin, G4double cut,
106 G4double tmax);
107
108 G4double ComputeMicroscopicCrossSection(G4double tkin,
109 G4double Z,
110 G4double cut);
111
112 G4double FindScaledEnergy(G4int Z, G4double rand, G4double logTkin,
113 G4double yymin, G4double yymax);
114
115 inline G4double MaxSecondaryEnergyForElement(G4double kineticEnergy,
116 G4double Z);
117
118 void MakeSamplingTables();
119
120 void StoreTables() const;
121
122 G4bool RetrieveTables();
123
124 virtual void DataCorrupted(G4int Z, G4double logTkin) const;
125
126 G4ParticleChangeForLoss* fParticleChange = nullptr;
127 const G4ParticleDefinition* particle = nullptr;
128 G4NistManager* nist = nullptr;
129
130 G4double factorForCross;
131 G4double sqrte;
132 G4double particleMass = 0.0;
133 G4double z13 = 0.0;
134 G4double z23 = 0.0;
135 G4double lnZ = 0.0;
136
137 G4double minPairEnergy;
138 G4double lowestKinEnergy;
139
140 G4double emin;
141 G4double emax;
142 G4double ymin = -5.0;
143 G4double dy = 0.005;
144
145 // Random numbers for sampling
146 G4double randNumbs[9];
147
148 G4int currentZ = 0;
149 G4int nYBinPerDecade = 4;
150 std::size_t nbiny = 1000;
151 std::size_t nbine = 0;
152
153 G4bool fTableToFile = false;
154
155 // static members
156 static const G4int NZDATPAIR = 5;
157 static const G4int NINTPAIR = 8;
158 static const G4int ZDATPAIR[NZDATPAIR];
159 static const G4double xgi[NINTPAIR];
160 static const G4double wgi[NINTPAIR];
161
162 private:
163
164 G4RiGeAngularGenerator* fAngularGenerator;
165 G4ParticleDefinition* theElectron;
166 G4ParticleDefinition* thePositron;
167 G4String dataName{""};
168 };
169
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171
172 inline void G4RiGeMuPairProductionModel::SetLowestKineticEnergy(G4double e)
173 {
174 lowestKinEnergy = e;
175 }
176
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178
179 inline
180 void G4RiGeMuPairProductionModel::SetParticle(const G4ParticleDefinition* p)
181 {
182 if(nullptr == particle) {
183 particle = p;
184 particleMass = particle->GetPDGMass();
185 }
186 }
187
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
189
190 inline G4double
191 G4RiGeMuPairProductionModel::MaxSecondaryEnergyForElement(G4double kineticEnergy,
192 G4double ZZ)
193 {
194 G4int Z = G4lrint(ZZ);
195 if(Z != currentZ) {
196 currentZ = Z;
197 z13 = nist->GetZ13(Z);
198 z23 = z13*z13;
199 lnZ = nist->GetLOGZ(Z);
200 }
201 return kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
202 }
203
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
205
206 #endif
207