| 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: G4RiGeAngularGenerator
33 //
34 // Authors: Girardo Depaola & Ricardo Pa
35 //
36 // Creation date: 27 October 2024
37 //
38 // -------------------------------------------
39 //
40
41 #include "G4RiGeAngularGenerator.hh"
42 #include "Randomize.hh"
43 #include "G4Log.hh"
44 #include "G4Exp.hh"
45 #include <CLHEP/Units/PhysicalConstants.h>
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 G4RiGeAngularGenerator::G4RiGeAngularGenerator
50 : G4VEmAngularDistribution("RiGeAngularGen")
51 {}
52
53 //....oooOO0OOooo........oooOO0OOooo........oo
54
55 G4ThreeVector&
56 G4RiGeAngularGenerator::SampleDirection(const
57 G4doub
58 const
59 {
60 // Sample gamma angle (Z - axis along the pa
61 G4double cost = SampleCosTheta(dp->GetKineti
62 dp->GetDefini
63 G4double sint = std::sqrt((1.0 - cost)*(1.0
64 G4double phi = CLHEP::twopi*G4UniformRand()
65
66 fLocalDirection.set(sint*std::cos(phi), sint
67 fLocalDirection.rotateUz(dp->GetMomentumDire
68
69 return fLocalDirection;
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oo
73
74 G4double G4RiGeAngularGenerator::SampleCosThet
75 G4double gEnergy,
76 G4double mass)
77 {
78 G4double gam = 1.0 + primKinEnergy/mass;
79 G4double rmax = gam*CLHEP::halfpi*std::min(1
80 G4double rmax2= rmax*rmax;
81 G4double x = G4UniformRand()*rmax2/(1.0 + rm
82
83 return std::cos(std::sqrt(x/(1.0 - x))/gam);
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oo
87
88 G4LorentzVector
89 G4RiGeAngularGenerator::Sample5DPairDirections
90 G4ThreeVector& dirElectron,
91 G4ThreeVector& dirPositron,
92 const G4double gEnergy, const
93 const G4double gMomentum,
94 const G4double muFinalMomentu
95 const G4double muFinalEnergy,
96 const G4double* randNumbs,
97 const G4double* W)
98 {
99 G4double muEnergy = dp->GetKineticEnergy();
100 G4ThreeVector muMomentumVector = dp->GetMome
101 G4double muMomentum = muMomentumVector.mag()
102 G4LorentzVector muFinalFourMomentum(muEnergy
103
104 // Electron mass
105 G4double eMass = CLHEP::electron_mass_c2;
106 G4double eMass2 = eMass*eMass;
107
108 // Muon mass
109 G4double muMass = dp->GetDefinition()->GetPD
110
111 G4double mint3 = 0.;
112 G4double maxt3 = CLHEP::pi;
113 G4double Cmin = std::cos(maxt3);
114 G4double Cmax = std::cos(mint3);
115
116 if (randNumbs[7] < W[0]) {
117 G4double A1 = -(q2 - 2.*muEnergy*gEnergy);
118 G4double B1 = -(2.*gMomentum*muMomentum);
119 G4double tginterval = G4Log((A1 + B1)/(A1
120
121 G4double costg = (-A1 + (A1 - B1)*G4Exp(B1
122 G4double sintg = std::sqrt((1.0 - costg)*(
123 G4double phig = CLHEP::twopi*randNumbs[2]
124 G4double sinpg = std::sin(phig);
125 G4double cospg = std::cos(phig);
126
127 G4ThreeVector dirGamma;
128 dirGamma.set(sintg*cospg, sintg*sinpg, cos
129 G4LorentzVector gFourMomentum(gEnergy, dir
130
131 G4double Ap = muMomentum*muMomentum + muFi
132 G4double A = Ap - 2.*muMomentum*gMomentum*
133 G4double B = 2.*muFinalMomentum*gMomentum*
134 G4double C = 2.*muFinalMomentum*gMomentum*
135 G4double absB = std::abs(B);
136 G4double t1interval = (1./(A + C + absB*mi
137 G4double t1 = (-(A + C) + 1./(1./(A + C +
138 G4double sint1 = std::sin(t1);
139 G4double cost1 = std::cos(t1);
140
141 G4ThreeVector dirMuon;
142 dirMuon.set(sint1, 0., cost1);
143
144 G4double cost5 = -1. + 2.*randNumbs[6];
145 G4double phi5 = CLHEP::twopi*randNumbs[8];
146
147 G4LorentzVector eFourMomentumMQ = eDP2(q2,
148 G4LorentzVector pFourMomentumMQ = pDP2(eMa
149
150 G4LorentzVector eFourMomentum = eFourMomen
151 G4LorentzVector pFourMomentum = pFourMomen
152
153 dirElectron = eFourMomentum.vect().unit();
154 dirPositron = pFourMomentum.vect().unit();
155
156 G4double phi = CLHEP::twopi*randNumbs[3];
157 PhiRotation(dirElectron, phi);
158 PhiRotation(dirPositron, phi);
159
160 } else if (randNumbs[7] >= W[0] && randNumbs
161 G4double A3 = q2 + 2.*gEnergy*muFinalEnerg
162 G4double B3 = -2.*gMomentum*muFinalMomentu
163
164 G4double tQ3interval = G4Log((A3 + B3)/(A3
165 G4double tQMG = (-A3 + (A3 - B3)*G4Exp(B3*
166 G4double phiQP = CLHEP::twopi*randNumbs[2]
167
168 G4double sintQ3 = std::sqrt(1. - tQMG*tQMG
169 G4double cospQP = std::cos(phiQP);
170 G4double sinpQP = std::sin(phiQP);
171
172 G4double Ap = muMomentum*muMomentum + muFi
173 G4double A = Ap + 2.*muFinalMomentum*gMome
174 G4double B = -2.*muMomentum*gMomentum*sint
175 G4double C = -2.*muMomentum*gMomentum*tQMG
176
177 G4double absB = std::abs(B);
178 G4double t3interval = (1./(A + C + absB*mi
179 G4double t3 = (-(A + C) + 1./(1./(A + C +
180 G4double sint3 = std::sin(t3);
181 G4double cost3 = std::cos(t3);
182
183 G4double cost = -sint3*sintQ3*cospQP + cos
184 G4double sint = std::sqrt((1. + cost)*(1.
185 G4double cosp = (sintQ3*cospQP*cost3 + sin
186 G4double sinp = sintQ3*sinpQP/sint;
187
188 G4ThreeVector dirGamma;
189 dirGamma.set(sint*cosp, sint*sinp, cost);
190 G4LorentzVector gFourMomentum(gEnergy, dir
191
192 G4ThreeVector dirMuon;
193 dirMuon.set(sint3, 0., cost3);
194
195 G4double cost5 = -1. + 2.*randNumbs[6];
196 G4double phi5 = CLHEP::twopi*randNumbs[8];
197
198 G4LorentzVector eFourMomentumMQ = eDP2(q2,
199 G4LorentzVector pFourMomentumMQ = pDP2(eMa
200
201 G4LorentzVector eFourMomentum = eFourMomen
202 G4LorentzVector pFourMomentum = pFourMomen
203
204 dirElectron = eFourMomentum.vect().unit();
205 dirPositron = pFourMomentum.vect().unit();
206
207 G4double phi = CLHEP::twopi*randNumbs[3];
208 PhiRotation(dirElectron, phi);
209 PhiRotation(dirPositron, phi);
210
211 } else if (randNumbs[7] >= W[1] && randNumbs
212 G4double phi3 = CLHEP::twopi*randNumbs[0];
213 G4double phi5 = CLHEP::twopi*randNumbs[1];
214 G4double phi6 = CLHEP::twopi*randNumbs[2];
215 G4double minmuFinalEnergy = muMass;
216 G4double muEnergyInterval = muEnergy - 2.*
217 G4double muFEnergy = minmuFinalEnergy + mu
218
219 G4double mineEnergy = eMass;
220 G4double maxeEnergy = muEnergy - muFEnergy
221 G4double eEnergyInterval = maxeEnergy - mi
222 G4double eEnergy = mineEnergy + eEnergyInt
223
224 G4double cosp3 = 1.;
225 G4double sinp3 = 0.;
226 G4double cosp5 = std::cos(phi5);
227 G4double sinp5 = std::sin(phi5);
228 G4double cosp6 = std::cos(phi6);
229 G4double sinp6 = std::sin(phi6);
230
231 G4double muFMomentum = std::sqrt(muFEnergy
232 G4double eMomentum = std::sqrt(eEnergy*eEn
233 G4double pEnergy = muEnergy - muFEnergy -
234 G4double pMomentum = std::sqrt(pEnergy*pEn
235
236 G4double A3 = -2.*muMass*muMass + 2.*muEne
237 G4double B3 = -2.*muMomentum*muFinalMoment
238 G4double cost3interval = G4Log((A3 + B3*Cm
239 G4double expanCost3r6 = G4Exp(B3*cost3inte
240 G4double cost3 = A3*(expanCost3r6 - 1.)/B3
241 G4double sint3 = std::sqrt((1. - cost3)*(1
242
243 G4ThreeVector muFinalMomentumVector(muFMom
244
245 G4LorentzVector muFourMomentum(muMomentum,
246 muFinalFourMomentum.set(muFEnergy, muFinal
247 G4LorentzVector auxVec1 = muFourMomentum -
248 G4double A5 = auxVec1.mag2() - 2.*eEnergy*
249 2.*muMomentumVector[2]*eMomentum - 2.*mu
250 G4double B5 = -2.*muFMomentum*eMomentum*(s
251 G4double absA5 = std::abs(A5);
252 G4double absB5 = std::abs(B5);
253 G4double mint5 = 0.;
254 G4double maxt5 = CLHEP::pi;
255 G4double t5interval = G4Log((absA5 + absB5
256 G4double argexp = absB5*t5interval*randNum
257 G4double t5 = -absA5/absB5 + G4Exp(argexp)
258 G4double sint5 = std::sin(t5);
259 G4double cost5 = std::cos(t5);
260
261 dirElectron.set(sint5*cosp5, sint5*sinp5,
262 G4ThreeVector eMomentumVector = eMomentum*
263
264 G4ThreeVector auxVec2 = muMomentumVector -
265 G4double p1mp3mp52 = auxVec2.dot(auxVec2);
266 G4double Bp = muFinalMomentum*(sint3*cosp3
267 eMomentum*(sint5*cosp5*cosp6 + sint5*sin
268 G4double Cp = -muMomentum + muFMomentum*co
269 G4double A6 = p1mp3mp52 + pMomentum*pMomen
270 G4double B6 = 2.*pMomentum*Bp;
271 G4double C6 = 2.*pMomentum*Cp;
272 G4double mint6 = 0.;
273 G4double maxt6 = CLHEP::pi;
274 G4double absA6C6 = std::abs(A6 + C6);
275 G4double absB6 = std::abs(B6);
276 G4double t6interval = (1./(absA6C6 + absB6
277 G4double t6 = (-absA6C6 + 1./(1./(absA6C6
278 G4double sint6 = std::sin(t6);
279 G4double cost6 = std::cos(t6);
280
281 dirPositron.set(sint6*cosp6, sint6*sinp6,
282
283 PhiRotation(dirElectron, phi3);
284 PhiRotation(dirPositron, phi3);
285
286 } else if (randNumbs[7] >= W[2]) {
287 G4double phi3 = CLHEP::twopi*randNumbs[0];
288 G4double phi6 = CLHEP::twopi*randNumbs[1];
289 G4double phi5 = CLHEP::twopi*randNumbs[2];
290 G4double minmuFinalEnergy = muMass;
291 G4double muFinalEnergyinterval = muEnergy
292 G4double muFEnergy = minmuFinalEnergy + mu
293
294 G4double minpEnergy = eMass;
295 G4double maxpEnergy = muEnergy - muFEnergy
296 G4double pEnergyinterval = maxpEnergy - mi
297 G4double pEnergy = minpEnergy + pEnergyint
298
299 G4double cosp3 = 1.;
300 G4double sinp3 = 0.;
301 G4double cosp5 = std::cos(phi5);
302 G4double sinp5 = std::sin(phi5);
303 G4double cosp6 = std::cos(phi6);
304 G4double sinp6 = std::sin(phi6);
305
306 G4double muFMomentum = std::sqrt(muFEnergy
307 G4double pMomentum = std::sqrt(pEnergy*pEn
308 G4double eEnergy = muEnergy - muFEnergy -
309 G4double eMomentum = std::sqrt(eEnergy*eEn
310
311 G4double A3 = -2.*muMass*muMass + 2.*muEne
312 G4double B3 = -2.*muMomentum*muFMomentum;
313 G4double cost3interval = G4Log((A3 + B3*Cm
314 G4double expanCost3r6 = G4Exp(B3*cost3inte
315 G4double cost3 = A3*(expanCost3r6 - 1.)/B3
316 G4double sint3 = std::sqrt((1. - cost3)*(1
317
318 G4ThreeVector muFinalMomentumVector;
319 muFinalMomentumVector.set(muFMomentum*sint
320 muFMomentum*cost3);
321
322 G4LorentzVector muFourMomentum(muMomentum,
323 muFinalFourMomentum.set(muFEnergy, muFinal
324 G4LorentzVector auxVec1 = muFourMomentum -
325 G4double A6 = auxVec1.mag2() - 2.*pEnergy*
326 2.*muMomentumVector[2]*pMomentum - 2.*mu
327 G4double B6 = -2.*muFMomentum*pMomentum*(s
328 G4double absA6 = std::abs(A6);
329 G4double absB6 = std::abs(B6);
330 G4double mint6 = 0.;
331 G4double maxt6 = CLHEP::pi;
332 G4double t6interval = G4Log((absA6 + absB6
333 G4double argexp = absB6*t6interval*randNum
334 G4double t6 = -absA6/absB6 + G4Exp(argexp)
335 G4double sint6 = std::sin(t6);
336 G4double cost6 = std::cos(t6);
337
338 dirPositron.set(sint6*cosp6, sint6*sinp6,
339 G4ThreeVector pMomentumVector = pMomentum*
340
341 G4ThreeVector auxVec2 = muMomentumVector -
342 G4double p1mp3mp62 = auxVec2.dot(auxVec2);
343 G4double Bp = muFMomentum*(sint3*cosp3*cos
344 pMomentum*(sint6*cosp6*cosp5 + sint6*sin
345 G4double Cp = -muMomentum + muFMomentum*co
346 G4double A5 = p1mp3mp62 + eMomentum*eMomen
347 G4double B5 = 2.*eMomentum*Bp;
348 G4double C5 = 2.*eMomentum*Cp;
349 G4double mint5 = 0.;
350 G4double maxt5 = CLHEP::pi;
351 G4double absA5C5 = std::abs(A5 + C5);
352 G4double absB5 = std::abs(B5);
353 G4double t5interval = (1./(absA5C5 + absB5
354 G4double t5 = (-absA5C5 + 1./(1./(absA5C5
355 G4double sint5 = std::sin(t5);
356 G4double cost5 = std::cos(t5);
357
358 dirElectron.set(sint5*cosp5, sint5*sinp5,
359
360 PhiRotation(dirElectron, phi3);
361 PhiRotation(dirPositron, phi3);
362 }
363 return muFinalFourMomentum;
364 }
365
366 //....oooOO0OOooo........oooOO0OOooo........oo
367 //....oooOO0OOooo........oooOO0OOooo........oo
368
369 void G4RiGeAngularGenerator::PhiRotation(G4Thr
370 {
371 G4double sinp = std::sin(phi);
372 G4double cosp = std::cos(phi);
373
374 G4double newX = dir.x()*cosp + dir.y()*sinp;
375 G4double newY = -dir.x()*sinp + dir.y()*cosp
376
377 dir.setX(newX);
378 dir.setY(newY);
379 }
380
381 //....oooOO0OOooo........oooOO0OOooo........oo
382
383 G4LorentzVector G4RiGeAngularGenerator::eDP2(G
384 G
385 G
386 {
387 G4double sint = std::sqrt((1.0 - x4)*(1.0 +
388 G4double cosp = std::cos(x5);
389 G4double sinp = std::sin(x5);
390
391 G4double QJM2 = (x1 + x3 - x2)*(x1 + x3 - x2
392
393 if (QJM2 < 0.) {
394 QJM2 = 1.e-13;
395 }
396
397 G4double QJM = std::sqrt(QJM2);
398 G4LorentzVector x6(std::sqrt(x2 + QJM2), QJM
399
400 return x6;
401 }
402
403 //....oooOO0OOooo........oooOO0OOooo........oo
404
405 G4LorentzVector G4RiGeAngularGenerator::pDP2(G
406 {
407 G4LorentzVector x7(x3 + x6.vect().dot(x6.vec
408 return x7;
409 }
410
411 //....oooOO0OOooo........oooOO0OOooo........oo
412
413 void G4RiGeAngularGenerator::PrintGeneratorInf
414 {
415 G4cout << "\n" << G4endl;
416 G4cout << "Angular Generator by RiGe algorit
417 }
418
419 //....oooOO0OOooo........oooOO0OOooo........oo
420