Factories based on python decorations of the trees

.gitignore

@@ -1,2 +1 @@
-python/__init__.py
 *.pyc

common/include/IndexRangeIterator.h

@@ -0,0 +1,28 @@
+#pragma once
+
+#include <boost/iterator/iterator_facade.hpp>
+
+/*
+ * Provide an iterator interface over an index that runs from 0 to N
+ */
+template<typename IDX>
+class IndexRangeIterator
+  : public boost::iterator_facade<IndexRangeIterator<IDX>,
+      IDX, // value
+      boost::random_access_traversal_tag,
+      IDX  // ref
+    >
+{
+public:
+  IndexRangeIterator(IDX idx, IDX max)
+    : m_idx{idx}, m_max{max} {}
+
+  IDX dereference() const { return m_idx; }
+  bool equal( const IndexRangeIterator& other ) const { return ( m_max == other.m_max ) && ( m_idx == other.m_idx ); }
+  void increment() { ++m_idx; }
+  void decrement() { --m_idx; }
+  void advance(std::ptrdiff_t d) { m_idx += d; }
+  std::ptrdiff_t distance_to(const IndexRangeIterator& other) const { return m_idx-other.m_idx; }
+private:
+  IDX m_idx, m_max;
+};

common/include/ScaleFactors.h

