| 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 // G4MuonRadiativeDecayChannelWithSpin class implementation
27 //
28 // References:
29 // - TRIUMF/TWIST Technote TN-55:
30 // "Radiative muon decay" by P. Depommier and A. Vacheret
31 // - Yoshitaka Kuno and Yasuhiro Okada
32 // "Muon Decays and Physics Beyond the Standard Model"
33 // Rev. Mod. Phys. 73, 151 (2001)
34 //
35 // Author: P.Gumplinger - Triumf, 25 July 2007
36 // Revision: D.Mingming - Center for HEP, Tsinghua Univ., 10 August 2011
37 // --------------------------------------------------------------------
38
39 #include "G4MuonRadiativeDecayChannelWithSpin.hh"
40
41 #include "G4DecayProducts.hh"
42 #include "G4LorentzVector.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "Randomize.hh"
46
47 G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin(
48 const G4String& theParentName, G4double theBR)
49 : G4VDecayChannel("Radiative Muon Decay", 1)
50 {
51 // set names for daughter particles
52 if (theParentName == "mu+") {
53 SetBR(theBR);
54 SetParent("mu+");
55 SetNumberOfDaughters(4);
56 SetDaughter(0, "e+");
57 SetDaughter(1, "gamma");
58 SetDaughter(2, "nu_e");
59 SetDaughter(3, "anti_nu_mu");
60 }
61 else if (theParentName == "mu-") {
62 SetBR(theBR);
63 SetParent("mu-");
64 SetNumberOfDaughters(4);
65 SetDaughter(0, "e-");
66 SetDaughter(1, "gamma");
67 SetDaughter(2, "anti_nu_e");
68 SetDaughter(3, "nu_mu");
69 }
70 else {
71 #ifdef G4VERBOSE
72 if (GetVerboseLevel() > 0) {
73 G4cout << "G4RadiativeMuonDecayChannel::G4RadiativeMuonDecayChannel():";
74 G4cout << " parent particle is not muon but ";
75 G4cout << theParentName << G4endl;
76 }
77 #endif
78 }
79 }
80
81 G4MuonRadiativeDecayChannelWithSpin&
82 G4MuonRadiativeDecayChannelWithSpin::operator=(const G4MuonRadiativeDecayChannelWithSpin& right)
83 {
84 if (this != &right) {
85 kinematics_name = right.kinematics_name;
86 verboseLevel = right.verboseLevel;
87 rbranch = right.rbranch;
88
89 // copy parent name
90 parent_name = new G4String(*right.parent_name);
91
92 // clear daughters_name array
93 ClearDaughtersName();
94
95 // recreate array
96 numberOfDaughters = right.numberOfDaughters;
97 if (numberOfDaughters > 0) {
98 if (daughters_name != nullptr) ClearDaughtersName();
99 daughters_name = new G4String*[numberOfDaughters];
100 // copy daughters name
101 for (G4int index = 0; index < numberOfDaughters; ++index) {
102 daughters_name[index] = new G4String(*right.daughters_name[index]);
103 }
104 }
105 parent_polarization = right.parent_polarization;
106 }
107 return *this;
108 }
109
110 G4DecayProducts* G4MuonRadiativeDecayChannelWithSpin::DecayIt(G4double)
111 {
112 #ifdef G4VERBOSE
113 if (GetVerboseLevel() > 1) G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt()";
114 #endif
115
116 CheckAndFillParent();
117 CheckAndFillDaughters();
118
119 // parent mass
120 G4double parentmass = G4MT_parent->GetPDGMass();
121
122 G4double EMMU = parentmass;
123
124 // daughters'mass
125 G4double daughtermass[4];
126 // G4double sumofdaughtermass = 0.0;
127 for (G4int index = 0; index < 4; ++index) {
128 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
129 // sumofdaughtermass += daughtermass[index];
130 }
131
132 G4double EMASS = daughtermass[0];
133
134 // create parent G4DynamicParticle at rest
135 G4ThreeVector dummy;
136 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
137
138 // create G4Decayproducts
139 auto products = new G4DecayProducts(*parentparticle);
140 delete parentparticle;
141
142 G4double eps = EMASS / EMMU;
143
144 G4double som0, x, y, xx, yy, zz;
145 G4double cthetaE, cthetaG, cthetaGE, phiE, phiG;
146 const std::size_t MAX_LOOP = 10000;
147
148 for (std::size_t loop_counter1 = 0; loop_counter1 < MAX_LOOP; ++loop_counter1) {
149 for (std::size_t loop_counter2 = 0; loop_counter2 < MAX_LOOP; ++loop_counter2) {
150 // -------------------------------------------------------------------
151 // Build two vectors of random length and random direction, for the
152 // positron and the photon.
153 // x/y is the length of the vector, xx, yy and zz the components,
154 // phi is the azimutal angle, theta the polar angle.
155 // -------------------------------------------------------------------
156
157 // For the positron
158 //
159 x = G4UniformRand();
160
161 rn3dim(xx, yy, zz, x);
162
163 if (std::fabs((xx * xx) + (yy * yy) + (zz * zz) - (x * x)) > 0.001) {
164 G4cout << "Norm of x not correct" << G4endl;
165 }
166
167 phiE = atan4(xx, yy);
168 cthetaE = zz / x;
169 G4double sthetaE = std::sqrt((xx * xx) + (yy * yy)) / x;
170
171 // What you get:
172 //
173 // x = positron energy
174 // phiE = azimutal angle of positron momentum
175 // cthetaE = cosine of polar angle of positron momentum
176 // sthetaE = sine of polar angle of positron momentum
177 //
178 //// G4cout << " x, xx, yy, zz " << x << " " << xx << " "
179 //// << yy << " " << zz << G4endl;
180 //// G4cout << " phiE, cthetaE, sthetaE " << phiE << " "
181 //// << cthetaE << " "
182 //// << sthetaE << " " << G4endl;
183
184 // For the photon
185 //
186 y = G4UniformRand();
187
188 rn3dim(xx, yy, zz, y);
189
190 if (std::fabs((xx * xx) + (yy * yy) + (zz * zz) - (y * y)) > 0.001) {
191 G4cout << " Norm of y not correct " << G4endl;
192 }
193
194 phiG = atan4(xx, yy);
195 cthetaG = zz / y;
196 G4double sthetaG = std::sqrt((xx * xx) + (yy * yy)) / y;
197
198 // What you get:
199 //
200 // y = photon energy
201 // phiG = azimutal angle of photon momentum
202 // cthetaG = cosine of polar angle of photon momentum
203 // sthetaG = sine of polar angle of photon momentum
204 //
205 //// G4cout << " y, xx, yy, zz " << y << " " << xx << " "
206 //// << yy << " " << zz << G4endl;
207 //// G4cout << " phiG, cthetaG, sthetaG " << phiG << " "
208 //// << cthetaG << " "
209 //// << sthetaG << " " << G4endl;
210
211 // Calculate the angle between positron and photon (cosine)
212 //
213 cthetaGE = cthetaE * cthetaG + sthetaE * sthetaG * std::cos(phiE - phiG);
214
215 //// G4cout << x << " " << cthetaE << " " << sthetaE << " "
216 //// << y << " " << cthetaG << " " << sthetaG << " "
217 //// << cthetaGE
218
219 G4double term0 = eps * eps;
220 G4double term1 = x * ((1.0 - eps) * (1.0 - eps)) + 2.0 * eps;
221 G4double beta =
222 std::sqrt(x * ((1.0 - eps) * (1.0 - eps)) * (x * ((1.0 - eps) * (1.0 - eps)) + 4.0 * eps))
223 / term1;
224 G4double delta = 1.0 - beta * cthetaGE;
225
226 G4double term3 = y * (1.0 - (eps * eps));
227 G4double term6 = term1 * delta * term3;
228
229 G4double Qsqr = (1.0 - term1 - term3 + term0 + 0.5 * term6) / ((1.0 - eps) * (1.0 - eps));
230
231 // Check the kinematics.
232 //
233 if (Qsqr >= 0.0 && Qsqr <= 1.0) break;
234
235 } // end loop count
236
237 // Do the calculation for -1 muon polarization (i.e. mu+)
238 //
239 G4double Pmu = -1.0;
240 if (GetParentName() == "mu-") {
241 Pmu = +1.0;
242 }
243
244 som0 = fron(Pmu, x, y, cthetaE, cthetaG, cthetaGE);
245
246 // Sample the decay rate
247 //
248 if (G4UniformRand() * 177.0 <= som0) break;
249 }
250
251 G4double E = EMMU / 2. * (x * ((1. - eps) * (1. - eps)) + 2. * eps);
252 G4double G = EMMU / 2. * y * (1. - eps * eps);
253
254 if (E < EMASS) E = EMASS;
255
256 // calculate daughter momentum
257 G4double daughtermomentum[4];
258
259 daughtermomentum[0] = std::sqrt(E * E - EMASS * EMASS);
260
261 G4double sthetaE = std::sqrt(1. - cthetaE * cthetaE);
262 G4double cphiE = std::cos(phiE);
263 G4double sphiE = std::sin(phiE);
264
265 // Coordinates of the decay positron with respect to the muon spin
266
267 G4double px = sthetaE * cphiE;
268 G4double py = sthetaE * sphiE;
269 G4double pz = cthetaE;
270
271 G4ThreeVector direction0(px, py, pz);
272
273 direction0.rotateUz(parent_polarization);
274
275 auto daughterparticle0 =
276 new G4DynamicParticle(G4MT_daughters[0], daughtermomentum[0] * direction0);
277
278 products->PushProducts(daughterparticle0);
279
280 daughtermomentum[1] = G;
281
282 G4double sthetaG = std::sqrt(1. - cthetaG * cthetaG);
283 G4double cphiG = std::cos(phiG);
284 G4double sphiG = std::sin(phiG);
285
286 // Coordinates of the decay gamma with respect to the muon spin
287
288 px = sthetaG * cphiG;
289 py = sthetaG * sphiG;
290 pz = cthetaG;
291
292 G4ThreeVector direction1(px, py, pz);
293
294 direction1.rotateUz(parent_polarization);
295
296 auto daughterparticle1 =
297 new G4DynamicParticle(G4MT_daughters[1], daughtermomentum[1] * direction1);
298
299 products->PushProducts(daughterparticle1);
300
301 // daughter 3 ,4 (neutrinos)
302 // create neutrinos in the C.M frame of two neutrinos
303
304 G4double energy2 = parentmass - E - G;
305
306 G4ThreeVector P34 = -1. * (daughtermomentum[0] * direction0 + daughtermomentum[1] * direction1);
307 G4double vmass2 = energy2 * energy2 - P34.mag2();
308 G4double vmass = std::sqrt(vmass2);
309
310 G4double costhetan = 2. * G4UniformRand() - 1.0;
311 G4double sinthetan = std::sqrt((1.0 - costhetan) * (1.0 + costhetan));
312 G4double phin = twopi * G4UniformRand() * rad;
313 G4double sinphin = std::sin(phin);
314 G4double cosphin = std::cos(phin);
315
316 G4ThreeVector direction2(sinthetan * cosphin, sinthetan * sinphin, costhetan);
317
318 auto daughterparticle2 = new G4DynamicParticle(G4MT_daughters[2], direction2 * (vmass / 2.));
319 auto daughterparticle3 =
320 new G4DynamicParticle(G4MT_daughters[3], direction2 * (-1.0 * vmass / 2.));
321
322 // boost to the muon rest frame
323 G4double beta = P34.mag() / energy2;
324 G4ThreeVector direction34 = P34.unit();
325
326 G4LorentzVector p4 = daughterparticle2->Get4Momentum();
327 p4.boost(direction34.x() * beta, direction34.y() * beta, direction34.z() * beta);
328 daughterparticle2->Set4Momentum(p4);
329
330 p4 = daughterparticle3->Get4Momentum();
331 p4.boost(direction34.x() * beta, direction34.y() * beta, direction34.z() * beta);
332 daughterparticle3->Set4Momentum(p4);
333
334 products->PushProducts(daughterparticle2);
335 products->PushProducts(daughterparticle3);
336
337 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
338 daughtermomentum[3] = daughterparticle3->GetTotalMomentum();
339
340 // output message
341 #ifdef G4VERBOSE
342 if (GetVerboseLevel() > 1) {
343 G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt() -";
344 G4cout << " create decay products in rest frame " << G4endl;
345 G4double TT = daughterparticle0->GetTotalEnergy() + daughterparticle1->GetTotalEnergy()
346 + daughterparticle2->GetTotalEnergy() + daughterparticle3->GetTotalEnergy();
347 G4cout << "e :" << daughterparticle0->GetTotalEnergy() / MeV << G4endl;
348 G4cout << "gamma:" << daughterparticle1->GetTotalEnergy() / MeV << G4endl;
349 G4cout << "nu2 :" << daughterparticle2->GetTotalEnergy() / MeV << G4endl;
350 G4cout << "nu2 :" << daughterparticle3->GetTotalEnergy() / MeV << G4endl;
351 G4cout << "total:" << (TT - parentmass) / keV << G4endl;
352 if (GetVerboseLevel() > 1) {
353 products->DumpInfo();
354 }
355 }
356 #endif
357
358 return products;
359 }
360
361 G4double G4MuonRadiativeDecayChannelWithSpin::fron(G4double Pmu, G4double x, G4double y,
362 G4double cthetaE, G4double cthetaG,
363 G4double cthetaGE)
364 {
365 G4double mu = 105.65;
366 G4double me = 0.511;
367 G4double rho = 0.75;
368 G4double del = 0.75;
369 G4double eps = 0.0;
370 G4double kap = 0.0;
371 G4double ksi = 1.0;
372
373 G4double delta = 1 - cthetaGE;
374
375 // Calculation of the functions f(x,y)
376
377 G4double f_1s = 12.0
378 * ((y * y) * (1.0 - y) + x * y * (2.0 - 3.0 * y) + 2.0 * (x * x) * (1.0 - 2.0 * y)
379 - 2.0 * (x * x * x));
380 G4double f0s = 6.0
381 * (-x * y * (2.0 - 3.0 * (y * y)) - 2.0 * (x * x) * (1.0 - y - 3.0 * (y * y))
382 + 2.0 * (x * x * x) * (1.0 + 2.0 * y));
383 G4double f1s =
384 3.0 * ((x * x) * y * (2.0 - 3.0 * y - 3.0 * (y * y)) - (x * x * x) * y * (4.0 + 3.0 * y));
385 G4double f2s = 1.5 * ((x * x * x) * (y * y) * (2.0 + y));
386
387 G4double f_1se = 12.0 * (x * y * (1.0 - y) + (x * x) * (2.0 - 3.0 * y) - 2.0 * (x * x * x));
388 G4double f0se = 6.0 * (-(x * x) * (2.0 - y - 2.0 * (y * y)) + (x * x * x) * (2.0 + 3.0 * y));
389 G4double f1se = -3.0 * (x * x * x) * y * (2.0 + y);
390 G4double f2se = 0.0;
391
392 G4double f_1sg = 12.0 * ((y * y) * (1.0 - y) + x * y * (1.0 - 2.0 * y) - (x * x) * y);
393 G4double f0sg =
394 6.0 * (-x * (y * y) * (2.0 - 3.0 * y) - (x * x) * y * (1.0 - 4.0 * y) + (x * x * x) * y);
395 G4double f1sg = 3.0 * ((x * x) * (y * y) * (1.0 - 3.0 * y) - 2.0 * (x * x * x) * (y * y));
396 G4double f2sg = 1.5 * (x * x * x) * (y * y * y);
397
398 G4double f_1v = 8.0
399 * ((y * y) * (3.0 - 2.0 * y) + 6.0 * x * y * (1.0 - y)
400 + 2.0 * (x * x) * (3.0 - 4.0 * y) - 4.0 * (x * x * x));
401 G4double f0v = 8.0
402 * (-x * y * (3.0 - y - (y * y)) - (x * x) * (3.0 - y - 4.0 * (y * y))
403 + 2.0 * (x * x * x) * (1.0 + 2.0 * y));
404 G4double f1v =
405 2.0 * ((x * x) * y * (6.0 - 5.0 * y - 2.0 * (y * y)) - 2.0 * (x * x * x) * y * (4.0 + 3.0 * y));
406 G4double f2v = 2.0 * (x * x * x) * (y * y) * (2.0 + y);
407
408 G4double f_1ve =
409 8.0 * (x * y * (1.0 - 2.0 * y) + 2.0 * (x * x) * (1.0 - 3.0 * y) - 4.0 * (x * x * x));
410 G4double f0ve =
411 4.0 * (-(x * x) * (2.0 - 3.0 * y - 4.0 * (y * y)) + 2.0 * (x * x * x) * (2.0 + 3.0 * y));
412 G4double f1ve = -4.0 * (x * x * x) * y * (2.0 + y);
413 G4double f2ve = 0.0;
414
415 G4double f_1vg = 8.0 * ((y * y) * (1.0 - 2.0 * y) + x * y * (1.0 - 4.0 * y) - 2.0 * (x * x) * y);
416 G4double f0vg =
417 4.0 * (2.0 * x * (y * y) * (1.0 + y) - (x * x) * y * (1.0 - 4.0 * y) + 2.0 * (x * x * x) * y);
418 G4double f1vg = 2.0 * ((x * x) * (y * y) * (1.0 - 2.0 * y) - 4.0 * (x * x * x) * (y * y));
419 G4double f2vg = 2.0 * (x * x * x) * (y * y * y);
420
421 G4double f_1t = 8.0
422 * ((y * y) * (3.0 - y) + 3.0 * x * y * (2.0 - y) + 2.0 * (x * x) * (3.0 - 2.0 * y)
423 - 2.0 * (x * x * x));
424 G4double f0t = 4.0
425 * (-x * y * (6.0 + (y * y)) - 2.0 * (x * x) * (3.0 + y - 3.0 * (y * y))
426 + 2.0 * (x * x * x) * (1.0 + 2.0 * y));
427 G4double f1t =
428 2.0 * ((x * x) * y * (6.0 - 5.0 * y + (y * y)) - (x * x * x) * y * (4.0 + 3.0 * y));
429 G4double f2t = (x * x * x) * (y * y) * (2.0 + y);
430
431 G4double f_1te = -8.0 * (x * y * (1.0 + 3.0 * y) + (x * x) * (2.0 + 3.0 * y) + 2.0 * (x * x * x));
432 G4double f0te = 4.0 * ((x * x) * (2.0 + 3.0 * y + 4.0 * (y * y)) + (x * x * x) * (2.0 + 3.0 * y));
433 G4double f1te = -2.0 * (x * x * x) * y * (2.0 + y);
434 G4double f2te = 0.0;
435
436 G4double f_1tg = -8.0 * ((y * y) * (1.0 + y) + x * y + (x * x) * y);
437 G4double f0tg = 4.0 * (x * (y * y) * (2.0 - y) + (x * x) * y * (1.0 + 2.0 * y) + (x * x * x) * y);
438 G4double f1tg = -2.0 * ((x * x) * (y * y) * (1.0 - y) + 2.0 * (x * x * x) * y);
439 G4double f2tg = (x * x * x) * (y * y * y);
440
441 G4double term = delta + 2.0 * (me * me) / ((mu * mu) * (x * x));
442 term = 1.0 / term;
443
444 G4double nss = term * f_1s + f0s + delta * f1s + (delta * delta) * f2s;
445 G4double nv = term * f_1v + f0v + delta * f1v + (delta * delta) * f2v;
446 G4double nt = term * f_1t + f0t + delta * f1t + (delta * delta) * f2t;
447
448 G4double nse = term * f_1se + f0se + delta * f1se + (delta * delta) * f2se;
449 G4double nve = term * f_1ve + f0ve + delta * f1ve + (delta * delta) * f2ve;
450 G4double nte = term * f_1te + f0te + delta * f1te + (delta * delta) * f2te;
451
452 G4double nsg = term * f_1sg + f0sg + delta * f1sg + (delta * delta) * f2sg;
453 G4double nvg = term * f_1vg + f0vg + delta * f1vg + (delta * delta) * f2vg;
454 G4double ntg = term * f_1tg + f0tg + delta * f1tg + (delta * delta) * f2tg;
455
456 G4double term1 = nv;
457 G4double term2 = 2.0 * nss + nv - nt;
458 G4double term3 = 2.0 * nss - 2.0 * nv + nt;
459
460 G4double term1e = 1.0 / 3.0 * (1.0 - 4.0 / 3.0 * del);
461 G4double term2e = 2.0 * nse + 5.0 * nve - nte;
462 G4double term3e = 2.0 * nse - 2.0 * nve + nte;
463
464 G4double term1g = 1.0 / 3.0 * (1.0 - 4.0 / 3.0 * del);
465 G4double term2g = 2.0 * nsg + 5.0 * nvg - ntg;
466 G4double term3g = 2.0 * nsg - 2.0 * nvg + ntg;
467
468 G4double som00 = term1 + (1.0 - 4.0 / 3.0 * rho) * term2 + eps * term3;
469 G4double som01 = Pmu * ksi
470 * (cthetaE * (nve - term1e * term2e + kap * term3e)
471 + cthetaG * (nvg - term1g * term2g + kap * term3g));
472
473 G4double som0 = (som00 + som01) / y;
474 som0 = fine_structure_const / 8. / (twopi * twopi * twopi) * som0;
475
476 return som0;
477 }
478