HigPhotonJetHLTOfflineSource.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     HigPhotonJetHLTOfflineSource
// Class:       HigPhotonJetHLTOfflineSource
// 

//
// Author: Xin Shi <Xin.Shi@cern.ch> 
// Created: 2014.07.22 
//

// system include files
#include <memory>
#include <iostream>

// user include files
#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/MonitorElement.h"
#include "DQMServices/Core/interface/DQMEDAnalyzer.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Common/interface/TriggerNames.h"

#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"

#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"
#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
#include "DataFormats/METReco/interface/PFMETCollection.h"
#include "DataFormats/METReco/interface/PFMET.h"
#include "DataFormats/JetReco/interface/PFJetCollection.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include "DataFormats/Math/interface/deltaPhi.h"

#include <TLorentzVector.h>
#include <TH2F.h>

//  Define the interface
class HigPhotonJetHLTOfflineSource : public DQMEDAnalyzer {

public:

  explicit HigPhotonJetHLTOfflineSource(const edm::ParameterSet&);

private:

  // Analyzer Methods
  void dqmBeginRun(const edm::Run &,
               const edm::EventSetup &) override;
  void bookHistograms(DQMStore::IBooker &,
                  edm::Run const &,
                  edm::EventSetup const &) override;  
  void analyze(const edm::Event &,
               const edm::EventSetup &) override;
  void endRun(const edm::Run &,
              const edm::EventSetup &) override;
  bool isMonitoredTriggerAccepted(const edm::TriggerNames&,
                  const edm::Handle<edm::TriggerResults>&); 

  // Input from Configuration File
  edm::ParameterSet pset_;
  std::string hltProcessName_;
  std::vector<std::string> hltPathsToCheck_;
  std::string dirname_; 
  bool verbose_;
  bool triggerAccept_; 
    
  edm::EDGetTokenT <edm::TriggerResults> triggerResultsToken_;
  edm::EDGetTokenT<reco::VertexCollection> pvToken_;
  edm::EDGetTokenT<reco::PhotonCollection> photonsToken_;
  edm::EDGetTokenT<reco::PFMETCollection> pfMetToken_;
  edm::EDGetTokenT<reco::PFJetCollection> pfJetsToken_;

  double pfjetMinPt_;  
  double photonMinPt_;  

  // Member Variables

  MonitorElement*  nvertices_reco_;
  MonitorElement*  nvertices_;
  MonitorElement*  nphotons_reco_;
  MonitorElement*  nphotons_;
  MonitorElement*  photonpt_reco_;
  MonitorElement*  photonpt_;
  MonitorElement*  photonrapidity_reco_;
  MonitorElement*  photonrapidity_;
  MonitorElement*  pfmet_reco_;
  MonitorElement*  pfmet_;
  MonitorElement*  pfmetphi_reco_;
  MonitorElement*  pfmetphi_;
  MonitorElement*  npfjets_reco_;
  MonitorElement*  npfjets_;
  MonitorElement*  delphiphomet_reco_;
  MonitorElement*  delphiphomet_;
  MonitorElement*  delphijetmet_reco_;
  MonitorElement*  delphijetmet_;
  MonitorElement*  invmassjj_reco_;
  MonitorElement*  invmassjj_;
  MonitorElement*  deletajj_reco_;
  MonitorElement*  deletajj_;
  MonitorElement*  triggers_reco_;
  MonitorElement*  triggers_;
  MonitorElement*  trigvsnvtx_reco_;
  MonitorElement*  trigvsnvtx_;
  
  double evtsrun_; 
  
};


// Class Methods 

