compareDigi.py

import ROOT
import random
import os
import sys
import numpy as np
import array
import math

ROOT.gROOT.SetBatch(0)
ROOT.gStyle.SetStatW(0.07)
ROOT.gStyle.SetStatH(0.06)

#ROOT.gStyle.SetOptStat(111110)

#ROOT.gStyle.SetErrorX(0)
#ROOT.gStyle.SetErrorY(0)

ROOT.gStyle.SetTitleStyle(0)
ROOT.gStyle.SetTitleAlign(13) ## coord in top left
ROOT.gStyle.SetTitleX(0.)
ROOT.gStyle.SetTitleY(1.)
ROOT.gStyle.SetTitleW(1)
ROOT.gStyle.SetTitleH(0.058)
ROOT.gStyle.SetTitleBorderSize(0)

ROOT.gStyle.SetPadLeftMargin(0.126)
ROOT.gStyle.SetPadRightMargin(0.10)
ROOT.gStyle.SetPadTopMargin(0.06)
ROOT.gStyle.SetPadBottomMargin(0.13)

ROOT.gStyle.SetMarkerStyle(1)


xmin=0
xmax=60
xbins=30
BINM=22
offsetLUT = {}
#binLow = [0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,12.0,14.0,16.0,18.0,20.0,24.0,28.0,32.0,36.0,42.0,50.0,60.0]
binLow = [0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,12.0,14.0,16.0,18.0,20.0,24.0,28.0,32.0,36.0,42.0,50.0]
ptbins = np.asarray(binLow)
	
def gethist1D(chain, den, todraw, x_bins):
	
    xBins = int(x_bins[1:-1].split(',')[0])
    xminBin = float(x_bins[1:-1].split(',')[1])
    xmaxBin = float(x_bins[1:-1].split(',')[2])
    h1 = ROOT.TH1F("h1","h1",xBins,xminBin,xmaxBin)
    chain.Draw("fabs(%s) >> h1"%todraw,den)
    #print "gethist1D, den ",den
    #h1.Print("ALL")
    return h1
		



def draw2Dplots_1(ch, xaxis, yaxis, x_bins, y_bins, xtitle, ytitle,cuts, text, picname):

    xBins = int(x_bins[1:-1].split(',')[0])
    xminBin = float(x_bins[1:-1].split(',')[1])
    xmaxBin = float(x_bins[1:-1].split(',')[2])
    yBins = int(y_bins[1:-1].split(',')[0])
    yminBin = float(y_bins[1:-1].split(',')[1])
    ymaxBin = float(y_bins[1:-1].split(',')[2])

    todrawb0 = "%s"%yaxis+":"+"%s>>b0"%xaxis
    
    c0 = ROOT.TCanvas("c0","c0")
    c0.SetGridx()
    c0.SetGridy()
    c0.SetTickx()
    c0.SetTicky()
    b0 = ROOT.TH2F("b0","b0",xBins,xminBin,xmaxBin,yBins,yminBin,ymaxBin)
    #b0.SetTitle("#scale[1.4]{#font[61]{CMS}} #font[52]{ preliminary}"+"  "*15+"Run=281976")
    b0.SetTitle("#scale[1.4]{#font[61]{CMS}} #font[52]{preliminary}"+"  "*15+"2017, RunH")
    #b0.SetTitle("%s Vs %s,%s"%(ytitle0, xtitle, st_title)) 
    #b1.SetTitleFont(62)
    b0.SetTitleSize(0.05)
    b0.GetXaxis().SetTitle("%s"%xtitle)
    b0.GetYaxis().SetTitle("%s"%ytitle)
    b0.GetXaxis().SetTitleSize(.05)
    b0.GetYaxis().SetTitleSize(.05)
    #b1.SetMaximum(30)
    b0.SetStats(0)
    print "todraw ",todrawb0, " cut ",cuts
    #hasxy = "&& fabs(%s)>0 && fabs(%s)>0"%(xaxis, yaxis0)
    drawoption = "colz"
    if xBins < 20 and yBins < 20:
        drawoption = "colztext"
    ch.Draw(todrawb0,cuts, drawoption)

    #tex1 = TLatex(0.35,.8,"#splitline{%s}{%.1f<|#eta|<%.1f,20<p_{T}<50}"%(st_title,etamin,etamax))
    #tex1 = ROOT.TLatex(0.15,.8,"%s"%(text))
    tex1 = ROOT.TLatex(0.2,.86,"%s"%(text))
    tex1.SetTextSize(0.05)
    tex1.SetTextFont(62)
    tex1.SetNDC()
    tex1.Draw("same")
    #c0.SaveAs("%s"%picname+".png")
    c0.SaveAs("%s"%picname+".pdf")
    c0.SaveAs("%s"%picname+".C")

