Reference Manual

Go to the documentation of this file.
00001 #include "HLTriggerOffline/Tau/interface/HLTTauMCProducer.h"
00002 
00003 using namespace edm;
00004 using namespace std;
00005 using namespace reco;
00006 
00007 HLTTauMCProducer::HLTTauMCProducer(const edm::ParameterSet& mc)
00008 {
00009  
00010   //One Parameter Set per Collection
00011 
00012   MC_      = mc.getUntrackedParameter<edm::InputTag>("GenParticles");
00013   ptMinMCTau_ = mc.getUntrackedParameter<double>("ptMinTau",5.);
00014   ptMinMCMuon_ = mc.getUntrackedParameter<double>("ptMinMuon",2.);
00015   ptMinMCElectron_ = mc.getUntrackedParameter<double>("ptMinElectron",5.);
00016   m_PDG_   = mc.getUntrackedParameter<std::vector<int> >("BosonID");
00017   etaMax = mc.getUntrackedParameter<double>("EtaMax",2.5);
00018 
00019   produces<LorentzVectorCollection>("LeptonicTauLeptons");
00020   produces<LorentzVectorCollection>("LeptonicTauElectrons");
00021   produces<LorentzVectorCollection>("LeptonicTauMuons");
00022   produces<LorentzVectorCollection>("HadronicTauOneProng");
00023   produces<LorentzVectorCollection>("HadronicTauThreeProng");
00024   produces<LorentzVectorCollection>("HadronicTauOneAndThreeProng");
00025   produces<LorentzVectorCollection>("TauOther");
00026   produces<LorentzVectorCollection>("Neutrina");
00027   produces<std::vector<int> >("Mothers");
00028 
00029 }
00030 
00031 HLTTauMCProducer::~HLTTauMCProducer(){ }
00032 
00033 void HLTTauMCProducer::produce(edm::Event& iEvent, const edm::EventSetup& iES)
00034 {
00035   //All the code from HLTTauMCInfo is here :-) 
00036   
00037   auto_ptr<LorentzVectorCollection> product_Electrons(new LorentzVectorCollection);
00038   auto_ptr<LorentzVectorCollection> product_Muons(new LorentzVectorCollection);
00039   auto_ptr<LorentzVectorCollection> product_Leptons(new LorentzVectorCollection);
00040   auto_ptr<LorentzVectorCollection> product_OneProng(new LorentzVectorCollection);
00041   auto_ptr<LorentzVectorCollection> product_ThreeProng(new LorentzVectorCollection);
00042   auto_ptr<LorentzVectorCollection> product_OneAndThreeProng(new LorentzVectorCollection);
00043   auto_ptr<LorentzVectorCollection> product_Other(new LorentzVectorCollection);
00044   auto_ptr<LorentzVectorCollection> product_Neutrina(new LorentzVectorCollection);
00045   auto_ptr<std::vector<int> > product_Mothers(new std::vector<int>);
00046   
00047   edm::Handle<GenParticleCollection> genParticles;
00048   iEvent.getByLabel(MC_, genParticles);
00049 
00050   GenParticleCollection::const_iterator p = genParticles->begin();
00051   
00052   for (;p != genParticles->end(); ++p ) {
00053     //Check the PDG ID
00054     bool pdg_ok = true;
00055     for(size_t pi =0;pi<m_PDG_.size();++pi)
00056       {
00057         if(abs((*p).pdgId())== m_PDG_[pi] && abs((*p).status()) == 3 ){
00058           pdg_ok = true;
00059           //      cout<<" Bsoson particles: "<< (*p).pdgId()<< " " <<(*p).status() << " "<< pdg_ok<<endl;
00060         }
00061       }
00062     
00063     // Check if the boson is one of interest and if there is a valid vertex
00064     if(  pdg_ok )
00065       {
00066         product_Mothers->push_back((*p).pdgId());
00067 
00068         std::vector<GenParticle*> decayProducts;
00069         
00070         TLorentzVector Boson((*p).px(),(*p).py(),(*p).pz(),(*p).energy());      
00071         
00072         for (GenParticle::const_iterator BosonIt=(*p).begin(); BosonIt != (*p).end(); BosonIt++){
00073           // cout<<" Dparticles: "<< (*BosonIt).pdgId() << " "<< (*BosonIt).status()<<endl;
00074           
00075           if (abs((*BosonIt).pdgId()) == 15 && ((*BosonIt).status()==3)) //if it is a Tau and decayed
00076             {
00077               decayProducts.clear();      
00078               //              cout<<" Boson daugther particles: "<< (*BosonIt).pdgId() << " "<< (*BosonIt).status()<< endl;           
00079               for (GenParticle::const_iterator TauIt = (*BosonIt).begin(); TauIt != (*BosonIt).end(); TauIt++) {
00080                 //      cout<<" Tau daughter particles: "<< (*TauIt).pdgId() << " "<< (*TauIt).status()<<endl;
00081 
00082                 if (abs((*TauIt).pdgId()) == 15 && ((*TauIt).status()==2)) //if it is a Tau and decayed
00083                   {                 
00084                     decayProducts = getGenStableDecayProducts((reco::GenParticle*) & (*TauIt));  
00085                     //  for (GenParticle::const_iterator TauIt2 = (*TauIt).begin(); TauIt2 != (*TauIt).end(); TauIt2++) {
00086                     //                cout<<" Real Tau particles: "<< (*TauIt2).pdgId() << " "<< (*TauIt2).status()<< " mother: "<< (*TauIt2).mother()->pdgId() << endl;
00087                     // }
00088                   }
00089               }
00090               
00091               
00092               if ( !decayProducts.empty() )
00093                 {
00094                   
00095                   LorentzVector Visible_Taus(0.,0.,0.,0.);
00096                   LorentzVector TauDecayProduct(0.,0.,0.,0.);
00097                   LorentzVector Neutrino(0.,0.,0.,0.);
00098                   
00099                   int numElectrons      = 0;
00100                   int numMuons          = 0;
00101                   int numChargedPions   = 0;
00102                   int numNeutralPions   = 0;
00103                   int numNeutrinos      = 0;
00104                   int numOtherParticles = 0;
00105                   
00106                   
00107                   for (std::vector<GenParticle*>::iterator pit = decayProducts.begin(); pit != decayProducts.end(); pit++)
00108                     {
00109                       int pdg_id = abs((*pit)->pdgId());
00110                       if (pdg_id == 11) numElectrons++;
00111                       else if (pdg_id == 13) numMuons++;
00112                       else if (pdg_id == 211) numChargedPions++;
00113                       else if (pdg_id == 111) numNeutralPions++;
00114                       else if (pdg_id == 12 || 
00115                                pdg_id == 14 || 
00116                                pdg_id == 16) {
00117                         numNeutrinos++;
00118                         if (pdg_id == 16) {
00119                           Neutrino.SetPxPyPzE((*pit)->px(),(*pit)->py(),(*pit)->pz(),(*pit)->energy());
00120                         }
00121                       }
00122                       else if (pdg_id != 22) {
00123                         numOtherParticles++;
00124                       }
00125                       
00126                       if (pdg_id != 12 &&
00127                           pdg_id != 14 && 
00128                           pdg_id != 16){
00129                         TauDecayProduct.SetPxPyPzE((*pit)->px(),(*pit)->py(),(*pit)->pz(),(*pit)->energy());
00130                         Visible_Taus+=TauDecayProduct;
00131                       } 
00132                       //                  cout<< "This has to be the same: " << (*pit)->pdgId() << " "<< (*pit)->status()<< " mother: "<< (*pit)->mother()->pdgId() << endl;
00133                     }
00134                   
00135                   int tauDecayMode = kOther;
00136                   
00137                   if ( numOtherParticles == 0 ){
00138                     if ( numElectrons == 1 ){
00139                       //--- tau decays into electrons
00140                       tauDecayMode = kElectron;
00141                     } else if ( numMuons == 1 ){
00142                       //--- tau decays into muons
00143                       tauDecayMode = kMuon;
00144                     } else {
00145                       //--- hadronic tau decays
00146                       switch ( numChargedPions ){
00147                       case 1 : 
00148                         switch ( numNeutralPions ){
00149                         case 0 : 
00150                           tauDecayMode = kOneProng0pi0;
00151                           break;
00152                         case 1 : 
00153                           tauDecayMode = kOneProng1pi0;
00154                           break;
00155                         case 2 : 
00156                           tauDecayMode = kOneProng2pi0;
00157                           break;
00158                         }
00159                         break;
00160                       case 3 : 
00161                         switch ( numNeutralPions ){
00162                         case 0 : 
00163                           tauDecayMode = kThreeProng0pi0;
00164                           break;
00165                         case 1 : 
00166                           tauDecayMode = kThreeProng1pi0;
00167                           break;
00168                         }
00169                         break;
00170                       }             
00171                     }
00172                   }
00173                   
00174                   
00175                   //              cout<< "So we have a: " << tauDecayMode <<endl;
00176                   
00177                   if(tauDecayMode == kElectron)
00178                     {                     
00179                       if((abs(Visible_Taus.eta())<etaMax)&&(Visible_Taus.pt()>ptMinMCElectron_)){
00180                         product_Electrons->push_back(Visible_Taus);
00181                         product_Leptons->push_back(Visible_Taus);
00182                       }
00183                     }
00184                   else if (tauDecayMode == kMuon)
00185                     {
00186                       
00187                       if((abs(Visible_Taus.eta())<etaMax)&&(Visible_Taus.pt()>ptMinMCMuon_)){
00188                         product_Muons->push_back(Visible_Taus);
00189                         product_Leptons->push_back(Visible_Taus);
00190                       }
00191                     }
00192                   else if(tauDecayMode == kOneProng0pi0 || 
00193                           tauDecayMode == kOneProng1pi0 || 
00194                           tauDecayMode == kOneProng2pi0 ) 
00195                     {
00196                       if ((abs(Visible_Taus.eta()) < etaMax) && (Visible_Taus.pt() > ptMinMCTau_)){
00197                         product_OneProng->push_back(Visible_Taus);
00198                         product_OneAndThreeProng->push_back(Visible_Taus);
00199                         product_Neutrina->push_back(Neutrino);
00200                       }
00201                     }
00202                   else if (tauDecayMode == kThreeProng0pi0 || 
00203                            tauDecayMode == kThreeProng1pi0 )
00204                     {
00205                       if((abs(Visible_Taus.eta())<etaMax)&&(Visible_Taus.pt()>ptMinMCTau_))  {
00206                         product_ThreeProng->push_back(Visible_Taus);
00207                         product_OneAndThreeProng->push_back(Visible_Taus);
00208                         product_Neutrina->push_back(Neutrino);
00209                       }                                                         
00210                     }
00211                   else if (tauDecayMode == kOther)
00212                     {
00213                       if((abs(Visible_Taus.eta())<etaMax)&&(Visible_Taus.pt()>ptMinMCTau_))  {
00214                         product_Other->push_back(Visible_Taus);
00215                       }
00216                     }                                          
00217                 }
00218             }
00219         }                              
00220         //  
00221       }
00222   }
00223   iEvent.put(product_Leptons,"LeptonicTauLeptons");
00224   iEvent.put(product_Electrons,"LeptonicTauElectrons");
00225   iEvent.put(product_Muons,"LeptonicTauMuons");
00226   iEvent.put(product_OneProng,"HadronicTauOneProng");
00227   iEvent.put(product_ThreeProng,"HadronicTauThreeProng");
00228   iEvent.put(product_OneAndThreeProng,"HadronicTauOneAndThreeProng");
00229   iEvent.put(product_Other, "TauOther");
00230   iEvent.put(product_Neutrina,"Neutrina"); 
00231   iEvent.put(product_Mothers,"Mothers"); 
00232   
00233                                                        
00234 }
00235 // Helper Function
00236 
00237 std::vector<reco::GenParticle*> HLTTauMCProducer::getGenStableDecayProducts(const reco::GenParticle* particle)
00238 {
00239   std::vector<GenParticle*> decayProducts;
00240   decayProducts.clear();
00241 
00242   //  std::cout << " Are we ever here?: "<< (*particle).numberOfDaughters() << std::endl;
00243   for ( GenParticle::const_iterator daughter_particle = (*particle).begin();daughter_particle != (*particle).end(); ++daughter_particle ){   
00244 
00245     int pdg_id = abs((*daughter_particle).pdgId());
00246 
00247 //    // check if particle is stable
00248     if ( pdg_id == 11 || pdg_id == 12 || pdg_id == 13 || pdg_id == 14 || pdg_id == 16 ||  pdg_id == 111 || pdg_id == 211 ){
00249       // stable particle, identified by pdg code
00250       decayProducts.push_back((reco::GenParticle*) &(* daughter_particle));
00251     } 
00252     else {
00253 //      // unstable particle, identified by non-zero decay vertex
00254       std::vector<GenParticle*> addDecayProducts = getGenStableDecayProducts((reco::GenParticle*) &(*daughter_particle));
00255       for ( std::vector<GenParticle*>::const_iterator adddaughter_particle = addDecayProducts.begin(); adddaughter_particle != addDecayProducts.end(); ++adddaughter_particle ){
00256         decayProducts.push_back((*adddaughter_particle));
00257       }
00258     }
00259   }
00260   return decayProducts;
00261 }
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/HLTriggerOffline/Tau/src/HLTTauMCProducer.cc
