analysis_INTT/plot/ClusADCCutComp.py

0001 #! /usr/bin/env python
0002 from optparse import OptionParser
0003 import sys
0004 import os
0005 import datetime
0006 from array import *
0007 from ROOT import *
0008 import numpy
0009 import math
0010 import glob
0011 from plotUtil import *
0012
0013 gROOT.SetBatch(True)
0014
0015 TickSize = 0.03
0016 AxisTitleSize = 0.05
0017 AxisLabelSize = 0.04
0018 LeftMargin = 0.15
0019 RightMargin = 0.08
0020 TopMargin = 0.08
0021 BottomMargin = 0.13
0022
0023 def Draw_1Dhist_datasimcomp(hdatas, hsims, gpadmargin, norm, logy, ymaxscale, XaxisName, Ytitle_unit, prelim, simlegtex, evtseltexts, outname):
0024     if len(hdatas) != len(hsims):
0025         print('Number of data histograms and simulation histograms do not match')
0026         sys.exit(1)
0027
0028     hsimcolor = colorset(len(hsims))
0029
0030     for hdata in hdatas:
0031         hdata.Sumw2()
0032     for hsim in hsims:
0033         hsim.Sumw2()
0034
0035     binwidth = hdatas[0].GetXaxis().GetBinWidth(1)
0036
0037     if norm == 'unity':
0038         for hdata in hdatas:
0039             hdata.Scale(1. / hdata.Integral(-1, -1))
0040         for hsim in hsims:
0041             hsim.Scale(1. / hsim.Integral(-1, -1))
0042     elif norm == 'data':
0043         for i, hsim in enumerate(hsims):
0044             hsim.Scale(hdatas[i].Integral(-1, -1) / hsim.Integral(-1, -1))
0045     else:
0046         if norm != 'none':
0047             print('Invalid normalization option: {}'.format(norm))
0048             sys.exit(1)
0049
0050     # Get the maximum bin content
0051     maxbincontent = max([hdata.GetMaximum() for hdata in hdatas] + [hsim.GetMaximum() for hsim in hsims])
0052
0053     c = TCanvas('c', 'c', 800, 700)
0054     pad1 = TPad( 'pad1', ' ', 0, 0.3, 1, 1)
0055     pad2 = TPad( 'pad2', ' ', 0, 0, 1, 0.3)
0056     pad1.SetRightMargin(gpadmargin[0])
0057     pad1.SetTopMargin(gpadmargin[1])
0058     pad1.SetLeftMargin(gpadmargin[2])
0059     pad1.SetBottomMargin(0.0)
0060     pad1.Draw()
0061     pad2.SetRightMargin(gpadmargin[0])
0062     pad2.SetTopMargin(0.0)
0063     pad2.SetLeftMargin(gpadmargin[2])
0064     pad2.SetBottomMargin(gpadmargin[3])
0065     pad2.Draw() # Draw the TPad on the TCanvas before plotting something on TPad
0066     # cd to the pad1
0067     pad1.cd()
0068     if logy:
0069         pad1.SetLogy()
0070
0071     for i, hsim in enumerate(hsims):
0072         if i == 0:
0073             if norm == 'unity':
0074                 if Ytitle_unit == '':
0075                     hsim.GetYaxis().SetTitle(
0076                         'Normalized entries / ({:g})'.format(binwidth))
0077                 else:
0078                     hsim.GetYaxis().SetTitle(
0079                         'Normalized entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0080             else:
0081                 if Ytitle_unit == '':
0082                     hsim.GetYaxis().SetTitle('Entries / ({:g})'.format(binwidth))
0083                 else:
0084                     hsim.GetYaxis().SetTitle(
0085                         'Entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0086
0087             if logy:
0088                 hsim.GetYaxis().SetRangeUser(hsim.GetMinimum(0)*0.5, (hsim.GetMaximum()) * ymaxscale)
0089             else:
0090                 hsim.GetYaxis().SetRangeUser(1E-3, (hsim.GetMaximum()) * ymaxscale)
0091
0092             hsim.GetXaxis().SetTitle(XaxisName)
0093             hsim.GetXaxis().SetTitleOffset(1.1)
0094             hsim.GetYaxis().SetTitleOffset(1.35)
0095             hsim.SetLineColorAlpha(TColor.GetColor(hsimcolor[i]), 0.95)
0096             hsim.SetLineWidth(2)
0097             hsim.SetMarkerSize(0)
0098             hsim.Draw('histe')
0099         else:
0100             hsim.SetLineColorAlpha(TColor.GetColor(hsimcolor[i]), 0.95)
0101             hsim.SetLineWidth(2)
0102             hsim.SetMarkerSize(0)
0103             hsim.Draw('histe same')
0104
0105     for i, hdata in enumerate(hdatas):
0106         hdata.SetMarkerStyle(20)
0107         hdata.SetMarkerSize(0.7)
0108         hdata.SetMarkerColorAlpha(TColor.GetColor(hsimcolor[i]), 0.95)
0109         hdata.SetLineColorAlpha(TColor.GetColor(hsimcolor[i]), 0.95)
0110         hdata.SetLineWidth(2)
0111         hdata.Draw('same PE1')
0112     shift = 0.45 if prelim else 0.75
0113     legylow = 0.2 + 0.04 * (len(hsims) - 1)
0114     leg = TLegend((1-RightMargin)-shift, (1-TopMargin)-legylow,
0115                   (1-RightMargin)-0.1, (1-TopMargin)-0.03)
0116     leg.SetTextSize(0.04)
0117     leg.SetFillStyle(0)
0118     # prelimtext = 'Preliminary' if prelim else 'Work-in-progress'
0119     prelimtext = 'Preliminary' if prelim else 'Internal'
0120     leg.AddEntry('', '#it{#bf{sPHENIX}} '+prelimtext, '')
0121     leg.AddEntry('', 'Au+Au #sqrt{s_{NN}}=200 GeV', '')
0122     for i, lt in enumerate(simlegtex):
0123         leg.AddEntry(hdatas[i], lt, "lpe")
0124     leg.Draw('same')
0125
0126     # event selection text
0127     legylow_evtselshift = 0.04 * len(evtseltexts)
0128     leg2 = TLegend((1-RightMargin)-shift, (1-TopMargin)-legylow-legylow_evtselshift,
0129                     (1-RightMargin)-0.1, (1-TopMargin)-legylow)
0130     leg2.SetTextSize(0.035)
0131     leg2.SetFillStyle(0)
0132     for i, evtseltext in enumerate(evtseltexts):
0133         leg2.AddEntry('', evtseltext, '')
0134     leg2.Draw('same')
0135     c.RedrawAxis()
0136
0137     c.Update()
0138     # cd to the pad2
0139     textscale = 2.7
0140     pad2.cd()
0141     pad2.SetLogy(0)
0142     pad2.SetGridy(0)
0143     # take the ratio = data/simulation
0144     l_hM_ratio = []
0145     for i, hsim in enumerate(hsims):
0146         hM_ratio = hdatas[i].Clone('hM_ratio'+str(i))
0147         hM_ratio.Divide(hsim)
0148         hM_ratio.SetMarkerStyle(20)
0149         hM_ratio.SetMarkerSize(0.7)
0150         hM_ratio.SetMarkerColorAlpha(TColor.GetColor(hsimcolor[i]), 0.95)
0151         hM_ratio.SetLineColorAlpha(TColor.GetColor(hsimcolor[i]), 0.95)
0152         hM_ratio.SetLineWidth(2)
0153         if i == 0:
0154             hM_ratio.GetYaxis().SetRangeUser(0.8+1E-3, 1.2-1E-3)
0155             hM_ratio.GetYaxis().SetNdivisions(505)
0156             hM_ratio.GetYaxis().SetTitle('Data/Sim.')
0157             hM_ratio.GetYaxis().SetTitleOffset(0.5)
0158             hM_ratio.GetYaxis().SetTickSize(TickSize)
0159             hM_ratio.GetYaxis().SetTitleSize(AxisTitleSize*textscale)
0160             hM_ratio.GetYaxis().SetLabelSize(AxisLabelSize*textscale)
0161             hM_ratio.GetXaxis().SetTickSize(TickSize*textscale)
0162             hM_ratio.GetXaxis().SetTitleSize(AxisTitleSize*textscale)
0163             hM_ratio.GetXaxis().SetLabelSize(AxisLabelSize*textscale)
0164             hM_ratio.GetXaxis().SetTitleOffset(1.1)
0165             hM_ratio.GetXaxis().SetTitle(XaxisName)
0166             hM_ratio.Draw('PE1X0')
0167         else:
0168             hM_ratio.Draw('PE1X0 same')
0169
0170         l_hM_ratio.append(hM_ratio)
0171
0172     # Draw a line at 1
0173     line = TLine(hM_ratio.GetXaxis().GetXmin(), 1, hM_ratio.GetXaxis().GetXmax(), 1)
0174     line.SetLineColor(kBlack)
0175     line.SetLineStyle(2)
0176     line.Draw()
0177     c.Draw()
0178     c.SaveAs(outname+'.pdf')
0179     c.SaveAs(outname+'.png')
0180     if(c):
0181         c.Close()
0182         gSystem.ProcessEvents()
0183         del c
0184         c = 0
0185
0186
0187
0188 if __name__ == '__main__':
0189     plotpath = './ClusADCCutComp/'
0190     os.makedirs(plotpath, exist_ok=True)
0191
0192     histbasepath = '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists'
0193     datahistspath = ['Data_CombinedNtuple_Run20869_HotDead_BCO_ADC_Survey/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet{}'.format(i) for i in range(3)]
0194     simhistspath = ['Sim_Ntuple_HIJING_new_20240424/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet{}'.format(i) for i in range(3)]
0195
0196     for i, datahistpath in enumerate(datahistspath):
0197         if os.path.isfile("{}/{}/hists_merged.root".format(histbasepath, datahistpath)):
0198             os.system("rm {}/{}/hists_merged.root".format(histbasepath, datahistpath))
0199             os.system("hadd -f -j 20 {}/{}/hists_merged.root {}/{}/hists_*.root".format(histbasepath, datahistpath, histbasepath, datahistpath))
0200         else:
0201             os.system("hadd -f -j 20 {}/{}/hists_merged.root {}/{}/hists_*.root".format(histbasepath, datahistpath, histbasepath, datahistpath))
0202
0203     for i, simhistpath in enumerate(simhistspath):
0204         if os.path.isfile("{}/{}/hists_merged.root".format(histbasepath, simhistpath)):
0205             os.system("rm {}/{}/hists_merged.root".format(histbasepath, simhistpath))
0206             os.system("hadd -f -j 20 {}/{}/hists_merged.root {}/{}/hists_*.root".format(histbasepath, simhistpath, histbasepath, simhistpath))
0207         else:
0208             os.system("hadd -f -j 20 {}/{}/hists_merged.root {}/{}/hists_*.root".format(histbasepath, simhistpath, histbasepath, simhistpath))
0209
0210     l_hM_Eta_reco_Centrality0to70_MBDAsymLe0p75_VtxZm30tom10_data = [GetHistogram(histbasepath+'/'+datahistpath+'/hists_merged.root', 'hM_Eta_reco_Centrality0to70_MBDAsymLe0p75_VtxZm30tom10') for datahistpath in datahistspath]
0211     l_hM_Eta_reco_Centrality0to70_MBDAsymLe0p75_VtxZm30tom10_sim = [GetHistogram(histbasepath+'/'+simhistpath+'/hists_merged.root', 'hM_Eta_reco_Centrality0to70_MBDAsymLe0p75_VtxZm30tom10') for simhistpath in simhistspath]
0212
0213     Draw_1Dhist_datasimcomp(l_hM_Eta_reco_Centrality0to70_MBDAsymLe0p75_VtxZm30tom10_data, l_hM_Eta_reco_Centrality0to70_MBDAsymLe0p75_VtxZm30tom10_sim, [0.05,0.06,0.15,0.32], 'data', False, 2.0, 'Reco-tracklet #eta', '', False, ['Cluster ADC > 35', 'Cluster ADC > 50', 'Cluster ADC > 65'], ['Centrality 0-70%, |Asymm._{MBD}|#leq0.75, -30#leqZ_{vtx}#leq-10[cm]'], plotpath+'Eta_reco_Centrality0to70_MBDAsymLe0p75_VtxZm30tom10_clusAdcCutComp')