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

Last change on this file since 1355 was 1340, checked in by garnier, 14 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.