def draw1D_compare(filelist, chs, xaxis_list, x_bins, xtitle, legs, cuts, text, picname):

    xBins = int(x_bins[1:-1].split(',')[0])
    xminBin = float(x_bins[1:-1].split(',')[1])
    xmaxBin = float(x_bins[1:-1].split(',')[2])
    tfilelist = []
    for f in filelist:
	tfilelist.append(ROOT.TFile(f, "READ"))
    color = [ROOT.kBlue, ROOT.kRed, ROOT.kGreen+2,ROOT.kMagenta+2, ROOT.kCyan+2]
    maker = [20,21,22,23,33]
    #title = "#scale[1.4]{#font[61]{CMS}} #font[52]{ preliminary}"+"  "*15+"Run=281976"
    title = "#scale[1.4]{#font[61]{CMS}} #font[52]{ Simulation preliminary}"+"  "*15
    hs1 = ROOT.THStack("hs1","%s"%title)
    hs2 = ROOT.THStack("hs2","%s"%title)
    legend = ROOT.TLegend(0.15,0.55,0.5,0.7)
    legend.SetFillColor(ROOT.kWhite)
    #legend.SetTextFont(42)
    #legend.SetTextSize(.05)
    
    i=0
    meanlist = []
    for ch in chs:
	tfilelist[i].cd()
	print "todraw ","%s>>hist%d"%(xaxis_list[i], i)," cut ",cuts[i]
	hist = ROOT.TH1F("hist%d"%i,"hist%d"%i, xBins, xminBin, xmaxBin)
	ch.Draw("%s>>hist%d"%(xaxis_list[i], i),cuts[i])
	ROOT.SetOwnership(hist, False)
	hist.SetLineColor(color[i])
	hist.SetLineWidth(2)
	hist.Sumw2()
	hs1.Add(hist)
	hs2.Add(hist.Scale(1.0/hist.Integral()))
	mean = hist.GetMean()
	meanlist.append(mean)
	rms = hist.GetRMS()
	#legend.AddEntry(hist, "%s, Mean: %.4f, RMS: %.4f"%(legs[i], mean, rms),"l")
	legend.AddEntry(hist, "%s, Mean: %.4f"%(legs[i], mean),"l")
	i +=1

    c0 = ROOT.TCanvas("c0","c0")
    c0.SetGridx()
    c0.SetGridy()
    c0.SetTickx()
    c0.SetTicky()
    c0.SetLogy()
    hs1.SetMinimum(.001)
    hs1.Draw("nostacke") 
    hs1.GetHistogram().GetXaxis().SetTitle("%s"%xtitle)
    hs1.GetHistogram().GetXaxis().SetTitleSize(.05)
    hs1.GetHistogram().GetYaxis().SetTitle("Normalized to unity")
    hs1.GetHistogram().GetYaxis().SetTitleSize(.05)
    #hs1.GetHistogram().GetYaxis().SetRangeUser(.001,2)
    legend.Draw("same")
    tex1 = ROOT.TLatex(0.2,.8,"%s"%(text))
    tex1.SetTextSize(0.05)
    tex1.SetTextFont(62)
    tex1.SetNDC()
    tex1.Draw("same")
    c0.SaveAs("%s"%picname+".C")
    #c0.SaveAs("%s"%picname+".png")
    c0.SaveAs("%s"%picname+".pdf")
#    c0.SaveAs("%s"%picname+".C")

    return meanlist
    """
    c1 = ROOT.TCanvas("c1","c1")
    c1.SetGridx()
    c1.SetGridy()
    c1.SetTickx()
    c1.SetTicky()
    hs2.Draw("nostack") 
    hs2.GetHistogram().GetXaxis().SetTitle("%s"%xtitle)
    hs2.GetHistogram().GetYaxis().SetTitle("Normalized to unity")
    legend.Draw("same")
    tex1.SetNDC()
    tex1.Draw("same")
    c1.SaveAs("%s"%picname+"_normalized.png")
    c1.SaveAs("%s"%picname+"_normalized.C")
    """

#######################################################


cscstations = [ [0,0], 
                [1,1], [1,2], [1,3],[1,4],
                [2,1], [2,2],
                [3,1], [3,2],
                [4,1], [4,2],]
chambernames = ["all",
		"ME1/b","ME1/2","ME1/3","ME1/a",
		"ME2/1","ME2/2",
		"ME3/1","ME3/2",
		"ME4/1","ME4/2",]
