HLTriggerOffline/JetMET/src/JetMETDQMPostProcessor.cc - master - githubmirror/cmssw

HLTriggerOffline/JetMET/src/JetMETDQMPostProcessor.cc
221 строка8 KB
Bruno Alves
Add turn-on plots.
24 сен 2025, 19:55
24 сен 2025, 19:55d7ec10f
// Numerator and denominator histograms filled in HLTriggerOffline/JetMET/src/HLTJetMETValidation.cc

#include "HLTriggerOffline/JetMET/interface/JetMETDQMPostProcessor.h"

#include "DQMServices/Core/interface/DQMStore.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include <cmath>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>

JetMETDQMPostProcessor::JetMETDQMPostProcessor(const edm::ParameterSet &pset) {
  subDir_ = pset.getUntrackedParameter<std::string>("subDir");
  patternJetTrg_ = pset.getUntrackedParameter<std::string>("PatternJetTrg", "");
  patternMetTrg_ = pset.getUntrackedParameter<std::string>("PatternMetTrg", "");
}

void JetMETDQMPostProcessor::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
  //////////////////////////////////
  // setup DQM stor               //
  //////////////////////////////////

  bool isJetDir = false;
  bool isMetDir = false;

  TPRegexp patternJet(patternJetTrg_);
  TPRegexp patternMet(patternMetTrg_);

  // go to the directory to be processed
  if (igetter.dirExists(subDir_))
    ibooker.cd(subDir_);
  else {
    edm::LogWarning("JetMETDQMPostProcessor") << "cannot find directory: " << subDir_ << " , skipping";
    return;
  }

  std::vector<std::string> subdirectories = igetter.getSubdirs();
  for (std::vector<std::string>::iterator dir = subdirectories.begin(); dir != subdirectories.end(); dir++) {
    ibooker.cd(*dir);

    isJetDir = false;
    isMetDir = false;

    if (TString(*dir).Contains(patternJet))
      isJetDir = true;
    if (TString(*dir).Contains(patternMet))
      isMetDir = true;

    if (isMetDir) {
      // GenMET
      dividehistos(
          ibooker, igetter, "GenMETTrgMC", "GenMET", "TurnOngMET", "Gen Missing ET", "Gen Missing ET Turn-On RelVal");
      dividehistos(ibooker,
                   igetter,
                   "GenMETTrg",
                   "GenMETTrgLow",
                   "TurnOngMETLow",
                   "Gen Missing ETLow",
                   "Gen Missing ET Turn-On Data");

      // HLTMET
      dividehistos(
          ibooker, igetter, "HLTMETTrgMC", "HLTMET", "TurnOnhMET", "HLT Missing ET", "HLT Missing ET Turn-On RelVal");
      dividehistos(ibooker,
                   igetter,
                   "HLTMETTrg",
                   "HLTMETTrgLow",
                   "TurnOnhMETLow",
                   "HLT Missing ETLow",
                   "HLT Missing ET Turn-On Data");
    }

    if (isJetDir) {
      // std::cout << "JetMETDQMPostProcessor - Jet paths: " << ibooker.pwd() <<
      // " " << *dir << std::endl;

      // GenJets
      dividehistos(
          ibooker, igetter, "GenJetPtTrgMC", "GenJetPt", "TurnOngJetPt", "Gen Jet Pt", "Gen Jet Pt Turn-On RelVal");
      dividehistos(ibooker,
                   igetter,
                   "GenJetPtTrg",
                   "GenJetPtTrgLow",
                   "TurnOngJetPtLow",
                   "Gen Jet PtLow",
                   "Gen Jet Pt Turn-On Data");
      dividehistos(
          ibooker, igetter, "GenJetEtaTrgMC", "GenJetEta", "TurnOngJetEta", "Gen Jet Eta", "Gen Jet Eta Turn-On RelVal");
      dividehistos(ibooker,
                   igetter,
                   "GenJetEtaTrg",
                   "GenJetEtaTrgLow",
                   "TurnOngJetEtaLow",
                   "Gen Jet EtaLow",
                   "Gen Jet Eta Turn-On Data");
      dividehistos(
          ibooker, igetter, "GenJetPhiTrgMC", "GenJetPhi", "TurnOngJetPhi", "Gen Jet Phi", "Gen Jet Phi Turn-On RelVal");
      dividehistos(ibooker,
                   igetter,
                   "GenJetPhiTrg",
                   "GenJetPhiTrgLow",
                   "TurnOngJetPhiLow",
                   "Gen Jet PhiLow",
                   "Gen Jet Phi Turn-On Data");

      // HLTJets
      dividehistos(
          ibooker, igetter, "HLTJetPtTrgMC", "HLTJetPt", "TurnOnhJetPt", "HLT Jet Pt", "HLT Jet Pt Turn-On RelVal");
      dividehistos(ibooker,
                   igetter,
                   "HLTJetPtTrg",
                   "HLTJetPtTrgLow",
                   "TurnOnhJetPtLow",
                   "HLT Jet PtLow",
                   "HLT Jet Pt Turn-On Data");
      dividehistos(
          ibooker, igetter, "HLTJetEtaTrgMC", "HLTJetEta", "TurnOnhJetEta", "HLT Jet Eta", "HLT Jet Eta Turn-On RelVal");
      dividehistos(ibooker,
                   igetter,
                   "HLTJetEtaTrg",
                   "HLTJetEtaTrgLow",
                   "TurnOnhJetEtaLow",
                   "HLT Jet EtaLow",
                   "HLT Jet Eta Turn-On Data");
      dividehistos(
          ibooker, igetter, "HLTJetPhiTrgMC", "HLTJetPhi", "TurnOnhJetPhi", "HLT Jet Phi", "HLT Jet Phi Turn-On RelVal");
      dividehistos(ibooker,
                   igetter,
                   "HLTJetPhiTrg",
                   "HLTJetPhiTrgLow",
                   "TurnOnhJetPhiLow",
                   "HLT Jet PhiLow",
                   "HLT Jet Phi Turn-On Data");
    }

    ibooker.goUp();
  }
}

