analysis_INTT/plot/closure.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, TFile, TCanvas, TLegend, TColor, gROOT, gSystem, gPad, kGreen, kBlack, kOrange, TGraphAsymmErrors
0008 import numpy
0009 import math
0010 import glob
0011 from plotUtil import GetHistogram, colorset
0012
0013 TickSize = 0.03
0014 AxisTitleSize = 0.05
0015 AxisLabelSize = 0.04
0016 LeftMargin = 0.15
0017 RightMargin = 0.08
0018 TopMargin = 0.08
0019 BottomMargin = 0.13
0020
0021 f_phobosAuAu200GeV = '/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/measurements/PHOBOS-PhysRevC.83.024913/auau_200GeV.root'
0022
0023 gROOT.SetBatch(True)
0024
0025 def Draw_1Dhist_datasimcomp(hdatas, hsims, gpadmargin, logy, ymaxscale, XaxisName, YaxisName, Ytitle_unit, prelim, addstr, datalegtex, simlegtex, outname):
0026     nhists = len(hdatas) + len(hsims)
0027     hsimcolor = colorset(len(hsims))
0028
0029     # get the centrality low and high from addstr: Centrality [0-5]%
0030     centstr = addstr.split(',')[0].split('[')[1].split(']')[0].split('-')
0031     centlow = int(centstr[0])
0032     centhigh = int(centstr[1])
0033     # get phobos measurement: the results are TGraphAsymmErrors
0034     f_phobos = TFile(f_phobosAuAu200GeV, 'READ')
0035     # get the TGraphAsymmErrors, if not found, return None and skip
0036     g_phobos = f_phobos.Get('AuAu_200GeV_Centrality_{}to{}'.format(centlow, centhigh))
0037     if g_phobos:
0038         nhists+=1
0039     else:
0040         print('PHOBOS measurement not found for centrality {}-{}%'.format(centlow, centhigh))
0041
0042     maxbincontent = -1E6
0043     xrange_low = 999
0044     xrange_high = -999
0045
0046     gdatas = []
0047     for hdata in hdatas:
0048         hdata.Sumw2()
0049         # maxbincontent = max(hdata.GetMaximum()+hdata.GetBinError(hdata.GetMaximumBin()), maxbincontent)
0050         tgae = TGraphAsymmErrors(hdata)
0051         # remove the points with zero bin content
0052         for i in range(tgae.GetN()-1, -1, -1):
0053             if tgae.GetY()[i] == 0:
0054                 tgae.RemovePoint(i)
0055         gdatas.append(tgae)
0056
0057     xrange_low = hdatas[0].GetXaxis().GetXmin()
0058     xrange_high = hdatas[0].GetXaxis().GetXmax()
0059
0060     for hsim in hsims:
0061         hsim.Sumw2()
0062         # Get the maximum bin content
0063         # maxbincontent = max(hsim.GetMaximum()+hsim.GetBinError(hsim.GetMaximumBin()), maxbincontent)
0064         # Get the x-axis range
0065         if hsim.GetXaxis().GetXmin() < xrange_low:
0066             xrange_low = hsim.GetXaxis().GetXmin()
0067         if hsim.GetXaxis().GetXmax() > xrange_high:
0068             xrange_high = hsim.GetXaxis().GetXmax()
0069
0070     binwidth = hdatas[0].GetXaxis().GetBinWidth(1)
0071
0072     # get the maximum bin content from HIJING truth, the first histogram in hsims
0073     maxbincontent = hsims[0].GetMaximum()
0074     print ("maxbincontent: ", maxbincontent)
0075     c = TCanvas('c', 'c', 800, 700)
0076     if logy:
0077         c.SetLogy()
0078     c.cd()
0079     gPad.SetRightMargin(gpadmargin[0])
0080     gPad.SetTopMargin(gpadmargin[1])
0081     gPad.SetLeftMargin(gpadmargin[2])
0082     gPad.SetBottomMargin(gpadmargin[3])
0083
0084     for i, hsim in enumerate(hsims):
0085         hsim.SetLineColor(TColor.GetColor(hsimcolor[i]))
0086         hsim.SetLineWidth(2)
0087         if 'generator' in simlegtex[i]: # HIING truth always the first
0088             hsim.GetYaxis().SetTitle(YaxisName)
0089             hsim.GetYaxis().SetRangeUser(0, maxbincontent * ymaxscale)
0090             # if logy:
0091             #     hsim.GetYaxis().SetRangeUser(hdata.GetMinimum(0)*0.5, maxbincontent * ymaxscale)
0092             # else:
0093             #     hsim.GetYaxis().SetRangeUser(0., maxbincontent * ymaxscale)
0094
0095             hsim.GetXaxis().SetTitle(XaxisName)
0096             # hsim.GetXaxis().SetRangeUser(xrange_low, xrange_high)
0097             hsim.GetXaxis().SetRange(hsim.FindFirstBinAbove(0),hsim.FindLastBinAbove(0))
0098             hsim.GetXaxis().SetTitleOffset(1.1)
0099             hsim.GetYaxis().SetTitleOffset(1.5)
0100             hsim.SetMarkerSize(0)
0101             hsim.Draw('histe')
0102         else:
0103             hsim.SetFillColorAlpha(TColor.GetColor(hsimcolor[i]), 0.3)
0104             # hsim.SetMarkerSize(0.6)
0105             hsim.SetMarkerColor(TColor.GetColor(hsimcolor[i]))
0106             hsim.Draw('E5 same')
0107
0108     for i, gdata in enumerate(gdatas): # the first is the CMS approach, the second is the PHOBOS approach
0109         # if hdata.Integral(-1,-1) == 0:
0110         #     continiue
0111
0112         gdata.SetMarkerStyle(34 if i==0 else 47)
0113         gdata.SetMarkerSize(1)
0114         gdata.SetMarkerColor(TColor.GetColor('#FF6700') if i == 0 else TColor.GetColor('#228B22'))
0115         gdata.SetLineColor(TColor.GetColor('#FF6700') if i == 0 else TColor.GetColor('#228B22'))
0116         gdata.SetLineWidth(1)
0117         gdata.SetFillColorAlpha(TColor.GetColor('#FF6700') if i == 0 else TColor.GetColor('#228B22'), 0.3)
0118         gdata.GetXaxis().SetRange(hdata.FindFirstBinAbove(0),hdata.FindLastBinAbove(0))
0119         gdata.Draw('5p same')
0120
0121     # Draw the PHOBOS measurement
0122     if g_phobos:
0123         g_phobos.SetMarkerStyle(25)
0124         g_phobos.SetMarkerSize(0.8)
0125         g_phobos.SetMarkerColor(TColor.GetColor('#035397'))
0126         g_phobos.SetLineColor(TColor.GetColor('#035397'))
0127         g_phobos.SetLineWidth(0)
0128         g_phobos.SetFillColorAlpha(TColor.GetColor('#035397'), 0.3)
0129         g_phobos.Draw('p same')
0130         g_phobos.Draw('3 same')
0131
0132
0133     shift = 0.45 if prelim else 0.75
0134     legyspace = 0.08
0135     if not g_phobos:
0136         legyspace = 0.085
0137     leg = TLegend((1-RightMargin)-shift, BottomMargin + 0.03,
0138                   (1-RightMargin)-0.1, BottomMargin + 0.03 + legyspace * nhists)
0139     leg.SetTextSize(0.035)
0140     leg.SetFillStyle(0)
0141     prelimtext = 'Preliminary' if prelim else 'Internal'
0142     leg.AddEntry('', '#it{#bf{sPHENIX}} '+prelimtext, '')
0143     leg.AddEntry('', 'Au+Au #sqrt{s_{NN}}=200 GeV', '')
0144     leg.AddEntry('', addstr, '')
0145     for i, lt in enumerate(datalegtex):
0146         leg.AddEntry(gdatas[i], lt, "pf");
0147     for i, lt in enumerate(simlegtex):
0148         if 'generator' in lt:
0149             leg.AddEntry(hsims[i], lt, "le");
0150         else:
0151             leg.AddEntry(hsims[i], lt, "pl");
0152     if g_phobos:
0153         leg.AddEntry(g_phobos, 'PHOBOS [Phys. Rev. C 83, 024913]', 'pf')
0154     leg.Draw()
0155     c.RedrawAxis()
0156     c.Draw()
0157     c.SaveAs(outname+'.pdf')
0158     c.SaveAs(outname+'.png')
0159     if(c):
0160         c.Close()
0161         gSystem.ProcessEvents()
0162         del c
0163         c = 0
0164
0165
0166 def GetMbinNum(fstr):
0167     mbinnum = -1
0168     if 'Centrality0to3' in fstr:
0169         mbinnum = 0
0170     elif 'Centrality3to6' in fstr:
0171         mbinnum = 1
0172     elif 'Centrality6to10' in fstr:
0173         mbinnum = 2
0174     elif 'Centrality10to15' in fstr:
0175         mbinnum = 3
0176     elif 'Centrality15to20' in fstr:
0177         mbinnum = 4
0178     elif 'Centrality20to25' in fstr:
0179         mbinnum = 5
0180     elif 'Centrality25to30' in fstr:
0181         mbinnum = 6
0182     elif 'Centrality30to35' in fstr:
0183         mbinnum = 7
0184     elif 'Centrality35to40' in fstr:
0185         mbinnum = 8
0186     elif 'Centrality40to45' in fstr:
0187         mbinnum = 9
0188     elif 'Centrality45to50' in fstr:
0189         mbinnum = 10
0190     elif 'Centrality50to55' in fstr:
0191         mbinnum = 11
0192     elif 'Centrality55to60' in fstr:
0193         mbinnum = 12
0194     elif 'Centrality60to65' in fstr:
0195         mbinnum = 13
0196     elif 'Centrality65to70' in fstr:
0197         mbinnum = 14
0198     elif 'Centrality0to70' in fstr:
0199         mbinnum = 70
0200     return mbinnum
0201
0202
0203 if __name__ == '__main__':
0204     parser = OptionParser(usage='usage: %prog ver [options -h]')
0205     parser.add_option('--datahistdir', dest='datahistdir', type='string', default='/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/systematics/', help='Histogram file name (data)')
0206     parser.add_option('--simhistdir', action='append', dest='simhistdir', type='string', help='Histogram file name (simulation). Example: /sphenix/user/hjheng/TrackletAna/analysis_INTT/plot/corrections/closure_HIJING_set2')
0207     parser.add_option('--filedesc', dest='filedesc', type='string', default='Centrality0to5_Zvtxm30p0tom10p0_noasel', help='String for input file and output plot description (cuts, binnings, etc...)')
0208     # parser.add_option('--simlegtext', action='append', dest='simlegtext', type='string', help='Legend text for simulation. Example: HIJING/EPOS/AMPT)')
0209     parser.add_option('--docompare', action='store_true', dest='docompare', default=False, help='Compare two analyses')
0210     parser.add_option('--plotdir', dest='plotdir', type='string', default='cross-checks', help='Plot directory')
0211
0212     (opt, args) = parser.parse_args()
0213     print('opt: {}'.format(opt))
0214
0215     datahistdir = opt.datahistdir
0216     simhistdir = opt.simhistdir
0217     # simlegtext = opt.simlegtext
0218     plotdir = opt.plotdir
0219     filedesc = opt.filedesc
0220     docompare = opt.docompare
0221
0222     os.makedirs('./corrections/{}'.format(plotdir), exist_ok=True)
0223
0224     MbinNum = GetMbinNum(filedesc)
0225     simlegtext = []
0226
0227     h1WEfinal_data = GetHistogram("{}/{}/finalhists_systematics_{}.root".format(datahistdir,filedesc,filedesc), 'hM_final')
0228     # loop over the bins, for |eta| > 1.1, set the bin content to 0
0229     for i in range(1, h1WEfinal_data.GetNbinsX()+1):
0230         if abs(h1WEfinal_data.GetBinCenter(i)) > 1.1:
0231             h1WEfinal_data.SetBinContent(i, 0)
0232             h1WEfinal_data.SetBinError(i, 0)
0233
0234     l_h1WEfinal_sim = []
0235     # l_h1WGhadron_sim = []
0236     for i, simhistd in enumerate(simhistdir):
0237         # tmplegtxt = simlegtext[i]
0238
0239         l_h1WEfinal_sim.append(GetHistogram("{}/{}/correction_hists.root".format(simhistd,filedesc), 'h1WGhadron'))
0240         simlegtext.append('HIJING (generator)')
0241         # From the CMS-style analysis
0242         # l_h1WEfinal_sim.append(GetHistogram("{}/{}/correction_hists.root".format(simhistd,filedesc), 'h1WEfinal'))
0243         # simlegtext.append('Simulation closure [The closest-match method]')
0244         # From the PHOBOS-style analysis
0245         # if docompare:
0246         #     l_h1WEfinal_sim.append(GetHistogram("/sphenix/user/ChengWei/INTT/INTT_commissioning/ZeroField/F4A_20869/2024_05_07/folder_Data_CombinedNtuple_Run20869_HotDead_BCO_ADC_Survey/evt_tracklet/complete_file/merged_file_folder/eta_closer_out_test_MultiZBin/hist_for_AN.root", 'dNdeta_1D_MCreco_loose_PostCorrection_corr_Mbin{}'.format(MbinNum)))
0247         #     simlegtext.append('Simulation closure [The combinatoric method]')
0248         # Raw generated hadron distribution
0249
0250     if docompare:
0251         hM_dNdeta_1D_Datareco_tight_PostCorrection = GetHistogram("/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/Run6_EvtZFitWidthChange/EvtVtxZ/FinalResult_10cm_Pol2BkgFit/completed/vtxZ_-10_10cm_MBin{}/Final_Mbin{}_00054280/Final_Mbin{}_00054280.root".format(MbinNum, MbinNum, MbinNum), 'h1D_dNdEta_reco')
0252         for i in range(1, hM_dNdeta_1D_Datareco_tight_PostCorrection.GetNbinsX()+1):
0253             if abs(hM_dNdeta_1D_Datareco_tight_PostCorrection.GetBinCenter(i)) > 1.1:
0254                 hM_dNdeta_1D_Datareco_tight_PostCorrection.SetBinContent(i, 0)
0255                 hM_dNdeta_1D_Datareco_tight_PostCorrection.SetBinError(i, 0)
0256
0257     # split the string by '_' and replace some characters
0258     descstr = filedesc.split('_')
0259     centstr = descstr[0]
0260     zvtxstr = descstr[1]
0261     aselstr = descstr[2]
0262
0263     centstr = centstr.replace('Centrality', 'Centrality [')
0264     centstr = centstr.replace('to', '-')
0265     centstr = centstr + ']%'
0266
0267     zvtxstr = zvtxstr.replace('m', '-')
0268     zvtxstr = zvtxstr.replace('p', '.')
0269     zvtxstr = zvtxstr.replace('to', ', ')
0270     zvtxstr = zvtxstr.replace('Zvtx', 'vtx_{Z} [')
0271     zvtxstr = zvtxstr + '] cm'
0272
0273     aselstr = '' if aselstr == 'noasel' else ', (with additional selection)'
0274
0275     # legstr = centstr + ', ' + zvtxstr + aselstr
0276     legstr = centstr
0277
0278     # Draw_1Dhist_datasimcomp(hdatas, hsims, gpadmargin, logy, ymaxscale, XaxisName, YaxisName, Ytitle_unit, prelim, addstr, datalegtex, simlegtex, outname):
0279     if docompare:
0280         Draw_1Dhist_datasimcomp([h1WEfinal_data, hM_dNdeta_1D_Datareco_tight_PostCorrection], l_h1WEfinal_sim, [0.06,0.06,0.15,0.13], False, 1.5, 'Pseudorapidity #eta', 'dN_{ch}/d#eta', '', False, legstr, ['Data - The closest-match method', 'Data - The combinatoric method'], simlegtext, './corrections/{}/dNdeta_{}_crosscheck'.format(plotdir,filedesc))
0281     else:
0282         Draw_1Dhist_datasimcomp([h1WEfinal_data], l_h1WEfinal_sim, [0.06,0.06,0.15,0.13], False, 1.5, 'Pseudorapidity #eta', 'dN_{ch}/d#eta', '', False, legstr, ['Data - The closest-match method'], simlegtext, './corrections/{}/dNdeta_{}_crosscheck'.format(plotdir,filedesc))