HigPhotonJetHLTOfflineSource::HigPhotonJetHLTOfflineSource(const edm::ParameterSet& pset) :
  pset_(pset)
{
  hltProcessName_ = pset.getParameter<std::string>("hltProcessName"); 
  hltPathsToCheck_ = pset.getParameter<std::vector<std::string>>("hltPathsToCheck"); 
  verbose_ = pset.getUntrackedParameter<bool>("verbose", false);
  triggerAccept_ = pset.getUntrackedParameter<bool>("triggerAccept", true);
  triggerResultsToken_ = consumes <edm::TriggerResults> (pset.getParameter<edm::InputTag>("triggerResultsToken"));
  dirname_ = pset.getUntrackedParameter<std::string>("dirname", std::string("HLT/Higgs/PhotonJet/"));
  pvToken_ = consumes<reco::VertexCollection> (pset.getParameter<edm::InputTag>("pvToken"));
  photonsToken_ = consumes<reco::PhotonCollection> (pset.getParameter<edm::InputTag>("photonsToken"));
  pfMetToken_ = consumes<reco::PFMETCollection> (pset.getParameter<edm::InputTag>("pfMetToken"));
  pfJetsToken_ = consumes<reco::PFJetCollection> (pset.getParameter<edm::InputTag>("pfJetsToken"));
  pfjetMinPt_ = pset.getUntrackedParameter<double>("pfjetMinPt", 0.0);
  photonMinPt_ = pset.getUntrackedParameter<double>("photonMinPt", 0.0);
}

void 
HigPhotonJetHLTOfflineSource::dqmBeginRun(const edm::Run & iRun, 
                      const edm::EventSetup & iSetup) 
{ // Initialize hltConfig
  HLTConfigProvider hltConfig;
  bool changedConfig;
  if (!hltConfig.init(iRun, iSetup, hltProcessName_, changedConfig)) {
    edm::LogError("HLTPhotonJetVal") << "Initialization of HLTConfigProvider failed!!"; 
    return;
  }
 
  evtsrun_ = 0; 
}


void 
HigPhotonJetHLTOfflineSource::bookHistograms(DQMStore::IBooker & iBooker, 
                         edm::Run const & iRun,
                         edm::EventSetup const & iSetup)
{
  iBooker.setCurrentFolder(dirname_);
  nvertices_reco_ = iBooker.book1D("nvertices_reco", "Reco: Number of vertices", 100, 0, 100); 
  nvertices_ = iBooker.book1D("nvertices", "Number of vertices", 100, 0, 100); 
  nphotons_reco_ = iBooker.book1D("nphotons_reco", "Reco: Number of photons", 100, 0, 10); 
  nphotons_ = iBooker.book1D("nphotons", "Number of photons", 100, 0, 10); 
  photonpt_reco_ = iBooker.book1D("photonpt_reco", "Reco: Photons pT", 100, 0, 500); 
  photonpt_ = iBooker.book1D("photonpt", "Photons pT", 100, 0, 500); 
  photonrapidity_reco_ = iBooker.book1D("photonrapidity_reco", "Reco: Photons rapidity;y_{#gamma}", 100, -2.5, 2.5); 
  photonrapidity_ = iBooker.book1D("photonrapidity", "Photons rapidity;y_{#gamma}", 100, -2.5, 2.5); 
  pfmet_reco_ = iBooker.book1D("pfmet_reco", "Reco: PF MET", 100, 0, 250); 
  pfmet_ = iBooker.book1D("pfmet", "PF MET", 100, 0, 250); 
  pfmetphi_reco_ = iBooker.book1D("pfmetphi_reco", "Reco: PF MET phi;#phi_{PFMET}", 100, -4, 4); 
  pfmetphi_ = iBooker.book1D("pfmetphi", "PF MET phi;#phi_{PFMET}", 100, -4, 4); 
  delphiphomet_reco_ = iBooker.book1D("delphiphomet_reco", "Reco: #Delta#phi(photon, MET);#Delta#phi(#gamma,MET)", 100, 0, 4); 
  delphiphomet_ = iBooker.book1D("delphiphomet", "#Delta#phi(photon, MET);#Delta#phi(#gamma,MET)", 100, 0, 4); 
  npfjets_reco_ = iBooker.book1D("npfjets_reco", "Reco: Number of PF Jets", 100, 0, 20); 
  npfjets_ = iBooker.book1D("npfjets", "Number of PF Jets", 100, 0, 20); 
  delphijetmet_reco_ = iBooker.book1D("delphijetmet_reco", "Reco: #Delta#phi(PFJet, MET);#Delta#phi(Jet,MET)", 100, 0, 4); 
  delphijetmet_ = iBooker.book1D("delphijetmet", "#Delta#phi(PFJet, MET);#Delta#phi(Jet,MET)", 100, 0, 4); 
  invmassjj_reco_ = iBooker.book1D("invmassjj_reco", "Reco: Inv mass two leading jets;M_{jj}[GeV]", 100, 0, 2000); 
  invmassjj_ = iBooker.book1D("invmassjj", "Inv mass two leading jets;M_{jj}[GeV]", 100, 0, 2000); 
  deletajj_reco_ = iBooker.book1D("deletajj_reco", "Reco: #Delta#eta(jj);|#Delta#eta_{jj}|", 100, 0, 6); 
  deletajj_ = iBooker.book1D("deletajj", "#Delta#eta(jj);|#Delta#eta_{jj}|", 100, 0, 6); 
  triggers_reco_ = iBooker.book1D("triggers_reco", "Reco: Triggers", hltPathsToCheck_.size(), 0, hltPathsToCheck_.size());
  triggers_ = iBooker.book1D("triggers", "Triggers", hltPathsToCheck_.size(), 0, hltPathsToCheck_.size());
  trigvsnvtx_reco_ = iBooker.book2D("trigvsnvtx_reco", "Reco: Trigger vs. # vertices;N_{vertices};Trigger", 100, 0, 100, hltPathsToCheck_.size(), 0, hltPathsToCheck_.size()); 
  trigvsnvtx_ = iBooker.book2D("trigvsnvtx", "Trigger vs. # vertices;N_{vertices};Trigger", 100, 0, 100, hltPathsToCheck_.size(), 0, hltPathsToCheck_.size()); 
}


void
HigPhotonJetHLTOfflineSource::analyze(const edm::Event& iEvent, 
                      const edm::EventSetup& iSetup)
{
  // Count total number of events in one run
  evtsrun_++; 

  edm::Handle<edm::TriggerResults> triggerResults;
  iEvent.getByToken(triggerResultsToken_, triggerResults);
  if(!triggerResults.isValid()) {
      edm::LogError("HigPhotonJetHLT")<<"Missing triggerResults collection" << std::endl;
      return;
  }

  // Check whether contains monitored trigger and accepted
  const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults); 
  bool triggered = isMonitoredTriggerAccepted(triggerNames, triggerResults); 

  // if (!triggered) return; 

  // Test scale 
  // if (evtsrun_ > 10) return;

  // N Vertices 
  edm::Handle<reco::VertexCollection> vertices;
  iEvent.getByToken(pvToken_, vertices);
  if(!vertices.isValid()) return;  
  if (verbose_)
    std::cout << "xshi:: N vertices : " << vertices->size() << std::endl;
  
  // Set trigger name labels
  for (size_t i = 0; i < hltPathsToCheck_.size(); i++) {
    triggers_->setBinLabel(i+1, hltPathsToCheck_[i]); 
  }

  // Fill trigger info
  for (unsigned int itrig = 0; itrig < triggerResults->size(); itrig++){
    const std::string& triggername = triggerNames.triggerName(itrig);
    for (size_t i = 0; i < hltPathsToCheck_.size(); i++) {
      if ( triggername.find(hltPathsToCheck_[i]) != std::string::npos) {
     triggers_reco_->Fill(i);
     trigvsnvtx_reco_->Fill(vertices->size(), i); 
     if (triggered) triggers_->Fill(i);
     if (triggered) trigvsnvtx_->Fill(vertices->size(), i);
      }
    }
  }

  nvertices_reco_->Fill(vertices->size());
  if (triggered) nvertices_->Fill(vertices->size());

  // PF MET
  edm::Handle<reco::PFMETCollection> pfmets;
  iEvent.getByToken(pfMetToken_, pfmets);
  if (!pfmets.isValid()) return;
  const reco::PFMET pfmet = pfmets->front();
  pfmet_reco_->Fill(pfmet.et()); 
  if (triggered) pfmet_->Fill(pfmet.et()); 
  if (verbose_)
    std::cout << "xshi:: number of pfmets: " << pfmets->size() << std::endl;

  pfmetphi_reco_->Fill(pfmet.phi()); 
  if (triggered) pfmetphi_->Fill(pfmet.phi()); 
  
  // Photons
  edm::Handle<reco::PhotonCollection> photons;
  iEvent.getByToken(photonsToken_, photons);
  if(!photons.isValid()) return;
  int nphotons = 0; 
  for(auto const & phoIter : *photons){
    if (phoIter.pt() < photonMinPt_ )  continue;
    nphotons++;
    photonpt_reco_->Fill(phoIter.pt());
    photonrapidity_reco_->Fill(phoIter.rapidity()); 
    if (triggered) photonpt_->Fill(phoIter.pt()); 
    if (triggered) photonrapidity_->Fill(phoIter.rapidity()); 
    double tmp_delphiphomet = fabs(deltaPhi(phoIter.phi(), pfmet.phi())); 
    delphiphomet_reco_->Fill(tmp_delphiphomet); 
    if (triggered) delphiphomet_->Fill(tmp_delphiphomet); 
  }
  nphotons_reco_->Fill(nphotons);
  if (triggered)  nphotons_->Fill(nphotons);
  
  // PF Jet
  edm::Handle<reco::PFJetCollection> pfjets;
  iEvent.getByToken(pfJetsToken_, pfjets);
  if(!pfjets.isValid()) return;
  if (verbose_)
    std::cout << "xshi:: N pfjets : " << pfjets->size() << std::endl;

  double min_delphijetmet = 6.0;
  TLorentzVector p4jet1, p4jet2, p4jj;
  // Two leading jets eta
  double etajet1(0), etajet2(0);
  int njet = 0;  
  for(auto const & jetIter : *pfjets){
    if (jetIter.pt() < pfjetMinPt_ ) continue; 
    njet++;

    double tmp_delphijetmet = fabs(deltaPhi(jetIter.phi(), pfmet.phi())); 
    if (tmp_delphijetmet < min_delphijetmet)
      min_delphijetmet = tmp_delphijetmet;

    if (njet == 1) {
      p4jet1.SetXYZM(jetIter.px(), jetIter.py(), jetIter.pz(), jetIter.mass()); 
      etajet1 = jetIter.eta(); 
    }
    if (njet == 2){
      p4jet2.SetXYZM(jetIter.px(), jetIter.py(), jetIter.pz(), jetIter.mass()); 
      etajet2 = jetIter.eta(); 
    }
  }
  npfjets_reco_->Fill(njet);   
  if (triggered) npfjets_->Fill(njet);   
  
  delphijetmet_reco_->Fill(min_delphijetmet); 
  if (triggered) delphijetmet_->Fill(min_delphijetmet); 
  p4jj = p4jet1 + p4jet2; 
  double deletajj = etajet1 - etajet2 ; 
  if (verbose_) 
    std::cout << "xshi:: invmass jj " << p4jj.M() << std::endl;
  
  invmassjj_reco_->Fill(p4jj.M());
  deletajj_reco_->Fill(deletajj); 
  if (triggered) invmassjj_->Fill(p4jj.M());
  if (triggered) deletajj_->Fill(deletajj); 
}


void 
HigPhotonJetHLTOfflineSource::endRun(const edm::Run & iRun, 
                     const edm::EventSetup& iSetup)
{
  // Normalize to the total number of events in the run
  TH2F* h = trigvsnvtx_->getTH2F();
  double integral = h->Integral();
  double norm = (integral > 0.) ? evtsrun_*hltPathsToCheck_.size()/integral : 1.;
  h->Scale(norm);
  if (verbose_) {
    std::cout << "xshi:: endRun total number of events: " << evtsrun_
          << ", integral = " << h->Integral()
          << ", norm = " << norm << std::endl;
  }
}

bool
HigPhotonJetHLTOfflineSource::isMonitoredTriggerAccepted(const edm::TriggerNames& triggerNames,
                             const edm::Handle<edm::TriggerResults>& triggerResults )
{
  for (unsigned int itrig = 0; itrig < triggerResults->size(); itrig++){
    // Only consider the triggered case.
    if ( triggerAccept_ && ( (*triggerResults)[itrig].accept() != 1) ) continue; 
    const std::string& triggername = triggerNames.triggerName(itrig);
    for (auto const & i : hltPathsToCheck_) {
      if ( triggername.find(i) != std::string::npos) {
    return true;
      }
    }
  }

  return false; 
}

//define this as a plug-in
DEFINE_FWK_MODULE(HigPhotonJetHLTOfflineSource);