Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Protected Types | Protected Member Functions | List of all members
G4KL3DecayChannel Class Reference

#include <G4KL3DecayChannel.hh>

Inheritance diagram for G4KL3DecayChannel:
Inheritance graph
[legend]
Collaboration diagram for G4KL3DecayChannel:
Collaboration graph
[legend]

Public Member Functions

 G4KL3DecayChannel (const G4String &theParentName, G4double theBR, const G4String &thePionName, const G4String &theLeptonName, const G4String &theNutrinoName)
 
virtual ~G4KL3DecayChannel ()
 
virtual G4DecayProducts * DecayIt (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 G4String & GetKinematicsName () const
 
G4double GetBR () const
 
G4int GetNumberOfDaughters () const
 
G4ParticleDefinition * GetParent ()
 
G4ParticleDefinition * GetDaughter (G4int anIndex)
 
G4int GetAngularMomentum ()
 
const G4String & GetParentName () const
 
const G4String & GetDaughterName (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 G4ThreeVector & GetPolarization () const
 

Protected Types

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

Protected Member Functions

 G4KL3DecayChannel (const G4KL3DecayChannel &)
 
G4KL3DecayChannel & operator= (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 &)
 
G4VDecayChannel & operator= (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
 
G4String * parent_name
 
G4String ** daughters_name
 
G4double rangeMass
 
G4ThreeVector parent_polarization
 
G4ParticleTable * particletable
 
G4int verboseLevel
 
G4ParticleDefinition * G4MT_parent
 
G4ParticleDefinition ** G4MT_daughters
 
G4double G4MT_parent_mass
 
G4double * G4MT_daughters_mass
 
G4double * G4MT_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.

68 {idPi=0, idLepton=1, idNutrino=2};

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 }
G4cout
G4GLOB_DLL std::ostream G4cout
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:303
G4endl
#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 }
G4cout
G4GLOB_DLL std::ostream G4cout
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:303
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
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
160  CheckAndFillParent();
161  G4double massK = G4MT_parent->GetPDGMass();
162 
163  // fill daughter particles and their mass
164  CheckAndFillDaughters();
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 }
G4VDecayChannel::CheckAndFillDaughters
void CheckAndFillDaughters()
Definition: G4VDecayChannel.hh:324
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4DecayProducts::PushProducts
G4int PushProducts(G4DynamicParticle *aParticle)
Definition: G4DecayProducts.cc:153
G4VDecayChannel::G4MT_parent
G4ParticleDefinition * G4MT_parent
Definition: G4VDecayChannel.hh:192
G4VDecayChannel::G4MT_daughters
G4ParticleDefinition ** G4MT_daughters
Definition: G4VDecayChannel.hh:193
G4KL3DecayChannel::DalitzDensity
G4double DalitzDensity(G4double parentmass, G4double Epi, G4double El, G4double Enu, G4double massPi, G4double massL, G4double massNu)
Definition: G4KL3DecayChannel.cc:315
rad
static constexpr double rad
Definition: G4SIunits.hh:149
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::setX
void setX(double)
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:311
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:303
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:173
G4KL3DecayChannel::PhaseSpace
void PhaseSpace(G4double Mparent, const G4double *Mdaughter, G4double *Edaughter, G4double *Pdaughter)
Definition: G4KL3DecayChannel.cc:249

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) {
122  kinematics_name = right.kinematics_name;
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
130  ClearDaughtersName();
131 
132  // recreate array
133  numberOfDaughters = right.numberOfDaughters;
134  if ( numberOfDaughters >0 ) {
135  if (daughters_name !=0) ClearDaughtersName();
136  daughters_name = new G4String*[numberOfDaughters];
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 }
G4int
int G4int
Definition: G4Types.hh:78
G4VDecayChannel::kinematics_name
G4String kinematics_name
Definition: G4VDecayChannel.hh:165
G4VDecayChannel::parent_name
G4String * parent_name
Definition: G4VDecayChannel.hh:171
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:173

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;
269  G4double energy;
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 }
G4int
int G4int
Definition: G4Types.hh:78
P
static double P[]
Definition: nf_gammaFunctions.cc:89
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:303
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
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: