Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KL3DecayChannel Class Reference

#include <G4KL3DecayChannel.hh>

Inheritance diagram for G4KL3DecayChannel:
Collaboration diagram for G4KL3DecayChannel:

Public Member Functions

 G4KL3DecayChannel (const G4String &theParentName, G4double theBR, const G4String &thePionName, const G4String &theLeptonName, const G4String &theNutrinoName)
 
virtual ~G4KL3DecayChannel ()
 
virtual G4DecayProductsDecayIt (G4double)
 
void SetDalitzParameter (G4double aLambda, G4double aXi)
 
G4double GetDalitzParameterLambda () const
 
G4double GetDalitzParameterXi () const
 
- Public Member Functions inherited from G4VDecayChannel
 G4VDecayChannel (const G4String &aName, G4int Verbose=1)
 
 G4VDecayChannel (const G4String &aName, const G4String &theParentName, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="", const G4String &theDaughterName4="")
 
virtual ~G4VDecayChannel ()
 
G4int operator== (const G4VDecayChannel &right) const
 
G4int operator!= (const G4VDecayChannel &right) const
 
G4int operator< (const G4VDecayChannel &right) const
 
const G4StringGetKinematicsName () const
 
G4double GetBR () const
 
G4int GetNumberOfDaughters () const
 
G4ParticleDefinitionGetParent ()
 
G4ParticleDefinitionGetDaughter (G4int anIndex)
 
G4int GetAngularMomentum ()
 
const G4StringGetParentName () const
 
const G4StringGetDaughterName (G4int anIndex) const
 
G4double GetParentMass () const
 
G4double GetDaughterMass (G4int anIndex) const
 
void SetParent (const G4ParticleDefinition *particle_type)
 
void SetParent (const G4String &particle_name)
 
void SetBR (G4double value)
 
void SetNumberOfDaughters (G4int value)
 
void SetDaughter (G4int anIndex, const G4ParticleDefinition *particle_type)
 
void SetDaughter (G4int anIndex, const G4String &particle_name)
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 
void DumpInfo ()
 
G4double GetRangeMass () const
 
void SetRangeMass (G4double val)
 
virtual G4bool IsOKWithParentMass (G4double parentMass)
 
void SetPolarization (const G4ThreeVector &)
 
const G4ThreeVectorGetPolarization () const
 

Protected Types

enum  { idPi =0, idLepton =1, idNutrino =2 }
 

Protected Member Functions

 G4KL3DecayChannel (const G4KL3DecayChannel &)
 
G4KL3DecayChanneloperator= (const G4KL3DecayChannel &)
 
void PhaseSpace (G4double Mparent, const G4double *Mdaughter, G4double *Edaughter, G4double *Pdaughter)
 
G4double DalitzDensity (G4double parentmass, G4double Epi, G4double El, G4double Enu, G4double massPi, G4double massL, G4double massNu)
 
- Protected Member Functions inherited from G4VDecayChannel
void ClearDaughtersName ()
 
void CheckAndFillDaughters ()
 
void CheckAndFillParent ()
 
 G4VDecayChannel ()
 
 G4VDecayChannel (const G4VDecayChannel &)
 
G4VDecayChanneloperator= (const G4VDecayChannel &)
 
G4double DynamicalMass (G4double massPDG, G4double width, G4double maxDev=+1.) const
 

Additional Inherited Members

- Protected Attributes inherited from G4VDecayChannel
G4String kinematics_name
 
G4double rbranch
 
G4int numberOfDaughters
 
G4Stringparent_name
 
G4String ** daughters_name
 
G4double rangeMass
 
G4ThreeVector parent_polarization
 
G4ParticleTableparticletable
 
G4int verboseLevel
 
G4ParticleDefinitionG4MT_parent
 
G4ParticleDefinition ** G4MT_daughters
 
G4double G4MT_parent_mass
 
G4doubleG4MT_daughters_mass
 
G4doubleG4MT_daughters_width
 
G4Mutex daughtersMutex
 
G4Mutex parentMutex
 
- Static Protected Attributes inherited from G4VDecayChannel
static const G4String noName = " "
 

Detailed Description

Definition at line 43 of file G4KL3DecayChannel.hh.

Member Enumeration Documentation

anonymous enum
protected
Enumerator
idPi 
idLepton 
idNutrino 

Definition at line 68 of file G4KL3DecayChannel.hh.

Constructor & Destructor Documentation

G4KL3DecayChannel::G4KL3DecayChannel ( const G4String theParentName,
G4double  theBR,
const G4String thePionName,
const G4String theLeptonName,
const G4String theNutrinoName 
)

Definition at line 54 of file G4KL3DecayChannel.cc.

