Add scanning part of the analysis

README.md

@@ -8,8 +8,8 @@ For recent instruction please visit: https://twiki.cern.ch/twiki/bin/view/CMS/Mi
   git cms-init
   git cms-addpkg RecoMET/METFilters 
 
-  git clone https://github.com/didukhle/MetScanning.git
-  scram b -j9
+  git clone https://github.com/cms-met/MetScanning
+  scram b -j $(nproc)
   ```
   You might need to run the following command if you want to access files via XROOT:
 ```

scan/README.md

@@ -0,0 +1,65 @@
+# METScanning
+
+This framework aims to centralize the MET scanners analysis in a more efficient framework.
+For more details see the [Wiki page](https://gitlab.cern.ch/mmahdavi/metscanning/-/wikis/home).
+
+
+# Mannual
+
+1. Set the environment and bBuild the executable file if needed.
+```
+cmsenv
+scram b -j $(nproc)
+```
+<br/>
+
+2. Run the analysis
+```
+runAnalysis INPUT_FILE1 [INPUT_FILE2 [...]] --output output.root --dataset DATASET --year YYYY --pass-rec [ON/off] --cross-check [on/OFF]
+```
+***`--cross-check ON` option will allow you to have cross-check histograms in NFR***
+
+- One can also run ```runAnalysis --help``` or ```runAnalysis -h``` for more detailed help message.
+- One may need to modify the configuration files stored in [`config`](https://github.com/cms-met/MetScanning/tree/master/scan/config) directory based on her/his analysis, such as different years, scenarios and etc. 
+- If one is interested in old-analysis strategy i.e. without BER and NFR control region, it is possible by passing the `--old-analysis` option to the program. Hence the obtained efficencies are with out any cut or ID applied.
+<br/>
+
+# Plotting Mannual
+For running plot scripts, you need to logout and login again if you have set `CMSSW` environment i.e. `cmsenv`. Then you need to set the proper environment,
+```
+cd plot
+. ./env
+cd -
+```
+There are two different python scripts which provide `efficiencies` and `occupancy maps` plots such as [plot_efficencies](https://github.com/cms-met/MetScanning/tree/master/scan/plot/plot_efficiencies) and [plot_occupancy_maps](https://github.com/cms-met/MetScanning/tree/master/scan/plot/plot_occupancy_maps) respectively. They are mentioned below separately.
+
+Before using the scripts, please make sure that the environment has been set up as mention above.
+
+- **Occupancy maps:**
+```
+plot_occupancy_maps [list of output root-files separated by space] --dataset [list of datasets w.r.t provided output root-files separated by space] --year YYYY
+```
+<br/>
+
+- **Efficiency plots:**
+
+For full plots, i.e all filters plots individually and merged plots for each provided dataset, `--all` option is needed like below.
+```
+plot_efficencies [list of output root-files separated by space] --dataset [list of datasets w.r.t provided output root-files separated by space] --year 2018 --all
+```
+
+If one needs also to obtain the actual final merged plot, all datasets must be provided. In this case the `--all` option is optional, for example, below command will plot all kind of plots and the actual final merged plot as well.
+```
+plot_efficencies jetht_output.root singlemuon_output.root egamma_output.root --dataset jetht singlemuon egamma --year 2018 --all
+```
+<br/>
+
+- **Plots Configurations:**
+
+The `plot_efficencies` script uses `yaml` configuration files which contain `binning` settings according to regions (`BER` and `NFR`) variables and/or datasets,  for each filter separately.
+
+The configuration files are provided in the [config](https://github.com/cms-met/MetScanning/tree/master/scan/config) directory. One may need to modify theses configuration files to obtian the proper plots.
+<br/>
+<br/>
+
+- *One can also run* ```plot_efficencies [--help, -h]``` or ```plot_occupancy_maps [--help, -h]``` for more detailed help messages.

scan/bin/BERegion.cc

@@ -0,0 +1,178 @@
+#include "BERegion.h"
+
+#include <TVector2.h>
+
+
+BERegion::BERegion(TTreeReader &reader, const std::string berCat)
+  : Tree(reader),
+    berCat_{berCat} {
+  cfg_ = YAML::LoadFile("config/BER.yaml");
+  minPfMet_ = Get<float>("min_pfmet");
+
+  if (berCat == "BadChCand") {
+    isBER = &BERegion::isBadChCandBER;
+    minDPhi_ = Get<float>("min_dphi");
+    pfCanMinPt_ = Get<float>("min_pfcandid_pt");
+    pfCanPdgId_ = Get<int>("abs_pfcand_pdgid");
+  }
+
+  else if (berCat == "GlbSTightHalo") {
+    isBER = &BERegion::isGlbSTHaloBER;
+    maxPfMetPhi_ = Get<float>("max_abs_pfmet_phi");
+    minPfMetPhi_ = Get<float>("min_abs_pfmet_phi");
+    minJetPt_ = Get<float>("min_jet_pt");
+    maxJetEta_ = Get<float>("max_abs_jet_eta");
+    maxJetChef_ = Get<float>("max_jet_chef");
+    jetPtGt25Num_ = Get<int>("num_jets_ptgt25");
+    jetPtGt200Num_ = Get<int>("num_jets_ptgt200");
+  }
+
+  else if (berCat == "BadPFMuon") {
+    isBER = &BERegion::isBadPFMuonBER;
+    minDPhi_ = Get<float>("min_dphi");
+    minMuPt_ = Get<float>("min_muon_pt");
+    maxAddMuPt_ = Get<float>("max_additional_muon_pt");
+    minDeltaR_ = Get<float>("min_delta_r");
+    jetPtGt25Num_ = Get<int>("num_jets_ptgt25");
+    jetPtGt200Num_ = Get<int>("num_jets_ptgt200");
+    minMuNum_ = Get<int>("min_num_muon");
+  }
+
+  else if (berCat == "EcalBadCal_U") {
+    isBER = &BERegion::isEcalBadCalBER;
+    minDPhi_ = Get<float>("min_dphi");
+    minJetPt_ = Get<float>("min_jet_pt");
+    maxJetEta_ = Get<float>("max_abs_jet_eta");
+    minJetEta_ = Get<float>("min_abs_jet_eta");
+    minJetNeef_ = Get<float>("min_jet_neef");
+    jetPtGt25Num_ = Get<int>("num_jets_ptgt25");
+    jetPtGt200Num_ = Get<int>("num_jets_ptgt200");
+  }
+
+  else if (berCat == "HBHENoise") {
+    isBER = &BERegion::isHBHENoiseBER;
+    minDPhi_ = Get<float>("min_dphi");
+    minJetPt_ = Get<float>("min_jet_pt");
+    maxJetEta_ = Get<float>("max_abs_jet_eta");
+    minJetNhef_ = Get<float>("min_jet_nhef");
+    jetPtGt25Num_ = Get<int>("num_jets_ptgt25");
+    jetPtGt200Num_ = Get<int>("num_jets_ptgt200");
+  }
+
+  else
+    isBER = &BERegion::isOthersBER;
+
+}
+
+
+bool BERegion::operator()() {
+  return (this->*isBER)();
+}
+
+
+bool BERegion::isBadChCandBER() {
+  int i = 0;
+  if (*pfmet > minPfMet_) {
+    for (long unsigned int j = 0; j < (*pfCandPt).size(); j++)
+      if ((*pfCandPt)[j] > pfCanMinPt_ and
+          std::abs((*pfCandPdgId)[j]) == pfCanPdgId_ and
+          std::fabs(
+            TVector2::Phi_mpi_pi(*pfmetPhi-(*pfCandPhi)[j])) > minDPhi_)
+        i++;
+  }
+
+  return (i >= 1);
+}
+
+
+bool BERegion::isGlbSTHaloBER() {
+  if (*pfmet > minPfMet_ and
+      (*jetPt).size() == 1 and
+      (*numJetsPt25) == 1 and
+      (std::fabs(TVector2::Phi_mpi_pi(*pfmetPhi)) < maxPfMetPhi_ or
+       std::fabs(TVector2::Phi_mpi_pi(*pfmetPhi)) > minPfMetPhi_) and
+      (*jetPt)[0] > minJetPt_ and
+      std::fabs((*jetEta)[0]) < maxJetEta_ and
+      (*jetCHEF)[0] < maxJetChef_)
+    return true;
+
+  return false;
+}
+
+
+bool BERegion::isBadPFMuonBER() {
+  if (*pfmet > minPfMet_ and
+      (*jetPt).size() == jetPtGt200Num_ and
+      (*numJetsPt25) == jetPtGt25Num_ and
+      (*muonPt).size() >= minMuNum_ and
+      (*muonPt)[0] > minMuPt_ and
+      std::fabs(
+        TVector2::Phi_mpi_pi(*pfmetPhi - (*muonPhi)[0])) >= minDPhi_ and
+      (((*muonPt).size() > minMuNum_) ? ((*muonPt)[1] < maxAddMuPt_) : true)
+      ) {
+    int k = 0;
+    for (long unsigned int j = 0; j < (*jetPt).size(); j++) {
+      // To be sure that the jets are not comming from muon jet cone
+      // by requiring DR < 0.4
+      float const deltaR = std::sqrt(
+          std::pow((*muonEta)[0] - (*jetEta)[j], 2) +
+          std::pow((*muonPhi)[0] - (*jetPhi)[j], 2));
+      if (deltaR < minDeltaR_)
+        continue;
+      k++;
+      break;
+    }
+
+    return (k == 0);
+  }
+
+  return false;
+}
+
+
+bool BERegion::isEcalBadCalBER() {
+  if (*pfmet > minPfMet_ and
+      (*jetPt).size() == jetPtGt200Num_ and
+      (*numJetsPt25) == jetPtGt25Num_ and
+      (*jetPt)[0] > minJetPt_ and
+      std::fabs((*jetEta)[0]) > minJetEta_ and
+      std::fabs((*jetEta)[0]) < maxJetEta_ and
+      (*jetNEEF)[0] > minJetNeef_ and
+      std::fabs(TVector2::Phi_mpi_pi(*pfmetPhi - (*jetPhi)[0])) > minDPhi_)
+    return true;
+
+  return false;
+}
+
+
+bool BERegion::isHBHENoiseBER() {
+  if (*pfmet > minPfMet_ and
+      (*jetPt).size() == jetPtGt200Num_ and
+      (*numJetsPt25) == jetPtGt25Num_ and
+      (*jetPt)[0] > minJetPt_ and
+      std::fabs((*jetEta)[0]) < maxJetEta_ and
+      (*jetNHEF)[0] > minJetNhef_ and
+      std::fabs(TVector2::Phi_mpi_pi(*pfmetPhi - (*jetPhi)[0])) > minDPhi_)
+    return true;
+
+  return false;
+}
+
+
+bool BERegion::isOthersBER() {
+  return (*pfmet > minPfMet_);
+}
+
+
+template <typename T>
+T BERegion::Get(std::string const field) {
+  if (cfg_[berCat_][field])
+    return cfg_[berCat_][field].as<T>();
+  else {
+    std::ostringstream message;
+    message << "BER configuration of " << berCat_
+      << " doesn't contain [" << field << "] field.";
+    throw std::runtime_error(message.str());
+  }
+}
+

scan/bin/BERegion.h

@@ -0,0 +1,56 @@
+#ifndef BEREGION_H
+#define BEREGION_H
+
+#include <yaml-cpp/yaml.h>
+
+#include "Tree.h"
+
+
+class BERegion final : private Tree {
+ public:
+   // Constructor
+   BERegion(TTreeReader &reader, const std::string="Others");
+   bool operator()();
+
+ private:
+   std::string berCat_;
+
+   // BER configuration.
+   YAML::Node cfg_;
+
+   // All BER parameters
+   float minPfMet_, minDPhi_, maxPfMetPhi_, minPfMetPhi_, minJetPt_, maxJetEta_,
+         minJetEta_, maxJetChef_, minJetNeef_, minJetNhef_, minMuPt_,
+         maxAddMuPt_, minDeltaR_, pfCanMinPt_;
+
+   int pfCanPdgId_, jetPtGt25Num_;
+
+   long unsigned int jetPtGt200Num_, minMuNum_; 
+
+   // Check whether the event is in BER of BadChCand filters.
+   bool isBadChCandBER();
+
+   // Check whether the event is in BER of GlbSTHalo and GlbTHalo filters.
+   bool isGlbSTHaloBER();
+
+   // Check whether the event is in BER of BadPFMuon filters.
+   bool isBadPFMuonBER();
+
+   // Check whether the event is in BER of EcalBadCal filters.
+   bool isEcalBadCalBER();
+
+   // Check whether the event is in BER of HBHENoise and HBHEIsoNoise filters.
+   bool isHBHENoiseBER();
+
+   // Check whether the event is in BER of other filters.
+   bool isOthersBER();
+
+   // A pointer to one of functions checking whether the event is in BER or not.
+   bool (BERegion::*isBER)();
+
+   // Return the value of given field from configuration file.
+   template <typename T>
+   T Get(std::string const field);
+};
+
+#endif

scan/bin/BuildFile.xml

@@ -0,0 +1,10 @@
+<flags CXXFLAGS="-std=c++17 -Wall -Wextra -O3" />
+
+<bin name="runAnalysis" file="*.cc">
+  <use name="boost_filesystem" />
+  <use name="boost_program_options" />
+  <use name="root" />
+  <use name="rootgraphics" />
+  <use name="yaml-cpp" />
+</bin>
+