| 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 file
30 //
31 //
32 // File name: G4MuonToMuonPairProductionModel
33 //
34 // Author: Siddharth Yajaman on the base of Vladimir Ivantchenko code
35 //
36 // Creation date: 12.07.2022
37 //
38 //
39 // -------------------------------------------------------------------
40 //
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
60
61 G4MuonToMuonPairProductionModel::G4MuonToMuonPairProductionModel(
62 const G4ParticleDefinition* p,
63 const G4String& nam)
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::classic_electr_radius/mueRatio, 2);
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76
77 G4double G4MuonToMuonPairProductionModel::ComputeDMicroscopicCrossSection(
78 G4double tkin,
79 G4double Z,
80 G4double pairEnergy)
81 // Calculates the differential (D) microscopic cross section
82 // using the cross section formula of Kelner, Kokoulin and Petrukhin (1999)
83 // Code written by Siddharth Yajaman (12/07/2022)
84 {
85 if (pairEnergy <= minPairEnergy)
86 return 0.0;
87
88 G4double totalEnergy = tkin + particleMass;
89 G4double residEnergy = totalEnergy - pairEnergy;
90
91 if (residEnergy <= muonMass) { return 0.0; }
92
93 G4double a0 = 1.0 / (totalEnergy * residEnergy);
94 G4double rhomax = 1.0 - 2*muonMass/pairEnergy;
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 NINTPAIR points)
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 = -asymmetry
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/totalEnergy);
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, pairEnergy);
129 }
130
131 G4double UMax = U_func(Z, rhomax*rhomax, ximax, Y, pairEnergy);
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) + xi[i]*(3 + rho2[i]))*
139 G4Log(1 + xii[i]) - 1 - 3*rho2[i] + beta*(1 - 2*rho2[i])
140 + ((1 + rho2[i])*(1 + 1.5*beta) - xii[i]*(1 + 2*beta)
141 *(1 - rho2[i]))*G4Log(xi1[i]);
142 sum += wgi[i]*(1.0 + rho[i])*phi*lnX;
143 }
144
145 return -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
146
147 }
148
149 G4double G4MuonToMuonPairProductionModel::U_func(G4double ZZ, G4double rho2,
150 G4double xi, G4double Y,
151 G4double pairEnergy,
152 const G4double B)
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(1.0));
158 G4double res = (0.65 * B / (A27*Z13) * mueRatio)/
159 (1 + (2*sqe*muonMass*muonMass*(B/Z13)*(1 + xi)*(1 + Y))
160 /(CLHEP::electron_mass_c2*pairEnergy*(1 - rho2)));
161 return res;
162 }
163
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165
166 void G4MuonToMuonPairProductionModel::SampleSecondaries(
167 std::vector<G4DynamicParticle*>* vdp,
168 const G4MaterialCutsCouple* couple,
169 const G4DynamicParticle* aDynamicParticle,
170 G4double tmin,
171 G4double tmax)
172 {
173 G4double kinEnergy = aDynamicParticle->GetKineticEnergy();
174 //G4cout << "--- G4MuonToMuonPairProductionModel::SampleSecondaries E(MeV)= "
175 // << kinEnergy << " "
176 // << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
177 G4double totalMomentum = std::sqrt(kinEnergy*(kinEnergy + 2.0*muonMass));
178
179 G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
180
181 // select randomly one element constituing the material
182 const G4Element* anElement = SelectRandomAtom(couple,particle,kinEnergy);
183
184 // define interval of energy transfer
185 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy,
186 anElement->GetZ());
187 G4double maxEnergy = std::min(tmax, maxPairEnergy);
188 G4double minEnergy = std::max(tmin, minPairEnergy);
189
190 if(minEnergy >= maxEnergy) { return; }
191 //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
192 // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
193 // << " ymin= " << ymin << " dy= " << dy << G4endl;
194
195 G4double coeff = G4Log(minPairEnergy/kinEnergy)/ymin;
196
197 // compute limits
198 G4double yymin = G4Log(minEnergy/kinEnergy)/coeff;
199 G4double yymax = G4Log(maxEnergy/kinEnergy)/coeff;
200
201 //G4cout << "yymin= " << yymin << " yymax= " << yymax << G4endl;
202
203 // units should not be used, bacause table was built without
204 G4double logTkin = G4Log(kinEnergy/CLHEP::MeV);
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= " << iz2 << G4endl;
226 do {
227 ++count;
228 // sampling using only one random number
229 G4double rand = G4UniformRand();
230
231 G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
232 if(iz1 != iz2) {
233 G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
234 G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1]);
235 G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2]);
236 //G4cout << count << ". x= " << x << " x2= " << x2
237 // << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
238 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
239 }
240 //G4cout << "x= " << x << " coeff= " << coeff << G4endl;
241 pairEnergy = kinEnergy*G4Exp(x*coeff);
242
243 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
244 } while((pairEnergy < minEnergy || pairEnergy > maxEnergy) && 10 > count);
245
246 //G4cout << "## pairEnergy(GeV)= " << pairEnergy/GeV
247 // << " Etot(GeV)= " << totalEnergy/GeV << G4endl;
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.5;
255 G4double muPlusEnergy = pairEnergy - muMinusEnergy;
256
257 // Sample angles
258 G4ThreeVector muMinusDirection, muPlusDirection;
259 //
260 GetAngularDistribution()->SamplePairDirections(aDynamicParticle,
261 muMinusEnergy, muPlusEnergy,
262 muMinusDirection, muPlusDirection);
263 // create G4DynamicParticle object for mu+mu-
264 muMinusEnergy = std::max(muMinusEnergy - muonMass, 0.0);
265 muPlusEnergy = std::max(muPlusEnergy - muonMass, 0.0);
266 G4DynamicParticle* aParticle1 =
267 new G4DynamicParticle(theMuonMinus,muMinusDirection,muMinusEnergy);
268 G4DynamicParticle* aParticle2 =
269 new G4DynamicParticle(theMuonPlus,muPlusDirection,muPlusEnergy);
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() + aParticle2->GetMomentum());
278 partDirection = partDirection.unit();
279
280 // if energy transfer is higher than threshold (very high by default)
281 // then stop tracking the primary particle and create a new secondary
282 if (pairEnergy > SecondaryThreshold()) {
283 fParticleChange->ProposeTrackStatus(fStopAndKill);
284 fParticleChange->SetProposedKineticEnergy(0.0);
285 G4DynamicParticle* newdp =
286 new G4DynamicParticle(particle, partDirection, kinEnergy);
287 vdp->push_back(newdp);
288 } else { // continue tracking the primary e-/e+ otherwise
289 fParticleChange->SetProposedMomentumDirection(partDirection);
290 fParticleChange->SetProposedKineticEnergy(kinEnergy);
291 }
292 //G4cout << "--- G4MuonToMuonPairProductionModel::SampleSecondaries done" << G4endl;
293 }
294
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
296
297 void G4MuonToMuonPairProductionModel::DataCorrupted(G4int Z, G4double logTkin) const
298 {
299 G4ExceptionDescription ed;
300 ed << "G4ElementData is not properly initialized Z= " << Z
301 << " Ekin(MeV)= " << G4Exp(logTkin)
302 << " IsMasterThread= " << IsMaster()
303 << " Model " << GetName();
304 G4Exception("G4MuonToMuonPairProductionModel","em0033",FatalException,ed,"");
305 }
306
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308