ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TauLeptonicDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TauLeptonicDecayChannel.cc
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 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class header file
31 //
32 // History: first implementation, based on object model of
33 // 30 May 1997 H.Kurashige
34 //
35 // Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige
36 // ------------------------------------------------------------
37 
38 #include "G4ParticleDefinition.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4DecayProducts.hh"
42 #include "G4VDecayChannel.hh"
44 #include "Randomize.hh"
45 #include "G4LorentzVector.hh"
46 #include "G4LorentzRotation.hh"
47 
48 
51 {
52 }
53 
54 
56  const G4String& theParentName,
57  G4double theBR,
58  const G4String& theLeptonName)
59  :G4VDecayChannel("Tau Leptonic Decay",1)
60 {
61  // set names for daughter particles
62  if (theParentName == "tau+") {
63  SetBR(theBR);
64  SetParent("tau+");
66  if ((theLeptonName=="e-"||theLeptonName=="e+")){
67  SetDaughter(0, "e+");
68  SetDaughter(1, "nu_e");
69  SetDaughter(2, "anti_nu_tau");
70  } else {
71  SetDaughter(0, "mu+");
72  SetDaughter(1, "nu_mu");
73  SetDaughter(2, "anti_nu_tau");
74  }
75  } else if (theParentName == "tau-") {
76  SetBR(theBR);
77  SetParent("tau-");
79  if ((theLeptonName=="e-"||theLeptonName=="e+")){
80  SetDaughter(0, "e-");
81  SetDaughter(1, "anti_nu_e");
82  SetDaughter(2, "nu_tau");
83  } else {
84  SetDaughter(0, "mu-");
85  SetDaughter(1, "anti_nu_mu");
86  SetDaughter(2, "nu_tau");
87  }
88  } else {
89 #ifdef G4VERBOSE
90  if (GetVerboseLevel()>0) {
91  G4cout << "G4TauLeptonicDecayChannel:: constructor :";
92  G4cout << " parent particle is not tau but ";
93  G4cout << theParentName << G4endl;
94  }
95 #endif
96  }
97 }
98 
100 {
101 }
102 
104  G4VDecayChannel(right)
105 {
106 }
107 
109 {
110  if (this != &right) {
112  verboseLevel = right.verboseLevel;
113  rbranch = right.rbranch;
114 
115  // copy parent name
116  parent_name = new G4String(*right.parent_name);
117 
118  // clear daughters_name array
120 
121  // recreate array
123  if ( numberOfDaughters >0 ) {
126  //copy daughters name
127  for (G4int index=0; index < numberOfDaughters; index++) {
128  daughters_name[index] = new G4String(*right.daughters_name[index]);
129  }
130  }
131  }
132  return *this;
133 }
134 
136 {
137  // this version neglects muon polarization
138  // assumes the pure V-A coupling
139  // gives incorrect energy spectrum for neutrinos
140 #ifdef G4VERBOSE
141  if (GetVerboseLevel()>1) G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
142 #endif
143 
146 
147  // parent mass
148  G4double parentmass = G4MT_parent->GetPDGMass();
149 
150  //daughters'mass
151  const G4int N_DAUGHTER=3;
152  G4double daughtermass[N_DAUGHTER];
153  for (G4int index=0; index<N_DAUGHTER; index++){
154  daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
155  }
156 
157  //create parent G4DynamicParticle at rest
158  G4ThreeVector dummy;
159  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
160  //create G4Decayproducts
161  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
162  delete parentparticle;
163 
164  // calculate daughter momentum
165  G4double daughtermomentum[N_DAUGHTER];
166 
167  // calcurate lepton momentum
168  G4double pmax = (parentmass*parentmass-daughtermass[0]*daughtermass[0])/2./parentmass;
169  G4double p, e;
170  G4double r;
171  const size_t MAX_LOOP=10000;
172  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
173  // determine momentum/energy
174  r = G4UniformRand();
175  p = pmax*G4UniformRand();
176  e = std::sqrt(p*p + daughtermass[0]*daughtermass[0]);
177  if (r < spectrum(p,e,parentmass,daughtermass[0]) ) break;
178  }
179 
180  //create daughter G4DynamicParticle
181  // daughter 0 (lepton)
182  daughtermomentum[0] = p;
183  G4double costheta, sintheta, phi, sinphi, cosphi;
184  costheta = 2.*G4UniformRand()-1.0;
185  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
186  phi = twopi*G4UniformRand()*rad;
187  sinphi = std::sin(phi);
188  cosphi = std::cos(phi);
189  G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
190  G4DynamicParticle * daughterparticle
191  = new G4DynamicParticle( G4MT_daughters[0], direction0*daughtermomentum[0]);
192  products->PushProducts(daughterparticle);
193 
194  // daughter 1 ,2 (nutrinos)
195  // create neutrinos in the C.M frame of two neutrinos
196  G4double energy2 = parentmass-e;
197  G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
198  G4double beta = -1.0*daughtermomentum[0]/energy2;
199  G4double costhetan = 2.*G4UniformRand()-1.0;
200  G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
201  G4double phin = twopi*G4UniformRand()*rad;
202  G4double sinphin = std::sin(phin);
203  G4double cosphin = std::cos(phin);
204 
205  G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
206  G4DynamicParticle * daughterparticle1
207  = new G4DynamicParticle( G4MT_daughters[1], direction1*(vmass/2.));
208  G4DynamicParticle * daughterparticle2
209  = new G4DynamicParticle( G4MT_daughters[2], direction1*(-1.0*vmass/2.));
210 
211  // boost to the muon rest frame
212  G4LorentzVector p4;
213  p4 = daughterparticle1->Get4Momentum();
214  p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
215  daughterparticle1->Set4Momentum(p4);
216  p4 = daughterparticle2->Get4Momentum();
217  p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
218  daughterparticle2->Set4Momentum(p4);
219  products->PushProducts(daughterparticle1);
220  products->PushProducts(daughterparticle2);
221  daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
222  daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
223 
224 
225  // output message
226 #ifdef G4VERBOSE
227  if (GetVerboseLevel()>1) {
228  G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
229  G4cout << " create decay products in rest frame " <<G4endl;
230  products->DumpInfo();
231  }
232 #endif
233  return products;
234 }
235 
236 
237 
238 
240  G4double e,
241  G4double mtau,
242  G4double ml)
243 {
244  G4double f1;
245  f1 = 3.0*e*(mtau*mtau+ml*ml)-4.0*mtau*e*e-2.0*mtau*ml*ml;
246  return p*(f1)/(mtau*mtau*mtau*mtau)/(0.6);
247 }
248 
249 
250