Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuonDecayChannelWithSpin.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 // 17 August 2004 P.Gumplinger and T.MacPhail
31 // samples Michel spectrum including 1st order
32 // radiative corrections
33 // Reference: Florian Scheck "Muon Physics", in Physics Reports
34 // (Review Section of Physics Letters) 44, No. 4 (1978)
35 // 187-248. North-Holland Publishing Company, Amsterdam
36 // at page 210 cc.
37 //
38 // W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25.
39 //
40 // ------------------------------------------------------------
41 //
43 
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "Randomize.hh"
47 
48 #include "G4DecayProducts.hh"
49 #include "G4LorentzVector.hh"
50 
53  parent_polarization(),
54  EMMU( 0.*MeV),
55  EMASS( 0.*MeV)
56 {
57 }
58 
60  G4double theBR)
61  : G4MuonDecayChannel(theParentName,theBR),
62  parent_polarization(),
63  EMMU( 0.*MeV),
64  EMASS( 0.*MeV)
65 {
66 }
67 
69 {
70 }
71 
73  G4MuonDecayChannel(right)
74 {
75  parent_polarization = right.parent_polarization;
76  EMMU = right.EMMU;
77  EMASS = right.EMASS;
78 }
79 
81 {
82  if (this != &right) {
84  verboseLevel = right.verboseLevel;
85  rbranch = right.rbranch;
86 
87  // copy parent name
88  parent_name = new G4String(*right.parent_name);
89 
90  // clear daughters_name array
92 
93  // recreate array
95  if ( numberOfDaughters >0 ) {
98  //copy daughters name
99  for (G4int index=0; index < numberOfDaughters; index++) {
101  }
102  }
103  parent_polarization = right.parent_polarization;
104  EMMU = right.EMMU;
105  EMASS = right.EMASS;
106  }
107  return *this;
108 }
109 
110 
112 {
113  // This version assumes V-A coupling with 1st order radiative correctons,
114  // the standard model Michel parameter values, but
115  // gives incorrect energy spectrum for neutrinos
116 
117 #ifdef G4VERBOSE
118  if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
119 #endif
120 
121  if (parent == 0) FillParent();
122  if (daughters == 0) FillDaughters();
123 
124  // parent mass
125  G4double parentmass = parent->GetPDGMass();
126 
127  EMMU = parentmass;
128 
129  //daughters'mass
130  G4double daughtermass[3];
131  G4double sumofdaughtermass = 0.0;
132  for (G4int index=0; index<3; index++){
133  daughtermass[index] = daughters[index]->GetPDGMass();
134  sumofdaughtermass += daughtermass[index];
135  }
136 
137  EMASS = daughtermass[0];
138 
139  //create parent G4DynamicParticle at rest
140  G4ThreeVector dummy;
141  G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
142  //create G4Decayproducts
143  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
144  delete parentparticle;
145 
146  // calcurate electron energy
147 
148  G4double michel_rho = 0.75; //Standard Model Michel rho
149  G4double michel_delta = 0.75; //Standard Model Michel delta
150  G4double michel_xsi = 1.00; //Standard Model Michel xsi
151  G4double michel_eta = 0.00; //Standard Model eta
152 
153  G4double rndm, x, ctheta;
154 
155  G4double FG;
156  G4double FG_max = 2.00;
157 
158  G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
159  G4double x0 = EMASS/W_mue;
160 
161  G4double x0_squared = x0*x0;
162 
163  // ***************************************************
164  // x0 <= x <= 1. and -1 <= y <= 1
165  //
166  // F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g(x)*y
167  // ***************************************************
168 
169  // ***** sampling F(x,y) directly (brute force) *****
170 
171  do{
172 
173  // Sample the positron energy by sampling from F
174 
175  rndm = G4UniformRand();
176 
177  x = x0 + rndm*(1.-x0);
178 
179  G4double x_squared = x*x;
180 
181  G4double F_IS, F_AS, G_IS, G_AS;
182 
183  F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
184  F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared));
185 
186  G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
187  G_IS = G_IS + michel_eta*(1.-x)*x0;
188 
189  G_AS = 3.*(michel_xsi-1.)*(1.-x);
190  G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared));
191  G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
192 
193  F_IS = F_IS + G_IS;
194  F_AS = F_AS + G_AS;
195 
196  // *** Radiative Corrections ***
197 
198  G4double R_IS = F_c(x,x0);
199 
200  G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
201 
202  // *** Radiative Corrections ***
203 
204  G4double R_AS = F_theta(x,x0);
205 
206  rndm = G4UniformRand();
207 
208  ctheta = 2.*rndm-1.;
209 
210  G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
211 
212  FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
213 
214  if(FG>FG_max){
215  G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<<G4endl;
216  FG_max = FG;
217  }
218 
219  rndm = G4UniformRand();
220 
221  }while(FG<rndm*FG_max);
222 
223  G4double energy = x * W_mue;
224 
225  rndm = G4UniformRand();
226 
227  G4double phi = twopi * rndm;
228 
229  if(energy < EMASS) energy = EMASS;
230 
231  // calculate daughter momentum
232  G4double daughtermomentum[3];
233 
234  daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
235 
236  G4double stheta = std::sqrt(1.-ctheta*ctheta);
237  G4double cphi = std::cos(phi);
238  G4double sphi = std::sin(phi);
239 
240  //Coordinates of the decay positron with respect to the muon spin
241 
242  G4double px = stheta*cphi;
243  G4double py = stheta*sphi;
244  G4double pz = ctheta;
245 
246  G4ThreeVector direction0(px,py,pz);
247 
248  direction0.rotateUz(parent_polarization);
249 
250  G4DynamicParticle * daughterparticle0
251  = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
252 
253  products->PushProducts(daughterparticle0);
254 
255 
256  // daughter 1 ,2 (neutrinos)
257  // create neutrinos in the C.M frame of two neutrinos
258  G4double energy2 = parentmass*(1.0 - x/2.0);
259  G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
260  G4double beta = -1.0*daughtermomentum[0]/energy2;
261  G4double costhetan = 2.*G4UniformRand()-1.0;
262  G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
263  G4double phin = twopi*G4UniformRand()*rad;
264  G4double sinphin = std::sin(phin);
265  G4double cosphin = std::cos(phin);
266 
267  G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
268  G4DynamicParticle * daughterparticle1
269  = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
270  G4DynamicParticle * daughterparticle2
271  = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
272 
273  // boost to the muon rest frame
274  G4LorentzVector p4;
275  p4 = daughterparticle1->Get4Momentum();
276  p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
277  daughterparticle1->Set4Momentum(p4);
278  p4 = daughterparticle2->Get4Momentum();
279  p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
280  daughterparticle2->Set4Momentum(p4);
281  products->PushProducts(daughterparticle1);
282  products->PushProducts(daughterparticle2);
283  daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
284  daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
285 
286  // output message
287 #ifdef G4VERBOSE
288  if (GetVerboseLevel()>1) {
289  G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
290  G4cout << " create decay products in rest frame " <<G4endl;
291  products->DumpInfo();
292  }
293 #endif
294  return products;
295 }
296 
297 G4double G4MuonDecayChannelWithSpin::R_c(G4double x){
298 
299  G4int n_max = (int)(100.*x);
300 
301  if(n_max<10)n_max=10;
302 
303  G4double L2 = 0.0;
304 
305  for(G4int n=1; n<=n_max; n++){
306  L2 += std::pow(x,n)/(n*n);
307  }
308 
309  G4double omega = std::log(EMMU/EMASS);
310 
311  G4double r_c;
312 
313  r_c = 2.*L2-(pi*pi/3.)-2.;
314  r_c = r_c + omega * (1.5+2.*std::log((1.-x)/x));
315  r_c = r_c - std::log(x)*(2.*std::log(x)-1.);
316  r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x);
317 
318  return r_c;
319 }