Added new tuple structure along with new methods of computing Tau Trigger SFs

.gitignore

@@ -3,8 +3,10 @@
 *.so
 *.pdf
 *.zip
-fit.exe
 *~
-*#*
-*.d
 __init__.py
+output
+*.code-workspace
+/.python
+/.vscode
+.DS_Store

Common/BuildFile.xml

@@ -1,4 +1,7 @@
 <use name="root"/>
 <use name="boost"/>
+<use name="CLHEP"/>
 <use name="DataFormats/HepMCCandidate"/>
+<use name="DataFormats/PatCandidates"/>
+<use name="FWCore/ParameterSet"/>
 <export> <lib name="1"/> </export>

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

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 <exception>
+#include <stdexcept>
+#include <string>
+#include <sstream>
+#include <vector>
+#include <set>
+#include <iostream>
+
+#include <TH1D.h>
+#include <Rtypes.h>
+
+#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<typename ValueType, typename Histogram>
+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<typename ObjectType, typename Selector, typename Comparitor>
+    std::vector<ObjectType> collect_objects(double weight, size_t n_objects, const Selector& selector,
+                                            const Comparitor& comparitor)
+    {
+        std::vector<ObjectType> 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<unsigned> counters;
+    std::vector<double> selections;
+    std::vector<double> selectionsSquaredErros;
+    std::vector<std::string> labels;
+    std::set<std::string> label_set;
+};
+
+namespace detail {
+struct DefaultSelectionManager{
+    template<typename ValueType>
+    void FillHistogram(ValueType value, const std::string& histogram_name) {}
+};
+} // namespace detail
+
+template<typename SelectionManager = detail::DefaultSelectionManager>
+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<typename ValueType>
+    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<cuts::ObjectSelector> : 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<TH1D> 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};
+};
+}

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<int>& pdg_to_exclude = {},
+                        const std::set<const reco::GenParticle*>& particles_to_exclude = {});
+
+    const std::set<const reco::GenParticle*>& 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<int>& pdg_to_exclude,
+                                 const std::set<const reco::GenParticle*>& particles_to_exclude);
+    void addParticle(const reco::GenParticle& particle);
+
+private:
+    std::map<ParticleType, std::set<const reco::GenParticle*>> particles;
+    std::map<ParticleType, LorentzVectorXYZ> momentum;
+};
+
 struct LeptonMatchResult {
     GenLeptonMatch match{GenLeptonMatch::NoMatch};
-    const reco::GenParticle* gen_particle{nullptr};
-    std::vector<const reco::GenParticle*> 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<const reco::GenParticle*> 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<const reco::GenParticle*>& daughters,
-                             const std::set<int>& pdg_to_exclude);
-
-LorentzVectorXYZ GetFinalStateMomentum(const reco::GenParticle& particle, std::vector<const reco::GenParticle*>& 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<LeptonMatchResult> CollectGenLeptons(const reco::GenParticleCollection& genParticles);
+LeptonMatchResult LeptonGenMatch(const LorentzVectorM& p4, const reco::GenParticleCollection& genParticles);
+LeptonMatchResult LeptonGenMatch(const LorentzVectorM& p4, const std::vector<LeptonMatchResult>& genLeptons);
 
 float GetNumberOfPileUpInteractions(edm::Handle<std::vector<PileupSummaryInfo>>& pu_infos);