| 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 // G4MuonDecayChannel class implementation <<
27 // 23 //
28 // Author: H.Kurashige, 30 May 1997 << 24 // $Id: G4MuonDecayChannel.cc,v 1.9 2001/07/11 10:02:00 gunter Exp $
29 // Contributions: << 25 // GEANT4 tag $Name: geant4-05-00 $
30 // - 2005 - M.Melissas, J.Brunner CPPM/IN2P3 << 26 //
31 // Added V-A fluxes for neutrinos using a ne << 27 //
32 // ------------------------------------------- << 28 // ------------------------------------------------------------
33 << 29 // GEANT 4 class header file
34 #include "G4MuonDecayChannel.hh" << 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 // ------------------------------------------------------------
35 36
36 #include "G4DecayProducts.hh" <<
37 #include "G4LorentzRotation.hh" <<
38 #include "G4LorentzVector.hh" <<
39 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh"
40 #include "G4PhysicalConstants.hh" << 38 #include "G4DecayProducts.hh"
41 #include "G4RotationMatrix.hh" <<
42 #include "G4SystemOfUnits.hh" <<
43 #include "G4VDecayChannel.hh" 39 #include "G4VDecayChannel.hh"
>> 40 #include "G4MuonDecayChannel.hh"
44 #include "Randomize.hh" 41 #include "Randomize.hh"
>> 42 #include "G4LorentzVector.hh"
>> 43 #include "G4LorentzRotation.hh"
45 44
46 G4MuonDecayChannel::G4MuonDecayChannel(const G << 45
47 : G4VDecayChannel("Muon Decay", 1) << 46 G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName,
>> 47 G4double theBR)
>> 48 :G4VDecayChannel("Muon Decay",1)
48 { 49 {
49 // set names for daughter particles 50 // set names for daughter particles
50 if (theParentName == "mu+") { 51 if (theParentName == "mu+") {
51 SetBR(theBR); 52 SetBR(theBR);
52 SetParent("mu+"); 53 SetParent("mu+");
53 SetNumberOfDaughters(3); 54 SetNumberOfDaughters(3);
54 SetDaughter(0, "e+"); 55 SetDaughter(0, "e+");
55 SetDaughter(1, "nu_e"); 56 SetDaughter(1, "nu_e");
56 SetDaughter(2, "anti_nu_mu"); 57 SetDaughter(2, "anti_nu_mu");
57 } << 58 } else if (theParentName == "mu-") {
58 else if (theParentName == "mu-") { <<
59 SetBR(theBR); 59 SetBR(theBR);
60 SetParent("mu-"); 60 SetParent("mu-");
61 SetNumberOfDaughters(3); 61 SetNumberOfDaughters(3);
62 SetDaughter(0, "e-"); 62 SetDaughter(0, "e-");
63 SetDaughter(1, "anti_nu_e"); 63 SetDaughter(1, "anti_nu_e");
64 SetDaughter(2, "nu_mu"); 64 SetDaughter(2, "nu_mu");
65 } << 65 } else {
66 else { <<
67 #ifdef G4VERBOSE 66 #ifdef G4VERBOSE
68 if (GetVerboseLevel() > 0) { << 67 if (GetVerboseLevel()>0) {
69 G4cout << "G4MuonDecayChannel:: construc 68 G4cout << "G4MuonDecayChannel:: constructor :";
70 G4cout << " parent particle is not muon 69 G4cout << " parent particle is not muon but ";
71 G4cout << theParentName << G4endl; 70 G4cout << theParentName << G4endl;
72 } 71 }
73 #endif 72 #endif
74 } 73 }
75 } 74 }
76 75
77 G4MuonDecayChannel& G4MuonDecayChannel::operat << 76 G4MuonDecayChannel::~G4MuonDecayChannel()
78 { 77 {
79 if (this != &right) { <<
80 kinematics_name = right.kinematics_name; <<
81 verboseLevel = right.verboseLevel; <<
82 rbranch = right.rbranch; <<
83 <<
84 // copy parent name <<
85 parent_name = new G4String(*right.parent_n <<
86 <<
87 // clear daughters_name array <<
88 ClearDaughtersName(); <<
89 <<
90 // recreate array <<
91 numberOfDaughters = right.numberOfDaughter <<
92 if (numberOfDaughters > 0) { <<
93 if (daughters_name != nullptr) ClearDaug <<
94 daughters_name = new G4String*[numberOfD <<
95 // copy daughters name <<
96 for (G4int index = 0; index < numberOfDa <<
97 daughters_name[index] = new G4String(* <<
98 } <<
99 } <<
100 } <<
101 return *this; <<
102 } 78 }
103 79
104 G4DecayProducts* G4MuonDecayChannel::DecayIt(G << 80 G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double)
105 { 81 {
106 // this version neglects muon polarization,a << 82 // this version neglects muon polarization
107 // assumes the pure V-A couplin 83 // assumes the pure V-A coupling
108 // the Neutrinos are correctly << 84 // gives incorrect energy spectrum for neutrinos
109 <<
110 #ifdef G4VERBOSE 85 #ifdef G4VERBOSE
111 if (GetVerboseLevel() > 1) G4cout << "G4Muon << 86 if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt ";
112 #endif 87 #endif
113 88
114 CheckAndFillParent(); << 89 if (parent == 0) FillParent();
115 CheckAndFillDaughters(); << 90 if (daughters == 0) FillDaughters();
116 << 91
117 // parent mass 92 // parent mass
118 G4double parentmass = G4MT_parent->GetPDGMas << 93 G4double parentmass = parent->GetPDGMass();
119 const G4int N_DAUGHTER = 3; <<
120 94
121 // daughters'mass << 95 //daughters'mass
122 G4double daughtermass[N_DAUGHTER]; << 96 G4double daughtermass[3];
123 // G4double sumofdaughtermass = 0.0; << 97 G4double sumofdaughtermass = 0.0;
124 for (G4int index = 0; index < N_DAUGHTER; ++ << 98 for (G4int index=0; index<3; index++){
125 daughtermass[index] = G4MT_daughters[index << 99 daughtermass[index] = daughters[index]->GetPDGMass();
126 // sumofdaughtermass += daughtermass[index << 100 sumofdaughtermass += daughtermass[index];
127 } 101 }
128 102
129 // create parent G4DynamicParticle at rest << 103 //create parent G4DynamicParticle at rest
130 G4ThreeVector dummy; 104 G4ThreeVector dummy;
131 auto parentparticle = new G4DynamicParticle( << 105 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
132 // create G4Decayproducts << 106 //create G4Decayproducts
133 auto products = new G4DecayProducts(*parentp << 107 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
134 delete parentparticle; 108 delete parentparticle;
135 109
136 // calculate daughter momentum 110 // calculate daughter momentum
137 G4double daughtermomentum[N_DAUGHTER]; << 111 G4double daughtermomentum[3];
138 // calculate electron energy << 112 G4double energy;
139 G4double xmax = (1.0 + daughtermass[0] * dau << 113 // calcurate electron energy
>> 114 G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass);
140 G4double x; 115 G4double x;
141 << 116 G4double r;
142 G4double Ee, Ene; << 117 do {
143 << 118 do {
144 G4double gam; << 119 r = G4UniformRand();
145 G4double EMax = parentmass / 2 - daughtermas << 120 x = xmax*G4UniformRand();
146 << 121 } while (r > (3.0 - 2.0*x)*x*x);
147 const std::size_t MAX_LOOP = 1000; << 122 energy = x*parentmass/2.0 - daughtermass[0];
148 // Generating Random Energy << 123 } while (energy <0.0);
149 for (std::size_t loop1 = 0; loop1 < MAX_LOOP << 124 //create daughter G4DynamicParticle
150 Ee = G4UniformRand(); << 125 // daughter 0 (electron)
151 for (std::size_t loop2 = 0; loop2 < MAX_LO << 126 daughtermomentum[0] = sqrt(energy*energy + 2.0*energy* daughtermass[0]);
152 x = xmax * G4UniformRand(); << 127 G4double costheta, sintheta, phi, sinphi, cosphi;
153 gam = G4UniformRand(); << 128 costheta = 2.*G4UniformRand()-1.0;
154 if (gam <= x * (1. - x)) break; << 129 sintheta = sqrt((1.0-costheta)*(1.0+costheta));
155 x = xmax; << 130 phi = 2.0*M_PI*G4UniformRand()*rad;
156 } << 131 sinphi = sin(phi);
157 Ene = x; << 132 cosphi = cos(phi);
158 if (Ene >= (1. - Ee)) break; << 133 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
159 Ene = 1. - Ee; << 134 G4DynamicParticle * daughterparticle
160 } << 135 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
161 G4double Enm = (2. - Ee - Ene); <<
162 <<
163 // initialisation of rotation parameters <<
164 <<
165 G4double costheta, sintheta, rphi, rtheta, r <<
166 costheta = 1. - 2. / Ee - 2. / Ene + 2. / En <<
167 sintheta = std::sqrt(1. - costheta * costhet <<
168 <<
169 rphi = twopi * G4UniformRand() * rad; <<
170 rtheta = (std::acos(2. * G4UniformRand() - 1 <<
171 rpsi = twopi * G4UniformRand() * rad; <<
172 <<
173 G4RotationMatrix rot; <<
174 rot.set(rphi, rtheta, rpsi); <<
175 <<
176 // electron 0 <<
177 daughtermomentum[0] = std::sqrt(Ee * Ee * EM <<
178 G4ThreeVector direction0(0.0, 0.0, 1.0); <<
179 <<
180 direction0 *= rot; <<
181 <<
182 auto daughterparticle = <<
183 new G4DynamicParticle(G4MT_daughters[0], d <<
184 <<
185 products->PushProducts(daughterparticle); 136 products->PushProducts(daughterparticle);
186 137
187 // electronic neutrino 1 << 138 // daughter 1 ,2 (nutrinos)
188 << 139 // create neutrinos in the C.M frame of two neutrinos
189 daughtermomentum[1] = std::sqrt(Ene * Ene * << 140 G4double energy2 = parentmass*(1.0 - x/2.0);
190 G4ThreeVector direction1(sintheta, 0.0, cost << 141 G4double vmass = sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
191 << 142 G4double beta = -1.0*daughtermomentum[0]/energy2;
192 direction1 *= rot; << 143 G4double costhetan = 2.*G4UniformRand()-1.0;
193 << 144 G4double sinthetan = sqrt((1.0-costhetan)*(1.0+costhetan));
194 auto daughterparticle1 = << 145 G4double phin = 2.0*M_PI*G4UniformRand()*rad;
195 new G4DynamicParticle(G4MT_daughters[1], d << 146 G4double sinphin = sin(phin);
>> 147 G4double cosphin = cos(phin);
>> 148
>> 149 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
>> 150 G4DynamicParticle * daughterparticle1
>> 151 = new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
>> 152 G4DynamicParticle * daughterparticle2
>> 153 = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
>> 154
>> 155 // boost to the muon rest frame
>> 156 G4LorentzVector p4;
>> 157 p4 = daughterparticle1->Get4Momentum();
>> 158 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
>> 159 daughterparticle1->Set4Momentum(p4);
>> 160 p4 = daughterparticle2->Get4Momentum();
>> 161 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
>> 162 daughterparticle2->Set4Momentum(p4);
196 products->PushProducts(daughterparticle1); 163 products->PushProducts(daughterparticle1);
197 <<
198 // muonnic neutrino 2 <<
199 <<
200 daughtermomentum[2] = std::sqrt(Enm * Enm * <<
201 G4ThreeVector direction2(-Ene / Enm * sinthe <<
202 <<
203 direction2 *= rot; <<
204 <<
205 auto daughterparticle2 = <<
206 new G4DynamicParticle(G4MT_daughters[2], d <<
207 products->PushProducts(daughterparticle2); 164 products->PushProducts(daughterparticle2);
>> 165 daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
>> 166 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
208 167
209 // output message << 168 // output message
210 #ifdef G4VERBOSE 169 #ifdef G4VERBOSE
211 if (GetVerboseLevel() > 1) { << 170 if (GetVerboseLevel()>1) {
212 G4cout << "G4MuonDecayChannel::DecayIt()"; << 171 G4cout << "G4MuonDecayChannel::DecayIt ";
213 G4cout << " create decay products in rest << 172 G4cout << " create decay products in rest frame " <<G4endl;
214 products->DumpInfo(); 173 products->DumpInfo();
215 } 174 }
216 #endif 175 #endif
217 return products; 176 return products;
218 } 177 }
>> 178
>> 179
>> 180
>> 181
>> 182
>> 183
219 184