Geant4: G4KL3DecayChannel Class Reference

anonymous enum [protected]

Enumerator:

idPi
idLepton
idNutrino

Definition at line 68 of file G4KL3DecayChannel.hh.

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

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

Definition at line 57 of file G4KL3DecayChannel.cc.

References daughterM, G4VDecayChannel::DumpInfo(), G4cout, G4endl, G4KL3DecayChannel(), G4VDecayChannel::GetVerboseLevel(), idLepton, idNutrino, and idPi.

Referenced by G4KL3DecayChannel().

00063                    :G4VDecayChannel("KL3 Decay",theParentName,
00064                                    theBR,  3,
00065                                    thePionName,theLeptonName,theNutrinoName)
00066 {
00067   static const G4String K_plus("kaon+");
00068   static const G4String K_minus("kaon-");
00069   static const G4String K_L("kaon0L");
00070   static const G4String Mu_plus("mu+");
00071   static const G4String Mu_minus("mu-");
00072   static const G4String E_plus("e+");
00073   static const G4String E_minus("e-");
00074   
00075   massK = 0.0;
00076   daughterM[idPi] = 0.0;
00077   daughterM[idLepton] = 0.0;
00078   daughterM[idNutrino] = 0.0;
00079 
00080   // check modes
00081   if ( ((theParentName == K_plus)&&(theLeptonName == E_plus)) ||
00082        ((theParentName == K_minus)&&(theLeptonName == E_minus))   ) {
00083     // K+- (Ke3)
00084     pLambda = 0.0286;
00085     pXi0    = -0.35;
00086    } else if ( ((theParentName == K_plus)&&(theLeptonName == Mu_plus)) ||
00087        ((theParentName == K_minus)&&(theLeptonName == Mu_minus))   ) {
00088     // K+- (Kmu3)
00089     pLambda = 0.033;
00090     pXi0    = -0.35;
00091   } else if ( (theParentName == K_L) && 
00092               ((theLeptonName == E_plus) ||(theLeptonName == E_minus))  ){
00093     // K0L (Ke3)
00094     pLambda = 0.0300;
00095     pXi0    = -0.11;
00096   } else if ( (theParentName == K_L) && 
00097               ((theLeptonName == Mu_plus) ||(theLeptonName == Mu_minus))  ){
00098     // K0L (Kmu3)
00099     pLambda = 0.034;
00100     pXi0    = -0.11;
00101   } else {
00102 #ifdef G4VERBOSE
00103     if (GetVerboseLevel()>2) {
00104       G4cout << "G4KL3DecayChannel:: constructor :";
00105       G4cout << "illegal arguments " << G4endl;;
00106       DumpInfo();
00107     }
00108 #endif
00109     // set values for K0L (Ke3) temporarily
00110     pLambda = 0.0300;
00111     pXi0    = -0.11;
00112   }
00113 }

G4KL3DecayChannel::~G4KL3DecayChannel ( ) [virtual]

Definition at line 115 of file G4KL3DecayChannel.cc.

00116 {
00117 }

G4KL3DecayChannel::G4KL3DecayChannel ( const G4KL3DecayChannel & ) [protected]

Definition at line 119 of file G4KL3DecayChannel.cc.

References daughterM, G4KL3DecayChannel(), idLepton, idNutrino, and idPi.

00119                                                                   :
00120   G4VDecayChannel(right),
00121   massK(right.massK), 
00122   pLambda(right.pLambda), 
00123   pXi0(right.pXi0)
00124 {
00125   daughterM[idPi] = right.daughterM[idPi];
00126   daughterM[idLepton] = right.daughterM[idLepton];
00127   daughterM[idNutrino] = right.daughterM[idNutrino];
00128 }

G4double G4KL3DecayChannel::DalitzDensity	(	G4double	Epi,
		G4double	El,
		G4double	Enu
	)			`[protected]`

Definition at line 327 of file G4KL3DecayChannel.cc.

References daughterM, G4cout, G4endl, G4VDecayChannel::GetVerboseLevel(), idLepton, idNutrino, and idPi.

Referenced by DecayIt().

