diff --git a/.gitignore b/.gitignore
index b06920b275f..836aadc9444 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,8 +3,10 @@
*.so
*.pdf
*.zip
-fit.exe
*~
-*#*
-*.d
__init__.py
+output
+*.code-workspace
+/.python
+/.vscode
+.DS_Store
\ No newline at end of file
diff --git a/Common/BuildFile.xml b/Common/BuildFile.xml
index 798808346e1..8416ef8eeba 100644
--- a/Common/BuildFile.xml
+++ b/Common/BuildFile.xml
@@ -1,4 +1,7 @@
+
+
+
diff --git a/Common/interface/AnalysisTypes.h b/Common/interface/AnalysisTypes.h
index c1c577ffe57..194eea4a1b5 100644
--- a/Common/interface/AnalysisTypes.h
+++ b/Common/interface/AnalysisTypes.h
@@ -54,17 +54,23 @@ ENUM_NAMES(GenQcdMatch) = {
{ GenQcdMatch::Gluon, "gen_gluon" }
};
-enum class TauType { e = 0, mu = 1, tau = 2, jet = 3 };
-ENUM_NAMES(TauType) = {
- { TauType::e, "e" }, { TauType::mu, "mu" }, { TauType::tau, "tau" }, { TauType::jet, "jet" }
+enum class LegType { e = 1, mu = 2, tau = 4, jet = 8 };
+ENUM_NAMES(LegType) = {
+ { LegType::e, "e" }, { LegType::mu, "mu" }, { LegType::tau, "tau" }, { LegType::jet, "jet" }
};
-inline constexpr TauType GenMatchToTauType(GenLeptonMatch gen_match)
+inline constexpr LegType GenMatchToLegType(GenLeptonMatch gen_match)
{
- if(gen_match == GenLeptonMatch::Electron || gen_match == GenLeptonMatch::TauElectron) return TauType::e;
- if(gen_match == GenLeptonMatch::Muon || gen_match == GenLeptonMatch::TauMuon) return TauType::mu;
- if(gen_match == GenLeptonMatch::Tau) return TauType::tau;
- return TauType::jet;
+ if(gen_match == GenLeptonMatch::Electron || gen_match == GenLeptonMatch::TauElectron) return LegType::e;
+ if(gen_match == GenLeptonMatch::Muon || gen_match == GenLeptonMatch::TauMuon) return LegType::mu;
+ if(gen_match == GenLeptonMatch::Tau) return LegType::tau;
+ return LegType::jet;
}
+enum class TauSelection { gen = 1, pt = 2, MVA = 4, DeepTau = 8 };
+ENUM_NAMES(TauSelection) = {
+ { TauSelection::gen, "gen" }, { TauSelection::pt, "pt" }, { TauSelection::MVA, "MVA"},
+ { TauSelection::DeepTau, "DeepTau" }
+};
+
} // namespace analysis
diff --git a/Common/interface/CutTools.h b/Common/interface/CutTools.h
new file mode 100644
index 00000000000..53fc38b3767
--- /dev/null
+++ b/Common/interface/CutTools.h
@@ -0,0 +1,204 @@
+/*! Common tools and definitions to apply cuts.
+This file is part of https://github.com/cms-tau-pog/TauTriggerTools. */
+
+#pragma once
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+#include "SmartHistogram.h"
+
+namespace cuts {
+
+class cut_failed : public std::exception {
+public:
+ cut_failed(size_t parameter_id) noexcept
+ : _param_id(parameter_id)
+ {
+ std::ostringstream ss;
+ ss << "Cut requirements are not fulfilled for parameter id = " << _param_id << ".";
+ message = ss.str();
+ }
+
+ ~cut_failed() noexcept {}
+
+ virtual const char* what() const noexcept { return message.c_str(); }
+ size_t param_id() const noexcept { return _param_id; }
+
+private:
+ size_t _param_id;
+ std::string message;
+};
+
+template
+ValueType fill_histogram(ValueType value, Histogram& histogram, double weight)
+{
+ histogram.Fill(value,weight);
+ return value;
+}
+
+
+class ObjectSelector{
+public:
+
+ virtual ~ObjectSelector(){}
+
+ void incrementCounter(size_t param_id, const std::string& param_label)
+ {
+ if (counters.size() < param_id)
+ throw std::runtime_error("counters out of range");
+ if (counters.size() == param_id){ //counters and selections filled at least once
+ counters.push_back(0);
+ selections.push_back(0);
+ selectionsSquaredErros.push_back(0);
+ const std::string label = make_unique_label(param_label);
+ labels.push_back(label);
+ label_set.insert(label);
+ }
+ counters.at(param_id)++;
+ }
+
+ void fill_selection(double weight = 1.0){
+ for (unsigned n = 0; n < counters.size(); ++n){
+ if(counters.at(n) > 0) {
+ selections.at(n) += weight;
+ selectionsSquaredErros.at(n) += weight * weight;
+ }
+ counters.at(n) = 0;
+ }
+ }
+
+ template
+ std::vector collect_objects(double weight, size_t n_objects, const Selector& selector,
+ const Comparitor& comparitor)
+ {
+ std::vector selected;
+ for (size_t n = 0; n < n_objects; ++n) {
+ try {
+ const ObjectType selectedCandidate = selector(n);
+ selected.push_back(selectedCandidate);
+ } catch(cuts::cut_failed&) {}
+ }
+
+ fill_selection(weight);
+ std::sort(selected.begin(), selected.end(), comparitor);
+
+ return selected;
+ }
+
+private:
+ std::string make_unique_label(const std::string& label)
+ {
+ if(!label_set.count(label)) return label;
+ for(size_t n = 2; ; ++n) {
+ std::ostringstream ss;
+ ss << label << "_" << n;
+ if(!label_set.count(ss.str())) return ss.str();
+ }
+ }
+
+protected:
+ std::vector counters;
+ std::vector selections;
+ std::vector selectionsSquaredErros;
+ std::vector labels;
+ std::set label_set;
+};
+
+namespace detail {
+struct DefaultSelectionManager{
+ template
+ void FillHistogram(ValueType value, const std::string& histogram_name) {}
+};
+} // namespace detail
+
+template
+class Cutter {
+public:
+ explicit Cutter(ObjectSelector* _objectSelector, SelectionManager* _selectionManager = nullptr)
+ : objectSelector(_objectSelector), selectionManager(_selectionManager), param_id(0) {}
+
+ bool Enabled() const { return objectSelector != nullptr; }
+ int CurrentParamId() const { return param_id; }
+
+ void operator()(bool expected, const std::string& label)
+ {
+ (*this)(expected, label, expected);
+ }
+
+ template
+ void operator()(bool expected, const std::string& label, const ValueType& value)
+ {
+ if(selectionManager) {
+ try {
+ selectionManager->FillHistogram(value, label);
+ }catch(std::exception& e) {
+ std::cout << "ERROR: " << e.what() << std::endl;
+ }
+ }
+ if(Enabled()) {
+ ++param_id;
+ if(!expected)
+ throw cut_failed(param_id -1);
+ objectSelector->incrementCounter(param_id - 1, label);
+ }
+ }
+
+ bool test(bool expected, const std::string& label)
+ {
+ try {
+ (*this)(expected, label);
+ return true;
+ } catch(cut_failed&) {}
+ return false;
+ }
+
+private:
+ ObjectSelector* objectSelector;
+ SelectionManager* selectionManager;
+ size_t param_id;
+};
+
+} // cuts
+
+namespace root_ext {
+
+template<>
+class SmartHistogram : public cuts::ObjectSelector, public AbstractHistogram {
+public:
+ using RootContainer = TH1D;
+
+ SmartHistogram(const std::string& name) : AbstractHistogram(name) {}
+
+ void SetSave(bool _save)
+ {
+ save = _save;
+ }
+
+ virtual void WriteRootObject()
+ {
+ if(!save || !selections.size() || !GetOutputDirectory() )
+ return;
+ std::unique_ptr selection_histogram(
+ new TH1D(Name().c_str(), Name().c_str(),selections.size(),-0.5,-0.5+selections.size()));
+ for (unsigned n = 0; n < selections.size(); ++n){
+ const std::string label = labels.at(n);
+ selection_histogram->GetXaxis()->SetBinLabel(n+1, label.c_str());
+ selection_histogram->SetBinContent(n+1,selections.at(n));
+ selection_histogram->SetBinError(n+1,std::sqrt(selectionsSquaredErros.at(n)));
+ }
+ root_ext::WriteObject(*selection_histogram, GetOutputDirectory());
+ }
+
+private:
+ bool save{true};
+};
+}
diff --git a/Common/interface/GenTruthTools.h b/Common/interface/GenTruthTools.h
index 8b77672ed64..7974fad1be6 100644
--- a/Common/interface/GenTruthTools.h
+++ b/Common/interface/GenTruthTools.h
@@ -13,24 +13,42 @@ namespace analysis {
namespace gen_truth {
+struct FinalState {
+public:
+ enum class ParticleType { visible, light_lepton, neutrino, gamma, charged_hadron, neutral_hadron };
+
+ explicit FinalState(const reco::GenParticle& particle, const std::set& pdg_to_exclude = {},
+ const std::set& particles_to_exclude = {});
+
+ const std::set& getParticles(ParticleType type) { return particles[type]; }
+ const LorentzVectorXYZ& getMomentum(ParticleType type) { return momentum[type]; }
+ size_t count(ParticleType type) { return getParticles(type).size(); }
+
+private:
+ void findFinalStateParticles(const reco::GenParticle& particle, const std::set& pdg_to_exclude,
+ const std::set& particles_to_exclude);
+ void addParticle(const reco::GenParticle& particle);
+
+private:
+ std::map> particles;
+ std::map momentum;
+};
+
struct LeptonMatchResult {
GenLeptonMatch match{GenLeptonMatch::NoMatch};
- const reco::GenParticle* gen_particle{nullptr};
- std::vector visible_daughters;
- LorentzVectorXYZ visible_daughters_p4;
- int n_chargedParticles;
- int n_neutralParticles;
+ const reco::GenParticle *gen_particle_firstCopy{nullptr}, *gen_particle_lastCopy{nullptr};
+ std::set visible_daughters, visible_rad;
+ LorentzVectorXYZ visible_p4, visible_rad_p4;
+ unsigned n_charged_hadrons{0}, n_neutral_hadrons{0}, n_gammas{0}, n_gammas_rad{0};
};
-void FindFinalStateDaughters(const reco::GenParticle& particle, std::set& daughters,
- const std::set& pdg_to_exclude);
-
-LorentzVectorXYZ GetFinalStateMomentum(const reco::GenParticle& particle, std::vector& visible_daughters,
- bool excludeInvisible, bool excludeLightLeptons);
-
-LeptonMatchResult LeptonGenMatch(const LorentzVectorM& p4,
- const reco::GenParticleCollection& genParticles);
+const reco::GenParticle* FindTerminalCopy(const reco::GenParticle& genParticle, bool first);
+bool FindLeptonGenMatch(const reco::GenParticle& particle, LeptonMatchResult& result,
+ const LorentzVectorM* ref_p4 = nullptr, double* best_match_dr2 = nullptr);
+std::vector CollectGenLeptons(const reco::GenParticleCollection& genParticles);
+LeptonMatchResult LeptonGenMatch(const LorentzVectorM& p4, const reco::GenParticleCollection& genParticles);
+LeptonMatchResult LeptonGenMatch(const LorentzVectorM& p4, const std::vector& genLeptons);
float GetNumberOfPileUpInteractions(edm::Handle>& pu_infos);
diff --git a/Common/interface/NumericPrimitives.h b/Common/interface/NumericPrimitives.h
new file mode 100644
index 00000000000..5c22bc867bf
--- /dev/null
+++ b/Common/interface/NumericPrimitives.h
@@ -0,0 +1,844 @@
+/*! Definition of the primitives that extend CERN ROOT functionality.
+This file is part of https://github.com/cms-tau-pog/TauTriggerTools. */
+
+#pragma once
+
+#include
+#include
+#include
+#include "TextIO.h"
+
+namespace analysis {
+
+enum class RangeBoundaries { Open, MinIncluded, MaxIncluded, BothIncluded };
+
+namespace detail {
+template::value>
+struct RangeSize {
+ static T size(const T& min, const T& max, RangeBoundaries) { return max - min; }
+};
+
+template
+struct RangeSize {
+ static T size(T min, T max, RangeBoundaries boundaries)
+ {
+ T delta = max - min;
+ if(boundaries == RangeBoundaries::BothIncluded)
+ delta += T(1);
+ else if(boundaries == RangeBoundaries::Open && delta != T(0))
+ delta -= T(1);
+ return delta;
+ }
+};
+
+template::value>
+struct Abs {
+ static T abs(T value) { return std::abs(value); }
+};
+
+template
+struct Abs {
+ static T abs(T value) { return value; }
+};
+
+template
+inline T FloatRound(T x, T /*ref*/)
+{
+ return x;
+}
+
+template<>
+inline float FloatRound(float x, float ref)
+{
+ static constexpr float precision = 1e-6f;
+ if(ref == 0) return x;
+ return std::round(x / ref / precision) * precision * ref;
+}
+
+template<>
+inline double FloatRound(double x, double ref)
+{
+ static constexpr double precision = 1e-8;
+ if(ref == 0) return x;
+ return std::round(x / ref / precision) * precision * ref;
+}
+
+
+template
+inline size_t GetNumberOfGridPoints(T min, T max, T step)
+{
+ return static_cast((max - min) / step);
+}
+
+template<>
+inline size_t GetNumberOfGridPoints(float min, float max, float step)
+{
+ static constexpr float precision = 1e-6f;
+ return static_cast(std::round((max - min) / step / precision) * precision);
+}
+
+template<>
+inline size_t GetNumberOfGridPoints(double min, double max, double step)
+{
+ static constexpr double precision = 1e-8;
+ return static_cast(std::round((max - min) / step / precision) * precision);
+}
+
+} // namespace detail
+
+template
+struct Range {
+ using ValueType = T;
+ using ConstRefType =
+ typename std::conditional::value, ValueType, const ValueType&>::type;
+
+ static const std::pair GetBoundariesSymbols(RangeBoundaries b)
+ {
+ static const std::map> symbols = {
+ { RangeBoundaries::Open, { '(', ')' } },
+ { RangeBoundaries::MinIncluded, { '[', ')' } },
+ { RangeBoundaries::MaxIncluded, { '(', ']' } },
+ { RangeBoundaries::BothIncluded, { '[', ']' } },
+ };
+ return symbols.at(b);
+ }
+
+ static RangeBoundaries CreateBoundaries(bool include_min, bool include_max)
+ {
+ if(include_min && !include_max)
+ return RangeBoundaries::MinIncluded;
+ if(!include_min && include_max)
+ return RangeBoundaries::MaxIncluded;
+ if(include_min && include_max)
+ return RangeBoundaries::BothIncluded;
+ return RangeBoundaries::Open;
+ }
+
+ Range() : _min(0), _max(0) {}
+ Range(ConstRefType min, ConstRefType max, RangeBoundaries boundaries = RangeBoundaries::BothIncluded) :
+ _min(min), _max(max), _boundaries(boundaries)
+ {
+ if(!IsValid(min, max))
+ throw exception("Invalid range [%1%, %2%].") % min % max;
+ }
+ Range(const Range& other) : _min(other._min), _max(other._max), _boundaries(other._boundaries) {}
+ Range(const Range& other, RangeBoundaries boundaries) :
+ _min(other._min), _max(other._max), _boundaries(boundaries) {}
+ virtual ~Range() {}
+ Range& operator=(const Range& other)
+ {
+ _min = other._min;
+ _max = other._max;
+ _boundaries = other._boundaries;
+ return *this;
+ }
+
+ ConstRefType min() const { return _min; }
+ ConstRefType max() const { return _max; }
+ T size() const { return detail::RangeSize::size(min(), max(), boundaries()); }
+
+ RangeBoundaries boundaries() const { return _boundaries; }
+ bool min_included() const
+ {
+ return boundaries() == RangeBoundaries::MinIncluded || boundaries() == RangeBoundaries::BothIncluded;
+ }
+ bool max_included() const
+ {
+ return boundaries() == RangeBoundaries::MaxIncluded || boundaries() == RangeBoundaries::BothIncluded;
+ }
+
+ bool Contains(ConstRefType v) const
+ {
+ if(min() == max())
+ return (min_included() || max_included()) && v == min();
+ const bool min_cond = (min_included() && v >= min()) || v > min();
+ const bool max_cond = (max_included() && v <= max()) || v < max();
+ return min_cond && max_cond;
+ }
+ static bool IsValid(ConstRefType min, ConstRefType max) { return min <= max; }
+
+ Range Extend(ConstRefType v, bool include = true) const
+ {
+ if(Contains(v))
+ return *this;
+ const auto new_min = std::min(min(), v);
+ const auto new_max = std::max(max(), v);
+ RangeBoundaries b;
+ if(new_min == v) {
+ if(include)
+ b = max_included() ? RangeBoundaries::BothIncluded : RangeBoundaries::MinIncluded;
+ else
+ b = max_included() ? RangeBoundaries::MaxIncluded : RangeBoundaries::Open;
+ } else {
+ if(include)
+ b = min_included() ? RangeBoundaries::BothIncluded : RangeBoundaries::MaxIncluded;
+ else
+ b = min_included() ? RangeBoundaries::MinIncluded : RangeBoundaries::Open;
+ }
+ return Range(new_min, new_max, b);
+ }
+
+ bool operator ==(const Range& other) const
+ {
+ return min() == other.min() && max() == other.max() && boundaries() == other.boundaries();
+ }
+ bool operator !=(const Range& other) const { return !(*this == other); }
+
+ bool Includes(const Range& other) const
+ {
+ const bool min_cond = min() == other.min() ? min_included() || !other.min_included() : min() < other.min();
+ const bool max_cond = max() == other.max() ? max_included() || !other.max_included() : max() > other.max();
+ return min_cond && max_cond;
+ }
+ bool Overlaps(const Range& other) const
+ {
+ if(min() == max())
+ return (min_included() || max_included()) && other.Contains(min());
+ const bool min_cond = min() == other.max() ? min_included() && other.max_included() : min() < other.max();
+ const bool max_cond = max() == other.min() ? max_included() && other.min_included() : max() > other.min();
+ return min_cond && max_cond;
+ }
+ Range Combine(const Range& other) const
+ {
+ if(!Overlaps(other))
+ throw exception("Unable to combine non overlapping ranges.");
+ const auto new_min = std::min(min(), other.min());
+ const auto new_max = std::max(max(), other.max());
+ const bool include_min = (new_min == min() && min_included())
+ || (new_min == other.min() && other.min_included());
+ const bool include_max = (new_max == max() && max_included())
+ || (new_max == other.max() && other.max_included());
+ const RangeBoundaries b = CreateBoundaries(include_min, include_max);
+ return Range(new_min, new_max, b);
+ }
+
+ std::string ToString(char sep = ':') const
+ {
+ std::ostringstream ss;
+ const auto b_sym = GetBoundariesSymbols(boundaries());
+ if(boundaries() != RangeBoundaries::BothIncluded)
+ ss << b_sym.first;
+ ss << min() << sep << max();
+ if(boundaries() != RangeBoundaries::BothIncluded)
+ ss << b_sym.second;
+ return ss.str();
+ }
+
+ static Range Parse(const std::string& str, const std::string& separators=": \t")
+ {
+ const auto values = SplitValueList(str, true, separators, true);
+ if(values.size() != 2)
+ throw exception("Invalid range '%1%'.") % str;
+ return Make(values);
+ }
+
+ static Range Read(std::istream& stream, const std::string& separators=": \t")
+ {
+ const auto values = ReadValueList(stream, 2, true, separators, true);
+ return Make(values);
+ }
+
+private:
+ static Range Make(const std::vector& values)
+ {
+ static const auto opened_b_symbols = GetBoundariesSymbols(RangeBoundaries::Open);
+ static const auto closed_b_symbols = GetBoundariesSymbols(RangeBoundaries::BothIncluded);
+ bool include_min = true, include_max = true;
+ std::string min_str = values.at(0), max_str = values.at(1);
+ if(min_str.size() && (min_str.front() == opened_b_symbols.first || min_str.front() == closed_b_symbols.first)) {
+ include_min = min_str.front() == closed_b_symbols.first;
+ min_str.erase(0, 1);
+ }
+ if(max_str.size() && (max_str.back() == opened_b_symbols.second || max_str.back() == closed_b_symbols.second)) {
+ include_max = max_str.back() == closed_b_symbols.second;
+ max_str.erase(max_str.size() - 1, 1);
+ }
+ const T min = ::analysis::Parse(min_str);
+ const T max = ::analysis::Parse(max_str);
+ const RangeBoundaries b = CreateBoundaries(include_min, include_max);
+ return Range(min, max, b);
+ }
+
+private:
+ T _min, _max;
+ RangeBoundaries _boundaries;
+};
+
+template
+std::ostream& operator<<(std::ostream& s, const Range& r)
+{
+ s << r.ToString(':');
+ return s;
+}
+
+template
+std::istream& operator>>(std::istream& s, Range& r)
+{
+ r = Range::Read(s);
+ return s;
+}
+
+template
+struct RelativeRange {
+ using ValueType = T;
+ using ConstRefType = typename Range::ConstRefType;
+ RelativeRange() : _down(0), _up(0) {}
+ RelativeRange(ConstRefType down, ConstRefType up) : _down(down), _up(up)
+ {
+ if(!IsValid(down, up))
+ throw exception("Invalid relative range [%1%, %2%].") % down % up;
+ }
+
+ ConstRefType down() const { return _down; }
+ ConstRefType up() const { return _up; }
+ Range ToAbsoluteRange(ConstRefType v) const { return Range(v + down(), v + up()); }
+ static bool IsValid(ConstRefType down, ConstRefType up) { return down <= 0 && up >= 0; }
+
+ std::string ToString(char sep = ' ') const
+ {
+ std::ostringstream ss;
+ ss << down() << sep << up();
+ return ss.str();
+ }
+
+ static RelativeRange Parse(const std::string& str, const std::string& separators=": \t")
+ {
+ const auto values = SplitValueList(str, true, separators, true);
+ if(values.size() != 2)
+ throw exception("Invalid relative range '%1%'.") % str;
+ return Make(values);
+ }
+
+ static RelativeRange Read(std::istream& stream, const std::string& separators=": \t")
+ {
+ const auto values = ReadValueList(stream, 2, true, separators, true);
+ return Make(values);
+ }
+
+private:
+ static RelativeRange Make(const std::vector& values)
+ {
+ const T down = ::analysis::Parse(values.at(0));
+ const T up = ::analysis::Parse(values.at(1));
+ return RelativeRange(down, up);
+ }
+
+private:
+ T _down, _up;
+};
+
+template
+std::ostream& operator<<(std::ostream& s, const RelativeRange& r)
+{
+ s << r.ToString(':');
+ return s;
+}
+
+template
+std::istream& operator>>(std::istream& s, RelativeRange& r)
+{
+ r = RelativeRange::Read(s);
+ return s;
+}
+
+template
+struct RangeWithStep : public Range {
+ using ValueType = typename Range::ValueType;
+ using ConstRefType = typename Range::ConstRefType;
+
+
+ enum class PrintMode { Step = 0, NGridPoints = 1, NBins = 2 };
+ struct iterator {
+ iterator(const RangeWithStep& _range, size_t _pos) : range(&_range), pos(_pos) {}
+ iterator& operator++() { ++pos; return *this; }
+ iterator operator++(int) { iterator iter(*this); operator++(); return iter; }
+ bool operator==(const iterator& other) { return range == other.range && pos == other.pos;}
+ bool operator!=(const iterator& other) { return !(*this == other); }
+ T operator*() { return range->grid_point_value(pos); }
+ private:
+ const RangeWithStep *range;
+ size_t pos;
+ };
+
+ RangeWithStep() : _step(0) {}
+ RangeWithStep(ConstRefType min, ConstRefType max, ConstRefType step) :
+ Range(min, max, RangeBoundaries::BothIncluded), _step(step)
+ {
+ }
+
+ ConstRefType step() const { return _step; }
+ T grid_point_value(size_t index) const
+ {
+ const T ref = std::max(detail::Abs::abs(this->min()), detail::Abs::abs(this->max()));
+ return detail::FloatRound(this->min() + T(index) * step(), ref);
+ }
+ size_t n_grid_points() const
+ {
+ if(this->max() == this->min()) return 1;
+ if(step() == 0)
+ throw exception("Number of grid points is not defined for a non-point range with the step = 0.");
+ size_t n_points = detail::GetNumberOfGridPoints(this->min(), this->max(), step());
+ if(this->Contains(grid_point_value(n_points)))
+ ++n_points;
+ return n_points;
+ }
+ size_t n_bins() const { return n_grid_points() - 1; }
+
+ size_t find_bin(T value) const
+ {
+ if(!this->Contains(value))
+ throw exception("find_bin: value is out of range.");
+ if(n_bins() == 0)
+ throw exception("find_bin: number of bins is 0.");
+ size_t bin_id = static_cast((value - this->min()) / step());
+ if(bin_id == n_bins())
+ --bin_id;
+ return bin_id;
+ }
+
+ iterator begin() const { return iterator(*this, 0); }
+ iterator end() const { return iterator(*this, n_grid_points()); }
+
+ std::string ToString(PrintMode mode = PrintMode::Step) const
+ {
+ std::ostringstream ss;
+ ss << this->min() << Separators().at(0) << this->max() << Separators().at(static_cast(mode));
+ if(mode == PrintMode::Step)
+ ss << step();
+ else if(mode == PrintMode::NGridPoints)
+ ss << n_grid_points();
+ else if(mode == PrintMode::NBins)
+ ss << n_bins();
+ else
+ throw exception("Unsupported RangeWithStep::PrintMode = %1%.") % static_cast(mode);
+ return ss.str();
+ }
+
+ static RangeWithStep Parse(const std::string& str)
+ {
+ const size_t first_split_pos = str.find_first_of(Separators());
+ if(first_split_pos != std::string::npos) {
+ const size_t last_split_pos = str.find_first_of(Separators(), first_split_pos + 1);
+ if(last_split_pos != std::string::npos) {
+ const size_t end_split_pos = str.find_last_of(Separators());
+ if(last_split_pos == end_split_pos) {
+ const size_t sep_pos = Separators().find(str.at(last_split_pos));
+ const PrintMode mode = static_cast(sep_pos);
+ std::vector values;
+ values.push_back(str.substr(0, first_split_pos));
+ values.push_back(str.substr(first_split_pos + 1, last_split_pos - first_split_pos - 1));
+ values.push_back(str.substr(last_split_pos + 1));
+ return Make(values, mode);
+ }
+ }
+ }
+ throw exception("Invalid range with step '%1%'.") % str;
+ }
+
+private:
+ static RangeWithStep Make(const std::vector& values, PrintMode mode)
+ {
+ const T min = ::analysis::Parse(values.at(0));
+ const T max = ::analysis::Parse(values.at(1));
+ T step(0);
+ if(mode == PrintMode::Step) {
+ step = ::analysis::Parse(values.at(2));
+ } else if(mode == PrintMode::NGridPoints) {
+ size_t n = ::analysis::Parse(values.at(2));
+ if(n == 0 || (n == 1 && max != min) || (n != 1 && max == min))
+ throw exception("Invalid number of grid points.");
+ if(max != min)
+ step = (max - min) / T(n - 1);
+ } else if(mode == PrintMode::NBins) {
+ size_t n = ::analysis::Parse(values.at(2));
+ if((n == 0 && max != min) || (n != 0 && max == min))
+ throw exception("Invalid number of bins.");
+ if(max != min)
+ step = (max - min) / T(n);
+ } else {
+ throw exception("Unsupported RangeWithStep::PrintMode = %1%.") % static_cast(mode);
+ }
+ return RangeWithStep(min, max, step);
+ }
+
+ static const std::string& Separators() { static const std::string sep = ":|/"; return sep; }
+
+private:
+ T _step;
+};
+
+template
+std::ostream& operator<<(std::ostream& s, const RangeWithStep& r)
+{
+ s << r.ToString();
+ return s;
+}
+
+template
+std::istream& operator>>(std::istream& s, RangeWithStep& r)
+{
+ std::string str;
+ s >> str;
+ r = RangeWithStep::Parse(str);
+ return s;
+}
+
+template
+struct Angle {
+ enum class Interval { Symmetric, Positive };
+ static constexpr double Pi() { return boost::math::constants::pi(); }
+ static constexpr double NumberOfPiPerPeriod() { return double(n_pi_per_period_num) / n_pi_per_period_denom; }
+ static constexpr double FullPeriod() { return n_pi_per_period_num * Pi() / n_pi_per_period_denom; }
+ static constexpr double HalfPeriod() { return FullPeriod() / 2; }
+ static constexpr double RadiansToDegreesFactor() { return 180. / Pi(); }
+
+ Angle() : _value(0), _interval(Interval::Symmetric) {}
+ Angle(double value, Interval interval = Interval::Symmetric)
+ : _value(AdjustValue(value, interval)), _interval(interval) {}
+
+ double value() const { return _value; }
+ double value_degrees() const { return value() * RadiansToDegreesFactor(); }
+ Interval interval() const { return _interval; }
+
+ Angle operator+(const Angle& other) const { return Angle(value() + other.value(), interval()); }
+ Angle operator-(const Angle& other) const { return Angle(value() - other.value(), interval()); }
+
+ static const Range& AngleValuesRange(Interval interval)
+ {
+ static const std::map> range_map = {
+ { Interval::Symmetric, { -HalfPeriod(), HalfPeriod() } },
+ { Interval::Positive, { 0, FullPeriod() } }
+ };
+ return range_map.at(interval);
+ }
+
+ static double AdjustValue(double value, Interval interval)
+ {
+ const Range& range = AngleValuesRange(interval);
+ value -= FullPeriod() * std::floor(value/FullPeriod());
+ while(value < range.min() || value >= range.max())
+ value += value < range.min() ? FullPeriod() : -FullPeriod();
+ return value;
+ }
+
+private:
+ double _value;
+ Interval _interval;
+};
+
+template
+std::ostream& operator<<(std::ostream& s, const Angle& a)
+{
+ s << a.value();
+ return s;
+}
+
+template
+std::istream& operator>>(std::istream& s, Angle& a)
+{
+ double value;
+ s >> value;
+ a = Angle(value);
+ return s;
+}
+
+template
+struct Range> {
+ using A = Angle;
+ using ValueType = A;
+
+ Range() : _min(0), _max(0) {}
+ Range(const A& min, const A& max) : _min(min), _max(max.value(), min.interval()) {}
+ virtual ~Range() {}
+
+ const A& min() const { return _min; }
+ const A& max() const { return _max; }
+ A size() const { return A(_max.value() - _min.value(), A::Interval::Positive); }
+ Range ToValueRange() const
+ {
+ const double min_value = min().value();
+ double max_value = max().value();
+ if(max_value < min_value)
+ max_value += A::FullPeriod();
+ return Range(min_value, max_value);
+ }
+
+ bool Contains(const A& a) const
+ {
+ const Range min_a_value_range = Range(min(), a).ToValueRange();
+ return ToValueRange().Contains(min_a_value_range.max());
+ }
+
+ static bool IsValid(const A& /*min*/, const A& /*max*/) { return true; }
+
+ Range Extend(const A& a) const
+ {
+ if(Contains(a))
+ return *this;
+ const A a_fixed(a.value(), min().interval());
+ const Range extend_min(a_fixed, max()), extend_max(min(), a_fixed);
+ return extend_max.size().value() < extend_min.size().value() ? extend_max : extend_min;
+ }
+
+ bool Includes(const Range& other) const
+ {
+ return Contains(other.min()) && Contains(other.max());
+ }
+
+ bool Overlaps(const Range& other) const
+ {
+ return Contains(other.min()) || Contains(other.max()) || other.Contains(min());
+ }
+
+ Range Combine(const Range& other) const
+ {
+ if(!Overlaps(other))
+ throw exception("Unable to combine non overlapping ranges.");
+ if(Includes(other))
+ return *this;
+ if(other.Includes(*this))
+ return other;
+ if(Contains(other.min()))
+ return Range(min(), other.max());
+ return Range(other.min(), max());
+ }
+ std::string ToString(char sep = ' ') const
+ {
+ std::ostringstream ss;
+ ss << min() << sep << max();
+ return ss.str();
+ }
+
+ static Range Parse(const std::string& str, const std::string& separators=": \t")
+ {
+ const auto values = SplitValueList(str, true, separators, true);
+ if(values.size() != 2)
+ throw exception("Invalid angle range '%1%'.") % str;
+ return Make(values);
+ }
+
+ static Range Read(std::istream& stream, const std::string& separators=": \t")
+ {
+ const auto values = ReadValueList(stream, 2, true, separators, true);
+ return Make(values);
+ }
+
+private:
+ static Range Make(const std::vector& values)
+ {
+ const A min = ::analysis::Parse(values.at(0));
+ const A max = ::analysis::Parse(values.at(1));
+ return Range(min, max);
+ }
+
+private:
+ A _min, _max;
+};
+
+template
+struct RangeMultiD {
+public:
+ using ValueType = typename Range::ValueType;
+ explicit RangeMultiD(size_t n_dim) : ranges(n_dim) {}
+ explicit RangeMultiD(const std::vector& _ranges) : ranges(_ranges) {}
+
+ size_t GetNumberOfDimensions() const { return ranges.size(); }
+ const Range& GetRange(size_t dim_id) const { Check(dim_id); return ranges.at(dim_id - 1); }
+ Range& GetRange(size_t dim_id) { Check(dim_id); return ranges.at(dim_id - 1); }
+
+ bool Contains(const std::vector& point) const
+ {
+ if(point.size() != GetNumberOfDimensions())
+ throw exception("Invalid number of dimensions.");
+ for(size_t n = 0; n < ranges.size(); ++n)
+ if(!ranges.at(n).Contains(point.at(n))) return false;
+ return true;
+ }
+
+private:
+ void Check(size_t dim_id) const
+ {
+ if(!dim_id || dim_id > GetNumberOfDimensions())
+ throw exception("Wrong dimension id = %1%") % dim_id;
+ }
+
+private:
+ std::vector ranges;
+};
+
+template
+struct MultiRange {
+ using ValueType = typename Range::ValueType;
+ using ConstRefType = typename Range::ConstRefType;
+ using RangeVec = std::vector;
+
+ static const std::string& Separator() { static const std::string sep = ", "; return sep; }
+
+ MultiRange() {}
+ explicit MultiRange(const RangeVec& _ranges) : ranges(_ranges) {}
+
+ bool Contains(const ValueType& point) const
+ {
+ for(const auto& range : ranges) {
+ if(range.Contains(point))
+ return true;
+ }
+ return false;
+ }
+
+ bool Overlaps(const Range& other) const
+ {
+ for(const auto& range : ranges) {
+ if(range.Overlaps(other))
+ return true;
+ }
+ return false;
+ }
+
+ std::string ToString() const
+ {
+ std::ostringstream ss;
+ for(const auto& range : ranges)
+ ss << range << Separator();
+ std::string str = ss.str();
+ if(str.size())
+ str.erase(str.size() - Separator().size());
+ return str;
+ }
+
+ static MultiRange Parse(const std::string& str)
+ {
+ const auto range_strs = SplitValueList(str, true, Separator(), true);
+ RangeVec ranges;
+ for(const auto& range_str : range_strs)
+ ranges.push_back(::analysis::Parse(range_str));
+ return MultiRange(ranges);
+ }
+
+private:
+ RangeVec ranges;
+};
+
+template
+std::ostream& operator<<(std::ostream& s, const MultiRange& r)
+{
+ s << r.ToString();
+ return s;
+}
+
+template
+std::istream& operator>>(std::istream& s, MultiRange& r)
+{
+ std::string str;
+ std::getline(s, str);
+ r = MultiRange::Parse(str);
+ return s;
+}
+
+
+struct NumericalExpression {
+ NumericalExpression() : _value(0) {}
+ NumericalExpression(const std::string& expression)
+ : _expression(expression)
+ {
+ const std::string formula = boost::str(boost::format("x*(%1%)") % expression);
+ TF1 fn("", formula.c_str(), 0, 1);
+// if(!fn.IsValid())
+// throw exception("Invalid numerical expression '%1%'") % expression;
+ _value = fn.Eval(1);
+ }
+
+ const std::string& expression() const { return _expression; }
+ double value() const { return _value; }
+ operator double() const { return _value; }
+
+private:
+ std::string _expression;
+ double _value;
+};
+
+inline std::ostream& operator<<(std::ostream& s, const NumericalExpression& e)
+{
+ s << e.expression();
+ return s;
+}
+
+inline std::istream& operator>>(std::istream& s, NumericalExpression& e)
+{
+ std::string line;
+ std::getline(s, line);
+ e = NumericalExpression(line);
+ return s;
+}
+
+struct Grid_ND {
+ using Position = std::vector;
+
+ struct iterator {
+ iterator(const Position& _pos, const Position& _limits) : pos(_pos), limits(&_limits) {}
+
+ bool operator==(const iterator& other) const {
+ if(pos.size() != other.pos.size()) return false;
+ for(size_t n = 0; n < pos.size(); ++n) {
+ if(pos.at(n) != other.pos.at(n)) return false;
+ }
+ return true;
+ }
+
+ bool operator!=(const iterator& other) { return !(*this == other); }
+
+ iterator& operator++()
+ {
+ ++pos.at(0);
+ for(size_t n = 0; n < pos.size() - 1 && pos.at(n) >= limits->at(n); ++n) {
+ ++pos.at(n+1);
+ pos.at(n) = 0;
+ }
+ return *this;
+ }
+
+ iterator operator++(int)
+ {
+ iterator cp(*this);
+ ++(*this);
+ return cp;
+ }
+
+ const Position& operator*() const { return pos; }
+ const Position* operator->() const { return &pos; }
+
+ private:
+ Position pos;
+ const Position* limits;
+ };
+
+ explicit Grid_ND(const Position& _limits) : limits(_limits)
+ {
+ if(!limits.size())
+ throw exception("Grid dimensions should be > 0");
+ for(size_t limit : limits) {
+ if(!limit)
+ throw exception("Grid range limit should be > 0.");
+ }
+ }
+
+ iterator begin() const
+ {
+ Position pos;
+ pos.assign(limits.size(), 0);
+ return iterator(pos, limits);
+ }
+
+ iterator end() const
+ {
+ Position pos;
+ pos.assign(limits.size(), 0);
+ pos.back() = limits.back();
+ return iterator(pos, limits);
+ }
+
+private:
+ Position limits;
+};
+
+} // namespace analysis
diff --git a/Common/interface/PatHelpers.h b/Common/interface/PatHelpers.h
new file mode 100644
index 00000000000..5f104ec3859
--- /dev/null
+++ b/Common/interface/PatHelpers.h
@@ -0,0 +1,39 @@
+/*! Various utility functions.
+This file is part of https://github.com/cms-tau-pog/TauTriggerTools. */
+
+#pragma once
+
+#include