Geant4  10.03.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 
46 
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "Randomize.hh"
50 #include "G4DecayProducts.hh"
51 #include "G4LorentzVector.hh"
52 
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+");
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-");
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 
92 {
93 }
94 
96  G4VDecayChannel(right)
97 {
98 }
99 
101 {
102  if (this != &right) {
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
112 
113  // recreate array
115  if ( numberOfDaughters >0 ) {
118  //copy daughters name
119  for (G4int index=0; index < numberOfDaughters; index++) {
120  daughters_name[index] = new G4String(*right.daughters_name[index]);
121  }
122  }
124  }
125  return *this;
126 }
127 
128 
130 {
131 
132 #ifdef G4VERBOSE
133  if (GetVerboseLevel()>1)
134  G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
135 #endif
136 
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*(1.0 - (x+y)/2.0);
377 
378  G4double vmass2 = energy2*energy2 -
379  (daughtermomentum[0]*direction0+daughtermomentum[1]*direction1)*
380  (daughtermomentum[0]*direction0+daughtermomentum[1]*direction1);
381  G4double vmass = std::sqrt(vmass2);
382 
383  G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2;
384  beta = -1.0 * std::min(beta,0.99);
385 
386  G4double costhetan = 2.*G4UniformRand()-1.0;
387  G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
388  G4double phin = twopi*G4UniformRand()*rad;
389  G4double sinphin = std::sin(phin);
390  G4double cosphin = std::cos(phin);
391 
392  G4ThreeVector direction2(sinthetan*cosphin,sinthetan*sinphin,costhetan);
393 
394  G4DynamicParticle * daughterparticle2
395  = new G4DynamicParticle( G4MT_daughters[2], direction2*(vmass/2.));
396  G4DynamicParticle * daughterparticle3
397  = new G4DynamicParticle( G4MT_daughters[3], direction2*(-1.0*vmass/2.));
398 
399  // boost to the muon rest frame
400 
401  G4ThreeVector direction34(direction0.x()+direction1.x(),
402  direction0.y()+direction1.y(),
403  direction0.z()+direction1.z());
404  direction34 = direction34.unit();
405 
406  G4LorentzVector p4 = daughterparticle2->Get4Momentum();
407  p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
408  daughterparticle2->Set4Momentum(p4);
409 
410  p4 = daughterparticle3->Get4Momentum();
411  p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
412  daughterparticle3->Set4Momentum(p4);
413 
414  products->PushProducts(daughterparticle2);
415  products->PushProducts(daughterparticle3);
416 
417  daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
418  daughtermomentum[3] = daughterparticle3->GetTotalMomentum();
419 
420 // output message
421 #ifdef G4VERBOSE
422  if (GetVerboseLevel()>1) {
423  G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
424  G4cout << " create decay products in rest frame " <<G4endl;
425  products->DumpInfo();
426  }
427 #endif
428  return products;
429 }
430 
431 G4double G4MuonRadiativeDecayChannelWithSpin::fron(G4double Pmu,
432  G4double x,
433  G4double y,
434  G4double cthetaE,
435  G4double cthetaG,
436  G4double cthetaGE)
437 {
438  G4double mu = 105.65;
439  G4double me = 0.511;
440  G4double rho = 0.75;
441  G4double del = 0.75;
442  G4double eps = 0.0;
443  G4double kap = 0.0;
444  G4double ksi = 1.0;
445 
446  G4double delta = 1-cthetaGE;
447 
448 // Calculation of the functions f(x,y)
449 
450  G4double f_1s = 12.0*((y*y)*(1.0-y)+x*y*(2.0-3.0*y)
451  +2.0*(x*x)*(1.0-2.0*y)-2.0*(x*x*x));
452  G4double f0s = 6.0*(-x*y*(2.0-3.0*(y*y))
453  -2.0*(x*x)*(1.0-y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y));
454  G4double f1s = 3.0*((x*x)*y*(2.0-3.0*y-3.0*(y*y))
455  -(x*x*x)*y*(4.0+3.0*y));
456  G4double f2s = 1.5*((x*x*x)*(y*y)*(2.0+y));
457 
458  G4double f_1se = 12.0*(x*y*(1.0-y)+(x*x)*(2.0-3.0*y)
459  -2.0*(x*x*x));
460  G4double f0se = 6.0*(-(x*x)*(2.0-y-2.0*(y*y))
461  +(x*x*x)*(2.0+3.0*y));
462  G4double f1se = -3.0*(x*x*x)*y*(2.0+y);
463  G4double f2se = 0.0;
464 
465  G4double f_1sg = 12.0*((y*y)*(1.0-y)+x*y*(1.0-2.0*y)
466  -(x*x)*y);
467  G4double f0sg = 6.0*(-x*(y*y)*(2.0-3.0*y)-(x*x)*y*(1.0-4.0*y)
468  +(x*x*x)*y);
469  G4double f1sg = 3.0*((x*x)*(y*y)*(1.0-3.0*y)
470  -2.0*(x*x*x)*(y*y));
471  G4double f2sg = 1.5*(x*x*x)*(y*y*y);
472 
473  G4double f_1v = 8.0*((y*y)*(3.0-2.0*y)+6.0*x*y*(1.0-y)
474  +2.0*(x*x)*(3.0-4.0*y)-4.0*(x*x*x));
475  G4double f0v = 8.0*(-x*y*(3.0-y-(y*y))-(x*x)*(3.0-y-4.0*(y*y))
476  +2.0*(x*x*x)*(1.0+2.0*y));
477  G4double f1v = 2.0*((x*x)*y*(6.0-5.0*y-2.0*(y*y))
478  -2.0*(x*x*x)*y*(4.0+3.0*y));
479  G4double f2v = 2.0*(x*x*x)*(y*y)*(2.0+y);
480 
481  G4double f_1ve = 8.0*(x*y*(1.0-2.0*y)
482  +2.0*(x*x)*(1.0-3.0*y)-4.0*(x*x*x));
483  G4double f0ve = 4.0*(-(x*x)*(2.0-3.0*y-4.0*(y*y))
484  +2.0*(x*x*x)*(2.0+3.0*y));
485  G4double f1ve = -4.0*(x*x*x)*y*(2.0+y);
486  G4double f2ve = 0.0;
487 
488  G4double f_1vg = 8.0*((y*y)*(1.0-2.0*y)+x*y*(1.0-4.0*y)
489  -2.0*(x*x)*y);
490  G4double f0vg = 4.0*(2.0*x*(y*y)*(1.0+y)-(x*x)*y*(1.0-4.0*y)
491  +2.0*(x*x*x)*y);
492  G4double f1vg = 2.0*((x*x)*(y*y)*(1.0-2.0*y)
493  -4.0*(x*x*x)*(y*y));
494  G4double f2vg = 2.0*(x*x*x)*(y*y*y);
495 
496  G4double f_1t = 8.0*((y*y)*(3.0-y)+3.0*x*y*(2.0-y)
497  +2.0*(x*x)*(3.0-2.0*y)-2.0*(x*x*x));
498  G4double f0t = 4.0*(-x*y*(6.0+(y*y))
499  -2.0*(x*x)*(3.0+y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y));
500  G4double f1t = 2.0*((x*x)*y*(6.0-5.0*y+(y*y))
501  -(x*x*x)*y*(4.0+3.0*y));
502  G4double f2t = (x*x*x)*(y*y)*(2.0+y);
503 
504  G4double f_1te = -8.0*(x*y*(1.0+3.0*y)+(x*x)*(2.0+3.0*y)
505  +2.0*(x*x*x));
506  G4double f0te = 4.0*((x*x)*(2.0+3.0*y+4.0*(y*y))
507  +(x*x*x)*(2.0+3.0*y));
508  G4double f1te = -2.0*(x*x*x)*y*(2.0+y);
509  G4double f2te = 0.0;
510 
511  G4double f_1tg = -8.0*((y*y)*(1.0+y)+x*y+(x*x)*y);
512  G4double f0tg = 4.0*(x*(y*y)*(2.0-y)+(x*x)*y*(1.0+2.0*y)
513  +(x*x*x)*y);
514  G4double f1tg = -2.0*((x*x)*(y*y)*(1.0-y)+2.0*(x*x*x)*y);
515  G4double f2tg = (x*x*x)*(y*y*y);
516 
517  G4double term = delta+2.0*(me*me)/((mu*mu)*(x*x));
518  term = 1.0/term;
519 
520  G4double nss = term*f_1s+f0s+delta*f1s+(delta*delta)*f2s;
521  G4double nv = term*f_1v+f0v+delta*f1v+(delta*delta)*f2v;
522  G4double nt = term*f_1t+f0t+delta*f1t+(delta*delta)*f2t;
523 
524  G4double nse = term*f_1se+f0se+delta*f1se+(delta*delta)*f2se;
525  G4double nve = term*f_1ve+f0ve+delta*f1ve+(delta*delta)*f2ve;
526  G4double nte = term*f_1te+f0te+delta*f1te+(delta*delta)*f2te;
527 
528  G4double nsg = term*f_1sg+f0sg+delta*f1sg+(delta*delta)*f2sg;
529  G4double nvg = term*f_1vg+f0vg+delta*f1vg+(delta*delta)*f2vg;
530  G4double ntg = term*f_1tg+f0tg+delta*f1tg+(delta*delta)*f2tg;
531 
532  G4double term1 = nv;
533  G4double term2 = 2.0*nss+nv-nt;
534  G4double term3 = 2.0*nss-2.0*nv+nt;
535 
536  G4double term1e = 1.0/3.0*(1.0-4.0/3.0*del);
537  G4double term2e = 2.0*nse+5.0*nve-nte;
538  G4double term3e = 2.0*nse-2.0*nve+nte;
539 
540  G4double term1g = 1.0/3.0*(1.0-4.0/3.0*del);
541  G4double term2g = 2.0*nsg+5.0*nvg-ntg;
542  G4double term3g = 2.0*nsg-2.0*nvg+ntg;
543 
544  G4double som00 = term1+(1.0-4.0/3.0*rho)*term2+eps*term3;
545  G4double som01 = Pmu*ksi*(cthetaE*(nve-term1e*term2e+kap*term3e)
546  +cthetaG*(nvg-term1g*term2g+kap*term3g));
547 
548  G4double som0 = (som00+som01)/y;
549  som0 = fine_structure_const/8./(twopi*twopi*twopi)*som0;
550 
551 // G4cout << x << " " << y << " " << som00 << " "
552 // << som01 << " " << som0 << G4endl;
553 
554  return som0;
555 }
void CheckAndFillDaughters()
void SetBR(G4double value)
G4MuonRadiativeDecayChannelWithSpin(const G4String &theParentName, G4double theBR)
const G4String & GetParentName() const
G4int PushProducts(G4DynamicParticle *aParticle)
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
static const G4double eps
tuple x
Definition: test.py:50
static constexpr double rad
Definition: G4SIunits.hh:149
int G4int
Definition: G4Types.hh:78
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double GetTotalMomentum() const
G4MuonRadiativeDecayChannelWithSpin & operator=(const G4MuonRadiativeDecayChannelWithSpin &)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4String kinematics_name
void DumpInfo() const
void SetNumberOfDaughters(G4int value)
HepLorentzVector & boost(double, double, double)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4String * parent_name
G4ThreeVector parent_polarization
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
G4int GetVerboseLevel() const
G4double GetPDGMass() const
Hep3Vector unit() const
void SetParent(const G4ParticleDefinition *particle_type)
void SetDaughter(G4int anIndex, const G4ParticleDefinition *particle_type)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4String ** daughters_name