source: trunk/source/particles/management/src/G4KL3DecayChannel.cc@ 1216

Last change on this file since 1216 was 1196, checked in by garnier, 16 years ago

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