diff --git a/apz/copytree.C b/apz/copytree.C new file mode 100644 index 0000000..402f72a --- /dev/null +++ b/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; +} diff --git a/apz/plotit.py b/apz/plotit.py new file mode 100755 index 0000000..118d1d2 --- /dev/null +++ b/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) + + diff --git a/apz/run.py b/apz/run.py new file mode 100755 index 0000000..883e472 --- /dev/null +++ b/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" diff --git a/apz/smallAna.C b/apz/smallAna.C new file mode 100644 index 0000000..39c8d77 --- /dev/null +++ b/apz/smallAna.C @@ -0,0 +1,206 @@ +#define smallAna_cxx + +#include "smallAna.h" +#include +#include +#include +#include +#include + +void smallAna::Begin(TTree * /*tree*/) +{ + TString option = GetOption(); + + +} + +void smallAna::SlaveBegin(TTree * /*tree*/) +{ + TString option = GetOption(); + cout<<"Options are: \n \t "<At(0))->GetString()); + _cut_gamma_eta = std::stod((string)((TObjString*)args->At(1))->GetString()); + _cut_lep1_pt = std::stod((string)((TObjString*)args->At(2))->GetString()); + _cut_lep2_pt = std::stod((string)((TObjString*)args->At(3))->GetString()); + + //histoFile = new TFile("my_higgsHistograms.root", "RECREATE"); + //fOutput->Add(histoFile); + fProofFile = new TProofOutputFile("my_higgsHistograms.root"); + fOutput->Add(fProofFile); + // Open the file + histoFile = fProofFile->OpenFile("RECREATE"); + + hists = new HistManager(histoFile); + for (size_t n=0; nFill(0); + + //if (nEvents[0]%10000 == 0){ + //Info("PROC", "Lpt pt = %.2f", lep1->Pt()); + + //std::cout<<"lep1 pt = "<Pt()<Pt() < _cut_lep1_pt) return kTRUE; // Original cut is 23 + if (lep2->Pt() < _cut_lep2_pt) return kTRUE; // Original cut is 4 + if (gamma->Pt() < _cut_gamma_pt) return kTRUE; // Original cut is 25 + if (abs(gamma->Eta()) > _cut_gamma_eta) return kTRUE; // Original cut is 1.4442 + //if (Mtri<110 || Mtri>170) return kTRUE; + if (lep1->DeltaR(*gamma)<1 || lep2->DeltaR(*gamma)<1) return kTRUE; // Original cut is 1.0 + + //hists->fill1DHist(lep1->Pt(), "lep1_pt_test",";p_T(l1)", 100, 0,100, 1,"DIR"); + + if (!trig1) return kTRUE; + + h1->Fill(1); + + nEvents[1]++; + makeHists("AfterCuts"); + + Float_t mll = (*lep1+*lep2).M(); + Float_t mllg = (*lep1+*lep2+*gamma).M(); + if (mll>3.0 && mll<3.2) { + nEvents[2]++; + makeHists("InJpsi"); + h1->Fill(2); + } + + if (mll>=2.76 && mll<=2.85) { + nEvents[3]++; + makeHists("InAPZ"); + h1->Fill(3); + } + + if ( ( mll>2.70 && mll<2.76) || (mll>2.85 && mll<2.9 ) ) { + nEvents[4]++; + makeHists("InSideBand"); + h1->Fill(4); + } + + const Float_t Mllg = (*lep1+*lep2+*gamma).M(); + + //if ( gamma->Pt()/Mllg > 0.35 && (*lep1+*lep2).Pt()/Mllg > 0.35 + // && Mllg>100 && Mllg<150) { // Optimal? + if ( gamma->Pt()/Mllg > 0.3 && (*lep1+*lep2).Pt()/Mllg > 0.30 + && Mllg>110 && Mllg<170) { // PAS + nEvents[5]++; + makeHists("M18"); + h1->Fill(5); + } + + h1->Fill(6); + + return kTRUE; +} + +void smallAna::makeHists(string tag){ + //void smallAna::makeHists(string tag, const TLorentzVector& l1, const TLorentzVector& l2, const TLorentzVector& g){ + const TLorentzVector diLep = (*lep1+*lep2); + const TLorentzVector tri = (*lep1+*lep2+*gamma); + const Float_t Mll = diLep.M(); + const Float_t Mjj = (*jet1+*jet2).M(); + const Float_t Mllg = tri.M(); + const Float_t Mllgj = (*lep1+*lep2+*gamma+*jet1).M(); + const Float_t Mgj = (*gamma+*jet1).M(); + + //const TVector2 diLep2(diLep.X(), diLep.Y()); + //const TVector2 gamma2(gamma->X(), gamma->Y()); + + hists->fill1DHist(lep1->Pt(), "lep1_pt_"+tag,";p_{T}(l1)", 60, 20,80, 1,tag); + hists->fill1DHist(lep2->Pt(), "lep2_pt_"+tag,";p_{T}(l2)", 50, 0,50, 1,tag); + hists->fill1DHist(gamma->Pt(),"gamma_pt_"+tag,";p_{T}(#gamma)", 100, 20,100, 1,tag); + + + hists->fill1DHist(lep1->Phi(), "lep1_phi_"+tag,";#phi(l1)", 50, -3.2, 3.2, 1,tag); + hists->fill1DHist(lep2->Phi(), "lep2_phi_"+tag,";#phi(l2)", 50, -3.2, 3.2, 1,tag); + hists->fill1DHist(gamma->Phi(), "gamma_phi_"+tag,";#phi(#gamma)", 50, -3.2, 3.2, 1,tag); + + + hists->fill1DHist(lep1->Eta(), "lep1_eta_"+tag,";#eta(l1)", 50, -2.4,2.4, 1,tag); + hists->fill1DHist(lep2->Eta(), "lep2_eta_"+tag,";#eta(l2)", 50, -2.4,2.4, 1,tag); + hists->fill1DHist(gamma->Eta(), "gamma_eta_"+tag,";#eta(#gamma)", 50, -2.4,2.4, 1,tag); + + + hists->fill1DHist(Mll, "Mll_50_"+tag,";m_{#mu#mu}", 100, 0,50, 1,tag); + hists->fill1DHist(Mll, "Mll_20_"+tag,";m_{#mu#mu}", 100, 0,20, 1,tag); + hists->fill1DHist(Mll, "Mll_20_to_50_"+tag,";m_{#mu#mu}", 50, 20,50, 1,tag); + hists->fill1DHist(Mll, "Mll_18_"+tag,";m_{#mu#mu}", 50, 10,30, 1,tag); + hists->fill1DHist(Mll, "Mll_Ups_"+tag,";m_{#mu#mu}", 60, 9,12, 1,tag); + hists->fill1DHist(Mll, "Mll_JPsi_"+tag,";m_{#mu#mu}", 50, 2,4, 1,tag); + hists->fill1DHist(Mll, "Mll_APZ_JPsi_"+tag,";m_{#mu#mu}", 70, 2.7, 3.4, 1,tag); + hists->fill1DHist(Mll, "Mll_APZ_"+tag,";m_{#mu#mu}", 30, 2.7,3.0, 1,tag); + + hists->fill1DHist(Mllg, "Mllg_full_"+tag,";m_{#mu#mu#gamma}", 100, 0,300, 1,tag); + hists->fill1DHist(Mllg, "Mllg_higg_"+tag,";m_{#mu#mu#gamma}", 50, 100,150, 1,tag); + hists->fill1DHist(Mllg, "Mllg_zed_"+tag,";m_{#mu#mu#gamma}", 40, 50,90, 1,tag); + + + hists->fill1DHist(lep1->DeltaR(*lep2), "dR_lep_"+tag,";#DeltaR(l1,l2)", 100, 0,1.5, 1,tag); + hists->fill1DHist(gamma->DeltaR(*lep1+*lep2), "dR_gamma_dilep_"+tag,";#DeltaR(#gamma,#mu#mu)", 50, 1, 4.5, 1,tag); + + hists->fill1DHist(diLep.Pt()/gamma->Pt(),"pt_ratio_dilep_gamma_"+tag,";p_{T}(#mu#mu)/p_{T}(#gamma)", 50, 0,2, 1,tag); + + hists->fill1DHist(met->Pt(), "met_Et_"+tag,";MET", 50, 0, 100, 1,tag); + hists->fill1DHist(met->Phi(), "met_phi_"+tag,";#phi(MET)", 50, -3.2, 3.2, 1,tag); + hists->fill1DHist(met->DeltaPhi(diLep), "met_dPhi_dilep_"+tag,";#Delta#phi(MET,#mu#mu)", 50, 0, 3.2, 1,tag); + hists->fill1DHist(met->DeltaPhi(*gamma), "met_dPhi_gamma_"+tag,";#Delta#phi(MET,#gamma)", 50, 0, 3.2, 1,tag); + + if (jet1->Pt()>0){ + hists->fill1DHist(jet1->Pt(), "jet1_pt_"+tag, ";p_{T}(j1)", 50, 20, 120, 1,tag); + hists->fill1DHist(jet1->Eta(), "jet1_eta_"+tag,";#eta(j1)", 50, -2.4, 2.4, 1,tag); + hists->fill1DHist(jet1->Phi(), "jet1_phi_"+tag,";#phi(j1)", 50, -3.2, 3.2, 1,tag); + + hists->fill1DHist(Mllgj, "Mllgj_full_"+tag,";m_{#mu#mu#gammaj}", 100, 50,500, 1,tag); + hists->fill1DHist(Mllgj, "Mllgj_higg_"+tag,";m_{#mu#mu#gammaj}", 20, 100,200, 1,tag); + hists->fill1DHist(Mllgj, "Mllgj_x200_"+tag,";m_{#mu#mu#gammaj}", 25, 180,230, 1,tag); + + hists->fill1DHist(Mgj, "Mgj_full_"+tag,";m_{#gammaj}", 50, 0,500, 1,tag); + hists->fill1DHist(Mgj, "Mgj_higg_"+tag,";m_{#gammaj}", 50, 100,200, 1,tag); + + hists->fill1DHist(gamma->DeltaR(*jet1), "dR_gamma_jet1_"+tag,";#DeltaR(#gamma,jet1)", 50, 0,5, 1,tag); + + if (jet2->Pt()>0){ + hists->fill1DHist(jet2->Pt(), "jet2_pt_"+tag, ";p_{T}(j2)", 50, 20, 100, 1,tag); + hists->fill1DHist(jet2->Eta(), "jet2_eta_"+tag,";#eta(j2)", 50, -2.4, 2.4, 1,tag); + hists->fill1DHist(jet2->Phi(), "jet2_phi_"+tag,";#phi(j2)", 50, -3.2, 3.2, 1,tag); + + hists->fill1DHist(Mjj, "Mjj_full_"+tag,";m_{jj}", 50, 0, 600, 1,tag); + hists->fill1DHist(Mjj, "Mjj_zoom_"+tag,";m_{jj}", 50, 100, 400, 1,tag); + + hists->fill1DHist(gamma->DeltaR(*jet2), "dR_gamma_jet2_"+tag,";#DeltaR(#gamma,jet2)", 50, 0,5, 1,tag); + + hists->fill1DHist(jet1->DeltaR(*jet2), "dRjj_"+tag,";#DeltaR(jj)", 50, 0,5, 1,tag); + } + } + +} + +void smallAna::SlaveTerminate() +{ + histoFile->cd(); + histoFile->Write(); + histoFile->Close(); + + cout<<" ** Slave is terminating, my master ...bow x3... **"< +#include +#include +#include +#include +// Header file for the classes stored in the TTree if any. +#include "TLorentzVector.h" + +#include "HistManager.h" + +#define nC 10 +class TFile; +class TProofOutputFile; + +class smallAna : public TSelector { + private: + + TFile* histoFile; + //TFile *fFile; + TProofOutputFile *fProofFile; // For optimized merging of the ntuple + HistManager *hists; + TH1F* h1; + UInt_t nEvents[nC]; + UInt_t totEvents; + + Float_t _cut_gamma_pt, _cut_gamma_eta, _cut_lep1_pt, _cut_lep2_pt; + + public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + + // Fixed size dimensions of array or collections stored in the TTree if any. + + // Declaration of leaf types + TLorentzVector *lep1; + TLorentzVector *lep2; + TLorentzVector *gamma; + TLorentzVector *jet1; + TLorentzVector *jet2; + TLorentzVector *met; + UInt_t run; + ULong64_t event; + + Bool_t trig1, trig2, trig3, trig4, trig5; + + // List of branches + TBranch *b_lep1; //! + TBranch *b_lep2; //! + TBranch *b_gamma; //! + TBranch *b_jet1; //! + TBranch *b_jet2; //! + TBranch *b_met; //! + TBranch *b_run; //! + TBranch *b_event; //! + + TBranch *b_trig1; + TBranch *b_trig2; + TBranch *b_trig3; + TBranch *b_trig4; + TBranch *b_trig5; + + smallAna(TTree * /*tree*/ =0) : fChain(0), fProofFile(0) { } + virtual ~smallAna() { } + virtual Int_t Version() const { return 2; } + virtual void Begin(TTree *tree); + virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); + virtual Bool_t Notify(); + virtual Bool_t Process(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; } + virtual void SetOption(const char *option) { fOption = option; } + virtual void SetObject(TObject *obj) { fObject = obj; } + virtual void SetInputList(TList *input) { fInput = input; } + virtual TList *GetOutputList() const { return fOutput; } + virtual void SlaveTerminate(); + virtual void Terminate(); + + + virtual void makeHists(string ); + + ClassDef(smallAna,0); +}; + +#endif + +#ifdef smallAna_cxx +void smallAna::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set object pointer + lep1 = 0; + lep2 = 0; + gamma = 0; + jet1 = 0; + jet2 = 0; + met = 0; + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("lep1", &lep1, &b_lep1); + fChain->SetBranchAddress("lep2", &lep2, &b_lep2); + fChain->SetBranchAddress("gamma", &gamma, &b_gamma); + fChain->SetBranchAddress("jet1", &jet1, &b_jet1); + fChain->SetBranchAddress("jet2", &jet2, &b_jet2); + fChain->SetBranchAddress("met", &met, &b_met); + fChain->SetBranchAddress("run", &run, &b_run); + fChain->SetBranchAddress("event", &event, &b_event); + + fChain->SetBranchAddress("trig1", &trig1, &b_trig1); + fChain->SetBranchAddress("trig2", &trig2, &b_trig2); + fChain->SetBranchAddress("trig3", &trig3, &b_trig3); + fChain->SetBranchAddress("trig4", &trig4, &b_trig4); + fChain->SetBranchAddress("trig5", &trig5, &b_trig5); +} + +Bool_t smallAna::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + //TFile *inFile = thisTree->GetCurrentFile(); + //cout<<" Opened file: \n"<GetName()<GetCurrentFile()->GetName()<