source: trunk/source/particles/management/src/G4MuonRadiativeDecayChannelWithSpin.cc @ 850

Last change on this file since 850 was 824, checked in by garnier, 16 years ago

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