/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/DQM/HLTEvF/src/HLTTauDQMCaloPlotter.cc

Go to the documentation of this file.

#include "DQM/HLTEvF/interface/HLTTauDQMCaloPlotter.h"
#include "Math/GenVector/VectorUtil.h"
#include <iostream>
#include <iomanip>
#include <fstream>

HLTTauDQMCaloPlotter::HLTTauDQMCaloPlotter(const edm::ParameterSet& iConfig,int etbins,int etabins,int phibins,double maxpt,bool ref,double dr):
 l2preJets_(iConfig.getParameter<std::vector<edm::InputTag> >("L2RegionalJets")),
 l2TauInfoAssoc_(iConfig.getParameter<edm::InputTag>("L2InfoAssociationInput")),
 doRef_(ref),
 matchDeltaRMC_(dr),
 triggerTag_(iConfig.getParameter<std::string>("DQMFolder")),
 l2Isolated_(iConfig.getParameter<edm::InputTag>("L2IsolatedJets")),
 EtMax_(maxpt),
 NPtBins_(etbins),
 NEtaBins_(etabins),
 NPhiBins_(phibins)
{
   store = &*edm::Service<DQMStore>();
 
  if(store)
    {

      //Create the histograms
      store->setCurrentFolder(triggerTag_);

      preJetEt= store->book1D("L2PreTauEt","L2 regional #tau E_{t}",NPtBins_,0,EtMax_);
       preJetEt->getTH1F()->GetXaxis()->SetTitle("L2 regional Jet E_{T}");
       preJetEt->getTH1F()->GetYaxis()->SetTitle("entries");

    preJetEta= store->book1D("L2PreTauEta","L2 regional #tau #eta",NEtaBins_,-2.5,2.5);
       preJetEta->getTH1F()->GetXaxis()->SetTitle("L2 regional Jet #eta");
       preJetEta->getTH1F()->GetYaxis()->SetTitle("entries");

    preJetPhi= store->book1D("L2PreTauPhi","L2 regional #tau #phi",NPhiBins_,-3.2,3.2);
       preJetPhi->getTH1F()->GetXaxis()->SetTitle("L2 regional Jet #phi");
       preJetPhi->getTH1F()->GetYaxis()->SetTitle("entries");

    jetEt= store->book1D("L2TauEt","L2 #tau E_{t}",NPtBins_,0,EtMax_);
       jetEt->getTH1F()->GetXaxis()->SetTitle("L2 selected Jet E_{T}");
       jetEt->getTH1F()->GetYaxis()->SetTitle("entries");

    jetEta= store->book1D("L2TauEta","L2 #tau #eta",NEtaBins_,-2.5,2.5);
       jetEta->getTH1F()->GetXaxis()->SetTitle("L2 selected Jet #eta");
       jetEta->getTH1F()->GetYaxis()->SetTitle("entries");

    jetPhi= store->book1D("L2TauPhi","L2 #tau #phi",NPhiBins_,-3.2,3.2);
       jetPhi->getTH1F()->GetXaxis()->SetTitle("L2 selected Jet #phi");
       jetPhi->getTH1F()->GetYaxis()->SetTitle("entries");

    jetEtRes= store->book1D("L2TauEtResol","L2 #tau E_{t} resolution",40,-2,2);
       jetEtRes->getTH1F()->GetXaxis()->SetTitle("L2 selected Jet #phi");
       jetEtRes->getTH1F()->GetYaxis()->SetTitle("entries");

    isoJetEt= store->book1D("L2IsoTauEt","L2 isolated #tau E_{t}",NPtBins_,0,EtMax_);
       isoJetEt->getTH1F()->GetXaxis()->SetTitle("L2 isolated Jet E_{T}");
       isoJetEt->getTH1F()->GetYaxis()->SetTitle("entries");

    isoJetEta= store->book1D("L2IsoTauEta","L2 isolated #tau #eta",NEtaBins_,-2.5,2.5);
       isoJetEta->getTH1F()->GetXaxis()->SetTitle("L2 isolated Jet #eta");
       isoJetEta->getTH1F()->GetYaxis()->SetTitle("entries");

    isoJetPhi= store->book1D("L2IsoTauPhi","L2 isolated #tau #phi",NPhiBins_,-3.2,3.2);
       isoJetPhi->getTH1F()->GetXaxis()->SetTitle("L2 isolated Jet #phi");
       isoJetPhi->getTH1F()->GetYaxis()->SetTitle("entries");

    ecalIsolEt=store->book1D("L2EcalIsolation","ECAL Isolation",40,0,20);
       ecalIsolEt->getTH1F()->GetXaxis()->SetTitle("L2 ECAL isolation E_{T}");
       ecalIsolEt->getTH1F()->GetYaxis()->SetTitle("entries");

    hcalIsolEt=store->book1D("L2HcalIsolation","HCAL Isolation",40,0,20);
       hcalIsolEt->getTH1F()->GetXaxis()->SetTitle("L2 HCAL isolation E_{T}");
       hcalIsolEt->getTH1F()->GetYaxis()->SetTitle("entries");

    seedHcalEt=store->book1D("L2HighestHCALCluster","Highest HCAL Cluster",40,0,80);
       seedHcalEt->getTH1F()->GetXaxis()->SetTitle("HCAL seed  E_{T}");
       seedHcalEt->getTH1F()->GetYaxis()->SetTitle("entries");

    seedEcalEt=store->book1D("L2HighestECALCluster","Highest ECAL Cluster",25,0,50);
       seedEcalEt->getTH1F()->GetXaxis()->SetTitle("ECAL seed  E_{T}");
       seedEcalEt->getTH1F()->GetYaxis()->SetTitle("entries");

    nEcalClusters=store->book1D("L2NEcalClusters","Nucmber of ECAL Clusters",20,0,20);
       nEcalClusters->getTH1F()->GetXaxis()->SetTitle("n. of ECAL Clusters");
       nEcalClusters->getTH1F()->GetYaxis()->SetTitle("entries");

    ecalClusterEtaRMS=store->book1D("L2EcalEtaRMS","ECAL Cluster #eta RMS",15,0,0.05);
       ecalClusterEtaRMS->getTH1F()->GetXaxis()->SetTitle("ECAL cluster #eta RMS");
       ecalClusterEtaRMS->getTH1F()->GetYaxis()->SetTitle("entries");

    ecalClusterPhiRMS=store->book1D("L2EcalPhiRMS","ECAL Cluster #phi RMS",30,0,0.1);
       ecalClusterPhiRMS->getTH1F()->GetXaxis()->SetTitle("ECAL cluster #phi RMS");
       ecalClusterPhiRMS->getTH1F()->GetYaxis()->SetTitle("entries");


    ecalClusterDeltaRRMS=store->book1D("L2EcalDeltaRRMS","ECAL Cluster #DeltaR RMS",30,0,0.1);
       ecalClusterDeltaRRMS->getTH1F()->GetXaxis()->SetTitle("ECAL cluster #DeltaR RMS");
       ecalClusterDeltaRRMS->getTH1F()->GetYaxis()->SetTitle("entries");

    nHcalClusters=store->book1D("L2NHcalClusters","Nucmber of HCAL Clusters",20,0,20);
       nHcalClusters->getTH1F()->GetXaxis()->SetTitle("n. of ECAL Clusters");
       nHcalClusters->getTH1F()->GetYaxis()->SetTitle("entries");

      hcalClusterEtaRMS=store->book1D("L2HcalEtaRMS","HCAL Cluster #eta RMS",15,0,0.05);
       hcalClusterEtaRMS->getTH1F()->GetXaxis()->SetTitle("HCAL cluster #eta RMS");
       hcalClusterEtaRMS->getTH1F()->GetYaxis()->SetTitle("entries");

      hcalClusterPhiRMS=store->book1D("L2HcalPhiRMS","HCAL Cluster #phi RMS",30,0,0.1);
       hcalClusterPhiRMS->getTH1F()->GetXaxis()->SetTitle("HCAL cluster #phi RMS");
       hcalClusterPhiRMS->getTH1F()->GetYaxis()->SetTitle("entries");

      hcalClusterDeltaRRMS=store->book1D("L2HcalDeltaRRMS","HCAL Cluster #DeltaR RMS",30,0,0.1);
       hcalClusterDeltaRRMS->getTH1F()->GetXaxis()->SetTitle("HCAL cluster #DeltaR RMS");
       hcalClusterDeltaRRMS->getTH1F()->GetYaxis()->SetTitle("entries");


      store->setCurrentFolder(triggerTag_+"/EfficiencyHelpers");

      recoEtEffNum=store->book1D("L2RecoTauEtEffNum","Efficiency vs E_{t}(Numerator)",NPtBins_,0,EtMax_);
      recoEtEffNum->getTH1F()->Sumw2();

      recoEtEffDenom=store->book1D("L2RecoTauEtEffDenom","Efficiency vs E_{t}(Denominator)",NPtBins_,0,EtMax_);
      recoEtEffDenom->getTH1F()->Sumw2();

      recoEtaEffNum=store->book1D("L2RecoTauEtaEffNum","Efficiency vs #eta (Numerator)",NEtaBins_,-2.5,2.5);
      recoEtaEffNum->getTH1F()->Sumw2();

      recoEtaEffDenom=store->book1D("L2RecoTauEtaEffDenom","Efficiency vs #eta(Denominator)",NEtaBins_,-2.5,2.5);
      recoEtaEffDenom->getTH1F()->Sumw2();

      recoPhiEffNum=store->book1D("L2RecoTauPhiEffNum","Efficiency vs #phi (Numerator)",NPhiBins_,-3.2,3.2);
      recoPhiEffNum->getTH1F()->Sumw2();

      recoPhiEffDenom=store->book1D("L2RecoTauPhiEffDenom","Efficiency vs #phi(Denominator)",NPhiBins_,-3.2,3.2);
      recoPhiEffDenom->getTH1F()->Sumw2();

      isoEtEffNum=store->book1D("L2IsoTauEtEffNum","Efficiency vs E_{t}(Numerator)",NPtBins_,0,EtMax_);
      isoEtEffNum->getTH1F()->Sumw2();

      isoEtEffDenom=store->book1D("L2IsoTauEtEffDenom","Efficiency vs E_{t}(Denominator)",NPtBins_,0,EtMax_);
      isoEtEffDenom->getTH1F()->Sumw2();

      isoEtaEffNum=store->book1D("L2IsoTauEtaEffNum","Efficiency vs #eta (Numerator)",NEtaBins_,-2.5,2.5);
      isoEtaEffNum->getTH1F()->Sumw2();

      isoEtaEffDenom=store->book1D("L2IsoTauEtaEffDenom","Efficiency vs #eta(Denominator)",NEtaBins_,-2.5,2.5);
      isoEtaEffDenom->getTH1F()->Sumw2();

      isoPhiEffNum=store->book1D("L2IsoTauPhiEffNum","Efficiency vs #phi (Numerator)",NPhiBins_,-3.2,3.2);
      isoPhiEffNum->getTH1F()->Sumw2();

      isoPhiEffDenom=store->book1D("L2IsoTauPhiEffDenom","Efficiency vs #phi(Denominator)",NPhiBins_,-3.2,3.2);
      isoPhiEffDenom->getTH1F()->Sumw2();
    }

}

