HLTTauRefProducer.cc

Go to the documentation of this file.

#include "DQMOffline/Trigger/interface/HLTTauRefProducer.h"
#include "TLorentzVector.h"
// TAU includes
#include "DataFormats/TauReco/interface/PFTau.h"
#include "DataFormats/TauReco/interface/PFTauDiscriminator.h"
#include "DataFormats/TauReco/interface/CaloTau.h"
#include "DataFormats/TauReco/interface/CaloTauDiscriminatorByIsolation.h"
// ELECTRON includes
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
#include "DataFormats/EgammaCandidates/interface/Electron.h"
#include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
#include "AnalysisDataFormats/Egamma/interface/ElectronIDAssociation.h"
#include "AnalysisDataFormats/Egamma/interface/ElectronID.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/TrackReco/interface/Track.h"
// MUON includes
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "TLorentzVector.h"
#include "DataFormats/Math/interface/deltaR.h"
//CaloTower includes
#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
#include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
#include "Math/GenVector/VectorUtil.h"


using namespace edm;
using namespace reco;
using namespace std;

HLTTauRefProducer::HLTTauRefProducer(const edm::ParameterSet& iConfig)
{

 
  //One Parameter Set per Collection

  ParameterSet pfTau = iConfig.getUntrackedParameter<edm::ParameterSet>("PFTaus");
  PFTaus_ = pfTau.getUntrackedParameter<InputTag>("PFTauProducer");
  PFTauDis_ = pfTau.getUntrackedParameter<std::vector<InputTag> >("PFTauDiscriminators");
  doPFTaus_ = pfTau.getUntrackedParameter<bool>("doPFTaus",false);
  ptMinPFTau_= pfTau.getUntrackedParameter<double>("ptMin",15.);

  ParameterSet  electrons = iConfig.getUntrackedParameter<edm::ParameterSet>("Electrons");
  Electrons_ = electrons.getUntrackedParameter<InputTag>("ElectronCollection");
  doElectrons_ = electrons.getUntrackedParameter<bool>("doElectrons",false);
  e_idAssocProd_ = electrons.getUntrackedParameter<InputTag>("IdCollection");
  e_ctfTrackCollection_= electrons.getUntrackedParameter<InputTag>("TrackCollection");
  ptMinElectron_= electrons.getUntrackedParameter<double>("ptMin",15.);
  e_doID_ = electrons.getUntrackedParameter<bool>("doID",false);
  e_doTrackIso_ = electrons.getUntrackedParameter<bool>("doTrackIso",false);
  e_trackMinPt_= electrons.getUntrackedParameter<double>("ptMinTrack",1.5);
  e_lipCut_= electrons.getUntrackedParameter<double>("lipMinTrack",1.5);
  e_minIsoDR_= electrons.getUntrackedParameter<double>("InnerConeDR",0.02);
  e_maxIsoDR_= electrons.getUntrackedParameter<double>("OuterConeDR",0.6);
  e_isoMaxSumPt_= electrons.getUntrackedParameter<double>("MaxIsoVar",0.02);

  ParameterSet  muons = iConfig.getUntrackedParameter<edm::ParameterSet>("Muons");
  Muons_ = muons.getUntrackedParameter<InputTag>("MuonCollection");
  doMuons_ = muons.getUntrackedParameter<bool>("doMuons",false);
  ptMinMuon_= muons.getUntrackedParameter<double>("ptMin",15.);

  ParameterSet  jets = iConfig.getUntrackedParameter<edm::ParameterSet>("Jets");
  Jets_ = jets.getUntrackedParameter<InputTag>("JetCollection");
  doJets_ = jets.getUntrackedParameter<bool>("doJets");
  ptMinJet_= jets.getUntrackedParameter<double>("etMin");

  ParameterSet  towers = iConfig.getUntrackedParameter<edm::ParameterSet>("Towers");
  Towers_ = towers.getUntrackedParameter<InputTag>("TowerCollection");
  doTowers_ = towers.getUntrackedParameter<bool>("doTowers");
  ptMinTower_= towers.getUntrackedParameter<double>("etMin");
  towerIsol_= towers.getUntrackedParameter<double>("towerIsolation");

  ParameterSet  photons = iConfig.getUntrackedParameter<edm::ParameterSet>("Photons");
  Photons_ = photons.getUntrackedParameter<InputTag>("PhotonCollection");
  doPhotons_ = photons.getUntrackedParameter<bool>("doPhotons");
  ptMinPhoton_= photons.getUntrackedParameter<double>("etMin");
  photonEcalIso_= photons.getUntrackedParameter<double>("ECALIso");


  etaMax = iConfig.getUntrackedParameter<double>("EtaMax",2.5);
  

  //recoCollections
  produces<LorentzVectorCollection>("PFTaus");
  produces<LorentzVectorCollection>("Electrons");
  produces<LorentzVectorCollection>("Muons");
  produces<LorentzVectorCollection>("Jets");
  produces<LorentzVectorCollection>("Photons");
  produces<LorentzVectorCollection>("Towers");

}

