Another TSelector analyzer

apz/copytree.C

@@ -0,0 +1,27 @@
+void copytree() {
+
+  //Get old file, old tree and set top branch address
+  //TFile *oldfile = new TFile("APZ-tree-Data-8TeV-Oct27.root");
+  TFile *oldfile = new TFile("APZ-tree-Data-MuEG-8TeV-2018July25.root");
+  TTree *oldtree = (TTree*)oldfile->Get("apzTree/apzTree");
+  //TLorentzVector *lep1   = new TLorentzVector();
+  //oldtree->SetBranchAddress("lep1",&lep1);
+  oldtree->SetBranchStatus("*",0);
+  oldtree->SetBranchStatus("lep1",1);
+  oldtree->SetBranchStatus("lep2",1);
+  oldtree->SetBranchStatus("gamma",1);
+  oldtree->SetBranchStatus("jet1",1);
+  oldtree->SetBranchStatus("jet2",1);
+  oldtree->SetBranchStatus("met",1);
+  oldtree->SetBranchStatus("run",1);
+  oldtree->SetBranchStatus("event",1);
+
+  //Create a new file + a clone of old tree in new file
+  TFile *newfile = new TFile("xsmall.root","recreate");
+  TTree *newtree = oldtree->CloneTree();
+
+  newtree->Print();
+  newfile->Write();
+  delete oldfile;
+  delete newfile;
+}

apz/plotit.py

@@ -0,0 +1,102 @@
+#! /usr/bin/env python
+
+import sys,os
+from math import sqrt
+from ROOT import *
+
+
+gROOT.SetBatch()
+gStyle.SetOptStat(0)
+
+def handleOverflowBins(hist):
+  if hist == None:
+    return
+  nBins   = hist.GetNbinsX()
+  lastBin = hist.GetBinContent(nBins)
+  ovflBin = hist.GetBinContent(nBins+1);
+  lastBinErr = hist.GetBinError(nBins);
+  ovflBinErr = hist.GetBinError(nBins+1);
+  firstBin    = hist.GetBinContent(1);
+  undflBin    = hist.GetBinContent(0);
+  firstBinErr = hist.GetBinError(1);
+  undflBinErr = hist.GetBinError(0);
+  hist.SetBinError(nBins, sqrt(pow(lastBinErr,2) + pow(ovflBinErr,2)) );
+  hist.SetBinContent(1, firstBin+undflBin);
+  hist.SetBinError(1, sqrt(pow(firstBinErr,2) + pow(undflBinErr,2)) );
+  hist.SetBinContent(0,0);
+  hist.SetBinContent(nBins+1,0);
+
+if __name__ == "__main__":
+  print "This is the __main__ part"
+
+  f = TFile("my_higgsHistograms.root",'OPEN')
+
+  ks = f.Get("AfterCuts")
+  names = []
+  for key in ks.GetListOfKeys():
+    names.append(key.GetName())
+
+  print names
+
+  h = {}
+
+
+  primary = "InAPZ"
+
+  for n in ["lep1_pt", "lep2_pt", "gamma_pt",
+            "lep1_phi","lep2_phi","gamma_phi",
+            "lep1_eta","lep2_eta","gamma_eta",
+            "dR_gamma_jet1","dR_gamma_jet2","dR_lep", "dRjj","dR_gamma_dilep",
+            "Mllg_full", "Mllg_higg", "Mllg_zed", "Mll_APZ_JPsi","Mll_18",
+            "pt_ratio_dilep_gamma",
+            "Mgj_full", "Mgj_higg",
+            "Mjj_full","Mjj_zoom","Mllgj_full","Mllgj_higg", "Mllgj_x200",
+            "jet1_pt","jet1_eta","jet1_phi",
+            "jet2_pt","jet2_eta","jet2_phi",
+            "met_Et","met_phi","met_dPhi_dilep","met_dPhi_gamma"]:
+
+
+    print "Plotting", n
+    m = 0
+    for d in [primary, 'AfterCuts', 'InJpsi']:
+      h[d] = f.Get(d+'/'+n+'_'+d)
+      if h[d] == None:
+        print "Not a good histo:", h[d]
+        continue
+
+      if d==primary:
+        opt = "hist"
+        h[d].SetLineColor(kGreen+2)
+      else:
+        opt = "hist same"
+        if d=='InJpsi':
+          h[d].SetLineColor(kRed+1)
+        if d=='AfterCuts':
+          h[d].SetLineColor(kBlue+1)
+        if d=='M18':
+          h[d].SetLineColor(kOrange+1)
+
+
+      norm = h[primary].Integral()
+
+      if h[d].Integral()!=0:
+        h[d].Scale(norm/h[d].Integral())
+
+      if n in ["Mllg_full","gamma_pt","Mjj_full","Mllgj","Mgj_full","pt_ratio_dilep_gamma","met"]:
+        handleOverflowBins(h[d])
+
+      h[d].Draw(opt)
+      m = max(h[d].GetMaximum(), m)
+      h[primary].SetMaximum(1.1*m)
+
+    c1.SaveAs("figs/fig_"+n+".png")
+
+
+  ## < indent of for loop
+  evt = f.Get("evt")
+  SplusB = evt.GetBinContent(3)
+  B = (9./11)*evt.GetBinContent(4)
+  significance = 2*(sqrt(SplusB) - sqrt(B))
+  print "S+B=%i, S=%.1f, signif=%.3f" %(SplusB, SplusB-B, significance)
+
+

