analysis_INTT/plot/systematics.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, kBlack, kGreen, kRed, kBlue, kOrange, kViolet, kAzure, kMagenta, kCyan, kYellow, kGray, kWhite, gPad, TPad, TLine
0008 import numpy
0009 import math
0010 import glob
0011 from plotUtil import GetHistogram, markerset, colorset2
0012
0013 TickSize = 0.03
0014 AxisTitleSize = 0.05
0015 AxisLabelSize = 0.045
0016 LeftMargin = 0.15
0017 RightMargin = 0.05
0018 TopMargin = 0.08
0019 BottomMargin = 0.13
0020 textscale = 2.6
0021
0022 gROOT.SetBatch(True)
0023
0024 def getRelativeDiff(h_2diff, h_nominal):
0025     # check if the histograms have the same binning
0026     if h_2diff.GetNbinsX() != h_nominal.GetNbinsX():
0027         print('Error: histograms have different number of bins')
0028         return None
0029
0030     h_reldiff = h_2diff.Clone('h_reldiff')
0031     h_reldiff.Reset()
0032     for i in range(1, h_2diff.GetNbinsX() + 1):
0033         absdiff = abs(h_2diff.GetBinContent(i) - h_nominal.GetBinContent(i))
0034         nominal = h_nominal.GetBinContent(i)
0035         reldiff = (absdiff / nominal) if nominal != 0 else 0
0036         h_reldiff.SetBinContent(i, reldiff)
0037
0038     return h_reldiff
0039
0040
0041 def getMaxRelDiff(l_hm_reldiff):
0042     hm_maxRelDiff = l_hm_reldiff[0].Clone('hm_maxRelDiff')
0043     hm_maxRelDiff.Reset()
0044     for hm_reldiff in l_hm_reldiff:
0045         for i in range(1, hm_reldiff.GetNbinsX() + 1):
0046             if hm_reldiff.GetBinContent(i) > hm_maxRelDiff.GetBinContent(i):
0047                 hm_maxRelDiff.SetBinContent(i, hm_reldiff.GetBinContent(i))
0048
0049     return hm_maxRelDiff
0050
0051
0052 def getFinalUncertainty(l_hm):
0053     # sum in quadrature of the uncertainties
0054     hM_TotalRelUnc = l_hm[0].Clone('hM_TotalRelUnc')
0055     hM_TotalRelUnc.Reset()
0056     for i in range(1, hM_TotalRelUnc.GetNbinsX() + 1):
0057         totalUnc = 0
0058         for hm in l_hm:
0059             totalUnc += hm.GetBinContent(i)**2
0060
0061         totalUnc = math.sqrt(totalUnc)
0062         hM_TotalRelUnc.SetBinContent(i, totalUnc)
0063
0064     return hM_TotalRelUnc
0065
0066 # provide a list of histograms to draw and compare. The first histogram in the list is the baseline
0067 def drawcomp(l_hm, l_hmleg, legheader, plotname):
0068     lmarker = markerset(len(l_hm))
0069     lcolor = colorset2(len(l_hm) - 1)
0070     ms = 1.2
0071     alpha = 0.8
0072     ymin = 1E9
0073     ymax = -1
0074
0075     # take the ratio to the baseline; the baseline is always the first in the list
0076     l_ratio = []
0077     for i in range(0, len(l_hm)):
0078         h_ratio = l_hm[i].Clone('h_ratio_variation{}'.format(i))
0079         h_ratio.Divide(l_hm[0])
0080         l_ratio.append(h_ratio)
0081
0082     for hm in l_hm:
0083         if hm.GetMinimum(0) < ymin:
0084             ymin = hm.GetMinimum(0)
0085         if hm.GetMaximum() > ymax:
0086             ymax = hm.GetMaximum()
0087
0088     # create a dummy histogram for axis settings
0089     hdummy = TH1F('hdummy', 'hdummy', l_hm[0].GetNbinsX(), l_hm[0].GetXaxis().GetXmin(), l_hm[0].GetXaxis().GetXmax())
0090
0091     c = TCanvas('c', 'c', 800, 700)
0092     pad1 = TPad( 'pad1', ' ', 0, 0.3, 1, 1)
0093     pad2 = TPad( 'pad2', ' ', 0, 0, 1, 0.3)
0094     pad1.SetTopMargin(0.06*(len(l_hm)+2))
0095     pad1.SetBottomMargin(0.0)
0096     pad1.Draw()
0097     pad2.SetTopMargin(0.0)
0098     pad2.SetBottomMargin(0.4)
0099     pad2.Draw()
0100     pad1.cd()
0101     # set the axis titles and labels for the histograms
0102     # hdummy.GetXaxis().SetTitle('Pseudorapidity #eta')
0103     hdummy.GetYaxis().SetTitle('dN_{ch}/d#eta')
0104     hdummy.GetYaxis().SetTitleOffset(1.6)
0105     hdummy.GetYaxis().SetNdivisions(505)
0106     hdummy.GetYaxis().SetRangeUser(ymin*0.975, ymax*1.025)
0107     hdummy.Draw()
0108     l_hm[0].SetMarkerStyle(20)
0109     l_hm[0].SetMarkerColor(kBlack)
0110     l_hm[0].SetMarkerSize(ms)
0111     l_hm[0].SetLineColor(kBlack)
0112     l_hm[0].SetLineWidth(0)
0113     l_hm[0].Draw('pex0 same')
0114     for i in range(1, len(l_hm)):
0115         l_hm[i].SetMarkerStyle(lmarker[i])
0116         l_hm[i].SetMarkerColorAlpha(TColor.GetColor(lcolor[i - 1]), alpha)
0117         l_hm[i].SetMarkerSize(ms if lmarker[i] != 33 else ms * 1.4)
0118         l_hm[i].SetLineColorAlpha(TColor.GetColor(lcolor[i - 1]), alpha)
0119         l_hm[i].SetLineWidth(0)
0120         l_hm[i].Draw('pex0 same')
0121     leg = TLegend(gPad.GetLeftMargin(), 1 - gPad.GetTopMargin() + 0.02, gPad.GetLeftMargin() + 0.2, 0.98)
0122     leg.SetHeader(legheader)
0123     leg.SetTextSize(0.045)
0124     leg.SetBorderSize(0)
0125     leg.SetFillStyle(0)
0126     for i in range(len(l_hm)):
0127         leg.AddEntry(l_hm[i], '{}'.format(l_hmleg[i]), 'p')
0128     leg.Draw()
0129     c.RedrawAxis()
0130     c.Update()
0131     # cd to the pad2 and draw the ratio
0132     pad2.cd()
0133     for i in range(1, len(l_ratio)):
0134         if i == 1:
0135             l_ratio[i].GetYaxis().SetTitle('Ratio')
0136             l_ratio[i].GetYaxis().SetTitleOffset(0.6)
0137             l_ratio[i].GetYaxis().SetRangeUser(0.9, 1.1)
0138             l_ratio[i].GetYaxis().SetNdivisions(502)
0139             l_ratio[i].GetYaxis().SetTitleSize(AxisTitleSize * textscale)
0140             l_ratio[i].GetYaxis().SetLabelSize(AxisLabelSize * textscale)
0141             l_ratio[i].GetXaxis().SetTitle('Pseudorapidity #eta')
0142             l_ratio[i].GetXaxis().SetTitleSize(AxisTitleSize * textscale)
0143             l_ratio[i].GetXaxis().SetLabelSize(AxisLabelSize * textscale)
0144             l_ratio[i].GetXaxis().SetTickLength(0.03*textscale)
0145             l_ratio[i].SetMarkerStyle(lmarker[i])
0146             l_ratio[i].SetMarkerColor(TColor.GetColor(lcolor[i - 1]))
0147             l_ratio[i].SetMarkerSize(ms if lmarker[i] != 33 else ms * 1.4)
0148             l_ratio[i].SetLineColor(TColor.GetColor(lcolor[i - 1]))
0149             l_ratio[i].SetLineWidth(0)
0150             l_ratio[i].Draw('pex0')
0151         else:
0152             l_ratio[i].SetMarkerStyle(lmarker[i])
0153             l_ratio[i].SetMarkerColor(TColor.GetColor(lcolor[i - 1]))
0154             l_ratio[i].SetMarkerSize(ms if lmarker[i] != 33 else ms * 1.4)
0155             l_ratio[i].SetLineColor(TColor.GetColor(lcolor[i - 1]))
0156             l_ratio[i].SetLineWidth(0)
0157             l_ratio[i].Draw('pex0 same')
0158     # draw a horizontal line at 1
0159     hline = TLine(l_ratio[1].GetXaxis().GetXmin(), 1, l_ratio[1].GetXaxis().GetXmax(), 1)
0160     hline.SetLineColor(kBlack)
0161     hline.SetLineWidth(1)
0162     hline.SetLineStyle(2)
0163     hline.Draw('same')
0164     c.Update()
0165     c.SaveAs(plotname+'.pdf')
0166     c.SaveAs(plotname+'.png')
0167     if (c):
0168         c.Close()
0169         gSystem.ProcessEvents()
0170         del c
0171         c = 0
0172
0173 if __name__ == '__main__':
0174     parser = OptionParser(usage='usage: %prog ver [options -h]')
0175     parser.add_option('--desc', dest='desc', type='string', default='Centrality0to70_Zvtxm30p0tom10p0_noasel', help='String for input file and output plot description (cuts, binnings, etc...)')
0176
0177     (opt, args) = parser.parse_args()
0178     print('opt: {}'.format(opt))
0179
0180     desc = opt.desc
0181
0182     # deltaR selection - baseline 0.5, variation  0.4, 0.6
0183     # cluster ADC cut - baseline 35, variation 50 and 0
0184     # cluster phi-size cut - baseline 40, variation none
0185     # random noise cluster
0186     # Statistical uncertainties in the geometric correction and alpha corrections as propagated to the final dNdEta
0187
0188     str_centrality = desc.split('_')[0]
0189     centrality_low = str_centrality.split('to')[0].replace('Centrality', '')
0190     centrality_high = str_centrality.split('to')[1]
0191     str_zvtx = desc.split('_')[1]
0192     zvtx_low = str_zvtx.split('to')[0].replace('Zvtx', '').replace('m', '-').replace('p', '.') if 'm' in str_zvtx else str_zvtx.split('to')[0].replace('Zvtx', '').replace('p', '.')
0193     zvtx_high = str_zvtx.split('to')[1].replace('m', '-').replace('p', '.') if 'm' in str_zvtx else str_zvtx.split('to')[1].replace('p', '.')
0194     str_asel = desc.split('_')[2]
0195     print (str_centrality, str_zvtx, str_asel)
0196     print ('centrality_low: {}, centrality_high: {}, zvtx_low: {}, zvtx_high: {}, str_asel: {}'.format(centrality_low, centrality_high, zvtx_low, zvtx_high, str_asel))
0197     alt_zvtx1 = 'Zvtxm10p0to0p0'
0198     alt_zvtx2 = 'Zvtx0p0to10p0'
0199
0200     os.makedirs('./systematics/{}'.format(desc), exist_ok=True)
0201
0202     h1WEfinal_nominal = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0203     # tracklet delta-R cut
0204     h1WEfinal_dRcutSet1 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p4_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0205     h1WEfinal_dRcutSet2 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p6_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0206     # cluster ADC cut
0207     h1WEfinal_clusAdcCutSet1 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet1_clusPhiSizeCutSet0_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0208     h1WEfinal_clusAdcCutSet2 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet2_clusPhiSizeCutSet0_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0209     # cluster phi-size cut
0210     h1WEfinal_clusphisizeSet1 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet1_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0211     # strangeness enhanced by 40% in HIJING
0212     h1WEfinal_strangeness = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_enhancedstrangeness_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0213     # event generator
0214     h1WEfinal_EPOS = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_EPOS_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0215     h1WEfinal_AMPT = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_AMPT_segment1/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0216     # different segments
0217     h1WEfinal_segment2 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_segment2/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0218     h1WEfinal_segment3 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_segment3/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0219     h1WEfinal_segment4 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_segment4/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0220     h1WEfinal_segment5 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_segment5/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0221     h1WEfinal_segment6 = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/corrections/Data_Run54280_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0_segment6/{}/correction_hists.root'.format(desc), 'h1WEfinal')
0222
0223     # get the uncertainty from the nominal
0224     hM_statunc_nominal = h1WEfinal_nominal.Clone('hM_statunc_nominal')
0225     hM_statunc_nominal.Reset()
0226     for i in range(1, h1WEfinal_nominal.GetNbinsX() + 1):
0227         statunc = h1WEfinal_nominal.GetBinError(i) / h1WEfinal_nominal.GetBinContent(i) if h1WEfinal_nominal.GetBinContent(i) != 0 else 0
0228         hM_statunc_nominal.SetBinContent(i, statunc)
0229
0230     # calculate the relative difference (absolute value of difference divided by the nominal value)
0231     # tracklet delta-R cut
0232     hM_reldiff_h1WEfinal_dRcutSet1 = getRelativeDiff(h1WEfinal_dRcutSet1, h1WEfinal_nominal)
0233     hM_reldiff_h1WEfinal_dRcutSet2 = getRelativeDiff(h1WEfinal_dRcutSet2, h1WEfinal_nominal)
0234     hM_maxreldiff_dRcut = getMaxRelDiff([hM_reldiff_h1WEfinal_dRcutSet1, hM_reldiff_h1WEfinal_dRcutSet2])
0235     # hM_maxreldiff_dRcut = getMaxRelDiff([hM_reldiff_h1WEfinal_dRcutSet1, hM_reldiff_h1WEfinal_dRcutSet2, hM_reldiff_h1WEfinal_dRcutSet3, hM_reldiff_h1WEfinal_dRcutSet4])
0236     # cluster ADC cut
0237     hM_reldiff_h1WEfinal_clusAdcCutSet1 = getRelativeDiff(h1WEfinal_clusAdcCutSet1, h1WEfinal_nominal)
0238     hM_reldiff_h1WEfinal_clusAdcCutSet2 = getRelativeDiff(h1WEfinal_clusAdcCutSet2, h1WEfinal_nominal)
0239     # hM_maxreldiff_clusAdcCut = getMaxRelDiff([hM_reldiff_h1WEfinal_clusAdcCutSet1])
0240     hM_maxreldiff_clusAdcCut = getMaxRelDiff([hM_reldiff_h1WEfinal_clusAdcCutSet1, hM_reldiff_h1WEfinal_clusAdcCutSet2])
0241     # cluster phi-size cut
0242     hM_reldiff_h1WEfinal_clusPhiSizeCutSet1 = getRelativeDiff(h1WEfinal_clusphisizeSet1, h1WEfinal_nominal)
0243     hM_maxreldiff_clusPhiSizeCut = getMaxRelDiff([hM_reldiff_h1WEfinal_clusPhiSizeCutSet1])
0244     # strangeness enhanced by 40%
0245     hM_reldiff_h1WEfinal_strangeness = getRelativeDiff(h1WEfinal_strangeness, h1WEfinal_nominal)
0246     hM_maxreldiff_strangeness = getMaxRelDiff([hM_reldiff_h1WEfinal_strangeness])
0247     # event generator
0248     hM_reldiff_h1WEfinal_EPOS = getRelativeDiff(h1WEfinal_EPOS, h1WEfinal_nominal)
0249     hM_reldiff_h1WEfinal_AMPT = getRelativeDiff(h1WEfinal_AMPT, h1WEfinal_nominal)
0250     hM_maxreldiff_eventgen = getMaxRelDiff([hM_reldiff_h1WEfinal_EPOS, hM_reldiff_h1WEfinal_AMPT])
0251     # different segments
0252     hM_reldiff_h1WEfinal_segment2 = getRelativeDiff(h1WEfinal_segment2, h1WEfinal_nominal)
0253     hM_reldiff_h1WEfinal_segment3 = getRelativeDiff(h1WEfinal_segment3, h1WEfinal_nominal)
0254     hM_reldiff_h1WEfinal_segment4 = getRelativeDiff(h1WEfinal_segment4, h1WEfinal_nominal)
0255     hM_reldiff_h1WEfinal_segment5 = getRelativeDiff(h1WEfinal_segment5, h1WEfinal_nominal)
0256     hM_reldiff_h1WEfinal_segment6 = getRelativeDiff(h1WEfinal_segment6, h1WEfinal_nominal)
0257     hM_maxreldiff_segment = getMaxRelDiff([hM_reldiff_h1WEfinal_segment2, hM_reldiff_h1WEfinal_segment3, hM_reldiff_h1WEfinal_segment4, hM_reldiff_h1WEfinal_segment5, hM_reldiff_h1WEfinal_segment6])
0258
0259     hM_TotalRelUnc = getFinalUncertainty([hM_statunc_nominal,
0260                                           hM_maxreldiff_dRcut,
0261                                           hM_maxreldiff_clusAdcCut,
0262                                           hM_maxreldiff_clusPhiSizeCut,
0263                                           hM_maxreldiff_strangeness,
0264                                           hM_maxreldiff_eventgen,
0265                                           hM_maxreldiff_segment
0266                                         ])
0267
0268     # get the absolute uncertainty
0269     hM_final = h1WEfinal_nominal.Clone('hM_final')
0270     hM_TotalUnc = hM_TotalRelUnc.Clone('hM_TotalUnc')
0271     for i in range(1, hM_TotalUnc.GetNbinsX() + 1):
0272         hM_TotalUnc.SetBinContent(i, hM_TotalRelUnc.GetBinContent(i) * h1WEfinal_nominal.GetBinContent(i))
0273         hM_final.SetBinError(i, hM_TotalRelUnc.GetBinContent(i) * h1WEfinal_nominal.GetBinContent(i))
0274
0275     # only include bins with |eta|<=1.1
0276     for i in range(1, hM_TotalRelUnc.GetNbinsX() + 1):
0277         if abs(hM_TotalRelUnc.GetXaxis().GetBinCenter(i)) > 1.25:
0278             # print ('i: {}, eta: {}'.format(i, hM_TotalRelUnc.GetXaxis().GetBinCenter(i)))
0279             hM_statunc_nominal.SetBinContent(i, 0)
0280             hM_maxreldiff_dRcut.SetBinContent(i, 0)
0281             hM_maxreldiff_clusAdcCut.SetBinContent(i, 0)
0282             hM_maxreldiff_clusPhiSizeCut.SetBinContent(i, 0)
0283             hM_maxreldiff_strangeness.SetBinContent(i, 0)
0284             hM_maxreldiff_eventgen.SetBinContent(i, 0)
0285             hM_maxreldiff_segment.SetBinContent(i, 0)
0286             hM_TotalRelUnc.SetBinContent(i, 0)
0287             hM_TotalUnc.SetBinContent(i, 0)
0288             hM_final.SetBinContent(i, 0)
0289
0290
0291     # Root file with the histograms (with the final uncertainty)
0292     fout = TFile('./systematics/{}/finalhists_systematics_{}.root'.format(desc, desc), 'recreate')
0293     fout.cd()
0294     hM_statunc_nominal.Write()
0295     hM_maxreldiff_dRcut.SetName('hM_maxreldiff_dRcut')
0296     hM_maxreldiff_dRcut.Write()
0297     hM_maxreldiff_clusAdcCut.SetName('hM_maxreldiff_clusAdcCut')
0298     hM_maxreldiff_clusAdcCut.Write()
0299     hM_maxreldiff_clusPhiSizeCut.SetName('hM_maxreldiff_clusPhiSizeCut')
0300     hM_maxreldiff_clusPhiSizeCut.Write()
0301     hM_maxreldiff_strangeness.SetName('hM_maxreldiff_strangeness')
0302     hM_maxreldiff_strangeness.Write()
0303     hM_maxreldiff_eventgen.SetName('hM_maxreldiff_eventgen')
0304     hM_maxreldiff_eventgen.Write()
0305     hM_maxreldiff_segment.SetName('hM_maxreldiff_segment')
0306     hM_maxreldiff_segment.Write()
0307     hM_TotalRelUnc.Write()
0308     hM_TotalUnc.Write()
0309     hM_final.Write()
0310     fout.Close()
0311
0312     # plotting
0313     # split the string by '_' and replace some characters
0314     descstr = desc.split('_')
0315     centstr = descstr[0]
0316     zvtxstr = descstr[1]
0317     aselstr = descstr[2]
0318     centstr = centstr.replace('Centrality', 'Centrality [')
0319     centstr = centstr.replace('to', '-')
0320     centstr = centstr + ']%'
0321     zvtxstr = zvtxstr.replace('m', '-')
0322     zvtxstr = zvtxstr.replace('p', '.')
0323     zvtxstr = zvtxstr.replace('to', ', ')
0324     zvtxstr = zvtxstr.replace('Zvtx', 'Z_{vtx} [')
0325     zvtxstr = zvtxstr + '] cm'
0326     aselstr = '' if aselstr == 'noasel' else ', (with additional selection)'
0327     legstr = centstr + ', ' + zvtxstr + aselstr
0328
0329     # drawcomp(l_hm, l_hmleg, legheader, plotname):
0330     drawcomp([h1WEfinal_nominal, h1WEfinal_dRcutSet1, h1WEfinal_dRcutSet2],
0331              ['Nominal (#DeltaR<0.5)', '#DeltaR<0.4', '#DeltaR<0.6'],
0332              '#DeltaR cut variation, ' + centstr, './systematics/{}/h1WEfinal_dRcut_variation_{}'.format(desc, desc))
0333     drawcomp([h1WEfinal_nominal, h1WEfinal_clusAdcCutSet1, h1WEfinal_clusAdcCutSet2],
0334              ['Nominal (Cluster ADC>35)', 'w.o cluster ADC cut', 'Cluster ADC>50'],
0335              'Cluster ADC cut variation, ' + centstr, './systematics/{}/h1WEfinal_clusAdcCut_variation_{}'.format(desc, desc))
0336     drawcomp([h1WEfinal_nominal, h1WEfinal_clusphisizeSet1],
0337              ['Nominal (Cluster #phi-size cut#leq40)', 'w.o cluster #phi-size cut'],
0338              'Cluster #phi-size cut variation, ' + centstr, './systematics/{}/h1WEfinal_clusPhiSizeCut_variation_{}'.format(desc, desc))
0339     drawcomp([h1WEfinal_nominal, h1WEfinal_strangeness],
0340              ['Nominal', 'Enhanced strangeness by 40%'],
0341              'Enhanced strangeness, ' + centstr, './systematics/{}/h1WEfinal_strangeness_variation_{}'.format(desc, desc))
0342     drawcomp([h1WEfinal_nominal, h1WEfinal_EPOS, h1WEfinal_AMPT],
0343              ['Nominal (HIJING)', 'EPOS', 'AMPT'],
0344              'Event generator variation, ' + centstr, './systematics/{}/h1WEfinal_eventgen_variation_{}'.format(desc, desc))
0345     drawcomp([h1WEfinal_nominal, h1WEfinal_segment2, h1WEfinal_segment3, h1WEfinal_segment4, h1WEfinal_segment5, h1WEfinal_segment6],
0346              ['Nominal', 'Segment 2', 'Segment 3', 'Segment 4', 'Segment 5', 'Segment 6'],
0347              'Run segment variation, ' + centstr, './systematics/{}/h1WEfinal_segment_variation_{}'.format(desc, desc))
0348
0349     # scale the histograms by 100 to get the percentage
0350     hM_TotalRelUnc.Scale(100)
0351     hM_statunc_nominal.Scale(100)
0352     hM_maxreldiff_dRcut.Scale(100)
0353     hM_maxreldiff_clusAdcCut.Scale(100)
0354     hM_maxreldiff_clusPhiSizeCut.Scale(100)
0355     hM_maxreldiff_strangeness.Scale(100)
0356     hM_maxreldiff_eventgen.Scale(100)
0357     hM_maxreldiff_segment.Scale(100)
0358
0359     # plot the relative difference
0360     markersty = 20
0361     c = TCanvas('c', 'c', 800, 700)
0362     c.cd()
0363     hM_TotalRelUnc.GetXaxis().SetTitle('#eta')
0364     hM_TotalRelUnc.GetYaxis().SetTitle('Relative uncertainty [%]')
0365     hM_TotalRelUnc.GetYaxis().SetTitleOffset(1.6)
0366     hM_TotalRelUnc.SetLineColor(kBlack)
0367     hM_TotalRelUnc.SetLineWidth(2)
0368     hM_TotalRelUnc.GetYaxis().SetRangeUser(0, hM_TotalRelUnc.GetMaximum() * 2.0)
0369     hM_TotalRelUnc.Draw('hist')
0370     hM_statunc_nominal.SetMarkerStyle(20)
0371     hM_statunc_nominal.SetMarkerColor(TColor.GetColor('#0F4C75'))
0372     hM_statunc_nominal.SetLineColor(TColor.GetColor('#0F4C75'))
0373     hM_statunc_nominal.SetLineWidth(3)
0374     hM_statunc_nominal.Draw('P same')
0375     hM_maxreldiff_dRcut.SetMarkerStyle(21)
0376     hM_maxreldiff_dRcut.SetMarkerColor(TColor.GetColor('#99ddff'))
0377     hM_maxreldiff_dRcut.SetLineColor(TColor.GetColor('#99ddff'))
0378     hM_maxreldiff_dRcut.SetLineWidth(3)
0379     hM_maxreldiff_dRcut.Draw('P same')
0380     hM_maxreldiff_clusAdcCut.SetMarkerStyle(33)
0381     hM_maxreldiff_clusAdcCut.SetMarkerSize(1.6)
0382     hM_maxreldiff_clusAdcCut.SetMarkerColor(TColor.GetColor('#44bb99'))
0383     hM_maxreldiff_clusAdcCut.SetLineColor(TColor.GetColor('#44bb99'))
0384     hM_maxreldiff_clusAdcCut.SetLineWidth(3)
0385     hM_maxreldiff_clusAdcCut.Draw('P same')
0386     hM_maxreldiff_clusPhiSizeCut.SetMarkerStyle(markersty)
0387     hM_maxreldiff_clusPhiSizeCut.SetMarkerColor(TColor.GetColor('#aadd00'))
0388     hM_maxreldiff_clusPhiSizeCut.SetLineColor(TColor.GetColor('#aadd00'))
0389     hM_maxreldiff_clusPhiSizeCut.SetLineWidth(2)
0390     hM_maxreldiff_clusPhiSizeCut.Draw('P same')
0391     hM_maxreldiff_strangeness.SetMarkerStyle(34)
0392     hM_maxreldiff_strangeness.SetMarkerColor(TColor.GetColor('#e99960'))
0393     hM_maxreldiff_strangeness.SetLineColor(TColor.GetColor('#e99960'))
0394     hM_maxreldiff_strangeness.SetLineWidth(3)
0395     hM_maxreldiff_strangeness.Draw('P same')
0396     hM_maxreldiff_eventgen.SetMarkerStyle(22)
0397     hM_maxreldiff_eventgen.SetMarkerColor(TColor.GetColor('#DE7C7D'))
0398     hM_maxreldiff_eventgen.SetLineColor(TColor.GetColor('#DE7C7D'))
0399     hM_maxreldiff_eventgen.SetLineWidth(3)
0400     hM_maxreldiff_eventgen.Draw('P same')
0401     hM_maxreldiff_segment.SetMarkerStyle(47)
0402     hM_maxreldiff_segment.SetMarkerColor(TColor.GetColor('#AA60C8'))
0403     hM_maxreldiff_segment.SetLineColor(TColor.GetColor('#AA60C8'))
0404     hM_maxreldiff_segment.SetLineWidth(3)
0405     hM_maxreldiff_segment.Draw('P same')
0406     leg = TLegend(0.2, 0.6, 0.45, 0.91)
0407     leg.SetHeader(legstr)
0408     leg.SetTextSize(0.035)
0409     leg.SetBorderSize(0)
0410     leg.SetFillStyle(0)
0411     leg.AddEntry(hM_TotalRelUnc, 'Total uncertainty', 'l')
0412     leg.AddEntry(hM_statunc_nominal, 'Statistical uncertainty in corrections', 'p')
0413     leg.AddEntry(hM_maxreldiff_dRcut, 'Tracklet #DeltaR cut variation', 'p')
0414     leg.AddEntry(hM_maxreldiff_clusAdcCut, 'Cluster ADC cut variation', 'p')
0415     leg.AddEntry(hM_maxreldiff_clusPhiSizeCut, 'Cluster #phi-size cut variation', 'p')
0416     leg.AddEntry(hM_maxreldiff_strangeness, 'Strangeness variation', 'p')
0417     leg.AddEntry(hM_maxreldiff_eventgen, 'Event generator variation', 'p')
0418     leg.AddEntry(hM_maxreldiff_segment, 'Run segment variation', 'p')
0419     leg.Draw()
0420     c.RedrawAxis()
0421     c.Draw()
0422     c.SaveAs('./systematics/{}/hM_TotalRelUnc_{}.pdf'.format(desc, desc))
0423     c.SaveAs('./systematics/{}/hM_TotalRelUnc_{}.png'.format(desc, desc))
0424     if(c):
0425         c.Close()
0426         gSystem.ProcessEvents()
0427         del c
0428         c = 0