analysis_INTT/plot/INTTVtxZ_Sim.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 TH1F, TH2F, TCanvas, TFile, TGraphAsymmErrors, TLegend, TColor, gROOT, gStyle, gPad, gSystem, kTRUE, kFALSE, kBird, kThermometer, RDataFrame, TF1
0008 import numpy as np
0009 import math
0010 import glob
0011 from plotUtil import Draw_2Dhist
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
0026 def Draw_1Dhist_fitGaussian(hist, norm1, logy, ymaxscale, XaxisName, Ytitle_unit, outname):
0027     hist.Sumw2()
0028     binwidth = hist.GetXaxis().GetBinWidth(1)
0029     c = TCanvas('c', 'c', 800, 700)
0030     if norm1:
0031         hist.Scale(1. / hist.Integral(-1, -1))
0032     if logy:
0033         c.SetLogy()
0034     c.cd()
0035     gPad.SetRightMargin(RightMargin)
0036     gPad.SetTopMargin(TopMargin)
0037     gPad.SetLeftMargin(LeftMargin)
0038     gPad.SetBottomMargin(BottomMargin)
0039     if norm1:
0040         if Ytitle_unit == '':
0041             hist.GetYaxis().SetTitle(
0042                 'Normalized entries / ({:g})'.format(binwidth))
0043         else:
0044             hist.GetYaxis().SetTitle(
0045                 'Normalized entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0046     else:
0047         if Ytitle_unit == '':
0048             hist.GetYaxis().SetTitle('Entries / ({:g})'.format(binwidth))
0049         else:
0050             hist.GetYaxis().SetTitle(
0051                 'Entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0052
0053     # hist.GetXaxis().SetRangeUser(hist.GetBinLowEdge(1)-binwidth, hist.GetBinLowEdge(hist.GetNbinsX())+2*binwidth)
0054     if logy:
0055         hist.GetYaxis().SetRangeUser(hist.GetMinimum(0)*0.5, (hist.GetMaximum()) * ymaxscale)
0056     else:
0057         hist.GetYaxis().SetRangeUser(0., (hist.GetMaximum()) * ymaxscale)
0058     hist.GetXaxis().SetTitle(XaxisName)
0059     hist.GetXaxis().SetTickSize(TickSize)
0060     hist.GetXaxis().SetTitleSize(AxisTitleSize)
0061     hist.GetXaxis().SetLabelSize(AxisLabelSize)
0062     hist.GetYaxis().SetTickSize(TickSize)
0063     hist.GetYaxis().SetTitleSize(AxisTitleSize)
0064     hist.GetYaxis().SetLabelSize(AxisLabelSize)
0065     hist.GetXaxis().SetTitleOffset(1.1)
0066     hist.GetYaxis().SetTitleOffset(1.35)
0067     hist.SetLineColor(1)
0068     hist.SetLineWidth(2)
0069     hist.Draw('hist')
0070     f1 = TF1('f1', 'gaus', hist.GetXaxis().GetXmin(), hist.GetXaxis().GetXmax())
0071     f1.SetParameter(1,hist.GetMean())
0072     f1.SetParLimits(1,-1,1)
0073     f1.SetParameter(2,hist.GetRMS())
0074     f1.SetParLimits(2,0,100)
0075     f1.SetLineColorAlpha(TColor.GetColor('#F54748'), 0.9)
0076     hist.Fit(f1, 'B')
0077     f1.Draw('same')
0078     leg = TLegend((1-RightMargin)-0.45, (1-TopMargin)-0.13,
0079                   (1-RightMargin)-0.1, (1-TopMargin)-0.03)
0080     leg.SetTextSize(0.040)
0081     leg.SetFillStyle(0)
0082     leg.AddEntry('', '#it{#bf{sPHENIX}} Simulation', '')
0083     leg.AddEntry('','Au+Au #sqrt{s_{NN}}=200 GeV','')
0084     leg.Draw()
0085     leg2 = TLegend(0.54, 0.65, 0.88, 0.78)
0086     leg2.SetTextSize(0.033)
0087     leg2.SetFillStyle(0)
0088     leg2.AddEntry(f1,'Gaussian fit','l')
0089     leg2.AddEntry('', '#mu={0:.2e}#pm{1:.2e}'.format(f1.GetParameter(1), f1.GetParError(1)), '')
0090     leg2.AddEntry('', '#sigma={0:.2e}#pm{1:.2e}'.format(f1.GetParameter(2), f1.GetParError(2)), '')
0091     leg2.Draw()
0092     c.RedrawAxis()
0093     c.Draw()
0094     c.SaveAs(outname+'.pdf')
0095     c.SaveAs(outname+'.png')
0096     if(c):
0097         c.Close()
0098         gSystem.ProcessEvents()
0099         del c
0100         c = 0
0101     return f1.GetParameter(1), f1.GetParError(1), f1.GetParameter(2), f1.GetParError(2)
0102
0103
0104 def Draw_HistGraph(h, xaxistitle, yaxistitle, yaxisrange, xaxisbinlabel, outname):
0105     c = TCanvas('c1', 'c1', 800, 700)
0106     gPad.SetRightMargin(0.08)
0107     c.cd()
0108     h.GetXaxis().SetTitle(xaxistitle)
0109     h.GetYaxis().SetTitle(yaxistitle)
0110     h.SetLineWidth(2)
0111     h.SetMarkerSize(1.5)
0112     h.GetXaxis().SetTitleOffset(1.2)
0113     h.GetXaxis().SetNdivisions(10, 0, 0, kTRUE)
0114     h.GetYaxis().SetTitleOffset(1.5)
0115     h.GetYaxis().SetRangeUser(yaxisrange[0], yaxisrange[1])
0116     # Set bin labels
0117     for i in range(len(xaxisbinlabel)):
0118         h.GetXaxis().SetBinLabel(i+1, xaxisbinlabel[i])
0119
0120     h.Draw('PEX0')
0121
0122     leg = TLegend(1 - RightMargin - 0.45, 1-TopMargin-0.15, 1 - RightMargin - 0.08, 1-TopMargin-0.035)
0123     leg.SetTextSize(0.040)
0124     leg.SetFillStyle(0)
0125     leg.AddEntry('', '#it{#bf{sPHENIX}} Simulation', '')
0126     leg.AddEntry('','Au+Au #sqrt{s_{NN}}=200 GeV','')
0127     leg.Draw()
0128     c.Draw()
0129     c.SaveAs(outname+'.pdf')
0130     c.SaveAs(outname+'.png')
0131     if(c):
0132         c.Close()
0133         gSystem.ProcessEvents()
0134         del c
0135         c = 0
0136
0137 def Draw_1DEff(EffErr, logx, XaxisName, axesrange, outname):
0138     c = TCanvas('c', 'c', 800, 700)
0139     if logx:
0140         c.SetLogx()
0141     c.cd()
0142     gPad.SetRightMargin(RightMargin)
0143     gPad.SetTopMargin(TopMargin)
0144     gPad.SetLeftMargin(LeftMargin)
0145     gPad.SetBottomMargin(BottomMargin)
0146     EffErr.GetXaxis().SetTitle(XaxisName)
0147     EffErr.GetYaxis().SetTitle('Efficiency')
0148     EffErr.GetXaxis().SetRangeUser(axesrange[0], axesrange[1])
0149     # EffErr.GetXaxis().SetMoreLogLabels()
0150     EffErr.GetYaxis().SetRangeUser(axesrange[2], axesrange[3])
0151     EffErr.SetMarkerColor(1)
0152     EffErr.SetMarkerSize(1.5)
0153     EffErr.SetMarkerStyle(20)
0154     EffErr.GetXaxis().SetTitleOffset(1.2)
0155     # EffErr.GetYaxis().SetTitleOffset(ytitleoffset)
0156     # EffErr.GetXaxis().SetNdivisions(10, 0, 0, kTRUE)
0157     EffErr.Draw('AP')
0158     leg = TLegend(LeftMargin+0.03, 1-TopMargin-0.15, LeftMargin+0.1, 1-TopMargin-0.035)
0159     leg.SetTextSize(0.040)
0160     leg.SetFillStyle(0)
0161     leg.AddEntry('', '#it{#bf{sPHENIX}} Simulation', '')
0162     leg.AddEntry('','Au+Au #sqrt{s_{NN}}=200 GeV','')
0163     leg.Draw()
0164     c.Draw()
0165     c.SaveAs(outname+'.pdf')
0166     c.SaveAs(outname+'.png')
0167     if(c):
0168         c.Close()
0169         gSystem.ProcessEvents()
0170         del c
0171         c = 0
0172
0173
0174 def Draw_1DEffComp(leff, lcolor, logx, XaxisName, NdivisionArg, legtext, axesrange, outname):
0175     c = TCanvas('c', 'c', 800, 700)
0176     if logx:
0177         c.SetLogx()
0178     c.cd()
0179     gPad.SetRightMargin(RightMargin)
0180     gPad.SetTopMargin(TopMargin)
0181     gPad.SetLeftMargin(LeftMargin)
0182     gPad.SetBottomMargin(BottomMargin)
0183     for i,eff in enumerate(leff):
0184         if i == 0:
0185             leff[i].GetXaxis().SetTitle(XaxisName)
0186             leff[i].GetXaxis().SetTitleOffset(1.2)
0187             leff[i].GetYaxis().SetTitle('Efficiency')
0188             leff[i].GetXaxis().SetRangeUser(axesrange[0], axesrange[1])
0189             leff[i].GetYaxis().SetRangeUser(axesrange[2], axesrange[3])
0190             leff[i].SetLineColor(TColor.GetColor(lcolor[i]))
0191             leff[i].SetMarkerColor(TColor.GetColor(lcolor[i]))
0192             leff[i].SetMarkerSize(1.3)
0193             leff[i].SetMarkerStyle(20)
0194             if len(NdivisionArg) != 0:
0195                 leff[i].GetXaxis().SetNdivisions(NdivisionArg[0], NdivisionArg[1], NdivisionArg[2], NdivisionArg[3])
0196             leff[i].Draw('AP')
0197         else:
0198             leff[i].SetLineColor(TColor.GetColor(lcolor[i]))
0199             leff[i].SetMarkerColor(TColor.GetColor(lcolor[i]))
0200             leff[i].SetMarkerSize(1.3)
0201             leff[i].SetMarkerStyle(20)
0202             leff[i].Draw('Psame')
0203
0204     leg = TLegend(LeftMargin, (1-TopMargin)-0.15,
0205                   LeftMargin+0.25, (1-TopMargin)-0.03)
0206     leg.SetTextSize(0.04)
0207     leg.SetFillStyle(0)
0208     leg.AddEntry('', '#it{#bf{sPHENIX}} Simulation', '')
0209     leg.AddEntry('','Au+Au #sqrt{s_{NN}}=200 GeV','')
0210     leg.Draw()
0211     leg1 = TLegend((1-gPad.GetRightMargin())-0.47, gPad.GetBottomMargin()+0.04,
0212                    (1-gPad.GetRightMargin())-0.05, gPad.GetBottomMargin()+0.06*len(legtext))
0213     leg1.SetTextSize(0.035)
0214     leg1.SetFillStyle(0)
0215     for i, eff in enumerate(leff):
0216         leg1.AddEntry(eff, legtext[i], "pe")
0217     leg1.Draw()
0218     c.Draw()
0219     c.SaveAs(outname+'.pdf')
0220     c.SaveAs(outname+'.png')
0221     if(c):
0222         c.Close()
0223         gSystem.ProcessEvents()
0224         del c
0225         c = 0
0226
0227 def Draw_2Dhist_eff(hist, logx, logy, logz, norm1, rmargin, XaxisName, YaxisName, drawopt, outname):
0228     # gStyle.SetPalette(kThermometer)
0229     gStyle.SetPaintTextFormat('1.3g')
0230     c = TCanvas('c', 'c', 800, 700)
0231     if logx:
0232         c.SetLogx()
0233     if logy:
0234         c.SetLogy()
0235     if logz:
0236         c.SetLogz()
0237     c.cd()
0238     gPad.SetRightMargin(rmargin)
0239     gPad.SetTopMargin(TopMargin)
0240     gPad.SetLeftMargin(LeftMargin)
0241     gPad.SetBottomMargin(BottomMargin)
0242     if norm1:
0243         hist.Scale(1. / hist.Integral(-1, -1, -1, -1))
0244     hist.GetXaxis().SetTitle(XaxisName)
0245     hist.GetYaxis().SetTitle(YaxisName)
0246     hist.GetZaxis().SetTitle('Efficiency')
0247     hist.GetXaxis().SetTickSize(TickSize)
0248     hist.GetYaxis().SetTickSize(TickSize)
0249     hist.GetXaxis().SetTitleSize(AxisTitleSize)
0250     hist.GetYaxis().SetTitleSize(AxisTitleSize)
0251     hist.GetXaxis().SetLabelSize(AxisLabelSize)
0252     hist.GetYaxis().SetLabelSize(AxisLabelSize)
0253     hist.GetXaxis().SetTitleOffset(1.1)
0254     hist.GetYaxis().SetTitleOffset(1.3)
0255     hist.GetZaxis().SetTitleOffset(1.2)
0256     hist.GetZaxis().SetLabelSize(AxisLabelSize)
0257     hist.GetZaxis().SetRangeUser(0, 1)
0258     hist.SetContour(1000)
0259     hist.SetMarkerSize(0.7)
0260     hist.Draw(drawopt)
0261     c.RedrawAxis()
0262     c.Draw()
0263     c.SaveAs(outname+'.pdf')
0264     c.SaveAs(outname+'.png')
0265     if(c):
0266         c.Close()
0267         gSystem.ProcessEvents()
0268         del c
0269         c = 0
0270
0271
0272
0273 if __name__ == '__main__':
0274     parser = OptionParser(usage="usage: %prog ver [options -h]")
0275     parser.add_option("-i", "--infiledir", dest="infiledir", default='/sphenix/user/hjheng/TrackletAna/minitree/INTT/VtxEvtMap_ana382_zvtx-20cm_dummyAlignParams', help="Input file directory")
0276     parser.add_option("-o", "--outdirprefix", dest="outdirprefix", default='ana382_zvtx-20cm_dummyAlignParams', help="Output file directory")
0277
0278     gStyle.SetPalette(kBird)
0279
0280     (opt, args) = parser.parse_args()
0281     print('opt: {}'.format(opt))
0282
0283     infiledir = opt.infiledir
0284     outdirprefix = opt.outdirprefix
0285
0286     plotpath = './RecoPV_ana/'
0287     os.makedirs('{}/{}'.format(plotpath,outdirprefix), exist_ok=True)
0288
0289     # Input file
0290     fname = '{}/minitree_merged.root'.format(infiledir)
0291     # os.system('hadd -f -j 20 {}/minitree_merged.root {}/minitree_00*.root'.format(infiledir, infiledir))
0292
0293
0294     # of clusters in inner layer, percentile
0295     df = RDataFrame('minitree', fname).Filter('is_min_bias')
0296     nparr_NClusInner = df.AsNumpy(columns=['NClusLayer1'])
0297     NclusInner_percentile = []
0298     Binedge_NclusInner_percentile = [1]
0299     NclusInner_percentile_cut = [0]
0300     NpercentileDiv = 20
0301     for i in range(NpercentileDiv-1):
0302         # print('percentile={}%, Nhits={}'.format((i+1)*5, np.percentile(np_NhitsL1['NClusLayer1'], (i+1)*5)))
0303         NclusInner_percentile.append(np.percentile(nparr_NClusInner['NClusLayer1'], (i+1)*5))
0304         Binedge_NclusInner_percentile.append(np.percentile(nparr_NClusInner['NClusLayer1'], (i+1)*5))
0305         if i % 2 == 1:
0306             NclusInner_percentile_cut.append(np.percentile(nparr_NClusInner['NClusLayer1'], (i+1)*5))
0307     NclusInner_percentile_cut.append(10000)
0308     Binedge_NclusInner_percentile.append(10000)
0309
0310     # Prepare histograms
0311     l_hM_DiffVtxZ = []
0312     for i in range(10):
0313         if i < 7:
0314             l_hM_DiffVtxZ.append(TH1F('hM_DiffVtxZ_Cent{}'.format(5*(2*i+1)), 'hM_DiffVtxZ_Cent{}'.format(5*(2*i+1)), 100, -2, 2))
0315         elif i >= 7:
0316             l_hM_DiffVtxZ.append(TH1F('hM_DiffVtxZ_Cent{}'.format(5*(2*i+1)), 'hM_DiffVtxZ_Cent{}'.format(5*(2*i+1)), 100, -5, 5))
0317         else:
0318             print('ERROR: i={}'.format(i))
0319             sys.exit()
0320
0321     hM_NClusInner_NtruthVtx1 = TH1F('hM_NClusInner_NtruthVtx1', 'hM_NClusInner_NtruthVtx1', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile))
0322     hM_NClusInner_NtruthVtx1_absdiffle120cm = TH1F('hM_NClusInner_NtruthVtx1_absdiffle120cm', 'hM_NClusInner_NtruthVtx1_absdiffle120cm', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile)) # very loose cut
0323     hM_NClusInner_NtruthVtx1_absdiffle5cm = TH1F('hM_NClusInner_NtruthVtx1_absdiffle5cm', 'hM_NClusInner_NtruthVtx1_absdiffle5cm', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile))
0324     hM_NClusInner_NtruthVtx1_absdiffle2cm = TH1F('hM_NClusInner_NtruthVtx1_absdiffle2cm', 'hM_NClusInner_NtruthVtx1_absdiffle2cm', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile))
0325     hM_NClusInner_NtruthVtx1_absdiffle1cm = TH1F('hM_NClusInner_NtruthVtx1_absdiffle1cm', 'hM_NClusInner_NtruthVtx1_absdiffle1cm', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile))
0326     hM_Centrality_NtruthVtx1 = TH1F('hM_Centrality_NtruthVtx1', 'hM_Centrality_NtruthVtx1', 10, 0, 100)
0327     hM_Centrality_NtruthVtx1_absdiffle120cm = TH1F('hM_Centrality_NtruthVtx1_absdiffle120cm', 'hM_Centrality_NtruthVtx1_absdiffle120cm', 10, 0, 100) # very loose cut
0328     hM_Centrality_NtruthVtx1_absdiffle5cm = TH1F('hM_Centrality_NtruthVtx1_absdiffle5cm', 'hM_Centrality_NtruthVtx1_absdiffle5cm', 10, 0, 100)
0329     hM_Centrality_NtruthVtx1_absdiffle2cm = TH1F('hM_Centrality_NtruthVtx1_absdiffle2cm', 'hM_Centrality_NtruthVtx1_absdiffle2cm', 10, 0, 100)
0330     hM_Centrality_NtruthVtx1_absdiffle1cm = TH1F('hM_Centrality_NtruthVtx1_absdiffle1cm', 'hM_Centrality_NtruthVtx1_absdiffle1cm', 10, 0, 100)
0331     hM_NClusInner_TruthPVz_NtruthVtx1 = TH2F('hM_NClusInner_TruthPVz_NtruthVtx1', 'hM_NClusInner_TruthPVz_NtruthVtx1', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile), 15, -30, 30)
0332     hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle120cm = TH2F('hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle120cm', 'hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle120cm', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile), 15, -30, 30) # very loose cut
0333     hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle5cm = TH2F('hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle5cm', 'hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle5cm', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile), 15, -30, 30)
0334     hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle2cm = TH2F('hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle2cm', 'hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle2cm', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile), 15, -30, 30)
0335     hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle1cm = TH2F('hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle1cm', 'hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle1cm', NpercentileDiv, np.asarray(Binedge_NclusInner_percentile), 15, -30, 30)
0336     hM_Centrality_TruthPVz_NtruthVtx1 = TH2F('hM_Centrality_TruthPVz_NtruthVtx1', 'hM_Centrality_TruthPVz_NtruthVtx1', 10, 0, 100, 15, -30, 30)
0337     hM_Centrality_TruthPVz_NtruthVtx1_absdiffle120cm = TH2F('hM_Centrality_TruthPVz_NtruthVtx1_absdiffle120cm', 'hM_Centrality_TruthPVz_NtruthVtx1_absdiffle120cm', 10, 0, 100, 15, -30, 30) # very loose cut
0338     hM_Centrality_TruthPVz_NtruthVtx1_absdiffle5cm = TH2F('hM_Centrality_TruthPVz_NtruthVtx1_absdiffle5cm', 'hM_Centrality_TruthPVz_NtruthVtx1_absdiffle5cm', 10, 0, 100, 15, -30, 30)
0339     hM_Centrality_TruthPVz_NtruthVtx1_absdiffle2cm = TH2F('hM_Centrality_TruthPVz_NtruthVtx1_absdiffle2cm', 'hM_Centrality_TruthPVz_NtruthVtx1_absdiffle2cm', 10, 0, 100, 15, -30, 30)
0340     hM_Centrality_TruthPVz_NtruthVtx1_absdiffle1cm = TH2F('hM_Centrality_TruthPVz_NtruthVtx1_absdiffle1cm', 'hM_Centrality_TruthPVz_NtruthVtx1_absdiffle1cm', 10, 0, 100, 15, -30, 30)
0341
0342     # Loop events and fill histograms
0343     f = TFile(fname, 'r')
0344     tree = f.Get('minitree')
0345     for idx in range(tree.GetEntries()):
0346         tree.GetEntry(idx)
0347         # skip MBD_centrality = nan
0348         if math.isnan(tree.MBD_centrality):
0349             continue
0350
0351         if tree.is_min_bias == False:
0352             continue
0353
0354         if tree.NTruthVtx == 1:
0355             idx = int(((tree.MBD_centrality / 5) - 1) / 2)
0356             l_hM_DiffVtxZ[idx].Fill(tree.PV_z - tree.TruthPV_z)
0357             hM_NClusInner_NtruthVtx1.Fill(tree.NClusLayer1)
0358             if abs(tree.PV_z - tree.TruthPV_z) <= 120:
0359                 hM_NClusInner_NtruthVtx1_absdiffle120cm.Fill(tree.NClusLayer1)
0360             if abs(tree.PV_z - tree.TruthPV_z) <= 5:
0361                 hM_NClusInner_NtruthVtx1_absdiffle5cm.Fill(tree.NClusLayer1)
0362             if abs(tree.PV_z - tree.TruthPV_z) <= 2:
0363                 hM_NClusInner_NtruthVtx1_absdiffle2cm.Fill(tree.NClusLayer1)
0364             if abs(tree.PV_z - tree.TruthPV_z) <= 1:
0365                 hM_NClusInner_NtruthVtx1_absdiffle1cm.Fill(tree.NClusLayer1)
0366
0367             # centbin = int(((tree.MBD_centrality / 5) - 1) / 2)
0368             hM_Centrality_NtruthVtx1.Fill(tree.MBD_centrality)
0369             if abs(tree.PV_z - tree.TruthPV_z) <= 120:
0370                 hM_Centrality_NtruthVtx1_absdiffle120cm.Fill(tree.MBD_centrality)
0371             if abs(tree.PV_z - tree.TruthPV_z) <= 5:
0372                 hM_Centrality_NtruthVtx1_absdiffle5cm.Fill(tree.MBD_centrality)
0373             if abs(tree.PV_z - tree.TruthPV_z) <= 2:
0374                 hM_Centrality_NtruthVtx1_absdiffle2cm.Fill(tree.MBD_centrality)
0375             if abs(tree.PV_z - tree.TruthPV_z) <= 1:
0376                 hM_Centrality_NtruthVtx1_absdiffle1cm.Fill(tree.MBD_centrality)
0377
0378             hM_NClusInner_TruthPVz_NtruthVtx1.Fill(tree.NClusLayer1, tree.TruthPV_z)
0379             if abs(tree.PV_z - tree.TruthPV_z) <= 120:
0380                 hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle120cm.Fill(tree.NClusLayer1, tree.TruthPV_z)
0381             if abs(tree.PV_z - tree.TruthPV_z) <= 5:
0382                 hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle5cm.Fill(tree.NClusLayer1, tree.TruthPV_z)
0383             if abs(tree.PV_z - tree.TruthPV_z) <= 2:
0384                 hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle2cm.Fill(tree.NClusLayer1, tree.TruthPV_z)
0385             if abs(tree.PV_z - tree.TruthPV_z) <= 1:
0386                 hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle1cm.Fill(tree.NClusLayer1, tree.TruthPV_z)
0387
0388             hM_Centrality_TruthPVz_NtruthVtx1.Fill(tree.MBD_centrality, tree.TruthPV_z)
0389             if abs(tree.PV_z - tree.TruthPV_z) <= 120:
0390                 hM_Centrality_TruthPVz_NtruthVtx1_absdiffle120cm.Fill(tree.MBD_centrality, tree.TruthPV_z)
0391             if abs(tree.PV_z - tree.TruthPV_z) <= 5:
0392                 hM_Centrality_TruthPVz_NtruthVtx1_absdiffle5cm.Fill(tree.MBD_centrality, tree.TruthPV_z)
0393             if abs(tree.PV_z - tree.TruthPV_z) <= 2:
0394                 hM_Centrality_TruthPVz_NtruthVtx1_absdiffle2cm.Fill(tree.MBD_centrality, tree.TruthPV_z)
0395             if abs(tree.PV_z - tree.TruthPV_z) <= 1:
0396                 hM_Centrality_TruthPVz_NtruthVtx1_absdiffle1cm.Fill(tree.MBD_centrality, tree.TruthPV_z)
0397
0398
0399     l_res_vtxZ = []
0400     l_reserr_vtxZ = []
0401     h_resolution_cent = TH1F('h_resolution_cent', 'h_resolution_cent', 10, 0, 10)
0402     h_bias_cent = TH1F('h_bias_cent', 'h_bias_cent', 10, 0, 10)
0403     for i, hist in enumerate(l_hM_DiffVtxZ):
0404         mean, meanerr, sigma, sigmaerr = Draw_1Dhist_fitGaussian(hist, False, True, 50, '#Deltaz(vtx_{Truth}, vtx_{Reco}) [cm]', 'cm', '{}/{}/DiffVtxZ_Cent{}'.format(plotpath,outdirprefix,5*(2*i+1)))
0405         l_res_vtxZ.append(sigma)
0406         l_reserr_vtxZ.append(sigmaerr)
0407         h_bias_cent.SetBinContent(i+1, mean)
0408         h_bias_cent.SetBinError(i+1, meanerr)
0409         h_resolution_cent.SetBinContent(i+1, sigma)
0410         h_resolution_cent.SetBinError(i+1, sigmaerr)
0411
0412     Draw_HistGraph(h_bias_cent, 'Centrality', '#mu_{#Deltaz}^{Gaussian} [cm]', [-0.1,0.1], ['{}-{}%'.format(10*i, 10*(i+1)) for i in range(len(l_reserr_vtxZ))], '{}/{}/VtxZBias_Cent'.format(plotpath,outdirprefix))
0413     Draw_HistGraph(h_resolution_cent, 'Centrality', '#sigma_{#Deltaz}^{Gaussian} [cm]', [0, h_resolution_cent.GetMaximum()*1.5], ['{}-{}%'.format(10*i, 10*(i+1)) for i in range(len(l_reserr_vtxZ))], '{}/{}/VtxZReolution_Cent'.format(plotpath,outdirprefix))
0414
0415     # Efficiency
0416     err_RecoPVEff_NClusInner_absdiffle60cm = TGraphAsymmErrors()
0417     err_RecoPVEff_NClusInner_absdiffle60cm.Divide(hM_NClusInner_NtruthVtx1_absdiffle120cm, hM_NClusInner_NtruthVtx1)
0418     err_RecoPVEff_NClusInner_absdiffle5cm = TGraphAsymmErrors()
0419     err_RecoPVEff_NClusInner_absdiffle5cm.Divide(hM_NClusInner_NtruthVtx1_absdiffle5cm, hM_NClusInner_NtruthVtx1)
0420     err_RecoPVEff_NClusInner_absdiffle2cm = TGraphAsymmErrors()
0421     err_RecoPVEff_NClusInner_absdiffle2cm.Divide(hM_NClusInner_NtruthVtx1_absdiffle2cm, hM_NClusInner_NtruthVtx1)
0422     err_RecoPVEff_NClusInner_absdiffle1cm = TGraphAsymmErrors()
0423     err_RecoPVEff_NClusInner_absdiffle1cm.Divide(hM_NClusInner_NtruthVtx1_absdiffle1cm, hM_NClusInner_NtruthVtx1)
0424     list_eff_NClusInner = [err_RecoPVEff_NClusInner_absdiffle60cm, err_RecoPVEff_NClusInner_absdiffle5cm, err_RecoPVEff_NClusInner_absdiffle2cm, err_RecoPVEff_NClusInner_absdiffle1cm]
0425     list_color = ['#03001C', '#9B0000', '#035397', '#347928']
0426     list_leg = ['|#Deltaz(vtx_{Reco},vtx_{Truth})|#leq120cm', '|#Deltaz(vtx_{Reco},vtx_{Truth})|#leq5cm', '|#Deltaz(vtx_{Reco},vtx_{Truth})|#leq2cm', '|#Deltaz(vtx_{Reco},vtx_{Truth})|#leq1cm']
0427     Draw_1DEffComp(list_eff_NClusInner, list_color, True, 'Number of clusters (inner)', [], list_leg, [0,10100,0,1.4], '{}/{}/RecoPVEff_NClusInner_EffComp'.format(plotpath,outdirprefix))
0428
0429     err_RecoPVEff_Centrality_absdiffle120cm = TGraphAsymmErrors()
0430     err_RecoPVEff_Centrality_absdiffle120cm.Divide(hM_Centrality_NtruthVtx1_absdiffle120cm, hM_Centrality_NtruthVtx1)
0431     # err_RecoPVEff_Centrality_absdiffle120cm.GetXaxis().SetMaxDigits(2)
0432     err_RecoPVEff_Centrality_absdiffle5cm = TGraphAsymmErrors()
0433     err_RecoPVEff_Centrality_absdiffle5cm.Divide(hM_Centrality_NtruthVtx1_absdiffle5cm, hM_Centrality_NtruthVtx1)
0434     # err_RecoPVEff_Centrality_absdiffle5cm.GetXaxis().SetMaxDigits(2)
0435     err_RecoPVEff_Centrality_absdiffle2cm = TGraphAsymmErrors()
0436     err_RecoPVEff_Centrality_absdiffle2cm.Divide(hM_Centrality_NtruthVtx1_absdiffle2cm, hM_Centrality_NtruthVtx1)
0437     err_RecoPVEff_Centrality_absdiffle1cm = TGraphAsymmErrors()
0438     err_RecoPVEff_Centrality_absdiffle1cm.Divide(hM_Centrality_NtruthVtx1_absdiffle1cm, hM_Centrality_NtruthVtx1)
0439     list_eff_Centrality = [err_RecoPVEff_Centrality_absdiffle120cm, err_RecoPVEff_Centrality_absdiffle5cm, err_RecoPVEff_Centrality_absdiffle2cm, err_RecoPVEff_Centrality_absdiffle1cm]
0440     Draw_1DEffComp(list_eff_Centrality, list_color, False, 'Centrality', [10,0,0,kTRUE], list_leg, [0,100,0,1.4], '{}/{}/RecoPVEff_CentralityBin_EffComp'.format(plotpath,outdirprefix))
0441
0442     # 2D efficiency
0443     eff_RecoPVEff_NClusInner_TruthPVz_absdiffle120cm = hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle120cm.Clone()
0444     eff_RecoPVEff_NClusInner_TruthPVz_absdiffle120cm.Divide(hM_NClusInner_TruthPVz_NtruthVtx1)
0445     eff_RecoPVEff_NClusInner_TruthPVz_absdiffle5cm = hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle5cm.Clone()
0446     eff_RecoPVEff_NClusInner_TruthPVz_absdiffle5cm.Divide(hM_NClusInner_TruthPVz_NtruthVtx1)
0447     eff_RecoPVEff_NClusInner_TruthPVz_absdiffle2cm = hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle2cm.Clone()
0448     eff_RecoPVEff_NClusInner_TruthPVz_absdiffle2cm.Divide(hM_NClusInner_TruthPVz_NtruthVtx1)
0449     eff_RecoPVEff_NClusInner_TruthPVz_absdiffle1cm = hM_NClusInner_TruthPVz_NtruthVtx1_absdiffle1cm.Clone()
0450     eff_RecoPVEff_NClusInner_TruthPVz_absdiffle1cm.Divide(hM_NClusInner_TruthPVz_NtruthVtx1)
0451     # Draw_2Dhist_eff(hist, logx, logy, logz, norm1, rmargin, XaxisName, YaxisName, drawopt, outname)
0452     Draw_2Dhist_eff(eff_RecoPVEff_NClusInner_TruthPVz_absdiffle120cm, True, False, False, False, 0.18, 'Number of clusters (inner)', 'Truth vtx_{z} [cm]', 'colz text45', '{}/{}/RecoPVEff_NClusInner_TruthPVz_absdiffle120cm'.format(plotpath,outdirprefix))
0453     Draw_2Dhist_eff(eff_RecoPVEff_NClusInner_TruthPVz_absdiffle5cm, True, False, False, False, 0.18, 'Number of clusters (inner)', 'Truth vtx_{z} [cm]', 'colz text45', '{}/{}/RecoPVEff_NClusInner_TruthPVz_absdiffle5cm'.format(plotpath,outdirprefix))
0454     Draw_2Dhist_eff(eff_RecoPVEff_NClusInner_TruthPVz_absdiffle2cm, True, False, False, False, 0.18, 'Number of clusters (inner)', 'Truth vtx_{z} [cm]', 'colz text45', '{}/{}/RecoPVEff_NClusInner_TruthPVz_absdiffle2cm'.format(plotpath,outdirprefix))
0455     Draw_2Dhist_eff(eff_RecoPVEff_NClusInner_TruthPVz_absdiffle1cm, True, False, False, False, 0.18, 'Number of clusters (inner)', 'Truth vtx_{z} [cm]', 'colz text45', '{}/{}/RecoPVEff_NClusInner_TruthPVz_absdiffle1cm'.format(plotpath,outdirprefix))
0456
0457     eff_RecoPVEff_Centrality_TruthPVz_absdiffle120cm = hM_Centrality_TruthPVz_NtruthVtx1_absdiffle120cm.Clone()
0458     eff_RecoPVEff_Centrality_TruthPVz_absdiffle120cm.Divide(hM_Centrality_TruthPVz_NtruthVtx1)
0459     eff_RecoPVEff_Centrality_TruthPVz_absdiffle120cm.GetXaxis().SetNdivisions(10, 0, 0, kTRUE)
0460     eff_RecoPVEff_Centrality_TruthPVz_absdiffle5cm = hM_Centrality_TruthPVz_NtruthVtx1_absdiffle5cm.Clone()
0461     eff_RecoPVEff_Centrality_TruthPVz_absdiffle5cm.Divide(hM_Centrality_TruthPVz_NtruthVtx1)
0462     eff_RecoPVEff_Centrality_TruthPVz_absdiffle5cm.GetXaxis().SetNdivisions(10, 0, 0, kTRUE)
0463     eff_RecoPVEff_Centrality_TruthPVz_absdiffle2cm = hM_Centrality_TruthPVz_NtruthVtx1_absdiffle2cm.Clone()
0464     eff_RecoPVEff_Centrality_TruthPVz_absdiffle2cm.Divide(hM_Centrality_TruthPVz_NtruthVtx1)
0465     eff_RecoPVEff_Centrality_TruthPVz_absdiffle2cm.GetXaxis().SetNdivisions(10, 0, 0, kTRUE)
0466     eff_RecoPVEff_Centrality_TruthPVz_absdiffle1cm = hM_Centrality_TruthPVz_NtruthVtx1_absdiffle1cm.Clone()
0467     eff_RecoPVEff_Centrality_TruthPVz_absdiffle1cm.Divide(hM_Centrality_TruthPVz_NtruthVtx1)
0468     eff_RecoPVEff_Centrality_TruthPVz_absdiffle1cm.GetXaxis().SetNdivisions(10, 0, 0, kTRUE)
0469     Draw_2Dhist_eff(eff_RecoPVEff_Centrality_TruthPVz_absdiffle120cm, False, False, False, False, 0.18, 'Centrality', 'Truth vtx_{z} [cm]', 'colz text45', '{}/{}/RecoPVEff_Centrality_TruthPVz_absdiffle120cm'.format(plotpath,outdirprefix))
0470     Draw_2Dhist_eff(eff_RecoPVEff_Centrality_TruthPVz_absdiffle5cm, False, False, False, False, 0.18, 'Centrality', 'Truth vtx_{z} [cm]', 'colz text45', '{}/{}/RecoPVEff_Centrality_TruthPVz_absdiffle5cm'.format(plotpath,outdirprefix))
0471     Draw_2Dhist_eff(eff_RecoPVEff_Centrality_TruthPVz_absdiffle2cm, False, False, False, False, 0.18, 'Centrality', 'Truth vtx_{z} [cm]', 'colz text45', '{}/{}/RecoPVEff_Centrality_TruthPVz_absdiffle2cm'.format(plotpath,outdirprefix))
0472     Draw_2Dhist_eff(eff_RecoPVEff_Centrality_TruthPVz_absdiffle1cm, False, False, False, False, 0.18, 'Centrality', 'Truth vtx_{z} [cm]', 'colz text45', '{}/{}/RecoPVEff_Centrality_TruthPVz_absdiffle1cm'.format(plotpath,outdirprefix))
0473