| 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 // G4KL3DecayChannel class implementation
27 //
28 // Author: H.Kurashige, 30 May 1997
29 // --------------------------------------------------------------------
30
31 #include "G4KL3DecayChannel.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 G4KL3DecayChannel::G4KL3DecayChannel(const G4String& theParentName, G4double theBR,
43 const G4String& thePionName, const G4String& theLeptonName,
44 const G4String& theNutrinoName)
45 : G4VDecayChannel("KL3 Decay", theParentName, theBR, 3, thePionName, theLeptonName,
46 theNutrinoName)
47 {
48 static const G4String K_plus("kaon+");
49 static const G4String K_minus("kaon-");
50 static const G4String K_L("kaon0L");
51 static const G4String Mu_plus("mu+");
52 static const G4String Mu_minus("mu-");
53 static const G4String E_plus("e+");
54 static const G4String E_minus("e-");
55
56 // check modes
57 if (((theParentName == K_plus) && (theLeptonName == E_plus))
58 || ((theParentName == K_minus) && (theLeptonName == E_minus)))
59 {
60 // K+- (Ke3)
61 pLambda = 0.0286;
62 pXi0 = -0.35;
63 }
64 else if (((theParentName == K_plus) && (theLeptonName == Mu_plus))
65 || ((theParentName == K_minus) && (theLeptonName == Mu_minus)))
66 {
67 // K+- (Kmu3)
68 pLambda = 0.033;
69 pXi0 = -0.35;
70 }
71 else if ((theParentName == K_L) && ((theLeptonName == E_plus) || (theLeptonName == E_minus))) {
72 // K0L (Ke3)
73 pLambda = 0.0300;
74 pXi0 = -0.11;
75 }
76 else if ((theParentName == K_L) && ((theLeptonName == Mu_plus) || (theLeptonName == Mu_minus))) {
77 // K0L (Kmu3)
78 pLambda = 0.034;
79 pXi0 = -0.11;
80 }
81 else {
82 #ifdef G4VERBOSE
83 if (GetVerboseLevel() > 2) {
84 G4cout << "G4KL3DecayChannel:: constructor :";
85 G4cout << "illegal arguments " << G4endl;
86 ;
87 DumpInfo();
88 }
89 #endif
90 // set values for K0L (Ke3) temporarily
91 pLambda = 0.0300;
92 pXi0 = -0.11;
93 }
94 }
95
96 G4KL3DecayChannel& G4KL3DecayChannel::operator=(const G4KL3DecayChannel& right)
97 {
98 if (this != &right) {
99 kinematics_name = right.kinematics_name;
100 verboseLevel = right.verboseLevel;
101 rbranch = right.rbranch;
102
103 // copy parent name
104 parent_name = new G4String(*right.parent_name);
105
106 // clear daughters_name array
107 ClearDaughtersName();
108
109 // recreate array
110 numberOfDaughters = right.numberOfDaughters;
111 if (numberOfDaughters > 0) {
112 if (daughters_name != nullptr) ClearDaughtersName();
113 daughters_name = new G4String*[numberOfDaughters];
114 // copy daughters name
115 for (G4int index = 0; index < numberOfDaughters; ++index) {
116 daughters_name[index] = new G4String(*right.daughters_name[index]);
117 }
118 }
119 pLambda = right.pLambda;
120 pXi0 = right.pXi0;
121 }
122 return *this;
123 }
124
125 G4DecayProducts* G4KL3DecayChannel::DecayIt(G4double)
126 {
127 // this version neglects muon polarization
128 // assumes the pure V-A coupling
129 // gives incorrect energy spectrum for Neutrinos
130 #ifdef G4VERBOSE
131 if (GetVerboseLevel() > 1) G4cout << "G4KL3DecayChannel::DecayIt " << G4endl;
132 #endif
133
134 // fill parent particle and its mass
135 CheckAndFillParent();
136 G4double massK = G4MT_parent->GetPDGMass();
137
138 // fill daughter particles and their mass
139 CheckAndFillDaughters();
140 G4double daughterM[3];
141 daughterM[idPi] = G4MT_daughters[idPi]->GetPDGMass();
142 daughterM[idLepton] = G4MT_daughters[idLepton]->GetPDGMass();
143 daughterM[idNutrino] = G4MT_daughters[idNutrino]->GetPDGMass();
144
145 // determine momentum/energy of daughters according to DalitzDensity
146 G4double daughterP[3], daughterE[3];
147 G4double w;
148 G4double r;
149 const size_t MAX_LOOP = 10000;
150 for (std::size_t loop_counter = 0; loop_counter < MAX_LOOP; ++loop_counter) {
151 r = G4UniformRand();
152 PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
153 w = DalitzDensity(massK, daughterE[idPi], daughterE[idLepton], daughterE[idNutrino],
154 daughterM[idPi], daughterM[idLepton], daughterM[idNutrino]);
155 if (r <= w) break;
156 }
157
158 // output message
159 #ifdef G4VERBOSE
160 if (GetVerboseLevel() > 1) {
161 G4cout << *daughters_name[0] << ":" << daughterP[0] / GeV << "[GeV/c]" << G4endl;
162 G4cout << *daughters_name[1] << ":" << daughterP[1] / GeV << "[GeV/c]" << G4endl;
163 G4cout << *daughters_name[2] << ":" << daughterP[2] / GeV << "[GeV/c]" << G4endl;
164 }
165 #endif
166
167 // create parent G4DynamicParticle at rest
168 auto direction = new G4ThreeVector(1.0, 0.0, 0.0);
169 auto parentparticle = new G4DynamicParticle(G4MT_parent, *direction, 0.0);
170 delete direction;
171
172 // create G4Decayproducts
173 auto products = new G4DecayProducts(*parentparticle);
174 delete parentparticle;
175
176 // create daughter G4DynamicParticle
177 G4double costheta, sintheta, phi, sinphi, cosphi;
178 G4double costhetan, sinthetan, phin, sinphin, cosphin;
179
180 // pion
181 costheta = 2. * G4UniformRand() - 1.0;
182 sintheta = std::sqrt((1.0 - costheta) * (1.0 + costheta));
183 phi = twopi * G4UniformRand() * rad;
184 sinphi = std::sin(phi);
185 cosphi = std::cos(phi);
186 direction = new G4ThreeVector(sintheta * cosphi, sintheta * sinphi, costheta);
187 G4ThreeVector momentum0 = (*direction) * daughterP[0];
188 auto daughterparticle = new G4DynamicParticle(G4MT_daughters[0], momentum0);
189 products->PushProducts(daughterparticle);
190
191 // neutrino
192 costhetan =
193 (daughterP[1] * daughterP[1] - daughterP[2] * daughterP[2] - daughterP[0] * daughterP[0])
194 / (2.0 * daughterP[2] * daughterP[0]);
195 sinthetan = std::sqrt((1.0 - costhetan) * (1.0 + costhetan));
196 phin = twopi * G4UniformRand() * rad;
197 sinphin = std::sin(phin);
198 cosphin = std::cos(phin);
199 direction->setX(sinthetan * cosphin * costheta * cosphi - sinthetan * sinphin * sinphi
200 + costhetan * sintheta * cosphi);
201 direction->setY(sinthetan * cosphin * costheta * sinphi + sinthetan * sinphin * cosphi
202 + costhetan * sintheta * sinphi);
203 direction->setZ(-sinthetan * cosphin * sintheta + costhetan * costheta);
204
205 G4ThreeVector momentum2 = (*direction) * daughterP[2];
206 daughterparticle = new G4DynamicParticle(G4MT_daughters[2], momentum2);
207 products->PushProducts(daughterparticle);
208
209 // lepton
210 G4ThreeVector momentum1 = (momentum0 + momentum2) * (-1.0);
211 daughterparticle = new G4DynamicParticle(G4MT_daughters[1], momentum1);
212 products->PushProducts(daughterparticle);
213
214 #ifdef G4VERBOSE
215 if (GetVerboseLevel() > 1) {
216 G4cout << "G4KL3DecayChannel::DecayIt ";
217 G4cout << " create decay products in rest frame " << G4endl;
218 G4cout << " decay products address=" << products << G4endl;
219 products->DumpInfo();
220 }
221 #endif
222 delete direction;
223 return products;
224 }
225
226 void G4KL3DecayChannel::PhaseSpace(G4double parentM, const G4double* M, G4double* E, G4double* P)
227 {
228 // Algorithm in this code was originally written in GDECA3 in GEANT3
229
230 // sum of daughters'mass
231 G4double sumofdaughtermass = 0.0;
232 G4int index;
233 const G4int N_DAUGHTER = 3;
234
235 for (index = 0; index < N_DAUGHTER; ++index) {
236 sumofdaughtermass += M[index];
237 }
238
239 // calculate daughter momentum. Generate two
240 G4double rd1, rd2, rd;
241 G4double momentummax = 0.0, momentumsum = 0.0;
242 G4double energy;
243 const size_t MAX_LOOP = 10000;
244 for (std::size_t loop_counter = 0; loop_counter < MAX_LOOP; ++loop_counter) {
245 rd1 = G4UniformRand();
246 rd2 = G4UniformRand();
247 if (rd2 > rd1) {
248 rd = rd1;
249 rd1 = rd2;
250 rd2 = rd;
251 }
252 momentummax = 0.0;
253 momentumsum = 0.0;
254 // daughter 0
255 energy = rd2 * (parentM - sumofdaughtermass);
256 P[0] = std::sqrt(energy * energy + 2.0 * energy * M[0]);
257 E[0] = energy;
258 if (P[0] > momentummax) momentummax = P[0];
259 momentumsum += P[0];
260 // daughter 1
261 energy = (1. - rd1) * (parentM - sumofdaughtermass);
262 P[1] = std::sqrt(energy * energy + 2.0 * energy * M[1]);
263 E[1] = energy;
264 if (P[1] > momentummax) momentummax = P[1];
265 momentumsum += P[1];
266 // daughter 2
267 energy = (rd1 - rd2) * (parentM - sumofdaughtermass);
268 P[2] = std::sqrt(energy * energy + 2.0 * energy * M[2]);
269 E[2] = energy;
270 if (P[2] > momentummax) momentummax = P[2];
271 momentumsum += P[2];
272 if (momentummax <= momentumsum - momentummax) break;
273 }
274 #ifdef G4VERBOSE
275 if (GetVerboseLevel() > 2) {
276 G4cout << "G4KL3DecayChannel::PhaseSpace ";
277 G4cout << "Kon mass:" << parentM / GeV << "GeV/c/c" << G4endl;
278 for (index = 0; index < 3; ++index) {
279 G4cout << index << " : " << M[index] / GeV << "GeV/c/c ";
280 G4cout << " : " << E[index] / GeV << "GeV ";
281 G4cout << " : " << P[index] / GeV << "GeV/c " << G4endl;
282 }
283 }
284 #endif
285 }
286
287 G4double G4KL3DecayChannel::DalitzDensity(G4double massK, G4double Epi, G4double El, G4double Enu,
288 G4double massPi, G4double massL, G4double massNu)
289 {
290 // KL3 decay - Dalitz Plot Density, see Chounet et al Phys. Rep. 4, 201
291 // Arguments
292 // Epi: kinetic enregy of pion
293 // El: kinetic enregy of lepton (e or mu)
294 // Enu: kinetic energy of nutrino
295 // Constants
296 // pLambda : linear energy dependence of f+
297 // pXi0 : = f+(0)/f-
298 // pNorm : normalization factor
299 // Variables
300 // Epi: total energy of pion
301 // El: total energy of lepton (e or mu)
302 // Enu: total energy of nutrino
303
304 // calculate total energy
305 Epi = Epi + massPi;
306 El = El + massL;
307 Enu = Enu + massNu;
308
309 G4double Epi_max = (massK * massK + massPi * massPi - massL * massL) / 2.0 / massK;
310 G4double E = Epi_max - Epi;
311 G4double q2 = massK * massK + massPi * massPi - 2.0 * massK * Epi;
312
313 G4double F = 1.0 + pLambda * q2 / massPi / massPi;
314 G4double Fmax = 1.0;
315 if (pLambda > 0.0) Fmax = (1.0 + pLambda * (massK * massK / massPi / massPi + 1.0));
316
317 G4double Xi = pXi0 * (1.0 + pLambda * q2 / massPi / massPi);
318
319 G4double coeffA = massK * (2.0 * El * Enu - massK * E) + massL * massL * (E / 4.0 - Enu);
320 G4double coeffB = massL * massL * (Enu - E / 2.0);
321 G4double coeffC = massL * massL * E / 4.0;
322
323 G4double RhoMax = (Fmax * Fmax) * (massK * massK * massK / 8.0);
324
325 G4double Rho = (F * F) * (coeffA + coeffB * Xi + coeffC * Xi * Xi);
326
327 #ifdef G4VERBOSE
328 if (GetVerboseLevel() > 2) {
329 G4cout << "G4KL3DecayChannel::DalitzDensity " << G4endl;
330 G4cout << " Pi[" << massPi / GeV << "GeV/c/c] :" << Epi / GeV << "GeV" << G4endl;
331 G4cout << " L[" << massL / GeV << "GeV/c/c] :" << El / GeV << "GeV" << G4endl;
332 G4cout << " Nu[" << massNu / GeV << "GeV/c/c] :" << Enu / GeV << "GeV" << G4endl;
333 G4cout << " F :" << F << " Fmax :" << Fmax << " Xi :" << Xi << G4endl;
334 G4cout << " A :" << coeffA << " B :" << coeffB << " C :" << coeffC << G4endl;
335 G4cout << " Rho :" << Rho << " RhoMax :" << RhoMax << G4endl;
336 }
337 #endif
338 return (Rho / RhoMax);
339 }
340