update for RunIISummer20UL top tagging scale factors for AK8 (with and without subjet b tagging) and HOTVR

RunIISummer20UL/README.md

@@ -0,0 +1,7 @@
+# TopTaggingScaleFactors (Run II Ultra Legacy)
+
+Please refer to this TWiki for documentation on the scale factors: https://twiki.cern.ch/twiki/bin/save/CMS/JetTopTagging
+
+You can have a look at the provided ROOT macro `readScaleFactor.cxx` to learn how to extract the scale factors from the ROOT files. Within the files, they are given in the form of `TGraphAsymmErrors`. Execute the macro with `root -l -b -q readScaleFactor.cxx`.
+
+Additionally, you can have a look at `identifyMergeCategory_PseudoCode.cxx` to learn how to get the correct merge category for each tagged jet. Please be aware that this is only pseudocode and not executable.

RunIISummer20UL/identifyMergeCategory_PseudoCode.cxx

@@ -0,0 +1,59 @@
+// Pseudocode! This does not run!
+
+enum class EMergeCategory {
+  isFullyMerged,
+  isSemiMerged,
+  isNotMerged,
+};
+
+EMergeCategory identifyMergeCategory_AK8(const vector<GenTop> & tops, const Jet & tagged_jet) {
+  if( tops.size() == 0 ) {
+    throw std::domain_error("The provided scale factors are not supposed to be applied on MC not containing top quarks");
+  }
+  EMergeCategory result = EMergeCategory::isNotMerged;
+  for(const auto top : tops) {
+    if( ! top->IsHadronicDecay() ) continue; // only generated top quarks decaying hadronically are considered for determining the merge category
+    Int_t n_matched = 0;
+    if( deltaR( top.bQuark(), tagged_jet ) < 0.8) n_matched++;
+    if( deltaR( top.WDaughter1(), tagged_jet ) < 0.8) n_matched++;
+    if( deltaR( top.WDaughter2(), tagged_jet ) < 0.8) n_matched++;
+    if( n_matched == 3 ) {
+      result = EMergeCategory::isFullyMerged;
+      break;
+    }
+    else if( n_matched == 2 ) {
+      result = EMergeCategory::isSemiMerged;
+      break;
+    }
+  }
+  return result;
+}
+
+EMergeCategory identifyMergeCategory_HOTVR(const vector<GenTop> & tops, const Jet & tagged_jet) {
+  if( tops.size() == 0 ) {
+    throw std::domain_error("The provided scale factors are not supposed to be applied on MC not containing top quarks");
+  }
+  EMergeCategory result = EMergeCategory::isNotMerged;
+  // HOTVR has a variable jet distance parameter; it is calculated using the raw (i.e. uncorrected) jet pT
+  // It is important to use the raw pT, since there are no energy corrections applied to the constituents during the clustering
+  const Double_t radius_min = 0.1;
+  const Double_t radius_max = 1.5;
+  const Double_t hotvr_rho = 600.;
+  const Double_t radius = min( radius_max, max( radius_min, hotvr_rho / ( tagged_jet.pt() * tagged_jet.JEC_factor_raw() ) ) );
+  for(const auto top : tops) {
+    if( ! top->IsHadronicDecay() ) continue; // only generated top quarks decaying hadronically are considered for determining the merge category
+    Int_t n_matched = 0;
+    if( deltaR( top.bQuark(), tagged_jet ) < radius) n_matched++;
+    if( deltaR( top.WDaughter1(), tagged_jet ) < radius) n_matched++;
+    if( deltaR( top.WDaughter2(), tagged_jet ) < radius) n_matched++;
+    if( n_matched == 3 ) {
+      result = EMergeCategory::isFullyMerged;
+      break;
+    }
+    else if( n_matched == 2 ) {
+      result = EMergeCategory::isSemiMerged;
+      break;
+    }
+  }
+  return result;
+}

RunIISummer20UL/readScaleFactor.cxx

