L1Muon2RecoTreeProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    L1Trigger/L1TNtuples
// Class:      L1JetRecoTreeProducer
//
// system include files
#include <memory>

// framework
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Common/interface/TriggerNames.h"

// cond formats
#include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"

// data formats
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/JetReco/interface/PFJetCollection.h"
#include "DataFormats/JetReco/interface/JetID.h"
#include "DataFormats/METReco/interface/PFMET.h"
#include "DataFormats/METReco/interface/PFMETCollection.h"

// RECO TRIGGER MATCHING:
#include "DataFormats/Math/interface/deltaR.h"
#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/HLTReco/interface/TriggerEvent.h"
#include "DataFormats/HLTReco/interface/TriggerObject.h"
#include "TString.h"
#include "TRegexp.h"
#include <utility>

//muons
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/GeometrySurface/interface/Cylinder.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "DataFormats/MuonReco/interface/MuonEnergy.h"
#include "DataFormats/MuonReco/interface/MuonTime.h"
#include "CondFormats/AlignmentRecord/interface/TrackerSurfaceDeformationRcd.h"

//taus
#include "DataFormats/TauReco/interface/PFTau.h"
#include "DataFormats/TauReco/interface/PFTauFwd.h"
#include "DataFormats/TauReco/interface/PFTauDiscriminator.h"

// ROOT output stuff
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "TH1.h"
#include "TTree.h"
#include "TF1.h"

//local  data formats
#include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoMuon2.h"
#include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoMet.h"

//vertices
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoVertexDataFormat.h"

using namespace std;

//
// class declaration
//

class L1Muon2RecoTreeProducer : public edm::one::EDAnalyzer<edm::one::SharedResources, edm::one::WatchRuns> {
public:
  explicit L1Muon2RecoTreeProducer(const edm::ParameterSet &);
  ~L1Muon2RecoTreeProducer() override;

private:
  void beginJob() override;
  void analyze(const edm::Event &, const edm::EventSetup &) override;
  void beginRun(const edm::Run &, const edm::EventSetup &) override;
  void endRun(edm::Run const &, edm::EventSetup const &) override {}
  void endJob() override;

public:
  L1Analysis::L1AnalysisRecoMuon2 *muon;

  L1Analysis::L1AnalysisRecoMuon2DataFormat *muon_data;

  double match_trigger(std::vector<int> &trigIndices,
                       const trigger::TriggerObjectCollection &trigObjs,
                       edm::Handle<trigger::TriggerEvent> &triggerEvent,
                       const reco::Muon &mu);
  void empty_hlt();

private:
  // output file
  edm::Service<TFileService> fs_;

  // tree
  TTree *tree_;

  // EDM input tags
  edm::EDGetTokenT<reco::MuonCollection> MuonToken_;
  edm::EDGetTokenT<edm::TriggerResults> TriggerResultsToken_;
  edm::EDGetTokenT<trigger::TriggerEvent> triggerSummaryLabelToken_;
  edm::EDGetTokenT<reco::VertexCollection> VtxToken_;
  edm::EDGetTokenT<reco::PFMETCollection> metToken_;

  // bool triggerMatching_;
  // edm::InputTag triggerSummaryLabel_;
  // std::string triggerProcessLabel_;
  // std::vector<std::string> isoTriggerNames_;
  // std::vector<std::string> triggerNames_;

  edm::Handle<edm::TriggerResults> isoTriggerToken_;
  edm::Handle<std::vector<std::string>> isoTriggerNamesToken_;
  // edm::Handle<edm::TriggerResults>     TriggerToken_;

  HLTConfigProvider hltConfig_;

  // debug stuff
  unsigned int maxMuon_;
  double triggerMaxDeltaR_;
  bool triggerMatching_;
  std::string triggerProcessLabel_;
  std::vector<std::string> isoTriggerNames_;
  std::vector<std::string> triggerNames_;
  std::vector<int> isoTriggerIndices_;
  std::vector<int> triggerIndices_;
};

L1Muon2RecoTreeProducer::L1Muon2RecoTreeProducer(const edm::ParameterSet &iConfig) {
  usesResource(TFileService::kSharedResource);

  maxMuon_ = iConfig.getParameter<unsigned int>("maxMuon");
  isoTriggerNames_ = iConfig.getParameter<std::vector<std::string>>("isoTriggerNames");
  // isoTriggerToken_         = iConfig.getParameter<std::vector<std::string>>("isoTriggerNames");
  // TriggerToken_         = iConfig.getParameter<std::vector<std::string>>("triggerNames");
  MuonToken_ = consumes<reco::MuonCollection>(iConfig.getUntrackedParameter("MuonToken", edm::InputTag("muons")));
  VtxToken_ = consumes<reco::VertexCollection>(
      iConfig.getUntrackedParameter("VertexToken", edm::InputTag("offlinePrimaryVertices")));

  TriggerResultsToken_ = consumes<edm::TriggerResults>(
      iConfig.getUntrackedParameter("TriggerResultsToken", edm::InputTag("TriggerResults", "", "HLT")));
  triggerSummaryLabelToken_ = consumes<trigger::TriggerEvent>(
      iConfig.getUntrackedParameter("triggerSummaryLabelToken", edm::InputTag("hltTriggerSummaryAOD", "", "HLT")));
  // TriggerResultsToken_      = consumes<edm::TriggerResults>(iConfig.getUntrackedParameter("TriggerResultsToken",edm::InputTag("triggerSummaryLabel")));
  // triggerSummaryLabelToken_ = consumes<trigger::TriggerEvent>(iConfig.getUntrackedParameter("triggerSummaryLabelToken",edm::InputTag("triggerSummaryLabel")));
  //iConfig.getParameter<edm::InputTag>("triggerSummaryLabel");
  metToken_ = consumes<reco::PFMETCollection>(iConfig.getUntrackedParameter("metToken", edm::InputTag("pfMet")));

  muon = new L1Analysis::L1AnalysisRecoMuon2(iConfig, consumesCollector());
  muon_data = muon->getData();

  tree_ = fs_->make<TTree>("Muon2RecoTree", "Muon2RecoTree");
  tree_->Branch("Muon", "L1Analysis::L1AnalysisRecoMuon2DataFormat", &muon_data, 32000, 3);

  triggerMaxDeltaR_ = iConfig.getParameter<double>("triggerMaxDeltaR");
  triggerMatching_ = iConfig.getUntrackedParameter<bool>("triggerMatching");
  triggerProcessLabel_ = iConfig.getUntrackedParameter<std::string>("triggerProcessLabel");
  isoTriggerNames_ = iConfig.getParameter<std::vector<std::string>>("isoTriggerNames");
  triggerNames_ = iConfig.getParameter<std::vector<std::string>>("triggerNames");

  // triggerMatching_     = iConfig.getUntrackedParameter<bool>("triggerMatching");
  // _ = iConfig.getParameter<edm::InputTag>("triggerSummaryLabel");
  // triggerProcessLabel_ = iConfig.getUntrackedParameter<std::string>("triggerProcessLabel");
  // triggerNames_        = iConfig.getParameter<std::vector<std::string> > ("triggerNames");
  // isoTriggerNames_     = iConfig.getParameter<std::vector<std::string> > ("isoTriggerNames");
  // triggerMaxDeltaR_    = iConfig.getParameter<double> ("triggerMaxDeltaR");
}

L1Muon2RecoTreeProducer::~L1Muon2RecoTreeProducer() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}

//
// member functions
//