60  :G4VDecayChannel("KL3 Decay",theParentName,
61  theBR, 3,
62  thePionName,theLeptonName,theNutrinoName)
63 {
64  static const G4String K_plus("kaon+");
65  static const G4String K_minus("kaon-");
66  static const G4String K_L("kaon0L");
67  static const G4String Mu_plus("mu+");
68  static const G4String Mu_minus("mu-");
69  static const G4String E_plus("e+");
70  static const G4String E_minus("e-");
71 
72  // check modes
73  if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
74  ((theParentName == K_minus)&&(theLeptonName == E_minus)) ) {
75  // K+- (Ke3)
76  pLambda = 0.0286;
77  pXi0 = -0.35;
78  } else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
79  ((theParentName == K_minus)&&(theLeptonName == Mu_minus)) ) {
80  // K+- (Kmu3)
81  pLambda = 0.033;
82  pXi0 = -0.35;
83  } else if ( (theParentName == K_L) &&
84  ((theLeptonName == E_plus) ||(theLeptonName == E_minus)) ){
85  // K0L (Ke3)
86  pLambda = 0.0300;
87  pXi0 = -0.11;
88  } else if ( (theParentName == K_L) &&
89  ((theLeptonName == Mu_plus) ||(theLeptonName == Mu_minus)) ){
90  // K0L (Kmu3)
91  pLambda = 0.034;
92  pXi0 = -0.11;
93  } else {
94 #ifdef G4VERBOSE
95  if (GetVerboseLevel()>2) {
96  G4cout << "G4KL3DecayChannel:: constructor :";
97  G4cout << "illegal arguments " << G4endl;;
98  DumpInfo();
99  }
100 #endif
101  // set values for K0L (Ke3) temporarily
102  pLambda = 0.0300;
103  pXi0 = -0.11;
104  }
105 }
G4GLOB_DLL std::ostream G4cout
G4int GetVerboseLevel() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

G4KL3DecayChannel::~G4KL3DecayChannel ( )
virtual

Definition at line 107 of file G4KL3DecayChannel.cc.

108 {
109 }
G4KL3DecayChannel::G4KL3DecayChannel ( const G4KL3DecayChannel right)
protected

Definition at line 111 of file G4KL3DecayChannel.cc.

111  :
112  G4VDecayChannel(right),
113  //massK(right.massK),
114  pLambda(right.pLambda),
115  pXi0(right.pXi0)
116 {
117 }

Here is the call graph for this function:

Member Function Documentation

G4double G4KL3DecayChannel::DalitzDensity ( G4double  parentmass,
G4double  Epi,
G4double  El,
G4double  Enu,
G4double  massPi,
G4double  massL,
G4double  massNu 
)
protected

Definition at line 315 of file G4KL3DecayChannel.cc.

317 {
318  // KL3 decay Dalitz Plot Density
319  // see Chounet et al Phys. Rep. 4, 201
320  // arguments
321  // Epi: kinetic enregy of pion
322  // El: kinetic enregy of lepton (e or mu)
323  // Enu: kinetic energy of nutrino
324  // constants
325  // pLambda : linear energy dependence of f+
326  // pXi0 : = f+(0)/f-
327  // pNorm : normalization factor
328  // variables
329  // Epi: total energy of pion
330  // El: total energy of lepton (e or mu)
331  // Enu: total energy of nutrino
332 
333  // calcurate total energy
334  Epi = Epi + massPi;
335  El = El + massL;
336  Enu = Enu + massNu;
337 
338  G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;
339  G4double E = Epi_max - Epi;
340  G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
341 
342  G4double F = 1.0 + pLambda*q2/massPi/massPi;
343  G4double Fmax = 1.0;
344  if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0));
345 
346  G4double Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi);
347 
348  G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
349  G4double coeffB = massL*massL*(Enu-E/2.0);
350  G4double coeffC = massL*massL*E/4.0;
351 
352  G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
353 
354  G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
355 
356 #ifdef G4VERBOSE
357  if (GetVerboseLevel()>2) {
358  G4cout << "G4KL3DecayChannel::DalitzDensity " <<G4endl;
359  G4cout << " Pi[" << massPi/GeV <<"GeV/c/c] :" << Epi/GeV << "GeV" <<G4endl;
360  G4cout << " L[" << massL/GeV <<"GeV/c/c] :" << El/GeV << "GeV" <<G4endl;
361  G4cout << " Nu[" << massNu/GeV <<"GeV/c/c] :" << Enu/GeV << "GeV" <<G4endl;
362  G4cout << " F :" << F << " Fmax :" << Fmax << " Xi :" << Xi << G4endl;
363  G4cout << " A :" << coeffA << " B :" << coeffB << " C :"<< coeffC <<G4endl;
364  G4cout << " Rho :" << Rho << " RhoMax :" << RhoMax << G4endl;
365  }
366 #endif
367  return (Rho/RhoMax);
368 }
G4GLOB_DLL std::ostream G4cout
G4int GetVerboseLevel() const
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4DecayProducts * G4KL3DecayChannel::DecayIt ( G4double  )
virtual

