CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/Validation/EventGenerator/src/TauDecay.cc

Go to the documentation of this file.
00001 #include "Validation/EventGenerator/interface/TauDecay.h"
00002 #include "Validation/EventGenerator/interface/PdtPdgMini.h"
00003 
00004 
00005 #include <iomanip>
00006 #include <cstdlib> 
00007 #include <iostream>
00008  
00009 TauDecay::TauDecay(){
00010   Reset();
00011 }
00012 
00013 TauDecay::~TauDecay(){
00014 
00015 }
00016 
00017 void TauDecay::Reset(){
00018   n_pi=0;
00019   n_pi0=0;
00020   n_K=0;
00021   n_K0L=0;
00022   n_K0S=0;
00023   n_gamma=0;
00024   n_nu=0;
00025   n_e=0;
00026   n_mu=0;
00027   n_a1=0;
00028   n_a10=0;
00029   n_rho=0;
00030   n_rho0=0;
00031   n_eta=0;
00032   n_omega=0;
00033   n_Kstar0=0;
00034   n_Kstar=0;
00035   unknown=0;
00036 }
00037 
00038 bool TauDecay::isTauFinalStateParticle(int pdgid){
00039   int id=abs(pdgid);
00040   if(id==PdtPdgMini::e_minus)   return true;  // e+-
00041   if(id==PdtPdgMini::nu_e)      return true;  // nu_e
00042   if(id==PdtPdgMini::mu_minus)  return true;  // mu+-
00043   if(id==PdtPdgMini::nu_mu)     return true;  // nu_mu
00044   if(id==PdtPdgMini::nu_tau)    return true;  // nu_tau
00045   if(id==PdtPdgMini::gamma)     return true;  // gamma happends in generator
00046   if(id==PdtPdgMini::pi0)       return true;  // pi0
00047   if(id==PdtPdgMini::pi_plus)   return true;  // pi+-
00048   if(id==PdtPdgMini::K_L0)      return true;  // K0L
00049   if(id==PdtPdgMini::K_S0)      return true;  // KS
00050   if(id==PdtPdgMini::K_plus)    return true;  // K+-
00051   return false;
00052 }
00053 
00054 bool TauDecay::isTauParticleCounter(int pdgid){
00055   int id=abs(pdgid);
00056   //count particles
00057   if(id==PdtPdgMini::pi_plus) { n_pi++;       return true;}
00058   if(id==PdtPdgMini::pi0)     { n_pi0++;      return true;}
00059   if(id==PdtPdgMini::K_plus)  { n_K++;        return true;}
00060   if(id==PdtPdgMini::K_L0)    { n_K0L++;      return true;}
00061   if(id==PdtPdgMini::K_S0)    { n_K0S++;      return true;}
00062   if(id==PdtPdgMini::gamma)   { n_gamma++;    return true;}
00063   if(id==PdtPdgMini::nu_tau ||
00064      id==PdtPdgMini::nu_e   ||
00065      id==PdtPdgMini::nu_mu)    { n_nu++;      return true;}
00066   if(id==PdtPdgMini::e_minus)  { n_e++;       return true;}
00067   if(id==PdtPdgMini::mu_minus) { n_mu++;      return true;}
00068   return false;
00069 }
00070 
00071 bool TauDecay::isTauResonanceCounter(int pdgid){
00072   int id=abs(pdgid);
00073   //count resonances
00074   if(id==PdtPdgMini::a_1_plus)   { n_a1++;      return true;}
00075   if(id==PdtPdgMini::a_10)       { n_a10++;     return true;}
00076   if(id==PdtPdgMini::rho_plus)   { n_rho++;     return true;}
00077   if(id==PdtPdgMini::rho0)       { n_rho0++;    return true;}
00078   if(id==PdtPdgMini::eta)        { n_eta++;     return true;}
00079   if(id==PdtPdgMini::omega)      { n_omega++;   return true;}
00080   //if(id==PdtPdgMini::K_S0)       { n_K0S++;     return true;}
00081   if(id==PdtPdgMini::K_star0)    { n_Kstar0++;  return true;}
00082   if(id==PdtPdgMini::K_star_plus){ n_Kstar++;   return true;}
00083   if(id==PdtPdgMini::W_plus)     { return true;}
00084   unknown++;
00085   return false;
00086 }
00087 
00088 void TauDecay::ClassifyDecayMode(unsigned int &JAK_ID,unsigned int &TauBitMask){
00089   //Reset Bits
00090   JAK_ID=0;
00091   TauBitMask=0;
00092   // Classify according to JAK and TauDecayStructure
00094   //
00095   // Exlusive modes remove first
00096   //
00097   if(n_pi>=1 && n_pi0>=1 && n_nu==1 && n_eta==1){ // eta modes
00098     JAK_ID=JAK_ETAPIPI0;
00099     TauBitMask=OneProng;
00100     if(n_pi0==1)TauBitMask+=OnePi0;
00101     if(n_pi0==2)TauBitMask+=TwoPi0;
00102     if(n_pi0==3)TauBitMask+=ThreePi0;
00103     ClassifyDecayResonance(TauBitMask);
00104     return;
00105   }
00106   //JAK_K0BK0PI
00107   if ((n_K0S+n_K0L==2 && n_pi==1 && n_pi0==0 && n_K==0 && n_nu==1) || (n_pi>=1 && n_pi0==0 && n_K0L+n_K0S==2 && n_K==0 && n_nu==1) ){
00108     JAK_ID=JAK_K0BK0PI;
00109     TauBitMask=OneProng;
00110     if(n_pi>1)TauBitMask+=KS0_to_pipi;
00111     ClassifyDecayResonance(TauBitMask);
00112     return;
00113   }
00114   //if(n_Kstar==1){ //JAK_KSTAR K0SPi && KPi0 
00115     if(n_e==0 && n_mu==0 && n_pi>=1 && n_pi0==0 && n_K==0 && n_K0L+n_K0S==1 && n_nu==1){
00116       JAK_ID=JAK_KSTAR;
00117       TauBitMask=OneProng;
00118       if(n_pi==3)TauBitMask+=KS0_to_pipi;
00119       ClassifyDecayResonance(TauBitMask);
00120       return;
00121     }
00122     if(n_e==0 && n_mu==0 && n_pi==0 && n_pi0==1 && n_K==1 && n_K0L==0 && n_K0S==0 && n_nu==1){
00123       JAK_ID=JAK_KSTAR;
00124       TauBitMask=OneProng;
00125       ClassifyDecayResonance(TauBitMask);
00126       return;
00127     }
00128     //}
00129   //JAK_PIK0PI0
00130   if(n_e==0 && n_mu==0 && n_pi>=1 && n_pi0==1 && n_K==0 && n_K0L+n_K0S==1  && n_nu==1){
00131     JAK_ID=JAK_PIK0PI0;
00132     TauBitMask=OneProng;
00133     TauBitMask+=OnePi0;
00134     if(n_pi==3)TauBitMask+=KS0_to_pipi;
00135     ClassifyDecayResonance(TauBitMask);
00136     return;
00137   }
00138   //JAK_KK0B
00139   if(n_e==0 && n_mu==0 && n_pi0==0 && n_K==1 && n_K0L+n_K0S==1 && n_nu==1){
00140     JAK_ID=JAK_KK0B;
00141     TauBitMask=OneProng;
00142     if(n_pi==2)TauBitMask+=KS0_to_pipi;
00143     return;
00144   }
00145   //JAK_ID=JAK_KK0BPI0
00146   if(n_e==0 && n_mu==0 && n_pi0==1 && n_K==1 && n_K0L+n_K0S==1  && n_nu==1){
00147     JAK_ID=JAK_KK0BPI0;
00148     TauBitMask=OneProng;
00149     TauBitMask+=OnePi0;
00150     if(n_pi==2)TauBitMask+=KS0_to_pipi;
00151     ClassifyDecayResonance(TauBitMask);
00152     return;
00153   }
00154   
00155   
00156   //Safty handelling for exlusive modes
00157   if (n_K0L!=0){
00158     std::cout << "Unknown mode with KL0: n_e " <<  n_e << " n_mu " << n_mu << " n_pi " << n_pi << " n_pi0 " << n_pi0 << " n_K " << n_K << "  n_K0L " << n_K0L << "  n_K0S " << n_K0S << " n_nu  " << n_nu << " n_gamma " << n_gamma 
00159               << std::endl;
00160     return;
00161   }
00162   if (n_K0S!=0){
00163     std::cout << "Unknown mode with KS0: n_e " <<  n_e << " n_mu " << n_mu << " n_pi " << n_pi << " n_pi0 " << n_pi0 << " n_K " << n_K << "  n_K0L " << n_K0L << "  n_K0S " << n_K0S << " n_nu  " << n_nu << " n_gamma " << n_gamma 
00164               << std::endl;
00165     return;
00166   }
00167   if (n_eta!=0){
00168     std::cout << "Unknown mode with eta: n_e " <<  n_e << " n_mu " << n_mu << " n_pi " << n_pi << " n_pi0 " << n_pi0 << " n_K " << n_K << "  n_K0L " << n_K0L << "  n_K0S " << n_K0S << " n_nu  " << n_nu << " n_gamma " << n_gamma 
00169               << std::endl;
00170     return;
00171   }
00172 
00173 
00174 
00175   if(n_pi+n_K+n_e+n_mu==1)TauBitMask=OneProng;
00176   if(n_pi+n_K==3)TauBitMask=ThreeProng;
00177   if(n_pi+n_K==5)TauBitMask=FiveProng;
00178   if(n_pi0==1)TauBitMask+=OnePi0;
00179   if(n_pi0==2)TauBitMask+=TwoPi0;
00180   if(n_pi0==3)TauBitMask+=ThreePi0;
00181   ClassifyDecayResonance(TauBitMask);
00183   //
00184   // Standard modes
00185   //
00186   if(n_e==1 && n_mu==0 && n_pi==0 && n_pi0==0 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==2){
00187     JAK_ID=JAK_ELECTRON;
00188       return;
00189   }
00190   if(n_e==0 && n_mu==1 && n_pi==0 && n_pi0==0 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==2){
00191     JAK_ID=JAK_MUON;
00192     return;
00193   }
00194   if(n_e==0 && n_mu==0 && n_pi==1 && n_pi0==0 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00195     JAK_ID=JAK_PION;
00196     return;
00197   }
00198   if(n_e==0 && n_mu==0 && n_pi==1 && n_pi0==1 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){// && n_rho==1){ removing intermediate resoance to be compatible with pythia8
00199     JAK_ID=JAK_RHO_PIPI0;
00200     return;
00201   }
00202   if(n_e==0 && n_mu==0 && n_pi==1 && n_pi0==2 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00203     JAK_ID=JAK_A1_3PI;
00204     return;
00205   }
00206   if(n_e==0 && n_mu==0 && n_pi==3 && n_pi0==0 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00207     JAK_ID=JAK_A1_3PI;
00208     return;
00209   }
00210   if(n_e==0 && n_mu==0 && n_pi==0 && n_pi0==0 && n_K==1 && n_K0L==0 && n_K0S==0 && n_nu==1){
00211     JAK_ID=JAK_KAON;
00212     return;
00213   }
00214   if(n_e==0 && n_mu==0 && n_pi==3 && n_pi0==1 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00215     JAK_ID=JAK_3PIPI0;
00216     return;
00217   }
00218   if(n_e==0 && n_mu==0 && n_pi==1 && n_pi0==3 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00219     JAK_ID=JAK_PI3PI0;
00220     return;
00221   }
00222   if(n_e==0 && n_mu==0 && n_pi==3 && n_pi0==2 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00223     JAK_ID=JAK_3PI2PI0;
00224     return;
00225   }
00226   if(n_e==0 && n_mu==0 && n_pi==5 && n_pi0==0 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00227     JAK_ID=JAK_5PI;
00228     return;
00229   }
00230   if(n_e==0 && n_mu==0 && n_pi==5 && n_pi0==1 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00231     JAK_ID=JAK_5PIPI0;
00232     return;
00233   }
00234   if(n_e==0 && n_mu==0 && n_pi==3 && n_pi0==3 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1){
00235     JAK_ID=JAK_3PI3PI0;
00236     return;
00237   }
00238   if(n_e==0 && n_mu==0 && n_pi==1 && n_pi0==0 && n_K==2 && n_K0L==0 && n_K0S==0 && n_nu==1){
00239     JAK_ID=JAK_KPIK;
00240     return;
00241   }
00242   if(n_e==0 && n_mu==0 && n_pi==0 && n_pi0==2 && n_K==1 && n_K0L==0 && n_K0S==0 && n_nu==1){
00243     JAK_ID=JAK_K2PI0;
00244     return;
00245   }
00246   if(n_e==0 && n_mu==0 && n_pi==2 && n_pi0==0 && n_K==1 && n_K0L==0 && n_K0S==0 && n_nu==1){
00247     JAK_ID=JAK_KPIPI;
00248     return;
00249   }
00250   // removing JAKID 21 to allow for compatibility with Pythia8 
00251   /*  if(n_e==0 && n_mu==0 && n_pi==1 && n_pi0==1 && n_K==0 && n_K0L==0 && n_K0S==0 && n_nu==1 && n_gamma>=1 && n_rho==0){
00252     JAK_ID=JAK_PIPI0GAM;
00253     return;
00254     }*/
00255   std::cout << "Tau Mode not found: n_e " <<  n_e << " n_mu " << n_mu << " n_pi " << n_pi << " n_pi0 " << n_pi0 << " n_K " << n_K << "  n_K0L " << n_K0L << "  n_K0S " << n_K0S << " n_nu  " << n_nu << " n_gamma " << n_gamma << std::endl;
00256   JAK_ID=JAK_UNKNOWN;
00257 }
00258 
00259 void TauDecay::ClassifyDecayResonance(unsigned int &TauBitMask){
00260   // Add Resonance info to TauBitMask
00261   if(n_a1>0)     TauBitMask+=Res_a1_pm;
00262   if(n_a10>0)    TauBitMask+=Res_a1_0;
00263   if(n_rho>0)    TauBitMask+=Res_rho_pm;
00264   if(n_rho0>0)   TauBitMask+=Res_rho_0;
00265   if(n_eta>0)    TauBitMask+=Res_eta;
00266   if(n_omega>0)  TauBitMask+=Res_omega;
00267   if(n_Kstar>0)  TauBitMask+=Res_Kstar_pm;
00268   if(n_Kstar0>0) TauBitMask+=Res_Kstar_0;
00269 }