// ------------ method called to for each event  ------------
void L1Muon2RecoTreeProducer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
  muon->init(iSetup);

  muon->Reset();
  edm::Handle<reco::MuonCollection> recoMuons;
  iEvent.getByToken(MuonToken_, recoMuons);

  edm::Handle<reco::VertexCollection> vertices;
  iEvent.getByToken(VtxToken_, vertices);

  edm::Handle<edm::TriggerResults> triggerResults;
  iEvent.getByToken(TriggerResultsToken_, triggerResults);

  edm::Handle<trigger::TriggerEvent> triggerSummaryLabel_;
  iEvent.getByToken(triggerSummaryLabelToken_, triggerSummaryLabel_);

  edm::Handle<reco::PFMETCollection> metLabel_;
  iEvent.getByToken(metToken_, metLabel_);

  int counter_met = 0;

  double METx = 0.;
  double METy = 0.;

  for (reco::PFMETCollection::const_iterator imet = metLabel_->begin();
       imet != metLabel_->end() && (unsigned)counter_met < 1;
       imet++) {
    METx = imet->px();
    METy = imet->py();
  }

  if (recoMuons.isValid()) {
    muon->SetMuon(iEvent, iSetup, recoMuons, vertices, METx, METy, maxMuon_);
  } else {
  }

  int counter_mu = 0;
  for (reco::MuonCollection::const_iterator imu = recoMuons->begin();
       imu != recoMuons->end() && (unsigned)counter_mu < maxMuon_;
       imu++) {
    //---------------------------------------------------------------------
    // TRIGGER MATCHING:
    // if specified the reconstructed muons are matched to a trigger
    //---------------------------------------------------------------------
    if (triggerMatching_) {
      double isoMatchDeltaR = 9999.;
      double matchDeltaR = 9999.;
      int hasIsoTriggered = 0;
      int hasTriggered = 0;

      int passesSingleMuonFlag = 0;

      // first check if the trigger results are valid:
      if (triggerResults.isValid()) {
        if (triggerSummaryLabel_.isValid()) {
          const edm::TriggerNames &trigNames = iEvent.triggerNames(*triggerResults);
          // for(UInt_t iPath = 0 ; iPath < trigNames.size() ; ++iPath)
          //   {
          //     cout<<iPath<<": "<<trigNames.triggerName(iPath)<<endl;
          //   }

          for (UInt_t iPath = 0; iPath < isoTriggerNames_.size(); ++iPath) {
            if (passesSingleMuonFlag == 1)
              continue;
            std::string pathName = isoTriggerNames_.at(iPath);

            bool passTrig = false;
            //cout<<"testing pathName                         = "<<pathName<<endl;
            //cout<<"trigNames.triggerIndex(pathName) = "<<trigNames.triggerIndex(pathName)<<endl;
            // cout<<"trigNames.triggerIndex(pathName)<trigNames.size()= "<<(trigNames.triggerIndex(pathName)<trigNames.size())<<endl;

            if (trigNames.triggerIndex(pathName) < trigNames.size())
              passTrig = triggerResults->accept(trigNames.triggerIndex(pathName));
            // cout<<"pass = "<<passTrig<<endl;
            if (passTrig)
              passesSingleMuonFlag = 1;
          }

          if (triggerSummaryLabel_.isValid()) {
            // get trigger objects:
            const trigger::TriggerObjectCollection triggerObjects = triggerSummaryLabel_->getObjects();

            matchDeltaR = match_trigger(triggerIndices_, triggerObjects, triggerSummaryLabel_, (*imu));
            if (matchDeltaR < triggerMaxDeltaR_)
              hasTriggered = 1;

            // cout<<"isoTriggerIndices_.size() = "<<isoTriggerIndices_.size()<<endl;
            // cout<<"triggerObjects.size() = "<<triggerObjects.size()<<endl;
            // cout<<"imu->eta() = "<<imu->eta()<<endl;

            isoMatchDeltaR = match_trigger(isoTriggerIndices_, triggerObjects, triggerSummaryLabel_, (*imu));
            if (isoMatchDeltaR < triggerMaxDeltaR_)
              hasIsoTriggered = 1;

          }  // end if (triggerEvent.isValid())

        }  // end if (triggerResults.isValid())

      }  // end if (triggerResults.isValid())

      // fill trigger matching variables:
      muon_data->hlt_isomu.push_back(hasIsoTriggered);
      muon_data->hlt_mu.push_back(hasTriggered);
      muon_data->hlt_isoDeltaR.push_back(isoMatchDeltaR);
      muon_data->hlt_deltaR.push_back(matchDeltaR);
      muon_data->passesSingleMuon.push_back(passesSingleMuonFlag);

    } else {
      empty_hlt();
    }  // end if (triggerMatching_)
  }

  tree_->Fill();
}

// ------------ method called once each job just before starting event loop  ------------
void L1Muon2RecoTreeProducer::beginJob(void) {}

// ------------ method called once each job just after ending the event loop  ------------
void L1Muon2RecoTreeProducer::endJob() {}

void L1Muon2RecoTreeProducer::empty_hlt() {
  muon_data->hlt_isomu.push_back(-9999);
  muon_data->hlt_mu.push_back(-9999);
  muon_data->hlt_isoDeltaR.push_back(-9999);
  muon_data->hlt_deltaR.push_back(-9999);
}

double L1Muon2RecoTreeProducer::match_trigger(std::vector<int> &trigIndices,
                                              const trigger::TriggerObjectCollection &trigObjs,
                                              edm::Handle<trigger::TriggerEvent> &triggerEvent,
                                              const reco::Muon &mu) {
  double matchDeltaR = 9999;

  for (size_t iTrigIndex = 0; iTrigIndex < trigIndices.size(); ++iTrigIndex) {
    int triggerIndex = trigIndices[iTrigIndex];
    const std::vector<std::string> moduleLabels(hltConfig_.moduleLabels(triggerIndex));
    // find index of the last module:
    const unsigned moduleIndex = hltConfig_.size(triggerIndex) - 2;
    // find index of HLT trigger name:
    const unsigned hltFilterIndex =
        triggerEvent->filterIndex(edm::InputTag(moduleLabels[moduleIndex], "", triggerProcessLabel_));

    if (hltFilterIndex < triggerEvent->sizeFilters()) {
      const trigger::Keys triggerKeys(triggerEvent->filterKeys(hltFilterIndex));
      const trigger::Vids triggerVids(triggerEvent->filterIds(hltFilterIndex));

      const unsigned nTriggers = triggerVids.size();
      for (size_t iTrig = 0; iTrig < nTriggers; ++iTrig) {
        // loop over all trigger objects:
        const trigger::TriggerObject trigObject = trigObjs[triggerKeys[iTrig]];

        double dRtmp = deltaR(mu, trigObject);

        if (dRtmp < matchDeltaR) {
          matchDeltaR = dRtmp;
        }

      }  // loop over different trigger objects
    }    // if trigger is in event (should apply hltFilter with used trigger...)
  }      // loop over muon candidates

  return matchDeltaR;
}

void L1Muon2RecoTreeProducer::beginRun(const edm::Run &run, const edm::EventSetup &eventSetup) {
  // Prepare for trigger matching for each new run:
  // Look up triggetIndices in the HLT config for the different paths
  if (triggerMatching_) {
    bool changed = true;
    if (!hltConfig_.init(run, eventSetup, triggerProcessLabel_, changed)) {
      // if you can't initialize hlt configuration, crash!
      std::cout << "Error: didn't find process" << triggerProcessLabel_ << std::endl;
      assert(false);
    }

    bool enableWildcard = true;
    for (size_t iTrig = 0; iTrig < triggerNames_.size(); ++iTrig) {
      // prepare for regular expression (with wildcards) functionality:
      TString tNameTmp = TString(triggerNames_[iTrig]);
      TRegexp tNamePattern = TRegexp(tNameTmp, enableWildcard);
      int tIndex = -1;
      // find the trigger index:
      for (unsigned ipath = 0; ipath < hltConfig_.size(); ++ipath) {
        // use TString since it provides reg exp functionality:
        TString tmpName = TString(hltConfig_.triggerName(ipath));
        if (tmpName.Contains(tNamePattern)) {
          tIndex = int(ipath);
          triggerIndices_.push_back(tIndex);
        }
      }
      if (tIndex < 0) {  // if can't find trigger path at all, give warning:
        std::cout << "Warning: Could not find trigger" << triggerNames_[iTrig] << std::endl;
        //assert(false);
      }
    }  // end for triggerNames
    for (size_t iTrig = 0; iTrig < isoTriggerNames_.size(); ++iTrig) {
      // prepare for regular expression functionality:
      TString tNameTmp = TString(isoTriggerNames_[iTrig]);
      TRegexp tNamePattern = TRegexp(tNameTmp, enableWildcard);
      int tIndex = -1;
      // find the trigger index:
      for (unsigned ipath = 0; ipath < hltConfig_.size(); ++ipath) {
        // use TString since it provides reg exp functionality:
        TString tmpName = TString(hltConfig_.triggerName(ipath));
        if (tmpName.Contains(tNamePattern)) {
          tIndex = int(ipath);
          isoTriggerIndices_.push_back(tIndex);
        }
      }
      if (tIndex < 0) {  // if can't find trigger path at all, give warning:
        std::cout << "Warning: Could not find trigger" << isoTriggerNames_[iTrig] << std::endl;
        //assert(false);
      }
    }  // end for isoTriggerNames
  }    // end if (triggerMatching_)
}

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