CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/DQMOffline/Trigger/src/HLTTauRefProducer.cc

Go to the documentation of this file.
00001 #include "DQMOffline/Trigger/interface/HLTTauRefProducer.h"
00002 #include "TLorentzVector.h"
00003 // TAU includes
00004 #include "DataFormats/TauReco/interface/PFTau.h"
00005 #include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
00006 #include "DataFormats/TauReco/interface/CaloTau.h"
00007 #include "DataFormats/TauReco/interface/CaloTauDiscriminatorByIsolation.h"
00008 // ELECTRON includes
00009 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00010 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00011 #include "DataFormats/EgammaCandidates/interface/Electron.h"
00012 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
00013 #include "AnalysisDataFormats/Egamma/interface/ElectronIDAssociation.h"
00014 #include "AnalysisDataFormats/Egamma/interface/ElectronID.h"
00015 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00016 #include "DataFormats/TrackReco/interface/Track.h"
00017 // MUON includes
00018 #include "DataFormats/MuonReco/interface/Muon.h"
00019 #include "DataFormats/JetReco/interface/CaloJet.h"
00020 #include "TLorentzVector.h"
00021 #include "DataFormats/Math/interface/deltaR.h"
00022 //CaloTower includes
00023 #include "DataFormats/CaloTowers/interface/CaloTower.h"
00024 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00025 #include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
00026 #include "Math/GenVector/VectorUtil.h"
00027 
00028 
00029 using namespace edm;
00030 using namespace reco;
00031 using namespace std;
00032 
00033 HLTTauRefProducer::HLTTauRefProducer(const edm::ParameterSet& iConfig)
00034 {
00035 
00036  
00037   //One Parameter Set per Collection
00038 
00039   ParameterSet pfTau = iConfig.getUntrackedParameter<edm::ParameterSet>("PFTaus");
00040   PFTaus_ = pfTau.getUntrackedParameter<InputTag>("PFTauProducer");
00041   PFTauDis_ = pfTau.getUntrackedParameter<std::vector<InputTag> >("PFTauDiscriminators");
00042   doPFTaus_ = pfTau.getUntrackedParameter<bool>("doPFTaus",false);
00043   ptMinPFTau_= pfTau.getUntrackedParameter<double>("ptMin",15.);
00044 
00045   ParameterSet  electrons = iConfig.getUntrackedParameter<edm::ParameterSet>("Electrons");
00046   Electrons_ = electrons.getUntrackedParameter<InputTag>("ElectronCollection");
00047   doElectrons_ = electrons.getUntrackedParameter<bool>("doElectrons",false);
00048   e_idAssocProd_ = electrons.getUntrackedParameter<InputTag>("IdCollection");
00049   e_ctfTrackCollection_= electrons.getUntrackedParameter<InputTag>("TrackCollection");
00050   ptMinElectron_= electrons.getUntrackedParameter<double>("ptMin",15.);
00051   e_doID_ = electrons.getUntrackedParameter<bool>("doID",false);
00052   e_doTrackIso_ = electrons.getUntrackedParameter<bool>("doTrackIso",false);
00053   e_trackMinPt_= electrons.getUntrackedParameter<double>("ptMinTrack",1.5);
00054   e_lipCut_= electrons.getUntrackedParameter<double>("lipMinTrack",1.5);
00055   e_minIsoDR_= electrons.getUntrackedParameter<double>("InnerConeDR",0.02);
00056   e_maxIsoDR_= electrons.getUntrackedParameter<double>("OuterConeDR",0.6);
00057   e_isoMaxSumPt_= electrons.getUntrackedParameter<double>("MaxIsoVar",0.02);
00058 
00059   ParameterSet  muons = iConfig.getUntrackedParameter<edm::ParameterSet>("Muons");
00060   Muons_ = muons.getUntrackedParameter<InputTag>("MuonCollection");
00061   doMuons_ = muons.getUntrackedParameter<bool>("doMuons",false);
00062   ptMinMuon_= muons.getUntrackedParameter<double>("ptMin",15.);
00063 
00064   ParameterSet  jets = iConfig.getUntrackedParameter<edm::ParameterSet>("Jets");
00065   Jets_ = jets.getUntrackedParameter<InputTag>("JetCollection");
00066   doJets_ = jets.getUntrackedParameter<bool>("doJets");
00067   ptMinJet_= jets.getUntrackedParameter<double>("etMin");
00068 
00069   ParameterSet  towers = iConfig.getUntrackedParameter<edm::ParameterSet>("Towers");
00070   Towers_ = towers.getUntrackedParameter<InputTag>("TowerCollection");
00071   doTowers_ = towers.getUntrackedParameter<bool>("doTowers");
00072   ptMinTower_= towers.getUntrackedParameter<double>("etMin");
00073   towerIsol_= towers.getUntrackedParameter<double>("towerIsolation");
00074 
00075   ParameterSet  photons = iConfig.getUntrackedParameter<edm::ParameterSet>("Photons");
00076   Photons_ = photons.getUntrackedParameter<InputTag>("PhotonCollection");
00077   doPhotons_ = photons.getUntrackedParameter<bool>("doPhotons");
00078   ptMinPhoton_= photons.getUntrackedParameter<double>("etMin");
00079   photonEcalIso_= photons.getUntrackedParameter<double>("ECALIso");
00080 
00081 
00082   etaMax = iConfig.getUntrackedParameter<double>("EtaMax",2.5);
00083   
00084 
00085   //recoCollections
00086   produces<LorentzVectorCollection>("PFTaus");
00087   produces<LorentzVectorCollection>("Electrons");
00088   produces<LorentzVectorCollection>("Muons");
00089   produces<LorentzVectorCollection>("Jets");
00090   produces<LorentzVectorCollection>("Photons");
00091   produces<LorentzVectorCollection>("Towers");
00092 
00093 }
00094 
00095 HLTTauRefProducer::~HLTTauRefProducer(){ }
00096 
00097 void HLTTauRefProducer::produce(edm::Event& iEvent, const edm::EventSetup& iES)
00098 {
00099   if(doPFTaus_)
00100     doPFTaus(iEvent,iES);
00101   if(doElectrons_)
00102     doElectrons(iEvent,iES);
00103   if(doMuons_)
00104     doMuons(iEvent,iES);
00105   if(doJets_)
00106     doJets(iEvent,iES);
00107   if(doPhotons_)
00108     doPhotons(iEvent,iES);
00109   if(doTowers_)
00110     doTowers(iEvent,iES);
00111 
00112 }
00113 
00114 void 
00115 HLTTauRefProducer::doPFTaus(edm::Event& iEvent,const edm::EventSetup& iES)
00116 {
00117       auto_ptr<LorentzVectorCollection> product_PFTaus(new LorentzVectorCollection);
00118       //Retrieve the collection
00119       edm::Handle<PFTauCollection> pftaus;
00120       if(iEvent.getByLabel(PFTaus_,pftaus))
00121         {
00122           for(unsigned int i=0;i<pftaus->size();++i)
00123             if((*pftaus)[i].pt()>ptMinPFTau_&&fabs((*pftaus)[i].eta())<etaMax)
00124               {
00125                 reco::PFTauRef thePFTau(pftaus,i);
00126                 for(unsigned int j=0;j<PFTauDis_.size();++j)
00127                 {
00128                   edm::Handle<PFTauDiscriminator> pftaudis;
00129                   if(iEvent.getByLabel(PFTauDis_[j],pftaudis))
00130                     {
00131                       if((*pftaudis)[thePFTau]>0.5)
00132                         {
00133                           LorentzVector vec((*pftaus)[i].px(),(*pftaus)[i].py(),(*pftaus)[i].pz(),(*pftaus)[i].energy());
00134                           product_PFTaus->push_back(vec);
00135                         }
00136                     }
00137                 }
00138               }
00139         }
00140 
00141       iEvent.put(product_PFTaus,"PFTaus");
00142       
00143 }
00144 
00145 
00146 void 
00147 HLTTauRefProducer::doElectrons(edm::Event& iEvent,const edm::EventSetup& iES)
00148 {
00149   auto_ptr<LorentzVectorCollection> product_Electrons(new LorentzVectorCollection);
00150   //Retrieve the collections
00151   
00152   edm::Handle<reco::ElectronIDAssociationCollection> pEleID;
00153   if(e_doID_){//UGLY HACK UNTIL GET ELETRON ID WORKING IN 210
00154    
00155       iEvent.getByLabel(e_idAssocProd_,pEleID);
00156     
00157     if (!pEleID.isValid()){
00158       edm::LogInfo("")<< "Error! Can't get electronIDAssocProducer by label. ";
00159       e_doID_ = false;
00160     }
00161   }
00162   edm::Handle<reco::TrackCollection> pCtfTracks;
00163   iEvent.getByLabel(e_ctfTrackCollection_, pCtfTracks);
00164   if (!pCtfTracks.isValid()) {
00165   edm::LogInfo("")<< "Error! Can't get " << e_ctfTrackCollection_.label() << " by label. ";
00166   iEvent.put(product_Electrons,"Electrons"); 
00167   return;
00168   }
00169   const reco::TrackCollection * ctfTracks = pCtfTracks.product();
00170   edm::Handle<GsfElectronCollection> electrons;
00171   if(iEvent.getByLabel(Electrons_,electrons))
00172     for(size_t i=0;i<electrons->size();++i)
00173       {
00174         edm::Ref<reco::GsfElectronCollection> electronRef(electrons,i);
00175         bool idDec=false;
00176         if(e_doID_){
00177           reco::ElectronIDAssociationCollection::const_iterator tagIDAssocItr;
00178           tagIDAssocItr = pEleID->find(electronRef);
00179           const reco::ElectronIDRef& id_tag = tagIDAssocItr->val;
00180           idDec=id_tag->cutBasedDecision();
00181         }else idDec=true;
00182         if((*electrons)[i].pt()>ptMinElectron_&&fabs((*electrons)[i].eta())<etaMax&&idDec)
00183           {
00184             if(e_doTrackIso_){
00185               reco::TrackCollection::const_iterator tr = ctfTracks->begin();
00186               double sum_of_pt_ele=0;
00187               for(;tr != ctfTracks->end();++tr)
00188                 {
00189                   double lip = (*electrons)[i].gsfTrack()->dz() - tr->dz();
00190                   if(tr->pt() > e_trackMinPt_ && fabs(lip) < e_lipCut_){
00191                     double dphi=fabs(tr->phi()-(*electrons)[i].trackMomentumAtVtx().phi());
00192                     if(dphi>acos(-1.))dphi=2*acos(-1.)-dphi;
00193                     double deta=fabs(tr->eta()-(*electrons)[i].trackMomentumAtVtx().eta());
00194                     double dr_ctf_ele = sqrt(deta*deta+dphi*dphi);
00195                     if((dr_ctf_ele>e_minIsoDR_) && (dr_ctf_ele<e_maxIsoDR_)){
00196                       double cft_pt_2 = (tr->pt())*(tr->pt());
00197                       sum_of_pt_ele += cft_pt_2;
00198                     }
00199                   }
00200                 }
00201               double isolation_value_ele = sum_of_pt_ele/((*electrons)[i].trackMomentumAtVtx().Rho()*(*electrons)[i].trackMomentumAtVtx().Rho());
00202               if(isolation_value_ele<e_isoMaxSumPt_){
00203                 LorentzVector vec((*electrons)[i].px(),(*electrons)[i].py(),(*electrons)[i].pz(),(*electrons)[i].energy());
00204                 product_Electrons->push_back(vec);
00205               } 
00206               
00207             }
00208             else{ 
00209               LorentzVector vec((*electrons)[i].px(),(*electrons)[i].py(),(*electrons)[i].pz(),(*electrons)[i].energy());
00210               product_Electrons->push_back(vec);
00211             }
00212           }
00213       }
00214   
00215   iEvent.put(product_Electrons,"Electrons"); 
00216 }
00217 
00218 void 
00219 HLTTauRefProducer::doMuons(edm::Event& iEvent,const edm::EventSetup& iES)
00220 {
00221   auto_ptr<LorentzVectorCollection> product_Muons(new LorentzVectorCollection);
00222   //Retrieve the collection
00223   edm::Handle<MuonCollection> muons;
00224       if(iEvent.getByLabel(Muons_,muons))
00225    
00226       for(size_t i = 0 ;i<muons->size();++i)
00227         {
00228          
00229             if((*muons)[i].pt()>ptMinMuon_&&fabs((*muons)[i].eta())<etaMax)
00230               {
00231                 LorentzVector vec((*muons)[i].px(),(*muons)[i].py(),(*muons)[i].pz(),(*muons)[i].energy());
00232                 product_Muons->push_back(vec);
00233               }
00234         }
00235 
00236 
00237       iEvent.put(product_Muons,"Muons");
00238  
00239 }
00240 
00241 
00242 void 
00243 HLTTauRefProducer::doJets(edm::Event& iEvent,const edm::EventSetup& iES)
00244 {
00245       auto_ptr<LorentzVectorCollection> product_Jets(new LorentzVectorCollection);
00246       //Retrieve the collection
00247       edm::Handle<CaloJetCollection> jets;
00248       if(iEvent.getByLabel(Jets_,jets))
00249       for(size_t i = 0 ;i<jets->size();++i)
00250         {
00251              if((*jets)[i].et()>ptMinJet_&&fabs((*jets)[i].eta())<etaMax)
00252               {
00253                 LorentzVector vec((*jets)[i].px(),(*jets)[i].py(),(*jets)[i].pz(),(*jets)[i].energy());
00254                 product_Jets->push_back(vec);
00255               }
00256         }
00257       iEvent.put(product_Jets,"Jets");
00258 }
00259 
00260 void 
00261 HLTTauRefProducer::doTowers(edm::Event& iEvent,const edm::EventSetup& iES)
00262 {
00263       auto_ptr<LorentzVectorCollection> product_Towers(new LorentzVectorCollection);
00264       //Retrieve the collection
00265       edm::Handle<CaloTowerCollection> towers;
00266       if(iEvent.getByLabel(Towers_,towers))
00267       for(size_t i = 0 ;i<towers->size();++i)
00268         {
00269              if((*towers)[i].pt()>ptMinTower_&&fabs((*towers)[i].eta())<etaMax)
00270               {
00271                 //calculate isolation
00272                 double isolET=0;
00273                 for(unsigned int j=0;j<towers->size();++j)
00274                   {
00275                     if(ROOT::Math::VectorUtil::DeltaR((*towers)[i].p4(),(*towers)[j].p4())<0.5)
00276                       isolET+=(*towers)[j].pt();
00277                   }
00278                 isolET-=(*towers)[i].pt();
00279                 if(isolET<towerIsol_)
00280                   {
00281                     LorentzVector vec((*towers)[i].px(),(*towers)[i].py(),(*towers)[i].pz(),(*towers)[i].energy());
00282                     product_Towers->push_back(vec);
00283                   }
00284               }
00285         }
00286       iEvent.put(product_Towers,"Towers");
00287 }
00288 
00289 
00290 void 
00291 HLTTauRefProducer::doPhotons(edm::Event& iEvent,const edm::EventSetup& iES)
00292 {
00293       auto_ptr<LorentzVectorCollection> product_Gammas(new LorentzVectorCollection);
00294       //Retrieve the collection
00295       edm::Handle<PhotonCollection> photons;
00296       if(iEvent.getByLabel(Photons_,photons))
00297       for(size_t i = 0 ;i<photons->size();++i)
00298         if((*photons)[i].ecalRecHitSumEtConeDR04()<photonEcalIso_)
00299         {
00300              if((*photons)[i].et()>ptMinPhoton_&&fabs((*photons)[i].eta())<etaMax)
00301               {
00302                 LorentzVector vec((*photons)[i].px(),(*photons)[i].py(),(*photons)[i].pz(),(*photons)[i].energy());
00303                 product_Gammas->push_back(vec);
00304               }
00305         }
00306       iEvent.put(product_Gammas,"Photons");
00307 }
00308