49 pLambda(0.0), pXi0(0.0)
62 thePionName,theLeptonName,theNutrinoName)
64 static const G4String K_plus(
"kaon+");
65 static const G4String K_minus(
"kaon-");
67 static const G4String Mu_plus(
"mu+");
68 static const G4String Mu_minus(
"mu-");
73 if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
74 ((theParentName == K_minus)&&(theLeptonName == E_minus)) ) {
78 }
else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
79 ((theParentName == K_minus)&&(theLeptonName == Mu_minus)) ) {
83 }
else if ( (theParentName == K_L) &&
84 ((theLeptonName == E_plus) ||(theLeptonName == E_minus)) ){
88 }
else if ( (theParentName == K_L) &&
89 ((theLeptonName == Mu_plus) ||(theLeptonName == Mu_minus)) ){
96 G4cout <<
"G4KL3DecayChannel:: constructor :";
114 pLambda(right.pLambda),
121 if (
this != &right) {
172 G4double daughterP[3], daughterE[3];
175 const size_t MAX_LOOP = 10000;
176 for (
size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
178 PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
180 daughterM[idPi],daughterM[idLepton],daughterM[idNutrino]);
199 delete parentparticle;
202 G4double costheta, sintheta, phi, sinphi, cosphi;
203 G4double costhetan, sinthetan, phin, sinphin, cosphin;
207 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
209 sinphi = std::sin(phi);
210 cosphi = std::cos(phi);
211 direction =
new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
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));
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);
239 G4cout <<
"G4KL3DecayChannel::DecayIt ";
240 G4cout <<
" create decay products in rest frame " <<
G4endl;
241 G4cout <<
" decay products address=" << products <<
G4endl;
259 const G4int N_DAUGHTER=3;
261 for (index=0; index<N_DAUGHTER; index++){
262 sumofdaughtermass += M[index];
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){
282 energy = rd2*(parentM - sumofdaughtermass);
283 P[0] = std::sqrt(energy*energy + 2.0*energy*M[0]);
285 if ( P[0] >momentummax )momentummax = P[0];
288 energy = (1.-rd1)*(parentM - sumofdaughtermass);
289 P[1] = std::sqrt(energy*energy + 2.0*energy*M[1]);
291 if ( P[1] >momentummax )momentummax = P[1];
294 energy = (rd1-rd2)*(parentM - sumofdaughtermass);
295 P[2] = std::sqrt(energy*energy + 2.0*energy*M[2]);
297 if ( P[2] >momentummax )momentummax = P[2];
299 if (momentummax <= momentumsum - momentummax )
break;
303 G4cout <<
"G4KL3DecayChannel::PhaseSpace ";
305 for (index=0; index<3; index++){
306 G4cout << index <<
" : " << M[index]/
GeV <<
"GeV/c/c ";
307 G4cout <<
" : " << E[index]/
GeV <<
"GeV ";
338 G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;
340 G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
344 if (
pLambda >0.0) Fmax = (1.0 +
pLambda*(massK*massK/massPi/massPi+1.0));
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;
352 G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
354 G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
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;
void CheckAndFillDaughters()
void ClearDaughtersName()
virtual G4DecayProducts * DecayIt(G4double)
CLHEP::Hep3Vector G4ThreeVector
G4int PushProducts(G4DynamicParticle *aParticle)
G4ParticleDefinition * G4MT_parent
virtual ~G4KL3DecayChannel()
G4ParticleDefinition ** G4MT_daughters
const G4double w[NPOINTSGL]
G4double DalitzDensity(G4double parentmass, G4double Epi, G4double El, G4double Enu, G4double massPi, G4double massL, G4double massNu)
G4GLOB_DLL std::ostream G4cout
static const double twopi
G4int GetVerboseLevel() const
G4double GetPDGMass() const
G4double energy(const ThreeVector &p, const G4double m)
G4KL3DecayChannel & operator=(const G4KL3DecayChannel &)
G4String ** daughters_name
void PhaseSpace(G4double Mparent, const G4double *Mdaughter, G4double *Edaughter, G4double *Pdaughter)
void CheckAndFillParent()