49 massK(0.0), pLambda(0.0), pXi0(0.0)
51 daughterM[idPi] = 0.0;
52 daughterM[idLepton] = 0.0;
53 daughterM[idNutrino] = 0.0;
65 thePionName,theLeptonName,theNutrinoName)
67 static const G4String K_plus(
"kaon+");
68 static const G4String K_minus(
"kaon-");
70 static const G4String Mu_plus(
"mu+");
71 static const G4String Mu_minus(
"mu-");
81 if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
82 ((theParentName == K_minus)&&(theLeptonName == E_minus)) ) {
86 }
else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
87 ((theParentName == K_minus)&&(theLeptonName == Mu_minus)) ) {
91 }
else if ( (theParentName == K_L) &&
92 ((theLeptonName == E_plus) ||(theLeptonName == E_minus)) ){
96 }
else if ( (theParentName == K_L) &&
97 ((theLeptonName == Mu_plus) ||(theLeptonName == Mu_minus)) ){
104 G4cout <<
"G4KL3DecayChannel:: constructor :";
122 pLambda(right.pLambda),
132 if (
this != &right) {
154 pLambda = right.pLambda;
189 G4double daughterP[3], daughterE[3];
213 delete parentparticle;
216 G4double costheta, sintheta, phi, sinphi, cosphi;
217 G4double costhetan, sinthetan, phin, sinphin, cosphin;
221 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
223 sinphi = std::sin(phi);
224 cosphi = std::cos(phi);
225 direction =
new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
232 costhetan = (daughterP[1]*daughterP[1]-daughterP[2]*daughterP[2]-daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]);
233 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
235 sinphin = std::sin(phin);
236 cosphin = std::cos(phin);
237 direction->
setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
238 direction->
setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
239 direction->
setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
253 G4cout <<
"G4KL3DecayChannel::DecayIt ";
254 G4cout <<
" create decay products in rest frame " <<
G4endl;
255 G4cout <<
" decay products address=" << products <<
G4endl;
273 for (index=0; index<3; index++){
274 sumofdaughtermass += M[index];
280 G4double momentummax=0.0, momentumsum = 0.0;
294 energy = rd2*(parentM - sumofdaughtermass);
295 P[0] = std::sqrt(energy*energy + 2.0*energy*M[0]);
297 if ( P[0] >momentummax )momentummax = P[0];
300 energy = (1.-rd1)*(parentM - sumofdaughtermass);
301 P[1] = std::sqrt(energy*energy + 2.0*energy*M[1]);
303 if ( P[1] >momentummax )momentummax = P[1];
306 energy = (rd1-rd2)*(parentM - sumofdaughtermass);
307 P[2] = std::sqrt(energy*energy + 2.0*energy*M[2]);
309 if ( P[2] >momentummax )momentummax = P[2];
311 }
while (momentummax > momentumsum - momentummax );
315 G4cout <<
"G4KL3DecayChannel::PhaseSpace ";
317 for (index=0; index<3; index++){
318 G4cout << index <<
" : " << M[index]/
GeV <<
"GeV/c/c ";
319 G4cout <<
" : " << E[index]/
GeV <<
"GeV ";
354 G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;
356 G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
358 G4double F = 1.0 + pLambda*q2/massPi/massPi;
360 if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0));
362 G4double Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi);
364 G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
365 G4double coeffB = massL*massL*(Enu-E/2.0);
366 G4double coeffC = massL*massL*E/4.0;
368 G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
370 G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
378 G4cout <<
" F :" << F <<
" Fmax :" << Fmax <<
" Xi :" << Xi <<
G4endl;
379 G4cout <<
" A :" << coeffA <<
" B :" << coeffB <<
" C :"<< coeffC <<
G4endl;
380 G4cout <<
" Rho :" << Rho <<
" RhoMax :" << RhoMax <<
G4endl;
void ClearDaughtersName()
virtual G4DecayProducts * DecayIt(G4double)
CLHEP::Hep3Vector G4ThreeVector
G4int PushProducts(G4DynamicParticle *aParticle)
G4ParticleDefinition * G4MT_parent
virtual ~G4KL3DecayChannel()
G4ParticleDefinition ** G4MT_daughters
double precision function energy(A, Z)
G4GLOB_DLL std::ostream G4cout
G4double DalitzDensity(G4double Epi, G4double El, G4double Enu)
G4int GetVerboseLevel() const
G4double GetPDGMass() const
G4KL3DecayChannel & operator=(const G4KL3DecayChannel &)
G4String ** daughters_name
void PhaseSpace(G4double Mparent, const G4double *Mdaughter, G4double *Edaughter, G4double *Pdaughter)
G4KL3DecayChannel(const G4String &theParentName, G4double theBR, const G4String &thePionName, const G4String &theLeptonName, const G4String &theNutrinoName)