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