apz/run.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+import sys,os
+from math import sqrt
+import argparse
+parser =  argparse.ArgumentParser(description='Running the tree selection', usage="./run.py")
+parser.add_argument("--proof",  dest="proof", action="store_true", default=False,  help="Run this thing with PROOF!")
+parser.add_argument("-c",  dest="cuts", type=str, default='default',
+                    choices=['default','optimal','scan'],help="Run this thing with PROOF!")
+
+opt = parser.parse_args()
+
+from ROOT import *
+
+workingPath = os.getcwd()
+parentDir = os.path.abspath(os.path.join(workingPath, os.pardir))
+SandM_path = parentDir+"/sugar-n-milk/"
+
+gSystem.SetBuildDir("buildDir", kTRUE)
+gSystem.AddIncludePath(" -I"+SandM_path)
+gROOT.LoadMacro(SandM_path+"/HistManager.cc+")
+
+#fChain = TChain("apzTree")
+#fChain.Add('APZ_smalltree_25July2018.root')
+
+fChain = TChain("apzTree/apzTree")
+fChain.Add('APZ-tree-Data-MuEG-8TeV-2018July25.root')
+
+timer = TStopwatch()
+timer.Start()
+
+optima = []
+
+# Default cuts:
+gamma_pt_cuts = [25]
+gamma_eta_cuts = [1.4442]
+lep1_pt_cuts = [23]
+lep2_pt_cuts = [4]
+
+if opt.cuts=='scan':
+  gamma_pt_cuts = range(26,34)
+  gamma_eta_cuts = [1.2, 1.3, 1.4, 1.4442, 1.5, 1.6, 1.8]
+  lep1_pt_cuts = range(23,26)
+  lep2_pt_cuts = [3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0]
+if opt.cuts=='optimal':
+  # Optimal cuts using 9 bins in signal region:
+  gamma_pt_cuts = [30]
+  gamma_eta_cuts = [1.5]
+  lep1_pt_cuts = [24]
+  lep2_pt_cuts = [6]
+  # Optimal cuts using 5 bins in signal region:
+  #gamma_pt_cuts = [30]
+  #gamma_eta_cuts = [1.4442]
+  #lep1_pt_cuts = [23]
+  #lep2_pt_cuts = [6]
+
+for cut_gamma_pt in gamma_pt_cuts:
+  for cut_gamma_eta in gamma_eta_cuts:
+    for cut_lep1_pt in lep1_pt_cuts:
+      for cut_lep2_pt in lep2_pt_cuts: 
+        fChain.Process("smallAna.C+", "%f %f %f %f" % (cut_gamma_pt, cut_gamma_eta, cut_lep1_pt, cut_lep2_pt))
+
+        f = TFile("my_higgsHistograms.root",'OPEN')  
+        evt = f.Get("evt")
+        SplusB = evt.GetBinContent(3)
+        B = (9./11)*evt.GetBinContent(4)
+        S = SplusB - B;
+        significance = 2*(sqrt(SplusB) - sqrt(B))
+        print "S+B=%i, S=%.1f, signif=%.3f" %(SplusB, S, significance)
+
+        info = [cut_gamma_pt, cut_gamma_eta, cut_lep1_pt, cut_lep2_pt, SplusB, S, significance]
+        optima.append(info)
+
+print optima
+sort_opt = sorted(optima, key=lambda x: float(x[6]))
+print 
+print 
+for l in sort_opt:
+  print l
+
+if opt.cuts!='scan':
+  print "Made histos, now let's make the plots"
+  execfile('plotit.py')
+
+print sort_opt[-1]
+print "Done!", "CPU Time: ", timer.CpuTime(), "RealTime : ", timer.RealTime()
+
+
+
+print "End of running"