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

Last change on this file since 835 was 824, checked in by garnier, 16 years ago

import all except CVS
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:  $
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.