| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4TauLeptonicDecayChannel class implementat <<
27 // 23 //
28 // Author: H.Kurashige, 30 May 1997 << 24 // $Id: G4TauLeptonicDecayChannel.cc,v 1.4 2004/12/10 18:02:04 gcosmo Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-08-00-patch-01 $
30 << 26 //
31 #include "G4TauLeptonicDecayChannel.hh" << 27 //
>> 28 // ------------------------------------------------------------
>> 29 // GEANT 4 class header file
>> 30 //
>> 31 // History: first implementation, based on object model of
>> 32 // 30 May 1997 H.Kurashige
>> 33 //
>> 34 // Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige
>> 35 // ------------------------------------------------------------
32 36
33 #include "G4DecayProducts.hh" <<
34 #include "G4LorentzRotation.hh" <<
35 #include "G4LorentzVector.hh" <<
36 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh"
37 #include "G4PhysicalConstants.hh" << 38 #include "G4DecayProducts.hh"
38 #include "G4SystemOfUnits.hh" <<
39 #include "G4VDecayChannel.hh" 39 #include "G4VDecayChannel.hh"
>> 40 #include "G4TauLeptonicDecayChannel.hh"
40 #include "Randomize.hh" 41 #include "Randomize.hh"
>> 42 #include "G4LorentzVector.hh"
>> 43 #include "G4LorentzRotation.hh"
41 44
42 G4TauLeptonicDecayChannel::G4TauLeptonicDecayC << 45
43 << 46 G4TauLeptonicDecayChannel::G4TauLeptonicDecayChannel(
44 : G4VDecayChannel("Tau Leptonic Decay", 1) << 47 const G4String& theParentName,
>> 48 G4double theBR,
>> 49 const G4String& theLeptonName)
>> 50 :G4VDecayChannel("Tau Leptonic Decay",1)
45 { 51 {
46 // set names for daughter particles 52 // set names for daughter particles
47 if (theParentName == "tau+") { 53 if (theParentName == "tau+") {
48 SetBR(theBR); 54 SetBR(theBR);
49 SetParent("tau+"); 55 SetParent("tau+");
50 SetNumberOfDaughters(3); 56 SetNumberOfDaughters(3);
51 if ((theLeptonName == "e-" || theLeptonNam << 57 if ((theLeptonName=="e-"||theLeptonName=="e+")){
52 SetDaughter(0, "e+"); 58 SetDaughter(0, "e+");
53 SetDaughter(1, "nu_e"); 59 SetDaughter(1, "nu_e");
54 SetDaughter(2, "anti_nu_tau"); 60 SetDaughter(2, "anti_nu_tau");
55 } << 61 } else {
56 else { <<
57 SetDaughter(0, "mu+"); 62 SetDaughter(0, "mu+");
58 SetDaughter(1, "nu_mu"); 63 SetDaughter(1, "nu_mu");
59 SetDaughter(2, "anti_nu_tau"); 64 SetDaughter(2, "anti_nu_tau");
60 } << 65 }
61 } << 66 } else if (theParentName == "tau-") {
62 else if (theParentName == "tau-") { <<
63 SetBR(theBR); 67 SetBR(theBR);
64 SetParent("tau-"); 68 SetParent("tau-");
65 SetNumberOfDaughters(3); 69 SetNumberOfDaughters(3);
66 if ((theLeptonName == "e-" || theLeptonNam << 70 if ((theLeptonName=="e-"||theLeptonName=="e+")){
67 SetDaughter(0, "e-"); 71 SetDaughter(0, "e-");
68 SetDaughter(1, "anti_nu_e"); 72 SetDaughter(1, "anti_nu_e");
69 SetDaughter(2, "nu_tau"); 73 SetDaughter(2, "nu_tau");
70 } << 74 } else {
71 else { <<
72 SetDaughter(0, "mu-"); 75 SetDaughter(0, "mu-");
73 SetDaughter(1, "anti_nu_mu"); 76 SetDaughter(1, "anti_nu_mu");
74 SetDaughter(2, "nu_tau"); 77 SetDaughter(2, "nu_tau");
75 } << 78 }
76 } << 79 } else {
77 else { <<
78 #ifdef G4VERBOSE 80 #ifdef G4VERBOSE
79 if (GetVerboseLevel() > 0) { << 81 if (GetVerboseLevel()>0) {
80 G4cout << "G4TauLeptonicDecayChannel:: c 82 G4cout << "G4TauLeptonicDecayChannel:: constructor :";
81 G4cout << " parent particle is not tau b 83 G4cout << " parent particle is not tau but ";
82 G4cout << theParentName << G4endl; 84 G4cout << theParentName << G4endl;
83 } 85 }
84 #endif 86 #endif
85 } 87 }
86 } 88 }
87 89
88 G4TauLeptonicDecayChannel& << 90 G4TauLeptonicDecayChannel::~G4TauLeptonicDecayChannel()
89 G4TauLeptonicDecayChannel::operator=(const G4T <<
90 { 91 {
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_n <<
98 <<
99 // clear daughters_name array <<
100 ClearDaughtersName(); <<
101 <<
102 // recreate array <<
103 numberOfDaughters = right.numberOfDaughter <<
104 if (numberOfDaughters > 0) { <<
105 if (daughters_name != nullptr) ClearDaug <<
106 daughters_name = new G4String*[numberOfD <<
107 // copy daughters name <<
108 for (G4int index = 0; index < numberOfDa <<
109 daughters_name[index] = new G4String(* <<
110 } <<
111 } <<
112 } <<
113 return *this; <<
114 } 92 }
115 93
116 G4DecayProducts* G4TauLeptonicDecayChannel::De << 94 G4DecayProducts *G4TauLeptonicDecayChannel::DecayIt(G4double)
117 { 95 {
118 // this version neglects muon polarization << 96 // this version neglects muon polarization
119 // assumes the pure V-A couplin 97 // assumes the pure V-A coupling
120 // gives incorrect energy spect 98 // gives incorrect energy spectrum for neutrinos
121 <<
122 #ifdef G4VERBOSE 99 #ifdef G4VERBOSE
123 if (GetVerboseLevel() > 1) G4cout << "G4TauL << 100 if (GetVerboseLevel()>1) G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
124 #endif 101 #endif
125 102
126 CheckAndFillParent(); << 103 if (parent == 0) FillParent();
127 CheckAndFillDaughters(); << 104 if (daughters == 0) FillDaughters();
128 << 105
129 // parent mass 106 // parent mass
130 G4double parentmass = G4MT_parent->GetPDGMas << 107 G4double parentmass = parent->GetPDGMass();
131 108
132 // daughters'mass << 109 //daughters'mass
133 const G4int N_DAUGHTER = 3; << 110 G4double daughtermass[3];
134 G4double daughtermass[N_DAUGHTER]; << 111 for (G4int index=0; index<3; index++){
135 for (G4int index = 0; index < N_DAUGHTER; ++ << 112 daughtermass[index] = daughters[index]->GetPDGMass();
136 daughtermass[index] = G4MT_daughters[index <<
137 } 113 }
138 114
139 // create parent G4DynamicParticle at rest << 115 //create parent G4DynamicParticle at rest
140 G4ThreeVector dummy; 116 G4ThreeVector dummy;
141 auto parentparticle = new G4DynamicParticle( << 117 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
142 // create G4Decayproducts << 118 //create G4Decayproducts
143 auto products = new G4DecayProducts(*parentp << 119 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
144 delete parentparticle; 120 delete parentparticle;
145 121
146 // calculate daughter momentum 122 // calculate daughter momentum
147 G4double daughtermomentum[N_DAUGHTER]; << 123 G4double daughtermomentum[3];
>> 124 G4double energy;
148 125
149 // calculate lepton momentum << 126 // calcurate lepton momentum
150 G4double pmax = (parentmass * parentmass - d << 127 G4double pmax = (parentmass*parentmass-daughtermass[0]*daughtermass[0])/2./parentmass;
151 G4double p, e; 128 G4double p, e;
152 G4double r; 129 G4double r;
153 const std::size_t MAX_LOOP = 10000; << 130 do {
154 for (std::size_t loop_counter = 0; loop_coun <<
155 // determine momentum/energy 131 // determine momentum/energy
156 r = G4UniformRand(); 132 r = G4UniformRand();
157 p = pmax * G4UniformRand(); << 133 p = pmax*G4UniformRand();
158 e = std::sqrt(p * p + daughtermass[0] * da << 134 e = std::sqrt(p*p + daughtermass[0]*daughtermass[0]);
159 if (r < spectrum(p, e, parentmass, daughte << 135 } while (r > spectrum(p,e,parentmass,daughtermass[0]) );
160 } << 136 energy = e- daughtermass[0];
161 137
162 // create daughter G4DynamicParticle << 138 //create daughter G4DynamicParticle
163 // daughter 0 (lepton) 139 // daughter 0 (lepton)
164 daughtermomentum[0] = p; 140 daughtermomentum[0] = p;
165 G4double costheta, sintheta, phi, sinphi, co << 141 G4double costheta, sintheta, phi, sinphi, cosphi;
166 costheta = 2. * G4UniformRand() - 1.0; << 142 costheta = 2.*G4UniformRand()-1.0;
167 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 143 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
168 phi = twopi * G4UniformRand() * rad; << 144 phi = twopi*G4UniformRand()*rad;
169 sinphi = std::sin(phi); 145 sinphi = std::sin(phi);
170 cosphi = std::cos(phi); 146 cosphi = std::cos(phi);
171 G4ThreeVector direction0(sintheta * cosphi, << 147 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
172 auto daughterparticle = << 148 G4DynamicParticle * daughterparticle
173 new G4DynamicParticle(G4MT_daughters[0], d << 149 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
174 products->PushProducts(daughterparticle); 150 products->PushProducts(daughterparticle);
175 151
176 // daughter 1 ,2 (nutrinos) 152 // daughter 1 ,2 (nutrinos)
177 // create neutrinos in the C.M frame of two 153 // create neutrinos in the C.M frame of two neutrinos
178 G4double energy2 = parentmass - e; << 154 G4double energy2 = parentmass-e;
179 G4double vmass = std::sqrt((energy2 - daught << 155 G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
180 G4double beta = -1.0 * daughtermomentum[0] / << 156 G4double beta = -1.0*daughtermomentum[0]/energy2;
181 G4double costhetan = 2. * G4UniformRand() - << 157 G4double costhetan = 2.*G4UniformRand()-1.0;
182 G4double sinthetan = std::sqrt((1.0 - costhe << 158 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
183 G4double phin = twopi * G4UniformRand() * ra << 159 G4double phin = twopi*G4UniformRand()*rad;
184 G4double sinphin = std::sin(phin); 160 G4double sinphin = std::sin(phin);
185 G4double cosphin = std::cos(phin); 161 G4double cosphin = std::cos(phin);
186 162
187 G4ThreeVector direction1(sinthetan * cosphin << 163 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
188 auto daughterparticle1 = new G4DynamicPartic << 164 G4DynamicParticle * daughterparticle1
189 auto daughterparticle2 = << 165 = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
190 new G4DynamicParticle(G4MT_daughters[2], d << 166 G4DynamicParticle * daughterparticle2
>> 167 = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
191 168
192 // boost to the muon rest frame 169 // boost to the muon rest frame
193 G4LorentzVector p4; 170 G4LorentzVector p4;
194 p4 = daughterparticle1->Get4Momentum(); 171 p4 = daughterparticle1->Get4Momentum();
195 p4.boost(direction0.x() * beta, direction0.y << 172 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
196 daughterparticle1->Set4Momentum(p4); 173 daughterparticle1->Set4Momentum(p4);
197 p4 = daughterparticle2->Get4Momentum(); 174 p4 = daughterparticle2->Get4Momentum();
198 p4.boost(direction0.x() * beta, direction0.y << 175 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
199 daughterparticle2->Set4Momentum(p4); 176 daughterparticle2->Set4Momentum(p4);
200 products->PushProducts(daughterparticle1); 177 products->PushProducts(daughterparticle1);
201 products->PushProducts(daughterparticle2); 178 products->PushProducts(daughterparticle2);
202 daughtermomentum[1] = daughterparticle1->Get 179 daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
203 daughtermomentum[2] = daughterparticle2->Get 180 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
204 181
205 // output message << 182
>> 183 // output message
206 #ifdef G4VERBOSE 184 #ifdef G4VERBOSE
207 if (GetVerboseLevel() > 1) { << 185 if (GetVerboseLevel()>1) {
208 G4cout << "G4TauLeptonicDecayChannel::Deca 186 G4cout << "G4TauLeptonicDecayChannel::DecayIt ";
209 G4cout << " create decay products in rest << 187 G4cout << " create decay products in rest frame " <<G4endl;
210 products->DumpInfo(); 188 products->DumpInfo();
211 } 189 }
212 #endif 190 #endif
213 return products; 191 return products;
214 } 192 }
215 193
216 G4double G4TauLeptonicDecayChannel::spectrum(G << 194
>> 195
>> 196
>> 197 G4double G4TauLeptonicDecayChannel::spectrum(G4double p,
>> 198 G4double e,
>> 199 G4double mtau,
>> 200 G4double ml)
217 { 201 {
218 G4double f1; 202 G4double f1;
219 f1 = 3.0 * e * (mtau * mtau + ml * ml) - 4.0 << 203 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 << 204 return p*(f1)/(mtau*mtau*mtau*mtau)/(0.6);
221 } 205 }
>> 206
>> 207
>> 208
222 209