Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/particles/management/src/G4MuonRadiativeDecayChannelWithSpin.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /particles/management/src/G4MuonRadiativeDecayChannelWithSpin.cc (Version 11.3.0) and /particles/management/src/G4MuonRadiativeDecayChannelWithSpin.cc (Version 8.1.p2)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // G4MuonRadiativeDecayChannelWithSpin class i    
 27 //                                                
 28 // References:                                    
 29 // - TRIUMF/TWIST Technote TN-55:                 
 30 //   "Radiative muon decay" by P. Depommier an    
 31 // - Yoshitaka Kuno and Yasuhiro Okada            
 32 //   "Muon Decays and Physics Beyond the Stand    
 33 //   Rev. Mod. Phys. 73, 151 (2001)               
 34 //                                                
 35 // Author: P.Gumplinger - Triumf, 25 July 2007    
 36 // Revision: D.Mingming - Center for HEP, Tsin    
 37 // -------------------------------------------    
 38                                                   
 39 #include "G4MuonRadiativeDecayChannelWithSpin.    
 40                                                   
 41 #include "G4DecayProducts.hh"                     
 42 #include "G4LorentzVector.hh"                     
 43 #include "G4PhysicalConstants.hh"                 
 44 #include "G4SystemOfUnits.hh"                     
 45 #include "Randomize.hh"                           
 46                                                   
 47 G4MuonRadiativeDecayChannelWithSpin::G4MuonRad    
 48   const G4String& theParentName, G4double theB    
 49   : G4VDecayChannel("Radiative Muon Decay", 1)    
 50 {                                                 
 51   // set names for daughter particles             
 52   if (theParentName == "mu+") {                   
 53     SetBR(theBR);                                 
 54     SetParent("mu+");                             
 55     SetNumberOfDaughters(4);                      
 56     SetDaughter(0, "e+");                         
 57     SetDaughter(1, "gamma");                      
 58     SetDaughter(2, "nu_e");                       
 59     SetDaughter(3, "anti_nu_mu");                 
 60   }                                               
 61   else if (theParentName == "mu-") {              
 62     SetBR(theBR);                                 
 63     SetParent("mu-");                             
 64     SetNumberOfDaughters(4);                      
 65     SetDaughter(0, "e-");                         
 66     SetDaughter(1, "gamma");                      
 67     SetDaughter(2, "anti_nu_e");                  
 68     SetDaughter(3, "nu_mu");                      
 69   }                                               
 70   else {                                          
 71 #ifdef G4VERBOSE                                  
 72     if (GetVerboseLevel() > 0) {                  
 73       G4cout << "G4RadiativeMuonDecayChannel::    
 74       G4cout << " parent particle is not muon     
 75       G4cout << theParentName << G4endl;          
 76     }                                             
 77 #endif                                            
 78   }                                               
 79 }                                                 
 80                                                   
 81 G4MuonRadiativeDecayChannelWithSpin&              
 82 G4MuonRadiativeDecayChannelWithSpin::operator=    
 83 {                                                 
 84   if (this != &right) {                           
 85     kinematics_name = right.kinematics_name;      
 86     verboseLevel = right.verboseLevel;            
 87     rbranch = right.rbranch;                      
 88                                                   
 89     // copy parent name                           
 90     parent_name = new G4String(*right.parent_n    
 91                                                   
 92     // clear daughters_name array                 
 93     ClearDaughtersName();                         
 94                                                   
 95     // recreate array                             
 96     numberOfDaughters = right.numberOfDaughter    
 97     if (numberOfDaughters > 0) {                  
 98       if (daughters_name != nullptr) ClearDaug    
 99       daughters_name = new G4String*[numberOfD    
100       // copy daughters name                      
101       for (G4int index = 0; index < numberOfDa    
102         daughters_name[index] = new G4String(*    
103       }                                           
104     }                                             
105     parent_polarization = right.parent_polariz    
106   }                                               
107   return *this;                                   
108 }                                                 
109                                                   
110 G4DecayProducts* G4MuonRadiativeDecayChannelWi    
111 {                                                 
112 #ifdef G4VERBOSE                                  
113   if (GetVerboseLevel() > 1) G4cout << "G4Muon    
114 #endif                                            
115                                                   
116   CheckAndFillParent();                           
117   CheckAndFillDaughters();                        
118                                                   
119   // parent mass                                  
120   G4double parentmass = G4MT_parent->GetPDGMas    
121                                                   
122   G4double EMMU = parentmass;                     
123                                                   
124   // daughters'mass                               
125   G4double daughtermass[4];                       
126   // G4double sumofdaughtermass = 0.0;            
127   for (G4int index = 0; index < 4; ++index) {     
128     daughtermass[index] = G4MT_daughters[index    
129     // sumofdaughtermass += daughtermass[index    
130   }                                               
131                                                   
132   G4double EMASS = daughtermass[0];               
133                                                   
134   // create parent G4DynamicParticle at rest      
135   G4ThreeVector dummy;                            
136   auto parentparticle = new G4DynamicParticle(    
137                                                   
138   // create G4Decayproducts                       
139   auto products = new G4DecayProducts(*parentp    
140   delete parentparticle;                          
141                                                   
142   G4double eps = EMASS / EMMU;                    
143                                                   
144   G4double som0, x, y, xx, yy, zz;                
145   G4double cthetaE, cthetaG, cthetaGE, phiE, p    
146   const std::size_t MAX_LOOP = 10000;             
147                                                   
148   for (std::size_t loop_counter1 = 0; loop_cou    
149     for (std::size_t loop_counter2 = 0; loop_c    
150       // -------------------------------------    
151       // Build two vectors of random length an    
152       // positron and the photon.                 
153       // x/y is the length of the vector, xx,     
154       // phi is the azimutal angle, theta the     
155       // -------------------------------------    
156                                                   
157       // For the positron                         
158       //                                          
159       x = G4UniformRand();                        
160                                                   
161       rn3dim(xx, yy, zz, x);                      
162                                                   
163       if (std::fabs((xx * xx) + (yy * yy) + (z    
164         G4cout << "Norm of x not correct" << G    
165       }                                           
166                                                   
167       phiE = atan4(xx, yy);                       
168       cthetaE = zz / x;                           
169       G4double sthetaE = std::sqrt((xx * xx) +    
170                                                   
171       // What you get:                            
172       //                                          
173       // x       = positron energy                
174       // phiE    = azimutal angle of positron     
175       // cthetaE = cosine of polar angle of po    
176       // sthetaE = sine of polar angle of posi    
177       //                                          
178       //// G4cout << " x, xx, yy, zz " << x  <    
179       ////                             << yy <    
180       //// G4cout << " phiE, cthetaE, sthetaE     
181       ////                                        
182       ////                                        
183                                                   
184       // For the photon                           
185       //                                          
186       y = G4UniformRand();                        
187                                                   
188       rn3dim(xx, yy, zz, y);                      
189                                                   
190       if (std::fabs((xx * xx) + (yy * yy) + (z    
191         G4cout << " Norm of y not correct " <<    
192       }                                           
193                                                   
194       phiG = atan4(xx, yy);                       
195       cthetaG = zz / y;                           
196       G4double sthetaG = std::sqrt((xx * xx) +    
197                                                   
198       // What you get:                            
199       //                                          
200       // y       = photon energy                  
201       // phiG    = azimutal angle of photon mo    
202       // cthetaG = cosine of polar angle of ph    
203       // sthetaG = sine of polar angle of phot    
204       //                                          
205       //// G4cout << " y, xx, yy, zz " << y  <    
206       ////                             << yy <    
207       //// G4cout << " phiG, cthetaG, sthetaG     
208       ////                                        
209       ////                                        
210                                                   
211       //      Calculate the angle between posi    
212       //                                          
213       cthetaGE = cthetaE * cthetaG + sthetaE *    
214                                                   
215       //// G4cout << x << " " << cthetaE << "     
216       ////        << y << " " << cthetaG << "     
217       ////        << cthetaGE                     
218                                                   
219       G4double term0 = eps * eps;                 
220       G4double term1 = x * ((1.0 - eps) * (1.0    
221       G4double beta =                             
222         std::sqrt(x * ((1.0 - eps) * (1.0 - ep    
223         / term1;                                  
224       G4double delta = 1.0 - beta * cthetaGE;     
225                                                   
226       G4double term3 = y * (1.0 - (eps * eps))    
227       G4double term6 = term1 * delta * term3;     
228                                                   
229       G4double Qsqr = (1.0 - term1 - term3 + t    
230                                                   
231       // Check the kinematics.                    
232       //                                          
233       if (Qsqr >= 0.0 && Qsqr <= 1.0) break;      
234                                                   
235     }  // end loop count                          
236                                                   
237     // Do the calculation for -1 muon polariza    
238     //                                            
239     G4double Pmu = -1.0;                          
240     if (GetParentName() == "mu-") {               
241       Pmu = +1.0;                                 
242     }                                             
243                                                   
244     som0 = fron(Pmu, x, y, cthetaE, cthetaG, c    
245                                                   
246     // Sample the decay rate                      
247     //                                            
248     if (G4UniformRand() * 177.0 <= som0) break    
249   }                                               
250                                                   
251   G4double E = EMMU / 2. * (x * ((1. - eps) *     
252   G4double G = EMMU / 2. * y * (1. - eps * eps    
253                                                   
254   if (E < EMASS) E = EMASS;                       
255                                                   
256   // calculate daughter momentum                  
257   G4double daughtermomentum[4];                   
258                                                   
259   daughtermomentum[0] = std::sqrt(E * E - EMAS    
260                                                   
261   G4double sthetaE = std::sqrt(1. - cthetaE *     
262   G4double cphiE = std::cos(phiE);                
263   G4double sphiE = std::sin(phiE);                
264                                                   
265   // Coordinates of the decay positron with re    
266                                                   
267   G4double px = sthetaE * cphiE;                  
268   G4double py = sthetaE * sphiE;                  
269   G4double pz = cthetaE;                          
270                                                   
271   G4ThreeVector direction0(px, py, pz);           
272                                                   
273   direction0.rotateUz(parent_polarization);       
274                                                   
275   auto daughterparticle0 =                        
276     new G4DynamicParticle(G4MT_daughters[0], d    
277                                                   
278   products->PushProducts(daughterparticle0);      
279                                                   
280   daughtermomentum[1] = G;                        
281                                                   
282   G4double sthetaG = std::sqrt(1. - cthetaG *     
283   G4double cphiG = std::cos(phiG);                
284   G4double sphiG = std::sin(phiG);                
285                                                   
286   // Coordinates of the decay gamma with respe    
287                                                   
288   px = sthetaG * cphiG;                           
289   py = sthetaG * sphiG;                           
290   pz = cthetaG;                                   
291                                                   
292   G4ThreeVector direction1(px, py, pz);           
293                                                   
294   direction1.rotateUz(parent_polarization);       
295                                                   
296   auto daughterparticle1 =                        
297     new G4DynamicParticle(G4MT_daughters[1], d    
298                                                   
299   products->PushProducts(daughterparticle1);      
300                                                   
301   // daughter 3 ,4 (neutrinos)                    
302   // create neutrinos in the C.M frame of two     
303                                                   
304   G4double energy2 = parentmass - E - G;          
305                                                   
306   G4ThreeVector P34 = -1. * (daughtermomentum[    
307   G4double vmass2 = energy2 * energy2 - P34.ma    
308   G4double vmass = std::sqrt(vmass2);             
309                                                   
310   G4double costhetan = 2. * G4UniformRand() -     
311   G4double sinthetan = std::sqrt((1.0 - costhe    
312   G4double phin = twopi * G4UniformRand() * ra    
313   G4double sinphin = std::sin(phin);              
314   G4double cosphin = std::cos(phin);              
315                                                   
316   G4ThreeVector direction2(sinthetan * cosphin    
317                                                   
318   auto daughterparticle2 = new G4DynamicPartic    
319   auto daughterparticle3 =                        
320     new G4DynamicParticle(G4MT_daughters[3], d    
321                                                   
322   // boost to the muon rest frame                 
323   G4double beta = P34.mag() / energy2;            
324   G4ThreeVector direction34 = P34.unit();         
325                                                   
326   G4LorentzVector p4 = daughterparticle2->Get4    
327   p4.boost(direction34.x() * beta, direction34    
328   daughterparticle2->Set4Momentum(p4);            
329                                                   
330   p4 = daughterparticle3->Get4Momentum();         
331   p4.boost(direction34.x() * beta, direction34    
332   daughterparticle3->Set4Momentum(p4);            
333                                                   
334   products->PushProducts(daughterparticle2);      
335   products->PushProducts(daughterparticle3);      
336                                                   
337   daughtermomentum[2] = daughterparticle2->Get    
338   daughtermomentum[3] = daughterparticle3->Get    
339                                                   
340   // output message                               
341 #ifdef G4VERBOSE                                  
342   if (GetVerboseLevel() > 1) {                    
343     G4cout << "G4MuonRadiativeDecayChannelWith    
344     G4cout << " create decay products in rest     
345     G4double TT = daughterparticle0->GetTotalE    
346                   + daughterparticle2->GetTota    
347     G4cout << "e    :" << daughterparticle0->G    
348     G4cout << "gamma:" << daughterparticle1->G    
349     G4cout << "nu2  :" << daughterparticle2->G    
350     G4cout << "nu2  :" << daughterparticle3->G    
351     G4cout << "total:" << (TT - parentmass) /     
352     if (GetVerboseLevel() > 1) {                  
353       products->DumpInfo();                       
354     }                                             
355   }                                               
356 #endif                                            
357                                                   
358   return products;                                
359 }                                                 
360                                                   
361 G4double G4MuonRadiativeDecayChannelWithSpin::    
362                                                   
363                                                   
364 {                                                 
365   G4double mu = 105.65;                           
366   G4double me = 0.511;                            
367   G4double rho = 0.75;                            
368   G4double del = 0.75;                            
369   G4double eps = 0.0;                             
370   G4double kap = 0.0;                             
371   G4double ksi = 1.0;                             
372                                                   
373   G4double delta = 1 - cthetaGE;                  
374                                                   
375   // Calculation of the functions f(x,y)          
376                                                   
377   G4double f_1s = 12.0                            
378                   * ((y * y) * (1.0 - y) + x *    
379                      - 2.0 * (x * x * x));        
380   G4double f0s = 6.0                              
381                  * (-x * y * (2.0 - 3.0 * (y *    
382                     + 2.0 * (x * x * x) * (1.0    
383   G4double f1s =                                  
384     3.0 * ((x * x) * y * (2.0 - 3.0 * y - 3.0     
385   G4double f2s = 1.5 * ((x * x * x) * (y * y)     
386                                                   
387   G4double f_1se = 12.0 * (x * y * (1.0 - y) +    
388   G4double f0se = 6.0 * (-(x * x) * (2.0 - y -    
389   G4double f1se = -3.0 * (x * x * x) * y * (2.    
390   G4double f2se = 0.0;                            
391                                                   
392   G4double f_1sg = 12.0 * ((y * y) * (1.0 - y)    
393   G4double f0sg =                                 
394     6.0 * (-x * (y * y) * (2.0 - 3.0 * y) - (x    
395   G4double f1sg = 3.0 * ((x * x) * (y * y) * (    
396   G4double f2sg = 1.5 * (x * x * x) * (y * y *    
397                                                   
398   G4double f_1v = 8.0                             
399                   * ((y * y) * (3.0 - 2.0 * y)    
400                      + 2.0 * (x * x) * (3.0 -     
401   G4double f0v = 8.0                              
402                  * (-x * y * (3.0 - y - (y * y    
403                     + 2.0 * (x * x * x) * (1.0    
404   G4double f1v =                                  
405     2.0 * ((x * x) * y * (6.0 - 5.0 * y - 2.0     
406   G4double f2v = 2.0 * (x * x * x) * (y * y) *    
407                                                   
408   G4double f_1ve =                                
409     8.0 * (x * y * (1.0 - 2.0 * y) + 2.0 * (x     
410   G4double f0ve =                                 
411     4.0 * (-(x * x) * (2.0 - 3.0 * y - 4.0 * (    
412   G4double f1ve = -4.0 * (x * x * x) * y * (2.    
413   G4double f2ve = 0.0;                            
414                                                   
415   G4double f_1vg = 8.0 * ((y * y) * (1.0 - 2.0    
416   G4double f0vg =                                 
417     4.0 * (2.0 * x * (y * y) * (1.0 + y) - (x     
418   G4double f1vg = 2.0 * ((x * x) * (y * y) * (    
419   G4double f2vg = 2.0 * (x * x * x) * (y * y *    
420                                                   
421   G4double f_1t = 8.0                             
422                   * ((y * y) * (3.0 - y) + 3.0    
423                      - 2.0 * (x * x * x));        
424   G4double f0t = 4.0                              
425                  * (-x * y * (6.0 + (y * y)) -    
426                     + 2.0 * (x * x * x) * (1.0    
427   G4double f1t =                                  
428     2.0 * ((x * x) * y * (6.0 - 5.0 * y + (y *    
429   G4double f2t = (x * x * x) * (y * y) * (2.0     
430                                                   
431   G4double f_1te = -8.0 * (x * y * (1.0 + 3.0     
432   G4double f0te = 4.0 * ((x * x) * (2.0 + 3.0     
433   G4double f1te = -2.0 * (x * x * x) * y * (2.    
434   G4double f2te = 0.0;                            
435                                                   
436   G4double f_1tg = -8.0 * ((y * y) * (1.0 + y)    
437   G4double f0tg = 4.0 * (x * (y * y) * (2.0 -     
438   G4double f1tg = -2.0 * ((x * x) * (y * y) *     
439   G4double f2tg = (x * x * x) * (y * y * y);      
440                                                   
441   G4double term = delta + 2.0 * (me * me) / ((    
442   term = 1.0 / term;                              
443                                                   
444   G4double nss = term * f_1s + f0s + delta * f    
445   G4double nv = term * f_1v + f0v + delta * f1    
446   G4double nt = term * f_1t + f0t + delta * f1    
447                                                   
448   G4double nse = term * f_1se + f0se + delta *    
449   G4double nve = term * f_1ve + f0ve + delta *    
450   G4double nte = term * f_1te + f0te + delta *    
451                                                   
452   G4double nsg = term * f_1sg + f0sg + delta *    
453   G4double nvg = term * f_1vg + f0vg + delta *    
454   G4double ntg = term * f_1tg + f0tg + delta *    
455                                                   
456   G4double term1 = nv;                            
457   G4double term2 = 2.0 * nss + nv - nt;           
458   G4double term3 = 2.0 * nss - 2.0 * nv + nt;     
459                                                   
460   G4double term1e = 1.0 / 3.0 * (1.0 - 4.0 / 3    
461   G4double term2e = 2.0 * nse + 5.0 * nve - nt    
462   G4double term3e = 2.0 * nse - 2.0 * nve + nt    
463                                                   
464   G4double term1g = 1.0 / 3.0 * (1.0 - 4.0 / 3    
465   G4double term2g = 2.0 * nsg + 5.0 * nvg - nt    
466   G4double term3g = 2.0 * nsg - 2.0 * nvg + nt    
467                                                   
468   G4double som00 = term1 + (1.0 - 4.0 / 3.0 *     
469   G4double som01 = Pmu * ksi                      
470                    * (cthetaE * (nve - term1e     
471                       + cthetaG * (nvg - term1    
472                                                   
473   G4double som0 = (som00 + som01) / y;            
474   som0 = fine_structure_const / 8. / (twopi *     
475                                                   
476   return som0;                                    
477 }                                                 
478