begins re-implementing old ipynb analysis/plots using new framework

.gitignore

@@ -2,6 +2,7 @@
 .ipynb_checkpoints/
 __pycache__/
 build/
+data/
 
 # C++ Build Files
 *.d

analysis/MiniTree.hpp

@@ -2176,7 +2176,7 @@ public :
    TBranch        *b_PhoGood_chHadIsoRC;   //!
    TBranch        *b_PhoGood_drMinParton;   //!
 
-   MiniTree(TTree *tree=0);
+   MiniTree();
    virtual ~MiniTree();
    virtual Int_t    Cut(Long64_t entry);
    virtual Int_t    GetEntry(Long64_t entry);
@@ -2188,26 +2188,9 @@ public :
 };
 
 
-MiniTree::MiniTree(TTree *tree) : fChain(0) 
-{
-// if parameter tree is not specified (or zero), connect the file
-// used to generate this class and read the Tree.
-   if (tree == 0) {
-      TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("TTTT_ext_treeProducerSusyMultilepton_tree.root");
-      if (!f || !f->IsOpen()) {
-         f = new TFile("TTTT_ext_treeProducerSusyMultilepton_tree.root");
-      }
-      f->GetObject("tree",tree);
-
-   }
-   Init(tree);
-}
+MiniTree::MiniTree() : fChain(0) { }
 
