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

Last change on this file since 1357 was 1340, checked in by garnier, 15 years ago

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