HLTTauDQMCaloPlotter::~HLTTauDQMCaloPlotter()
{
}

void
HLTTauDQMCaloPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup,const LVColl& McInfo)
{
   edm::Handle<reco::L2TauInfoAssociation> l2TauInfoAssoc;
   edm::Handle<reco::CaloJetCollection> l2Isolated;
   edm::Handle<reco::CaloJetCollection> l2Regional;
   reco::CaloJetCollection l2RegionalJets;


   
   //Merge the L2 Regional Collections!
   reco::CaloJetCollection l2MergedJets;

   for(unsigned int j=0;j<l2preJets_.size();++j) {

     bool  gotPreJets= iEvent.getByLabel(l2preJets_[j],l2Regional) &&l2Regional.isValid();

     if(gotPreJets)
       if((!l2Regional.failedToGet())) {
         for(unsigned int i=0;i<l2Regional->size();++i) 
           l2MergedJets.push_back(l2Regional->at(i));

       }
   }

   //Sort
   SorterByPt sorter;
   std::sort(l2MergedJets.begin(),l2MergedJets.end(),sorter);

   //Remove Collinear Jets
   reco::CaloJetCollection l2CleanJets;
   while(l2MergedJets.size()>0) {
     l2CleanJets.push_back(l2MergedJets.at(0));
     reco::CaloJetCollection tmp;
     for(unsigned int i=1 ;i<l2MergedJets.size();++i) {
       double DR = ROOT::Math::VectorUtil::DeltaR(l2MergedJets.at(0).p4(),l2MergedJets.at(i).p4());
       if(DR>0.1) 
         tmp.push_back(l2MergedJets.at(i));
     }

     l2MergedJets.swap(tmp);
     tmp.clear();
   }


   //Now fill the regional jet plots by ref if you do ref to 
   //avoid double counting!

   if(doRef_)     {
       for(unsigned int i=0;i<McInfo.size();++i) {
         std::pair<bool,reco::CaloJet> m = inverseMatch(McInfo.at(i),l2CleanJets);
         if(m.first) {
           preJetEt->Fill(m.second.pt());
           preJetEta->Fill(m.second.eta());
           preJetPhi->Fill(m.second.phi());
           recoEtEffDenom->Fill(McInfo.at(i).pt());
           recoEtaEffDenom->Fill(McInfo.at(i).eta());
           recoPhiEffDenom->Fill(McInfo.at(i).phi());
           l2RegionalJets.push_back(m.second);
         }
       }
   }
   else {
     for(unsigned int i=0;i<l2CleanJets.size();++i) {
       reco::CaloJet jet = l2CleanJets.at(i);
           preJetEt->Fill(jet.pt());
           preJetEta->Fill(jet.eta());
           preJetPhi->Fill(jet.phi());
           recoEtEffDenom->Fill(jet.pt());
           recoEtaEffDenom->Fill(jet.eta());
           recoPhiEffDenom->Fill(jet.phi());
           l2RegionalJets.push_back(jet);
     }
   }
      




   bool gotL2=iEvent.getByLabel(l2TauInfoAssoc_,l2TauInfoAssoc) && l2TauInfoAssoc.isValid();



     if(gotL2)     {
       //If the Collection exists do work
       if(l2TauInfoAssoc->size()>0)
         for(reco::L2TauInfoAssociation::const_iterator p = l2TauInfoAssoc->begin();p!=l2TauInfoAssoc->end();++p)
           {
             //Retrieve The L2TauIsolationInfo Class from the AssociationMap
             const reco::L2TauIsolationInfo l2info = p->val;
       
             //Retrieve the Jet From the AssociationMap
             const reco::Jet& jet =*(p->key);
             
             std::pair<bool,LV> m =match(jet,McInfo); 
             
             if((doRef_&&m.first)||(!doRef_))
                 {
                   ecalIsolEt->Fill(l2info.ecalIsolEt());
                   hcalIsolEt->Fill(l2info.hcalIsolEt());
                   seedEcalEt->Fill(l2info.seedEcalHitEt());
                   seedHcalEt->Fill(l2info.seedHcalHitEt());

                   nEcalClusters->Fill(l2info.nEcalHits());
                   ecalClusterEtaRMS->Fill(l2info.ecalClusterShape()[0]);
                   ecalClusterPhiRMS->Fill(l2info.ecalClusterShape()[1]);
                   ecalClusterDeltaRRMS->Fill(l2info.ecalClusterShape()[2]);

                   nHcalClusters->Fill(l2info.nHcalHits());
                   hcalClusterEtaRMS->Fill(l2info.hcalClusterShape()[0]);
                   hcalClusterPhiRMS->Fill(l2info.hcalClusterShape()[1]);
                   hcalClusterDeltaRRMS->Fill(l2info.hcalClusterShape()[2]);

                   jetEt->Fill(jet.et());
                   jetEta->Fill(jet.eta());
                   jetPhi->Fill(jet.phi());

                   LV refLV;
                   if(doRef_)
                     refLV = m.second;
                   else
                     refLV = jet.p4();
                   
                   if(doRef_)
                     jetEtRes->Fill((jet.pt()-refLV.pt())/refLV.pt());

                   if(matchJet(jet,l2RegionalJets))
                     {
                       recoEtEffNum->Fill(refLV.pt());
                       recoEtaEffNum->Fill(refLV.eta());
                       recoPhiEffNum->Fill(refLV.phi());
                     }

                   isoEtEffDenom->Fill(refLV.pt());
                   isoEtaEffDenom->Fill(refLV.eta());
                   isoPhiEffDenom->Fill(refLV.phi());


                   bool  gotIsoL2=iEvent.getByLabel(l2Isolated_,l2Isolated) &&l2Isolated.isValid();


        

                   if(gotIsoL2)
                     if(l2Isolated.isValid());
                     {
                       if(matchJet(jet,*l2Isolated))
                         { 
                           isoJetEt->Fill(jet.et());
                           isoJetEta->Fill(jet.eta());
                           isoJetPhi->Fill(jet.phi());

                           isoEtEffNum->Fill(refLV.pt());
                           isoEtaEffNum->Fill(refLV.eta());
                           isoPhiEffNum->Fill(refLV.phi());
                         }
                     }
                 }
           } 
               
     }

}




