File indexing completed on 2025-08-05 08:11:16
0001
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 as np
0009 import math
0010 import glob
0011 from plotUtil import Draw_1Dhist
0012
0013 gROOT.LoadMacro('./sPHENIXStyle/sPhenixStyle.C')
0014 gROOT.ProcessLine('SetsPhenixStyle()')
0015 gROOT.SetBatch(True)
0016
0017 TickSize = 0.03
0018 AxisTitleSize = 0.05
0019 AxisLabelSize = 0.04
0020 LeftMargin = 0.15
0021 RightMargin = 0.08
0022 TopMargin = 0.08
0023 BottomMargin = 0.13
0024
0025 NpercentileDiv = 10
0026 interval = 100/NpercentileDiv
0027
0028 def Draw_1Dhist_wPercentile(hist, l_percentile, norm1, logx, logy, ymaxscale, XaxisName, Ytitle_unit, outname):
0029 hist.Sumw2()
0030 binwidth = hist.GetXaxis().GetBinWidth(1)
0031 binwidth2 = hist.GetXaxis().GetBinWidth(2)
0032 printbinwidth = True
0033 if binwidth != binwidth2:
0034 printbinwidth = False
0035
0036 histmax = hist.GetMaximum()
0037
0038 c = TCanvas('c', 'c', 800, 700)
0039 if norm1:
0040 hist.Scale(1. / hist.Integral(-1, -1))
0041 if logy:
0042 c.SetLogy()
0043 if logx:
0044 c.SetLogx()
0045 c.cd()
0046 gPad.SetRightMargin(RightMargin)
0047 gPad.SetTopMargin(TopMargin)
0048 gPad.SetLeftMargin(LeftMargin)
0049 gPad.SetBottomMargin(BottomMargin)
0050 if printbinwidth:
0051 if norm1:
0052 if Ytitle_unit == '':
0053 hist.GetYaxis().SetTitle('Normalized entries / ({:g})'.format(binwidth))
0054 else:
0055 hist.GetYaxis().SetTitle('Normalized entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0056 else:
0057 if Ytitle_unit == '':
0058 hist.GetYaxis().SetTitle('Entries / ({:g})'.format(binwidth))
0059 else:
0060 hist.GetYaxis().SetTitle('Entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0061 else:
0062 if norm1:
0063 if Ytitle_unit == '':
0064 hist.GetYaxis().SetTitle('Normalized entries')
0065 else:
0066 hist.GetYaxis().SetTitle('Normalized entries {unit})'.format(unit=Ytitle_unit))
0067 else:
0068 if Ytitle_unit == '':
0069 hist.GetYaxis().SetTitle('Entries')
0070 else:
0071 hist.GetYaxis().SetTitle('Entries {unit}'.format(unit=Ytitle_unit))
0072
0073
0074 if logy:
0075 hist.GetYaxis().SetRangeUser(hist.GetMinimum(0)*0.5, (hist.GetMaximum()) * ymaxscale)
0076 else:
0077 hist.GetYaxis().SetRangeUser(0., (hist.GetMaximum()) * ymaxscale)
0078 hist.GetXaxis().SetTitle(XaxisName)
0079 hist.GetXaxis().SetTickSize(TickSize)
0080 hist.GetXaxis().SetTitleSize(AxisTitleSize)
0081 hist.GetXaxis().SetLabelSize(AxisLabelSize)
0082 hist.GetYaxis().SetTickSize(TickSize)
0083 hist.GetYaxis().SetTitleSize(AxisTitleSize)
0084 hist.GetYaxis().SetLabelSize(AxisLabelSize)
0085 hist.GetXaxis().SetTitleOffset(1.1)
0086 hist.GetYaxis().SetTitleOffset(1.35)
0087 hist.SetLineColor(1)
0088 hist.SetLineWidth(2)
0089 hist.Draw('hist')
0090 leg = TLegend(1 - RightMargin - 0.45, 1-TopMargin-0.15, 1 - RightMargin - 0.08, 1-TopMargin-0.035)
0091 leg.SetTextSize(0.040)
0092 leg.SetFillStyle(0)
0093 leg.AddEntry('', '#it{#bf{sPHENIX}} Simulation', '')
0094 leg.AddEntry('','Au+Au #sqrt{s_{NN}}=200 GeV','')
0095 leg.Draw()
0096 linestorage = []
0097 textstorage = []
0098 for i,p in enumerate(l_percentile):
0099
0100
0101 pline = TLine(p, 0, p, histmax*0.6)
0102 pline.SetLineWidth(1)
0103 pline.SetLineStyle(kDashed)
0104 pline.SetLineColor(2)
0105 if len(linestorage)==0:
0106 pline.Draw()
0107 else:
0108 pline.Draw('same')
0109 gPad.Update()
0110 linestorage.append(pline)
0111
0112 if (logx == False and i > 2) or logx == True:
0113 print (i, int(interval*((len(l_percentile))-(i+1))), int(interval*((len(l_percentile))-(i))))
0114 ptext = TText(p, histmax*0.25, '{:d}-{:d}%'.format(int(interval*((len(l_percentile)+1)-(i+2))), int(interval*((len(l_percentile)+1)-(i+1)))))
0115 ptext.SetTextAlign(13)
0116 ptext.SetTextSize(0.02)
0117 ptext.SetTextColor(2)
0118 ptext.SetTextAngle(90)
0119 if len(textstorage)==0:
0120 ptext.Draw()
0121 else:
0122 ptext.Draw('same')
0123 gPad.Update()
0124 textstorage.append(ptext)
0125
0126 c.RedrawAxis()
0127 c.Draw()
0128 c.SaveAs(outname+'.pdf')
0129 c.SaveAs(outname+'.png')
0130 if(c):
0131 c.Close()
0132 gSystem.ProcessEvents()
0133 del c
0134 c = 0
0135
0136 def Draw_2Dhist_wPercentile(hist, centvar, l_percentile, logz, norm1, rmargin, XaxisName, YaxisName, drawopt, outname):
0137 c = TCanvas('c', 'c', 800, 700)
0138 if logz:
0139 c.SetLogz()
0140 c.cd()
0141 gPad.SetRightMargin(rmargin)
0142 gPad.SetTopMargin(TopMargin)
0143 gPad.SetLeftMargin(LeftMargin)
0144 gPad.SetBottomMargin(BottomMargin)
0145 if norm1:
0146 hist.Scale(1. / hist.Integral(-1, -1, -1, -1))
0147 hist.GetXaxis().SetTitle(XaxisName)
0148 hist.GetYaxis().SetTitle(YaxisName)
0149 hist.GetXaxis().SetTickSize(TickSize)
0150 hist.GetYaxis().SetTickSize(TickSize)
0151 hist.GetXaxis().SetTitleSize(AxisTitleSize)
0152 hist.GetYaxis().SetTitleSize(AxisTitleSize)
0153 hist.GetXaxis().SetLabelSize(AxisLabelSize)
0154 hist.GetYaxis().SetLabelSize(AxisLabelSize)
0155 hist.GetXaxis().SetTitleOffset(1.1)
0156 hist.GetYaxis().SetTitleOffset(1.3)
0157 hist.GetZaxis().SetLabelSize(AxisLabelSize)
0158 hist.SetContour(1000)
0159 hist.Draw(drawopt)
0160
0161 linestorage = []
0162 for i,p in enumerate(l_percentile):
0163 if centvar == 'Centrality_mbdquantity':
0164 pline = TLine(0, p, hist.GetXaxis().GetXmax(), p)
0165 elif centvar == 'NClusLayer1':
0166 pline = TLine(p, 0, p, hist.GetYaxis().GetXmax())
0167 pline.SetLineWidth(1)
0168 pline.SetLineStyle(kDashed)
0169 pline.SetLineColor(2)
0170 if len(linestorage)==0:
0171 pline.Draw()
0172 else:
0173 pline.Draw('same')
0174 gPad.Update()
0175 linestorage.append(pline)
0176 gPad.Update()
0177
0178 c.RedrawAxis()
0179 c.Draw()
0180 c.SaveAs(outname+'.pdf')
0181 c.SaveAs(outname+'.png')
0182 if(c):
0183 c.Close()
0184 gSystem.ProcessEvents()
0185 del c
0186 c = 0
0187
0188
0189 if __name__ == '__main__':
0190 parser = OptionParser(usage="usage: %prog ver [options -n]")
0191 parser.add_option("-f", "--inputfile", dest="inputfile", type="string", default='/sphenix/user/hjheng/TrackletAna/minitree/INTT/VtxEvtMap_ana382_zvtx-20cm_dummyAlignParams/INTTVtxZ.root', help="Input ntuple file name")
0192 parser.add_option("-d", "--plotdir", dest="plotdir", type="string", default='./centProxy/ana382_zvtx-20cm_dummyAlignParams', help="Plot directory")
0193 parser.add_option("-c", "--centralityvar", dest="centralityvar", type="string", default='Centrality_mbdquantity', help="Centrality variable name [Centrality_mbdquantity or NClusLayer1]")
0194
0195 (opt, args) = parser.parse_args()
0196
0197 inputfile = opt.inputfile
0198 plotdir = opt.plotdir
0199 centralityvar = opt.centralityvar
0200 maxcent = 50 if centralityvar == 'Centrality_mbdquantity' else 5000
0201
0202 os.makedirs(plotdir, exist_ok=True)
0203
0204 df = ROOT.RDataFrame('minitree', inputfile)
0205 np_centvar = df.AsNumpy(columns=[centralityvar])
0206 CentVar_percentile = []
0207 Binedge_CentVar_percentile = [0]
0208 CentVar_percentile_cut = [0]
0209 for i in range(NpercentileDiv-1):
0210 print('percentile={}-{}%, Centrality quantity({})={}'.format(i*interval, (i+1)*interval, centralityvar, np.percentile(np_centvar[centralityvar], (i+1)*interval)))
0211 CentVar_percentile.append(np.percentile(np_centvar[centralityvar], (i+1)*interval))
0212 Binedge_CentVar_percentile.append(np.percentile(np_centvar[centralityvar], (i+1)*interval))
0213 if i % 2 == 1:
0214 CentVar_percentile_cut.append(np.percentile(np_centvar[centralityvar], (i+1)*interval))
0215 CentVar_percentile_cut.append(maxcent)
0216 Binedge_CentVar_percentile.append(maxcent)
0217
0218 with open('{}/Centrality_{}_bin.txt'.format(plotdir, centralityvar), 'w') as f:
0219 for i in Binedge_CentVar_percentile:
0220 print('{:g}'.format(i), file=f)
0221
0222 hM_NClusLayer1 = TH1F('hM_NClusLayer1', 'hM_NClusLayer1', 200, 0, 4000)
0223 hM_MBDquantity = TH1F('hM_MBDquantity', 'hM_MBDquantity', 200, 0, 25)
0224 hM_NClusLayer1_MBDquantity = TH2F('hM_NClusLayer1_MBDquantity', 'hM_NClusLayer1_MBDquantity', 200, 0, 4000, 200, 0, 25)
0225 f = TFile(inputfile, 'r')
0226 tree = f.Get('minitree')
0227 for idx in range(tree.GetEntries()):
0228 tree.GetEntry(idx)
0229 hM_NClusLayer1.Fill(tree.NClusLayer1)
0230 hM_MBDquantity.Fill(tree.Centrality_mbdquantity)
0231 hM_NClusLayer1_MBDquantity.Fill(tree.NClusLayer1, tree.Centrality_mbdquantity)
0232
0233 if centralityvar == 'Centrality_mbdquantity':
0234 Draw_1Dhist_wPercentile(hM_MBDquantity, CentVar_percentile, False, False, True, 5, 'MBD charge sum (N+S)', '', '{}/MBDquantity_wPercentile'.format(plotdir))
0235 Draw_1Dhist(hM_NClusLayer1, False, False, True, 5, 'Number of clusters (Layer 3+4)', '', '{}/NClusLayer1'.format(plotdir))
0236 elif centralityvar == 'NClusLayer1':
0237 Draw_1Dhist_wPercentile(hM_NClusLayer1, CentVar_percentile, False, False, True, 5, 'Number of clusters (Layer 3+4)', '', '{}/NClusLayer1_wPercentile'.format(plotdir))
0238 Draw_1Dhist(hM_MBDquantity, False, False, True, 5, 'MBD charge sum (N+S)', '', '{}/MBDquantity'.format(plotdir))
0239
0240 Draw_2Dhist_wPercentile(hM_NClusLayer1_MBDquantity, centralityvar, CentVar_percentile, True, False, 0.13, 'Number of clusters (Layer 3+4)', 'MBD charge sum (N+S)', 'colz', '{}/NClusLayer1_MBDquantity'.format(plotdir))
