Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
centProxy.py
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file centProxy.py
1 #! /usr/bin/env python
2 from optparse import OptionParser
3 import sys
4 import os
5 import datetime
6 from array import *
7 from ROOT import *
8 import numpy as np
9 import math
10 import glob
11 from plotUtil import Draw_1Dhist
12 
13 gROOT.LoadMacro('./sPHENIXStyle/sPhenixStyle.C')
14 gROOT.ProcessLine('SetsPhenixStyle()')
15 gROOT.SetBatch(True)
16 
17 TickSize = 0.03
18 AxisTitleSize = 0.05
19 AxisLabelSize = 0.04
20 LeftMargin = 0.15
21 RightMargin = 0.08
22 TopMargin = 0.08
23 BottomMargin = 0.13
24 
25 NpercentileDiv = 10
26 interval = 100/NpercentileDiv
27 
28 def Draw_1Dhist_wPercentile(hist, l_percentile, norm1, logx, logy, ymaxscale, XaxisName, Ytitle_unit, outname):
29  hist.Sumw2()
30  binwidth = hist.GetXaxis().GetBinWidth(1)
31  binwidth2 = hist.GetXaxis().GetBinWidth(2)
32  printbinwidth = True
33  if binwidth != binwidth2:
34  printbinwidth = False
35 
36  histmax = hist.GetMaximum()
37 
38  c = TCanvas('c', 'c', 800, 700)
39  if norm1:
40  hist.Scale(1. / hist.Integral(-1, -1))
41  if logy:
42  c.SetLogy()
43  if logx:
44  c.SetLogx()
45  c.cd()
46  gPad.SetRightMargin(RightMargin)
47  gPad.SetTopMargin(TopMargin)
48  gPad.SetLeftMargin(LeftMargin)
49  gPad.SetBottomMargin(BottomMargin)
50  if printbinwidth:
51  if norm1:
52  if Ytitle_unit == '':
53  hist.GetYaxis().SetTitle('Normalized entries / ({:g})'.format(binwidth))
54  else:
55  hist.GetYaxis().SetTitle('Normalized entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
56  else:
57  if Ytitle_unit == '':
58  hist.GetYaxis().SetTitle('Entries / ({:g})'.format(binwidth))
59  else:
60  hist.GetYaxis().SetTitle('Entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
61  else:
62  if norm1:
63  if Ytitle_unit == '':
64  hist.GetYaxis().SetTitle('Normalized entries')
65  else:
66  hist.GetYaxis().SetTitle('Normalized entries {unit})'.format(unit=Ytitle_unit))
67  else:
68  if Ytitle_unit == '':
69  hist.GetYaxis().SetTitle('Entries')
70  else:
71  hist.GetYaxis().SetTitle('Entries {unit}'.format(unit=Ytitle_unit))
72 
73  # hist.GetXaxis().SetRangeUser(hist.GetBinLowEdge(1)-binwidth, hist.GetBinLowEdge(hist.GetNbinsX())+2*binwidth)
74  if logy:
75  hist.GetYaxis().SetRangeUser(hist.GetMinimum(0)*0.5, (hist.GetMaximum()) * ymaxscale)
76  else:
77  hist.GetYaxis().SetRangeUser(0., (hist.GetMaximum()) * ymaxscale)
78  hist.GetXaxis().SetTitle(XaxisName)
79  hist.GetXaxis().SetTickSize(TickSize)
80  hist.GetXaxis().SetTitleSize(AxisTitleSize)
81  hist.GetXaxis().SetLabelSize(AxisLabelSize)
82  hist.GetYaxis().SetTickSize(TickSize)
83  hist.GetYaxis().SetTitleSize(AxisTitleSize)
84  hist.GetYaxis().SetLabelSize(AxisLabelSize)
85  hist.GetXaxis().SetTitleOffset(1.1)
86  hist.GetYaxis().SetTitleOffset(1.35)
87  hist.SetLineColor(1)
88  hist.SetLineWidth(2)
89  hist.Draw('hist')
90  leg = TLegend(1 - RightMargin - 0.45, 1-TopMargin-0.15, 1 - RightMargin - 0.08, 1-TopMargin-0.035)
91  leg.SetTextSize(0.040)
92  leg.SetFillStyle(0)
93  leg.AddEntry('', '#it{#bf{sPHENIX}} Simulation', '')
94  leg.AddEntry('','Au+Au #sqrt{s_{NN}}=200 GeV','')
95  leg.Draw()
96  linestorage = []
97  textstorage = []
98  for i,p in enumerate(l_percentile):
99  # if (logx == False and i < 3):
100  # continue
101  pline = TLine(p, 0, p, histmax*0.6)
102  pline.SetLineWidth(1)
103  pline.SetLineStyle(kDashed)
104  pline.SetLineColor(2)
105  if len(linestorage)==0:
106  pline.Draw()
107  else:
108  pline.Draw('same')
109  gPad.Update()
110  linestorage.append(pline)
111 
112  if (logx == False and i > 2) or logx == True:
113  print (i, int(interval*((len(l_percentile))-(i+1))), int(interval*((len(l_percentile))-(i))))
114  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)))))
115  ptext.SetTextAlign(13)
116  ptext.SetTextSize(0.02)
117  ptext.SetTextColor(2)
118  ptext.SetTextAngle(90)
119  if len(textstorage)==0:
120  ptext.Draw()
121  else:
122  ptext.Draw('same')
123  gPad.Update()
124  textstorage.append(ptext)
125 
126  c.RedrawAxis()
127  c.Draw()
128  c.SaveAs(outname+'.pdf')
129  c.SaveAs(outname+'.png')
130  if(c):
131  c.Close()
132  gSystem.ProcessEvents()
133  del c
134  c = 0
135 
136 def Draw_2Dhist_wPercentile(hist, centvar, l_percentile, logz, norm1, rmargin, XaxisName, YaxisName, drawopt, outname):
137  c = TCanvas('c', 'c', 800, 700)
138  if logz:
139  c.SetLogz()
140  c.cd()
141  gPad.SetRightMargin(rmargin)
142  gPad.SetTopMargin(TopMargin)
143  gPad.SetLeftMargin(LeftMargin)
144  gPad.SetBottomMargin(BottomMargin)
145  if norm1:
146  hist.Scale(1. / hist.Integral(-1, -1, -1, -1))
147  hist.GetXaxis().SetTitle(XaxisName)
148  hist.GetYaxis().SetTitle(YaxisName)
149  hist.GetXaxis().SetTickSize(TickSize)
150  hist.GetYaxis().SetTickSize(TickSize)
151  hist.GetXaxis().SetTitleSize(AxisTitleSize)
152  hist.GetYaxis().SetTitleSize(AxisTitleSize)
153  hist.GetXaxis().SetLabelSize(AxisLabelSize)
154  hist.GetYaxis().SetLabelSize(AxisLabelSize)
155  hist.GetXaxis().SetTitleOffset(1.1)
156  hist.GetYaxis().SetTitleOffset(1.3)
157  hist.GetZaxis().SetLabelSize(AxisLabelSize)
158  hist.SetContour(1000)
159  hist.Draw(drawopt)
160 
161  linestorage = []
162  for i,p in enumerate(l_percentile):
163  if centvar == 'Centrality_mbdquantity':
164  pline = TLine(0, p, hist.GetXaxis().GetXmax(), p)
165  elif centvar == 'NClusLayer1':
166  pline = TLine(p, 0, p, hist.GetYaxis().GetXmax())
167  pline.SetLineWidth(1)
168  pline.SetLineStyle(kDashed)
169  pline.SetLineColor(2)
170  if len(linestorage)==0:
171  pline.Draw()
172  else:
173  pline.Draw('same')
174  gPad.Update()
175  linestorage.append(pline)
176  gPad.Update()
177 
178  c.RedrawAxis()
179  c.Draw()
180  c.SaveAs(outname+'.pdf')
181  c.SaveAs(outname+'.png')
182  if(c):
183  c.Close()
184  gSystem.ProcessEvents()
185  del c
186  c = 0
187 
188 
189 if __name__ == '__main__':
190  parser = OptionParser(usage="usage: %prog ver [options -n]")
191  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")
192  parser.add_option("-d", "--plotdir", dest="plotdir", type="string", default='./centProxy/ana382_zvtx-20cm_dummyAlignParams', help="Plot directory")
193  parser.add_option("-c", "--centralityvar", dest="centralityvar", type="string", default='Centrality_mbdquantity', help="Centrality variable name [Centrality_mbdquantity or NClusLayer1]")
194 
195  (opt, args) = parser.parse_args()
196 
197  inputfile = opt.inputfile
198  plotdir = opt.plotdir
199  centralityvar = opt.centralityvar
200  maxcent = 50 if centralityvar == 'Centrality_mbdquantity' else 5000
201 
202  os.makedirs(plotdir, exist_ok=True)
203 
204  df = ROOT.RDataFrame('minitree', inputfile)
205  np_centvar = df.AsNumpy(columns=[centralityvar])
206  CentVar_percentile = []
207  Binedge_CentVar_percentile = [0]
208  CentVar_percentile_cut = [0]
209  for i in range(NpercentileDiv-1):
210  print('percentile={}-{}%, Centrality quantity({})={}'.format(i*interval, (i+1)*interval, centralityvar, np.percentile(np_centvar[centralityvar], (i+1)*interval)))
211  CentVar_percentile.append(np.percentile(np_centvar[centralityvar], (i+1)*interval))
212  Binedge_CentVar_percentile.append(np.percentile(np_centvar[centralityvar], (i+1)*interval))
213  if i % 2 == 1:
214  CentVar_percentile_cut.append(np.percentile(np_centvar[centralityvar], (i+1)*interval))
215  CentVar_percentile_cut.append(maxcent)
216  Binedge_CentVar_percentile.append(maxcent)
217 
218  with open('{}/Centrality_{}_bin.txt'.format(plotdir, centralityvar), 'w') as f:
219  for i in Binedge_CentVar_percentile:
220  print('{:g}'.format(i), file=f)
221 
222  hM_NClusLayer1 = TH1F('hM_NClusLayer1', 'hM_NClusLayer1', 200, 0, 4000)
223  hM_MBDquantity = TH1F('hM_MBDquantity', 'hM_MBDquantity', 200, 0, 25)
224  hM_NClusLayer1_MBDquantity = TH2F('hM_NClusLayer1_MBDquantity', 'hM_NClusLayer1_MBDquantity', 200, 0, 4000, 200, 0, 25)
225  f = TFile(inputfile, 'r')
226  tree = f.Get('minitree')
227  for idx in range(tree.GetEntries()):
228  tree.GetEntry(idx)
229  hM_NClusLayer1.Fill(tree.NClusLayer1)
230  hM_MBDquantity.Fill(tree.Centrality_mbdquantity)
231  hM_NClusLayer1_MBDquantity.Fill(tree.NClusLayer1, tree.Centrality_mbdquantity)
232 
233  if centralityvar == 'Centrality_mbdquantity':
234  Draw_1Dhist_wPercentile(hM_MBDquantity, CentVar_percentile, False, False, True, 5, 'MBD charge sum (N+S)', '', '{}/MBDquantity_wPercentile'.format(plotdir))
235  Draw_1Dhist(hM_NClusLayer1, False, False, True, 5, 'Number of clusters (Layer 3+4)', '', '{}/NClusLayer1'.format(plotdir))
236  elif centralityvar == 'NClusLayer1':
237  Draw_1Dhist_wPercentile(hM_NClusLayer1, CentVar_percentile, False, False, True, 5, 'Number of clusters (Layer 3+4)', '', '{}/NClusLayer1_wPercentile'.format(plotdir))
238  Draw_1Dhist(hM_MBDquantity, False, False, True, 5, 'MBD charge sum (N+S)', '', '{}/MBDquantity'.format(plotdir))
239 
240  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))