Introduces the filval analysis fromwork wip

.gitignore

@@ -6,3 +6,8 @@ __pycache__/
 *.d
 *.so
 *.pcm
+
+# output data files
+*.root
+*.png
+*.pdf

MiniTree.cpp

@@ -1,159 +0,0 @@
-#ifndef minitree_cxx
-#define minitree_cxx
-
-#include <algorithm>
-#include <functional>
-#include <cstdio>
-#include <string>
-#include <vector>
-#include <map>
-
-#include <TStyle.h>
-
-#include "MiniTree.h"
-#include "HistCollection.h"
-using namespace std;
-
-
-HistCollection* build_histograms(const char* sample_name, const char* filename){
-    HistCollection *hc = new HistCollection(sample_name);
-    TFile *f = TFile::Open(filename);
-    TTree *tree = (TTree*) f->Get("tree");
-    MiniTree minitree(tree);
-    minitree.Loop(hc);
-    delete tree;
-    f->Close();
-    delete f;
-    printf("Finished filling histograms for file %s\n", filename);
-    return hc;
-}
-
-void set_branch_statuses(MiniTree *minitree){
-    auto on = [minitree](const char* bname){
-        minitree->fChain->SetBranchStatus(bname, true);
-    };
-    auto off = [minitree](const char* bname){
-        minitree->fChain->SetBranchStatus(bname, false);
-    };
-    off("*");
-    on("nLepGood");
-    on("LepGood_pt");
-    on("LepGood_eta");
-    on("LepGood_phi");
-    on("LepGood_mcPt");
-    on("LepGood_charge");
-    on("nJet");
-    on("Jet_btagCMVA");
-    on("Jet_pt");
-    on("nGenTop");
-}
-
-double delta_r2(double eta1, double eta2, double phi1, double phi2){
-    return pow(eta1 - eta2, 2) + pow(phi1 - phi2, 2);
-}
-
-double isolation(int n, int i, function<double(int, int)> delta_r2_calculator){
-    double min_iso = -1;
-    for (int j=0; j<n; j++){
-        if (i == j) continue;
-        double curr_iso = delta_r2_calculator(i, j);
-        if (curr_iso < min_iso || min_iso == -1){
-            min_iso = curr_iso;
-        }
-    }
-    if (min_iso == -1)
-        return min_iso;
-    else
-        return sqrt(min_iso);
-}
-
-void MiniTree::Loop(HistCollection* hc){
-
-    /* Define first a few little subroutines that will be used later in the
-     * main loop.
-     */
-    auto lepton_delta_r2 = [this](int i, int j){
-        return delta_r2(this->LepGood_eta[i], this->LepGood_eta[j],
-                        this->LepGood_phi[i], this->LepGood_phi[j]);
-    };
-    auto lepton_isolation = [this, lepton_delta_r2](int i){
-        return isolation(this->nLepGood, i, function<double(int, int)>(lepton_delta_r2));
-    };
-    auto mini_isolation_cut = [this, lepton_isolation](int i){
-        /* Implement the mini-isolation cut which utilizes a dynamic isolation cone size
-         * based on lepton pt
-         */
-        double r0 = 0.5;
-        double pt0 = 7.5; // GeV
-
-        double iso_cut = min(r0, pt0/this->LepGood_pt[i]);
-        double isolation = lepton_isolation(i);
-        return isolation > iso_cut;
-    };
-
-    if (fChain == 0) return;
-    /* Main analysis code begins here
-     *
-     */
-    set_branch_statuses(this);
-
-    Long64_t nentries = fChain->GetEntriesFast();
-
-    Long64_t nbytes = 0, nb = 0;
-    for (Long64_t jentry=0; jentry<nentries;jentry++) {
-        Long64_t ientry = LoadTree(jentry);
-        if (ientry < 0) break;
-        nb = fChain->GetEntry(jentry);   nbytes += nb;
-
-        /* Main analysis code begins here
-         *
-         */
-        hc->lepton_count->Fill(nLepGood);
-        hc->top_quark_count->Fill(nGenTop);
-        for(int i=0; i<nLepGood; i++){
-            double delta_pt_val = LepGood_pt[i] - LepGood_mcPt[i];
-            hc->delta_pt->Fill(delta_pt_val);
-            hc->lepton_count_vs_delta_pt->Fill(nLepGood, delta_pt_val);
-        }
-        hc->jet_count->Fill(nJet);
-        int b_jet_count = 0;
-        for(int i=0; i<nJet; i++){
-            hc->b_jet_discriminator->Fill(Jet_btagCMVA[i]);
-            b_jet_count += (Jet_btagCMVA[i] > 0);
-        }
-        hc->b_jet_count->Fill(b_jet_count);
-
-        for(int i=0; i<nJet; i++){
-            hc->b_jet_pt_all->Fill(Jet_pt[i]);
-            if (b_jet_count >= 3){
-                hc->b_jet_pt_3_or_more->Fill(Jet_pt[i]);
-            }
-        }
-        int lepton_count_pass_miniiso= 0;
-        for(int i=0; i<nLepGood; i++){
-            lepton_count_pass_miniiso += mini_isolation_cut(i);
-        }
-        hc->lepton_count_pass_miniiso->Fill(lepton_count_pass_miniiso);
-
-        if (nLepGood == 3){
-            hc->jet_count_trilepton->Fill(nJet);
-        }
-        if (nLepGood == 2 && (LepGood_charge[0] * LepGood_charge[1]) == 1){
-            hc->jet_count_ss_dilepton->Fill(nJet);
-        }
-        if (nLepGood == 2 && (LepGood_charge[0] * LepGood_charge[1]) == -1){
-            hc->jet_count_os_dilepton->Fill(nJet);
-        }
-
-        if (nLepGood == 3){
-            vector<double> isos = {lepton_isolation(0),
-                                   lepton_isolation(1),
-                                   lepton_isolation(2)};
-            sort(isos.begin(), isos.end()); //Sorts ascending by default
-            hc->trilepton_iso_1->Fill(isos[2]);
-            hc->trilepton_iso_2->Fill(isos[1]);
-            hc->trilepton_iso_3->Fill(isos[0]);
-        }
-    }
-}
-#endif // tree_cxx

