Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuonRadiativeDecayChannelWithSpin.cc
Go to the documentation of this file.
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  parent_polarization()
56 {
57 }
58 
59 G4MuonRadiativeDecayChannelWithSpin::
60  G4MuonRadiativeDecayChannelWithSpin(const G4String& theParentName,
61  G4double theBR)
62  : G4VDecayChannel("Radiative Muon Decay",1),
63  parent_polarization()
64 {
65  // set names for daughter particles
66  if (theParentName == "mu+") {
67  SetBR(theBR);
68  SetParent("mu+");
69  SetNumberOfDaughters(4);
70  SetDaughter(0, "e+");
71  SetDaughter(1, "gamma");
72  SetDaughter(2, "nu_e");
73  SetDaughter(3, "anti_nu_mu");
74  } else if (theParentName == "mu-") {
75  SetBR(theBR);
76  SetParent("mu-");
77  SetNumberOfDaughters(4);
78  SetDaughter(0, "e-");
79  SetDaughter(1, "gamma");
80  SetDaughter(2, "anti_nu_e");
81  SetDaughter(3, "nu_mu");
82  } else {
83 #ifdef G4VERBOSE
84  if (GetVerboseLevel()>0) {
85  G4cout << "G4RadiativeMuonDecayChannel:: constructor :";
86  G4cout << " parent particle is not muon but ";
87  G4cout << theParentName << G4endl;
88  }
89 #endif
90  }
91 }
92 
93 G4MuonRadiativeDecayChannelWithSpin::~G4MuonRadiativeDecayChannelWithSpin()
94 {
95 }
96 
97 G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin(const G4MuonRadiativeDecayChannelWithSpin &right):
98  G4VDecayChannel(right)
99 {
100  parent_polarization = right.parent_polarization;
101 }
102 
103 G4MuonRadiativeDecayChannelWithSpin & G4MuonRadiativeDecayChannelWithSpin::operator=(const G4MuonRadiativeDecayChannelWithSpin & right)
104 {
105  if (this != &right) {
106  kinematics_name = right.kinematics_name;
107  verboseLevel = right.verboseLevel;
108  rbranch = right.rbranch;
109 
110  // copy parent name
111  parent_name = new G4String(*right.parent_name);
112 
113  // clear daughters_name array
114  ClearDaughtersName();
115 
116  // recreate array
117  numberOfDaughters = right.numberOfDaughters;
118  if ( numberOfDaughters >0 ) {
119  if (daughters_name !=0) ClearDaughtersName();
120  daughters_name = new G4String*[numberOfDaughters];
121  //copy daughters name
122  for (G4int index=0; index < numberOfDaughters; index++) {
123  daughters_name[index] = new G4String(*right.daughters_name[index]);
124  }
125  }
126  parent_polarization = right.parent_polarization;
127  }
128  return *this;
129 }
130 
131 
132 G4DecayProducts *G4MuonRadiativeDecayChannelWithSpin::DecayIt(G4double)
133 {
134 
135 #ifdef G4VERBOSE
136  if (GetVerboseLevel()>1)
137  G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
138 #endif
139 
140  if (parent == 0) FillParent();
141  if (daughters == 0) FillDaughters();
142 
143  // parent mass
144  G4double parentmass = parent->GetPDGMass();
145 
146  G4double EMMU = parentmass;
147 
148  //daughters'mass
149  G4double daughtermass[4];
150  G4double sumofdaughtermass = 0.0;
151  for (G4int index=0; index<4; index++){
152  daughtermass[index] = daughters[index]->GetPDGMass();
153  sumofdaughtermass += daughtermass[index];
154  }
155 
156  G4double EMASS = daughtermass[0];
157 
158  //create parent G4DynamicParticle at rest
159  G4ThreeVector dummy;
160  G4DynamicParticle * parentparticle =
161  new G4DynamicParticle( parent, dummy, 0.0);
162  //create G4Decayproducts
163  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
164  delete parentparticle;
165 
166  G4int i = 0;
167 
168  G4double eps = EMASS/EMMU;
169 
170  G4double som0, Qsqr, x, y, xx, yy, zz;
171  G4double cthetaE, cthetaG, cthetaGE, phiE, phiG;
172 
173  do {
174 
175 // leap1:
176 
177  i++;
178 
179 // leap2:
180 
181  do {
182 //
183 //--------------------------------------------------------------------------
184 // Build two vectors of random length and random direction, for the
185 // positron and the photon.
186 // x/y is the length of the vector, xx, yy and zz the components,
187 // phi is the azimutal angle, theta the polar angle.
188 //--------------------------------------------------------------------------
189 //
190 // For the positron
191 //
192  x = G4UniformRand();
193 
194  rn3dim(xx,yy,zz,x);
195 
196  if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(x*x))>0.001){
197  G4cout << "Norm of x not correct" << G4endl;
198  }
199 
200  phiE = atan4(xx,yy);
201  cthetaE = zz/x;
202  G4double sthetaE = std::sqrt((xx*xx)+(yy*yy))/x;
203 //
204 // What you get:
205 //
206 // x = positron energy
207 // phiE = azimutal angle of positron momentum
208 // cthetaE = cosine of polar angle of positron momentum
209 // sthetaE = sine of polar angle of positron momentum
210 //
216 //
217 //-----------------------------------------------------------------------
218 //
219 // For the photon
220 //
221  y = G4UniformRand();
222 
223  rn3dim(xx,yy,zz,y);
224 
225  if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(y*y))>0.001){
226  G4cout << " Norm of y not correct " << G4endl;
227  }
228 
229  phiG = atan4(xx,yy);
230  cthetaG = zz/y;
231  G4double sthetaG = std::sqrt((xx*xx)+(yy*yy))/y;
232 //
233 // What you get:
234 //
235 // y = photon energy
236 // phiG = azimutal angle of photon momentum
237 // cthetaG = cosine of polar angle of photon momentum
238 // sthetaG = sine of polar angle of photon momentum
239 //
245 //
246 //-----------------------------------------------------------------------
247 //
248 // Maybe certain restrictions on the kinematical variables:
249 //
254 //
255 //-----------------------------------------------------------------------
256 //
257 // Calculate the angle between positron and photon (cosine)
258 //
259  cthetaGE = cthetaE*cthetaG+sthetaE*sthetaG*std::cos(phiE-phiG);
260 //
264 //
265 //-----------------------------------------------------------------------
266 //
267  G4double term0 = eps*eps;
268  G4double term1 = x*((1.0-eps)*(1.0-eps))+2.0*eps;
269  G4double beta = std::sqrt( x*((1.0-eps)*(1.0-eps))*
270  (x*((1.0-eps)*(1.0-eps))+4.0*eps))/term1;
271  G4double delta = 1.0-beta*cthetaGE;
272 
273  G4double term3 = y*(1.0-(eps*eps));
274  G4double term6 = term1*delta*term3;
275 
276  Qsqr = (1.0-term1-term3+term0+0.5*term6)/((1.0-eps)*(1.0-eps));
277 //
278 //-----------------------------------------------------------------------
279 //
280 // Check the kinematics.
281 //
282  } while ( Qsqr<0.0 || Qsqr>1.0 );
283 //
285 //
286 // Do the calculation for -1 muon polarization (i.e. mu+)
287 //
288  G4double Pmu = -1.0;
289  if (GetParentName() == "mu-")Pmu = +1.0;
290 //
291 // and for Fronsdal
292 //
293 //-----------------------------------------------------------------------
294 //
295  som0 = fron(Pmu,x,y,cthetaE,cthetaG,cthetaGE);
296 //
304 //
305 //-----------------------------------------------------------------------
306 //
308 //
309 //----------------------------------------------------------------------
310 //
311 // Sample the decay rate
312 //
313 
314  } while (G4UniformRand()*177.0 > som0);
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( 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( 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*(1.0 - (x+y)/2.0);
377 
378  G4double vmass = std::sqrt((energy2-
379  (daughtermomentum[0]+daughtermomentum[1]))*
380  (energy2+
381  (daughtermomentum[0]+daughtermomentum[1])));
382  G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2;
383  beta = -1.0 * std::min(beta,0.99);
384 
385  G4double costhetan = 2.*G4UniformRand()-1.0;
386  G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
387  G4double phin = twopi*G4UniformRand()*rad;
388  G4double sinphin = std::sin(phin);
389  G4double cosphin = std::cos(phin);
390 
391  G4ThreeVector direction2(sinthetan*cosphin,sinthetan*sinphin,costhetan);
392 
393  G4DynamicParticle * daughterparticle2
394  = new G4DynamicParticle( daughters[2], direction2*(vmass/2.));
395  G4DynamicParticle * daughterparticle3
396  = new G4DynamicParticle( daughters[3], direction2*(-1.0*vmass/2.));
397 
398  // boost to the muon rest frame
399 
400  G4ThreeVector direction34(direction0.x()+direction1.x(),
401  direction0.y()+direction1.y(),
402  direction0.z()+direction1.z());
403  direction34 = direction34.unit();
404 
405  G4LorentzVector p4 = daughterparticle2->Get4Momentum();
406  p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
407  daughterparticle2->Set4Momentum(p4);
408 
409  p4 = daughterparticle3->Get4Momentum();
410  p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
411  daughterparticle3->Set4Momentum(p4);
412 
413  products->PushProducts(daughterparticle2);
414  products->PushProducts(daughterparticle3);
415 
416  daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
417  daughtermomentum[3] = daughterparticle3->GetTotalMomentum();
418 
419 // output message
420 #ifdef G4VERBOSE
421  if (GetVerboseLevel()>1) {
422  G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
423  G4cout << " create decay products in rest frame " <<G4endl;
424  products->DumpInfo();
425  }
426 #endif
427  return products;
428 }
429 
430 G4double G4MuonRadiativeDecayChannelWithSpin::fron(G4double Pmu,
431  G4double x,
432  G4double y,
433  G4double cthetaE,
434  G4double cthetaG,
435  G4double cthetaGE)
436 {
437  G4double mu = 105.65;
438  G4double me = 0.511;
439  G4double rho = 0.75;
440  G4double del = 0.75;
441  G4double eps = 0.0;
442  G4double kap = 0.0;
443  G4double ksi = 1.0;
444 
445  G4double delta = 1-cthetaGE;
446 
447 // Calculation of the functions f(x,y)
448 
449  G4double f_1s = 12.0*((y*y)*(1.0-y)+x*y*(2.0-3.0*y)
450  +2.0*(x*x)*(1.0-2.0*y)-2.0*(x*x*x));
451  G4double f0s = 6.0*(-x*y*(2.0-3.0*(y*y))
452  -2.0*(x*x)*(1.0-y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y));
453  G4double f1s = 3.0*((x*x)*y*(2.0-3.0*y-3.0*(y*y))
454  -(x*x*x)*y*(4.0+3.0*y));
455  G4double f2s = 1.5*((x*x*x)*(y*y)*(2.0+y));
456 
457  G4double f_1se = 12.0*(x*y*(1.0-y)+(x*x)*(2.0-3.0*y)
458  -2.0*(x*x*x));
459  G4double f0se = 6.0*(-(x*x)*(2.0-y-2.0*(y*y))
460  +(x*x*x)*(2.0+3.0*y));
461  G4double f1se = -3.0*(x*x*x)*y*(2.0+y);
462  G4double f2se = 0.0;
463 
464  G4double f_1sg = 12.0*((y*y)*(1.0-y)+x*y*(1.0-2.0*y)
465  -(x*x)*y);
466  G4double f0sg = 6.0*(-x*(y*y)*(2.0-3.0*y)-(x*x)*y*(1.0-4.0*y)
467  +(x*x*x)*y);
468  G4double f1sg = 3.0*((x*x)*(y*y)*(1.0-3.0*y)
469  -2.0*(x*x*x)*(y*y));
470  G4double f2sg = 1.5*(x*x*x)*(y*y*y);
471 
472  G4double f_1v = 8.0*((y*y)*(3.0-2.0*y)+6.0*x*y*(1.0-y)
473  +2.0*(x*x)*(3.0-4.0*y)-4.0*(x*x*x));
474  G4double f0v = 8.0*(-x*y*(3.0-y-(y*y))-(x*x)*(3.0-y-4.0*(y*y))
475  +2.0*(x*x*x)*(1.0+2.0*y));
476  G4double f1v = 2.0*((x*x)*y*(6.0-5.0*y-2.0*(y*y))
477  -2.0*(x*x*x)*y*(4.0+3.0*y));
478  G4double f2v = 2.0*(x*x*x)*(y*y)*(2.0+y);
479 
480  G4double f_1ve = 8.0*(x*y*(1.0-2.0*y)
481  +2.0*(x*x)*(1.0-3.0*y)-4.0*(x*x*x));
482  G4double f0ve = 4.0*(-(x*x)*(2.0-3.0*y-4.0*(y*y))
483  +2.0*(x*x*x)*(2.0+3.0*y));
484  G4double f1ve = -4.0*(x*x*x)*y*(2.0+y);
485  G4double f2ve = 0.0;
486 
487  G4double f_1vg = 8.0*((y*y)*(1.0-2.0*y)+x*y*(1.0-4.0*y)
488  -2.0*(x*x)*y);
489  G4double f0vg = 4.0*(2.0*x*(y*y)*(1.0+y)-(x*x)*y*(1.0-4.0*y)
490  +2.0*(x*x*x)*y);
491  G4double f1vg = 2.0*((x*x)*(y*y)*(1.0-2.0*y)
492  -4.0*(x*x*x)*(y*y));
493  G4double f2vg = 2.0*(x*x*x)*(y*y*y);
494 
495  G4double f_1t = 8.0*((y*y)*(3.0-y)+3.0*x*y*(2.0-y)
496  +2.0*(x*x)*(3.0-2.0*y)-2.0*(x*x*x));
497  G4double f0t = 4.0*(-x*y*(6.0+(y*y))
498  -2.0*(x*x)*(3.0+y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y));
499  G4double f1t = 2.0*((x*x)*y*(6.0-5.0*y+(y*y))
500  -(x*x*x)*y*(4.0+3.0*y));
501  G4double f2t = (x*x*x)*(y*y)*(2.0+y);
502 
503  G4double f_1te = -8.0*(x*y*(1.0+3.0*y)+(x*x)*(2.0+3.0*y)
504  +2.0*(x*x*x));
505  G4double f0te = 4.0*((x*x)*(2.0+3.0*y+4.0*(y*y))
506  +(x*x*x)*(2.0+3.0*y));
507  G4double f1te = -2.0*(x*x*x)*y*(2.0+y);
508  G4double f2te = 0.0;
509 
510  G4double f_1tg = -8.0*((y*y)*(1.0+y)+x*y+(x*x)*y);
511  G4double f0tg = 4.0*(x*(y*y)*(2.0-y)+(x*x)*y*(1.0+2.0*y)
512  +(x*x*x)*y);
513  G4double f1tg = -2.0*((x*x)*(y*y)*(1.0-y)+2.0*(x*x*x)*y);
514  G4double f2tg = (x*x*x)*(y*y*y);
515 
516  G4double term = delta+2.0*(me*me)/((mu*mu)*(x*x));
517  term = 1.0/term;
518 
519  G4double nss = term*f_1s+f0s+delta*f1s+(delta*delta)*f2s;
520  G4double nv = term*f_1v+f0v+delta*f1v+(delta*delta)*f2v;
521  G4double nt = term*f_1t+f0t+delta*f1t+(delta*delta)*f2t;
522 
523  G4double nse = term*f_1se+f0se+delta*f1se+(delta*delta)*f2se;
524  G4double nve = term*f_1ve+f0ve+delta*f1ve+(delta*delta)*f2ve;
525  G4double nte = term*f_1te+f0te+delta*f1te+(delta*delta)*f2te;
526 
527  G4double nsg = term*f_1sg+f0sg+delta*f1sg+(delta*delta)*f2sg;
528  G4double nvg = term*f_1vg+f0vg+delta*f1vg+(delta*delta)*f2vg;
529  G4double ntg = term*f_1tg+f0tg+delta*f1tg+(delta*delta)*f2tg;
530 
531  G4double term1 = nv;
532  G4double term2 = 2.0*nss+nv-nt;
533  G4double term3 = 2.0*nss-2.0*nv+nt;
534 
535  G4double term1e = 1.0/3.0*(1.0-4.0/3.0*del);
536  G4double term2e = 2.0*nse+5.0*nve-nte;
537  G4double term3e = 2.0*nse-2.0*nve+nte;
538 
539  G4double term1g = 1.0/3.0*(1.0-4.0/3.0*del);
540  G4double term2g = 2.0*nsg+5.0*nvg-ntg;
541  G4double term3g = 2.0*nsg-2.0*nvg+ntg;
542 
543  G4double som00 = term1+(1.0-4.0/3.0*rho)*term2+eps*term3;
544  G4double som01 = Pmu*ksi*(cthetaE*(nve-term1e*term2e+kap*term3e)
545  +cthetaG*(nvg-term1g*term2g+kap*term3g));
546 
547  G4double som0 = (som00+som01)/y;
548  som0 = fine_structure_const/8./(twopi*twopi*twopi)*som0;
549 
550 // G4cout << x << " " << y << " " << som00 << " "
551 // << som01 << " " << som0 << G4endl;
552 
553  return som0;
554 }
Generated on Sat May 25 2013 14:33:23 for Geant4 by   doxygen 1.8.4