| 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: G4RiGeAngularGenerator
33 //
34 // Authors: Girardo Depaola & Ricardo Pacheco
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........oooOO0OOooo........oooOO0OOooo....
48
49 G4RiGeAngularGenerator::G4RiGeAngularGenerator()
50 : G4VEmAngularDistribution("RiGeAngularGen")
51 {}
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55 G4ThreeVector&
56 G4RiGeAngularGenerator::SampleDirection(const G4DynamicParticle* dp,
57 G4double gEnergy, G4int,
58 const G4Material*)
59 {
60 // Sample gamma angle (Z - axis along the parent particle).
61 G4double cost = SampleCosTheta(dp->GetKineticEnergy(), gEnergy,
62 dp->GetDefinition()->GetPDGMass());
63 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
64 G4double phi = CLHEP::twopi*G4UniformRand();
65
66 fLocalDirection.set(sint*std::cos(phi), sint*std::sin(phi), cost);
67 fLocalDirection.rotateUz(dp->GetMomentumDirection());
68
69 return fLocalDirection;
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74 G4double G4RiGeAngularGenerator::SampleCosTheta(G4double primKinEnergy,
75 G4double gEnergy,
76 G4double mass)
77 {
78 G4double gam = 1.0 + primKinEnergy/mass;
79 G4double rmax = gam*CLHEP::halfpi*std::min(1.0, gam*mass/gEnergy - 1.0);
80 G4double rmax2= rmax*rmax;
81 G4double x = G4UniformRand()*rmax2/(1.0 + rmax2);
82
83 return std::cos(std::sqrt(x/(1.0 - x))/gam);
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87
88 G4LorentzVector
89 G4RiGeAngularGenerator::Sample5DPairDirections(const G4DynamicParticle* dp,
90 G4ThreeVector& dirElectron,
91 G4ThreeVector& dirPositron,
92 const G4double gEnergy, const G4double q2,
93 const G4double gMomentum,
94 const G4double muFinalMomentum,
95 const G4double muFinalEnergy,
96 const G4double* randNumbs,
97 const G4double* W)
98 {
99 G4double muEnergy = dp->GetKineticEnergy();
100 G4ThreeVector muMomentumVector = dp->GetMomentum();
101 G4double muMomentum = muMomentumVector.mag();
102 G4LorentzVector muFinalFourMomentum(muEnergy, muMomentumVector);
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()->GetPDGMass();
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 - B1))/B1;
120
121 G4double costg = (-A1 + (A1 - B1)*G4Exp(B1*tginterval*randNumbs[1]))/B1;
122 G4double sintg = std::sqrt((1.0 - costg)*(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, costg);
129 G4LorentzVector gFourMomentum(gEnergy, dirGamma*gMomentum);
130
131 G4double Ap = muMomentum*muMomentum + muFinalMomentum*muFinalMomentum + gMomentum*gMomentum;
132 G4double A = Ap - 2.*muMomentum*gMomentum*costg;
133 G4double B = 2.*muFinalMomentum*gMomentum*sintg*cospg;
134 G4double C = 2.*muFinalMomentum*gMomentum*costg - 2.*muMomentum*muFinalMomentum;
135 G4double absB = std::abs(B);
136 G4double t1interval = (1./(A + C + absB*mint3) - 1./(A + C + absB*maxt3))/absB;
137 G4double t1 = (-(A + C) + 1./(1./(A + C + absB*mint3) - absB*t1interval*randNumbs[0]))/absB;
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, eMass2, eMass2, cost5, phi5);
148 G4LorentzVector pFourMomentumMQ = pDP2(eMass2, eFourMomentumMQ);
149
150 G4LorentzVector eFourMomentum = eFourMomentumMQ.boost(gFourMomentum.boostVector());
151 G4LorentzVector pFourMomentum = pFourMomentumMQ.boost(gFourMomentum.boostVector());
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[7] < W[1]) {
161 G4double A3 = q2 + 2.*gEnergy*muFinalEnergy;
162 G4double B3 = -2.*gMomentum*muFinalMomentum;
163
164 G4double tQ3interval = G4Log((A3 + B3)/(A3 - B3))/B3;
165 G4double tQMG = (-A3 + (A3 - B3)*G4Exp(B3*tQ3interval*randNumbs[0]))/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 + muFinalMomentum*muFinalMomentum + gMomentum*gMomentum;
173 G4double A = Ap + 2.*muFinalMomentum*gMomentum*tQMG;
174 G4double B = -2.*muMomentum*gMomentum*sintQ3*cospQP;
175 G4double C = -2.*muMomentum*gMomentum*tQMG - 2.*muMomentum*muFinalMomentum;
176
177 G4double absB = std::abs(B);
178 G4double t3interval = (1./(A + C + absB*mint3) - 1./(A + C + absB*maxt3))/absB;
179 G4double t3 = (-(A + C) + 1./(1./(A + C + absB*mint3) - absB*t3interval*randNumbs[0]))/absB;
180 G4double sint3 = std::sin(t3);
181 G4double cost3 = std::cos(t3);
182
183 G4double cost = -sint3*sintQ3*cospQP + cost3*tQMG;
184 G4double sint = std::sqrt((1. + cost)*(1. - cost));
185 G4double cosp = (sintQ3*cospQP*cost3 + sint3*tQMG)/sint;
186 G4double sinp = sintQ3*sinpQP/sint;
187
188 G4ThreeVector dirGamma;
189 dirGamma.set(sint*cosp, sint*sinp, cost);
190 G4LorentzVector gFourMomentum(gEnergy, dirGamma*gMomentum);
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, eMass2, eMass2, cost5, phi5);
199 G4LorentzVector pFourMomentumMQ = pDP2(eMass2, eFourMomentumMQ);
200
201 G4LorentzVector eFourMomentum = eFourMomentumMQ.boost(gFourMomentum.boostVector());
202 G4LorentzVector pFourMomentum = pFourMomentumMQ.boost(gFourMomentum.boostVector());
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[7] < W[2]) {
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.*eMass - minmuFinalEnergy;
217 G4double muFEnergy = minmuFinalEnergy + muEnergyInterval*randNumbs[3];
218
219 G4double mineEnergy = eMass;
220 G4double maxeEnergy = muEnergy - muFEnergy - eMass;
221 G4double eEnergyInterval = maxeEnergy - mineEnergy;
222 G4double eEnergy = mineEnergy + eEnergyInterval*randNumbs[4];
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*muFEnergy - muMass*muMass);
232 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
233 G4double pEnergy = muEnergy - muFEnergy - eEnergy;
234 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
235
236 G4double A3 = -2.*muMass*muMass + 2.*muEnergy*muFinalEnergy;
237 G4double B3 = -2.*muMomentum*muFinalMomentum;
238 G4double cost3interval = G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
239 G4double expanCost3r6 = G4Exp(B3*cost3interval*randNumbs[5]);
240 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
241 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
242
243 G4ThreeVector muFinalMomentumVector(muFMomentum*sint3, 0., muFMomentum*cost3);
244
245 G4LorentzVector muFourMomentum(muMomentum, muMomentumVector);
246 muFinalFourMomentum.set(muFEnergy, muFinalMomentumVector);
247 G4LorentzVector auxVec1 = muFourMomentum - muFinalFourMomentum;
248 G4double A5 = auxVec1.mag2() - 2.*eEnergy*(muEnergy - muFEnergy) +
249 2.*muMomentumVector[2]*eMomentum - 2.*muFMomentum*eMomentum*cost3;
250 G4double B5 = -2.*muFMomentum*eMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5);
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*maxt5)/(absA5 + absB5*mint5))/absB5;
256 G4double argexp = absB5*t5interval*randNumbs[6] + G4Log(absA5 + absB5*mint5);
257 G4double t5 = -absA5/absB5 + G4Exp(argexp)/absB5;
258 G4double sint5 = std::sin(t5);
259 G4double cost5 = std::cos(t5);
260
261 dirElectron.set(sint5*cosp5, sint5*sinp5, cost5);
262 G4ThreeVector eMomentumVector = eMomentum*dirElectron;
263
264 G4ThreeVector auxVec2 = muMomentumVector - muFinalMomentumVector - eMomentumVector;
265 G4double p1mp3mp52 = auxVec2.dot(auxVec2);
266 G4double Bp = muFinalMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6) +
267 eMomentum*(sint5*cosp5*cosp6 + sint5*sinp5*sinp6);
268 G4double Cp = -muMomentum + muFMomentum*cost3 + eMomentum*cost5;
269 G4double A6 = p1mp3mp52 + pMomentum*pMomentum;
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*mint6) - 1./(absA6C6 + absB6*maxt6))/absB6;
277 G4double t6 = (-absA6C6 + 1./(1./(absA6C6 + absB6*mint6) - absB6*t6interval*randNumbs[8]))/absB6;
278 G4double sint6 = std::sin(t6);
279 G4double cost6 = std::cos(t6);
280
281 dirPositron.set(sint6*cosp6, sint6*sinp6, cost6);
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 - 2.*eMass - minmuFinalEnergy;
292 G4double muFEnergy = minmuFinalEnergy + muFinalEnergyinterval*randNumbs[3];
293
294 G4double minpEnergy = eMass;
295 G4double maxpEnergy = muEnergy - muFEnergy - eMass;
296 G4double pEnergyinterval = maxpEnergy - minpEnergy;
297 G4double pEnergy = minpEnergy + pEnergyinterval*randNumbs[4];
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*muFEnergy - muMass*muMass);
307 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
308 G4double eEnergy = muEnergy - muFEnergy - pEnergy;
309 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
310
311 G4double A3 = -2.*muMass*muMass + 2.*muEnergy*muFinalEnergy;
312 G4double B3 = -2.*muMomentum*muFMomentum;
313 G4double cost3interval = G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
314 G4double expanCost3r6 = G4Exp(B3*cost3interval*randNumbs[5]);
315 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
316 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
317
318 G4ThreeVector muFinalMomentumVector;
319 muFinalMomentumVector.set(muFMomentum*sint3*cosp3, muFMomentum*sint3*sinp3,
320 muFMomentum*cost3);
321
322 G4LorentzVector muFourMomentum(muMomentum, muMomentumVector);
323 muFinalFourMomentum.set(muFEnergy, muFinalMomentumVector);
324 G4LorentzVector auxVec1 = muFourMomentum - muFinalFourMomentum;
325 G4double A6 = auxVec1.mag2() - 2.*pEnergy*(muEnergy - muFEnergy) +
326 2.*muMomentumVector[2]*pMomentum - 2.*muFMomentum*pMomentum*cost3;
327 G4double B6 = -2.*muFMomentum*pMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6);
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*maxt6)/(absA6 + absB6*mint6))/absB6;
333 G4double argexp = absB6*t6interval*randNumbs[6] + G4Log(absA6 + absB6*mint6);
334 G4double t6 = -absA6/absB6 + G4Exp(argexp)/absB6;
335 G4double sint6 = std::sin(t6);
336 G4double cost6 = std::cos(t6);
337
338 dirPositron.set(sint6*cosp6, sint6*sinp6, cost6);
339 G4ThreeVector pMomentumVector = pMomentum*dirPositron;
340
341 G4ThreeVector auxVec2 = muMomentumVector - muFinalMomentumVector - pMomentumVector;
342 G4double p1mp3mp62 = auxVec2.dot(auxVec2);
343 G4double Bp = muFMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5) +
344 pMomentum*(sint6*cosp6*cosp5 + sint6*sinp6*sinp5);
345 G4double Cp = -muMomentum + muFMomentum*cost3 + pMomentum*cost6;
346 G4double A5 = p1mp3mp62 + eMomentum*eMomentum;
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*mint5) - 1./(absA5C5 + absB5*maxt5))/absB5;
354 G4double t5 = (-absA5C5 + 1./(1./(absA5C5 + absB5*mint5) - absB5*t5interval*randNumbs[8]))/absB5;
355 G4double sint5 = std::sin(t5);
356 G4double cost5 = std::cos(t5);
357
358 dirElectron.set(sint5*cosp5, sint5*sinp5, cost5);
359
360 PhiRotation(dirElectron, phi3);
361 PhiRotation(dirPositron, phi3);
362 }
363 return muFinalFourMomentum;
364 }
365
366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
368
369 void G4RiGeAngularGenerator::PhiRotation(G4ThreeVector& dir, G4double phi)
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........oooOO0OOooo........oooOO0OOooo....
382
383 G4LorentzVector G4RiGeAngularGenerator::eDP2(G4double x1, G4double x2,
384 G4double x3, G4double x4,
385 G4double x5)
386 {
387 G4double sint = std::sqrt((1.0 - x4)*(1.0 + x4));
388 G4double cosp = std::cos(x5);
389 G4double sinp = std::sin(x5);
390
391 G4double QJM2 = (x1 + x3 - x2)*(x1 + x3 - x2)/(4.*x1) - x3;
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*sint*cosp, QJM*sint*sinp, QJM*x4);
399
400 return x6;
401 }
402
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
404
405 G4LorentzVector G4RiGeAngularGenerator::pDP2(G4double x3, const G4LorentzVector& x6)
406 {
407 G4LorentzVector x7(x3 + x6.vect().dot(x6.vect()), -x6.vect());
408 return x7;
409 }
410
411 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
412
413 void G4RiGeAngularGenerator::PrintGeneratorInformation() const
414 {
415 G4cout << "\n" << G4endl;
416 G4cout << "Angular Generator by RiGe algorithm" << G4endl;
417 }
418
419 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
420