MiniTree.h

@@ -12,8 +12,6 @@
 #include <TChain.h>
 #include <TFile.h>
 
-#include "HistCollection.h"
-
 // Header file for the classes stored in the TTree if any.
 
 class MiniTree {
@@ -2183,14 +2181,12 @@ public :
    virtual Int_t    GetEntry(Long64_t entry);
    virtual Long64_t LoadTree(Long64_t entry);
    virtual void     Init(TTree *tree);
-   virtual void     Loop(HistCollection *hc);
+   virtual void     Loop();
    virtual Bool_t   Notify();
    virtual void     Show(Long64_t entry = -1);
 };
 
-#endif
 
-#ifdef minitree_cxx
 MiniTree::MiniTree(TTree *tree) : fChain(0) 
 {
 // if parameter tree is not specified (or zero), connect the file
@@ -3351,4 +3347,5 @@ Int_t MiniTree::Cut(Long64_t entry)
    return 1;
 }
 
-#endif // #ifdef minitree_cxx
+void MiniTree::Loop() { }
+#endif

MiniTreeDataSet.hpp

@@ -0,0 +1,56 @@
+#ifndef minitreedataset_h
+#define minitreedataset_h
+
+#include <string>
+
+#include "filval/filval.hpp"
+#include "filval_root/filval_root.hpp"
+#include "MiniTree.hpp"
+
+using namespace filval;
+using namespace filval::root;
+
+class MiniTreeDataSet : public filval::DataSet,
+                        public MiniTree{
+    private:
+        template <typename T>
+        void track_branch(const std::string bname, T *bval){
+            values[bname] = new filval::ObservedValue<T>(bval);
+            fChain->SetBranchStatus(bname.c_str(), true);
+        }
+        long next_entry;
+        long nentries;
+        bool load_next(){
+            if (next_entry >= nentries) return false;
+            fChain->GetEntry(next_entry);
+            ++next_entry;
+            return true;
+        }
+    public:
+        MiniTreeDataSet(TTree *tree)
+          :MiniTree(tree){
+            fChain->SetBranchStatus("*", false);
+
+            track_branch<int>("nLepGood", &nLepGood);
+            track_branch<int>("nJet", &nJet);
+
+            add_value(new FilterGreaterThan<int>(values["nLepGood"], 3), "nLepGoodCut");
+
+            add_container(new ContainerTH1I("nLepGood", "Lepton Multiplicity", 10, 0, 10, values["nLepGood"]));
+            add_container(new ContainerTH1I("nLepGood2", "Lepton Multiplicity", 10, 0, 10, values["nLepGood"]));
+            containers["nLepGood2"]->add_filter(values["nLepGoodCut"]);
+
+            add_value(new DerivedPair<int, int>(&values, "nLepGood", "nJet"), "nLepvsnJet");
+            add_container(new ContainerTGraph("nLepvsnJet", values["nLepvsnJet"]));
+
+
+            next_entry = 0;
+            nentries = fChain->GetEntriesFast();
+          }
+
+        void register_container(GenContainer* container){
+            containers[container->get_name()] = container;
+        }
+
+};
+#endif // minitreedataset_h

TTTT_Analysis.cpp

@@ -0,0 +1,94 @@
+#include <iostream>
+#include <vector>
+#include <utility>
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TCanvas.h"
+
+#include "filval/filval.hpp"
+#include "filval_root/filval_root.hpp"
+
+#include "MiniTreeDataSet.hpp"
+
+using namespace std;
+using namespace filval;
+using namespace filval::root;
+
+
+void print_pair(DerivedPair<double, double> dp){
+    pair<double, double> p = dp.get_value();
+    cout << "(" << p.first << ", " << p.second << ")\n";
+}
+
+void test1(){
+    double x = 12;
+    double y = 12;
+
+    ObservedValue<double> x_val(&x);
+    ObservedValue<double> y_val(&y);
+
+    DerivedPair<double, double> dp(&x_val, &y_val);
+    print_pair(dp);
+    x = 2;
+    y = 2;
+    print_pair(dp);
+
+    ContainerVector<double> cont(&x_val, "cont");
+    x = 12;
+    cont.fill();
+    x = 2;
+    cont.fill();
+    auto *container = cont.get_container();
+
+    for( auto v: *container )
+        cout << v << ", ";
+    cout << endl;
+}
+
+void test2(){
+    double x = 12;
+    ObservedValue<double> x_val(&x);
+    ContainerTH1D hist("h1", "Hist", 20, 0, 20, &x_val);
+    hist.fill();
+    hist.fill();
+    hist.fill();
+    x = 11;
+    hist.fill();
+    hist.fill();
+    hist.fill();
+    hist.fill();
+    hist.fill();
+    hist.fill();
+
+    TH1D* h = (TH1D*) hist.get_container();
+    TCanvas can("c1");
+    h->Draw();
+    can.Draw();
+    can.SaveAs("outfile.png");
+}
+
+void test3(){
+    TFile *f = TFile::Open("./data/TTTT_ext_treeProducerSusyMultilepton_tree.root");
+    TTree *tree = (TTree*) f->Get("tree");
+    MiniTreeDataSet mtds(tree);
+    mtds.process();
+    TH1* hist = ((ContainerTH1I*)mtds.get_container("nLepGood"))->get_container();
+    TCanvas can("c1");
+    hist->Draw();
+    can.Draw();
+    can.SaveAs("outfile.png");
+    can.Clear();
+    TGraph* graph= ((ContainerTGraph*)mtds.get_container("nLepvsnJet"))->get_container();
+    graph->Draw("A*");
+    can.Draw();
+    can.SaveAs("outfile3.png");
+
+    delete tree;
+    f->Close();
+    delete f;
+}
+
+int main(int argc, const char* argv[]){
+    test3();
+}