//----------------------------------------------------------------------
TProfile *JetMETDQMPostProcessor::dividehistos(DQMStore::IBooker &ibooker,
                                               DQMStore::IGetter &igetter,
                                               const std::string &numName,
                                               const std::string &denomName,
                                               const std::string &outName,
                                               const std::string &label,
                                               const std::string &titel) {
  TH1F *num = getHistogram(ibooker, igetter, ibooker.pwd() + "/" + numName);
  TH1F *denom = getHistogram(ibooker, igetter, ibooker.pwd() + "/" + denomName);

  if (num == nullptr)
    edm::LogWarning("JetMETDQMPostProcessor")
        << "numerator histogram " << ibooker.pwd() + "/" + numName << " does not exist";
  if (denom == nullptr)
    edm::LogWarning("JetMETDQMPostProcessor")
        << "denominator histogram " << ibooker.pwd() + "/" + denomName << " does not exist";

  // Check if histograms actually exist
  if (!num || !denom)
    return nullptr;

  MonitorElement *meOut = ibooker.bookProfile(
      outName, titel, num->GetXaxis()->GetNbins(), num->GetXaxis()->GetXmin(), num->GetXaxis()->GetXmax(), 0., 1.2);
  meOut->setEfficiencyFlag();
  TProfile *out = meOut->getTProfile();
  out->GetXaxis()->SetTitle(label.c_str());
  out->SetYTitle("Efficiency");
  out->SetOption("PE");
  out->SetLineColor(2);
  out->SetLineWidth(2);
  out->SetMarkerStyle(20);
  out->SetMarkerSize(0.8);
  out->SetStats(kFALSE);

  for (int i = 1; i <= num->GetNbinsX(); i++) {
    double e, low, high;
    Efficiency((int)num->GetBinContent(i), (int)denom->GetBinContent(i), 0.683, e, low, high);
    double err = e - low > high - e ? e - low : high - e;
    // here is the trick to store info in TProfile:
    out->SetBinContent(i, e);
    out->SetBinEntries(i, 1);
    out->SetBinError(i, sqrt(e * e + err * err));
  }
  return out;
}

//----------------------------------------------------------------------
TH1F *JetMETDQMPostProcessor::getHistogram(DQMStore::IBooker &ibooker,
                                           DQMStore::IGetter &igetter,
                                           const std::string &histoPath) {
  ibooker.pwd();
  MonitorElement *monElement = igetter.get(histoPath);
  if (monElement != nullptr)
    return monElement->getTH1F();
  else
    return nullptr;
}

//----------------------------------------------------------------------
void JetMETDQMPostProcessor::Efficiency(
    int passing, int total, double level, double &mode, double &lowerBound, double &upperBound) {
  // protection
  if (total == 0) {
    mode = 0.5;
    lowerBound = 0;
    upperBound = 1;
    return;
  }
  mode = passing / ((double)total);

  // see http://root.cern.ch/root/html/TEfficiency.html#compare
  lowerBound = TEfficiency::Wilson(total, passing, level, false);
  upperBound = TEfficiency::Wilson(total, passing, level, true);
}

//----------------------------------------------------------------------

DEFINE_FWK_MODULE(JetMETDQMPostProcessor);