00328 {
00329   // KL3 decay   Dalitz Plot Density
00330   //               see Chounet et al Phys. Rep. 4, 201
00331   //  arguments
00332   //    Epi: kinetic enregy of pion
00333   //    El:  kinetic enregy of lepton (e or mu)
00334   //    Enu: kinetic energy of nutrino
00335   //  constants
00336   //    pLambda : linear energy dependence of f+
00337   //    pXi0    : = f+(0)/f-
00338   //    pNorm   : normalization factor
00339   //  variables
00340   //    Epi: total energy of pion
00341   //    El:  total energy of lepton (e or mu)
00342   //    Enu: total energy of nutrino
00343 
00344   // mass of daughters
00345   G4double massPi = daughterM[idPi];
00346   G4double massL  = daughterM[idLepton]; 
00347   G4double massNu = daughterM[idNutrino];
00348 
00349   // calcurate total energy
00350   Epi = Epi + massPi;
00351   El  = El  + massL;
00352   Enu = Enu + massNu;
00353   
00354   G4double Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;
00355   G4double E  = Epi_max - Epi;
00356   G4double q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
00357 
00358   G4double F    = 1.0 + pLambda*q2/massPi/massPi;
00359   G4double Fmax = 1.0;
00360   if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0));
00361 
00362   G4double Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi);
00363 
00364   G4double coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
00365   G4double coeffB = massL*massL*(Enu-E/2.0);
00366   G4double coeffC = massL*massL*E/4.0;
00367 
00368   G4double RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
00369 
00370   G4double Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
00371  
00372 #ifdef G4VERBOSE
00373   if (GetVerboseLevel()>2) {
00374     G4cout << "G4KL3DecayChannel::DalitzDensity  " <<G4endl;
00375     G4cout << " Pi[" << massPi/GeV <<"GeV/c/c] :" << Epi/GeV << "GeV" <<G4endl;
00376     G4cout << " L[" << massL/GeV <<"GeV/c/c] :" << El/GeV << "GeV" <<G4endl;
00377     G4cout << " Nu[" << massNu/GeV <<"GeV/c/c] :" << Enu/GeV << "GeV" <<G4endl;
00378     G4cout << " F :" << F  << " Fmax :" << Fmax << "  Xi :" << Xi << G4endl;
00379     G4cout << " A :" << coeffA << "  B :" << coeffB << "  C :"<< coeffC <<G4endl; 
00380     G4cout << " Rho :" << Rho  << "   RhoMax :" << RhoMax << G4endl;
00381   }
00382 #endif
00383   return (Rho/RhoMax);
00384 }

G4DecayProducts * G4KL3DecayChannel::DecayIt ( G4double ) [virtual]

Implements G4VDecayChannel.

Definition at line 164 of file G4KL3DecayChannel.cc.

References DalitzDensity(), daughterM, G4VDecayChannel::daughters, G4VDecayChannel::daughters_name, G4DecayProducts::DumpInfo(), G4VDecayChannel::FillDaughters(), G4VDecayChannel::FillParent(), G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetPDGMass(), G4VDecayChannel::GetVerboseLevel(), idLepton, idNutrino, idPi, G4VDecayChannel::parent, PhaseSpace(), and G4DecayProducts::PushProducts().

