analysis_INTT/plot/Final.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, TPad, gPad, TLegend, TColor, gROOT, TGraphErrors, TGraphAsymmErrors, TGaxis, kHAlignLeft, kVAlignTop, kHAlignRight, kVAlignBottom
0008 import numpy
0009 import math
0010 import glob
0011 from array import array
0012 import json
0013 from plotUtil import GetHistogram
0014 from measurements.measurements import *
0015 import numpy as np
0016 from closure import GetMbinNum
0017
0018 gROOT.SetBatch(True)
0019
0020 TickSize = 0.03
0021 AxisTitleSize = 0.05
0022 AxisLabelSize = 0.04
0023 LeftMargin = 0.22
0024 RightMargin = 0.06
0025 TopMargin = 0.05
0026 BottomMargin = 0.13
0027 FillAlpha = 0.5
0028
0029 include_rawhijing = False
0030
0031 # compute uncertainty - weighted average (! valid only for uncorrelated uncertainties, for example RHIC combination !)
0032 def weighted_average(values, uncertainties):
0033     """
0034     Computes the weighted average and its uncertainty for a given set of values and uncertainties.
0035
0036     Parameters:
0037     values (list or array): Data points.
0038     uncertainties (list or array): Corresponding uncertainties of the data points.
0039
0040     Returns:
0041     tuple: Weighted average and its uncertainty.
0042     """
0043     values = np.array(values)
0044     uncertainties = np.array(uncertainties)
0045
0046     # Compute weights
0047     weights = 1 / uncertainties**2
0048
0049     # Compute weighted average
0050     weighted_avg = np.sum(weights * values) / np.sum(weights)
0051
0052     # Compute uncertainty in weighted average
0053     weighted_uncertainty = np.sqrt(1 / np.sum(weights))
0054
0055     return weighted_avg, weighted_uncertainty
0056
0057 def average(values, uncertainties):
0058     values = np.array(values)
0059     uncertainties = np.array(uncertainties)
0060
0061     avg = np.mean(values)
0062     avg_err = np.mean(uncertainties)
0063
0064     return avg, avg_err
0065
0066
0067 def msrmnt_sphenix_generator_auau_200gev(fnameprefix):
0068     centralitybin = sphenix_centrality_interval()
0069     centnpart, centnparterr = sphenix_centralitynpart()
0070
0071     dNdEta_eta0, dNdEta_eta0_divnpart2, Centrality, dNdEta_eta0_err, dNdEta_eta0_divnpart2_err, Centrality_err = [], [], [], [], [], []
0072
0073     for i in range(len(centralitybin) - 1):
0074         # hM_dNdEtafinal = GetHistogram('./corrections/HIJING_closure_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0/Centrality{}to{}_Zvtxm10p0to10p0_noasel/correction_hists.root'.format(centralitybin[i],centralitybin[i+1]), 'h1WGhadron')
0075         hM_dNdEtafinal = GetHistogram('./corrections/{}/Centrality{}to{}_Zvtxm10p0to10p0_noasel/correction_hists.root'.format(fnameprefix,centralitybin[i],centralitybin[i+1]), 'h1WGhadron')
0076
0077         dNdEta_eta0_centrality = (hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(0)) + hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(-0.2)) + hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(0.2))) / 3.
0078         list_dNdEta_eta0 = [hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(0)), hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(-0.2)), hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(0.2))]
0079         list_dNdEta_eta0_err = [hM_dNdEtafinal.GetBinError(hM_dNdEtafinal.FindBin(0)), hM_dNdEtafinal.GetBinError(hM_dNdEtafinal.FindBin(-0.2)), hM_dNdEtafinal.GetBinError(hM_dNdEtafinal.FindBin(0.2))]
0080         dNdEta_eta0_centrality, dNdEta_eta0_err_centrality = average(list_dNdEta_eta0, list_dNdEta_eta0_err)
0081         # print ('dNdEta_eta0_centrality={}, dNdta_eta0_err_centrality={}'.format(dNdEta_eta0_centrality, dNdEta_eta0_err_centrality))
0082         rel_err_dNdEta = dNdEta_eta0_err_centrality / dNdEta_eta0_centrality
0083         rel_err_npart = centnparterr[i] / centnpart[i]
0084         dNdEta_eta0_divnpart2_err_centrality = dNdEta_eta0_centrality / (centnpart[i] / 2.) * math.sqrt(rel_err_dNdEta**2 + rel_err_npart**2) # error propagation; TODO: include the uncertainty for <npart>
0085
0086         dNdEta_eta0.append(dNdEta_eta0_centrality)
0087         dNdEta_eta0_divnpart2.append(dNdEta_eta0_centrality / (centnpart[i] / 2.))
0088         Centrality.append((centralitybin[i] + centralitybin[i+1]) / 2.)
0089         dNdEta_eta0_err.append(dNdEta_eta0_err_centrality)
0090         dNdEta_eta0_divnpart2_err.append(0)
0091         Centrality_err.append(0)
0092
0093
0094     dNdEta_eta0 = array('f', dNdEta_eta0)
0095     dNdEta_eta0_divnpart2 = array('f', dNdEta_eta0_divnpart2)
0096     Centrality = array('f', Centrality)
0097     dNdEta_eta0_err = array('f', dNdEta_eta0_err)
0098     dNdEta_eta0_divnpart2_err = array('f', dNdEta_eta0_divnpart2_err)
0099     Centrality_err = array('f', Centrality_err)
0100     centnpart = array('f', centnpart)
0101     centnparterr = array('f', centnparterr)
0102     gsphenix_rawhijing_auau_200gev = TGraphErrors(len(Centrality), Centrality, dNdEta_eta0, Centrality_err, dNdEta_eta0_err)
0103     gsphenix_rawhijing_divnpart2_auau_200gev = TGraphErrors(len(centnpart), centnpart, dNdEta_eta0_divnpart2, centnparterr, dNdEta_eta0_divnpart2_err)
0104
0105     return gsphenix_rawhijing_auau_200gev, gsphenix_rawhijing_divnpart2_auau_200gev
0106
0107 def msrmnt_generator_dNdetaNpart2_auau_200gev(filepath):
0108     centnpart, centnparterr = sphenix_centralitynpart()
0109     npart_max = max(centnpart)
0110     npart_min = min(centnpart)
0111
0112     file = TFile(filepath, 'READ')
0113     g_truthdNdEta = file.Get('g_truthdNdEta')
0114
0115     # exclude points outside npart_min and npart_max
0116     # Iterate in reverse order to avoid index issues when removing points
0117     for i in range(g_truthdNdEta.GetN()-1, -1, -1):
0118         if g_truthdNdEta.GetX()[i] < npart_min or g_truthdNdEta.GetX()[i] > npart_max:
0119             g_truthdNdEta.RemovePoint(i)
0120
0121     return g_truthdNdEta
0122
0123
0124 def msrmnt_sphenix_data_auau_200gev(approach='cms'): # approach: 'cms', 'phobos', 'combined'
0125     centralitybin = sphenix_centrality_interval()
0126     centnpart, centnparterr = sphenix_centralitynpart()
0127
0128     dNdEta_eta0, dNdEta_eta0_divnpart2, Centrality, dNdEta_eta0_err, dNdEta_eta0_divnpart2_err, Centrality_err = [], [], [], [], [], []
0129
0130     for i in range(len(centralitybin) - 1):
0131         if approach == 'cms':
0132             hM_dNdEtafinal = GetHistogram('./systematics/Centrality{}to{}_Zvtxm10p0to10p0_noasel/finalhists_systematics_Centrality{}to{}_Zvtxm10p0to10p0_noasel.root'.format(centralitybin[i],centralitybin[i+1],centralitybin[i],centralitybin[i+1]), 'hM_final')
0133         elif approach == 'phobos':
0134             mbin = GetMbinNum('Centrality{}to{}'.format(centralitybin[i],centralitybin[i+1]))
0135             hM_dNdEtafinal = GetHistogram('/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/Run7/EvtVtxZ/FinalResult_10cm_Pol2BkgFit_DeltaPhi0p026/completed/vtxZ_-10_10cm_MBin{}/Final_Mbin{}_00054280/Final_Mbin{}_00054280.root'.format(mbin,mbin,mbin), 'h1D_dNdEta_reco')
0136         elif approach == 'combined':
0137             # hM_dNdEtafinal_cms = GetHistogram('./systematics/Centrality{}to{}_Zvtxm10p0to10p0_noasel/finalhists_systematics_Centrality{}to{}_Zvtxm10p0to10p0_noasel.root'.format(centralitybin[i],centralitybin[i+1],centralitybin[i],centralitybin[i+1]), 'hM_final')
0138             # mbin = GetMbinNum('Centrality{}to{}'.format(centralitybin[i],centralitybin[i+1]))
0139             # hM_dNdEtafinal_phobos = GetHistogram('/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/Run5/EvtVtxZ/FinalResult/completed/vtxZ_-10_10cm_MBin{}/Final_Mbin{}_00054280/Final_Mbin{}_00054280.root'.format(mbin,mbin,mbin), 'h1D_dNdEta_reco')
0140             fcombined = TFile('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/systematics/combined/combined.root')
0141             tgae_dNdEtafinal = fcombined.Get('tgae_combine_Centrality{}to{}'.format(centralitybin[i],centralitybin[i+1]))
0142             hM_dNdEtafinal = GetHistogram('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/systematics/combined/combined.root', 'h_combined_Centrality{}to{}'.format(centralitybin[i],centralitybin[i+1]))
0143         else:
0144             sys.exit('Approach {} is not supported. Use default approach of cms'.format(approach))
0145
0146         if approach == 'cms' or approach == 'phobos':
0147             dNdEta_eta0_centrality = (hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(0)) + hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(-0.2)) + hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(0.2))) / 3.
0148             list_dNdEta_eta0 = [hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(0)), hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(-0.2)), hM_dNdEtafinal.GetBinContent(hM_dNdEtafinal.FindBin(0.2))]
0149             list_dNdEta_eta0_err = [hM_dNdEtafinal.GetBinError(hM_dNdEtafinal.FindBin(0)), hM_dNdEtafinal.GetBinError(hM_dNdEtafinal.FindBin(-0.2)), hM_dNdEtafinal.GetBinError(hM_dNdEtafinal.FindBin(0.2))]
0150             dNdEta_eta0_centrality, dNdEta_eta0_err_centrality = average(list_dNdEta_eta0, list_dNdEta_eta0_err)
0151             # print ('dNdEta_eta0_centrality={}, dNdta_eta0_err_centrality={}'.format(dNdEta_eta0_centrality, dNdEta_eta0_err_centrality))
0152             rel_err_dNdEta = dNdEta_eta0_err_centrality / dNdEta_eta0_centrality
0153             rel_err_npart = centnparterr[i] / centnpart[i]
0154             dNdEta_eta0_divnpart2_err_centrality = dNdEta_eta0_centrality / (centnpart[i] / 2.) * math.sqrt(rel_err_dNdEta**2 + rel_err_npart**2)
0155
0156             # print('Approach {}: centrality {}-{}, dNdeta +- uncertainty = {:.3f} +- {:.3f} ({:.3f}%), dNdeta/(0.5 * <npart>) +- uncertainty = {:.3f} +- {:.3f} ({:.3f}%)'.format(approach,
0157             #                                                                                                                                                                      centralitybin[i],
0158             #                                                                                                                                                                      centralitybin[i+1],
0159             #                                                                                                                                                                      dNdEta_eta0_centrality,
0160             #                                                                                                                                                                      dNdEta_eta0_err_centrality,
0161             #                                                                                                                                                                      dNdEta_eta0_err_centrality / dNdEta_eta0_centrality * 100,
0162             #                                                                                                                                                                      dNdEta_eta0_centrality / (centnpart[i] / 2.),
0163             #                                                                                                                                                                      dNdEta_eta0_divnpart2_err_centrality,
0164             #                                                                                                                                                                      dNdEta_eta0_divnpart2_err_centrality / (dNdEta_eta0_centrality / (centnpart[i] / 2.)) * 100))
0165
0166             dNdEta_eta0.append(dNdEta_eta0_centrality)
0167             dNdEta_eta0_divnpart2.append(dNdEta_eta0_centrality / (centnpart[i] / 2.))
0168             Centrality.append((centralitybin[i] + centralitybin[i+1]) / 2.)
0169             dNdEta_eta0_err.append(dNdEta_eta0_err_centrality)
0170             dNdEta_eta0_divnpart2_err.append(dNdEta_eta0_divnpart2_err_centrality)
0171             Centrality_err.append(0)
0172         elif approach == 'combined':
0173             # use the histogram to get the indices of the bins, then use the indices to get the values and errors from the TGraphAsymmErrors
0174             # print ('Eta = 0, dNdEta combined = {:.3f} +- {:.3f}'.format(tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(0)-1], tgae_dNdEtafinal.GetEYlow()[hM_dNdEtafinal.GetXaxis().FindBin(0)-1]))
0175             # print ('Eta = -0.2, dNdEta combined = {:.3f} +- {:.3f}'.format(tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(-0.2)-1], tgae_dNdEtafinal.GetEYlow()[hM_dNdEtafinal.GetXaxis().FindBin(-0.2)-1]))
0176             # print ('Eta = 0.2, dNdEta combined = {:.3f} +- {:.3f}'.format(tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(0.2)-1], tgae_dNdEtafinal.GetEYlow()[hM_dNdEtafinal.GetXaxis().FindBin(0.2)-1]))
0177             dNdEta_eta0_centrality = (tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(0)-1] + tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(-0.2)-1] + tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(0.2)-1]) / 3.
0178             list_dNdEta_eta0 = [tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(0)-1], tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(-0.2)-1], tgae_dNdEtafinal.GetY()[hM_dNdEtafinal.GetXaxis().FindBin(0.2)-1]]
0179             list_dNdEta_eta0_err = [tgae_dNdEtafinal.GetEYlow()[hM_dNdEtafinal.GetXaxis().FindBin(0)-1], tgae_dNdEtafinal.GetEYlow()[hM_dNdEtafinal.GetXaxis().FindBin(-0.2)-1], tgae_dNdEtafinal.GetEYlow()[hM_dNdEtafinal.GetXaxis().FindBin(0.2)-1]]
0180             dNdEta_eta0_centrality, dNdEta_eta0_err_centrality = average(list_dNdEta_eta0, list_dNdEta_eta0_err)
0181             rel_err_dNdEta = dNdEta_eta0_err_centrality / dNdEta_eta0_centrality
0182             rel_err_npart = centnparterr[i] / centnpart[i]
0183             dNdEta_eta0_divnpart2_err_centrality = dNdEta_eta0_centrality / (centnpart[i] / 2.) * math.sqrt(rel_err_dNdEta**2 + rel_err_npart**2)
0184             # print in the Latex table format
0185             precision = '.2f'  # Default precision, can be changed to '.2f' or any other format
0186             # print ('{}\\%-{}\\% & ${:{}}\\pm{:{}}$ & ${:{}}\\pm{:{}}$ & ${:{}}\\pm{:{}}$ \\\\'.format(
0187             #     centralitybin[i], centralitybin[i+1],
0188             #     dNdEta_eta0_centrality, precision, dNdEta_eta0_err_centrality, precision,
0189             #     centnpart[i], precision, centnparterr[i], precision,
0190             #     dNdEta_eta0_centrality / (centnpart[i] / 2.), precision, dNdEta_eta0_divnpart2_err_centrality, precision))
0191
0192             dNdEta_eta0.append(dNdEta_eta0_centrality)
0193             dNdEta_eta0_divnpart2.append(dNdEta_eta0_centrality / (centnpart[i] / 2.))
0194             Centrality.append((centralitybin[i] + centralitybin[i+1]) / 2.)
0195             dNdEta_eta0_err.append(dNdEta_eta0_err_centrality)
0196             dNdEta_eta0_divnpart2_err.append(dNdEta_eta0_divnpart2_err_centrality)
0197             Centrality_err.append(0)
0198         else:
0199             sys.exit('Approach {} is not supported. Use default approach of cms'.format(approach))
0200
0201     dNdEta_eta0 = array('f', dNdEta_eta0)
0202     dNdEta_eta0_divnpart2 = array('f', dNdEta_eta0_divnpart2)
0203     Centrality = array('f', Centrality)
0204     dNdEta_eta0_err = array('f', dNdEta_eta0_err)
0205     dNdEta_eta0_divnpart2_err = array('f', dNdEta_eta0_divnpart2_err)
0206     Centrality_err = array('f', Centrality_err)
0207     centnpart = array('f', centnpart)
0208     centnparterr = array('f', centnparterr)
0209     gsphenix_data_auau_200gev = TGraphAsymmErrors(len(Centrality), Centrality, dNdEta_eta0, Centrality_err, Centrality_err, dNdEta_eta0_err, dNdEta_eta0_err)
0210     gsphenix_data_divnpart2_auau_200gev = TGraphAsymmErrors(len(centnpart), centnpart, dNdEta_eta0_divnpart2, centnparterr, centnparterr, dNdEta_eta0_divnpart2_err, dNdEta_eta0_divnpart2_err)
0211
0212     return gsphenix_data_auau_200gev, gsphenix_data_divnpart2_auau_200gev
0213
0214
0215 def msrmnt_dict():
0216     msrmntdict = {}
0217
0218     msrmntdict['cms_pbpb_2p76'] = cms_pbpb_2p76()
0219     msrmntdict['cms_xexe_5p44'] = cms_xexe_5p44()
0220     msrmntdict['atlas_pbpb_2p76'] = atlas_pbpb_2p76()
0221     msrmntdict['alice_pbpb_2p76'] = alice_pbpb_2p76()
0222     msrmntdict['alice_pbpb_5p02'] = alice_pbpb_5p02()
0223     msrmntdict['alice_xexe_5p44'] = alice_xexe_5p44()
0224     msrmntdict['phobos_auau_0p2'] = phobos_auau_0p2()
0225     msrmntdict['phobos_auau_0p13'] = phobos_auau_0p13()
0226     msrmntdict['phobos_auau_0p0624'] = phobos_auau_0p0624()
0227     msrmntdict['phobos_auau_0p0196'] = phobos_auau_0p0196()
0228     msrmntdict['phenix_auau_0p2'] = phenix_auau_0p2()
0229     msrmntdict['phenix_auau_0p13'] = phenix_auau_0p13()
0230     msrmntdict['phenix_auau_0p0624'] = phenix_auau_0p0624()
0231     msrmntdict['phenix_auau_0p039'] = phenix_auau_0p039()
0232     msrmntdict['phenix_auau_0p027'] = phenix_auau_0p027()
0233     msrmntdict['phenix_auau_0p0196'] = phenix_auau_0p0196()
0234     msrmntdict['phenix_auau_0p0145'] = phenix_auau_0p0145()
0235     msrmntdict['phenix_auau_0p0077'] = phenix_auau_0p0077()
0236     msrmntdict['brahms_auau_0p2'] = brahms_auau_0p2()
0237
0238     # sphenix_rawhijing_auau_200gev, sphenix_rawhijing_divnpart2_auau_200gev = msrmnt_sphenix_generator_auau_200gev('HIJING_closure_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0')
0239     # sphenix_rawepos_auau_200gev, sphenix_rawepos_divnpart2_auau_200gev = msrmnt_sphenix_generator_auau_200gev('EPOS_closure_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0')
0240     # sphenix_rawampt_auau_200gev, sphenix_rawampt_divnpart2_auau_200gev = msrmnt_sphenix_generator_auau_200gev('AMPT_closure_dRcut0p5_NominalVtxZ_RandomClusSet0_clusAdcCutSet0_clusPhiSizeCutSet0')
0241     sphenix_rawhijing_divnpart2_auau_200gev = msrmnt_generator_dNdetaNpart2_auau_200gev('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/generator_dNdEtaNpart2/truthdNdEta_HIJING.root')
0242     sphenix_rawepos_divnpart2_auau_200gev = msrmnt_generator_dNdetaNpart2_auau_200gev('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/generator_dNdEtaNpart2/truthdNdEta_EPOS.root')
0243     sphenix_rawampt_divnpart2_auau_200gev = msrmnt_generator_dNdetaNpart2_auau_200gev('/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/generator_dNdEtaNpart2/truthdNdEta_AMPT.root')
0244     # sphnx_sim_auau_200gev, sphnx_sim_divnpart2_auau_200gev = msrmnt_sphenix_sim_auau_200gev()
0245     sphnx_data_auau_200gev_cmsapp, sphnx_data_divnpart2_auau_200gev_cmsapp = msrmnt_sphenix_data_auau_200gev('cms')
0246     sphnx_data_auau_200gev_phobosapp, sphnx_data_divnpart2_auau_200gev_phobosapp = msrmnt_sphenix_data_auau_200gev('phobos')
0247     sphnx_data_auau_200gev_combined, sphnx_data_divnpart2_auau_200gev_combined = msrmnt_sphenix_data_auau_200gev('combined')
0248
0249     # msrmntdict['sphenix_rawhijing_auau_200gev'] = sphenix_rawhijing_auau_200gev
0250     # msrmntdict['sphenix_rawepos_auau_200gev'] = sphenix_rawepos_auau_200gev
0251     # msrmntdict['sphenix_rawampt_auau_200gev'] = sphenix_rawampt_auau_200gev
0252     # msrmntdict['sphenix_sim_auau_200gev'] = sphnx_sim_auau_200gev
0253     msrmntdict['sphenix_data_auau_200gev_cmsapp'] = sphnx_data_auau_200gev_cmsapp
0254     msrmntdict['sphenix_data_auau_200gev_cmsapp_redraw'] = sphnx_data_auau_200gev_cmsapp
0255     msrmntdict['sphenix_data_auau_200gev_phobosapp'] = sphnx_data_auau_200gev_phobosapp
0256     msrmntdict['sphenix_data_auau_200gev_phobosapp_redraw'] = sphnx_data_auau_200gev_phobosapp
0257     msrmntdict['sphenix_data_auau_200gev_combined'] = sphnx_data_auau_200gev_combined
0258     msrmntdict['sphenix_data_auau_200gev_combined_redraw'] = sphnx_data_auau_200gev_combined
0259
0260     # dndeta / (<npart> / 2)
0261     msrmntdict['phenix_auau_0p2_divnpart2'] = phenix_auau_0p2_divnpart2() # PHENIX only provide dndeta / (<npart>)
0262     msrmntdict['brahms_auau_0p2_divnpart2'] = brahms_auau_0p2_divnpart2()
0263     msrmntdict['phobos_auau_0p2_divnpart2'] = phobos_auau_0p2_divnpart2()
0264     msrmntdict['sphenix_rawhijing_auau_0p2_divnpart2'] = sphenix_rawhijing_divnpart2_auau_200gev
0265     msrmntdict['sphenix_rawepos_auau_0p2_divnpart2'] = sphenix_rawepos_divnpart2_auau_200gev
0266     msrmntdict['sphenix_rawampt_auau_0p2_divnpart2'] = sphenix_rawampt_divnpart2_auau_200gev
0267     # msrmntdict['sphenix_sim_auau_0p2_divnpart2'] = sphnx_sim_divnpart2_auau_200gev
0268     msrmntdict['sphenix_data_auau_0p2_divnpart2_cmsapp'] = sphnx_data_divnpart2_auau_200gev_cmsapp
0269     msrmntdict['sphenix_data_auau_0p2_divnpart2_cmsapp_redraw'] = sphnx_data_divnpart2_auau_200gev_cmsapp
0270     msrmntdict['sphenix_data_auau_0p2_divnpart2_phobosapp'] = sphnx_data_divnpart2_auau_200gev_phobosapp
0271     msrmntdict['sphenix_data_auau_0p2_divnpart2_phobosapp_redraw'] = sphnx_data_divnpart2_auau_200gev_phobosapp
0272     msrmntdict['sphenix_data_auau_0p2_divnpart2_combined'] = sphnx_data_divnpart2_auau_200gev_combined
0273     msrmntdict['sphenix_data_auau_0p2_divnpart2_combined_redraw'] = sphnx_data_divnpart2_auau_200gev_combined
0274
0275     msrmntdict['rhic_combine_dndetadivnpart2'] = combine_RHIC()
0276     msrmntdict['rhic_combine_dndetadivnpart2_redraw'] = combine_RHIC()
0277
0278     msrmntdict['empty'] = TGraphErrors()
0279     return msrmntdict
0280
0281 def combine_RHIC():
0282     # combine PHENIX, BRAHMS, PHOBOS, and sPHENIX for dNdeta / (<npart> / 2)
0283     # print out the dNdeta for each experiment, then the npart/2 for each experiment
0284     # print out the combined dNdeta and npart/2
0285     # start with PHENIX
0286
0287     centralitybin_phenix = [[0, 5], [5, 10], [10 ,15], [15, 20], [20, 25], [25, 30], [30, 35], [35, 40], [40, 45], [45, 50], [50, 55], [55, 60], [60, 65], [65, 70]]
0288     centralitybin_phobos = [[0, 3], [3, 6], [6,10], [10, 15],[15, 20], [20, 25], [25, 30], [30, 35], [35, 40], [40, 45], [45, 50], [50, 55], [55, 60], [60, 65], [65, 70]]
0289     centralitybin_brahms = [[0, 5], [5, 10], [10, 20], [20, 30], [30, 40], [40, 50]]
0290     centralitybin_sphenix = centralitybin_phobos
0291     sphenix_npart, sphenix_nparterr = sphenix_centralitynpart()
0292     sphnx_data_auau_200gev_combined, sphnx_data_divnpart2_auau_200gev_combined = msrmnt_sphenix_data_auau_200gev('combined')
0293
0294     # centralitybin_combine = [[0,3], [3,6], [6,10], [10,15], [15,20], [20,25], [25,30], [30,35], [35,40], [40,45], [45,50], [50,55], [55,60], [60,65], [65,70]]
0295     centralitybin_combine = [[0,3], [0,5], [3,6], [5,10], [6,10], [10,15], [15,20], [20,25], [25,30], [30,35], [35,40], [40,45], [45,50], [50,55], [55,60], [60,65], [65,70]]
0296
0297     gr_phenix_auau_0p2 = phenix_auau_0p2()
0298     phenix_dndeta, phenix_dndeta_err = [], []
0299     for i in range(gr_phenix_auau_0p2.GetN()):
0300         phenix_dndeta.append(gr_phenix_auau_0p2.GetY()[i])
0301         phenix_dndeta_err.append(gr_phenix_auau_0p2.GetErrorY(i))
0302
0303     # add 60-65% and 65-70% points for phenix that are not in the graph and HepData
0304     phenix_dndeta.append(37.5)
0305     phenix_dndeta.append(25.6)
0306     phenix_dndeta_err.append(5.4)
0307     phenix_dndeta_err.append(4.5)
0308
0309     gr_brahms_auau_0p2_divnpart2 = brahms_auau_0p2_divnpart2()
0310
0311     # print('------------------------------------')
0312     # for i in range(len(phenix_dndeta)):
0313     #     print ('PHENIX: Centrality = {}-{}%, dNdeta/(0.5npart) = {:.3f}+-{:.3f}, npart = {:.3f}+-{:.3f}'.format(centralitybin_phenix[i][0], centralitybin_phenix[i][1], phenix_auau_0p2_dndetadivnpart2[-(i+1)], phenix_auau_0p2_dndetadivnpart2_err[-(i+1)], phenix_auau_0p2_npart[-(i+1)],phenix_auau_0p2_nparterr[-(i+1)]))
0314
0315     # print('------------------------------------')
0316     # for i in range(len(phobos_auau_0p2_raw)):
0317     #     print ('PHOBOS: Centrality = {}-{}%, dNdeta/(0.5npart) = {:.3f}+-{:.3f}, npart = {:.3f}+-{:.3f}'.format(centralitybin_phobos[i][0], centralitybin_phobos[i][1], phobos_auau_0p2_dndetadivnpart2[-(i+1)], phobos_auau_0p2_dndetadivnpart2_err[-(i+1)], phobos_auau_0p2_npart[-(i+1)],phobos_auau_0p2_nparterr[-(i+1)]))
0318
0319     # print('------------------------------------')
0320     # for i in range(len(brahms_auau_0p2_raw)):
0321     #     print ('BRAHMS: Centrality = {}-{}%, dNdeta/(0.5npart) = {:.3f}+-{:.3f}, npart = {:.3f}+-{:.3f}'.format(centralitybin_brahms[i][0], centralitybin_brahms[i][1], gr_brahms_auau_0p2_divnpart2.GetY()[len(gr_brahms_auau_0p2_divnpart2.GetY())-i-1], gr_brahms_auau_0p2_divnpart2.GetErrorY(len(gr_brahms_auau_0p2_divnpart2.GetY())-i-1), brahms_auau_0p2_npart[-(i+1)], brahms_auau_0p2_nparterr[-(i+1)]))
0322
0323     # print('------------------------------------')
0324     # for i in range(len(sphnx_data_auau_200gev_combined.GetY())):
0325     #     print ('sPHENIX: Centrality = {}-{}%, dNdeta/(0.5npart) = {:.3f}+-{:.3f}, npart = {:.3f}+-{:.3f}'.format(centralitybin_sphenix[i][0], centralitybin_sphenix[i][1], sphnx_data_divnpart2_auau_200gev_combined.GetY()[i], sphnx_data_divnpart2_auau_200gev_combined.GetErrorY(i), sphenix_npart[i], sphenix_nparterr[i]))
0326
0327     # print('------------------------------------')
0328
0329     rhic_combine_dndetadivnpart2 = []
0330     rhic_combine_dndetadivnpart2_err = []
0331     rhic_combine_npart = []
0332     rhic_combine_nparterr = []
0333
0334     # combine the data points
0335     for i in range(len(centralitybin_combine)):
0336         # print('Centrality = {}-{} %'.format(centralitybin_combine[i][0], centralitybin_combine[i][1]))
0337         # check if centralitybin_combine[i] is in centralitybin_phenix, centralitybin_phobos, centralitybin_brahms, and centralitybin_sphenix
0338         # if it is, then combine the data points using the function weighted_average(values, uncertainties)
0339
0340         values_dndeta = []
0341         uncertainties_dndeta = []
0342         values_npart = []
0343         uncertainties_npart = []
0344
0345         if centralitybin_combine[i] in centralitybin_phenix:
0346             phenix_index = centralitybin_phenix.index(centralitybin_combine[i])
0347             values_dndeta.append(phenix_auau_0p2_dndetadivnpart2[-(phenix_index+1)])
0348             uncertainties_dndeta.append(phenix_auau_0p2_dndetadivnpart2_err[-(phenix_index+1)])
0349             values_npart.append(phenix_auau_0p2_npart[-(phenix_index+1)])
0350             uncertainties_npart.append(phenix_auau_0p2_nparterr[-(phenix_index+1)])
0351             # print ('PHENIX: dNdeta/(0.5npart) = {:.3f}+-{:.3f}, npart = {:.3f}+-{:.3f}'.format(phenix_auau_0p2_dndetadivnpart2[-(phenix_index+1)], phenix_auau_0p2_dndetadivnpart2_err[-(phenix_index+1)], phenix_auau_0p2_npart[-(phenix_index+1)], phenix_auau_0p2_nparterr[-(phenix_index+1)]))
0352
0353         if centralitybin_combine[i] in centralitybin_phobos:
0354             phobos_index = centralitybin_phobos.index(centralitybin_combine[i])
0355             values_dndeta.append(phobos_auau_0p2_dndetadivnpart2[-(phobos_index+1)])
0356             uncertainties_dndeta.append(phobos_auau_0p2_dndetadivnpart2_err[-(phobos_index+1)])
0357             values_npart.append(phobos_auau_0p2_npart[-(phobos_index+1)])
0358             uncertainties_npart.append(phobos_auau_0p2_nparterr[-(phobos_index+1)])
0359             # print ('PHOBOS: dNdeta/(0.5npart) = {:.3f}+-{:.3f}, npart = {:.3f}+-{:.3f}'.format(phobos_auau_0p2_dndetadivnpart2[-(phobos_index+1)], phobos_auau_0p2_dndetadivnpart2_err[-(phobos_index+1)], phobos_auau_0p2_npart[-(phobos_index+1)], phobos_auau_0p2_nparterr[-(phobos_index+1)]))
0360
0361         if centralitybin_combine[i] in centralitybin_brahms:
0362             brahms_index = centralitybin_brahms.index(centralitybin_combine[i])
0363             values_dndeta.append(gr_brahms_auau_0p2_divnpart2.GetY()[len(gr_brahms_auau_0p2_divnpart2.GetY())-i-1])
0364             uncertainties_dndeta.append(gr_brahms_auau_0p2_divnpart2.GetErrorY(len(gr_brahms_auau_0p2_divnpart2.GetY())-i-1))
0365             values_npart.append(brahms_auau_0p2_npart[-(brahms_index+1)])
0366             uncertainties_npart.append(brahms_auau_0p2_nparterr[-(brahms_index+1)])
0367             # print ('BRAHMS: dNdeta/(0.5npart) = {:.3f}+-{:.3f}, npart = {:.3f}+-{:.3f}'.format(gr_brahms_auau_0p2_divnpart2.GetY()[len(gr_brahms_auau_0p2_divnpart2.GetY())-i-1], gr_brahms_auau_0p2_divnpart2.GetErrorY(len(gr_brahms_auau_0p2_divnpart2.GetY())-i-1), brahms_auau_0p2_npart[-(brahms_index+1)], brahms_auau_0p2_nparterr[-(brahms_index+1)]))
0368
0369         if centralitybin_combine[i] in centralitybin_sphenix:
0370             sphenix_index = centralitybin_sphenix.index(centralitybin_combine[i])
0371             values_dndeta.append(sphnx_data_divnpart2_auau_200gev_combined.GetY()[sphenix_index])
0372             uncertainties_dndeta.append(sphnx_data_divnpart2_auau_200gev_combined.GetErrorY(sphenix_index))
0373             values_npart.append(sphenix_npart[sphenix_index])
0374             uncertainties_npart.append(sphenix_nparterr[sphenix_index])
0375             # print ('sPHENIX: dNdeta/(0.5npart) = {:.3f}+-{:.3f}, npart = {:.3f}+-{:.3f}'.format(sphnx_data_divnpart2_auau_200gev_combined.GetY()[sphenix_index], sphnx_data_divnpart2_auau_200gev_combined.GetErrorY(sphenix_index), sphenix_npart[sphenix_index], sphenix_nparterr[sphenix_index]))
0376
0377         dndeta_combined, dndeta_combined_err = weighted_average(values_dndeta, uncertainties_dndeta)
0378         npart_combined, npart_combined_err = weighted_average(values_npart, uncertainties_npart)
0379         # print ('Combined: dNdeta/(0.5npart) = ${:.2f}\pm{:.2f}~({:.2f}\%)$, npart = {:.3f}+-{:.3f}'.format(dndeta_combined, dndeta_combined_err, dndeta_combined_err/dndeta_combined*100, npart_combined, npart_combined_err))
0380         # print('------------------------------------')
0381
0382         rhic_combine_dndetadivnpart2.append(dndeta_combined)
0383         rhic_combine_dndetadivnpart2_err.append(dndeta_combined_err)
0384         rhic_combine_npart.append(npart_combined)
0385         rhic_combine_nparterr.append(npart_combined_err)
0386
0387     # make a TGraphAsymmErrors for the combined data points
0388     gr_rhic_combined = TGraphAsymmErrors(len(centralitybin_combine), array('f', rhic_combine_npart), array('f', rhic_combine_dndetadivnpart2), array('f', rhic_combine_nparterr), array('f', rhic_combine_nparterr), array('f', rhic_combine_dndetadivnpart2_err), array('f', rhic_combine_dndetadivnpart2_err))
0389
0390     return gr_rhic_combined
0391
0392
0393 def gstyle(g, mstyle, msize, color, lstyle, lwidth, alpha):
0394     g.SetMarkerStyle(mstyle)
0395     g.SetMarkerSize(msize)
0396     g.SetMarkerColor(color)
0397     g.SetLineColor(color)
0398     g.SetLineStyle(lstyle)
0399     g.SetLineWidth(lwidth)
0400     g.SetFillColorAlpha(color, alpha)
0401
0402
0403 def dNdeta_eta0_Summary(docombine, jsonpath, canvname, axisrange = [-1, 71, 20, 1100], canvsize = [800, 600], grlegpos = [0.35, 0.18, 0.9, 0.45], prelimtext=''):
0404     xmin = axisrange[0]
0405     xmax = axisrange[1]
0406     ymin = axisrange[2]
0407     ymax = axisrange[3]
0408     # Create a canvas
0409     c = TCanvas('c_'+canvname, 'c_'+canvname, canvsize[0], canvsize[1])
0410     c.cd()
0411     # Set left margin and logarithmic y-axis
0412     gPad.SetTopMargin(TopMargin)
0413     gPad.SetLeftMargin(LeftMargin)
0414     gPad.SetLogy()
0415     # Create a frame for the graph
0416     gframe = TH1F('gframe_' + canvname, 'gframe_' + canvname, 1, xmin, xmax)
0417     gframe.SetLabelOffset(999, 'X')
0418     gframe.SetTickLength(0, 'X')
0419     gframe.GetXaxis().SetTitle('Centrality [%]')
0420     gframe.GetYaxis().SetTitle('#LT#frac{dN_{ch}}{d#eta}#GT')
0421     gframe.GetYaxis().SetMoreLogLabels()
0422     gframe.GetYaxis().SetTitleOffset(2.1)
0423     gframe.SetAxisRange(ymin, ymax, 'Y')
0424     gframe.Draw()
0425     # Create and draw the x-axis on top
0426     axis = TGaxis(xmax, ymin, xmin, ymin, xmin, xmax, 511, '-')
0427     axis.SetLabelOffset(-0.032)
0428     axis.SetLabelFont(43)
0429     # print('gframe.GetYaxis().GetLabelSize()={}'.format(gframe.GetYaxis().GetLabelSize()))
0430     axis.SetLabelSize(30)
0431     axis.Draw()
0432     # Create and draw the y-axis on top
0433     axisup = TGaxis(xmax, ymax, xmin, ymax, xmin, xmax, 511, '+')
0434     axisup.SetLabelOffset(1)
0435     axisup.Draw()
0436     # Create the legend for the TGraphs
0437     # gr_leg = TLegend(grlegpos[0], grlegpos[1], grlegpos[2], grlegpos[3])
0438     # gr_leg = TLegend(gPad.GetLeftMargin()+0.05, (1 - TopMargin) - 0.35, gPad.GetLeftMargin()+0.35, (1 - TopMargin) - 0.18)
0439     gr_leg = TLegend(1 - gPad.GetRightMargin() - 0.55, gPad.GetBottomMargin() + 0.05, 1 - gPad.GetRightMargin() - 0.05, gPad.GetBottomMargin() + 0.2)
0440     gr_leg.SetNColumns(2)
0441     gr_leg.SetTextSize(0.033)
0442     gr_leg.SetFillStyle(0)
0443     # Load measurement dictionaries
0444     dict_msrmnt= msrmnt_dict()
0445     # Load the json file
0446     with open(jsonpath, 'r') as fp:
0447         grtags = json.load(fp)
0448         for tagName, par in grtags.items():
0449             if not include_rawhijing and 'rawhijing' in tagName:
0450                 continue
0451
0452             if docombine:
0453                 if 'phobosapp' in tagName or 'cmsapp' in tagName:
0454                     continue
0455             else:
0456                 if 'combined' in tagName:
0457                     continue
0458
0459             gr = dict_msrmnt[tagName]
0460             gstyle(gr, par['markerstyle'], par['markersize'], TColor.GetColor(par['markercolor']), par['linestyle'], par['linewidth'], par['alpha'])
0461             gr.Draw(par['DrawOption'][0])
0462             if par['DrawOption'][1] != '':
0463                 gr.Draw(par['DrawOption'][1])
0464             # add legend if par['Legendtext'] does not contain 'redraw'
0465             if 'redraw' not in tagName:
0466                 gr_leg.AddEntry(gr, par['Legendtext'], par['Legendstyle'])
0467
0468     gr_leg.Draw()
0469     sphnxleg = TLegend(gPad.GetLeftMargin()+0.05, (1 - TopMargin) - 0.18, gPad.GetLeftMargin()+0.1, (1 - TopMargin) - 0.06)
0470     sphnxleg.SetTextAlign(kHAlignLeft + kVAlignTop)
0471     sphnxleg.SetTextSize(0.045)
0472     sphnxleg.SetFillStyle(0)
0473     sphnxleg.AddEntry('', '#it{#bf{sPHENIX}} ' + prelimtext, '')
0474     sphnxleg.AddEntry('', 'Au+Au #sqrt{s_{NN}}=200 GeV', '')
0475     sphnxleg.Draw()
0476     if (docombine):
0477         c.SaveAs(plotpath + 'dNdEta_eta0_{}_combine.png'.format(canvname))
0478         c.SaveAs(plotpath + 'dNdEta_eta0_{}_combine.pdf'.format(canvname))
0479     else:
0480         c.SaveAs(plotpath + 'dNdEta_eta0_{}.png'.format(canvname));
0481         c.SaveAs(plotpath + 'dNdEta_eta0_{}.pdf'.format(canvname));
0482
0483
0484 def dNdeta_eta0_divnpart2_Summary(docombine, rhiccombine, plotstyle, jsonpath, canvname, axisrange = [0, 400, 1.49, 4.71], canvsize = [800, 600], grlegpos = [0.35, 0.18, 0.9, 0.45], prelimtext=''):
0485     # plotstyle:
0486     # 1 - all measurements and theories comparison
0487     # 2 - only measurement comparison
0488     # 3 - only sPHENIX measurement and theories
0489
0490     xmin = axisrange[0]
0491     xmax = axisrange[1]
0492     ymin = axisrange[2]
0493     ymax = axisrange[3]
0494     # Create a canvas
0495     c = TCanvas('c_'+canvname, 'c_'+canvname, canvsize[0], canvsize[1])
0496     c.cd()
0497     # Set left margin and logarithmic y-axis
0498     gPad.SetTopMargin(TopMargin)
0499     gPad.SetLeftMargin(LeftMargin - 0.05)
0500     gPad.SetLogy(0)
0501     # Create a frame for the graph
0502     gframe = TH1F('gframe_' + canvname, 'gframe_' + canvname, 1, axisrange[0], axisrange[1])
0503     # gframe.GetXaxis().SetTitle('#LTN_{part}#GT')
0504     gframe.GetXaxis().SetTitle('N_{part}')
0505     gframe.GetYaxis().SetTitle('#LT#frac{dN_{ch}}{d#eta}#GT/(#LTN_{part}#GT/2)')
0506     gframe.GetYaxis().SetTitleOffset(1.55)
0507     gframe.GetXaxis().SetRangeUser(axisrange[0], axisrange[1])
0508     gframe.SetAxisRange(axisrange[2], axisrange[3], 'Y')
0509     gframe.Draw()
0510     # Create the legend for the TGraphs
0511     # gr_leg = TLegend(grlegpos[0], grlegpos[1], grlegpos[2], grlegpos[3])
0512     # gr_leg = TLegend(gPad.GetLeftMargin(), 1 - gPad.GetTopMargin() + 0.02, gPad.GetLeftMargin() + 0.72, 0.97)
0513
0514     grleg_x1 = 1 - gPad.GetRightMargin() - 0.43
0515     grleg_y1 = gPad.GetBottomMargin() + 0.05 # default for non-RHIC combine
0516     grleg_x2 = 1 - gPad.GetRightMargin() - 0.04 # default for non-RHIC combine
0517     grleg_y2 = gPad.GetBottomMargin() + 0.24 # default for non-RHIC combine
0518     if rhiccombine:
0519         grleg_x1 = 1 - gPad.GetRightMargin() - 0.55
0520         grleg_y1 = gPad.GetBottomMargin() + 0.08
0521         grleg_x2 = 1 - gPad.GetRightMargin() - 0.12
0522         grleg_y2 = gPad.GetBottomMargin() + 0.24
0523     else:
0524         if plotstyle == 1:
0525             grleg_x1 = 1 - gPad.GetRightMargin() - 0.43
0526             grleg_y1 = gPad.GetBottomMargin() + 0.05
0527             grleg_x2 = 1 - gPad.GetRightMargin() - 0.04
0528             grleg_y2 = gPad.GetBottomMargin() + 0.24
0529         elif plotstyle == 2: # only measurement comparison
0530             grleg_x1 = 1 - gPad.GetRightMargin() - 0.31
0531             grleg_y1 = gPad.GetBottomMargin() + 0.05
0532             grleg_x2 = 1 - gPad.GetRightMargin() - 0.1
0533             grleg_y2 = gPad.GetBottomMargin() + 0.24
0534         elif plotstyle == 3:
0535             grleg_x1 = 1 - gPad.GetRightMargin() - 0.31
0536             grleg_y1 = gPad.GetBottomMargin() + 0.05
0537             grleg_x2 = 1 - gPad.GetRightMargin() - 0.1
0538             grleg_y2 = gPad.GetBottomMargin() + 0.24
0539         else:
0540             print ('Error: Invalid plotstyle provided. Must be 1, 2, or 3. Using default grleg_x1.')
0541             grleg_x1 = 1 - gPad.GetRightMargin() - 0.43
0542             grleg_y1 = gPad.GetBottomMargin() + 0.05
0543             grleg_x2 = 1 - gPad.GetRightMargin() - 0.04
0544             grleg_y2 = gPad.GetBottomMargin() + 0.24
0545
0546
0547     gr_leg = TLegend(grleg_x1, grleg_y1, grleg_x2, grleg_y2)
0548     gr_leg.SetNColumns(2) if plotstyle == 1 else gr_leg.SetNColumns(1)
0549     gr_leg.SetTextSize(0.033)
0550     gr_leg.SetFillStyle(0)
0551
0552     leg_generator = TLegend(1 - gPad.GetRightMargin() - 0.545,
0553                             gPad.GetBottomMargin() + 0.03,
0554                             1 - gPad.GetRightMargin() - 0.02,
0555                             gPad.GetBottomMargin() + 0.08)
0556     leg_generator.SetNColumns(3)
0557     leg_generator.SetTextSize(0.033)
0558     leg_generator.SetFillStyle(0)
0559
0560     print ('-----------------------------------------------')
0561     # Load measurement dictionaries
0562     dict_msrmnt= msrmnt_dict()
0563     # Load the json file
0564     with open(jsonpath, 'r') as fp:
0565         grtags = json.load(fp)
0566         for tagName, par in grtags.items():
0567             if docombine:
0568                 if 'phobosapp' in tagName or 'cmsapp' in tagName:
0569                     continue
0570             else:
0571                 if 'combined' in tagName:
0572                     continue
0573
0574             if plotstyle == 2: # only measurement comparison
0575                 if 'rawhijing' in tagName or 'rawepos' in tagName or 'rawampt' in tagName or 'empty' in tagName:
0576                     continue
0577
0578             if plotstyle == 3: # only sPHENIX measurement and theories
0579                 print (tagName)
0580                 if tagName.startswith('phenix_') or tagName.startswith('brahms_') or tagName.startswith('phobos_') or 'empty' in tagName:
0581                     continue
0582
0583             # print (tagName)
0584
0585             gr = dict_msrmnt[tagName]
0586             gstyle(gr, par['markerstyle'], par['markersize'], TColor.GetColor(par['markercolor']), par['linestyle'], par['linewidth'], par['alpha'])
0587             gr.Draw(par['DrawOption'][0])
0588             if par['DrawOption'][1] != '':
0589                 gr.Draw(par['DrawOption'][1])
0590
0591             # add legend if par['Legendtext'] does not contain 'redraw'
0592             if not rhiccombine:
0593                 if ('redraw' not in tagName):
0594                     gr_leg.AddEntry(gr, par['Legendtext'], par['Legendstyle'])
0595             else:
0596                 if ('redraw' not in tagName and 'rawhijing' not in tagName and 'rawepos' not in tagName and 'rawampt' not in tagName):
0597                     gr_leg.AddEntry(gr, par['Legendtext'], par['Legendstyle'])
0598
0599                 if ('rawhijing' in tagName or 'rawepos' in tagName or 'rawampt' in tagName):
0600                     leg_generator.AddEntry(gr, par['Legendtext'], par['Legendstyle'])
0601
0602     gr_leg.Draw()
0603     if rhiccombine:
0604         leg_generator.Draw()
0605
0606     sphnxleg = TLegend(gPad.GetLeftMargin()+0.05, (1 - TopMargin) - 0.18, gPad.GetLeftMargin()+0.1, (1 - TopMargin) - 0.06)
0607     sphnxleg.SetTextAlign(kHAlignLeft + kVAlignTop)
0608     sphnxleg.SetTextSize(0.045)
0609     sphnxleg.SetFillStyle(0)
0610     sphnxleg.AddEntry('', '#it{#bf{sPHENIX}} ' + prelimtext, '')
0611     sphnxleg.AddEntry('', 'Au+Au #sqrt{s_{NN}}=200 GeV', '')
0612     sphnxleg.Draw()
0613     if (docombine):
0614         if not rhiccombine:
0615             c.SaveAs(plotpath + 'dNdEta_eta0_divnpart2_{}_combine_style{}.png'.format(canvname, plotstyle))
0616             c.SaveAs(plotpath + 'dNdEta_eta0_divnpart2_{}_combine_style{}.pdf'.format(canvname, plotstyle))
0617         else:
0618             c.SaveAs(plotpath + 'dNdEta_eta0_divnpart2_{}_combine_rhiccombine_style{}.png'.format(canvname, plotstyle))
0619             c.SaveAs(plotpath + 'dNdEta_eta0_divnpart2_{}_combine_rhiccombine_style{}.pdf'.format(canvname, plotstyle))
0620     else:
0621         if not rhiccombine:
0622             c.SaveAs(plotpath + 'dNdEta_eta0_divnpart2_{}_style{}.png'.format(canvname, plotstyle))
0623             c.SaveAs(plotpath + 'dNdEta_eta0_divnpart2_{}_style{}.pdf'.format(canvname, plotstyle))
0624         else:
0625             c.SaveAs(plotpath + 'dNdEta_eta0_divnpart2_{}_rhiccombine_style{}.png'.format(canvname, plotstyle))
0626             c.SaveAs(plotpath + 'dNdEta_eta0_divnpart2_{}_rhiccombine_style{}.pdf'.format(canvname, plotstyle))
0627
0628
0629 if __name__ == '__main__':
0630     parser = OptionParser(usage='usage: %prog ver [options -h]')
0631     parser.add_option('--combine', action='store_true', dest='combine', default=False, help='Combine two analyses')
0632
0633     (opt, args) = parser.parse_args()
0634     print('opt: {}'.format(opt))
0635
0636     combine = opt.combine
0637
0638     plotpath = './dNdEtaFinal/'
0639     os.makedirs(plotpath, exist_ok=True)
0640
0641     dNdeta_eta0_Summary(combine, './measurements/RHIC_plotparam.json', 'RHIC')
0642     dNdeta_eta0_divnpart2_Summary(combine, False, 1, './measurements/RHIC_dndetadivnpart2_plotparam.json', 'RHIC')
0643     dNdeta_eta0_divnpart2_Summary(combine, False, 2, './measurements/RHIC_dndetadivnpart2_plotparam.json', 'RHIC')
0644     dNdeta_eta0_divnpart2_Summary(combine, False, 3, './measurements/RHIC_dndetadivnpart2_plotparam.json', 'RHIC')
0645     dNdeta_eta0_divnpart2_Summary(combine, True, 1, './measurements/RHIC_combine_dndetadivnpart2_plotparam.json', 'RHIC')
0646
0647     combine_RHIC()
0648