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

Last change on this file since 1192 was 992, checked in by garnier, 17 years ago

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