00165 {
00166   // this version neglects muon polarization 
00167   //              assumes the pure V-A coupling
00168   //              gives incorrect energy spectrum for Nutrinos
00169 #ifdef G4VERBOSE
00170   if (GetVerboseLevel()>1) G4cout << "G4KL3DecayChannel::DecayIt " << G4endl;
00171 #endif
00172 
00173   // fill parent particle and its mass
00174   if (parent == 0) {
00175     FillParent();
00176   }
00177   massK = parent->GetPDGMass();
00178 
00179   // fill daughter particles and their mass
00180   if (daughters == 0) {
00181     FillDaughters();
00182   }
00183   daughterM[idPi] = daughters[idPi]->GetPDGMass();
00184   daughterM[idLepton] = daughters[idLepton]->GetPDGMass();
00185   daughterM[idNutrino] = daughters[idNutrino]->GetPDGMass();
00186 
00187   // determine momentum/energy of daughters 
00188   //  according to DalitzDensity 
00189   G4double daughterP[3], daughterE[3];
00190   G4double w;
00191   G4double r;
00192   do {
00193     r = G4UniformRand();
00194     PhaseSpace(massK, &daughterM[0], &daughterE[0], &daughterP[0]);
00195     w = DalitzDensity(daughterE[idPi],daughterE[idLepton],daughterE[idNutrino]);
00196   } while ( r > w);
00197 
00198   // output message
00199 #ifdef G4VERBOSE
00200   if (GetVerboseLevel()>1) {
00201     G4cout << *daughters_name[0] << ":" << daughterP[0]/GeV << "[GeV/c]" <<G4endl;
00202     G4cout << *daughters_name[1] << ":" << daughterP[1]/GeV << "[GeV/c]" <<G4endl;
00203     G4cout << *daughters_name[2] << ":" << daughterP[2]/GeV << "[GeV/c]" <<G4endl;
00204   }
00205 #endif
00206    //create parent G4DynamicParticle at rest
00207   G4ThreeVector* direction = new G4ThreeVector(1.0,0.0,0.0);
00208   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, *direction, 0.0);
00209   delete direction;
00210 
00211   //create G4Decayproducts
00212   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00213   delete parentparticle;
00214 
00215   //create daughter G4DynamicParticle 
00216   G4double costheta, sintheta, phi, sinphi, cosphi; 
00217   G4double costhetan, sinthetan, phin, sinphin, cosphin;
00218  
00219   // pion
00220   costheta = 2.*G4UniformRand()-1.0;
00221   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00222   phi  = twopi*G4UniformRand()*rad;
00223   sinphi = std::sin(phi);
00224   cosphi = std::cos(phi);
00225   direction = new G4ThreeVector(sintheta*cosphi,sintheta*sinphi,costheta);
00226   G4ThreeVector momentum0 =  (*direction)*daughterP[0]; 
00227   G4DynamicParticle * daughterparticle 
00228        = new G4DynamicParticle( daughters[0], momentum0);
00229   products->PushProducts(daughterparticle);
00230 
00231   // neutrino
00232   costhetan = (daughterP[1]*daughterP[1]-daughterP[2]*daughterP[2]-daughterP[0]*daughterP[0])/(2.0*daughterP[2]*daughterP[0]);
00233   sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00234   phin  = twopi*G4UniformRand()*rad;
00235   sinphin = std::sin(phin);
00236   cosphin = std::cos(phin);
00237   direction->setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); 
00238   direction->setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); 
00239   direction->setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
00240 
00241   G4ThreeVector momentum2 =  (*direction)*daughterP[2]; 
00242   daughterparticle = new G4DynamicParticle( daughters[2], momentum2);
00243   products->PushProducts(daughterparticle);
00244 
00245   //lepton
00246   G4ThreeVector momentum1 = (momentum0 + momentum2) * (-1.0);
00247   daughterparticle = 
00248        new G4DynamicParticle( daughters[1], momentum1);
00249   products->PushProducts(daughterparticle);
00250 
00251 #ifdef G4VERBOSE
00252   if (GetVerboseLevel()>1) {
00253      G4cout << "G4KL3DecayChannel::DecayIt ";
00254      G4cout << "  create decay products in rest frame " <<G4endl;
00255      G4cout << "  decay products address=" << products << G4endl;
00256      products->DumpInfo();
00257   }
00258 #endif
00259   delete direction;
00260   return products;
00261 }

G4double G4KL3DecayChannel::GetDalitzParameterLambda ( ) const [inline]

Definition at line 114 of file G4KL3DecayChannel.hh.

00115 {
00116   return  pLambda;
00117 }

G4double G4KL3DecayChannel::GetDalitzParameterXi ( ) const [inline]

Definition at line 120 of file G4KL3DecayChannel.hh.

00121 {
00122   return  pXi0;
00123 }

G4KL3DecayChannel & G4KL3DecayChannel::operator= ( const G4KL3DecayChannel & ) [protected]

Definition at line 130 of file G4KL3DecayChannel.cc.

References G4VDecayChannel::ClearDaughtersName(), daughterM, G4VDecayChannel::daughters_name, idLepton, idNutrino, idPi, G4VDecayChannel::kinematics_name, massK, G4VDecayChannel::numberOfDaughters, G4VDecayChannel::parent_name, pLambda, pXi0, G4VDecayChannel::rbranch, and G4VDecayChannel::verboseLevel.

