| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4MuonDecayChannelWithSpin class implementa
27 //
28 // References:
29 // - Florian Scheck "Muon Physics", in Physics
30 // (Review Section of Physics Letters) 44, N
31 // 187-248. North-Holland Publishing Company
32 // - W.E. Fisher and F. Scheck, Nucl. Phys. B8
33
34 // Authors: P.Gumplinger and T.MacPhail, 17 Au
35 // -------------------------------------------
36
37 #include "G4MuonDecayChannelWithSpin.hh"
38
39 #include "G4DecayProducts.hh"
40 #include "G4LorentzVector.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh"
44
45 G4MuonDecayChannelWithSpin::G4MuonDecayChannel
46
47 : G4MuonDecayChannel(theParentName, theBR)
48 {}
49
50 G4MuonDecayChannelWithSpin&
51 G4MuonDecayChannelWithSpin::operator=(const G4
52 {
53 if (this != &right) {
54 kinematics_name = right.kinematics_name;
55 verboseLevel = right.verboseLevel;
56 rbranch = right.rbranch;
57
58 // copy parent name
59 delete parent_name;
60 parent_name = new G4String(*right.parent_n
61
62 // clear daughters_name array
63 ClearDaughtersName();
64
65 // recreate array
66 numberOfDaughters = right.numberOfDaughter
67 if (numberOfDaughters > 0) {
68 daughters_name = new G4String*[numberOfD
69 // copy daughters name
70 for (G4int index = 0; index < numberOfDa
71 daughters_name[index] = new G4String(*
72 }
73 }
74 }
75 return *this;
76 }
77
78 G4DecayProducts* G4MuonDecayChannelWithSpin::D
79 {
80 // This version assumes V-A coupling with 1s
81 // the standard model Michel pa
82 // gives incorrect energy spect
83
84 #ifdef G4VERBOSE
85 if (GetVerboseLevel() > 1) G4cout << "G4Muon
86 #endif
87
88 CheckAndFillParent();
89 CheckAndFillDaughters();
90
91 // parent mass
92 G4double parentmass = G4MT_parent->GetPDGMas
93
94 G4double EMMU = parentmass;
95
96 // daughters'mass
97 G4double daughtermass[3];
98 // G4double sumofdaughtermass = 0.0;
99 for (G4int index = 0; index < 3; ++index) {
100 daughtermass[index] = G4MT_daughters[index
101 // sumofdaughtermass += daughtermass[index
102 }
103
104 G4double EMASS = daughtermass[0];
105
106 // create parent G4DynamicParticle at rest
107 G4ThreeVector dummy;
108 auto parentparticle = new G4DynamicParticle(
109 // create G4Decayproducts
110 auto products = new G4DecayProducts(*parentp
111 delete parentparticle;
112
113 // calculate electron energy
114
115 G4double michel_rho = 0.75; // Standard Mod
116 G4double michel_delta = 0.75; // Standard M
117 G4double michel_xsi = 1.00; // Standard Mod
118 G4double michel_eta = 0.00; // Standard Mod
119
120 G4double rndm, x, ctheta;
121
122 G4double FG;
123 G4double FG_max = 2.00;
124
125 G4double W_mue = (EMMU * EMMU + EMASS * EMAS
126 G4double x0 = EMASS / W_mue;
127
128 G4double x0_squared = x0 * x0;
129
130 // *****************************************
131 // x0 <= x <= 1. and -1 <= y <= 1
132 //
133 // F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g
134 // *****************************************
135
136 // ***** sampling F(x,y) directly (brute for
137
138 const std::size_t MAX_LOOP = 10000;
139 for (std::size_t loop_count = 0; loop_count
140 // Sample the positron energy by sampling
141
142 rndm = G4UniformRand();
143
144 x = x0 + rndm * (1. - x0);
145
146 G4double x_squared = x * x;
147
148 G4double F_IS, F_AS, G_IS, G_AS;
149
150 F_IS = 1. / 6. * (-2. * x_squared + 3. * x
151 F_AS = 1. / 6. * std::sqrt(x_squared - x0_
152
153 G_IS = 2. / 9. * (michel_rho - 0.75) * (4.
154 G_IS = G_IS + michel_eta * (1. - x) * x0;
155
156 G_AS = 3. * (michel_xsi - 1.) * (1. - x);
157 G_AS =
158 G_AS + 2. * (michel_xsi * michel_delta -
159 G_AS = 1. / 9. * std::sqrt(x_squared - x0_
160
161 F_IS = F_IS + G_IS;
162 F_AS = F_AS + G_AS;
163
164 // *** Radiative Corrections ***
165
166 const G4double omega = std::log(EMMU / EMA
167 G4double R_IS = F_c(x, x0, omega);
168
169 G4double F = 6. * F_IS + R_IS / std::sqrt(
170
171 // *** Radiative Corrections ***
172
173 G4double R_AS = F_theta(x, x0, omega);
174
175 rndm = G4UniformRand();
176
177 ctheta = 2. * rndm - 1.;
178
179 G4double G = 6. * F_AS - R_AS / std::sqrt(
180
181 FG = std::sqrt(x_squared - x0_squared) * F
182
183 if (FG > FG_max) {
184 G4Exception("G4MuonDecayChannelWithSpin:
185 "Problem in Muon Decay: FG >
186 FG_max = FG;
187 }
188
189 rndm = G4UniformRand();
190
191 if (FG >= rndm * FG_max) break;
192 }
193
194 G4double energy = x * W_mue;
195
196 rndm = G4UniformRand();
197
198 G4double phi = twopi * rndm;
199
200 if (energy < EMASS) energy = EMASS;
201
202 // Calculate daughter momentum
203 G4double daughtermomentum[3];
204
205 daughtermomentum[0] = std::sqrt(energy * ene
206
207 G4double stheta = std::sqrt(1. - ctheta * ct
208 G4double cphi = std::cos(phi);
209 G4double sphi = std::sin(phi);
210
211 // Coordinates of the decay positron with re
212 G4double px = stheta * cphi;
213 G4double py = stheta * sphi;
214 G4double pz = ctheta;
215
216 G4ThreeVector direction0(px, py, pz);
217
218 direction0.rotateUz(parent_polarization);
219
220 auto daughterparticle0 =
221 new G4DynamicParticle(G4MT_daughters[0], d
222
223 products->PushProducts(daughterparticle0);
224
225 // daughter 1 ,2 (neutrinos)
226 // create neutrinos in the C.M frame of two
227 G4double energy2 = parentmass - energy;
228 G4double vmass = std::sqrt((energy2 - daught
229 G4double beta = -1.0 * daughtermomentum[0] /
230 G4double costhetan = 2. * G4UniformRand() -
231 G4double sinthetan = std::sqrt((1.0 - costhe
232 G4double phin = twopi * G4UniformRand() * ra
233 G4double sinphin = std::sin(phin);
234 G4double cosphin = std::cos(phin);
235
236 G4ThreeVector direction1(sinthetan * cosphin
237 auto daughterparticle1 = new G4DynamicPartic
238 auto daughterparticle2 =
239 new G4DynamicParticle(G4MT_daughters[2], d
240
241 // boost to the muon rest frame
242 G4LorentzVector p4;
243 p4 = daughterparticle1->Get4Momentum();
244 p4.boost(direction0.x() * beta, direction0.y
245 daughterparticle1->Set4Momentum(p4);
246 p4 = daughterparticle2->Get4Momentum();
247 p4.boost(direction0.x() * beta, direction0.y
248 daughterparticle2->Set4Momentum(p4);
249 products->PushProducts(daughterparticle1);
250 products->PushProducts(daughterparticle2);
251 daughtermomentum[1] = daughterparticle1->Get
252 daughtermomentum[2] = daughterparticle2->Get
253
254 // output message
255 #ifdef G4VERBOSE
256 if (GetVerboseLevel() > 1) {
257 G4cout << "G4MuonDecayChannelWithSpin::Dec
258 G4cout << " create decay products in rest
259 G4double TT = daughterparticle0->GetTotalE
260 + daughterparticle2->GetTota
261 G4cout << "e " << daughterparticle0->GetT
262 G4cout << "nu1" << daughterparticle1->GetT
263 G4cout << "nu2" << daughterparticle2->GetT
264 G4cout << "total" << (TT - parentmass) / k
265 if (GetVerboseLevel() > 2) {
266 products->DumpInfo();
267 }
268 }
269 #endif
270
271 return products;
272 }
273
274 G4double G4MuonDecayChannelWithSpin::R_c(G4dou
275 {
276 auto n_max = (G4int)(100. * x);
277
278 if (n_max < 10) n_max = 10;
279
280 G4double L2 = 0.0;
281
282 for (G4int n = 1; n <= n_max; ++n) {
283 L2 += std::pow(x, n) / (n * n);
284 }
285
286 G4double r_c;
287
288 r_c = 2. * L2 - (pi * pi / 3.) - 2.;
289 r_c = r_c + omega * (1.5 + 2. * std::log((1.
290 r_c = r_c - std::log(x) * (2. * std::log(x)
291 r_c = r_c + (3. * std::log(x) - 1. - 1. / x)
292
293 return r_c;
294 }
295