filval/README.md

@@ -0,0 +1,16 @@
+A FILter-VALue System
+=====================
+This is a header-only, generic, data analysis system that allows for creating performant generation of *Plots*. *Plots* contain *Values* and can make use of *Filters*.
+*Filters* can also depend of *Values*, and *Values* can depend on other *Values*. A *Dataset* is a generic object that contains a series of observations. The individual observations
+consist of a series of *Observed Values*. One can also define *Derived Values* which are calculated from *Observed Values* or other *Derived Values*. Care is taken automatically 
+so *Derived Values* are calculated at most once per observation.
+
+
+
+```C++
+MyDataSet myDataSet("somefile.root", "tree"); // MyDataSet subclasses DataSet
+TTreeValue<int> countmyDataSet;
+myDataSet.addValue
+
+Hist1D myplot(count, myDataSet, ); // Hist1D subclasses Plot
+```

filval/container.hpp

@@ -0,0 +1,72 @@
+#ifndef container_hpp
+#define container_hpp
+#include "value.hpp"
+#include "filter.hpp"
+#include <vector>
+
+namespace filval{
+class GenContainer{
+    private:
+        std::string name;
+        std::string desc;
+        std::vector<Filter*> filters;
+    protected:
+        virtual void _fill() = 0;
+    public:
+        GenContainer(const std::string name)
+          :name(name){ }
+        void add_filter(GenValue* filter){
+            filters.push_back(dynamic_cast<Filter*>(filter));
+        }
+        void fill(){
+            for (auto filter : filters){
+                if (!filter->get_value()){
+                    return;
+                }
+            }
+            _fill();
+        }
+        void set_description(const std::string& description){
+            desc = description;
+        }
+        const std::string& get_name(){
+            return name;
+        }
+};
+typedef std::map<std::string, GenContainer*> ContainerSet;
+
+template <typename H>
+class Container : public GenContainer{
+    protected:
+        H* container;
+    public:
+        Container(H* container, const std::string name)
+          :GenContainer(name),
+           container(container){ }
+        virtual H* get_container(){
+            return container;
+        }
+
+};
+
+template <typename T>
+class ContainerVector : public Container<std::vector<T> >{
+    private:
+        Value<T>* value;
+
+        void _fill(){
+            this->container->push_back(value->get_value());
+        }
+    public:
+        ContainerVector(std::vector<T> *container, Value<T>* value, const std::string name)
+          :Container<std::vector<T> >(container, name),
+           value(value){ }
+        ContainerVector(Value<T>* value, const std::string name)
+          :Container<std::vector<T> >(NULL, name),
+           value(value){
+            this->container = new std::vector<T>();
+        }
+};
+
+}
+#endif // container_hpp

filval/dataset.hpp

@@ -0,0 +1,36 @@
+#ifndef dataset_hpp
+#define dataset_hpp
+#include "value.hpp"
+#include "container.hpp"
+
+namespace filval{
+class DataSet{
+    protected:
+        ValueSet values;
+        ContainerSet containers;
+        virtual bool load_next() = 0;
+    public:
+        void process(){
+            while( load_next() ){
+                for(auto val : values)
+                    val.second->reset();
+                for(auto con : containers)
+                    con.second->fill();
+            }
+        }
+        void add_value(GenValue *value, const std::string& value_name ){
+            values[value_name] = value;
+        }
+        GenValue* get_value(std::string value_name){
+            return values[value_name];
+        }
+
+        void add_container(GenContainer *container){
+            containers[container->get_name()] = container;
+        }
+        GenContainer* get_container(std::string container_name){
+            return containers[container_name];
+        }
+};
+}
+#endif // dataset_hpp

filval/filter.hpp

@@ -0,0 +1,101 @@
+#ifndef filter_h
+#define filter_h
+#include <iostream>
+#include "value.hpp"
+/* A Filter is a special type of derived value that can only return a boolean.
+ * Container objects have a vector of filters that control if a "fill" call
+ * actually places data into the container or not.
+ */
+namespace filval {
+
+class Filter : public DerivedValue<bool>{ };
+
+
+class FilterAnd : public Filter {
+    protected:
+        Filter *filterA;
+        Filter *filterB;
+    void update_value(){
+        value = filterA->get_value() && filterB->get_value();
+    }
+    public:
+        FilterAnd(Filter *filterA, Filter *filterB)
+            :filterA(filterA), filterB(filterB){ }
+};
+
+class FilterOr : public Filter {
+    private:
+        Filter *filterA;
+        Filter *filterB;
+    void update_value(){
+        value = filterA->get_value() || filterB->get_value();
+    }
+    public:
+        FilterOr(Filter *filterA, Filter *filterB)
+            :filterA(filterA), filterB(filterB){ }
+};
+
+class FilterInv : public Filter {
+    private:
+        Filter *filter;
+    void update_value(){
+        value = !filter->get_value();
+    }
+    public:
+        FilterInv(Filter *filter)
+            :filter(filter){ }
+};
+
+template <typename T>
+class FilterGreaterThan : public Filter {
+    private:
+        Value<T> *filter_value;
+        Value<T> *range_low;
+        void update_value(){
+            value = filter_value->get_value() > range_low->get_value();
+        }
+    public:
+        FilterGreaterThan(GenValue* filter_value, GenValue* range_low)
+            :filter_value(dynamic_cast<Value<T>*>(filter_value)),
+             range_low(dynamic_cast<Value<T>*>(range_low)) { }
+
+        FilterGreaterThan(GenValue* filter_value, T range_low)
+            :filter_value(dynamic_cast<Value<T>*>(filter_value)){
+            this->range_low = new ConstantValue<T>(range_low);
+        }
+};
+
+template <typename T>
+class FilterLessThan : public Filter {
+    private:
+        Value<T> *filter_value;
+        Value<T> *range_high;
+        void update_value(){
+            value = filter_value->get_value() < range_high->get_value();
+        }
+    public:
+        FilterLessThan(GenValue* filter_value, GenValue* range_high)
+            :filter_value(dynamic_cast<Value<T>*>(filter_value)),
+             range_high(dynamic_cast<Value<T>*>(range_high)) { }
+
+        FilterLessThan(GenValue* filter_value, T range_high)
+            :filter_value(dynamic_cast<Value<T>*>(filter_value)){
+            this->range_high = new ConstantValue<T>(range_high);
+        }
+};
+
+
+template <typename T>
+class RangeFilter : public FilterAnd {
+    public:
+        RangeFilter(Value<T> *filter_value, Value<T> *range_low, Value<T> *range_high){
+            this->filterA = new FilterLessThan<T>(filter_value, range_high);
+            this->filterB = new FilterGreaterThan<T>(filter_value, range_low);
+         }
+        RangeFilter(Value<T> *filter_value, T range_low, T range_high){
+            this->filterA = new FilterLessThan<T>(filter_value, range_high);
+            this->filterB = new FilterGreaterThan<T>(filter_value, range_low);
+        }
+};
+}
+#endif // filter_h

