Adds basic event selection

analysis/TTTT_Analysis.cpp

@@ -70,10 +70,6 @@ void enable_extra_branches(MiniTreeDataSet& mt){
 void declare_values(MiniTreeDataSet& mt){
 
     // Define a couple selections to be used in the top-mass reconstruction.
-    auto& b_mva_filter = GenFunction::register_function<bool(Particle)>("b_mva_filter",
-        FUNC(([cut=B_JET_WP](const Particle& j){
-                return j.jet.b_cmva > cut;
-        })));
     auto& b_pdgid_filter = GenFunction::register_function<bool(Particle)>("b_pdgid_filter",
         FUNC(([](const Particle& j){
                 return j.genpart.pdgId == 5 || j.genpart.pdgId==-5;
@@ -102,7 +98,7 @@ void declare_values(MiniTreeDataSet& mt){
     // Here is the calculation of the Top Reconstructed Mass from Particle
     auto jets = lookup<std::vector<Particle>>("jets");
 
-    auto b_jets = filter(b_mva_filter, jets, "reco_b_jets");
+    auto b_jets = lookup<std::vector<Particle>>("b_jets");
     auto w_dijets = tup_filter<Particle,Particle>(w_mass_filter, combinations<Particle,2>(jets, "reco_dijets"));
 
     auto top_cands = cart_product<std::tuple<Particle,Particle>, Particle>(w_dijets, b_jets);
@@ -301,6 +297,10 @@ void declare_containers(MiniTreeDataSet& mt){
     mt.register_container<Vector<std::vector<float>>>("Jet_phi",             lookup<std::vector<float>>("Jet_phi"));
     mt.register_container<Vector<std::vector<float>>>("Jet_mass",            lookup<std::vector<float>>("Jet_mass"));
 
+    mt.register_container<PassCount>("SR4j_count", SR4j);
+    mt.register_container<PassCount>("SR5j_count", SR5j);
+    mt.register_container<PassCount>("SR6j_count", SR6j);
+
 }
 
 
@@ -318,6 +318,7 @@ void run_analysis(const std::string& input_filename, bool silent){
     create_all_common_values(mt);
     enable_extra_branches(mt);
     declare_values(mt);
+    init_selection();
     declare_containers(mt);
 
     mt.process(silent);

analysis/common/obj_types.hpp

@@ -129,6 +129,16 @@ construct_jets_value(MiniTreeDataSet& mt){
     return jets;
 }
 
+decltype(auto)
+construct_b_jets_value(MiniTreeDataSet& mt){
+    auto b_jets = fv::filter(GenFunction::register_function<bool(Particle)>("is_b_jet",
+        FUNC(([](const Particle& p){
+            if(p.tag != Particle::JET) return false;
+            return p.jet.b_cmva > 0;
+        }))), lookup<std::vector<Particle>>("jets"), "b_jets");
+    return b_jets;
+}
+
 
 decltype(auto)
 construct_leptons_value(MiniTreeDataSet& mt){
@@ -205,6 +215,7 @@ construct_mc_jets_value(MiniTreeDataSet& mt){
 
 void create_all_common_values(MiniTreeDataSet& mt){
     construct_jets_value(mt);
+    construct_b_jets_value(mt);
     construct_mc_jets_value(mt);
     construct_leptons_value(mt);
 }

analysis/selection.hpp

@@ -34,54 +34,77 @@
 #ifndef SELECTION_HPP
 #define SELECTION_HPP
 
-#include <iostream>
 #include <vector>
-#include <utility>
-#include <numeric>
-#include <limits>
 
 #include "filval/filval.hpp"
 #include "filval_root/filval_root.hpp"
 
-#include <TSystem.h>
+#include "analysis/common/obj_types.hpp"
 
-#define JET_REQUIREMENT   6
-#define B_JET_REQUIREMENT 3
-#define B_JET_WP          0.  // lower bound on CMVA value
-using namespace fv;
-struct Selection {
-public:
+ObsFilter* SR4j;
+ObsFilter* SR5j;
+ObsFilter* SR6j;
 
-    static ObsFilter* trilepton_filter;
+void init_selection(){
+    auto leptons = lookup<std::vector<Particle>>("leptons");
+    auto jets    = lookup<std::vector<Particle>>("jets");
 
-    static ObsFilter* jet_multiplicity_filter;
-    static ObsFilter* b_jet_multiplicity_filter;
+    // Require *exactly* three charged leptons
+    auto trilepton = obs_filter("trilepton",GenFunction::register_function<bool()>("trilepton",
+        FUNC(([leptons](){
+            return leptons->get_value().size() == 3;
+        }))));
 
-    static void init(){
-        auto nLep = lookup<int>("nLepGood");
-        trilepton_filter = obs_filter("trilepton_filter",GenFunction::register_function<bool()>("trilepton_filter",
-            FUNC(([nLep=nLep](){
-                return nLep->get_value() == 3;
-            }))));
+    // Require at least three charged b-jets
+    auto b_jet = obs_filter("b_jet",GenFunction::register_function<bool()>("b_jet",
+        FUNC(([jets](){
+            int n_b_jets = 0;
+            for(auto jet : jets->get_value()){
+                if(jet.tag == Particle::JET && jet.jet.b_cmva > 0)
+                    n_b_jets++;
+            }
+            return n_b_jets >= 3;
+        }))));
 
+    // Require all same-flavour OS dilepton combinations are outside the Z mass
+    // window (70,105)
+    auto z_mass_window = obs_filter("z_mass_window",GenFunction::register_function<bool()>("z_mass_window",
+        FUNC(([leptons](){
+            auto& leps = leptons->get_value();
+            int n = leps.size();
+            for(int i = 0; i < n; i++){
+                for(int j = i+1; j < n; j++){
+                    const Particle& p1 = leps[i];
+                    const Particle& p2 = leps[j];
+                    if(abs(p1.lepton.pdg_id) != abs(p2.lepton.pdg_id)) continue;
+                    if(p1.lepton.charge == p2.lepton.charge) continue;
+                    double m = (p1.v + p2.v).M();
+                    if(70 < m && m < 105)
+                        return false;
+                }
+            }
+            return true;
+        }))));
 
-        auto nJet = lookup<int>("nJet");
-        jet_multiplicity_filter = obs_filter("jet_multiplicity_filter",GenFunction::register_function<bool()>("jet_multiplicity_filter",
-            FUNC(([nJet=nJet](){
-                return nJet->get_value() >= JET_REQUIREMENT;
-            }))));
 
-        auto bJet_MVA = lookup<std::vector<float>>("Jet_btagCMVA");
-        b_jet_multiplicity_filter = obs_filter("b_jet_multiplicity_filter",GenFunction::register_function<bool()>("b_jet_multiplicity_filter",
-            FUNC(([bJet_MVA=bJet_MVA](){
-                int n_b_jet = 0;
-                for(auto j : bJet_MVA->get_value()){
-                    if(j > B_JET_WP){
-                        n_b_jet++;
-                    }
-                }
-                return n_b_jet >= B_JET_REQUIREMENT;
-            }))));
-    }
-};
+    auto J4 = obs_filter("4jet_selection",GenFunction::register_function<bool()>("4jet_selection",
+        FUNC(([jets](){
+            return jets->get_value().size() >= 4;
+        }))));
+
+    auto J5 = obs_filter("5jet_selection",GenFunction::register_function<bool()>("5jet_selection",
+        FUNC(([jets](){
+            return jets->get_value().size() >= 5;
+        }))));
+
+    auto J6 = obs_filter("6jet_selection",GenFunction::register_function<bool()>("6jet_selection",
+        FUNC(([jets](){
+            return jets->get_value().size() >= 6;
+        }))));
+
+    auto base_sel = ObsFilter::conj(z_mass_window, ObsFilter::conj(trilepton, b_jet));
+    SR4j = ObsFilter::conj(base_sel, J4);
+    SR5j = ObsFilter::conj(base_sel, J5);
+    SR6j = ObsFilter::conj(base_sel, J6);
+}
 #endif // SELECTION_HPP

filval/filter.hpp

@@ -56,8 +56,21 @@ class ObsFilter : public DerivedValue<bool>{
           :DerivedValue<bool>(name),
            filter_function(GenFunction::register_function<bool()>("filter::"+name, filter_function, impl)){ }
 
-        /** Return a new filter that is the conjuction of the two source
-         * filters.
+        /** Return a new filter that is the conjuction of the two source filters.
+         */
+        static ObsFilter* conj(ObsFilter *f1, ObsFilter *f2){
+            auto new_name = f1->get_name() + "&&" + f2->get_name();
+            return new ObsFilter(new_name, [f1,f2](){return f1->get_value() && f2->get_value();});
+        }
+
+        /** Return a new filter that is the conjuction of the two source filters.
+         */
+        static ObsFilter* disj(ObsFilter *f1, ObsFilter *f2){
+            auto new_name = f1->get_name() + "||" + f2->get_name();
+            return new ObsFilter(new_name, [f1, f2](){return f1->get_value() && f2->get_value();});
+        }
+
+        /** Return a new filter that is the conjuction of the two source filters.
          */
         ObsFilter* operator*(ObsFilter *f){
             auto new_name = this->get_name() + "&&" + f->get_name();

filval_root/container.hpp

@@ -73,6 +73,7 @@ namespace fv::root::util{
                 INFO("Cannot save STL container " << type_name <<" as pdf");
                 break;
             case ROOT:
+                /* DEBUG("Writing object \"" << obj_name << "\" of type \"" << type_name << "\"\n"); */
                 gDirectory->WriteObjectAny(container, type_name.c_str(), obj_name.c_str());
                 break;
             default:
@@ -276,7 +277,6 @@ class _Counter : public Container<std::map<D,int>,V>{
             std::string type_name = "std::map<"+fv::util::get_type_name(typeid(D))+",int>";
             util::save_as_stl(this->get_container(), type_name, this->get_name(), option);
         }
-
 };
 
 
@@ -304,6 +304,28 @@ class CounterMany : public _Counter<std::vector<V>,V>{
         }
 };
 
+class PassCount : public Container<int,bool>{
+    private:
+
+        void _fill(){
+            if(this->value->get_value()){
+                (*this->container)++;
+            }
+        }
+    public:
+        PassCount(const std::string& name, Value<bool>* value)
+          :Container<int,bool>(name, value){
+            this->container = new int(0);
+        }
+
+        void save_as(const std::string& fname, const SaveOption& option = SaveOption::PNG) {
+            //ROOT(hilariously) cannot serialize basic data types, we wrap this
+            //in a vector.
+            std::vector<int> v({*this->get_container()});
+            util::save_as_stl(&v, "std::vector<int>", this->get_name(), option);
+        }
+};
+
 
 template <typename... ArgTypes>
 class MVA : public Container<TMVA::DataLoader,typename MVAData<ArgTypes...>::type>{

python/utils.py

@@ -278,7 +278,7 @@ class ResultSet:
             CANVAS.SetCanvasSize(SINGLE_PLOT_SIZE[0],
                                  SINGLE_PLOT_SIZE[1])
 
-        colors = it.cycle([ROOT.kRed, ROOT.kBlue, ROOT.kGreen])
+        colors = it.cycle([ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kYellow])
         stack = ROOT.THStack(hist_name+"_stack", title)
         if labels is None:
             labels = [hist.GetName() for hist in hists]