Back to home page

sPhenix code displayed by LXR
 
 

    

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

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, TTree, TChain, TLegend, TLatex, TF1, gROOT, gPad, gSystem, kHAlignRight, kVAlignBottom
0008 import numpy as np
0009 import math
0010 import glob
0011 from plotUtil import *
0012 
0013 # gROOT.LoadMacro('./sPHENIXStyle/sPhenixStyle.C')
0014 # gROOT.ProcessLine('SetsPhenixStyle()')
0015 gROOT.SetBatch(True)
0016 
0017 def drawcomparison_data(hists, normhist, logy, histcolor, legentry, legposition, xtitle, ytitle, yscale, outname):
0018     
0019     hmax = -999
0020     for i, h in enumerate(hists):
0021         h.Scale(normhist.Integral(-1,-1)/h.Integral(-1,-1))
0022         if h.GetMaximum() > hmax:
0023             hmax = h.GetMaximum()
0024     
0025     c = TCanvas('c', 'c', 800, 700)
0026     c.cd()
0027     if logy: 
0028         c.SetLogy()
0029     # c.SetLeftMargin(LeftMargin)
0030     # c.SetRightMargin(RightMargin)
0031     c.SetTopMargin(0.18)
0032     # c.SetBottomMargin(BottomMargin)
0033     
0034     # leg = TLegend(legposition[0], legposition[1], legposition[2], legposition[3])
0035     leg = TLegend(gPad.GetLeftMargin(), 1-gPad.GetTopMargin()+0.01, 1-gPad.GetRightMargin(), 0.98)
0036     leg.SetNColumns(2)
0037     leg.SetBorderSize(0)
0038     leg.SetFillStyle(0)
0039     leg.SetTextSize(0.03)
0040     
0041     for i, h in enumerate(hists):
0042         h.SetLineColor(TColor.GetColor(histcolor[i]))
0043         h.SetMarkerColor(TColor.GetColor(histcolor[i]))
0044         h.SetMarkerStyle(20)
0045         h.SetMarkerSize(0.6)
0046         h.SetLineWidth(2)
0047         
0048         if i == 0:
0049             h.GetXaxis().SetTitle(xtitle)
0050             h.GetYaxis().SetTitle(ytitle)
0051             h.GetXaxis().SetTitleOffset(1.2)
0052             h.GetYaxis().SetTitleOffset(1.6)
0053             if logy:
0054                 h.GetYaxis().SetRangeUser(0.5, hmax * yscale)
0055             else:
0056                 h.GetYaxis().SetRangeUser(0., hmax * yscale)
0057                 
0058             h.Draw('pe1')
0059         else:
0060             h.Draw('same pe1')
0061             
0062         leg.AddEntry(h, legentry[i], 'pe1')
0063             
0064         
0065     leg.Draw()
0066     legsphnx = TLegend((1-gPad.GetRightMargin())-0.43, (1-gPad.GetTopMargin())-0.15, (1-gPad.GetRightMargin())-0.11, (1-gPad.GetTopMargin())-0.05)
0067     legsphnx.SetTextSize(0.04)
0068     legsphnx.SetFillStyle(0)
0069     legsphnx.AddEntry("", "#it{#bf{sPHENIX}} Internal", "")
0070     legsphnx.AddEntry("", "Au+Au #sqrt{s_{NN}}=200 GeV", "")
0071     legsphnx.Draw()
0072     c.SaveAs('{}.pdf'.format(outname))
0073     c.SaveAs('{}.png'.format(outname))
0074     if(c):
0075         c.Close()
0076         gSystem.ProcessEvents()
0077         del c
0078         c = 0
0079     
0080     
0081 def Draw_1Dhist_fitGaussian(hist, norm1, logy, ymaxscale, XaxisName, Ytitle_unit, outname):
0082     hist.Sumw2()
0083     binwidth = hist.GetXaxis().GetBinWidth(1)
0084     c = TCanvas('c', 'c', 800, 700)
0085     if norm1:
0086         hist.Scale(1. / hist.Integral(-1, -1))
0087     if logy:
0088         c.SetLogy()
0089     c.cd()
0090     gPad.SetRightMargin(RightMargin)
0091     gPad.SetTopMargin(TopMargin)
0092     gPad.SetLeftMargin(LeftMargin)
0093     gPad.SetBottomMargin(BottomMargin)
0094     if norm1:
0095         if Ytitle_unit == '':
0096             hist.GetYaxis().SetTitle(
0097                 'Normalized entries / ({:g})'.format(binwidth))
0098         else:
0099             hist.GetYaxis().SetTitle(
0100                 'Normalized entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0101     else:
0102         if Ytitle_unit == '':
0103             hist.GetYaxis().SetTitle('Entries / ({:g})'.format(binwidth))
0104         else:
0105             hist.GetYaxis().SetTitle(
0106                 'Entries / ({:g} {unit})'.format(binwidth, unit=Ytitle_unit))
0107 
0108     # hist.GetXaxis().SetRangeUser(hist.GetBinLowEdge(1)-binwidth, hist.GetBinLowEdge(hist.GetNbinsX())+2*binwidth)
0109     if logy:
0110         hist.GetYaxis().SetRangeUser(hist.GetMinimum(0)*0.5, (hist.GetMaximum()) * ymaxscale)
0111     else:
0112         hist.GetYaxis().SetRangeUser(0., (hist.GetMaximum()) * ymaxscale)
0113     hist.GetXaxis().SetTitle(XaxisName)
0114     hist.GetXaxis().SetTickSize(TickSize)
0115     hist.GetXaxis().SetTitleSize(AxisTitleSize)
0116     hist.GetXaxis().SetLabelSize(AxisLabelSize)
0117     hist.GetYaxis().SetTickSize(TickSize)
0118     hist.GetYaxis().SetTitleSize(AxisTitleSize)
0119     hist.GetYaxis().SetLabelSize(AxisLabelSize)
0120     hist.GetXaxis().SetTitleOffset(1.1)
0121     hist.GetYaxis().SetTitleOffset(1.35)
0122     hist.SetLineColor(1)
0123     hist.SetLineWidth(2)
0124     hist.Draw('hist')
0125     histxmin = hist.GetBinLowEdge(hist.FindFirstBinAbove(0))
0126     histxmax = hist.GetBinLowEdge(hist.FindLastBinAbove(0))+binwidth 
0127     f1 = TF1('f1', 'gaus', histxmin, histxmax)
0128     f1.SetParameter(1,hist.GetMean())
0129     f1.SetParLimits(1,-50,0)
0130     f1.SetParameter(2,hist.GetRMS())
0131     f1.SetParLimits(2,0,100)
0132     f1.SetLineColorAlpha(TColor.GetColor('#F54748'), 0.9)
0133     hist.Fit(f1, 'B')
0134     f1.Draw('same')
0135     leg = TLegend((1-RightMargin)-0.45, (1-TopMargin)-0.13,
0136                   (1-RightMargin)-0.1, (1-TopMargin)-0.03)
0137     leg.SetTextSize(0.040)
0138     leg.SetFillStyle(0)
0139     leg.AddEntry('', '#it{#bf{sPHENIX}} Internal', '')
0140     leg.AddEntry('','Au+Au #sqrt{s_{NN}}=200 GeV','')
0141     leg.Draw()
0142     leg2 = TLegend(0.52, 0.65, 0.88, 0.78)
0143     leg2.SetTextSize(0.033)
0144     leg2.SetFillStyle(0)
0145     leg2.AddEntry(f1,'Gaussian fit','l')
0146     leg2.AddEntry('', '#mu={0:.3e}#pm{1:.3e}'.format(f1.GetParameter(1), f1.GetParError(1)), '')
0147     leg2.AddEntry('', '#sigma={0:.3e}#pm{1:.3e}'.format(f1.GetParameter(2), f1.GetParError(2)), '')
0148     leg2.Draw()
0149     c.RedrawAxis()
0150     c.Draw()
0151     c.SaveAs(outname+'.pdf')
0152     c.SaveAs(outname+'.png')
0153     if(c):
0154         c.Close()
0155         gSystem.ProcessEvents()
0156         del c
0157         c = 0
0158     return f1.GetParameter(1), f1.GetParError(1), f1.GetParameter(2), f1.GetParError(2)
0159 
0160 def Draw_2Dhist_wtext(hist, IsData, logz, norm1, rmargin, XaxisName, YaxisName, ZaxisName, addtext, drawopt, outname):
0161     c = TCanvas('c', 'c', 800, 700)
0162     if logz:
0163         c.SetLogz()
0164     c.cd()
0165     if ZaxisName == '':
0166         gPad.SetRightMargin(rmargin)
0167     else:
0168         gPad.SetRightMargin(rmargin+0.03)
0169 
0170     gPad.SetTopMargin(TopMargin)
0171     gPad.SetLeftMargin(LeftMargin)
0172     gPad.SetBottomMargin(BottomMargin)
0173     if norm1:
0174         hist.Scale(1. / hist.Integral(-1, -1, -1, -1))
0175     hist.GetXaxis().SetTitle(XaxisName)
0176     hist.GetYaxis().SetTitle(YaxisName)
0177     hist.GetXaxis().SetTickSize(TickSize)
0178     hist.GetYaxis().SetTickSize(TickSize)
0179     hist.GetXaxis().SetTitleSize(AxisTitleSize)
0180     hist.GetYaxis().SetTitleSize(AxisTitleSize)
0181     hist.GetXaxis().SetLabelSize(AxisLabelSize)
0182     hist.GetYaxis().SetLabelSize(AxisLabelSize)
0183     if ZaxisName != '':
0184         hist.GetZaxis().SetTitle(ZaxisName)
0185         hist.GetZaxis().SetTitleSize(AxisTitleSize)
0186         hist.GetZaxis().SetTitleOffset(1.2*(rmargin/0.14))
0187         
0188     hist.GetXaxis().SetTitleOffset(1.1)
0189     hist.GetYaxis().SetTitleOffset(1.3)
0190     hist.GetZaxis().SetLabelSize(AxisLabelSize)
0191     hist.SetContour(1000)
0192     hist.Draw(drawopt)
0193 
0194     # rightshift = 0.09 if IsData else 0.1
0195     rightshift = 0
0196     leg = TLegend((1-RightMargin)-0.2, (1-TopMargin)+0.01, 1-gPad.GetRightMargin(), (1-TopMargin)+0.04)
0197     leg.SetTextAlign(kHAlignRight+kVAlignBottom)
0198     leg.SetTextSize(0.045)
0199     leg.SetFillStyle(0)
0200     if IsData:
0201         leg.AddEntry("", "#it{#bf{sPHENIX}} Internal", "")
0202         # leg.AddEntry("", "Au+Au #sqrt{s_{NN}}=200 GeV", "")
0203     else:
0204         leg.AddEntry("", "#it{#bf{sPHENIX}} Simulation", "")
0205         # leg.AddEntry("", "Au+Au #sqrt{s_{NN}}=200 GeV", "")
0206     leg.Draw()
0207     c.RedrawAxis()
0208     # add text
0209     if addtext != '':
0210         tex = TLatex()
0211         tex.SetNDC()
0212         tex.SetTextSize(0.04)
0213         tex.DrawLatex(0.2, 0.85, addtext)
0214     c.Draw()
0215     c.SaveAs(outname+'.pdf')
0216     c.SaveAs(outname+'.png')
0217     if(c):
0218         c.Close()
0219         gSystem.ProcessEvents()
0220         del c
0221         c = 0
0222         
0223         
0224 if __name__ == '__main__':
0225     parser = OptionParser(usage="usage: %prog ver [options -h]")
0226     parser.add_option("-i", "--infiledir", dest="infiledir", default='/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/hists/Data_CombinedNtuple_Run20869_HotDead_BCO_ADC_Survey/RecoVtx', help="Input file directory")
0227     parser.add_option("-o", "--outdirprefix", dest="outdirprefix", default='Data_Run20869_HotDead_BCO_ADC_Survey', help="Output file directory")
0228 
0229     (opt, args) = parser.parse_args()
0230     print('opt: {}'.format(opt)) 
0231     
0232     infiledir = opt.infiledir
0233     outdirprefix = opt.outdirprefix
0234     plotpath = './RecoPV_ana/'
0235     os.makedirs('{}/{}'.format(plotpath,outdirprefix), exist_ok=True)
0236     
0237     if os.path.isfile("{}/hists_merged.root".format(infiledir)):
0238         os.system("rm {}/hists_merged.root".format(infiledir))
0239         os.system("hadd -j 20 -f {}/hists_merged.root {}/hists_*.root".format(infiledir, infiledir))
0240     else:
0241         os.system("hadd -j 20 -f {}/hists_merged.root {}/hists_*.root".format(infiledir, infiledir))
0242     
0243     hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10 = GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10')
0244     hM_INTTVtxZ_MBDAsymm_Centrality0to70_Inclusive = GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_MBDAsymm_Centrality0to70_Inclusive')
0245     hM_INTTVtxZ_MBDAsymm_Centrality0to70_MBDAsymLe1_VtxZm10to10 = GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_MBDAsymm_Centrality0to70_MBDAsymLe1_VtxZm10to10')
0246     hM_INTTVtxZ_MBDVtxZ_Centrality0to70_MBDAsymLe1 = GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_MBDVtxZ_Centrality0to70_MBDAsymLe1')
0247     
0248     centrality_cut = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
0249     l_hM_INTTVtxZ_Centrality = [GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_Centrality_{:d}to{:d}'.format(centrality_cut[i], centrality_cut[i+1])) for i in range(len(centrality_cut)-1)]
0250     l_hM_INTTVtxZ_Centrality_MBDAsymLe1 = [GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_Centrality_{:d}to{:d}_MBDAsymLe1'.format(centrality_cut[i], centrality_cut[i+1])) for i in range(len(centrality_cut)-1)]
0251     l_hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10 = [GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_Centrality_{:d}to{:d}_MBDAsymLe1_VtxZm10to10'.format(centrality_cut[i], centrality_cut[i+1])) for i in range(len(centrality_cut)-1)]
0252     l_hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10_forfit = [GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_Centrality_{:d}to{:d}_MBDAsymLe1_VtxZm10to10'.format(centrality_cut[i], centrality_cut[i+1])) for i in range(len(centrality_cut)-1)]
0253     l_hM_INTTVtxZ_MBDAsymm_Centrality_MBDAsymLe1_VtxZm10to10 = [GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_MBDAsymm_Centrality_{:d}to{:d}_MBDAsymLe1_VtxZm10to10'.format(centrality_cut[i], centrality_cut[i+1])) for i in range(len(centrality_cut)-1)]
0254     
0255     colors = ['#001219','#013a63', '#005f73', '#0a9396', '#94d2bd', '#ee9b00', '#ca6702', '#bb3e03', '#9b2226', '#7b193a']
0256     # drawcomparison_data(hists, normhist, logy, histcolor, legentry, legposition, xtitle, ytitle, yscale, outname)
0257     drawcomparison_data(l_hM_INTTVtxZ_Centrality, l_hM_INTTVtxZ_Centrality[0], True, colors, ['Centrality 0-10%', 'Centrality 10-20%', 'Centrality 20-30%', 'Centrality 30-40%', 'Centrality 40-50%', 'Centrality 50-60%', 'Centrality 60-70%', 'Centrality 70-80%', 'Centrality 80-90%', 'Centrality 90-100%'], [0.17, 0.6, 0.42, 0.9], 'INTT vtx_{Z} [cm]', 'Normalized entries', 10, '{}/{}/InttVtxz_MBDCentrality_Comp'.format(plotpath, outdirprefix))
0258     drawcomparison_data(l_hM_INTTVtxZ_Centrality_MBDAsymLe1, l_hM_INTTVtxZ_Centrality_MBDAsymLe1[0], True, colors, ['Centrality 0-10%', 'Centrality 10-20%', 'Centrality 20-30%', 'Centrality 30-40%', 'Centrality 40-50%', 'Centrality 50-60%', 'Centrality 60-70%', 'Centrality 70-80%', 'Centrality 80-90%', 'Centrality 90-100%'], [0.17, 0.6, 0.42, 0.9], 'INTT vtx_{Z} [cm]', 'Normalized entries', 10, '{}/{}/InttVtxz_MBDCentrality_MBDAsymLe1_Comp'.format(plotpath, outdirprefix)) 
0259     drawcomparison_data(l_hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10, l_hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10[0], False, colors, ['Centrality 0-10%', 'Centrality 10-20%', 'Centrality 20-30%', 'Centrality 30-40%', 'Centrality 40-50%', 'Centrality 50-60%', 'Centrality 60-70%', 'Centrality 70-80%', 'Centrality 80-90%', 'Centrality 90-100%'], [0.17, 0.6, 0.42, 0.9], 'INTT vtx_{Z} [cm]', 'Normalized entries', 1.2, '{}/{}/InttVtxz_MBDCentrality_MBDAsymLe1_VtxZm10to10_Comp'.format(plotpath, outdirprefix))
0260     
0261     l_hM_INTTVtxZ_MBDVtxZ_Centrality = [GetHistogram('{}/hists_merged.root'.format(infiledir), 'hM_INTTVtxZ_MBDVtxZ_Centrality_{:d}to{:d}'.format(centrality_cut[i], centrality_cut[i+1])) for i in range(len(centrality_cut)-1)]
0262     # Draw_2Dhist(hist, IsData, logz, norm1, rmargin, XaxisName, YaxisName, ZaxisName, drawopt, outname):
0263     # Draw_2Dhist_wtext(hist, IsData, logz, norm1, rmargin, XaxisName, YaxisName, ZaxisName, addtext, drawopt, outname):
0264     Draw_2Dhist_wtext(hM_INTTVtxZ_MBDVtxZ_Centrality0to70_MBDAsymLe1, True, False, False, 0.18, 'INTT vtx_{Z} [cm]', 'MBD vtx_{Z} [cm]', 'Entries', "Centrality: 0-70%", 'colz', '{}/{}/InttVtxz_MbdVtxz_MBDCentrality0to70_MBDAsymLe1'.format(plotpath, outdirprefix))
0265     Draw_2Dhist_wtext(hM_INTTVtxZ_MBDAsymm_Centrality0to70_Inclusive, True, False, False, 0.15, 'INTT vtx_{Z} [cm]', 'MBD charge asymmetry', 'Entries', "Centrality: 0-70%", 'colz', '{}/{}/InttVtxz_MbdAsymm_MBDCentrality0to70'.format(plotpath, outdirprefix))
0266     Draw_2Dhist_wtext(hM_INTTVtxZ_MBDAsymm_Centrality0to70_MBDAsymLe1_VtxZm10to10, True, False, False, 0.15, 'INTT vtx_{Z} [cm]', 'MBD charge asymmetry', 'Entries', "Centrality: 0-70%", 'colz', '{}/{}/InttVtxz_MbdAsymm_MBDCentrality0to70_MBDAsymLe1_VtxZm10to10'.format(plotpath, outdirprefix))
0267     
0268     for i in range(len(centrality_cut)-1):
0269         # Draw_1Dhist_fitGaussian(hist, norm1, logy, ymaxscale, XaxisName, Ytitle_unit, outname)
0270         Draw_1Dhist_fitGaussian(l_hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10_forfit[i], False, False, 1.3, 'INTT vtx_{Z} [cm]', 'cm', '{}/{}/InttVtxz_MBDCentrality_MBDAsymLe1_VtxZm10to10_{:d}to{:d}_GaussianFit'.format(plotpath, outdirprefix, centrality_cut[i], centrality_cut[i+1]))
0271         Draw_2Dhist_wtext(l_hM_INTTVtxZ_MBDVtxZ_Centrality[i], True, False, False, 0.15, 'INTT vtx_{Z} [cm]', 'MBD vtx_{Z} [cm]', 'Entries', "Centrality: {}-{}%".format(centrality_cut[i], centrality_cut[i+1]), 'colz', '{}/{}/InttVtxz_MbdVtxz_MBDCentrality{:d}to{:d}'.format(plotpath, outdirprefix, centrality_cut[i], centrality_cut[i+1]))
0272     
0273     Draw_1Dhist_fitGaussian(hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10, False, False, 1.3, 'INTT vtx_{Z} [cm]', 'cm', '{}/{}/InttVtxz_MBDCentrality0to70_MBDAsymLe1_VtxZm10to10_GaussianFit'.format(plotpath, outdirprefix))
[ 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!