HLTTauRefProducer::~HLTTauRefProducer(){ }

void HLTTauRefProducer::produce(edm::Event& iEvent, const edm::EventSetup& iES)
{
  if(doPFTaus_)
    doPFTaus(iEvent,iES);
  if(doElectrons_)
    doElectrons(iEvent,iES);
  if(doMuons_)
    doMuons(iEvent,iES);
  if(doJets_)
    doJets(iEvent,iES);
  if(doPhotons_)
    doPhotons(iEvent,iES);
  if(doTowers_)
    doTowers(iEvent,iES);

}

void 
HLTTauRefProducer::doPFTaus(edm::Event& iEvent,const edm::EventSetup& iES)
{
      auto_ptr<LorentzVectorCollection> product_PFTaus(new LorentzVectorCollection);
      //Retrieve the collection
      edm::Handle<PFTauCollection> pftaus;
      if(iEvent.getByLabel(PFTaus_,pftaus))
    {
      for(unsigned int i=0;i<pftaus->size();++i)
        if((*pftaus)[i].pt()>ptMinPFTau_&&fabs((*pftaus)[i].eta())<etaMax)
          {
        reco::PFTauRef thePFTau(pftaus,i);
        for(unsigned int j=0;j<PFTauDis_.size();++j)
        {
          edm::Handle<PFTauDiscriminator> pftaudis;
          if(iEvent.getByLabel(PFTauDis_[j],pftaudis))
            {
              if((*pftaudis)[thePFTau]>0.5)
            {
              LorentzVector vec((*pftaus)[i].px(),(*pftaus)[i].py(),(*pftaus)[i].pz(),(*pftaus)[i].energy());
              product_PFTaus->push_back(vec);
            }
            }
        }
          }
    }

      iEvent.put(product_PFTaus,"PFTaus");
      
}


void 
HLTTauRefProducer::doElectrons(edm::Event& iEvent,const edm::EventSetup& iES)
{
  auto_ptr<LorentzVectorCollection> product_Electrons(new LorentzVectorCollection);
  //Retrieve the collections
  
  edm::Handle<reco::ElectronIDAssociationCollection> pEleID;
  if(e_doID_){//UGLY HACK UNTIL GET ELETRON ID WORKING IN 210
   
      iEvent.getByLabel(e_idAssocProd_,pEleID);
    
    if (!pEleID.isValid()){
      edm::LogInfo("")<< "Error! Can't get electronIDAssocProducer by label. ";
      e_doID_ = false;
    }
  }
  edm::Handle<reco::TrackCollection> pCtfTracks;
  iEvent.getByLabel(e_ctfTrackCollection_, pCtfTracks);
  if (!pCtfTracks.isValid()) {
  edm::LogInfo("")<< "Error! Can't get " << e_ctfTrackCollection_.label() << " by label. ";
  iEvent.put(product_Electrons,"Electrons"); 
  return;
  }
  const reco::TrackCollection * ctfTracks = pCtfTracks.product();
  edm::Handle<GsfElectronCollection> electrons;
  if(iEvent.getByLabel(Electrons_,electrons))
    for(size_t i=0;i<electrons->size();++i)
      {
    edm::Ref<reco::GsfElectronCollection> electronRef(electrons,i);
    bool idDec=false;
    if(e_doID_){
      reco::ElectronIDAssociationCollection::const_iterator tagIDAssocItr;
      tagIDAssocItr = pEleID->find(electronRef);
      const reco::ElectronIDRef& id_tag = tagIDAssocItr->val;
      idDec=id_tag->cutBasedDecision();
    }else idDec=true;
    if((*electrons)[i].pt()>ptMinElectron_&&fabs((*electrons)[i].eta())<etaMax&&idDec)
      {
        if(e_doTrackIso_){
          reco::TrackCollection::const_iterator tr = ctfTracks->begin();
          double sum_of_pt_ele=0;
          for(;tr != ctfTracks->end();++tr)
        {
          double lip = (*electrons)[i].gsfTrack()->dz() - tr->dz();
          if(tr->pt() > e_trackMinPt_ && fabs(lip) < e_lipCut_){
            double dphi=fabs(tr->phi()-(*electrons)[i].trackMomentumAtVtx().phi());
            if(dphi>acos(-1.))dphi=2*acos(-1.)-dphi;
            double deta=fabs(tr->eta()-(*electrons)[i].trackMomentumAtVtx().eta());
            double dr_ctf_ele = sqrt(deta*deta+dphi*dphi);
            if((dr_ctf_ele>e_minIsoDR_) && (dr_ctf_ele<e_maxIsoDR_)){
              double cft_pt_2 = (tr->pt())*(tr->pt());
              sum_of_pt_ele += cft_pt_2;
            }
          }
        }
          double isolation_value_ele = sum_of_pt_ele/((*electrons)[i].trackMomentumAtVtx().Rho()*(*electrons)[i].trackMomentumAtVtx().Rho());
          if(isolation_value_ele<e_isoMaxSumPt_){
        LorentzVector vec((*electrons)[i].px(),(*electrons)[i].py(),(*electrons)[i].pz(),(*electrons)[i].energy());
        product_Electrons->push_back(vec);
          } 
          
        }
        else{ 
          LorentzVector vec((*electrons)[i].px(),(*electrons)[i].py(),(*electrons)[i].pz(),(*electrons)[i].energy());
          product_Electrons->push_back(vec);
        }
      }
      }
  
  iEvent.put(product_Electrons,"Electrons"); 
}