00131 {
00132   if (this != &right) { 
00133     kinematics_name = right.kinematics_name;
00134     verboseLevel = right.verboseLevel;
00135     rbranch = right.rbranch;
00136 
00137     // copy parent name
00138     parent_name = new G4String(*right.parent_name);
00139 
00140     // clear daughters_name array
00141     ClearDaughtersName();
00142 
00143     // recreate array
00144     numberOfDaughters = right.numberOfDaughters;
00145     if ( numberOfDaughters >0 ) {
00146       if (daughters_name !=0) ClearDaughtersName();
00147       daughters_name = new G4String*[numberOfDaughters];
00148       //copy daughters name
00149       for (G4int index=0; index < numberOfDaughters; index++) {
00150           daughters_name[index] = new G4String(*right.daughters_name[index]);
00151       }
00152     }
00153     massK = right.massK; 
00154     pLambda = right.pLambda; 
00155     pXi0 = right.pXi0;
00156     daughterM[idPi] = right.daughterM[idPi];
00157     daughterM[idLepton] = right.daughterM[idLepton];
00158     daughterM[idNutrino] = right.daughterM[idNutrino];
00159   }
00160   return *this;
00161 }

void G4KL3DecayChannel::PhaseSpace	(	G4double	Mparent,
		const G4double *	Mdaughter,
		G4double *	Edaughter,
		G4double *	Pdaughter
	)			`[protected]`

Definition at line 263 of file G4KL3DecayChannel.cc.

References G4cout, G4endl, G4UniformRand, and G4VDecayChannel::GetVerboseLevel().

Referenced by DecayIt().

00268 {
00269   
00270   //sum of daughters'mass
00271   G4double sumofdaughtermass = 0.0;
00272   G4int index;
00273   for (index=0; index<3; index++){
00274     sumofdaughtermass += M[index];
00275   }
00276 
00277   //calculate daughter momentum
00278   //  Generate two 
00279   G4double rd1, rd2, rd;
00280   G4double momentummax=0.0, momentumsum = 0.0;
00281   G4double energy;
00282 
00283   do {
00284     rd1 = G4UniformRand();
00285     rd2 = G4UniformRand();
00286     if (rd2 > rd1) {
00287       rd  = rd1;
00288       rd1 = rd2;
00289       rd2 = rd;
00290     } 
00291     momentummax = 0.0;
00292     momentumsum = 0.0;
00293     // daughter 0
00294     energy = rd2*(parentM - sumofdaughtermass);
00295     P[0] = std::sqrt(energy*energy + 2.0*energy*M[0]);
00296     E[0] = energy;
00297     if ( P[0] >momentummax )momentummax =  P[0];
00298     momentumsum  +=  P[0];
00299     // daughter 1
00300     energy = (1.-rd1)*(parentM - sumofdaughtermass);
00301     P[1] = std::sqrt(energy*energy + 2.0*energy*M[1]);
00302     E[1] = energy;
00303     if ( P[1] >momentummax )momentummax =  P[1];
00304     momentumsum  +=  P[1];
00305     // daughter 2
00306     energy = (rd1-rd2)*(parentM - sumofdaughtermass);
00307     P[2] = std::sqrt(energy*energy + 2.0*energy*M[2]);
00308     E[2] = energy;
00309     if ( P[2] >momentummax )momentummax =  P[2];
00310     momentumsum  +=  P[2];
00311   } while (momentummax >  momentumsum - momentummax );
00312 
00313 #ifdef G4VERBOSE
00314   if (GetVerboseLevel()>2) {
00315      G4cout << "G4KL3DecayChannel::PhaseSpace    ";
00316      G4cout << "Kon mass:" << parentM/GeV << "GeV/c/c" << G4endl;
00317      for (index=0; index<3; index++){
00318        G4cout << index << " : " << M[index]/GeV << "GeV/c/c  ";
00319        G4cout << " : " << E[index]/GeV << "GeV  ";
00320        G4cout << " : " << P[index]/GeV << "GeV/c " << G4endl;
00321      }
00322   }
00323 #endif
00324 }

void G4KL3DecayChannel::SetDalitzParameter	(	G4double	aLambda,
		G4double	aXi
	)			`[inline]`

Definition at line 107 of file G4KL3DecayChannel.hh.

00108 {
00109    pLambda  = aLambda;
00110    pXi0      = aXi;
00111 }

G4double G4KL3DecayChannel::daughterM[3] [protected]

Definition at line 69 of file G4KL3DecayChannel.hh.

Referenced by DalitzDensity(), DecayIt(), G4KL3DecayChannel(), and operator=().


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
Protected Types
	idPi = 0
	idLepton = 1
	idNutrino = 2
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 Epi, G4double El, G4double Enu)
Protected Attributes
G4double	daughterM [3]

G4KL3DecayChannel Class Reference

Public Member Functions

Protected Types

Protected Member Functions

Protected Attributes

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation