L1ECALPrefiringWeightProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    ProdTutorial/L1ECALPrefiringWeightProducer
// Class:      L1ECALPrefiringWeightProducer
//
//
// Original Author:  localusers user
//         Created:  Thu, 08 Nov 2018 16:16:00 GMT
//
//
 
// system include files
#include <memory>
 
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
 
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
 
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/StreamID.h"
 
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/Photon.h"
#include "TH2.h"
 
#include <iostream>
enum fluctuations { central = 0, up, down };
 
class L1ECALPrefiringWeightProducer : public edm::global::EDProducer<> {
public:
  explicit L1ECALPrefiringWeightProducer(const edm::ParameterSet&);
  ~L1ECALPrefiringWeightProducer() override;
 
  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
 
private:
  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
 
  double getPrefiringRate(double eta, double pt, TH2F* h_prefmap, fluctuations fluctuation) const;
 
  edm::EDGetTokenT<std::vector<pat::Photon> > photons_token_;
  edm::EDGetTokenT<std::vector<pat::Jet> > jets_token_;
 
  TH2F* h_prefmap_photon;
  TH2F* h_prefmap_jet;
  std::string dataera_;
  bool useEMpt_;
  double prefiringRateSystUnc_;
  bool skipwarnings_;
};
 
L1ECALPrefiringWeightProducer::L1ECALPrefiringWeightProducer(const edm::ParameterSet& iConfig) {
  photons_token_ = consumes<std::vector<pat::Photon> >(iConfig.getParameter<edm::InputTag>("ThePhotons"));
  jets_token_ = consumes<std::vector<pat::Jet> >(iConfig.getParameter<edm::InputTag>("TheJets"));
 
  dataera_ = iConfig.getParameter<std::string>("DataEra");
  useEMpt_ = iConfig.getParameter<bool>("UseJetEMPt");
  prefiringRateSystUnc_ = iConfig.getParameter<double>("PrefiringRateSystematicUncty");
  skipwarnings_ = iConfig.getParameter<bool>("SkipWarnings");
 
  TFile* file_prefiringmaps_;
  std::string fname = iConfig.getParameter<std::string>("L1Maps");
  edm::FileInPath mapsfilepath("PhysicsTools/PatUtils/data/" + fname);
  file_prefiringmaps_ = new TFile(mapsfilepath.fullPath().c_str(), "read");
  if (file_prefiringmaps_ == nullptr && !skipwarnings_)
    std::cout << "File with maps not found. All prefiring weights set to 0. " << std::endl;
 
  TString mapphotonfullname = "L1prefiring_photonptvseta_" + dataera_;
  if (!file_prefiringmaps_->Get(mapphotonfullname) && !skipwarnings_)
    std::cout << "Photon map not found. All photons prefiring weights set to 0. " << std::endl;
  h_prefmap_photon = (TH2F*)file_prefiringmaps_->Get(mapphotonfullname);
 
  TString mapjetfullname = (useEMpt_) ? "L1prefiring_jetemptvseta_" + dataera_ : "L1prefiring_jetptvseta_" + dataera_;
  if (!file_prefiringmaps_->Get(mapjetfullname) && !skipwarnings_)
    std::cout << "Jet map not found. All jets prefiring weights set to 0. " << std::endl;
  h_prefmap_jet = (TH2F*)file_prefiringmaps_->Get(mapjetfullname);
  file_prefiringmaps_->Close();
  delete file_prefiringmaps_;
  produces<double>("nonPrefiringProb").setBranchAlias("nonPrefiringProb");
  produces<double>("nonPrefiringProbUp").setBranchAlias("nonPrefiringProbUp");
  produces<double>("nonPrefiringProbDown").setBranchAlias("nonPrefiringProbDown");
}
 
L1ECALPrefiringWeightProducer::~L1ECALPrefiringWeightProducer() {
  delete h_prefmap_photon;
  delete h_prefmap_jet;
}
 