@@ -0,0 +1,203 @@
+#pragma once
+
+#include "BinnedValuesJSONParser.h"
+#include <stdexcept>
+
+class ILeptonScaleFactor {
+public:
+  virtual ~ILeptonScaleFactor() {}
+  virtual float get(const Parameters& parameters, Variation variation) const = 0;
+};
+class IDiLeptonScaleFactor {
+public:
+  virtual ~IDiLeptonScaleFactor() {}
+  virtual float get(const Parameters& l1Param, const Parameters& l2Param, Variation variation) const = 0;
+};
+
+class IJetScaleFactor {
+public:
+  virtual ~IJetScaleFactor() {}
+
+  enum Flavour {
+    Light  = 0,
+    Charm  = 1,
+    Beauty = 2
+  };
+  static inline Flavour get_flavour(int hadron_flavour) {
+    switch(hadron_flavour) {
+      case 5:
+        return Flavour::Beauty;
+      case 4:
+        return Flavour::Charm;
+      default:
+        return Flavour::Light;
+    }
+  }
+  virtual float get(const Parameters& parameters, Flavour flavour, Variation variation) const = 0;
+};
+
+/**
+ * Simple case: one scale factor from file
+ */
+class ScaleFactor : public ILeptonScaleFactor {
+public:
+  explicit ScaleFactor(const std::string& file)
+  : m_values{BinnedValuesJSONParser(file).get_values()}
+  {}
+  virtual ~ScaleFactor() {}
+
+  virtual float get(const Parameters& parameters, Variation variation) const override final {
+    return m_values.get(parameters)[static_cast<std::size_t>(variation)];
+  }
+private:
+  BinnedValues m_values;
+};
+
+/**
+ * Dilepton scalefactor from two lepton scalefactors (takes a parameter set for each)
+ */
+class DiLeptonFromLegsScaleFactor : public IDiLeptonScaleFactor {
+public:
+  DiLeptonFromLegsScaleFactor( std::unique_ptr<ILeptonScaleFactor>&& l1_leg1
+                             , std::unique_ptr<ILeptonScaleFactor>&& l1_leg2
+                             , std::unique_ptr<ILeptonScaleFactor>&& l2_leg1
+                             , std::unique_ptr<ILeptonScaleFactor>&& l2_leg2 )
+  : m_l1_leg1{std::move(l1_leg1)}
+  , m_l1_leg2{std::move(l1_leg2)}
+  , m_l2_leg1{std::move(l2_leg1)}
+  , m_l2_leg2{std::move(l2_leg2)}
+  {}
+  virtual ~DiLeptonFromLegsScaleFactor() {}
+
+  virtual float get(const Parameters& l1Param, const Parameters& l2Param, Variation variation) const override final {
+    const float eff_lep1_leg1 = m_l1_leg1->get(l1Param, variation);
+    const float eff_lep1_leg2 = m_l1_leg2->get(l1Param, variation);
+    const float eff_lep2_leg1 = m_l2_leg1->get(l2Param, variation);
+    const float eff_lep2_leg2 = m_l2_leg2->get(l2Param, variation);
+
+    if ( variation == Nominal ) {
+      return -(eff_lep1_leg1 * eff_lep2_leg1) +
+          (1 - (1 - eff_lep1_leg2)) * eff_lep2_leg1 +
+          eff_lep1_leg1 * (1 - (1 - eff_lep2_leg2)) ;
+    } else {
+      const float nom_eff_lep1_leg1 = m_l1_leg1->get(l1Param, Nominal);
+      const float nom_eff_lep1_leg2 = m_l1_leg2->get(l1Param, Nominal);
+      const float nom_eff_lep2_leg1 = m_l2_leg1->get(l2Param, Nominal);
+      const float nom_eff_lep2_leg2 = m_l2_leg2->get(l2Param, Nominal);
+
+      const float nominal = -(eff_lep1_leg1 * eff_lep2_leg1) +
+          (1 - (1 - eff_lep1_leg2)) * eff_lep2_leg1 +
+          eff_lep1_leg1 * (1 - (1 - eff_lep2_leg2)) ;
+
+      const float error_squared =
+          std::pow(1 - nom_eff_lep2_leg1 - (1 - nom_eff_lep2_leg2), 2) *
+          std::pow(eff_lep1_leg1, 2) +
+          std::pow(nom_eff_lep2_leg1, 2) *
+          std::pow(eff_lep1_leg2, 2) +
+          std::pow(1 - nom_eff_lep1_leg1 - (1 - nom_eff_lep1_leg2), 2) *
+          std::pow(eff_lep2_leg1, 2) +
+          std::pow(nom_eff_lep1_leg1, 2) *
+          std::pow(eff_lep2_leg2, 2);
+
+      if ( variation == Up ) {
+        return std::min(nominal+std::sqrt(error_squared), 1.f);
+      } else if ( variation == Down ) {
+        return std::max(0.f, nominal-std::sqrt(error_squared));
+      }
+    }
+  }
+private:
+  std::unique_ptr<ILeptonScaleFactor> m_l1_leg1;
+  std::unique_ptr<ILeptonScaleFactor> m_l1_leg2;
+  std::unique_ptr<ILeptonScaleFactor> m_l2_leg1;
+  std::unique_ptr<ILeptonScaleFactor> m_l2_leg2;
+};
+
+/**
+ * B-tagging scale factor (depends on flavour)
+ */
+class BTaggingScaleFactor : public IJetScaleFactor {
+public:
+  BTaggingScaleFactor(const std::string& lightFile, const std::string& charmFile, const std::string& beautyFile)
+  : m_lightValues {BinnedValuesJSONParser(lightFile ).get_values()}
+  , m_charmValues {BinnedValuesJSONParser(charmFile ).get_values()}
+  , m_beautyValues{BinnedValuesJSONParser(beautyFile).get_values()}
+  {}
+  virtual ~BTaggingScaleFactor() {}
+
+  virtual float get(const Parameters& parameters, Flavour flavour, Variation variation) const override final {
+    switch(flavour) {
+      case Flavour::Beauty:
+        return m_beautyValues.get(parameters)[static_cast<std::size_t>(variation)];
+      case Flavour::Charm:
+        return m_charmValues .get(parameters)[static_cast<std::size_t>(variation)];
+      case Flavour::Light:
+        return m_lightValues .get(parameters)[static_cast<std::size_t>(variation)];
+      default:
+        throw std::invalid_argument("Unsupported flavour: "+flavour);
+    }
+  }
+private:
+  BinnedValues m_lightValues;
+  BinnedValues m_charmValues;
+  BinnedValues m_beautyValues;
+};
+
+#include <vector>
+#include <algorithm>
+#include "IndexRangeIterator.h"
+/**
+ * Weight between different scale factors
+ */
+template<class ISingle,typename... Args>
+class WeightedScaleFactor : public ISingle {
+public:
+  WeightedScaleFactor( const std::vector<float>& probs, std::vector<std::unique_ptr<ISingle>>&& sfs )
+  : m_terms{std::move(sfs)}
+  {
+    const double norm{1./std::accumulate(std::begin(probs), std::end(probs), 0.)};
+    std::transform(std::begin(probs), std::end(probs), std::back_inserter(m_probs), [norm] ( float p ) -> float { return norm*p; } );
+  }
+  virtual ~WeightedScaleFactor() {}
+
+  virtual float get(Args&&... args, Variation variation) const override final {
+    return std::accumulate(IndexRangeIterator<std::size_t>(0, m_terms.size()), IndexRangeIterator<std::size_t>(m_terms.size(), m_terms.size()), 0.,
+        [this,&args...,variation] ( float wsum, std::size_t i ) -> float {
+          return wsum + m_probs[i]*m_terms[i]->get(std::forward<Args>(args)..., variation);
+        }
+      );
+  }
+private:
+  std::vector<float> m_probs;
+  std::vector<std::unique_ptr<ISingle>> m_terms;
+};
+using WScaleFactor         = WeightedScaleFactor<ILeptonScaleFactor,const Parameters&>;
+using WLLScaleFactor       = WeightedScaleFactor<IDiLeptonScaleFactor,const Parameters&,const Parameters&>;
+using WBTaggingScaleFactor = WeightedScaleFactor<IJetScaleFactor,const Parameters&,IJetScaleFactor::Flavour>;
+
+#include <random>
+/**
+ * Sample according to fractions from different scale factors
+ */
+template<class ISingle,typename... Args>
+class SampledScaleFactor : public ISingle
+{
+public:
+  SampledScaleFactor( const std::vector<float> probs, std::vector<std::unique_ptr<ISingle>>&& sfs )
+  : m_rg{42}
+  , m_probs{std::discrete_distribution<std::size_t>(std::begin(probs), std::end(probs))}
+  , m_terms{std::move(sfs)}
+  {}
+  virtual ~SampledScaleFactor() {}
+
+  virtual float get(Args... args, Variation variation) const override final {
+    return m_terms[m_probs(m_rg)]->get(std::forward<Args>(args)..., variation);
+  }
+private:
+  mutable std::mt19937 m_rg;
+  mutable std::discrete_distribution<std::size_t> m_probs;
+  std::vector<std::unique_ptr<ISingle>> m_terms;
+};
+using SmpScaleFactor         = SampledScaleFactor<ILeptonScaleFactor,const Parameters&>;
+using SmpLLScaleFactor       = SampledScaleFactor<IDiLeptonScaleFactor,const Parameters&,const Parameters&>;
+using SmpBTaggingScaleFactor = SampledScaleFactor<IJetScaleFactor,const Parameters&,IJetScaleFactor::Flavour>;

common/include/kinematics.h

@@ -0,0 +1,30 @@
+#pragma once
+
+#include <Math/LorentzVector.h>
+
+namespace Kinematics {
+  using LorentzVector = ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiE4D<float> >;
+
+  float deltaPhi( const LorentzVector& pa, const LorentzVector& pb )
+  {
+    return std::acos(std::min<float>(std::max<float>(-1.,
+            (pa.Px()*pb.Px()+pa.Py()*pb.Py())/(pa.Pt()*pb.Pt())
+           ), 1.));
+  }
+  float deltaR( const LorentzVector& pa, const LorentzVector& pb )
+  {
+    return std::sqrt(std::pow(pb.Eta()-pa.Eta(), 2)+std::pow(deltaPhi(pa, pb), 2));
+  }
+
+  float signedDeltaPhi( const LorentzVector& pa, const LorentzVector& pb )
+  {
+    return ( ( pb.Py()*pa.Px()-pb.Px()*pa.Py() > 0. ) ? 1. : -1. )*deltaPhi(pa, pb);
+  }
+
+  float signedDeltaEta( const LorentzVector& pa, const LorentzVector& pb )
+  {
+    const float etaA{pa.Eta()};
+    const float etaB{pb.Eta()};
+    return ( etaA > 0. ? 1. : -1. )*( etaB - etaA ); // (small) positive and negative = more and less forward
+  }
+};

installpy_standalone.sh

@@ -0,0 +1,42 @@
+#!/bin/sh
+#
+# Creates a scram-like (symlinked) python install directory for *Tools packages
+# in the equivalent of ${CMSSW_BASE}/src/install_python (or ${PREFIX}/install_python, if PREFIX is defined)
+#
+thisscript="$(realpath ${0})"
+commontoolspath="$(dirname ${thisscript})"
+cp3llbbpath="$(dirname ${commontoolspath})"
+if [[ "$(basename ${cp3llbbpath})" != "cp3_llbb" ]]; then
+  echo "Directory above CommonTools is not called cp3_llbb, aborting"
+  exit 1
+fi
+basepath="$(dirname ${cp3llbbpath})"
+if [[ -z "${PREFIX}" ]]; then
+  PREFIX="${basepath}"
+fi
+installpath="${PREFIX}/install_python"
+if [[ -a "${installpath}" ]]; then
+  echo "Install path ${installpath} already exists, aborting"
+  exit 1
+fi
+echo "Installing into ${installpath}, make sure to add it to PYTHONPATH as well"
+mkdir -p "${installpath}"
+pushd "${basepath}" > /dev/null
+for pkgpy in */*Tools/python; do
+  pkgpath_py="${basepath}/${pkgpy}"
+  if [[ ! -d "${pkgpath_py}" ]]; then
+    echo "ASSERT FAILURE: ${pkgpath_py}"
+    exit 1
+  fi
+  pkgname_full="$(dirname ${pkgpy})"
+  pkgname="$(basename ${pkgname_full})"
+  hatname="$(dirname  ${pkgname_full})"
+  hatinstalldir="${installpath}/${hatname}"
+  mkdir -p "${hatinstalldir}"
+  hatinitpy="${hatinstalldir}/__init__.py"
+  if [[ ! -f "${hatinitpy}" ]]; then
+    echo "" > "${hatinitpy}"
+  fi
+  ln -s "${pkgpath_py}" "${hatinstalldir}/${pkgname}"
+  echo "Installed ${pkgname_full}"
+done

python/__init__.py

@@ -0,0 +1,13 @@
+__all__ = ("pathCommonTools", "pathCP3llbb", "pathCMS")
+import os.path
+pathCommonTools = os.path.dirname(os.path.dirname(os.path.realpath(__file__)))
+pathCP3llbb = os.path.dirname(os.path.abspath(pathCommonTools))
+_pathCMS_paths = os.path.dirname(os.path.dirname(pathCP3llbb))
+import os
+_pathCMS_env = os.getenv("CMSSW_BASE")
+pathCMS = _pathCMS_paths
+if _pathCMS_env == "":
+    print("Warning: Could not get CMSSW_BASE variable from the environment")
+elif _pathCMS_env != _pathCMS_paths:
+    print("Warning: CommonTools is not using the CMSSW release version from $CMSSW_BASE ({0} versus {1}), using the former".format(_pathCMS_paths, _pathCMS_env))
+    pathCMS = _pathCMS_paths