std::pair<bool,LV> 
HLTTauDQMCaloPlotter::match(const reco::Jet& jet,const LVColl& McInfo)
{

  //Loop On the Collection and see if your tau jet is matched to one there
 //Also find the nearest Matched MC Particle to your Jet (to be complete)
 
 bool matched=false;
 LV mLV;

 if(McInfo.size()>0)
  for(LVColl::const_iterator it = McInfo.begin();it!=McInfo.end();++it)
   {
          double delta = ROOT::Math::VectorUtil::DeltaR(jet.p4(),*it);
          if(delta<matchDeltaRMC_)
            {
              matched=true;
              mLV=*it;
            }
   }


 std::pair<bool,LV> p = std::make_pair(matched,mLV);
 return p;
}


std::pair<bool,reco::CaloJet> 
HLTTauDQMCaloPlotter::inverseMatch(const LV& jet,const reco::CaloJetCollection& jets)
{

  //Loop On the Collection and see if your tau jet is matched to one there
  
  //MATCH THE neartes energy jet in the delta R we want
 
 bool matched=false;
 reco::CaloJet mjet;
 double distance=100000;
 if(jets.size()>0)
   for(reco::CaloJetCollection::const_iterator it = jets.begin();it!=jets.end();++it)
     if(ROOT::Math::VectorUtil::DeltaR(it->p4(),jet)<matchDeltaRMC_)
       {
         matched=true;
         double delta=fabs(jet.pt()-it->pt());
         if(delta<distance)
           {
             distance=delta;
             mjet = *it;
           }
       }

 std::pair<bool,reco::CaloJet> p = std::make_pair(matched,mjet);
 return p;
}


bool 
HLTTauDQMCaloPlotter::matchJet(const reco::Jet& jet,const reco::CaloJetCollection& McInfo)
{
  //Loop On the Collection and see if your tau jet is matched to one there
 //Also find the nearest Matched MC Particle to your Jet (to be complete)
 
 bool matched=false;

 if(McInfo.size()>0)
  for(reco::CaloJetCollection::const_iterator it = McInfo.begin();it!=McInfo.end();++it)
   {
     //           double delta = ROOT::Math::VectorUtil::DeltaR(jet.p4(),it->p4());
     //   if(delta<matchDeltaRMC_)
     if(jet.p4()==it->p4())
       {
         matched=true;
         
       }
   }
 return matched;
}