void L1ECALPrefiringWeightProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
  using namespace edm;
 
  //Photons
  edm::Handle<std::vector<pat::Photon> > thePhotons;
  iEvent.getByToken(photons_token_, thePhotons);
 
  //Jets
  edm::Handle<std::vector<pat::Jet> > theJets;
  iEvent.getByToken(jets_token_, theJets);
 
  //Probability for the event NOT to prefire, computed with the prefiring maps per object.
  //Up and down values correspond to the resulting value when shifting up/down all prefiring rates in prefiring maps.
  double nonPrefiringProba[3] = {1., 1., 1.};  //0: central, 1: up, 2: down
 
  for (const auto fluct : {fluctuations::central, fluctuations::up, fluctuations::down}) {
    for (const auto& photon : *thePhotons) {
      double pt_gam = photon.pt();
      double eta_gam = photon.eta();
      if (pt_gam < 20.)
        continue;
      if (fabs(eta_gam) < 2.)
        continue;
      if (fabs(eta_gam) > 3.)
        continue;
      double prefiringprob_gam = getPrefiringRate(eta_gam, pt_gam, h_prefmap_photon, fluct);
      nonPrefiringProba[fluct] *= (1. - prefiringprob_gam);
    }
 
    //Now applying the prefiring maps to jets in the affected regions.
    for (const auto& jet : *theJets) {
      double pt_jet = jet.pt();
      double eta_jet = jet.eta();
      double phi_jet = jet.phi();
      if (pt_jet < 20.)
        continue;
      if (fabs(eta_jet) < 2.)
        continue;
      if (fabs(eta_jet) > 3.)
        continue;
 
      //Loop over photons to remove overlap
      double nonprefiringprobfromoverlappingphotons = 1.;
      for (const auto& photon : *thePhotons) {
        double pt_gam = photon.pt();
        double eta_gam = photon.eta();
        double phi_gam = photon.phi();
        if (pt_gam < 20.)
          continue;
        if (fabs(eta_gam) < 2.)
          continue;
        if (fabs(eta_gam) > 3.)
          continue;
        double dR = reco::deltaR(eta_jet, phi_jet, eta_gam, phi_gam);
        if (dR > 0.4)
          continue;
        double prefiringprob_gam = getPrefiringRate(eta_gam, pt_gam, h_prefmap_photon, fluct);
        nonprefiringprobfromoverlappingphotons *= (1. - prefiringprob_gam);
      }
      //useEMpt =true if one wants to use maps parametrized vs Jet EM pt instead of pt.
      if (useEMpt_)
        pt_jet *= (jet.neutralEmEnergyFraction() + jet.chargedEmEnergyFraction());
      double nonprefiringprobfromoverlappingjet = 1. - getPrefiringRate(eta_jet, pt_jet, h_prefmap_jet, fluct);
 
      if (nonprefiringprobfromoverlappingphotons == 1.) {
        nonPrefiringProba[fluct] *= nonprefiringprobfromoverlappingjet;
      }
      //If overlapping photons have a non prefiring rate larger than the jet, then replace these weights by the jet one
      else if (nonprefiringprobfromoverlappingphotons > nonprefiringprobfromoverlappingjet) {
        if (nonprefiringprobfromoverlappingphotons != 0.) {
          nonPrefiringProba[fluct] *= nonprefiringprobfromoverlappingjet / nonprefiringprobfromoverlappingphotons;
        } else {
          nonPrefiringProba[fluct] = 0.;
        }
      }
      //Last case: if overlapping photons have a non prefiring rate smaller than the jet, don't consider the jet in the event weight, and do nothing.
    }
  }
 
  auto nonPrefiringProb = std::make_unique<double>(nonPrefiringProba[0]);
  auto nonPrefiringProbUp = std::make_unique<double>(nonPrefiringProba[1]);
  auto nonPrefiringProbDown = std::make_unique<double>(nonPrefiringProba[2]);
  iEvent.put(std::move(nonPrefiringProb), "nonPrefiringProb");
  iEvent.put(std::move(nonPrefiringProbUp), "nonPrefiringProbUp");
  iEvent.put(std::move(nonPrefiringProbDown), "nonPrefiringProbDown");
}
 
double L1ECALPrefiringWeightProducer::getPrefiringRate(double eta,
                                                       double pt,
                                                       TH2F* h_prefmap,
                                                       fluctuations fluctuation) const {
  if (h_prefmap == nullptr && !skipwarnings_)
    std::cout << "Prefiring map not found, setting prefiring rate to 0 " << std::endl;
  if (h_prefmap == nullptr)
    return 0.;
  //Check pt is not above map overflow
  int nbinsy = h_prefmap->GetNbinsY();
  double maxy = h_prefmap->GetYaxis()->GetBinLowEdge(nbinsy + 1);
  if (pt >= maxy)
    pt = maxy - 0.01;
  int thebin = h_prefmap->FindBin(eta, pt);
 
  double prefrate = h_prefmap->GetBinContent(thebin);
 
  double statuncty = h_prefmap->GetBinError(thebin);
  double systuncty = prefiringRateSystUnc_ * prefrate;
 
  if (fluctuation == up)
    prefrate = std::min(1., prefrate + sqrt(pow(statuncty, 2) + pow(systuncty, 2)));
  if (fluctuation == down)
    prefrate = std::max(0., prefrate - sqrt(pow(statuncty, 2) + pow(systuncty, 2)));
  return prefrate;
}
 
// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void L1ECALPrefiringWeightProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
 
  desc.add<edm::InputTag>("ThePhotons", edm::InputTag("slimmedPhotons"));
  desc.add<edm::InputTag>("TheJets", edm::InputTag("slimmedJets"));
  desc.add<std::string>("L1Maps", "L1PrefiringMaps.root");
  desc.add<std::string>("DataEra", "2017BtoF");
  desc.add<bool>("UseJetEMPt", false);
  desc.add<double>("PrefiringRateSystematicUncty", 0.2);
  desc.add<bool>("SkipWarnings", true);
  descriptions.add("l1ECALPrefiringWeightProducer", desc);
}
 
//define this as a plug-in
DEFINE_FWK_MODULE(L1ECALPrefiringWeightProducer);