Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4TauLeptonicDecayChannel.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// $Id$
28//
29//
30// ------------------------------------------------------------
31// GEANT 4 class header file
32//
33// History: first implementation, based on object model of
34// 30 May 1997 H.Kurashige
35//
36// Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige
37// ------------------------------------------------------------
38
41#include "G4SystemOfUnits.hh"
42#include "G4DecayProducts.hh"
43#include "G4VDecayChannel.hh"
45#include "Randomize.hh"
46#include "G4LorentzVector.hh"
47#include "G4LorentzRotation.hh"
48
49
52{
53}
54
55
57 const G4String& theParentName,
58 G4double theBR,
59 const G4String& theLeptonName)
60 :G4VDecayChannel("Tau Leptonic Decay",1)
61{
62 // set names for daughter particles
63 if (theParentName == "tau+") {
64 SetBR(theBR);
65 SetParent("tau+");
67 if ((theLeptonName=="e-"||theLeptonName=="e+")){
68 SetDaughter(0, "e+");
69 SetDaughter(1, "nu_e");
70 SetDaughter(2, "anti_nu_tau");
71 } else {
72 SetDaughter(0, "mu+");
73 SetDaughter(1, "nu_mu");
74 SetDaughter(2, "anti_nu_tau");
75 }
76 } else if (theParentName == "tau-") {
77 SetBR(theBR);
78 SetParent("tau-");
80 if ((theLeptonName=="e-"||theLeptonName=="e+")){
81 SetDaughter(0, "e-");
82 SetDaughter(1, "anti_nu_e");
83 SetDaughter(2, "nu_tau");
84 } else {
85 SetDaughter(0, "mu-");
86 SetDaughter(1, "anti_nu_mu");
87 SetDaughter(2, "nu_tau");
88 }
89 } else {
90#ifdef G4VERBOSE
91 if (GetVerboseLevel()>0) {
92 G4cout << "G4TauLeptonicDecayChannel:: constructor :";
93 G4cout << " parent particle is not tau but ";
94 G4cout << theParentName << G4endl;
95 }
96#endif
97 }
98}
99
101{
102}
103
105 G4VDecayChannel(right)
106{
107}
108
110{
111 if (this != &right) {
114 rbranch = right.rbranch;
115
116 // copy parent name
117 parent_name = new G4String(*right.parent_name);
118
119 // clear daughters_name array
121
122 // recreate array
124 if ( numberOfDaughters >0 ) {
127 //copy daughters name
128 for (G4int index=0; index < numberOfDaughters; index++) {
129 daughters_name[index] = new G4String(*right.daughters_name[index]);
130 }
131 }
132 }
133 return *this;
134}
135
137{
138 // this version neglects muon polarization
139 // assumes the pure V-A coupling
140 // gives incorrect energy spectrum for neutrinos
141#ifdef G4VERBOSE
142 if (GetVerboseLevel()>1) G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
143#endif
144
145 if (parent == 0) FillParent();
146 if (daughters == 0) FillDaughters();
147
148 // parent mass
149 G4double parentmass = parent->GetPDGMass();
150
151 //daughters'mass
152 G4double daughtermass[3];
153 for (G4int index=0; index<3; index++){
154 daughtermass[index] = daughters[index]->GetPDGMass();
155 }
156
157 //create parent G4DynamicParticle at rest
158 G4ThreeVector dummy;
159 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
160 //create G4Decayproducts
161 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
162 delete parentparticle;
163
164 // calculate daughter momentum
165 G4double daughtermomentum[3];
166
167 // calcurate lepton momentum
168 G4double pmax = (parentmass*parentmass-daughtermass[0]*daughtermass[0])/2./parentmass;
169 G4double p, e;
170 G4double r;
171 do {
172 // determine momentum/energy
173 r = G4UniformRand();
174 p = pmax*G4UniformRand();
175 e = std::sqrt(p*p + daughtermass[0]*daughtermass[0]);
176 } while (r > spectrum(p,e,parentmass,daughtermass[0]) );
177
178 //create daughter G4DynamicParticle
179 // daughter 0 (lepton)
180 daughtermomentum[0] = p;
181 G4double costheta, sintheta, phi, sinphi, cosphi;
182 costheta = 2.*G4UniformRand()-1.0;
183 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
184 phi = twopi*G4UniformRand()*rad;
185 sinphi = std::sin(phi);
186 cosphi = std::cos(phi);
187 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
188 G4DynamicParticle * daughterparticle
189 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
190 products->PushProducts(daughterparticle);
191
192 // daughter 1 ,2 (nutrinos)
193 // create neutrinos in the C.M frame of two neutrinos
194 G4double energy2 = parentmass-e;
195 G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
196 G4double beta = -1.0*daughtermomentum[0]/energy2;
197 G4double costhetan = 2.*G4UniformRand()-1.0;
198 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
199 G4double phin = twopi*G4UniformRand()*rad;
200 G4double sinphin = std::sin(phin);
201 G4double cosphin = std::cos(phin);
202
203 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
204 G4DynamicParticle * daughterparticle1
205 = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
206 G4DynamicParticle * daughterparticle2
207 = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
208
209 // boost to the muon rest frame
211 p4 = daughterparticle1->Get4Momentum();
212 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
213 daughterparticle1->Set4Momentum(p4);
214 p4 = daughterparticle2->Get4Momentum();
215 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
216 daughterparticle2->Set4Momentum(p4);
217 products->PushProducts(daughterparticle1);
218 products->PushProducts(daughterparticle2);
219 daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
220 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
221
222
223 // output message
224#ifdef G4VERBOSE
225 if (GetVerboseLevel()>1) {
226 G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
227 G4cout << " create decay products in rest frame " <<G4endl;
228 products->DumpInfo();
229 }
230#endif
231 return products;
232}
233
234
235
236
237G4double G4TauLeptonicDecayChannel::spectrum(G4double p,
238 G4double e,
239 G4double mtau,
240 G4double ml)
241{
242 G4double f1;
243 f1 = 3.0*e*(mtau*mtau+ml*ml)-4.0*mtau*e*e-2.0*mtau*ml*ml;
244 return p*(f1)/(mtau*mtau*mtau*mtau)/(0.6);
245}
246
247
248
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
double z() const
double x() const
double y() const
HepLorentzVector & boost(double, double, double)
void DumpInfo() const
G4int PushProducts(G4DynamicParticle *aParticle)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
G4double GetTotalMomentum() const
G4TauLeptonicDecayChannel & operator=(const G4TauLeptonicDecayChannel &)
virtual G4DecayProducts * DecayIt(G4double)
G4String * parent_name
void SetBR(G4double value)
G4String ** daughters_name
G4int GetVerboseLevel() const
void SetNumberOfDaughters(G4int value)
void SetDaughter(G4int anIndex, const G4ParticleDefinition *particle_type)
G4ParticleDefinition * parent
G4ParticleDefinition ** daughters
G4String kinematics_name
void SetParent(const G4ParticleDefinition *particle_type)