filval/filval.hpp

@@ -0,0 +1,7 @@
+#ifndef filval_hpp
+#define filval_hpp
+#include "value.hpp"
+#include "filter.hpp"
+#include "container.hpp"
+#include "dataset.hpp"
+#endif // filval_hpp

filval/main.cpp

@@ -0,0 +1,45 @@
+#include <iostream>
+#include <vector>
+#include <utility>
+
+#include "value.hpp"
+
+using namespace std;
+
+
+void print_pair(DerivedPair<double, double> p){
+    pair<double, double> pr = p.get_value();
+    cout << "(" << pr.first << ", " << pr.second << ")\n";
+}
+
+int main(int argc, const char* argv[]){
+    /* Value<int, int> v; */
+    double x = 12;
+    double y = 12;
+
+    ObservedValue<double> x_val(&x);
+    ObservedValue<double> y_val(&y);
+
+    DerivedPair<double, double> dp(&x_val, &y_val);
+    print_pair(dp);
+    x = 2;
+    y = 2;
+    print_pair(dp);
+
+    ContainerVector<double> cont(&x_val);
+    x = 12;
+    cont.fill();
+    x = 2;
+    cont.fill();
+    auto *container = cont.get_container();
+
+    for( auto v: *container )
+        cout << v << ", ";
+    cout << endl;
+    /* ValA val; */
+    /* val.reset(); */
+    /* Address x = val.get_value(); */
+    /* cout << "Hello World " << x.street_address << "\n"; */
+    /* x = val.get_value(); */
+    /* cout << "Hello World " << x.street_address << "\n"; */
+}

filval/value.hpp

@@ -0,0 +1,115 @@
+#ifndef value_hpp
+#define value_hpp
+#include <iostream>
+#include <utility>
+#include <map>
+
+namespace filval{
+
+class GenValue{
+    public:
+        virtual void reset() = 0;
+};
+typedef std::map<std::string, GenValue*> ValueSet;
+
+/* class ValueSet: public std::map<std::string, GenValue*>{ */
+/*     public: */
+/*         GenValue*& operator[](const std::string& key){ */
+/*             GenValue*& value = (*this)[key]; */
+/*             if (value == NULL){ */
+/*                 std::cerr << "ERROR: key \""<<key */
+/*                     <<"\" not in valueset" << std::endl; */
+/*             } */
+/*             return value; */
+/*         } */
+/* }; */
+
+
+template <typename V>
+class Value : public GenValue{
+    public:
+        virtual V& get_value() = 0;
+};
+
+
+template <typename V>
+class ObservedValue : public Value<V>{
+    /* For "observed" values, there is nothing to calculate since this is
+     * merely a wrapper around a field in the observation. A pointer to the
+     * value is kept and it's value is read when requested.
+     */
+    private:
+        V *val_ref;
+    public:
+        ObservedValue(V *val_ref)
+            : val_ref(val_ref){}
+        V& get_value(){
+            return *val_ref;
+        }
+        void reset(){ }
+};
+
+
+template <typename V>
+class DerivedValue : public Value<V>{
+    /* A "derived" value is the result of some sort of calculation. Since it is
+     * desireable for the calculation to occur at most a single time for each
+     * observation, the result of the calculation is stored in the object. be
+     * sure that "reset" is called between processing observations to force a
+     * re-calculation.
+     */
+    protected:
+        V value; // The value to be calculated
+        bool value_valid;
+
+        virtual void update_value() = 0;
+    public:
+        V& get_value(){
+            if (!value_valid){
+                /* std::cout << "updating value!\n"; */
+                update_value();
+                value_valid = true;
+            }
+            return value;
+        }
+
+        void reset(){
+            value_valid = false;
+        }
+};
+
+
+template <typename T1, typename T2>
+class DerivedPair : public DerivedValue<std::pair<T1, T2> >{
+    protected:
+        std::pair<Value<T1>*, Value<T2>* > value_pair;
+        void update_value(){
+            this->value.first = value_pair.first->get_value();
+            this->value.second = value_pair.second->get_value();
+        }
+    public:
+        DerivedPair(ValueSet *values, const std::string &label1, const std::string &label2){
+            ValueSet &valueSet = *values;
+            value_pair.first  = (Value<T1>*) valueSet[label1];
+            value_pair.second = (Value<T2>*) valueSet[label2];
+         }
+        DerivedPair(Value<T1> *value1, Value<T2> *value2){
+            value_pair.first  = value1;
+            value_pair.second = value2;
+         }
+
+};
+
+template <typename V>
+class ConstantValue : public DerivedValue<V>{
+    protected:
+        V const_value;
+        void update_value(){
+            this->value = const_value;
+        }
+    public:
+        ConstantValue(V const_value)
+            :const_value(const_value) { }
+};
+}
+#endif // value_hpp