Implements G4VDecayChannel.

Definition at line 150 of file G4KL3DecayChannel.cc.

151 {
152  // this version neglects muon polarization
153  // assumes the pure V-A coupling
154  // gives incorrect energy spectrum for Nutrinos
155 #ifdef G4VERBOSE
156  if (GetVerboseLevel()>1) G4cout << "G4KL3DecayChannel::DecayIt " << G4endl;
157 #endif
158 
159  // fill parent particle and its mass
161  G4double massK = G4MT_parent->GetPDGMass();
162 
163  // fill daughter particles and their mass
165  G4double daughterM[3];
166  daughterM[idPi] = G4MT_daughters[idPi]->GetPDGMass();
167  daughterM[idLepton] = G4MT_daughters[idLepton]->GetPDGMass();
168  daughterM[idNutrino] = G4MT_daughters[idNutrino]->GetPDGMass();
169 
170  // determine momentum/energy of daughters
171  // according to DalitzDensity
172  G4double daughterP[3], daughterE[3];
173  G4double w;
174  G4double r;
175  const size_t MAX_LOOP = 10000;
176  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
177  r = G4UniformRand();
178  PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
179  w = DalitzDensity(massK,daughterE[idPi],daughterE[idLepton],daughterE[idNutrino],
180  daughterM[idPi],daughterM[idLepton],daughterM[idNutrino]);
181  if ( r <= w) break;
182  }
183 
184  // output message
185 #ifdef G4VERBOSE
186  if (GetVerboseLevel()>1) {
187  G4cout << *daughters_name[0] << ":" << daughterP[0]/GeV << "[GeV/c]" <<G4endl;
188  G4cout << *daughters_name[1] << ":" << daughterP[1]/GeV << "[GeV/c]" <<G4endl;
189  G4cout << *daughters_name[2] << ":" << daughterP[2]/GeV << "[GeV/c]" <<G4endl;
190  }
191 #endif
192  //create parent G4DynamicParticle at rest
193  G4ThreeVector* direction = new G4ThreeVector(1.0,0.0,0.0);
194  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, *direction, 0.0);
195  delete direction;
196 
197  //create G4Decayproducts
198  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
199  delete parentparticle;
200 
201  //create daughter G4DynamicParticle
202  G4double costheta, sintheta, phi, sinphi, cosphi;
203  G4double costhetan, sinthetan, phin, sinphin, cosphin;
204 
205  // pion
206  costheta = 2.*G4UniformRand()-1.0;
207  sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
208  phi = twopi*G4UniformRand()*rad;
209  sinphi = std::sin(phi);
210  cosphi = std::cos(phi);
211  direction = new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
212  G4ThreeVector momentum0 = (*direction)*daughterP[0];
213  G4DynamicParticle * daughterparticle
214  = new G4DynamicParticle( G4MT_daughters[0], momentum0);
215  products->PushProducts(daughterparticle);
216 
217  // neutrino
218  costhetan = (daughterP[1]*daughterP[1]-daughterP[2]*daughterP[2]-daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]);
219  sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
220  phin = twopi*G4UniformRand()*rad;
221  sinphin = std::sin(phin);
222  cosphin = std::cos(phin);
223  direction->setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
224  direction->setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
225  direction->setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
226 
227  G4ThreeVector momentum2 = (*direction)*daughterP[2];
228  daughterparticle = new G4DynamicParticle( G4MT_daughters[2], momentum2);
229  products->PushProducts(daughterparticle);
230 
231  //lepton
232  G4ThreeVector momentum1 = (momentum0 + momentum2) * (-1.0);
233  daughterparticle =
234  new G4DynamicParticle( G4MT_daughters[1], momentum1);
235  products->PushProducts(daughterparticle);
236 
237 #ifdef G4VERBOSE
238  if (GetVerboseLevel()>1) {
239  G4cout << "G4KL3DecayChannel::DecayIt ";
240  G4cout << " create decay products in rest frame " <<G4endl;
241  G4cout << " decay products address=" << products << G4endl;
242  products->DumpInfo();
243  }
244 #endif
245  delete direction;
246  return products;
247 }
void CheckAndFillDaughters()
CLHEP::Hep3Vector G4ThreeVector
G4int PushProducts(G4DynamicParticle *aParticle)
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
G4double DalitzDensity(G4double parentmass, G4double Epi, G4double El, G4double Enu, G4double massPi, G4double massL, G4double massNu)
static constexpr double rad
Definition: G4SIunits.hh:149
void setY(double)
void setZ(double)
void setX(double)
static constexpr double twopi
Definition: G4SIunits.hh:76
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void DumpInfo() const
G4int GetVerboseLevel() const
G4double GetPDGMass() const
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4String ** daughters_name
void PhaseSpace(G4double Mparent, const G4double *Mdaughter, G4double *Edaughter, G4double *Pdaughter)

