| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4MuonToMuonPairProductionMo
33 //
34 // Author: Siddharth Yajaman on the bas
35 //
36 // Creation date: 12.07.2022
37 //
38 //
39 // -------------------------------------------
40 //
41 //....oooOO0OOooo........oooOO0OOooo........oo
42 //....oooOO0OOooo........oooOO0OOooo........oo
43
44 #include "G4MuonToMuonPairProductionModel.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4EmParameters.hh"
48 #include "G4MuonMinus.hh"
49 #include "G4MuonPlus.hh"
50 #include "Randomize.hh"
51 #include "G4Material.hh"
52 #include "G4Element.hh"
53 #include "G4ParticleChangeForLoss.hh"
54 #include "G4Log.hh"
55 #include "G4Exp.hh"
56 #include <iostream>
57 #include <fstream>
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
60
61 G4MuonToMuonPairProductionModel::G4MuonToMuonP
62 const G4Parti
63 const G4Strin
64 : G4MuPairProductionModel(p, nam)
65 {
66 theMuonMinus = G4MuonMinus::MuonMinus();
67 theMuonPlus = G4MuonPlus::MuonPlus();
68 muonMass = theMuonPlus->GetPDGMass();
69 minPairEnergy = 2.*muonMass;
70 mueRatio = muonMass/CLHEP::electron_mass_c2;
71 factorForCross = 2./(3*CLHEP::pi)*
72 std::pow(CLHEP::fine_structure_const*CLHEP
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oo
76
77 G4double G4MuonToMuonPairProductionModel::Comp
78 G4d
79 G4d
80 G4d
81 // Calculates the differential (D) microscopi
82 // using the cross section formula of Kelner,
83 // Code written by Siddharth Yajaman (12/07/20
84 {
85 if (pairEnergy <= minPairEnergy)
86 return 0.0;
87
88 G4double totalEnergy = tkin + particleMass;
89 G4double residEnergy = totalEnergy - pairEne
90
91 if (residEnergy <= muonMass) { return 0.0; }
92
93 G4double a0 = 1.0 / (totalEnergy * residEner
94 G4double rhomax = 1.0 - 2*muonMass/pairEnerg
95 G4double tmnexp = 1. - rhomax;
96
97 if(tmnexp >= 1.0) { return 0.0; }
98
99 G4double tmn = G4Log(tmnexp);
100
101 G4double z2 = Z*Z;
102 G4double beta = 0.5*pairEnergy*pairEnergy*a0
103 G4double xi0 = 0.5*beta;
104
105 // Gaussian integration in ln(1-ro) ( with N
106 G4double rho[NINTPAIR];
107 G4double rho2[NINTPAIR];
108 G4double xi[NINTPAIR];
109 G4double xi1[NINTPAIR];
110 G4double xii[NINTPAIR];
111
112 for (G4int i = 0; i < NINTPAIR; ++i)
113 {
114 rho[i] = G4Exp(tmn*xgi[i]) - 1.0; // rho =
115 rho2[i] = rho[i] * rho[i];
116 xi[i] = xi0*(1.0-rho2[i]);
117 xi1[i] = 1.0 + xi[i];
118 xii[i] = 1.0 / xi[i];
119 }
120
121 G4double ximax = xi0*(1. - rhomax*rhomax);
122
123 G4double Y = 10 * std::sqrt(particleMass/tot
124 G4double U[8];
125
126 for (G4int i = 0; i < NINTPAIR; ++i)
127 {
128 U[i] = U_func(Z, rho2[i], xi[i], Y, pairEn
129 }
130
131 G4double UMax = U_func(Z, rhomax*rhomax, xim
132
133 G4double sum = 0.0;
134 for (G4int i = 0; i < NINTPAIR; ++i)
135 {
136 G4double X = 1 + U[i] - UMax;
137 G4double lnX = G4Log(X);
138 G4double phi = ((2 + rho2[i])*(1 + beta) +
139 G4Log(1 + xii[i]) - 1 - 3*
140 + ((1 + rho2[i])*(1 + 1.5*
141 *(1 - rho2[i]))*G4Log(xi1[
142 sum += wgi[i]*(1.0 + rho[i])*phi*lnX;
143 }
144
145 return -tmn*sum*factorForCross*z2*residEnerg
146
147 }
148
149 G4double G4MuonToMuonPairProductionModel::U_fu
150
151
152
153 {
154 G4int Z = G4lrint(ZZ);
155 G4double A27 = nist->GetA27(Z);
156 G4double Z13 = nist->GetZ13(Z);
157 static const G4double sqe = std::sqrt(G4Exp(
158 G4double res = (0.65 * B / (A27*Z13) * mueRa
159 (1 + (2*sqe*muonMass*muonMass*(B/Z13)*(1 +
160 /(CLHEP::electron_mass_c2*pairEnergy*(1 -
161 return res;
162 }
163
164 //....oooOO0OOooo........oooOO0OOooo........oo
165
166 void G4MuonToMuonPairProductionModel::SampleSe
167 std::vector<G4Dy
168 const G4Material
169 const G4DynamicP
170 G4double tmin,
171 G4double tmax)
172 {
173 G4double kinEnergy = aDynamicParticle->GetKi
174 //G4cout << "--- G4MuonToMuonPairProductionM
175 // << kinEnergy << " "
176 // << aDynamicParticle->GetDefinitio
177 G4double totalMomentum = std::sqrt(kinEnergy
178
179 G4ThreeVector partDirection = aDynamicPartic
180
181 // select randomly one element constituing t
182 const G4Element* anElement = SelectRandomAto
183
184 // define interval of energy transfer
185 G4double maxPairEnergy = MaxSecondaryEnergyF
186
187 G4double maxEnergy = std::min(tmax, maxPairE
188 G4double minEnergy = std::max(tmin, minPairE
189
190 if(minEnergy >= maxEnergy) { return; }
191 //G4cout << "emin= " << minEnergy << " emax=
192 // << " minPair= " << minPairEnergy << " max
193 // << " ymin= " << ymin << " dy= " << dy
194
195 G4double coeff = G4Log(minPairEnergy/kinEner
196
197 // compute limits
198 G4double yymin = G4Log(minEnergy/kinEnergy)/
199 G4double yymax = G4Log(maxEnergy/kinEnergy)/
200
201 //G4cout << "yymin= " << yymin << " yymax=
202
203 // units should not be used, bacause table w
204 G4double logTkin = G4Log(kinEnergy/CLHEP::Me
205
206 // sample mu-mu+ pair energy first
207
208 // select sample table via Z
209 G4int iz1(0), iz2(0);
210 for(G4int iz=0; iz<NZDATPAIR; ++iz) {
211 if(currentZ == ZDATPAIR[iz]) {
212 iz1 = iz2 = iz;
213 break;
214 } else if(currentZ < ZDATPAIR[iz]) {
215 iz2 = iz;
216 if(iz > 0) { iz1 = iz-1; }
217 else { iz1 = iz2; }
218 break;
219 }
220 }
221 if(0 == iz1) { iz1 = iz2 = NZDATPAIR-1; }
222
223 G4double pairEnergy = 0.0;
224 G4int count = 0;
225 //G4cout << "start loop Z1= " << iz1 << " Z2
226 do {
227 ++count;
228 // sampling using only one random number
229 G4double rand = G4UniformRand();
230
231 G4double x = FindScaledEnergy(iz1, rand, l
232 if(iz1 != iz2) {
233 G4double x2 = FindScaledEnergy(iz2, rand
234 G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1
235 G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2
236 //G4cout << count << ". x= " << x << "
237 // << " Z1= " << iz1 << " Z2
238 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
239 }
240 //G4cout << "x= " << x << " coeff= " << c
241 pairEnergy = kinEnergy*G4Exp(x*coeff);
242
243 // Loop checking, 03-Aug-2015, Vladimir Iv
244 } while((pairEnergy < minEnergy || pairEnerg
245
246 //G4cout << "## pairEnergy(GeV)= " << pairEn
247 // << " Etot(GeV)= " << totalEnergy/
248
249 // sample r=(mu+mu-)/pairEnergy ( uniformly
250 G4double rmax = 1 - 2*muonMass/(pairEnergy);
251 G4double r = rmax * (-1.+2.*G4UniformRand())
252
253 // compute energies from pairEnergy,r
254 G4double muMinusEnergy = (1.-r)*pairEnergy*0
255 G4double muPlusEnergy = pairEnergy - muMinus
256
257 // Sample angles
258 G4ThreeVector muMinusDirection, muPlusDirect
259 //
260 GetAngularDistribution()->SamplePairDirectio
261
262
263 // create G4DynamicParticle object for mu+mu
264 muMinusEnergy = std::max(muMinusEnergy - muo
265 muPlusEnergy = std::max(muPlusEnergy - muonM
266 G4DynamicParticle* aParticle1 =
267 new G4DynamicParticle(theMuonMinus,muMinus
268 G4DynamicParticle* aParticle2 =
269 new G4DynamicParticle(theMuonPlus,muPlusDi
270 // Fill output vector
271 vdp->push_back(aParticle1);
272 vdp->push_back(aParticle2);
273
274 // primary change
275 kinEnergy -= pairEnergy;
276 partDirection *= totalMomentum;
277 partDirection -= (aParticle1->GetMomentum()
278 partDirection = partDirection.unit();
279
280 // if energy transfer is higher than thresho
281 // then stop tracking the primary particle a
282 if (pairEnergy > SecondaryThreshold()) {
283 fParticleChange->ProposeTrackStatus(fStopA
284 fParticleChange->SetProposedKineticEnergy(
285 G4DynamicParticle* newdp =
286 new G4DynamicParticle(particle, partDire
287 vdp->push_back(newdp);
288 } else { // continue tracking the primary e-
289 fParticleChange->SetProposedMomentumDirect
290 fParticleChange->SetProposedKineticEnergy(
291 }
292 //G4cout << "--- G4MuonToMuonPairProductionM
293 }
294
295 //....oooOO0OOooo........oooOO0OOooo........oo
296
297 void G4MuonToMuonPairProductionModel::DataCorr
298 {
299 G4ExceptionDescription ed;
300 ed << "G4ElementData is not properly initial
301 << " Ekin(MeV)= " << G4Exp(logTkin)
302 << " IsMasterThread= " << IsMaster()
303 << " Model " << GetName();
304 G4Exception("G4MuonToMuonPairProductionModel
305 }
306
307 //....oooOO0OOooo........oooOO0OOooo........oo
308