filval_root/README.md

@@ -0,0 +1,2 @@
+ROOT compatability layer for FILter-VALue System
+================================================

filval_root/container.hpp

@@ -0,0 +1,110 @@
+#ifndef root_container_hpp
+#define root_container_hpp
+#include <utility>
+#include "../filval/filval.hpp"
+#include "TH1.h"
+#include "TH2.h"
+#include "TGraph.h"
+#include <iostream>
+
+namespace filval::root {
+
+template <typename T>
+class ContainerTH1 : public Container<TH1>{
+    private:
+        void _fill(){
+            container->Fill(value->get_value());
+        }
+    protected:
+        Value<T> *value;
+        ContainerTH1(TH1* container, const std::string &name, Value<T> *value)
+          :Container<TH1>(container, name),
+           value(value){ }
+};
+
+class ContainerTH1D : public ContainerTH1<double>{
+    public:
+        ContainerTH1D(const std::string& name, const std::string& title,
+                      int nbins, double low, double high, GenValue* value)
+          :ContainerTH1<double>(NULL, name, dynamic_cast<Value<double>*>(value)){
+               this->container = new TH1D(name.c_str(), title.c_str(), nbins, low, high);
+        }
+};
+
+class ContainerTH1I : public ContainerTH1<int>{
+    public:
+        ContainerTH1I(const std::string& name, const std::string& title,
+                      int nbins, int low, int high, GenValue* value)
+          :ContainerTH1<int>(NULL, name, dynamic_cast<Value<int>*>(value)){
+               this->container = new TH1I(name.c_str(), title.c_str(), nbins, low, high);
+        }
+};
+
+
+template <typename T>
+class ContainerTH2 : public Container<TH2>{
+    private:
+        void _fill(){
+            std::pair<T, T> val = value->get_value();
+            container->Fill(val.first, val.second);
+        }
+    protected:
+        Value<std::pair<T, T> > *value;
+        ContainerTH2(TH2* container, const std::string &name, Value<std::pair<T, T> > *value)
+          :Container<TH2>(container, name),
+           value(value){ }
+};
+
+class ContainerTH2D : public ContainerTH2<double>{
+    public:
+        ContainerTH2D(const std::string& name, const std::string& title,
+                      int nbins_x, double low_x, double high_x,
+                      int nbins_y, double low_y, double high_y,
+                      GenValue* value)
+          :ContainerTH2<double>(NULL, name, dynamic_cast<Value<std::pair<double, double> >*>(value)){
+               this->container = new TH2D(name.c_str(), title.c_str(), nbins_x, low_x, high_x, nbins_y, low_y, high_y);
+        }
+};
+
+class ContainerTH2I : public ContainerTH2<int>{
+    public:
+        ContainerTH2I(const std::string& name, const std::string& title,
+                      int nbins_x, int low_x, int high_x,
+                      int nbins_y, int low_y, int high_y,
+                      GenValue* value)
+          :ContainerTH2<int>(NULL, name, dynamic_cast<Value<std::pair<int, int> >*>(value)){
+               this->container = new TH2I(name.c_str(), title.c_str(), nbins_x, low_x, high_x, nbins_y, low_y, high_y);
+        }
+};
+
+class ContainerTGraph : public Container<TGraph>{
+    private:
+        Value<std::pair<int, int> > *value;
+        std::vector<int> x_data;
+        std::vector<int> y_data;
+        bool data_modified;
+        void _fill(){
+            auto val = value->get_value();
+            x_data.push_back(val.first);
+            y_data.push_back(val.second);
+            /* std::cout << x_data.size() << std::endl; */
+            data_modified = true;
+        }
+    public:
+        ContainerTGraph(const std::string &name, GenValue* value)
+          :Container<TGraph>(new TGraph(), name),
+           value(dynamic_cast<Value<std::pair<int, int> >*>(value)),
+           data_modified(false){ }
+
+        TGraph* get_container(){
+            if (data_modified){
+                delete container;
+                container = new TGraph(x_data.size(), x_data.data(), y_data.data());
+                data_modified = false;
+            }
+            return container;
+        }
+};
+
+}
+#endif // root_container_hpp

filval_root/filval_root.hpp

@@ -0,0 +1,6 @@
+#ifndef filval_root_hpp
+#define filval_root_hpp
+#include "value.hpp"
+#include "container.hpp"
+/* #include "dataset.hpp" */
+#endif // filval_root_hpp

filval_root/value.hpp

@@ -0,0 +1,32 @@
+#ifndef root_values_hpp
+#define root_values_hpp
+#include "value.hpp"
+#include "TLorentzVector.h"
+
+namespace filval::root{
+
+class DerivedLorentzVector : public DerivedValue<TLorentzVector>{
+    protected:
+        Value<double> *pt;
+        Value<double> *eta;
+        Value<double> *phi;
+        Value<double> *m;
+        void update_value(){
+            value.SetPtEtaPhiM(pt->get_value(), eta->get_value(), phi->get_value(), m->get_value());
+        }
+    public:
+        DerivedLorentzVector(ValueSet *values,
+                             const std::string &pt_label,
+                             const std::string &eta_label,
+                             const std::string &phi_label,
+                             const std::string &m_label){
+            ValueSet &valueSet = *values;
+            pt = (Value<double>*) valueSet[pt_label];
+            eta = (Value<double>*) valueSet[eta_label];
+            phi = (Value<double>*) valueSet[phi_label];
+            m = (Value<double>*) valueSet[m_label];
+         }
+
+};
+}
+#endif // root_values_hpp

