54 thePionName, theLeptonName, theNutrinoName)
56 static const G4String K_plus(
"kaon+");
57 static const G4String K_minus(
"kaon-");
59 static const G4String Mu_plus(
"mu+");
60 static const G4String Mu_minus(
"mu-");
65 if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
66 ((theParentName == K_minus)&&(theLeptonName == E_minus)) )
72 else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
73 ((theParentName == K_minus)&&(theLeptonName == Mu_minus)) )
79 else if ( (theParentName == K_L) &&
80 ((theLeptonName == E_plus) || (theLeptonName == E_minus)) )
86 else if ( (theParentName == K_L) &&
87 ((theLeptonName == Mu_plus) || (theLeptonName == Mu_minus)) )
98 G4cout <<
"G4KL3DecayChannel:: constructor :";
115 pLambda(right.pLambda),
173 G4double daughterP[3], daughterE[3];
176 const size_t MAX_LOOP = 10000;
177 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
180 PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
209 delete parentparticle;
212 G4double costheta, sintheta, phi, sinphi, cosphi;
213 G4double costhetan, sinthetan, phin, sinphin, cosphin;
217 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
219 sinphi = std::sin(phi);
220 cosphi = std::cos(phi);
221 direction =
new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
228 costhetan = (daughterP[1]*daughterP[1]
229 - daughterP[2]*daughterP[2]
230 - daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]);
231 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
233 sinphin = std::sin(phin);
234 cosphin = std::cos(phin);
235 direction->
setX( sinthetan*cosphin*costheta*cosphi
236 - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
237 direction->
setY( sinthetan*cosphin*costheta*sinphi
238 + 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 const G4int N_DAUGHTER=3;
275 for (index=0; index<N_DAUGHTER; ++index)
277 sumofdaughtermass +=
M[index];
282 G4double momentummax=0.0, momentumsum = 0.0;
284 const size_t MAX_LOOP=10000;
285 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
298 energy = rd2*(parentM - sumofdaughtermass);
301 if (
P[0] >momentummax )momentummax =
P[0];
304 energy = (1.-rd1)*(parentM - sumofdaughtermass);
307 if (
P[1] >momentummax )momentummax =
P[1];
310 energy = (rd1-rd2)*(parentM - sumofdaughtermass);
313 if (
P[2] >momentummax )momentummax =
P[2];
315 if (momentummax <= momentumsum - momentummax )
break;
320 G4cout <<
"G4KL3DecayChannel::PhaseSpace ";
322 for (index=0; index<3; ++index)
324 G4cout << index <<
" : " <<
M[index]/
GeV <<
"GeV/c/c ";
325 G4cout <<
" : " << E[index]/
GeV <<
"GeV ";
356 G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;
358 G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
362 if (
pLambda >0.0) Fmax = (1.0 +
pLambda*(massK*massK/massPi/massPi+1.0));
366 G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
367 G4double coeffB = massL*massL*(Enu-E/2.0);
368 G4double coeffC = massL*massL*E/4.0;
370 G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
372 G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
381 G4cout <<
" F :" << F <<
" Fmax :" << Fmax <<
" Xi :" << Xi <<
G4endl;
382 G4cout <<
" A :" << coeffA <<
" B :" << coeffB <<
" C :"<< coeffC
384 G4cout <<
" Rho :" << Rho <<
" RhoMax :" << RhoMax <<
G4endl;
static constexpr double twopi
static constexpr double rad
static constexpr double GeV
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
G4int PushProducts(G4DynamicParticle *aParticle)
void PhaseSpace(G4double Mparent, const G4double *Mdaughter, G4double *Edaughter, G4double *Pdaughter)
virtual G4DecayProducts * DecayIt(G4double)
G4double DalitzDensity(G4double parentmass, G4double Epi, G4double El, G4double Enu, G4double massPi, G4double massL, G4double massNu)
virtual ~G4KL3DecayChannel()
G4KL3DecayChannel & operator=(const G4KL3DecayChannel &)
G4double GetPDGMass() const
G4ParticleDefinition ** G4MT_daughters
void CheckAndFillParent()
G4String ** daughters_name
G4int GetVerboseLevel() const
G4ParticleDefinition * G4MT_parent
void CheckAndFillDaughters()
void ClearDaughtersName()
G4double energy(const ThreeVector &p, const G4double m)