ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KL3DecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4KL3DecayChannel.cc
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 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class header file
31 //
32 // History: first implementation, based on object model of
33 // 30 May 1997 H.Kurashige
34 // ------------------------------------------------------------
35 
36 #include "G4ParticleDefinition.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.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 G4KL3DecayChannel::G4KL3DecayChannel()
47  :G4VDecayChannel(),
48  pLambda(0.0), pXi0(0.0)
49 {
50 }
51 
52 
53 G4KL3DecayChannel::G4KL3DecayChannel(
54  const G4String& theParentName,
55  G4double theBR,
56  const G4String& thePionName,
57  const G4String& theLeptonName,
58  const G4String& theNutrinoName)
59  :G4VDecayChannel("KL3 Decay",theParentName,
60  theBR, 3,
61  thePionName,theLeptonName,theNutrinoName)
62 {
63  static const G4String K_plus("kaon+");
64  static const G4String K_minus("kaon-");
65  static const G4String K_L("kaon0L");
66  static const G4String Mu_plus("mu+");
67  static const G4String Mu_minus("mu-");
68  static const G4String E_plus("e+");
69  static const G4String E_minus("e-");
70 
71  // check modes
72  if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
73  ((theParentName == K_minus)&&(theLeptonName == E_minus)) ) {
74  // K+- (Ke3)
75  pLambda = 0.0286;
76  pXi0 = -0.35;
77  } else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
78  ((theParentName == K_minus)&&(theLeptonName == Mu_minus)) ) {
79  // K+- (Kmu3)
80  pLambda = 0.033;
81  pXi0 = -0.35;
82  } else if ( (theParentName == K_L) &&
83  ((theLeptonName == E_plus) ||(theLeptonName == E_minus)) ){
84  // K0L (Ke3)
85  pLambda = 0.0300;
86  pXi0 = -0.11;
87  } else if ( (theParentName == K_L) &&
88  ((theLeptonName == Mu_plus) ||(theLeptonName == Mu_minus)) ){
89  // K0L (Kmu3)
90  pLambda = 0.034;
91  pXi0 = -0.11;
92  } else {
93 #ifdef G4VERBOSE
94  if (GetVerboseLevel()>2) {
95  G4cout << "G4KL3DecayChannel:: constructor :";
96  G4cout << "illegal arguments " << G4endl;;
97  DumpInfo();
98  }
99 #endif
100  // set values for K0L (Ke3) temporarily
101  pLambda = 0.0300;
102  pXi0 = -0.11;
103  }
104 }
105 
106 G4KL3DecayChannel::~G4KL3DecayChannel()
107 {
108 }
109 
110 G4KL3DecayChannel::G4KL3DecayChannel(const G4KL3DecayChannel &right):
111  G4VDecayChannel(right),
112  //massK(right.massK),
113  pLambda(right.pLambda),
114  pXi0(right.pXi0)
115 {
116 }
117 
118 G4KL3DecayChannel & G4KL3DecayChannel::operator=(const G4KL3DecayChannel & right)
119 {
120  if (this != &right) {
121  kinematics_name = right.kinematics_name;
122  verboseLevel = right.verboseLevel;
123  rbranch = right.rbranch;
124 
125  // copy parent name
126  parent_name = new G4String(*right.parent_name);
127 
128  // clear daughters_name array
129  ClearDaughtersName();
130 
131  // recreate array
132  numberOfDaughters = right.numberOfDaughters;
133  if ( numberOfDaughters >0 ) {
134  if (daughters_name !=0) ClearDaughtersName();
135  daughters_name = new G4String*[numberOfDaughters];
136  //copy daughters name
137  for (G4int index=0; index < numberOfDaughters; index++) {
138  daughters_name[index] = new G4String(*right.daughters_name[index]);
139  }
140  }
141  //massK = right.massK;
142  pLambda = right.pLambda;
143  pXi0 = right.pXi0;
144  }
145  return *this;
146 }
147 
148 
149 G4DecayProducts* G4KL3DecayChannel::DecayIt(G4double)
150 {
151  // this version neglects muon polarization
152  // assumes the pure V-A coupling
153  // gives incorrect energy spectrum for Nutrinos
154 #ifdef G4VERBOSE
155  if (GetVerboseLevel()>1) G4cout << "G4KL3DecayChannel::DecayIt " << G4endl;
156 #endif
157 
158  // fill parent particle and its mass
159  CheckAndFillParent();
160  G4double massK = G4MT_parent->GetPDGMass();
161 
162  // fill daughter particles and their mass
163  CheckAndFillDaughters();
164  G4double daughterM[3];
165  daughterM[idPi] = G4MT_daughters[idPi]->GetPDGMass();
166  daughterM[idLepton] = G4MT_daughters[idLepton]->GetPDGMass();
167  daughterM[idNutrino] = G4MT_daughters[idNutrino]->GetPDGMass();
168 
169  // determine momentum/energy of daughters
170  // according to DalitzDensity
171  G4double daughterP[3], daughterE[3];
172  G4double w;
173  G4double r;
174  const size_t MAX_LOOP = 10000;
175  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
176  r = G4UniformRand();
177  PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
178  w = DalitzDensity(massK,daughterE[idPi],daughterE[idLepton],daughterE[idNutrino],
179  daughterM[idPi],daughterM[idLepton],daughterM[idNutrino]);
180  if ( r <= w) break;
181  }
182 
183  // output message
184 #ifdef G4VERBOSE
185  if (GetVerboseLevel()>1) {
186  G4cout << *daughters_name[0] << ":" << daughterP[0]/GeV << "[GeV/c]" <<G4endl;
187  G4cout << *daughters_name[1] << ":" << daughterP[1]/GeV << "[GeV/c]" <<G4endl;
188  G4cout << *daughters_name[2] << ":" << daughterP[2]/GeV << "[GeV/c]" <<G4endl;
189  }
190 #endif
191  //create parent G4DynamicParticle at rest
192  G4ThreeVector* direction = new G4ThreeVector(1.0,0.0,0.0);
193  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, *direction, 0.0);
194  delete direction;
195 
196  //create G4Decayproducts
197  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
198  delete parentparticle;
199 
200  //create daughter G4DynamicParticle
201  G4double costheta, sintheta, phi, sinphi, cosphi;
202  G4double costhetan, sinthetan, phin, sinphin, cosphin;
203 
204  // pion
205  costheta = 2.*G4UniformRand()-1.0;
206  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
207  phi = twopi*G4UniformRand()*rad;
208  sinphi = std::sin(phi);
209  cosphi = std::cos(phi);
210  direction = new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
211  G4ThreeVector momentum0 = (*direction)*daughterP[0];
212  G4DynamicParticle * daughterparticle
213  = new G4DynamicParticle( G4MT_daughters[0], momentum0);
214  products->PushProducts(daughterparticle);
215 
216  // neutrino
217  costhetan = (daughterP[1]*daughterP[1]-daughterP[2]*daughterP[2]-daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]);
218  sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
219  phin = twopi*G4UniformRand()*rad;
220  sinphin = std::sin(phin);
221  cosphin = std::cos(phin);
222  direction->setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
223  direction->setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
224  direction->setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
225 
226  G4ThreeVector momentum2 = (*direction)*daughterP[2];
227  daughterparticle = new G4DynamicParticle( G4MT_daughters[2], momentum2);
228  products->PushProducts(daughterparticle);
229 
230  //lepton
231  G4ThreeVector momentum1 = (momentum0 + momentum2) * (-1.0);
232  daughterparticle =
233  new G4DynamicParticle( G4MT_daughters[1], momentum1);
234  products->PushProducts(daughterparticle);
235 
236 #ifdef G4VERBOSE
237  if (GetVerboseLevel()>1) {
238  G4cout << "G4KL3DecayChannel::DecayIt ";
239  G4cout << " create decay products in rest frame " <<G4endl;
240  G4cout << " decay products address=" << products << G4endl;
241  products->DumpInfo();
242  }
243 #endif
244  delete direction;
245  return products;
246 }
247 
248 void G4KL3DecayChannel::PhaseSpace(G4double parentM,
249  const G4double* M,
250  G4double* E,
251  G4double* P )
252 // algorism of this code is originally written in GDECA3 of GEANT3
253 {
254 
255  //sum of daughters'mass
256  G4double sumofdaughtermass = 0.0;
257  G4int index;
258  const G4int N_DAUGHTER=3;
259 
260  for (index=0; index<N_DAUGHTER; index++){
261  sumofdaughtermass += M[index];
262  }
263 
264  //calculate daughter momentum
265  // Generate two
266  G4double rd1, rd2, rd;
267  G4double momentummax=0.0, momentumsum = 0.0;
268  G4double energy;
269  const size_t MAX_LOOP=10000;
270  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
271  rd1 = G4UniformRand();
272  rd2 = G4UniformRand();
273  if (rd2 > rd1) {
274  rd = rd1;
275  rd1 = rd2;
276  rd2 = rd;
277  }
278  momentummax = 0.0;
279  momentumsum = 0.0;
280  // daughter 0
281  energy = rd2*(parentM - sumofdaughtermass);
282  P[0] = std::sqrt(energy*energy + 2.0*energy*M[0]);
283  E[0] = energy;
284  if ( P[0] >momentummax )momentummax = P[0];
285  momentumsum += P[0];
286  // daughter 1
287  energy = (1.-rd1)*(parentM - sumofdaughtermass);
288  P[1] = std::sqrt(energy*energy + 2.0*energy*M[1]);
289  E[1] = energy;
290  if ( P[1] >momentummax )momentummax = P[1];
291  momentumsum += P[1];
292  // daughter 2
293  energy = (rd1-rd2)*(parentM - sumofdaughtermass);
294  P[2] = std::sqrt(energy*energy + 2.0*energy*M[2]);
295  E[2] = energy;
296  if ( P[2] >momentummax )momentummax = P[2];
297  momentumsum += P[2];
298  if (momentummax <= momentumsum - momentummax ) break;
299  }
300 #ifdef G4VERBOSE
301  if (GetVerboseLevel()>2) {
302  G4cout << "G4KL3DecayChannel::PhaseSpace ";
303  G4cout << "Kon mass:" << parentM/GeV << "GeV/c/c" << G4endl;
304  for (index=0; index<3; index++){
305  G4cout << index << " : " << M[index]/GeV << "GeV/c/c ";
306  G4cout << " : " << E[index]/GeV << "GeV ";
307  G4cout << " : " << P[index]/GeV << "GeV/c " << G4endl;
308  }
309  }
310 #endif
311 }
312 
313 
314 G4double G4KL3DecayChannel::DalitzDensity(G4double massK, G4double Epi, G4double El, G4double Enu,
315  G4double massPi, G4double massL , G4double massNu )
316 {
317  // KL3 decay Dalitz Plot Density
318  // see Chounet et al Phys. Rep. 4, 201
319  // arguments
320  // Epi: kinetic enregy of pion
321  // El: kinetic enregy of lepton (e or mu)
322  // Enu: kinetic energy of nutrino
323  // constants
324  // pLambda : linear energy dependence of f+
325  // pXi0 : = f+(0)/f-
326  // pNorm : normalization factor
327  // variables
328  // Epi: total energy of pion
329  // El: total energy of lepton (e or mu)
330  // Enu: total energy of nutrino
331 
332  // calcurate total energy
333  Epi = Epi + massPi;
334  El = El + massL;
335  Enu = Enu + massNu;
336 
337  G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;
338  G4double E = Epi_max - Epi;
339  G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
340 
341  G4double F = 1.0 + pLambda*q2/massPi/massPi;
342  G4double Fmax = 1.0;
343  if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0));
344 
345  G4double Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi);
346 
347  G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
348  G4double coeffB = massL*massL*(Enu-E/2.0);
349  G4double coeffC = massL*massL*E/4.0;
350 
351  G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
352 
353  G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
354 
355 #ifdef G4VERBOSE
356  if (GetVerboseLevel()>2) {
357  G4cout << "G4KL3DecayChannel::DalitzDensity " <<G4endl;
358  G4cout << " Pi[" << massPi/GeV <<"GeV/c/c] :" << Epi/GeV << "GeV" <<G4endl;
359  G4cout << " L[" << massL/GeV <<"GeV/c/c] :" << El/GeV << "GeV" <<G4endl;
360  G4cout << " Nu[" << massNu/GeV <<"GeV/c/c] :" << Enu/GeV << "GeV" <<G4endl;
361  G4cout << " F :" << F << " Fmax :" << Fmax << " Xi :" << Xi << G4endl;
362  G4cout << " A :" << coeffA << " B :" << coeffB << " C :"<< coeffC <<G4endl;
363  G4cout << " Rho :" << Rho << " RhoMax :" << RhoMax << G4endl;
364  }
365 #endif
366  return (Rho/RhoMax);
367 }
368 
369