HistCollection.h

@@ -30,6 +30,7 @@ class HistCollection{
   private:
     map<std::string, HistWithPlottingInfo> histograms;
     std::string sample_name;
+    unsigned int event_count;
 
     std::string hist_id(const std::string &id){
         return sample_name+"_"+id;
@@ -56,6 +57,8 @@ class HistCollection{
     // List of included histogram objects
     TH1D *lepton_count;
     TH1D *lepton_count_pass_miniiso;
+    TH1D *lepton_count_pass_miniiso_event;
+    TH1D *lepton_count_pass_miniiso_ratio;
     TH1D *delta_pt;
     TH2D *lepton_count_vs_delta_pt;
     TH1D *top_quark_count;
@@ -73,9 +76,22 @@ class HistCollection{
     TH1D *trilepton_iso_1;
     TH1D *trilepton_iso_2;
     TH1D *trilepton_iso_3;
+    TH1D *trilepton_jet_iso_1;
+    TH1D *trilepton_jet_iso_2;
+    TH1D *trilepton_jet_iso_3;
+
+    TH1D *dijet_mass;
+    TH1D *dijet_mass_no_b;
+    TH1D *dijet_mass_one_b;
+    TH1D *dijet_mass_only_b;
+
+    TH1D *trijet_mass;
+    TH1D *trijet_mass_no_b;
+    TH1D *trijet_mass_one_b;
+    TH1D *trijet_mass_only_b;
 
     HistCollection(const std::string &sample_name)
-      : sample_name(sample_name)
+      : sample_name(sample_name), event_count(0)
     {
         // ---------------------------------------------
         lepton_count = book_histogram_1d("lepton_count", "Lepton Multiplicity",
@@ -85,6 +101,14 @@ class HistCollection{
                                                       "Number of Leptons passing mini-isolation",
                                                       8, 0, 8);
         // ---------------------------------------------
+        lepton_count_pass_miniiso_event = book_histogram_1d("lepton_count_pass_miniiso_event",
+                                                            "Lepton Multiplicity - Mini-iso Cut",
+                                                            8, 0, 8);
+        // ---------------------------------------------
+        lepton_count_pass_miniiso_ratio = book_histogram_1d("lepton_count_pass_miniiso_ratio",
+                                                            "Lepton Multiplicity - Mini-iso pass ratio",
+                                                            8, 0, 8);
+        // ---------------------------------------------
         delta_pt = book_histogram_1d("delta_pt", "{\\Delta P_T}",
                                      100, -50, 50);
         // ---------------------------------------------
@@ -121,13 +145,25 @@ class HistCollection{
                                         50, 0, 2000);
         // ---------------------------------------------
         trilepton_iso_1 = book_histogram_1d("trilepton_iso_1", "Trilepton Isolation - First",
-                                            50, 0, 2*PI);
+                                            25, 0, 2*PI);
         // ---------------------------------------------
         trilepton_iso_2 = book_histogram_1d("trilepton_iso_2", "Trilepton Isolation - Second",
-                                            50, 0, 2*PI);
+                                            25, 0, 2*PI);
         // ---------------------------------------------
         trilepton_iso_3 = book_histogram_1d("trilepton_iso_3", "Trilepton Isolation - Third",
-                                            50, 0, 2*PI);
+                                            25, 0, 2*PI);
+        // ---------------------------------------------
+        trilepton_jet_iso_1 = book_histogram_1d("trilepton_jet_iso_1", "Trilepton Jet Isolation - First",
+                                            25, 0, 2*PI);
+        // ---------------------------------------------
+        trilepton_jet_iso_2 = book_histogram_1d("trilepton_jet_iso_2", "Trilepton Jet Isolation - Second",
+                                            25, 0, 2*PI);
+        // ---------------------------------------------
+        trilepton_jet_iso_3 = book_histogram_1d("trilepton_jet_iso_3", "Trilepton Jet Isolation - Third",
+                                            25, 0, 2*PI);
+        // ---------------------------------------------
+        dijet_mass = book_histogram_1d("dijet_mass", "Di-Jet Mass",
+                                       100, 0, 1500);
         // ---------------------------------------------
     }
 
@@ -160,6 +196,10 @@ class HistCollection{
         return sample_name;
     }
 
+    unsigned int get_event_count(){
+        return event_count;
+    }
+
     vector<std::string> get_ids(){
         vector<std::string> vec;
         for(auto& p: histograms)
@@ -171,6 +211,17 @@ class HistCollection{
         TH1* hist = histograms[id].hist;
         return hist;
     }
+
+    HistCollection& operator++(){
+        event_count++;
+        return *this;
+    }
+
+    HistCollection operator++(int){
+        HistCollection tmp(*this);
+        operator++();
+        return tmp;
+    }
 };
 
 

legacy/MiniTree.cpp

@@ -0,0 +1,288 @@
+#ifndef minitree_cxx
+#define minitree_cxx
+
+#include <algorithm>
+#include <functional>
+#include <cstdio>
+#include <string>
+#include <vector>
+#include <map>
+
+#include <TStyle.h>
+#include <TLorentzVector.h>
+
+#include "MiniTree.h"
+#include "HistCollection.h"
+using namespace std;
+
+
+HistCollection* build_histograms(const char* sample_name, const char* filename){
+    HistCollection *hc = new HistCollection(sample_name);
+    TFile *f = TFile::Open(filename);
+    TTree *tree = (TTree*) f->Get("tree");
+    MiniTree minitree(tree);
+    minitree.Loop(hc);
+    delete tree;
+    f->Close();
+    delete f;
+    printf("Finished filling histograms for file %s with %d events.\n",
+            filename, hc->get_event_count());
+    return hc;
+}
+
+/* HistCollection* load_histograms(const char* sample_name){ */
+/*     string filename = string(sample_name) + "_histcollection.root"; */
+/*     TFile *f = TFile::Open(filename.c_str()); */
+/*     HistCollection *hc = new HistCollection(sample_name); */
+/*     TFile *f = TFile::Open(filename); */
+/*     TTree *tree = (TTree*) f->Get("tree"); */
+/*     MiniTree minitree(tree); */
+/*     minitree.Loop(hc); */
+/*     delete tree; */
+/*     f->Close(); */
+/*     delete f; */
+/*     printf("Finished filling histograms for file %s with %d events.\n", */
+/*             filename, hc->get_event_count()); */
+/*     return hc; */
+/* } */
+
+void set_branch_statuses(MiniTree *minitree){
+    auto on = [minitree](const char* bname){
+        minitree->fChain->SetBranchStatus(bname, true);
+    };
+    auto off = [minitree](const char* bname){
+        minitree->fChain->SetBranchStatus(bname, false);
+    };
+    off("*");
+    on("nLepGood");
+    on("LepGood_pt");
+    on("LepGood_eta");
+    on("LepGood_phi");
+    on("LepGood_mcPt");
+    on("LepGood_charge");
+    on("nJet");
+    on("Jet_btagCMVA");
+    on("Jet_pt");
+    on("Jet_eta");
+    on("Jet_phi");
+    on("nGenTop");
+    on("ngenLep");
+    on("genLep_sourceId");
+    on("genLep_pt");
+}
+
+double delta_r2(double eta1, double phi1, double eta2, double phi2){
+    return pow(eta1 - eta2, 2) + pow(phi1 - phi2, 2);
+}
+
+double isolation(int i, int n, function<double(int, int)> delta_r2_calculator){
+    double min_iso = -1;
+    for (int j=0; j<n; j++){
+        if (i == j) continue;
+        double curr_iso = delta_r2_calculator(i, j);
+        if ((curr_iso < min_iso && curr_iso >= 0) || min_iso == -1){
+            min_iso = curr_iso;
+        }
+    }
+    if (min_iso == -1)
+        return min_iso;
+    else
+        return sqrt(min_iso);
+}
+
+TLorentzVector g_vec1;
+TLorentzVector g_vec2;
+double diobject_mass(double pt1, double eta1, double phi1, double mass1,
+                     double pt2, double eta2, double phi2, double mass2){
+    g_vec1.SetPtEtaPhiM(pt1, eta1, phi1, mass1);
+    g_vec2.SetPtEtaPhiM(pt2, eta2, phi2, mass2);
+    g_vec1 += g_vec2;
+    return g_vec1.M();
+}
+
+void MiniTree::Loop(HistCollection* hc){
+
+    /* Define first a few little subroutines that will be used later in the
+     * main loop.
+     */
+    // Functions to calculate the \delta R^2 between two objects
+    auto lepton_lepton_delta_r2 = [this](int i, int j){
+        return delta_r2(this->LepGood_eta[i], this->LepGood_phi[i],
+                        this->LepGood_eta[j], this->LepGood_phi[j]);
+    };
+    auto lepton_jet_delta_r2 = [this](int i, int j){
+        return delta_r2(this->LepGood_eta[i], this->LepGood_phi[i],
+                        this->Jet_eta[j], this->Jet_phi[j]);
+    };
+    auto jet_jet_delta_r2 = [this](int i, int j){
+        return delta_r2(this->Jet_eta[i], this->Jet_phi[i],
+                        this->Jet_eta[j], this->Jet_phi[j]);
+    };
+
+    // Funtions to calculate isolation of objects
+    auto lepton_lepton_isolation = [this, lepton_lepton_delta_r2](int i){
+        if (i >= this->nLepGood) return -1.0;
+        return isolation(i, this->nLepGood, function<double(int, int)>(lepton_lepton_delta_r2));
+    };
+    auto lepton_jet_isolation = [this, lepton_jet_delta_r2](int i){
+        if (i >= this->nLepGood) return -1.0;
+        return isolation(i, this->nJet, function<double(int, int)>(lepton_jet_delta_r2));
+    };
+    auto mini_isolation_cut = [this, lepton_jet_isolation](int i){
+        /* Implement the mini-isolation cut which utilizes a dynamic isolation cone size
+         * based on lepton pt
+         */
+        double r0 = 0.5;
+        double pt0 = 7.5; // GeV
+
+        double iso_cut = min(r0, pt0/this->LepGood_pt[i]);
+        double isolation = lepton_jet_isolation(i);
+        return isolation > iso_cut;
+    };
+    auto jet_jet_isolation = [this, jet_jet_delta_r2](int i){
+        if (i >= this->nJet) return -1.0;
+        return isolation(i, this->nJet, function<double(int, int)>(jet_jet_delta_r2));
+    };
+    auto dijet_mass = [this](int i, int j){
+        if (i >= this->nJet || j >= this->nJet) return -1.0;
+        return diobject_mass(this->Jet_pt[i], this->Jet_eta[i],
+                             this->Jet_phi[i], this->Jet_mass[i],
+                             this->Jet_pt[j], this->Jet_eta[j],
+                             this->Jet_phi[j], this->Jet_mass[j]);
+    };
+    auto dilepton_mass = [this](int i, int j){
+        if (i >= this->nLepGood || j >= this->nLepGood) return -1.0;
+        return diobject_mass(this->LepGood_pt[i], this->LepGood_eta[i],
+                             this->LepGood_phi[i], this->LepGood_mass[i],
+                             this->LepGood_pt[j], this->LepGood_eta[j],
+                             this->LepGood_phi[j], this->LepGood_mass[j]);
+    };
+    auto dijet_mass_set = [this, dijet_mass](){
+        auto *mass_set = new vector<double>();
+        for(int i=0; i<(this->nJet-1); i++){
+            for(int j=i+1; j<this->nJet; j++){
+                mass_set->push_back(dijet_mass(i, j));
+            }
+        }
+        return mass_set;
+    };
+
+    auto nBJets = [this](float mva_cut = 0){
+        int n = 0;
+        for(int i=0; i<this->nJet; ++i)
+            if(this->Jet_btagCMVA[i] > mva_cut) ++n;
+        return n;
+    };
+
+    auto lepton_mass_window_pass = [this, dilepton_mass](double low, double high){
+        for (int i=0; i<this->nLepGood, ++i){
+            for (int j=0; j<this->nLepGood, ++j){
+                double mass = dilepton_mass(i, j);
+                if (low < mass && mass < high){
+                    return false;
+                }
+            }
+        }
+        return true;
+    }
+    auto ss_dilepton_base_selection[this, nBJets, lepton_mass_window_pass](){
+        if (this->nLepGood != 2) return false;
+        if (this->LepGood_charge[0] != this->LepGood_charge[1]) return false;
+        if (this->nJet < 4) return false;
+        if (nBJets() < 3) return false;
+        if (!lepton_mass_window_pass(70, 105)) return false;
+        return true;
+    }
+
+    auto trilepton_base_selection[this](){
+        if (nLepGood != 3) return false;
+        if (nJet < 4) return false;
+        if (nBJets() < 3) return false;
+        if (!lepton_mass_window_pass(70, 105)) return false;
+        return true;
+    }
+
+    if (fChain == 0) return;
+    /* Set branch statuses so only used branches are loaded. This
+     * can greatly reduce execution time.
+     */
+    set_branch_statuses(this);
+
+    Long64_t nentries = fChain->GetEntriesFast();
+
+    Long64_t nbytes = 0, nb = 0;
+    for (Long64_t jentry=0; jentry<nentries;jentry++) {
+        Long64_t ientry = LoadTree(jentry);
+        if (ientry < 0) break;
+        nb = fChain->GetEntry(jentry);   nbytes += nb;
+
+        /* Main analysis code begins here
+         *
+         */
+        ++(*hc); // Increments event count
+        hc->lepton_count->Fill(nLepGood);
+        hc->top_quark_count->Fill(nGenTop);
+        for(int i=0; i<nLepGood; i++){
+            double delta_pt_val = LepGood_pt[i] - LepGood_mcPt[i];
+            hc->delta_pt->Fill(delta_pt_val);
+            hc->lepton_count_vs_delta_pt->Fill(nLepGood, delta_pt_val);
+        }
+        hc->jet_count->Fill(nJet);
+        int b_jet_count = 0;
+        for(int i=0; i<nJet; i++){
+            hc->b_jet_discriminator->Fill(Jet_btagCMVA[i]);
+            b_jet_count += (Jet_btagCMVA[i] > 0);
+        }
+        hc->b_jet_count->Fill(b_jet_count);
+
+        for(int i=0; i<nJet; i++){
+            hc->b_jet_pt_all->Fill(Jet_pt[i]);
+            if (b_jet_count >= 3){
+                hc->b_jet_pt_3_or_more->Fill(Jet_pt[i]);
+            }
+        }
+        int lepton_count_pass_miniiso= 0;
+        for(int i=0; i<nLepGood; i++){
+            lepton_count_pass_miniiso += mini_isolation_cut(i);
+        }
+        hc->lepton_count_pass_miniiso->Fill(lepton_count_pass_miniiso);
+        if (lepton_count_pass_miniiso == nLepGood){
+            hc->lepton_count_pass_miniiso_event->Fill(nLepGood);
+        }
+
+        if (nLepGood == 3){
+            hc->jet_count_trilepton->Fill(nJet);
+        }
+        if (nLepGood == 2 && (LepGood_charge[0] * LepGood_charge[1]) == 1){
+            hc->jet_count_ss_dilepton->Fill(nJet);
+        }
+        if (nLepGood == 2 && (LepGood_charge[0] * LepGood_charge[1]) == -1){
+            hc->jet_count_os_dilepton->Fill(nJet);
+        }
+
+        if (nLepGood == 3){
+            vector<double> isos = {lepton_jet_isolation(0),
+                                   lepton_jet_isolation(1),
+                                   lepton_jet_isolation(2)};
+            sort(isos.begin(), isos.end()); //Sorts ascending by default
+            hc->trilepton_iso_1->Fill(isos[2]);
+            hc->trilepton_iso_2->Fill(isos[1]);
+            hc->trilepton_iso_3->Fill(isos[0]);
+            vector<double> jet_isos = {jet_jet_isolation(0),
+                                       jet_jet_isolation(1),
+                                       jet_jet_isolation(2)};
+            sort(jet_isos.begin(), jet_isos.end());
+            hc->trilepton_jet_iso_1->Fill(jet_isos[2]);
+            hc->trilepton_jet_iso_2->Fill(jet_isos[1]);
+            hc->trilepton_jet_iso_3->Fill(jet_isos[0]);
+        }
+
+        vector<double> *mass_set = dijet_mass_set();
+        for( double &mass : *mass_set ){
+            hc->dijet_mass->Fill(mass);
+        }
+    }
+    hc->lepton_count_pass_miniiso_ratio->Divide(hc->lepton_count_pass_miniiso_event,
+                                                hc->lepton_count);
+}
+#endif // tree_cxx