allvars = ["quality","bend","pattern","key_WG","key_hs", "bx", "totStubs", "trknmb"]
xbinsall = ["(17,-1.5,15.5)","(3,-1.5,1.5)","(13,-1.5,11.5)","(122,-1.5,120.5)","(162,-1.5,160.5)", "(17,-1.5,15.5)", "(7, -0.5, 6.5)", "(4,-1.5, 2.5)"]
bx_bins = "(16, 0.5,16.5)"


xtitle = "timing of LCTs, BX"
xaxis = "bx"
xaxis_list = [xaxis, xaxis]

#draw1D_compare(chs, xaxis_list, bx_bins, xtitle, legs, cuts, text, outputdir+"stubscompare_bx_emualtor")
def compareDataAndEmulation(rootfile, outputdir):
    for chname in ['alcttree','clcttree','lcttree']:
    #for chname in ['lcttree']:
        obj = chname[:-4]
        ch = ROOT.TChain("lctreader/"+chname)
        ch.Add(rootfile)
        for ivar, var in enumerate(allvars):
            x_bins = xbinsall[ivar]
            xaxis = var+"_data"; yaxis = var+"_emul"
            #cuts = xaxis+">=0 && "+yaxis+">=0"
            #cuts =  xaxis+">=0"
            cuts = "(bx_corr_emul>=0 || bx_data>=0)"
            xtitle = obj+","+var +", Data"
            ytitle = obj+","+var +", Re-emulation"
            if var == "bx":
                yaxis = "bx_corr_emul"
                ytitle = obj+",converted bx"+", Re-emulation"
                if chname == 'alcttree':
                    x_bins = "(10, -1.5, 8.5)"
                elif chname == 'clcttree':
                    x_bins = "(6, -1.5, 4.5)"
                elif chname == "lcttree":
                    x_bins = "(4, -1.5, 2.5)"

            if var == "quality" and (chname == 'alcttree' or chname == 'clcttree'):
                x_bins = "(8,-1.5, 6.5)"
            text = "Data Vs Re-emulation"
            #print "text ",text," xaxis ",xaxis," yaxis ",yaxis," cuts ",cuts
            picname = os.path.join(outputdir, obj+"_"+var+"_data_reemulation")
            draw2Dplots_1(ch, xaxis, yaxis, x_bins, x_bins, xtitle, ytitle,cuts, text, picname)
            for ichambertype in range(1, len(cscstations)):
                chambername = chambernames[ichambertype]
                st = cscstations[ichambertype][0]
                ring = cscstations[ichambertype][1]
                if ring == 4:
                    continue
                cuts_ch = "station == %d && ring ==%d "%(st, ring)
                picname = os.path.join(outputdir, obj+"_"+var+"_data_reemulation_st%d_ring%d"%(st, ring))
                if ichambertype == 1:
                    chambername = "ME1/1"
                    #for k in range(1, 37):
                    #    text_pme11 = "Data Vs Re-emulation, p"+chambername+", chamber%d"%k
                    #    picname_pme11 = os.path.join(outputdir, obj+"_"+var+"_data_reemulation_st%d_ring%d_chamber%d_endcap1"%(st, ring, k))
                    #    cuts_pme11 = cuts+" && "+cuts_ch+" && chamber ==%d && endcap==1"%k
                    #    draw2Dplots_1(ch, xaxis, yaxis, x_bins, x_bins, xtitle, ytitle, cuts_pme11, text_pme11, picname_pme11)
                    #    text_mme11 = "Data Vs Re-emulation, m"+chambername+", chamber%d"%k
                    #    picname_mme11 = os.path.join(outputdir, obj+"_"+var+"_data_reemulation_st%d_ring%d_chamber%d_endcap2"%(st, ring, k))
                    #    cuts_mme11 = cuts+" && "+cuts_ch+" && chamber ==%d && endcap==2"%k
                    #    draw2Dplots_1(ch, xaxis, yaxis, x_bins, x_bins, xtitle, ytitle, cuts_mme11, text_mme11, picname_mme11)
                text = "Data Vs Re-emulation, "+chambername
                draw2Dplots_1(ch, xaxis, yaxis, x_bins, x_bins, xtitle, ytitle,cuts+" && "+cuts_ch, text, picname)




rootfile = "TPEHists_run2017H_20171213.root"
outputdir = "Run2017H_data_reemulation_20171213/"

os.system("mkdir -p "+outputdir)
compareDataAndEmulation(rootfile, outputdir)