| Geant4 Cross Reference |
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 // G4TauLeptonicDecayChannel class implementation
27 //
28 // Author: H.Kurashige, 30 May 1997
29 // --------------------------------------------------------------------
30
31 #include "G4TauLeptonicDecayChannel.hh"
32
33 #include "G4DecayProducts.hh"
34 #include "G4LorentzRotation.hh"
35 #include "G4LorentzVector.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4VDecayChannel.hh"
40 #include "Randomize.hh"
41
42 G4TauLeptonicDecayChannel::G4TauLeptonicDecayChannel(const G4String& theParentName, G4double theBR,
43 const G4String& theLeptonName)
44 : G4VDecayChannel("Tau Leptonic Decay", 1)
45 {
46 // set names for daughter particles
47 if (theParentName == "tau+") {
48 SetBR(theBR);
49 SetParent("tau+");
50 SetNumberOfDaughters(3);
51 if ((theLeptonName == "e-" || theLeptonName == "e+")) {
52 SetDaughter(0, "e+");
53 SetDaughter(1, "nu_e");
54 SetDaughter(2, "anti_nu_tau");
55 }
56 else {
57 SetDaughter(0, "mu+");
58 SetDaughter(1, "nu_mu");
59 SetDaughter(2, "anti_nu_tau");
60 }
61 }
62 else if (theParentName == "tau-") {
63 SetBR(theBR);
64 SetParent("tau-");
65 SetNumberOfDaughters(3);
66 if ((theLeptonName == "e-" || theLeptonName == "e+")) {
67 SetDaughter(0, "e-");
68 SetDaughter(1, "anti_nu_e");
69 SetDaughter(2, "nu_tau");
70 }
71 else {
72 SetDaughter(0, "mu-");
73 SetDaughter(1, "anti_nu_mu");
74 SetDaughter(2, "nu_tau");
75 }
76 }
77 else {
78 #ifdef G4VERBOSE
79 if (GetVerboseLevel() > 0) {
80 G4cout << "G4TauLeptonicDecayChannel:: constructor :";
81 G4cout << " parent particle is not tau but ";
82 G4cout << theParentName << G4endl;
83 }
84 #endif
85 }
86 }
87
88 G4TauLeptonicDecayChannel&
89 G4TauLeptonicDecayChannel::operator=(const G4TauLeptonicDecayChannel& right)
90 {
91 if (this != &right) {
92 kinematics_name = right.kinematics_name;
93 verboseLevel = right.verboseLevel;
94 rbranch = right.rbranch;
95
96 // copy parent name
97 parent_name = new G4String(*right.parent_name);
98
99 // clear daughters_name array
100 ClearDaughtersName();
101
102 // recreate array
103 numberOfDaughters = right.numberOfDaughters;
104 if (numberOfDaughters > 0) {
105 if (daughters_name != nullptr) ClearDaughtersName();
106 daughters_name = new G4String*[numberOfDaughters];
107 // copy daughters name
108 for (G4int index = 0; index < numberOfDaughters; ++index) {
109 daughters_name[index] = new G4String(*right.daughters_name[index]);
110 }
111 }
112 }
113 return *this;
114 }
115
116 G4DecayProducts* G4TauLeptonicDecayChannel::DecayIt(G4double)
117 {
118 // this version neglects muon polarization
119 // assumes the pure V-A coupling
120 // gives incorrect energy spectrum for neutrinos
121
122 #ifdef G4VERBOSE
123 if (GetVerboseLevel() > 1) G4cout << "G4TauLeptonicDecayChannel::DecayIt()";
124 #endif
125
126 CheckAndFillParent();
127 CheckAndFillDaughters();
128
129 // parent mass
130 G4double parentmass = G4MT_parent->GetPDGMass();
131
132 // daughters'mass
133 const G4int N_DAUGHTER = 3;
134 G4double daughtermass[N_DAUGHTER];
135 for (G4int index = 0; index < N_DAUGHTER; ++index) {
136 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
137 }
138
139 // create parent G4DynamicParticle at rest
140 G4ThreeVector dummy;
141 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
142 // create G4Decayproducts
143 auto products = new G4DecayProducts(*parentparticle);
144 delete parentparticle;
145
146 // calculate daughter momentum
147 G4double daughtermomentum[N_DAUGHTER];
148
149 // calculate lepton momentum
150 G4double pmax = (parentmass * parentmass - daughtermass[0] * daughtermass[0]) / 2. / parentmass;
151 G4double p, e;
152 G4double r;
153 const std::size_t MAX_LOOP = 10000;
154 for (std::size_t loop_counter = 0; loop_counter < MAX_LOOP; ++loop_counter) {
155 // determine momentum/energy
156 r = G4UniformRand();
157 p = pmax * G4UniformRand();
158 e = std::sqrt(p * p + daughtermass[0] * daughtermass[0]);
159 if (r < spectrum(p, e, parentmass, daughtermass[0])) break;
160 }
161
162 // create daughter G4DynamicParticle
163 // daughter 0 (lepton)
164 daughtermomentum[0] = p;
165 G4double costheta, sintheta, phi, sinphi, cosphi;
166 costheta = 2. * G4UniformRand() - 1.0;
167 sintheta = std::sqrt((1.0 - costheta) * (1.0 + costheta));
168 phi = twopi * G4UniformRand() * rad;
169 sinphi = std::sin(phi);
170 cosphi = std::cos(phi);
171 G4ThreeVector direction0(sintheta * cosphi, sintheta * sinphi, costheta);
172 auto daughterparticle =
173 new G4DynamicParticle(G4MT_daughters[0], direction0 * daughtermomentum[0]);
174 products->PushProducts(daughterparticle);
175
176 // daughter 1 ,2 (nutrinos)
177 // create neutrinos in the C.M frame of two neutrinos
178 G4double energy2 = parentmass - e;
179 G4double vmass = std::sqrt((energy2 - daughtermomentum[0]) * (energy2 + daughtermomentum[0]));
180 G4double beta = -1.0 * daughtermomentum[0] / energy2;
181 G4double costhetan = 2. * G4UniformRand() - 1.0;
182 G4double sinthetan = std::sqrt((1.0 - costhetan) * (1.0 + costhetan));
183 G4double phin = twopi * G4UniformRand() * rad;
184 G4double sinphin = std::sin(phin);
185 G4double cosphin = std::cos(phin);
186
187 G4ThreeVector direction1(sinthetan * cosphin, sinthetan * sinphin, costhetan);
188 auto daughterparticle1 = new G4DynamicParticle(G4MT_daughters[1], direction1 * (vmass / 2.));
189 auto daughterparticle2 =
190 new G4DynamicParticle(G4MT_daughters[2], direction1 * (-1.0 * vmass / 2.));
191
192 // boost to the muon rest frame
193 G4LorentzVector p4;
194 p4 = daughterparticle1->Get4Momentum();
195 p4.boost(direction0.x() * beta, direction0.y() * beta, direction0.z() * beta);
196 daughterparticle1->Set4Momentum(p4);
197 p4 = daughterparticle2->Get4Momentum();
198 p4.boost(direction0.x() * beta, direction0.y() * beta, direction0.z() * beta);
199 daughterparticle2->Set4Momentum(p4);
200 products->PushProducts(daughterparticle1);
201 products->PushProducts(daughterparticle2);
202 daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
203 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
204
205 // output message
206 #ifdef G4VERBOSE
207 if (GetVerboseLevel() > 1) {
208 G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
209 G4cout << " create decay products in rest frame " << G4endl;
210 products->DumpInfo();
211 }
212 #endif
213 return products;
214 }
215
216 G4double G4TauLeptonicDecayChannel::spectrum(G4double p, G4double e, G4double mtau, G4double ml)
217 {
218 G4double f1;
219 f1 = 3.0 * e * (mtau * mtau + ml * ml) - 4.0 * mtau * e * e - 2.0 * mtau * ml * ml;
220 return p * (f1) / (mtau * mtau * mtau * mtau) / (0.6);
221 }
222