@@ -0,0 +1,187 @@
+// includes only needed for nicely printing scale factors to terminal:
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+
+
+//______________________________
+// CONSTANTS AND TYPEDEFS
+
+enum class EYear {
+  isUL16preVFP,
+  isUL16postVFP,
+  isUL17,
+  isUL18,
+};
+
+const std::map<EYear, TString> kRootFiles = {
+  {EYear::isUL16preVFP, "scaleFactors/TopTaggingScaleFactors_RunIISummer20UL16preVFP_PUPPIv15.root"},
+  {EYear::isUL16postVFP, "scaleFactors/TopTaggingScaleFactors_RunIISummer20UL16postVFP_PUPPIv15.root"},
+  {EYear::isUL17, "scaleFactors/TopTaggingScaleFactors_RunIISummer20UL17_PUPPIv15.root"},
+  {EYear::isUL18, "scaleFactors/TopTaggingScaleFactors_RunIISummer20UL18_PUPPIv15.root"},
+};
+
+enum class ETagger {
+  isAK8,
+  isAK8_subjetbtag,
+  isHOTVR,
+};
+
+enum class EWorkingPoint {
+  isWP0p38, // for AK8 taggers
+  isWP0p47, // %
+  isWP0p52, // %
+  isWP0p61, // %
+  isWP0p69, // %
+  isDefault, // for HOTVR only
+};
+
+enum class EMergeCategory {
+  isFullyMerged,
+  isSemiMerged,
+  isNotMerged,
+};
+
+const std::map<ETagger, TString> kTaggerToString = {
+  {ETagger::isAK8, "AK8_PUPPI"},
+  {ETagger::isAK8_subjetbtag, "AK8_PUPPI_subjetbtag"},
+  {ETagger::isHOTVR, "HOTVR_PUPPI"},
+};
+
+const std::map<EWorkingPoint, TString> kWorkingPointToString = {
+  {EWorkingPoint::isWP0p38, "wp0p38_vt"},
+  {EWorkingPoint::isWP0p47, "wp0p47_t"},
+  {EWorkingPoint::isWP0p52, "wp0p52_m"},
+  {EWorkingPoint::isWP0p61, "wp0p61_l"},
+  {EWorkingPoint::isWP0p69, "wp0p69_vl"},
+  {EWorkingPoint::isDefault, "default"},
+};
+
+const std::map<EWorkingPoint, TString> kWorkingPointToMisTagEffString = {
+  {EWorkingPoint::isWP0p38, "mis0p001"},
+  {EWorkingPoint::isWP0p47, "mis0p005"},
+  {EWorkingPoint::isWP0p52, "mis0p010"},
+  {EWorkingPoint::isWP0p61, "mis0p025"},
+  {EWorkingPoint::isWP0p69, "mis0p050"},
+};
+
+const std::map<EMergeCategory, TString> kMergeCategoryToString = {
+  {EMergeCategory::isFullyMerged, "FullyMerged"},
+  {EMergeCategory::isSemiMerged, "SemiMerged"},
+  {EMergeCategory::isNotMerged, "NotMerged"},
+};
+
+typedef struct {
+  Double_t fCentral = 1; // central value
+  Double_t fErrTotUp = 0; // total uncertainty
+  Double_t fErrTotDown = 0; // %
+  Double_t fErrStatUp = 0; // statistical component
+  Double_t fErrStatDown = 0; // %
+  Double_t fErrSystUp = 0; // systematic component
+  Double_t fErrSystDown = 0; // %
+  Double_t fErrUncorrUp = 0; // era-uncorrelated component
+  Double_t fErrUncorrDown = 0; // %
+  Double_t fErrCorrUp = 0; // era-correlated component
+  Double_t fErrCorrDown = 0; // %
+} ScaleFactor_t;
+
+
+//______________________________
+// FUNCTIONS
+
+ScaleFactor_t loadScaleFactor(const Double_t jetPt, const ETagger tagger, const EWorkingPoint wp, const EMergeCategory mergeCat, const EYear year) {
+  if((wp == EWorkingPoint::isDefault && tagger != ETagger::isHOTVR) || (tagger == ETagger::isHOTVR && wp != EWorkingPoint::isDefault)) {
+    throw std::invalid_argument("Given combination of tagger/WP not valid");
+  }
+  TFile *rootFile = TFile::Open(kRootFiles.at(year).Data(), "READ");
+  const TString graph_basename = kTaggerToString.at(tagger)+"_"+kWorkingPointToString.at(wp)+(tagger == ETagger::isAK8 ? "_"+kWorkingPointToMisTagEffString.at(wp) : "")+"/"+kMergeCategoryToString.at(mergeCat)+"_";
+  const TGraphAsymmErrors *graph_tot = (TGraphAsymmErrors*)rootFile->Get((graph_basename+"tot").Data());
+  const TGraphAsymmErrors *graph_stat = (TGraphAsymmErrors*)rootFile->Get((graph_basename+"stat").Data());
+  const TGraphAsymmErrors *graph_syst = (TGraphAsymmErrors*)rootFile->Get((graph_basename+"syst").Data());
+  const TGraphAsymmErrors *graph_uncorr = (TGraphAsymmErrors*)rootFile->Get((graph_basename+"uncorr").Data());
+  const TGraphAsymmErrors *graph_corr = (TGraphAsymmErrors*)rootFile->Get((graph_basename+"corr").Data());
+
+  UInt_t ibin = 0;
+  Double_t x;
+  Double_t y;
+  graph_tot->GetPoint(ibin, x, y);
+
+  ScaleFactor_t sf;
+
+  const Double_t minimum_jetPt = x - graph_tot->GetErrorXlow(ibin);
+  if(jetPt < minimum_jetPt) {
+    std::cerr << "No scale factor for such low jet pT available. Will set scale factor to 1 with zero uncertainty" << std::endl;
+    return sf;
+  }
+
+  while(!(jetPt >= x - graph_tot->GetErrorXlow(ibin) && jetPt < x + graph_tot->GetErrorXhigh(ibin))) { // check if jetPt within pT interval of bin with bin number "ibin"
+    ++ibin;
+    graph_tot->GetPoint(ibin, x, y);
+    if(ibin + 1 >= graph_tot->GetN()) break; // use scale factor of last bin if given jetPt exceeds highest bin edge of the graph
+  }
+
+  sf = ScaleFactor_t{
+    .fCentral = y,
+    .fErrTotUp = graph_tot->GetErrorYhigh(ibin),
+    .fErrTotDown = graph_tot->GetErrorYlow(ibin),
+    .fErrStatUp = graph_stat->GetErrorYhigh(ibin),
+    .fErrStatDown = graph_stat->GetErrorYlow(ibin),
+    .fErrSystUp = graph_syst->GetErrorYhigh(ibin),
+    .fErrSystDown = graph_syst->GetErrorYlow(ibin),
+    .fErrUncorrUp = graph_uncorr->GetErrorYhigh(ibin),
+    .fErrUncorrDown = graph_uncorr->GetErrorYlow(ibin),
+    .fErrCorrUp = graph_corr->GetErrorYhigh(ibin),
+    .fErrCorrDown = graph_corr->GetErrorYlow(ibin),
+  };
+
+  rootFile->Close();
+
+  return sf;
+}
+
+void printScaleFactor(const ScaleFactor_t & sf, const Int_t precision = 3) {
+  std::stringstream ss;
+  ss << "SF: ";
+  ss << std::fixed << std::setprecision(precision) << sf.fCentral;
+  ss << " +" << std::fixed << std::setprecision(precision) << sf.fErrTotUp;
+  ss << "/";
+  ss << "-" << std::fixed << std::setprecision(precision) << sf.fErrTotDown;
+  ss << " (tot)";
+  //________________
+  ss << " [";
+  ss << " +" << std::fixed << std::setprecision(precision) << sf.fErrStatUp;
+  ss << "/";
+  ss << "-" << std::fixed << std::setprecision(precision) << sf.fErrStatDown;
+  ss << " (stat) |";
+  ss << " +" << std::fixed << std::setprecision(precision) << sf.fErrSystUp;
+  ss << "/";
+  ss << "-" << std::fixed << std::setprecision(precision) << sf.fErrSystDown;
+  ss << " (syst) ]";
+  //________________
+  ss << " [";
+  ss << " +" << std::fixed << std::setprecision(precision) << sf.fErrUncorrUp;
+  ss << "/";
+  ss << "-" << std::fixed << std::setprecision(precision) << sf.fErrUncorrDown;
+  ss << " (uncorr) |";
+  ss << " +" << std::fixed << std::setprecision(precision) << sf.fErrCorrUp;
+  ss << "/";
+  ss << "-" << std::fixed << std::setprecision(precision) << sf.fErrCorrDown;
+  ss << " (corr) ]";
+  std::cout << ss.str() << std::endl;
+}
+
+
+//______________________________
+// MACRO
+
+void readScaleFactor() {
+  // Exemplary configuration:
+  const Double_t jetPt = 345.6789; // in GeV
+  const ETagger tagger = ETagger::isAK8;
+  const EWorkingPoint wp = EWorkingPoint::isWP0p52;
+  const EMergeCategory mergeCat = EMergeCategory::isFullyMerged; // in practice, this needs to be determined for each individual jet; see TWiki how to do it
+  const EYear year = EYear::isUL17;
+
+  const ScaleFactor_t sf = loadScaleFactor(jetPt, tagger, wp, mergeCat, year);
+  printScaleFactor(sf);
+}