initial commit of rudimentary analysis system for TTTT

HistCollection.h

@@ -0,0 +1,73 @@
+#ifndef histcollection_h
+#define histcollection_h
+
+#include <string>
+#include <sstream>
+#include <map>
+
+#include <TH1.h>
+#include <TH2.h>
+
+
+
+class HistCollection{
+  private:
+    map<std::string, TH1*> histograms;
+    std::string prefix;
+    
+    TH1D* book_histogram_1d(const std::string &id, const std::string &title,
+                            int nbins, double min, double max){
+        const char* full_id = (prefix+id).c_str();
+        TH1D *hist = new TH1D(full_id, title.c_str(), nbins, min, max);
+        histograms[id] = hist;
+        return hist;
+    };
+
+    TH2D* book_histogram_2d(const std::string id, const std::string title,
+                            int nbins_x, double min_x, double max_x,
+                            int nbins_y, double min_y, double max_y){
+        const char* full_id = (prefix+id).c_str();
+        TH2D *hist = new TH2D(full_id, title.c_str(), nbins_x, min_x, max_x, nbins_y, min_y, max_y);
+        histograms[id] = hist;
+        return hist;
+    };
+  public:
+    /* List of included histogram objects
+    */
+    TH1D *lepton_count;
+    TH1D *delta_pt;
+    TH2D *lepton_count_vs_delta_pt;
+    
+    
+    HistCollection(const std::string &filename)
+      : prefix(filename.substr(0, filename.size()-5))
+    {
+//        this->prefix = prefix.substr(0, prefix.size()-5);
+        lepton_count = book_histogram_1d("lepton_count", "Good Lepton Count", 10, 0, 10);
+        delta_pt = book_histogram_1d("delta_pt", "{\\\\Delta P_T}", 100, -50, 50);
+        lepton_count_vs_delta_pt = book_histogram_2d("lepton_count_vs_delta_pt", "Good Lepton Count Vs {\\\\delta P_T}",
+                                                     10, 0, 10,
+                                                     100, -50, 50);
+    };
+    
+    HistCollection(): HistCollection("") {};
+    
+    ~HistCollection(){
+        for(auto& p: histograms)
+            delete p.second;
+    };
+    
+    vector<std::string> get_ids(){
+        vector<std::string> vec;
+        for(auto& p: histograms)
+            vec.push_back(p.first);
+        return vec;
+    }
+    TH1* get_by_id(const std::string &id){
+        TH1* hist = histograms[id];
+        return hist;
+    }
+};
+
+
+#endif // histcollection_h

MiniTree.cpp

@@ -0,0 +1,49 @@
+#ifndef minitree_cxx
+#define minitree_cxx
+#include <cstdio>
+#include <string>
+#include <map>
+
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+
+#include "MiniTree.h"
+#include "HistCollection.h"
+
+
+HistCollection* build_histograms(const char* filename){
+    HistCollection *hc = new HistCollection(filename);
+    TFile *f = TFile::Open(filename);
+    TTree *tree = (TTree*) f->Get("tree");
+    MiniTree minitree(tree);
+    minitree.Loop(hc);
+    delete tree;
+    f->Close();
+    delete f;
+    std::printf("Finished filling histograms for file %s\n", filename);
+    return hc;
+}
+
+
+void MiniTree::Loop(HistCollection* hc){
+    if (fChain == 0) return;
+
+    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;
+        // if (Cut(ientry) < 0) continue;
+        hc->lepton_count->Fill(nLepGood);
+        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);
+        }
+    }
+}
+#endif // tree_cxx

plotter.py

@@ -0,0 +1,7 @@
+import ROOT
+
+
+canvas = ROOT.TCanvas("c1")
+def do_plots():
+    ROOT.lepton_count.Draw()
+    canvas.Draw()

utils.py

@@ -0,0 +1,48 @@
+import io
+import sys
+import ROOT
+
+class OutputCapture:
+    def __init__(self):
+        self.my_stdout = io.StringIO()
+        self.my_stderr = io.StringIO()
+
+    def get_stdout(self):
+        self.my_stdout.seek(0)
+        return self.my_stdout.read()
+
+    def get_stderr(self):
+        self.my_stderr.seek(0)
+        return self.my_stderr.read()
+
+    def __enter__(self):
+        self.stdout = sys.stdout
+        self.stderr = sys.stderr
+        sys.stdout = self.my_stdout
+        sys.stderr = self.my_stderr
+
+    def __exit__(self, *args):
+        sys.stdout = self.stdout
+        sys.stderr = self.stderr
+        self.stdout = None
+        self.stderr = None
+
+def bin_range(n, end=None):
+    if end is None:
+        return range(1, n+1)
+    else:
+        return range(n+1, end+1)
+        
+def normalize_columns(hist2d):
+    normHist = ROOT.TH2D(hist2d)
+    cols, rows = hist2d.GetNbinsX(), hist2d.GetNbinsY()
+    for col in bin_range(cols):
+        sum_ = 0;
+        for row in bin_range(rows):
+            sum_ += hist2d.GetBinContent(col, row)
+        if sum_ == 0:
+            continue
+        for row in bin_range(rows):
+            norm = hist2d.GetBinContent(col, row) / sum_
+            normHist.SetBinContent(col, row, norm)
+    return normHist