Here is the call graph for this function:

G4double G4KL3DecayChannel::GetDalitzParameterLambda ( ) const
inline

Definition at line 113 of file G4KL3DecayChannel.hh.

114 {
115  return pLambda;
116 }
G4double G4KL3DecayChannel::GetDalitzParameterXi ( ) const
inline

Definition at line 119 of file G4KL3DecayChannel.hh.

120 {
121  return pXi0;
122 }
G4KL3DecayChannel & G4KL3DecayChannel::operator= ( const G4KL3DecayChannel right)
protected

Definition at line 119 of file G4KL3DecayChannel.cc.

120 {
121  if (this != &right) {
123  verboseLevel = right.verboseLevel;
124  rbranch = right.rbranch;
125 
126  // copy parent name
127  parent_name = new G4String(*right.parent_name);
128 
129  // clear daughters_name array
131 
132  // recreate array
134  if ( numberOfDaughters >0 ) {
137  //copy daughters name
138  for (G4int index=0; index < numberOfDaughters; index++) {
139  daughters_name[index] = new G4String(*right.daughters_name[index]);
140  }
141  }
142  //massK = right.massK;
143  pLambda = right.pLambda;
144  pXi0 = right.pXi0;
145  }
146  return *this;
147 }
int G4int
Definition: G4Types.hh:78
G4String kinematics_name
G4String * parent_name
G4String ** daughters_name

Here is the call graph for this function:

void G4KL3DecayChannel::PhaseSpace ( G4double  Mparent,
const G4double Mdaughter,
G4double Edaughter,
G4double Pdaughter 
)
protected

Definition at line 249 of file G4KL3DecayChannel.cc.

254 {
255 
256  //sum of daughters'mass
257  G4double sumofdaughtermass = 0.0;
258  G4int index;
259  const G4int N_DAUGHTER=3;
260 
261  for (index=0; index<N_DAUGHTER; index++){
262  sumofdaughtermass += M[index];
263  }
264 
265  //calculate daughter momentum
266  // Generate two
267  G4double rd1, rd2, rd;
268  G4double momentummax=0.0, momentumsum = 0.0;
270  const size_t MAX_LOOP=10000;
271  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
272  rd1 = G4UniformRand();
273  rd2 = G4UniformRand();
274  if (rd2 > rd1) {
275  rd = rd1;
276  rd1 = rd2;
277  rd2 = rd;
278  }
279  momentummax = 0.0;
280  momentumsum = 0.0;
281  // daughter 0
282  energy = rd2*(parentM - sumofdaughtermass);
283  P[0] = std::sqrt(energy*energy + 2.0*energy*M[0]);
284  E[0] = energy;
285  if ( P[0] >momentummax )momentummax = P[0];
286  momentumsum += P[0];
287  // daughter 1
288  energy = (1.-rd1)*(parentM - sumofdaughtermass);
289  P[1] = std::sqrt(energy*energy + 2.0*energy*M[1]);
290  E[1] = energy;
291  if ( P[1] >momentummax )momentummax = P[1];
292  momentumsum += P[1];
293  // daughter 2
294  energy = (rd1-rd2)*(parentM - sumofdaughtermass);
295  P[2] = std::sqrt(energy*energy + 2.0*energy*M[2]);
296  E[2] = energy;
297  if ( P[2] >momentummax )momentummax = P[2];
298  momentumsum += P[2];
299  if (momentummax <= momentumsum - momentummax ) break;
300  }
301 #ifdef G4VERBOSE
302  if (GetVerboseLevel()>2) {
303  G4cout << "G4KL3DecayChannel::PhaseSpace ";
304  G4cout << "Kon mass:" << parentM/GeV << "GeV/c/c" << G4endl;
305  for (index=0; index<3; index++){
306  G4cout << index << " : " << M[index]/GeV << "GeV/c/c ";
307  G4cout << " : " << E[index]/GeV << "GeV ";
308  G4cout << " : " << P[index]/GeV << "GeV/c " << G4endl;
309  }
310  }
311 #endif
312 }
int G4int
Definition: G4Types.hh:78
static double P[]
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4int GetVerboseLevel() const
G4double energy(const ThreeVector &p, const G4double m)
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4KL3DecayChannel::SetDalitzParameter ( G4double  aLambda,
G4double  aXi 
)
inline

Definition at line 106 of file G4KL3DecayChannel.hh.

107 {
108  pLambda = aLambda;
109  pXi0 = aXi;
110 }

The documentation for this class was generated from the following files: