SUSY_HLT_VBF_Mu.cc

Go to the documentation of this file.

#include "DataFormats/HLTReco/interface/TriggerObject.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "FWCore/Common/interface/TriggerNames.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "HLTriggerOffline/SUSYBSM/interface/SUSY_HLT_VBF_Mu.h"
#include "TLorentzVector.h"
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <regex.h>
#include <string>
#include <sys/types.h>

SUSY_HLT_VBF_Mu::SUSY_HLT_VBF_Mu(const edm::ParameterSet &ps) {
  edm::LogInfo("SUSY_HLT_VBF_Mu") << "Constructor SUSY_HLT_VBF_Mu::SUSY_HLT_VBF_Mu " << std::endl;
  // Get parameters from configuration file
  theTrigSummary_ = consumes<trigger::TriggerEvent>(ps.getParameter<edm::InputTag>("trigSummary"));
  theMuonCollection_ = consumes<reco::MuonCollection>(ps.getParameter<edm::InputTag>("MuonCollection"));
  thePfJetCollection_ = consumes<reco::PFJetCollection>(ps.getParameter<edm::InputTag>("pfJetCollection"));
  thePfMETCollection_ = consumes<reco::PFMETCollection>(ps.getParameter<edm::InputTag>("pfMETCollection"));
  theCaloJetCollection_ = consumes<reco::CaloJetCollection>(ps.getParameter<edm::InputTag>("caloJetCollection"));
  theCaloMETCollection_ = consumes<reco::CaloMETCollection>(ps.getParameter<edm::InputTag>("caloMETCollection"));
  triggerResults_ = consumes<edm::TriggerResults>(ps.getParameter<edm::InputTag>("TriggerResults"));
  HLTProcess_ = ps.getParameter<std::string>("HLTProcess");
  triggerPath_ = ps.getParameter<std::string>("TriggerPath");
  triggerMuFilter_ = ps.getParameter<edm::InputTag>("TriggerFilterMuon");
  triggerHTFilter_ = ps.getParameter<edm::InputTag>("TriggerFilterHT");
  triggerMetFilter_ = ps.getParameter<edm::InputTag>("TriggerFilterMET");
  triggerDiJetFilter_ = ps.getParameter<edm::InputTag>("TriggerFilterMJJ");
  triggerCaloMETFilter_ = ps.getParameter<edm::InputTag>("TriggerFilterCaloMET");
  ptThrJetTrig_ = ps.getUntrackedParameter<double>("PtThrJetTrig");
  etaThrJetTrig_ = ps.getUntrackedParameter<double>("EtaThrJetTrig");
  ptThrJet_ = ps.getUntrackedParameter<double>("PtThrJet");
  etaThrJet_ = ps.getUntrackedParameter<double>("EtaThrJet");
  deltaetaVBFJets = ps.getUntrackedParameter<double>("DeltaEtaVBFJets");
  pfmetOnlinethreshold = ps.getUntrackedParameter<double>("PFMetCutOnline");
  muonOnlinethreshold = ps.getUntrackedParameter<double>("MuonCutOnline");
  htOnlinethreshold = ps.getUntrackedParameter<double>("HTCutOnline");
  mjjOnlinethreshold = ps.getUntrackedParameter<double>("MJJCutOnline");
}

SUSY_HLT_VBF_Mu::~SUSY_HLT_VBF_Mu() {
  edm::LogInfo("SUSY_HLT_VBF_Mu") << "Destructor SUSY_HLT_VBF_Mu::~SUSY_HLT_VBF_Mu " << std::endl;
}

void SUSY_HLT_VBF_Mu::dqmBeginRun(edm::Run const &run, edm::EventSetup const &e) {
  bool changed;

  if (!fHltConfig.init(run, e, HLTProcess_, changed)) {
    edm::LogError("SUSY_HLT_VBF_Mu") << "Initialization of HLTConfigProvider failed!!";
    return;
  }

  bool pathFound = false;
  const std::vector<std::string> allTrigNames = fHltConfig.triggerNames();
  for (size_t j = 0; j < allTrigNames.size(); ++j) {
    if (allTrigNames[j].find(triggerPath_) != std::string::npos) {
      pathFound = true;
    }
  }

  if (!pathFound) {
    LogDebug("SUSY_HLT_VBF_Mu") << "Path not found"
                                << "\n";
    return;
  }

  edm::LogInfo("SUSY_HLT_VBF_Mu") << "SUSY_HLT_VBF_Mu::beginRun" << std::endl;
}

void SUSY_HLT_VBF_Mu::bookHistograms(DQMStore::IBooker &ibooker_, edm::Run const &, edm::EventSetup const &) {
  edm::LogInfo("SUSY_HLT_VBF_Mu") << "SUSY_HLT_VBF_Mu::bookHistograms" << std::endl;
  // book at beginRun
  bookHistos(ibooker_);
}

void SUSY_HLT_VBF_Mu::analyze(edm::Event const &e, edm::EventSetup const &eSetup) {
  edm::LogInfo("SUSY_HLT_VBF_Mu") << "SUSY_HLT_VBF_Mu::analyze" << std::endl;

  //-------------------------------
  //--- Jets
  //-------------------------------

  edm::Handle<reco::PFJetCollection> pfJetCollection;
  e.getByToken(thePfJetCollection_, pfJetCollection);
  if (!pfJetCollection.isValid()) {
    edm::LogError("SUSY_HLT_VBF_Mu") << "invalid collection: PFJets"
                                     << "\n";
    return;
  }

  edm::Handle<reco::CaloJetCollection> caloJetCollection;
  e.getByToken(theCaloJetCollection_, caloJetCollection);
  if (!caloJetCollection.isValid()) {
    edm::LogError("SUSY_HLT_VBF_Mu") << "invalid collection: CaloJets"
                                     << "\n";
    return;
  }

  //-------------------------------
  //--- Muon
  //-------------------------------
  edm::Handle<reco::MuonCollection> MuonCollection;
  e.getByToken(theMuonCollection_, MuonCollection);
  if (!MuonCollection.isValid()) {
    edm::LogError("SUSY_HLT_VBF_Mu") << "invalid collection: Muons "
                                     << "\n";
    return;
  }

  //-------------------------------
  //--- MET
  //-------------------------------

  edm::Handle<reco::CaloMETCollection> caloMETCollection;
  e.getByToken(theCaloMETCollection_, caloMETCollection);
  if (!caloMETCollection.isValid()) {
    edm::LogError("SUSY_HLT_VBF_Mu") << "invalid collection: CaloMET"
                                     << "\n";
    return;
  }

  edm::Handle<reco::PFMETCollection> pfMETCollection;
  e.getByToken(thePfMETCollection_, pfMETCollection);
  if (!pfMETCollection.isValid()) {
    edm::LogError("SUSY_HLT_VBF_Mu") << "invalid collection: PFMET"
                                     << "\n";
    return;
  }
  //

  //-------------------------------
  //--- Trigger
  //-------------------------------
  edm::Handle<edm::TriggerResults> hltresults;
  e.getByToken(triggerResults_, hltresults);
  if (!hltresults.isValid()) {
    edm::LogError("SUSY_HLT_VBF_Mu") << "invalid collection: TriggerResults"
                                     << "\n";
    return;
  }
  edm::Handle<trigger::TriggerEvent> triggerSummary;
  e.getByToken(theTrigSummary_, triggerSummary);
  if (!triggerSummary.isValid()) {
    edm::LogError("SUSY_HLT_VBF_Mu") << "invalid collection: TriggerSummary"
                                     << "\n";
    return;
  }

  // get online objects

  std::vector<float> ptMuon, etaMuon, phiMuon;
  size_t filterIndexMu = triggerSummary->filterIndex(triggerMuFilter_);
  trigger::TriggerObjectCollection triggerObjects = triggerSummary->getObjects();
  if (!(filterIndexMu >= triggerSummary->sizeFilters())) {
    const trigger::Keys &keys = triggerSummary->filterKeys(filterIndexMu);
    for (size_t j = 0; j < keys.size(); ++j) {
      trigger::TriggerObject foundObject = triggerObjects[keys[j]];
      if (fabs(foundObject.id()) == 13) {  // It's a muon

        bool same = false;
        for (unsigned int x = 0; x < ptMuon.size(); x++) {
          if (fabs(ptMuon[x] - foundObject.pt()) < 0.01)
            same = true;
        }

        if (!same) {
          h_triggerMuPt->Fill(foundObject.pt());
          h_triggerMuEta->Fill(foundObject.eta());
          h_triggerMuPhi->Fill(foundObject.phi());
          ptMuon.push_back(foundObject.pt());
          etaMuon.push_back(foundObject.eta());
          phiMuon.push_back(foundObject.phi());
        }
      }
    }
  }

  // get online objects

  size_t filterIndexMet = triggerSummary->filterIndex(triggerMetFilter_);
  size_t dijetFilterIndex = triggerSummary->filterIndex(triggerDiJetFilter_);

  if (!(filterIndexMet >= triggerSummary->sizeFilters())) {
    const trigger::Keys &keys = triggerSummary->filterKeys(filterIndexMet);
    for (size_t j = 0; j < keys.size(); ++j) {
      trigger::TriggerObject foundObject = triggerObjects[keys[j]];
      h_triggerMet->Fill(foundObject.pt());
      h_triggerMetPhi->Fill(foundObject.phi());
    }
  }

  dijet = -1;

  std::vector<float> ptJet, etaJet, phiJet;
  if (!(dijetFilterIndex >= triggerSummary->sizeFilters())) {
    const trigger::Keys &KEYS(triggerSummary->filterKeys(dijetFilterIndex));
    const size_t nK(KEYS.size());
    const trigger::TriggerObjectCollection &TOC(triggerSummary->getObjects());

    for (size_t i = 0; i < nK; i++) {
      const trigger::TriggerObject &TO1(TOC[KEYS[i]]);
      if (TO1.pt() > ptThrJetTrig_ && fabs(TO1.eta()) < etaThrJetTrig_) {
        // for dijet part
        for (size_t j = i; j < nK; j++) {
          if (i < j) {
            const trigger::TriggerObject &TO2(TOC[KEYS[j]]);
            if (TO2.pt() > ptThrJetTrig_ && fabs(TO2.eta()) < etaThrJetTrig_) {
              double tmpdeta = fabs(TO1.eta() - TO2.eta());
              double tmpopposite = TO1.eta() * TO2.eta();
              if (tmpdeta > deltaetaVBFJets && tmpopposite < 0) {
                TLorentzVector j1(TO1.px(), TO1.py(), TO1.pz(), TO1.energy());
                TLorentzVector j2(TO2.px(), TO2.py(), TO2.pz(), TO2.energy());
                double tmpmass = (j1 + j2).M();
                if (dijet < tmpmass) {
                  dijet = tmpmass;
                }
              }
            }
          }
        }
      }
    }
    h_DiJetMass->Fill(dijet);
  }

  size_t filterIndexCaloMET = triggerSummary->filterIndex(triggerCaloMETFilter_);
  if (filterIndexCaloMET < triggerSummary->sizeFilters()) {
    const trigger::Keys &keys = triggerSummary->filterKeys(filterIndexCaloMET);
    if (!keys.empty()) {
      float met_h = triggerObjects[keys[0]].pt();
      h_triggerCaloMet->Fill(met_h);
    }
  }

  size_t filterIndexHt = triggerSummary->filterIndex(triggerHTFilter_);
  if (filterIndexHt < triggerSummary->sizeFilters()) {
    const trigger::Keys &keys = triggerSummary->filterKeys(filterIndexHt);
    if (!keys.empty()) {
      float ht = triggerObjects[keys[0]].pt();
      h_ht->Fill(ht);
    }
  }

  bool hasFired = false;

  const edm::TriggerNames &trigNames = e.triggerNames(*hltresults);
  unsigned int numTriggers = trigNames.size();

  for (unsigned int hltIndex = 0; hltIndex < numTriggers; ++hltIndex) {
    if (trigNames.triggerName(hltIndex).find(triggerPath_) != std::string::npos && hltresults->wasrun(hltIndex) &&
        hltresults->accept(hltIndex))
      hasFired = true;
  }

  // Matching the muon
  int indexOfMatchedMuon = -1;
  int offlineCounter = 0;
  for (reco::MuonCollection::const_iterator muon = MuonCollection->begin(); muon != MuonCollection->end(); ++muon) {
    for (size_t off_i = 0; off_i < ptMuon.size(); ++off_i) {
      if (reco::deltaR(muon->eta(), muon->phi(), etaMuon[off_i], phiMuon[off_i]) < 0.5) {
        indexOfMatchedMuon = offlineCounter;
        break;
      }
    }
    offlineCounter++;
  }

  float pfHT = 0.0;

  for (reco::PFJetCollection::const_iterator i_pfjet = pfJetCollection->begin(); i_pfjet != pfJetCollection->end();
       ++i_pfjet) {
    if (i_pfjet->pt() < ptThrJet_)
      continue;
    if (fabs(i_pfjet->eta()) > etaThrJet_)
      continue;
    pfHT += i_pfjet->pt();
  }

  //

  dijetOff = -1;

  size_t jetCol = pfJetCollection->size();

  for (size_t i = 0; i < jetCol; i++) {
    if (pfJetCollection->at(i).pt() > ptThrJetTrig_ && fabs(pfJetCollection->at(i).eta()) < etaThrJetTrig_) {
      for (size_t j = i; j < jetCol; j++) {
        if (i < j) {
          if (pfJetCollection->at(j).pt() > ptThrJetTrig_ && fabs(pfJetCollection->at(j).eta()) < etaThrJetTrig_) {
            double tmpdetaOff = fabs(pfJetCollection->at(i).eta() - pfJetCollection->at(j).eta());
            double tmpoppositeOff = pfJetCollection->at(i).eta() * pfJetCollection->at(j).eta();
            if (tmpdetaOff > deltaetaVBFJets && tmpoppositeOff < 0) {
              TLorentzVector j1Off(pfJetCollection->at(i).px(),
                                   pfJetCollection->at(i).py(),
                                   pfJetCollection->at(i).pz(),
                                   pfJetCollection->at(i).energy());
              TLorentzVector j2Off(pfJetCollection->at(j).px(),
                                   pfJetCollection->at(j).py(),
                                   pfJetCollection->at(j).pz(),
                                   pfJetCollection->at(j).energy());
              double tmpmassOff = (j1Off + j2Off).M();
              if (dijetOff < tmpmassOff)
                dijetOff = tmpmassOff;
            }
          }
        }
      }
    }
  }

  //  For trigger turn on curves

  // for muon
  if (indexOfMatchedMuon > -1 && (dijetOff > mjjOnlinethreshold) &&
      (pfMETCollection->begin()->et() > pfmetOnlinethreshold) && (pfHT > htOnlinethreshold)) {
    h_den_muonpt->Fill(MuonCollection->at(indexOfMatchedMuon).pt());
    h_den_muoneta->Fill(MuonCollection->at(indexOfMatchedMuon).eta());
    if (hasFired) {
      h_num_muonpt->Fill(MuonCollection->at(indexOfMatchedMuon).pt());
      h_num_muoneta->Fill(MuonCollection->at(indexOfMatchedMuon).eta());
    }
  }

  // For MJJ
  if (indexOfMatchedMuon > -1 && (MuonCollection->at(indexOfMatchedMuon).pt() > muonOnlinethreshold) &&
      (pfMETCollection->begin()->et() > pfmetOnlinethreshold) && (pfHT > htOnlinethreshold)) {
    h_den_mjj->Fill(dijetOff);
    if (hasFired) {
      h_num_mjj->Fill(dijetOff);
    }
  }

  // For HT
  //
  if (indexOfMatchedMuon > -1 && (MuonCollection->at(indexOfMatchedMuon).pt() > muonOnlinethreshold) &&
      (pfMETCollection->begin()->et() > pfmetOnlinethreshold) && (dijetOff > mjjOnlinethreshold)) {
    h_den_ht->Fill(pfHT);
    if (hasFired) {
      h_num_ht->Fill(pfHT);
    }
  }

  if (indexOfMatchedMuon > -1 && (dijetOff > mjjOnlinethreshold) && (pfHT > htOnlinethreshold) &&
      (MuonCollection->at(indexOfMatchedMuon).pt() > muonOnlinethreshold)) {
    h_den_met->Fill(pfMETCollection->begin()->et());
    if (hasFired) {
      h_num_met->Fill(pfMETCollection->begin()->et());
    }
  }
}

void SUSY_HLT_VBF_Mu::endRun(edm::Run const &run, edm::EventSetup const &eSetup) {
  edm::LogInfo("SUSY_HLT_VBF_Mu") << "SUSY_HLT_VBF_Mu::endRun" << std::endl;
}

void SUSY_HLT_VBF_Mu::bookHistos(DQMStore::IBooker &ibooker_) {
  ibooker_.cd();
  ibooker_.setCurrentFolder("HLT/SUSYBSM/SUSY_HLT_VBF_Mu");

  // offline quantities

  // online quantities
  h_triggerMuPt = ibooker_.book1D("triggerMuPt", "Trigger Muon Pt; GeV", 50, 0.0, 500.0);
  h_triggerMuEta = ibooker_.book1D("triggerMuEta", "Trigger Muon Eta", 20, -3.0, 3.0);
  h_triggerMuPhi = ibooker_.book1D("triggerMuPhi", "Trigger Muon Phi", 20, -3.5, 3.5);

  h_triggerCaloMet = ibooker_.book1D("h_triggerCaloMet", "Trigger Calo MET; GeV", 20, 0.0, 500.0);

  h_ht = ibooker_.book1D("h_ht", "h_ht", 30, 0.0, 1500.0);

  h_triggerMet = ibooker_.book1D("triggerMet", "Trigger MET; GeV", 20, 0.0, 500.0);
  h_triggerMetPhi = ibooker_.book1D("triggerMetPhi", "Trigger MET Phi", 20, -3.5, 3.5);

  h_DiJetMass = ibooker_.book1D("h_DiJetMass", "h_DiJetMass", 500, 0.0, 5000.0);

  // num and den hists to be divided in harvesting step to make turn on curves

  h_den_muonpt = ibooker_.book1D("h_den_muonpt", "h_den_muonpt", 50, 0.0, 500.0);
  h_num_muonpt = ibooker_.book1D("h_num_muonpt", "h_num_muonpt", 50, 0.0, 500.0);

  h_den_muoneta = ibooker_.book1D("h_den_muoneta", "h_den_muoneta", 20, -3.0, 3.0);
  h_num_muoneta = ibooker_.book1D("h_num_muoneta", "h_num_muoneta", 20, -3.0, 3.0);

  h_den_mjj = ibooker_.book1D("h_den_mjj", "h_den_mjj", 500, 0.0, 5000.0);
  h_num_mjj = ibooker_.book1D("h_num_mjj", "h_num_mjj", 500, 0.0, 5000.0);

  h_den_met = ibooker_.book1D("h_den_met", "h_den_met", 20, 0.0, 500.0);
  h_num_met = ibooker_.book1D("h_num_met", "h_num_met", 20, 0.0, 500.0);

  h_den_ht = ibooker_.book1D("h_den_ht", "h_den_ht", 30, 0.0, 1500.0);
  h_num_ht = ibooker_.book1D("h_num_ht", "h_num_ht", 30, 0.0, 1500.0);

  ibooker_.cd();
}

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