-MiniTree::~MiniTree()
-{
-   if (!fChain) return;
-   delete fChain->GetCurrentFile();
-}
+MiniTree::~MiniTree() { }
 
 Int_t MiniTree::GetEntry(Long64_t entry)
 {

analysis/MiniTreeDataSet.hpp

@@ -15,6 +15,10 @@ using namespace fv::root;
 class MiniTreeDataSet : public DataSet,
                         public MiniTree{
     private:
+        std::string input_filename;
+        std::string output_filename;
+        TFile* input_file;
+        TFile* output_file;
         long next_entry;
         long nentries;
         bool load_next(){
@@ -31,12 +35,21 @@ class MiniTreeDataSet : public DataSet,
         }
 
     public:
-        MiniTreeDataSet(TTree *tree)
-          :MiniTree(tree){
-            next_entry = 0;
+        MiniTreeDataSet(const std::string& input_filename, const std::string output_filename)
+          :input_filename(input_filename),
+           output_filename(output_filename),
+           next_entry(0) {
+            input_file = TFile::Open(input_filename.c_str());
+            Init((TTree*) input_file->Get("tree"));
             nentries = fChain->GetEntriesFast();
+            output_file = TFile::Open(output_filename.c_str(), "UPDATE");
           }
 
+        ~MiniTreeDataSet(){
+            input_file->Close();
+            output_file->Close();
+        }
+
         template <typename T>
         Value<T>* track_branch(const std::string& bname){
             TBranch* branch = fChain->GetBranch(bname.c_str());
@@ -65,6 +78,7 @@ class MiniTreeDataSet : public DataSet,
             return new PointerValue<T>(bname, bref);
         }
         void save_all(){
+            output_file->cd();
             for(auto container : containers){
                 container.second->save_as("outfile", SaveOption::ROOT);
             }

analysis/TTTT_Analysis.cpp

@@ -94,6 +94,8 @@ void declare_values(MiniTreeDataSet& mt){
     new WrapperVector<float>("nLepGood", "LepGood_phi", "LepGood_phi");
     new WrapperVector<float>("nLepGood", "LepGood_mass", "LepGood_mass");
 
+    new WrapperVector<float>("nJet", "Jet_btagCMVA", "Jet_btagCMVA");
+
 
     new ZipMapFour<float, float>(get_energy, "LepGood_pt", "LepGood_eta", "LepGood_phi", "LepGood_mass",
                                  "lepton_energy");
@@ -106,12 +108,22 @@ void declare_values(MiniTreeDataSet& mt){
     new Mean<float>("lepton_energy", "lepton_energy_mean");
 
 
-    new Filter("nLepGood>=3", FUNC(([nLepGood=lookup("nLepGood")](){
-            return dynamic_cast<Value<int>*>(nLepGood)->get_value() >=3;})));
+    new Count<float>(GenFunction::register_function<bool(float)>("bJet_Selection", FUNC(([](float x){return x>0;}))),
+                     "Jet_btagCMVA",  "b_jet_count");
+
+    new Filter("trilepton", FUNC(([nLepGood=lookup("nLepGood")](){
+            return dynamic_cast<Value<int>*>(nLepGood)->get_value() ==3;})));
+
+    new Filter("os-dilepton", FUNC(([nLepGood=lookup("nLepGood"), LepGood_charge=lookup("LepGood_charge")](){
+                    bool is_dilepton = dynamic_cast<Value<int>*>(nLepGood)->get_value()==2;
+                    int* charge = dynamic_cast<Value<int*>*>(LepGood_charge)->get_value();
+                    return is_dilepton && (charge[0] != charge[1]);})));
+
+    new Filter("ss-dilepton", FUNC(([nLepGood=lookup("nLepGood"), LepGood_charge=lookup("LepGood_charge")](){
+                    bool is_dilepton = dynamic_cast<Value<int>*>(nLepGood)->get_value()==2;
+                    int* charge = dynamic_cast<Value<int*>*>(LepGood_charge)->get_value();
+                    return is_dilepton && (charge[0] == charge[1]);})));
 
-    new Filter("nLepGood<=4", FUNC(([nLepGood=lookup("nLepGood")](){
-            return dynamic_cast<Value<int>*>(nLepGood)->get_value() <=4;})));
-    new RangeFilter<int>("3<=nLepGood<5", dynamic_cast<Value<int>*>(lookup("nLepGood")), 3, 5);
 }
 
 void declare_containers(MiniTreeDataSet& mt){
@@ -123,38 +135,50 @@ void declare_containers(MiniTreeDataSet& mt){
     mt.register_container(new ContainerTH1F("lepton_energy_min", "Lepton Energy - Min", lookup("lepton_energy_min"), 50, 0, 500));
     mt.register_container(new ContainerTH1F("lepton_energy_rng", "Lepton Energy - Range", lookup("lepton_energy_range"), 50, 0, 500));
 
-    mt.register_container(new ContainerTH1I("nLepGood2", "Lepton Multiplicity", lookup("nLepGood"), 10, 0, 10));
-    mt.get_container("nLepGood2")->add_filter(lookup_filter("3<=nLepGood<5"));
-
-    mt.register_container(new ContainerTH1I("nLepGood3", "Lepton Multiplicity", lookup("nLepGood"), 10, 0, 10));
-    mt.get_container("nLepGood3")->add_filter(!(*lookup_filter("3<=nLepGood<5")));
 
     mt.register_container(new ContainerTGraph("nLepvsnJet", new Pair<int, int>("nLepGood", "nJet") ));
 
     mt.register_container(new ContainerTH2FMany("lepton_energy_vs_pt", "Lepton Energy - Range", lookup("lepton_energy_lepton_pt"),
                                                50, 0, 500, 50, 0, 500));
+
+    mt.register_container(new ContainerTH1I("b_jet_count", "B-Jet Multiplicity", lookup("b_jet_count"), 10, 0, 10));
+
+
+    mt.register_container(new ContainerTH1I("jet_count_os_dilepton", "Jet Multiplicity - OS Dilepton Events", lookup("nJet"), 14, 0, 14));
+    mt.get_container("jet_count_os_dilepton")->add_filter(lookup_filter("os-dilepton"));
+    mt.register_container(new ContainerTH1I("jet_count_ss_dilepton", "Jet Multiplicity - SS Dilepton Events", lookup("nJet"), 14, 0, 14));
+    mt.get_container("jet_count_ss_dilepton")->add_filter(lookup_filter("ss-dilepton"));
+    mt.register_container(new ContainerTH1I("jet_count_trilepton", "Jet Multiplicity - Trilepton Events", lookup("nJet"), 14, 0, 14));
+    mt.get_container("jet_count_trilepton")->add_filter(lookup_filter("trilepton"));
+}
+
+std::string replace_suffix(const std::string& input, const std::string& new_suffix){
+    return input.substr(0, input.find_last_of(".")) + new_suffix;
 }
 
-void run_analysis(const std::string& filename){
-    TFile *f = TFile::Open(filename.c_str());
-    TTree *tree = (TTree*) f->Get("tree");
-    MiniTreeDataSet mt(tree);
+void run_analysis(const std::string& input_filename, bool silent){
+    string log_filename = replace_suffix(input_filename, "_result.log");
+    Log::init_logger(log_filename, LogPriority::kLogDebug);
+
+    string output_filename = replace_suffix(input_filename, "_result.root");
+    MiniTreeDataSet mt(input_filename, output_filename);
+
     enable_branches(mt);
     declare_values(mt);
     declare_containers(mt);
-    mt.process();
+
+    mt.process(silent);
     mt.save_all();
 }
 
 int main(int argc, char * argv[])
 {
-    Log::init_logger("log.txt", LogPriority::kLogDebug);
     ArgParser args(argc, argv);
-    if(args.cmdOptionExists("-f")) {
-        run_analysis(args.getCmdOption("-f"));
-    }
-    else {
+    if(!args.cmdOptionExists("-f")) {
         cout << "Usage: ./main -f input_minitree.root" << endl;
+        return -1;
     }
-    return 0;
+    bool silent = args.cmdOptionExists("-s");
+    string input_filename = args.getCmdOption("-f");
+    run_analysis(input_filename, silent);
 }

filval/dataset.hpp

@@ -20,20 +20,20 @@ class DataSet{
         virtual int get_current_event() = 0;
 
     public:
-        void process(){
+        void process(bool silent=false){
             int events, current_event;
             summary();
             events = get_events();
-            std::cout << std::endl;
+            if (!silent) std::cout << std::endl;
             while( load_next() ){
                 current_event = get_current_event();
-                std::cout << "\rprocessing event: " << current_event+1 << "/" << events << std::flush;
+                if (!silent) std::cout << "\rprocessing event: " << current_event+1 << "/" << events << std::flush;
                 GenValue::reset();
                 for(auto con : containers){
                     con.second->fill();
                 }
             }
-            std::cout << " Finished!" << std::endl;
+            if (!silent) std::cout << " Finished!" << std::endl;
         }
 
         virtual void save_all(){

filval/value.hpp

@@ -486,6 +486,30 @@ class ZipMapFour : public DerivedValue<std::vector<R> >{
                       alias){ }
 };
 
+/**
+ *
+ */
+template<typename T>
+class Count : public DerivedValue<int>{
+    private:
+        Function<bool(T)>& selector;
+        Value<std::vector<T> >* v;
+        void update_value(){
+            value = 0;
+            for(auto val : v->get_value()){
+                if(selector(val))
+                    value++;
+            }
+        }
+    public:
+        Count(Function<bool(T)>& selector, Value<std::vector<T>>* v, const std::string alias="")
+          :DerivedValue<int>("count("+selector.get_name()+":"+v->get_name()+")", alias),
+           selector(selector), v(v) { }
+
+        Count(Function<bool(T)>& selector, const std::string& v_name, const std::string alias="")
+          :Count(selector, dynamic_cast<Value<std::vector<T> >*>(GenValue::get_value(v_name)), alias) { }
+};
+
 
 /**
  * Reduce a Value of type vector<T> to just a T.

filval_root/container.hpp

@@ -19,11 +19,11 @@ void _save_img(TObject* container, const std::string& fname){
     delete c1;
 }
 
-void _save_bin(TObject* container, const std::string& fname){
-    INFO("Saving object: " << container->GetName() << " into file " << fname);
-    TFile* f = TFile::Open(fname.c_str(), "UPDATE");
+void _save_bin(TObject* container){
+    INFO("Saving object: " << container->GetName() << " into file " << gDirectory->GetName());
+    /* TFile* f = TFile::Open(fname.c_str(), "UPDATE"); */
     container->Write(container->GetName(), TObject::kOverwrite);
-    f->Close();
+    /* f->Close(); */
 }
 
 void _save_as(TObject* container, const std::string& fname, const SaveOption& option = SaveOption::PNG) {
@@ -33,7 +33,7 @@ void _save_as(TObject* container, const std::string& fname, const SaveOption& op
         case PDF:
             _save_img(container, fname+".pdf"); break;
         case ROOT:
-            _save_bin(container, fname+".root"); break;
+            _save_bin(container); break;
         default:
             break;
     }

python/utils.py

@@ -1,10 +1,47 @@
 import io
 import sys
 from math import ceil, sqrt
+from subprocess import run
 import itertools as it
 import ROOT
 
 
+class HistCollection:
+    def __init__(self, sample_name, input_filename,
+                 exe_path="../build/main",
+                 rebuild_hists = False):
+        self.sample_name = sample_name
+        if rebuild_hists:
+            run([exe_path, "-s", "-f", input_filename])
+        output_filename = input_filename.replace(".root", "_result.root")
+        self._file = ROOT.TFile.Open(output_filename)
+        l = self._file.GetListOfKeys()
+        self.map = {}
+        for i in range(l.GetSize()):
+            name = l.At(i).GetName()
+            self.map[name] = self._file.Get(name)
+            setattr(self, name, self.map[name])
+
+    def draw(self, canvas, shape=None):
+        if shape is None:
+            n = int(ceil(sqrt(len(self.map))))
+            shape = (n, n)
+        print(shape)
+        canvas.Clear()
+        canvas.Divide(*shape)
+        for i, hist in enumerate(self.map.values()):
+            canvas.cd(i+1)
+            try:
+                hist.SetStats(False)
+            except AttributeError:
+                pass
+            print(i, hist, str(type(hist)))
+            draw_option = ""
+            if (type(hist) == ROOT.TH2F):
+                draw_option = "COLZ"
+            hist.Draw(draw_option)
+
+
 class OutputCapture:
     def __init__(self):
         self.my_stdout = io.StringIO()
@@ -93,7 +130,7 @@ def stack_hist(hists,
 def stack_hist_array(canvas, histcollections, fields, titles,
                      shape=None, **kwargs):
     def get_hist_set(attrname):
-        hists, labels = zip(*[(getattr(h, attrname), h.get_sample_name())
+        hists, labels = zip(*[(getattr(h, attrname), h.sample_name)
                               for h in histcollections])
         return hists, labels
     n_fields = len(fields)