Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:17 

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 def Draw_1Dhist_datasimcomp_evtsel(hdatas, hsims, gpadmargin, norm, logy, ymaxscale, XaxisName, Ytitle_unit, prelim, evtseltexts, outname):
0016     # make sure the number of histograms to compare are the same
0017     if len(hdatas) != len(hsims):
0018         print('Number of histograms to compare are not the same')
0019         sys.exit(1)
0020         
0021     nhists = len(hdatas)
0022     hcolor = ['#4a536b','#608BC1','#B03052']
0023 
0024     for i in range(nhists):
0025         hdatas[i].Sumw2()
0026         hsims[i].Sumw2()
0027 
0028     binwidth = hdatas[0].GetXaxis().GetBinWidth(1)
0029 
0030     # Get the maximum bin content 
0031     maxbincontent = -1
0032     minbincontent = 1E9
0033     if norm == 'unity':
0034         for i in range(nhists):
0035             hdatas[i].Scale(1.0 / hdatas[i].Integral(-1,-1))
0036             hsims[i].Scale(1.0 / hsims[i].Integral(-1,-1))
0037             maxbinc = max(hdatas[i].GetMaximum(), hsims[i].GetMaximum())
0038             minbinc = min(hdatas[i].GetMinimum(0), hsims[i].GetMinimum(0))
0039             if maxbinc > maxbincontent:
0040                 maxbincontent = maxbinc
0041             if minbinc < minbincontent:
0042                 minbincontent = minbinc
0043     elif norm == 'data': # normalized to sim to data 
0044         for i in range(nhists):
0045             hsims[i].Scale(hdatas[i].Integral(-1,-1) / hsims[i].Integral(-1,-1))
0046             maxbinc = max(hdatas[i].GetMaximum(), hsims[i].GetMaximum())
0047             minbinc = min(hdatas[i].GetMinimum(0), hsims[i].GetMinimum(0))
0048             if maxbinc > maxbincontent:
0049                 maxbincontent = maxbinc
0050             if minbinc < minbincontent:
0051                 minbincontent = minbinc
0052     else:
0053         if norm != 'none':
0054             print('Invalid normalization option: {}'.format(norm))
0055             sys.exit(1)
0056 
0057     c = TCanvas('c', 'c', 800, 700)
0058     pad1 = TPad( 'pad1', ' ', 0, 0.3, 1, 1)
0059     pad2 = TPad( 'pad2', ' ', 0, 0, 1, 0.3)
0060     pad1.SetRightMargin(gpadmargin[0])
0061     pad1.SetTopMargin(gpadmargin[1])
0062     pad1.SetLeftMargin(gpadmargin[2])
0063     pad1.SetBottomMargin(0.0)
0064     pad1.Draw()
0065     pad2.SetRightMargin(gpadmargin[0])
0066     pad2.SetTopMargin(0.0)
0067     pad2.SetLeftMargin(gpadmargin[2])
0068     pad2.SetBottomMargin(gpadmargin[3])
0069     pad2.Draw() # Draw the TPad on the TCanvas before plotting something on TPad
0070     # cd to the pad1
0071     pad1.cd()
0072     if logy:
0073         pad1.SetLogy()
0074     
0075     for i, hsim in enumerate(hsims):
0076         if i == 0:
0077             if norm == 'unity':
0078                 if Ytitle_unit == '':
0079                     hsim.GetYaxis().SetTitle(
0080                         'Normalized entries / ({:g})'.format(binwidth))
0081                 else:
0082                     hsim.GetYaxis().SetTitle(
0083                         'Normalized entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0084             else:
0085                 if Ytitle_unit == '':
0086                     hsim.GetYaxis().SetTitle('Entries / ({:g})'.format(binwidth))
0087                 else:
0088                     hsim.GetYaxis().SetTitle(
0089                         'Entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0090 
0091             if logy:
0092                 hsim.GetYaxis().SetRangeUser(minbincontent*0.5, maxbincontent * ymaxscale)
0093             else:
0094                 hsim.GetYaxis().SetRangeUser(1E-3, (hsim.GetMaximum()) * ymaxscale)
0095 
0096             # hsim.GetXaxis().SetTitle(XaxisName)
0097             # hsim.GetXaxis().SetTitleOffset(1.1)
0098             hsim.GetXaxis().SetLabelOffset(999)
0099             hsim.GetXaxis().SetLabelSize(0)
0100             hsim.GetYaxis().SetTitleOffset(1.35)
0101             hsim.SetLineColor(TColor.GetColor(hcolor[i]))
0102             hsim.SetLineWidth(2)
0103             hsim.SetMarkerSize(0)
0104             hsim.Draw('histe')
0105         else:
0106             hsim.SetLineColor(TColor.GetColor(hcolor[i]))
0107             hsim.SetLineWidth(2)
0108             hsim.SetMarkerSize(0)
0109             hsim.Draw('histe same')
0110 
0111     for i, hdata in enumerate(hdatas):
0112         hdata.SetMarkerStyle(20)
0113         hdata.SetMarkerSize(0.6)
0114         hdata.SetMarkerColor(TColor.GetColor(hcolor[i]))
0115         hdata.SetLineColor(TColor.GetColor(hcolor[i]))
0116         hdata.SetLineWidth(2)
0117         hdata.Draw('same PE1')
0118     
0119     # dummy hist for legend
0120     hdummy_data = hdatas[0].Clone('hdummy_data')
0121     hdummy_data.SetMarkerStyle(20)
0122     hdummy_data.SetMarkerSize(0.6)
0123     hdummy_data.SetLineWidth(2)
0124     hdummy_data.SetMarkerColor(kBlack)
0125     hdummy_data.SetLineColor(kBlack)
0126     hdummy_sim = hsims[0].Clone('hdummy_sim')
0127     hdummy_sim.SetLineWidth(2)
0128     hdummy_sim.SetLineColor(kBlack)
0129     # Legend
0130     shift = 0.43 if prelim else 0.73
0131     legylow = 0.15 + 0.0 * (len(hsims))
0132     leg = TLegend((1-RightMargin)-shift, (1-TopMargin)-legylow,
0133                   (1-RightMargin)-0.1, (1-TopMargin)-0.02)
0134     leg.SetTextSize(0.05)
0135     leg.SetFillStyle(0)
0136     # prelimtext = 'Preliminary' if prelim else 'Work-in-progress'
0137     prelimtext = 'Preliminary' if prelim else 'Internal'
0138     leg.AddEntry('', '#it{#bf{sPHENIX}} '+prelimtext, '')
0139     leg.AddEntry('', 'Au+Au #sqrt{s_{NN}}=200 GeV', '')
0140     leg.SetNColumns(2)
0141     leg.AddEntry(hdummy_data, 'Data', "pe")
0142     leg.AddEntry(hdummy_sim, 'HIJING', "le")
0143     leg.Draw('same')
0144     
0145     # event selection text
0146     legylow_evtselshift = 0.05 * len(evtseltexts)
0147     leg2 = TLegend((1-RightMargin)-shift, (1-TopMargin)-legylow-legylow_evtselshift,
0148                     (1-RightMargin)-0.3, (1-TopMargin)-legylow)
0149     leg2.SetTextSize(0.04)
0150     leg2.SetFillStyle(0)
0151     for i, evtseltext in enumerate(evtseltexts):
0152         leg2.AddEntry(hsims[i], evtseltext, 'l')
0153     leg2.Draw('same')
0154     c.RedrawAxis()
0155     
0156     c.Update()
0157     # cd to the pad2
0158     textscale = 2.7
0159     pad2.cd()
0160     pad2.SetLogy(0)
0161     pad2.SetGridy(0)
0162     # take the ratio = data/simulation
0163     l_hM_ratio = []
0164     for i, hsim in enumerate(hsims):
0165         hM_ratio = hdatas[i].Clone('hM_ratio')
0166         hM_ratio.Divide(hsim)
0167         # hM_ratio.SetMarkerStyle(20)
0168         # hM_ratio.SetMarkerSize(0.7)
0169         # hM_ratio.SetMarkerColor(TColor.GetColor(hcolor[i]))
0170         hM_ratio.SetLineColor(TColor.GetColor(hcolor[i]))
0171         hM_ratio.SetLineWidth(2)
0172         if i == 0:
0173             hM_ratio.GetYaxis().SetRangeUser(0.001, 1.999)
0174             # hM_ratio.GetYaxis().SetRangeUser(0.501, 1.499)
0175             hM_ratio.GetYaxis().SetNdivisions(505)
0176             hM_ratio.GetYaxis().SetTitle('Data/Sim.')
0177             hM_ratio.GetYaxis().SetTitleOffset(0.5)
0178             hM_ratio.GetYaxis().SetTickSize(TickSize)
0179             hM_ratio.GetYaxis().SetTitleSize(AxisTitleSize*textscale)
0180             hM_ratio.GetYaxis().SetLabelSize(AxisLabelSize*textscale)
0181             hM_ratio.GetXaxis().SetTickSize(TickSize*textscale)
0182             hM_ratio.GetXaxis().SetTitleSize(AxisTitleSize*textscale)
0183             hM_ratio.GetXaxis().SetLabelSize(AxisLabelSize*textscale)
0184             hM_ratio.GetXaxis().SetTitleOffset(1.1)
0185             hM_ratio.GetXaxis().SetTitle(XaxisName)
0186             hM_ratio.Draw('l hist')
0187         else:
0188             hM_ratio.Draw('l hist same')
0189             
0190         l_hM_ratio.append(hM_ratio)
0191     
0192     # Draw a line at 1
0193     line = TLine(hM_ratio.GetXaxis().GetXmin(), 1, hM_ratio.GetXaxis().GetXmax(), 1)
0194     line.SetLineColor(kBlack)
0195     line.SetLineStyle(2)
0196     line.Draw()
0197     c.Draw()
0198     c.SaveAs(outname+'.pdf')
0199     c.SaveAs(outname+'.png')
0200     if(c):
0201         c.Close()
0202         gSystem.ProcessEvents()
0203         del c
0204         c = 0
0205 
0206 if __name__ == '__main__':
0207     plotdir = './DataSimComp/EvtSelComp'
0208     if not os.path.exists(plotdir):
0209         os.makedirs(plotdir)
0210     
0211     datafilestocomp = [
0212         '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_CombinedNtuple_Run54280_HotChannel_BCO/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_TrigBit13InttZvtx30cm/hists_merged.root',
0213         '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_CombinedNtuple_Run54280_HotChannel_BCO/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_TrigBit12InttZvtx10cm/hists_merged.root',
0214         '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_CombinedNtuple_Run54280_HotChannel_BCO/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_TrigBit12InttZvtx10cm_clusPhisizele6/hists_merged.root'
0215     ]
0216     
0217     simfilestocomp = [
0218         '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_Ntuple_HIJING_ana443_20241102/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_TrigBit13InttZvtx30cm/hists_merged.root',
0219         '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_Ntuple_HIJING_ana443_20241102/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_TrigBit12InttZvtx10cm/hists_merged.root',
0220         '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Sim_Ntuple_HIJING_ana443_20241102/dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_TrigBit12InttZvtx10cm_clusPhisizele6/hists_merged.root'
0221     ]
0222     
0223     hdatas_NClusLayer1_zvtxwei = []
0224     hsims_NClusLayer1_zvtxwei = []
0225     hdatas_clusphi_zvtxwei = []
0226     hsims_clusphi_zvtxwei = []
0227     hdatas_cluseta_zvtxwei = []
0228     hsims_cluseta_zvtxwei = []
0229     hdatas_dEta_reco_altrange = []
0230     hsims_dEta_reco_altrange = []
0231     hdatas_dR_reco = []
0232     hsims_dR_reco = []
0233     hdatas_dPhi_reco_altrange = []
0234     hsims_dPhi_reco_altrange = []
0235     hdatas_Eta_reco = []
0236     hsims_Eta_reco = []
0237     hdatas_Phi_reco = []
0238     hsims_Phi_reco = []
0239     
0240     for datafile in datafilestocomp:
0241         hdatas_NClusLayer1_zvtxwei.append(GetHistogram(datafile, 'hM_NClusLayer1_zvtxwei'))
0242         hdatas_clusphi_zvtxwei.append(GetHistogram(datafile, 'hM_clusphi_zvtxwei'))
0243         hdatas_cluseta_zvtxwei.append(GetHistogram(datafile, 'hM_cluseta_zvtxwei'))
0244         hdatas_dEta_reco_altrange.append(GetHistogram(datafile, 'hM_dEta_reco_altrange'))
0245         hdatas_dR_reco.append(GetHistogram(datafile, 'hM_dR_reco'))
0246         hdatas_dPhi_reco_altrange.append(GetHistogram(datafile, 'hM_dPhi_reco_altrange'))
0247         hdatas_Eta_reco.append(GetHistogram(datafile, 'hM_Eta_reco'))
0248         hdatas_Phi_reco.append(GetHistogram(datafile, 'hM_Phi_reco'))
0249     for simfile in simfilestocomp:
0250         hsims_NClusLayer1_zvtxwei.append(GetHistogram(simfile, 'hM_NClusLayer1_zvtxwei'))
0251         hsims_clusphi_zvtxwei.append(GetHistogram(simfile, 'hM_clusphi_zvtxwei'))
0252         hsims_cluseta_zvtxwei.append(GetHistogram(simfile, 'hM_cluseta_zvtxwei'))
0253         hsims_dEta_reco_altrange.append(GetHistogram(simfile, 'hM_dEta_reco_altrange'))
0254         hsims_dR_reco.append(GetHistogram(simfile, 'hM_dR_reco'))
0255         hsims_dPhi_reco_altrange.append(GetHistogram(simfile, 'hM_dPhi_reco_altrange'))
0256         hsims_Eta_reco.append(GetHistogram(simfile, 'hM_Eta_reco'))
0257         hsims_Phi_reco.append(GetHistogram(simfile, 'hM_Phi_reco'))
0258             
0259     print (len(hdatas_NClusLayer1_zvtxwei), len(hsims_NClusLayer1_zvtxwei))
0260     padmargin = [0.08,0.06,0.15,0.32]
0261     # Draw_1Dhist_datasimcomp_evtsel(hdatas, hsims, gpadmargin, norm, logy, ymaxscale, XaxisName, Ytitle_unit, prelim, evtseltexts, outname)
0262     Draw_1Dhist_datasimcomp_evtsel(hdatas_NClusLayer1_zvtxwei, hsims_NClusLayer1_zvtxwei, padmargin, 'data', True, 100, 'Number of clusters (inner)', '', False, ['Trigger bit 13 (MBD S&N#geq2, |Z_{vtx}|<30cm)', 'Trigger bit 12 (MBD S&N#geq2, |Z_{vtx}|<10cm)', 'Trigger bit 12, cluster#phi size<6, ADC>35'], plotdir+'/EvtSelComp_NClusLayer1_zvtxwei')
0263     Draw_1Dhist_datasimcomp_evtsel(hdatas_clusphi_zvtxwei, hsims_clusphi_zvtxwei, padmargin, 'data', False, 1.6, 'Cluster #phi', '', False, ['Trigger bit 13 (MBD S&N#geq2, |Z_{vtx}|<30cm)', 'Trigger bit 12 (MBD S&N#geq2, |Z_{vtx}|<10cm)', 'Trigger bit 12, cluster#phi size<6, ADC>35'], plotdir+'/EvtSelComp_clusphi_zvtxwei')
0264     Draw_1Dhist_datasimcomp_evtsel(hdatas_cluseta_zvtxwei, hsims_cluseta_zvtxwei, padmargin, 'data', False, 1.7, 'Cluster #eta', '', False, ['Trigger bit 13 (MBD S&N#geq2, |Z_{vtx}|<30cm)', 'Trigger bit 12 (MBD S&N#geq2, |Z_{vtx}|<10cm)', 'Trigger bit 12, cluster#phi size<6, ADC>35'], plotdir+'/EvtSelComp_cluseta_zvtxwei')
0265     Draw_1Dhist_datasimcomp_evtsel(hdatas_dEta_reco_altrange, hsims_dEta_reco_altrange, padmargin, 'data', True, 100, '#Delta#eta', '', False, ['Trigger bit 13 (MBD S&N#geq2, |Z_{vtx}|<30cm)', 'Trigger bit 12 (MBD S&N#geq2, |Z_{vtx}|<10cm)', 'Trigger bit 12, cluster#phi size<6, ADC>35'], plotdir+'/EvtSelComp_dEta_reco_altrange')
0266     Draw_1Dhist_datasimcomp_evtsel(hdatas_dR_reco, hsims_dR_reco, padmargin, 'data', True, 100, '#DeltaR', '', False, ['Trigger bit 13 (MBD S&N#geq2, |Z_{vtx}|<30cm)', 'Trigger bit 12 (MBD S&N#geq2, |Z_{vtx}|<10cm)', 'Trigger bit 12, cluster#phi size<6, ADC>35'], plotdir+'/EvtSelComp_dR_reco')
0267     Draw_1Dhist_datasimcomp_evtsel(hdatas_dPhi_reco_altrange, hsims_dPhi_reco_altrange, padmargin, 'data', True, 100, '#Delta#phi', '', False, ['Trigger bit 13 (MBD S&N#geq2, |Z_{vtx}|<30cm)', 'Trigger bit 12 (MBD S&N#geq2, |Z_{vtx}|<10cm)', 'Trigger bit 12, cluster#phi size<6, ADC>35'], plotdir+'/EvtSelComp_dPhi_reco_altrange')
0268     Draw_1Dhist_datasimcomp_evtsel(hdatas_Eta_reco, hsims_Eta_reco, padmargin, 'data', False, 1.7, 'Reco-tracklet #eta', '', False, ['Trigger bit 13 (MBD S&N#geq2, |Z_{vtx}|<30cm)', 'Trigger bit 12 (MBD S&N#geq2, |Z_{vtx}|<10cm)', 'Trigger bit 12, cluster#phi size<6, ADC>35'], plotdir+'/EvtSelComp_Eta_reco')
0269     Draw_1Dhist_datasimcomp_evtsel(hdatas_Phi_reco, hsims_Phi_reco, padmargin, 'data', False, 1.6, 'Reco-tracklet #phi', '', False, ['Trigger bit 13 (MBD S&N#geq2, |Z_{vtx}|<30cm)', 'Trigger bit 12 (MBD S&N#geq2, |Z_{vtx}|<10cm)', 'Trigger bit 12, cluster#phi size<6, ADC>35'], plotdir+'/EvtSelComp_Phi_reco')
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!