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

Last change on this file since 1344 was 1340, checked in by garnier, 15 years ago

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