void 
HLTTauRefProducer::doMuons(edm::Event& iEvent,const edm::EventSetup& iES)
{
  auto_ptr<LorentzVectorCollection> product_Muons(new LorentzVectorCollection);
  //Retrieve the collection
  edm::Handle<MuonCollection> muons;
      if(iEvent.getByLabel(Muons_,muons))
   
      for(size_t i = 0 ;i<muons->size();++i)
    {
     
        if((*muons)[i].pt()>ptMinMuon_&&fabs((*muons)[i].eta())<etaMax)
          {
        LorentzVector vec((*muons)[i].px(),(*muons)[i].py(),(*muons)[i].pz(),(*muons)[i].energy());
        product_Muons->push_back(vec);
          }
    }


      iEvent.put(product_Muons,"Muons");
 
}


void 
HLTTauRefProducer::doJets(edm::Event& iEvent,const edm::EventSetup& iES)
{
      auto_ptr<LorentzVectorCollection> product_Jets(new LorentzVectorCollection);
      //Retrieve the collection
      edm::Handle<CaloJetCollection> jets;
      if(iEvent.getByLabel(Jets_,jets))
      for(size_t i = 0 ;i<jets->size();++i)
    {
         if((*jets)[i].et()>ptMinJet_&&fabs((*jets)[i].eta())<etaMax)
          {
        LorentzVector vec((*jets)[i].px(),(*jets)[i].py(),(*jets)[i].pz(),(*jets)[i].energy());
        product_Jets->push_back(vec);
          }
    }
      iEvent.put(product_Jets,"Jets");
}

void 
HLTTauRefProducer::doTowers(edm::Event& iEvent,const edm::EventSetup& iES)
{
      auto_ptr<LorentzVectorCollection> product_Towers(new LorentzVectorCollection);
      //Retrieve the collection
      edm::Handle<CaloTowerCollection> towers;
      if(iEvent.getByLabel(Towers_,towers))
      for(size_t i = 0 ;i<towers->size();++i)
    {
         if((*towers)[i].pt()>ptMinTower_&&fabs((*towers)[i].eta())<etaMax)
          {
        //calculate isolation
        double isolET=0;
        for(unsigned int j=0;j<towers->size();++j)
          {
            if(ROOT::Math::VectorUtil::DeltaR((*towers)[i].p4(),(*towers)[j].p4())<0.5)
              isolET+=(*towers)[j].pt();
          }
        isolET-=(*towers)[i].pt();
        if(isolET<towerIsol_)
          {
            LorentzVector vec((*towers)[i].px(),(*towers)[i].py(),(*towers)[i].pz(),(*towers)[i].energy());
            product_Towers->push_back(vec);
          }
          }
    }
      iEvent.put(product_Towers,"Towers");
}


void 
HLTTauRefProducer::doPhotons(edm::Event& iEvent,const edm::EventSetup& iES)
{
      auto_ptr<LorentzVectorCollection> product_Gammas(new LorentzVectorCollection);
      //Retrieve the collection
      edm::Handle<PhotonCollection> photons;
      if(iEvent.getByLabel(Photons_,photons))
      for(size_t i = 0 ;i<photons->size();++i)
    if((*photons)[i].ecalRecHitSumEtConeDR04()<photonEcalIso_)
    {
         if((*photons)[i].et()>ptMinPhoton_&&fabs((*photons)[i].eta())<etaMax)
          {
        LorentzVector vec((*photons)[i].px(),(*photons)[i].py(),(*photons)[i].pz(),(*photons)[i].energy());
        product_Gammas->push_back(vec);
          }
    }
      iEvent.put(product_Gammas,"Photons");
}