| 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 // G4MuonDecayChannel class implementation
27 //
28 // Author: H.Kurashige, 30 May 1997
29 // Contributions:
30 // - 2005 - M.Melissas, J.Brunner CPPM/IN2P3
31 // Added V-A fluxes for neutrinos using a new algorithm, 2005
32 // --------------------------------------------------------------------
33
34 #include "G4MuonDecayChannel.hh"
35
36 #include "G4DecayProducts.hh"
37 #include "G4LorentzRotation.hh"
38 #include "G4LorentzVector.hh"
39 #include "G4ParticleDefinition.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4RotationMatrix.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4VDecayChannel.hh"
44 #include "Randomize.hh"
45
46 G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName, G4double theBR)
47 : G4VDecayChannel("Muon Decay", 1)
48 {
49 // set names for daughter particles
50 if (theParentName == "mu+") {
51 SetBR(theBR);
52 SetParent("mu+");
53 SetNumberOfDaughters(3);
54 SetDaughter(0, "e+");
55 SetDaughter(1, "nu_e");
56 SetDaughter(2, "anti_nu_mu");
57 }
58 else if (theParentName == "mu-") {
59 SetBR(theBR);
60 SetParent("mu-");
61 SetNumberOfDaughters(3);
62 SetDaughter(0, "e-");
63 SetDaughter(1, "anti_nu_e");
64 SetDaughter(2, "nu_mu");
65 }
66 else {
67 #ifdef G4VERBOSE
68 if (GetVerboseLevel() > 0) {
69 G4cout << "G4MuonDecayChannel:: constructor :";
70 G4cout << " parent particle is not muon but ";
71 G4cout << theParentName << G4endl;
72 }
73 #endif
74 }
75 }
76
77 G4MuonDecayChannel& G4MuonDecayChannel::operator=(const G4MuonDecayChannel& right)
78 {
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_name);
86
87 // clear daughters_name array
88 ClearDaughtersName();
89
90 // recreate array
91 numberOfDaughters = right.numberOfDaughters;
92 if (numberOfDaughters > 0) {
93 if (daughters_name != nullptr) ClearDaughtersName();
94 daughters_name = new G4String*[numberOfDaughters];
95 // copy daughters name
96 for (G4int index = 0; index < numberOfDaughters; ++index) {
97 daughters_name[index] = new G4String(*right.daughters_name[index]);
98 }
99 }
100 }
101 return *this;
102 }
103
104 G4DecayProducts* G4MuonDecayChannel::DecayIt(G4double)
105 {
106 // this version neglects muon polarization,and electron mass
107 // assumes the pure V-A coupling
108 // the Neutrinos are correctly V-A
109
110 #ifdef G4VERBOSE
111 if (GetVerboseLevel() > 1) G4cout << "G4MuonDecayChannel::DecayIt ";
112 #endif
113
114 CheckAndFillParent();
115 CheckAndFillDaughters();
116
117 // parent mass
118 G4double parentmass = G4MT_parent->GetPDGMass();
119 const G4int N_DAUGHTER = 3;
120
121 // daughters'mass
122 G4double daughtermass[N_DAUGHTER];
123 // G4double sumofdaughtermass = 0.0;
124 for (G4int index = 0; index < N_DAUGHTER; ++index) {
125 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
126 // sumofdaughtermass += daughtermass[index];
127 }
128
129 // create parent G4DynamicParticle at rest
130 G4ThreeVector dummy;
131 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
132 // create G4Decayproducts
133 auto products = new G4DecayProducts(*parentparticle);
134 delete parentparticle;
135
136 // calculate daughter momentum
137 G4double daughtermomentum[N_DAUGHTER];
138 // calculate electron energy
139 G4double xmax = (1.0 + daughtermass[0] * daughtermass[0] / parentmass / parentmass);
140 G4double x;
141
142 G4double Ee, Ene;
143
144 G4double gam;
145 G4double EMax = parentmass / 2 - daughtermass[0];
146
147 const std::size_t MAX_LOOP = 1000;
148 // Generating Random Energy
149 for (std::size_t loop1 = 0; loop1 < MAX_LOOP; ++loop1) {
150 Ee = G4UniformRand();
151 for (std::size_t loop2 = 0; loop2 < MAX_LOOP; ++loop2) {
152 x = xmax * G4UniformRand();
153 gam = G4UniformRand();
154 if (gam <= x * (1. - x)) break;
155 x = xmax;
156 }
157 Ene = x;
158 if (Ene >= (1. - Ee)) break;
159 Ene = 1. - Ee;
160 }
161 G4double Enm = (2. - Ee - Ene);
162
163 // initialisation of rotation parameters
164
165 G4double costheta, sintheta, rphi, rtheta, rpsi;
166 costheta = 1. - 2. / Ee - 2. / Ene + 2. / Ene / Ee;
167 sintheta = std::sqrt(1. - costheta * costheta);
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 * EMax * EMax + 2.0 * Ee * EMax * daughtermass[0]);
178 G4ThreeVector direction0(0.0, 0.0, 1.0);
179
180 direction0 *= rot;
181
182 auto daughterparticle =
183 new G4DynamicParticle(G4MT_daughters[0], direction0 * daughtermomentum[0]);
184
185 products->PushProducts(daughterparticle);
186
187 // electronic neutrino 1
188
189 daughtermomentum[1] = std::sqrt(Ene * Ene * EMax * EMax + 2.0 * Ene * EMax * daughtermass[1]);
190 G4ThreeVector direction1(sintheta, 0.0, costheta);
191
192 direction1 *= rot;
193
194 auto daughterparticle1 =
195 new G4DynamicParticle(G4MT_daughters[1], direction1 * daughtermomentum[1]);
196 products->PushProducts(daughterparticle1);
197
198 // muonnic neutrino 2
199
200 daughtermomentum[2] = std::sqrt(Enm * Enm * EMax * EMax + 2.0 * Enm * EMax * daughtermass[2]);
201 G4ThreeVector direction2(-Ene / Enm * sintheta, 0, -Ee / Enm - Ene / Enm * costheta);
202
203 direction2 *= rot;
204
205 auto daughterparticle2 =
206 new G4DynamicParticle(G4MT_daughters[2], direction2 * daughtermomentum[2]);
207 products->PushProducts(daughterparticle2);
208
209 // output message
210 #ifdef G4VERBOSE
211 if (GetVerboseLevel() > 1) {
212 G4cout << "G4MuonDecayChannel::DecayIt()";
213 G4cout << " create decay products in rest frame " << G4endl;
214 products->DumpInfo();
215 }
216 #endif
217 return products;
218 }
219