analysis_INTT/plot/EtaDepClusADCCut.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 from plot_DataSimComp import Draw_1Dhist_datasimcomp
0013
0014 gROOT.SetBatch(True)
0015
0016 smin = 5
0017 smax = 55
0018 sstep = 1
0019
0020 esmin = 0.6
0021 esmax = 1.2
0022 esstep = 0.02
0023
0024 def theta2pseudorapidity(theta):
0025     theta = theta * TMath.Pi() / 180.0
0026     return -1.0 * TMath.Log(TMath.Tan(theta/2.0))
0027
0028
0029 def EtaDepADCCut_INTTPrivate():
0030     thetastep = [0.001, 15, 20, 25, 30, 35, 45, 55, 125, 135, 145, 150, 155, 160, 165, 179.999]
0031     # reverse thetastep
0032     thetastep = thetastep[::-1]
0033     adccut_theta = [225, 165, 135, 120, 105, 90, 75, 60, 75, 90, 105, 120, 135, 165, 225]
0034     etastep = []
0035     for i in range(0, len(thetastep)):
0036         etastep.append(theta2pseudorapidity(thetastep[i]))
0037
0038     hm = TH1F('hm', 'hm', len(etastep)-1, array('d', etastep))
0039     for j in range(hm.GetNbinsX()):
0040         hm.SetBinContent(j+1, adccut_theta[j])
0041
0042     print (etastep)
0043
0044     return hm
0045
0046
0047 def funclist(scalemin, scalemax, scalestep):
0048     l_func = []
0049     for i in range(scalemin, scalemax, scalestep):
0050         f = TF1('f', '{}*TMath::CosH(x)'.format(i), -5, 5)
0051         l_func.append(f)
0052     return l_func
0053
0054
0055 def Draw_2Dhist_wFunc(hist_wcut, hist_wocut, IsData, logz, norm1, rmargin, XaxisName, YaxisName, ZaxisName, drawopt_wcut, drawopt_wocut, outname):
0056     # fs = funclist(smin, smax, sstep)
0057     # nColors = gStyle.GetNumberOfColors()
0058     # nFuncs = len(fs)
0059
0060     # strip the outname to get the const scale: it should be at the end of the string, remove 'constscale'
0061     scale = -999
0062     etascale = -999
0063     if 'inttprivate' not in outname:
0064         scale = outname.split('_')[-2].replace('constscale','')
0065         etascale = outname.split('_')[-1].replace('etascale','').replace('p','.')
0066         print ('scale = {}, etascale = {}'.format(scale, etascale))
0067
0068     c = TCanvas('c', 'c', 800, 700)
0069     if logz:
0070         c.SetLogz()
0071     c.cd()
0072     if ZaxisName == '':
0073         gPad.SetRightMargin(rmargin)
0074     else:
0075         gPad.SetRightMargin(rmargin+0.03)
0076
0077     gPad.SetTopMargin(TopMargin)
0078     gPad.SetLeftMargin(LeftMargin)
0079     gPad.SetBottomMargin(BottomMargin)
0080     if norm1:
0081         hist_wcut.Scale(1. / hist_wcut.Integral(-1, -1, -1, -1))
0082     hist_wcut.GetXaxis().SetTitle(XaxisName)
0083     hist_wcut.GetYaxis().SetTitle(YaxisName)
0084     hist_wcut.GetXaxis().SetTickSize(TickSize)
0085     hist_wcut.GetYaxis().SetTickSize(TickSize)
0086     hist_wcut.GetXaxis().SetTitleSize(AxisTitleSize)
0087     hist_wcut.GetYaxis().SetTitleSize(AxisTitleSize)
0088     hist_wcut.GetXaxis().SetLabelSize(AxisLabelSize)
0089     hist_wcut.GetYaxis().SetLabelSize(AxisLabelSize)
0090     if ZaxisName != '':
0091         hist_wcut.GetZaxis().SetTitle(ZaxisName)
0092         hist_wcut.GetZaxis().SetTitleSize(AxisTitleSize)
0093         hist_wcut.GetZaxis().SetTitleOffset(1.1)
0094
0095     hist_wcut.GetXaxis().SetTitleOffset(1.1)
0096     hist_wcut.GetYaxis().SetTitleOffset(1.3)
0097     hist_wcut.GetZaxis().SetLabelSize(AxisLabelSize)
0098     hist_wcut.SetContour(1000)
0099     hist_wcut.Draw(drawopt_wcut)
0100
0101     hist_wocut.SetLineWidth(1)
0102     hist_wocut.Draw(drawopt_wocut)
0103
0104     if 'inttprivate' not in outname:
0105         adcstep = [20+i*30 for i in range(0, 20)] # quantized step
0106         etastep = []
0107         f = TF1('f', '{}*TMath::CosH({}*x)'.format(scale, etascale), -10, 10)
0108         for i in range(0, 20):
0109             etastep.insert(0, f.GetX(adcstep[i], -10, 0))
0110             etastep.append(f.GetX(adcstep[i], 0, 10))
0111
0112         # remove nan values
0113         etastep = [x for x in etastep if not math.isnan(x)]
0114         print(etastep)
0115
0116         hm_cut_v2 = TH1F('hm_cut_v2', 'hm_cut_v2', len(etastep)-1, array('d', etastep))
0117         for j in range(hm_cut_v2.GetNbinsX()):
0118             # Set bin content at the quantized step
0119             if hm_cut_v2.GetBinLowEdge(j+1) < 0:
0120                 hm_cut_v2.SetBinContent(j+1, f.Eval(hm_cut_v2.GetBinLowEdge(j+1)))
0121             else:
0122                 hm_cut_v2.SetBinContent(j+1, f.Eval(hm_cut_v2.GetBinCenter(j+1) + hm_cut_v2.GetBinWidth(j+1)/2))
0123             # hm_cut_v2.SetBinContent(j+1, f.Eval(hm_cut_v2.GetBinCenter(j+1)))
0124
0125         hm_cut_v2.SetLineColorAlpha(kRed, 1)
0126         hm_cut_v2.SetLineWidth(2)
0127         hm_cut_v2.Draw('hist same')
0128
0129         f.SetLineColorAlpha(kRed, 1)
0130         f.SetLineStyle(2)
0131         f.SetLineWidth(2)
0132         f.Draw('same')
0133     else:
0134         hm_cut = EtaDepADCCut_INTTPrivate()
0135         hm_cut.SetLineColorAlpha(kRed, 1)
0136         hm_cut.SetLineWidth(2)
0137         hm_cut.Draw('hist same')
0138
0139     # hm_cut = TH1F('hm_cut', 'hm_cut', 40, -10, 10)
0140     # for j in range(hm_cut.GetNbinsX()):
0141     #     hm_cut.SetBinContent(j+1, f.Eval(hm_cut.GetBinCenter(j+1)))
0142     # hm_cut.SetLineColorAlpha(kRed, 1)
0143     # hm_cut.SetLineWidth(2)
0144     # hm_cut.Draw('hist same')
0145
0146     # lhm_cut = []
0147     # for i, func in enumerate(fs):
0148     #     hm_cut = TH1F('hm_cut_{}'.format(i), 'hm_cut_{}'.format(i), 40, -10, 10)
0149     #     for j in range(hm_cut.GetNbinsX()):
0150     #         hm_cut.SetBinContent(j+1, func.Eval(hm_cut.GetBinCenter(j+1)))
0151
0152     #     lhm_cut.append(hm_cut)
0153     #     fcolor = int(float(nColors) / nFuncs * i)
0154     #     # func.SetLineColorAlpha(gStyle.GetColorPalette(fcolor), 0.5)
0155     #     # func.SetLineWidth(2)
0156     #     # func.Draw('same')
0157     #     lhm_cut[i].SetLineColorAlpha(gStyle.GetColorPalette(fcolor), 0.5)
0158     #     lhm_cut[i].SetLineWidth(2)
0159     #     lhm_cut[i].Draw('hist same')
0160
0161     # Add text for the function
0162     tex = TLatex()
0163     tex.SetTextSize(0.04)
0164     tex.SetTextAlign(11)
0165     tex.SetTextFont(42)
0166     tex.SetTextColor(kRed)
0167     tex.SetNDC()
0168     if 'inttprivate' not in outname:
0169         tex.DrawLatex(LeftMargin+0.05, 0.85, 'ADC cut = {:d}#timescosh({:.2f}#times#eta)'.format(int(scale), float(etascale)))
0170     else:
0171         tex.DrawLatex(LeftMargin+0.05, 0.85, 'INTT private #eta-dependent ADC cut')
0172
0173
0174     rightshift = 0.0 if IsData else 0.0
0175     leg = TLegend((1-RightMargin)-0.45, (1-TopMargin)+0.01, (1-RightMargin)-rightshift, (1-TopMargin)+0.04)
0176     leg.SetTextAlign(kHAlignRight+kVAlignBottom)
0177     leg.SetTextSize(0.045)
0178     leg.SetFillStyle(0)
0179     if IsData:
0180         leg.AddEntry("", "#it{#bf{sPHENIX}} Work-in-progress", "")
0181         # leg.AddEntry("", "Au+Au #sqrt{s_{NN}}=200 GeV", "")
0182     else:
0183         leg.AddEntry("", "#it{#bf{sPHENIX}} Simulation", "")
0184         # leg.AddEntry("", "Au+Au #sqrt{s_{NN}}=200 GeV", "")
0185     leg.Draw()
0186     c.RedrawAxis()
0187     c.Draw()
0188     c.SaveAs(outname+'.pdf')
0189     c.SaveAs(outname+'.png')
0190     if(c):
0191         c.Close()
0192         gSystem.ProcessEvents()
0193         del c
0194         c = 0
0195
0196
0197 if __name__ == '__main__':
0198     parser = OptionParser()
0199     parser.add_option("-d", "--plotdir", dest="plotdir", type="string", default='/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/RecoCluster', help="Plot directory")
0200
0201     (opt, args) = parser.parse_args()
0202     print('opt: {}'.format(opt))
0203
0204     plotdir = opt.plotdir
0205     os.makedirs(plotdir, exist_ok=True)
0206     os.makedirs(plotdir+'/ClusADCCutOptim', exist_ok=True)
0207     os.makedirs(plotdir+'/ClusADCCutOptim/Hist1D_ClusEtaPV', exist_ok=True)
0208     os.makedirs(plotdir+'/ClusADCCutOptim/Hist2D_ClusEtaPVADC_wStepFunc', exist_ok=True)
0209
0210     datahistfile = '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_CombinedNtuple_Run20869_HotDead_BCO_ADC_Survey/Cluster/hists_merged.root'
0211     simhistfile = '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_Ntuple_HIJING_new_20240424/Cluster/hists_merged.root'
0212
0213     if not os.path.isfile(datahistfile) or not os.path.isfile(simhistfile):
0214         print("Error: file for the merged histograms not found")
0215         sys.exit(1)
0216
0217     EtaDepADCCut_INTTPrivate()
0218
0219     hM_ClusEtaPV_ClusADC_all_zoomin_data = GetHistogram(datahistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin')
0220     hM_ClusEtaPV_ClusADC_all_zoomin_sim = GetHistogram(simhistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin')
0221     NClus_total_data = hM_ClusEtaPV_ClusADC_all_zoomin_data.Integral(-1,-1,-1,-1)
0222     NClus_total_sim = hM_ClusEtaPV_ClusADC_all_zoomin_sim.Integral(-1,-1,-1,-1)
0223     print ('Number of clusters in histogram (data) = ', hM_ClusEtaPV_ClusADC_all_zoomin_data.Integral(-1,-1,-1,-1))
0224     print ('Number of clusters in histogram (sim) = ', hM_ClusEtaPV_ClusADC_all_zoomin_sim.Integral(-1,-1,-1,-1))
0225
0226
0227     # Draw_2Dhist_wFunc(hist, IsData, logz, norm1, rmargin, XaxisName, YaxisName, ZaxisName, drawopt, outname)
0228     # Draw_2Dhist_wFunc(hM_ClusEtaPV_ClusADC_all_zoomin_data, True, False, False, 0.08, 'Cluster #eta', 'Cluster ADC', '', 'box', plotdir+'/ClusADCCutOptim/hM_ClusEtaPV_ClusADC_all_zoomin_data_wTF1Cosh')
0229     # Draw_2Dhist_wFunc(hM_ClusEtaPV_ClusADC_all_zoomin_sim, False, False, False, 0.08, 'Cluster #eta', 'Cluster ADC', '', 'box', plotdir+'/ClusADCCutOptim/hM_ClusEtaPV_ClusADC_all_zoomin_sim_wTF1Cosh')
0230
0231     fraclost_constscale_data = []
0232     fracpreserve_constscale_sim = []
0233     mindistto1 = 999.
0234     bestscale = 0
0235     bestfraclos_data = 0
0236     bestfracpreserve_sim = 0
0237
0238     l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut = []
0239     l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut = []
0240     for i in range(smin, smax, sstep):
0241         l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut.append(GetHistogram(datahistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin_constscale{:d}_etascale1p0'.format(i)))
0242         l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut.append(GetHistogram(simhistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin_constscale{:d}_etascale1p0'.format(i)))
0243
0244         Nclus_data = l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut[int((i-smin)/sstep)].Integral(-1,-1,-1,-1)
0245         Nclus_sim = l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut[int((i-smin)/sstep)].Integral(-1,-1,-1,-1)
0246         fl = (NClus_total_data - Nclus_data) / NClus_total_data
0247         fp = Nclus_sim / NClus_total_sim
0248         fraclost_constscale_data.append(fl)
0249         fracpreserve_constscale_sim.append(fp)
0250         # print ('Number of clusters in histogram (data) for ADC cut with const. scale {} = '.format(i), l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut[int((i-smin)/sstep)].Integral(-1,-1,-1,-1))
0251         # print ('Number of clusters in histogram (sim) = ', l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut[int((i-smin)/sstep)].Integral(-1,-1,-1,-1))
0252         if fp > 0.945 and fp < 0.955:
0253             print ('Scale = {}, fraction lost in data = {:.4f}, fraction preserved in sim = {:.4f}'.format(i, fl, fp))
0254
0255         if fp > 0.895 and fp < 0.905:
0256             print ('Scale = {}, fraction lost in data = {:.4f}, fraction preserved in sim = {:.4f}'.format(i, fl, fp))
0257
0258         if fp > 0.845 and fp < 0.855:
0259             print ('Scale = {}, fraction lost in data = {:.4f}, fraction preserved in sim = {:.4f}'.format(i, fl, fp))
0260
0261         if fp > 0.795 and fp < 0.805:
0262             print ('Scale = {}, fraction lost in data = {:.4f}, fraction preserved in sim = {:.4f}'.format(i, fl, fp))
0263
0264         # Calculate the distance to the point (1,1), the closer the better
0265         distto1 = math.sqrt((1-fl)**2 + (1-fp)**2)
0266         if distto1 < mindistto1:
0267             mindistto1 = distto1
0268             bestscale = i
0269             bestfraclos_data = fl
0270             bestfracpreserve_sim = fp
0271
0272
0273         # Draw_1Dhist_datasimcomp(hdata, hsims, gpadmargin, norm, logy, ymaxscale, XaxisName, Ytitle_unit, prelim, simlegtex, evtseltexts, outname):
0274         # plot the cluster eta distribution for each const scale
0275         hM_ClusEtaPV_EtaDepADCCut_constscale_data = GetHistogram(datahistfile, 'hM_ClusEtaPV_EtaDepADCCut_constscale{:d}_etascale1p0'.format(i))
0276         hM_ClusEtaPV_EtaDepADCCut_constscale_sim = GetHistogram(simhistfile, 'hM_ClusEtaPV_EtaDepADCCut_constscale{:d}_etascale1p0'.format(i))
0277         # Draw_1Dhist_datasimcomp(hM_cluseta_zvtxwei_data, l_hM_cluseta_zvtxwei_sim, [0.1,0.08,0.15,0.13], 'data', False, 1.8, 'Cluster #eta', '', False, simlegtext, ['Cluster ADC > 35'], './DataSimComp/{}/Cluster_Eta_zvtxwei'.format(plotdir))
0278         Draw_1Dhist_datasimcomp(hM_ClusEtaPV_EtaDepADCCut_constscale_data, [hM_ClusEtaPV_EtaDepADCCut_constscale_sim], [0.1,0.08,0.15,0.13], 'data', False, 1.8, 'Cluster #eta', '', False, ['Simulation'], ['Cluster ADC cut > {:d}#timescosh(1#times#eta)'.format(i)], plotdir+'/ClusADCCutOptim/Hist1D_ClusEtaPV/hM_ClusEtaPV_EtaDepADCCut_constscale{:d}_etascale1p0'.format(i))
0279         # Draw_2Dhist_wFunc(hist_wcut, hist_wocut, IsData, logz, norm1, rmargin, XaxisName, YaxisName, ZaxisName, drawopt_wcut, drawopt_wocut, outname):
0280         Draw_2Dhist_wFunc(l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut[int((i-smin)/sstep)], hM_ClusEtaPV_ClusADC_all_zoomin_data, True, False, False, 0.08, 'Cluster #eta', 'Cluster ADC', '', 'col', 'box same', plotdir+'/ClusADCCutOptim/Hist2D_ClusEtaPVADC_wStepFunc/hM_ClusEtaPV_ClusADC_all_zoomin_data_wTF1Cosh_constscale{:d}_etascale1p0'.format(i))
0281         Draw_2Dhist_wFunc(l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut[int((i-smin)/sstep)], hM_ClusEtaPV_ClusADC_all_zoomin_sim, False, False, False, 0.08, 'Cluster #eta', 'Cluster ADC', '', 'col', 'box same', plotdir+'/ClusADCCutOptim/Hist2D_ClusEtaPVADC_wStepFunc/hM_ClusEtaPV_ClusADC_all_zoomin_sim_wTF1Cosh_constscale{:d}_etascale1p0'.format(i))
0282
0283
0284     fraclost_etascale_data = []
0285     fracpreserve_etascale_sim = []
0286     l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut_escale = []
0287     l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut_escale = []
0288     for i in range(0, int(float((esmax-esmin)/esstep))):
0289         es = esmin + i*esstep
0290         esstr = '{:.2f}'.format(es).replace('.','p')
0291         print (es, esstr)
0292         l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut_escale.append(GetHistogram(datahistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin_constscale30_etascale{}'.format(esstr)))
0293         l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut_escale.append(GetHistogram(simhistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin_constscale30_etascale{}'.format(esstr)))
0294         hM_ClusEtaPV_EtaDepADCCut_etascale_data = GetHistogram(datahistfile, 'hM_ClusEtaPV_EtaDepADCCut_constscale30_etascale{}'.format(esstr))
0295         hM_ClusEtaPV_EtaDepADCCut_etascale_sim = GetHistogram(simhistfile, 'hM_ClusEtaPV_EtaDepADCCut_constscale30_etascale{}'.format(esstr))
0296
0297         Nclus_data = l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut_escale[i].Integral(-1,-1,-1,-1)
0298         Nclus_sim = l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut_escale[i].Integral(-1,-1,-1,-1)
0299         fl = (NClus_total_data - Nclus_data) / NClus_total_data
0300         fp = Nclus_sim / NClus_total_sim
0301         fraclost_etascale_data.append(fl)
0302         fracpreserve_etascale_sim.append(fp)
0303
0304         Draw_1Dhist_datasimcomp(hM_ClusEtaPV_EtaDepADCCut_etascale_data, [hM_ClusEtaPV_EtaDepADCCut_etascale_sim], [0.1,0.08,0.15,0.13], 'data', False, 1.8, 'Cluster #eta', '', False, ['Simulation'], ['Cluster ADC cut > 30#timescosh({:.2f}#times#eta)'.format(es)], plotdir+'/ClusADCCutOptim/Hist1D_ClusEtaPV/hM_ClusEtaPV_EtaDepADCCut_constscale30_etascale{}'.format(esstr))
0305         Draw_2Dhist_wFunc(l_hM_ClusEtaPV_ClusADC_all_zoomin_data_adccut_escale[i], hM_ClusEtaPV_ClusADC_all_zoomin_data, True, False, False, 0.08, 'Cluster #eta', 'Cluster ADC', '', 'col', 'box same', plotdir+'/ClusADCCutOptim/Hist2D_ClusEtaPVADC_wStepFunc/hM_ClusEtaPV_ClusADC_all_zoomin_data_wTF1Cosh_constscale30_etascale{}'.format(esstr))
0306         Draw_2Dhist_wFunc(l_hM_ClusEtaPV_ClusADC_all_zoomin_sim_adccut_escale[i], hM_ClusEtaPV_ClusADC_all_zoomin_sim, False, False, False, 0.08, 'Cluster #eta', 'Cluster ADC', '', 'col', 'box same', plotdir+'/ClusADCCutOptim/Hist2D_ClusEtaPVADC_wStepFunc/hM_ClusEtaPV_ClusADC_all_zoomin_sim_wTF1Cosh_constscale30_etascale{}'.format(esstr))
0307
0308
0309     fraclost_inttprivate_data = []
0310     fracpreserve_inttprivate_sim = []
0311     hM_ClusEtaPV_EtaDepADCCut_inttprivate_data = GetHistogram(datahistfile, 'hM_ClusEtaPV_EtaDepADCCut_inttprivate')
0312     hM_ClusEtaPV_EtaDepADCCut_inttprivate_sim = GetHistogram(simhistfile, 'hM_ClusEtaPV_EtaDepADCCut_inttprivate')
0313     hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_data = GetHistogram(datahistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate')
0314     hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_sim = GetHistogram(simhistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate')
0315     Nclus_data = hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_data.Integral(-1,-1,-1,-1)
0316     Nclus_sim = hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_sim.Integral(-1,-1,-1,-1)
0317     fl = (NClus_total_data - Nclus_data) / NClus_total_data
0318     fp = Nclus_sim / NClus_total_sim
0319     fraclost_inttprivate_data.append(fl)
0320     fracpreserve_inttprivate_sim.append(fp)
0321     Draw_1Dhist_datasimcomp(hM_ClusEtaPV_EtaDepADCCut_inttprivate_data, [hM_ClusEtaPV_EtaDepADCCut_inttprivate_sim], [0.1,0.08,0.15,0.13], 'data', False, 1.8, 'Cluster #eta', '', False, ['Simulation'], ['Cluster ADC cut (INTT private study)'], plotdir+'/ClusADCCutOptim/Hist1D_ClusEtaPV/hM_ClusEtaPV_EtaDepADCCut_inttprivate')
0322     hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_data = GetHistogram(datahistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate')
0323     hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_sim = GetHistogram(simhistfile, 'hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate')
0324     Draw_2Dhist_wFunc(hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_data, hM_ClusEtaPV_ClusADC_all_zoomin_data, True, False, False, 0.08, 'Cluster #eta', 'Cluster ADC', '', 'col', 'box same', plotdir+'/ClusADCCutOptim/Hist2D_ClusEtaPVADC_wStepFunc/hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_data_wTF1Cosh')
0325     Draw_2Dhist_wFunc(hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_sim, hM_ClusEtaPV_ClusADC_all_zoomin_sim, False, False, False, 0.08, 'Cluster #eta', 'Cluster ADC', '', 'col', 'box same', plotdir+'/ClusADCCutOptim/Hist2D_ClusEtaPVADC_wStepFunc/hM_ClusEtaPV_ClusADC_all_zoomin_inttprivate_sim_wTF1Cosh')
0326
0327
0328
0329
0330     # print('Best scale = {}, with minimum distance to (1,1) = {:.4f}, fraction lost in data = {:.4f}, fraction preserved in sim = {:.4f}'.format(bestscale, mindistto1, bestfraclos_data, bestfracpreserve_sim))
0331
0332     # Draw TGraph for the fraction of lost clusters, X-axis is data and Y-axis is sim
0333     # Get the maximum and minimum from the lists
0334     minfraclost = min(min(fraclost_constscale_data), min(fraclost_etascale_data), min(fraclost_inttprivate_data))
0335     maxfraclost = max(max(fraclost_constscale_data), max(fraclost_etascale_data), max(fraclost_inttprivate_data))
0336     minfracpreserve = min(min(fracpreserve_constscale_sim), min(fracpreserve_etascale_sim), min(fracpreserve_inttprivate_sim))
0337     maxfracpreserve = max(max(fracpreserve_constscale_sim), max(fracpreserve_etascale_sim), max(fracpreserve_inttprivate_sim))
0338
0339     g_fraclost_constscale = TGraph(len(fraclost_constscale_data), array('d', fraclost_constscale_data), array('d', fracpreserve_constscale_sim))
0340     g_fraclost_etascale = TGraph(len(fraclost_etascale_data), array('d', fraclost_etascale_data), array('d', fracpreserve_etascale_sim))
0341     g_fraclost_inttprivate = TGraph(len(fraclost_inttprivate_data), array('d', fraclost_inttprivate_data), array('d', fracpreserve_inttprivate_sim))
0342     c_fraclost = TCanvas('c_fraclost', 'c_fraclost', 800, 700)
0343     c_fraclost.cd()
0344     c_fraclost.SetTopMargin(TopMargin)
0345     c_fraclost.SetLeftMargin(LeftMargin)
0346     c_fraclost.SetBottomMargin(BottomMargin)
0347     c_fraclost.SetRightMargin(RightMargin)
0348     g_fraclost_inttprivate.SetMarkerStyle(20)
0349     g_fraclost_inttprivate.SetMarkerSize(1)
0350     g_fraclost_inttprivate.SetMarkerColor(kBlack)
0351     g_fraclost_inttprivate.SetLineColor(kBlack)
0352     g_fraclost_inttprivate.SetLineWidth(2)
0353     g_fraclost_inttprivate.GetXaxis().SetTitle('Fraction of clusters excluded in data')
0354     g_fraclost_inttprivate.GetYaxis().SetTitle('Fraction of clusters perserved in sim')
0355     g_fraclost_inttprivate.GetXaxis().SetTitleSize(AxisTitleSize)
0356     g_fraclost_inttprivate.GetYaxis().SetTitleSize(AxisTitleSize)
0357     g_fraclost_inttprivate.GetXaxis().SetLabelSize(AxisLabelSize)
0358     g_fraclost_inttprivate.GetYaxis().SetLabelSize(AxisLabelSize)
0359     g_fraclost_inttprivate.GetXaxis().SetTitleOffset(1.1)
0360     g_fraclost_inttprivate.GetYaxis().SetTitleOffset(1.3)
0361     g_fraclost_inttprivate.GetXaxis().SetTickSize(TickSize)
0362     g_fraclost_inttprivate.GetYaxis().SetTickSize(TickSize)
0363     g_fraclost_inttprivate.GetXaxis().SetLimits(0, maxfraclost+0.03)
0364     g_fraclost_inttprivate.GetHistogram().SetMinimum(minfracpreserve-0.03)
0365     g_fraclost_inttprivate.GetHistogram().SetMaximum(1.02)
0366     g_fraclost_inttprivate.Draw('AP')
0367     g_fraclost_constscale.SetMarkerStyle(20)
0368     g_fraclost_constscale.SetMarkerSize(1)
0369     g_fraclost_constscale.SetMarkerColor(TColor.GetColor('#0F4C75'))
0370     g_fraclost_constscale.SetLineColor(TColor.GetColor('#0F4C75'))
0371     g_fraclost_constscale.SetLineWidth(2)
0372     g_fraclost_constscale.Draw('P same')
0373     g_fraclost_etascale.SetMarkerStyle(20)
0374     g_fraclost_etascale.SetMarkerSize(1)
0375     g_fraclost_etascale.SetMarkerColor(TColor.GetColor('#810000'))
0376     g_fraclost_etascale.SetLineColor(TColor.GetColor('#810000'))
0377     g_fraclost_etascale.SetLineWidth(2)
0378     g_fraclost_etascale.Draw('P same')
0379     # legend
0380     leg = TLegend(0.18, 0.17, 0.4, 0.35)
0381     leg.SetHeader('Cluster ADC cut = c#timescosh(b#eta)')
0382     leg.SetTextSize(0.04)
0383     leg.SetFillStyle(0)
0384     leg.AddEntry(g_fraclost_constscale, 'Vary constant scale c (fix #eta scale)', 'p')
0385     leg.AddEntry(g_fraclost_etascale, 'Vary #eta scale b (fix constant scale)', 'p')
0386     leg.AddEntry(g_fraclost_inttprivate, 'INTT private study', 'p')
0387     leg.Draw()
0388     c_fraclost.Draw()
0389     c_fraclost.SaveAs(plotdir+'/ClusADCCutOptim/g_fraclost.pdf')
0390     c_fraclost.SaveAs(plotdir+'/ClusADCCutOptim/g_fraclost.png')
0391
0392
0393