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

Last change on this file since 1315 was 1196, checked in by garnier, 15 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.