#!/usr/bin/python import os, sys, math import ROOT from ROOT import TH1D, TH2D, TFile, TCanvas, TLegend from ROOT import TGraphErrors, TLine, TGaxis def systematics_hist(lg1, output, legend, max, plotpol = 0): # generate a horizontal histogram with alpha numeric bins hrho = TH1D('hrho', '; ; Systematic Uncertainties (#times 10^{4})', len(lg1), 0, len(lg1)) hdelta = TH1D('hdelta', '; ; Systematic Uncertainties (#times 10^{4})', len(lg1), 0, len(lg1)) hxi = TH1D('hxi', '; ; Systematic Uncertainties (#times 10^{4})', len(lg1), 0, len(lg1)) # assign each of the bins to a particular value i = len(lg1) for a in lg1: hrho.SetBinContent(i,a[1]) hrho.GetXaxis().SetBinLabel(i,a[0] + '\t \t \t %(i)0.3f \t %(j)0.3f \t %(k)0.3f \t' %{'i':a[1],'j':a[2],'k':a[3]}) hdelta.SetBinContent(i,a[2]) hdelta.GetXaxis().SetBinLabel(i,a[0]) hxi.SetBinContent(i,a[3]) hxi.GetXaxis().SetBinLabel(i,a[0]) i-=1 # h = [hrho, hdelta, hxi] vheight = (len(lg1) + 1)*50 if plotpol : vheight = 100 c1 = TCanvas('c1','systematics', 800, vheight) c1.SetFillColor(10) c1.SetLeftMargin(0.38) c1.SetTopMargin(0.01) c1.SetRightMargin(0.05) #c1.SetBottomMargin(0.15) c1.SetGridx() c1.SetGridy() leg = TLegend(0.75, 0.20, 0.94, 0.40) leg.AddEntry(hrho,legend[0],'f') leg.AddEntry(hdelta,legend[1],'f') leg.AddEntry(hxi,legend[2],'f') leg.SetFillColor(10) if not plotpol: c1.SetBottomMargin(0.15) hrho.SetBarWidth(0.25) hrho.SetBarOffset(0.70) hrho.SetFillColor(632) hrho.SetLineColor(1) hrho.SetStats(0) hrho.GetXaxis().SetLabelSize(0.07) hrho.GetYaxis().SetLabelSize(0.045) hrho.GetYaxis().SetTitleSize(0.045) hdelta.SetBarWidth(0.25) hdelta.SetBarOffset(0.40) hdelta.SetFillColor(600) hdelta.SetLineColor(1) hxi.SetBarWidth(0.25) hxi.SetBarOffset(0.10) hxi.SetFillColor(418) hxi.SetLineColor(1) hrho.SetMaximum(max) hrho.GetXaxis().SetLabelOffset(0.05) hrho.Draw('hbar2') hdelta.Draw('hbar2same') hxi.Draw('hbar2same') leg.Draw('same') else: hxi.SetBarWidth(0.5) hxi.SetFillColor(418) hxi.SetLineColor(1) c1.SetBottomMargin(0.45) c1.SetTopMargin(0.01) hxi.GetXaxis().SetLabelSize(0.3) hxi.GetYaxis().SetLabelSize(0.2) hxi.GetYaxis().SetTitleSize(0.2) hxi.SetStats(0) hxi.Draw('hbar2') c1.Print(output) def systematics_hist2(lg1, lg2, output, legend, max, plotpol = 0): if len(lg1) != len(lg2): print 'Systematics input has the wrong number of elements' # generate a horizontal histogram with alpha numeric bins hrho1 = TH1D('hrho1', '; ; Systematic Uncertainties (#times 10^{4})', len(lg1), 0, len(lg1)) hdelta1 = TH1D('hdelta1', '; ; Systematic Uncertainties (#times 10^{4})', len(lg1), 0, len(lg1)) hxi1 = TH1D('hxi1', '; ; Systematic Uncertainties (#times 10^{4})', len(lg1), 0, len(lg1)) hrho2 = TH1D('hrho2', '; ; Systematic Uncertainties (#times 10^{4})', len(lg2), 0, len(lg2)) hdelta2 = TH1D('hdelta2', '; ; Systematic Uncertainties (#times 10^{4})', len(lg2), 0, len(lg2)) hxi2 = TH1D('hxi2', '; ; Systematic Uncertainties (#times 10^{4})', len(lg2), 0, len(lg2)) # assign each of the bins to a particular value i = len(lg1) for a in lg1: hrho1.SetBinContent(i,a[1]) hrho1.GetXaxis().SetBinLabel(i,a[0]) hdelta1.SetBinContent(i,a[2]) hdelta1.GetXaxis().SetBinLabel(i,a[0]) hxi1.SetBinContent(i,a[3]) hxi1.GetXaxis().SetBinLabel(i,a[0]) i-=1 i = len(lg2) for a in lg2: hrho2.SetBinContent(i,a[1]) hrho2.GetXaxis().SetBinLabel(i,a[0]) hdelta2.SetBinContent(i,a[2]) hdelta2.GetXaxis().SetBinLabel(i,a[0]) hxi2.SetBinContent(i,a[3]) hxi2.GetXaxis().SetBinLabel(i,a[0]) i-=1 # h = [hrho, hdelta, hxi] vheight = (len(lg1) + 1)*50 if plotpol : vheight = 100 c1 = TCanvas('c1','systematics', 800, vheight) c1.SetFillColor(10) c1.SetLeftMargin(0.30) c1.SetTopMargin(0.01) c1.SetRightMargin(0.05) c1.SetGridx() c1.SetGridy() leg = TLegend(0.70, 0.11, 0.94, 0.45) leg.AddEntry(hrho1,legend[0],'f') leg.AddEntry(hdelta1,legend[1],'f') leg.AddEntry(hxi1,legend[2],'f') leg.AddEntry(hrho2,legend[3],'f') leg.AddEntry(hdelta2,legend[4],'f') leg.AddEntry(hxi2,legend[5],'f') leg.SetFillColor(10) if not plotpol: hrho1.SetBarWidth(0.25) hrho1.SetBarOffset(0.70) hrho1.SetFillColor(622) hrho1.SetLineColor(1) hrho1.SetStats(0) hrho1.GetXaxis().SetLabelSize(0.06) hrho1.GetYaxis().SetLabelSize(0.045) hrho1.GetYaxis().SetTitleSize(0.045) hrho2.SetBarWidth(0.25) hrho2.SetBarOffset(0.70) hrho2.SetFillColor(632) hrho2.SetLineColor(1) hrho2.SetStats(0) hrho2.GetXaxis().SetLabelSize(0.06) hdelta1.SetBarWidth(0.25) hdelta1.SetBarOffset(0.40) hdelta1.SetFillColor(590) hdelta1.SetLineColor(1) hdelta2.SetBarWidth(0.25) hdelta2.SetBarOffset(0.40) hdelta2.SetFillColor(600) hdelta2.SetLineColor(1) hxi1.SetBarWidth(0.25) hxi1.SetBarOffset(0.10) hxi1.SetFillColor(408) hxi1.SetLineColor(1) hxi2.SetBarWidth(0.25) hxi2.SetBarOffset(0.10) hxi2.SetFillColor(418) hxi2.SetLineColor(1) hrho1.SetMaximum(max) hrho1.Draw('hbar2') hdelta1.Draw('hbar2same') hxi1.Draw('hbar2same') hrho2.Draw('hbar2same') hdelta2.Draw('hbar2same') hxi2.Draw('hbar2same') leg.Draw('same') else: hxi.SetBarWidth(0.5) hxi.SetFillColor(418) hxi.SetLineColor(1) c1.SetBottomMargin(0.35) hxi.GetXaxis().SetLabelSize(0.30) hxi.GetYaxis().SetLabelSize(0.15) hxi.GetYaxis().SetTitleSize(0.15) hxi.SetStats(0) hxi.Draw('hbar2') c1.Print(output) def decay_par_graph(lpar, output, xmin, xmax, plotleg=1): ROOT.gStyle.SetOptStat(0) ROOT.gStyle.SetOptFit(111) ROOT.gStyle.SetPalette(1) ROOT.gStyle.SetCanvasColor(10) ROOT.gStyle.SetLabelFont(132) ROOT.gStyle.SetLabelFont(132,'X') ROOT.gStyle.SetLabelFont(132,'Y') ROOT.gStyle.SetLabelFont(132,'Z') ROOT.gStyle.SetTitleFont(132) ROOT.gStyle.SetTitleFont(132,'X') ROOT.gStyle.SetTitleFont(132,'Y') #ROOT.gStyle.SetTitleOffset(1.5,'X') ROOT.gStyle.SetTitleOffset(1.25,'Y') ROOT.gStyle.SetTitleFont(132,'Z') ROOT.gStyle.SetTextFont(132) ROOT.gStyle.SetStatFont(132) ROOT.gStyle.SetTitleFillColor(10) ROOT.gStyle.SetStatColor(10) ROOT.gStyle.SetPadLeftMargin(0.15) ROOT.gStyle.SetPadRightMargin(0.15) htot = TH2D('htot','',100, xmin, xmax, len(lpar), -len(lpar), 0) #hsys = TH1('hsys','',len(lpar), 0, len(lpar)) gtot = TGraphErrors(len(lpar)) gsys = TGraphErrors(len(lpar)) l1 = TLine(xmin, -len(lpar) + 1, xmax, -len(lpar) + 1) l1.SetLineStyle(2) l1.SetLineWidth(2) l1.SetLineColor(2) l2 = TLine(xmin, -1, xmax, -1) l2.SetLineStyle(2) l2.SetLineWidth(2) l2.SetLineColor(2) leg = TLegend(0.80, 0.16, 0.985, 0.6) leg.SetFillColor(10) i = 0 k = len(lpar) for a in lpar: gtot.SetPoint(i, a[1], -i-0.5) gtot.SetPointError(i, a[2], 0.3) gsys.SetPoint(i, a[1], -i-0.5) gsys.SetPointError(i, a[3], 0.2) htot.GetYaxis().SetBinLabel(k,a[0]) i += 1 k -= 1 #for i in range(0,len(lpar)): # temp = gtot.GetHistogram().GetYaxis().FindBin(-i) # gtot.GetHistogram().GetYaxis().SetBinLabel(temp, a[0]) c1 = TCanvas('c1','dk pars',800, 200) c1.SetFillColor(10) c1.SetLeftMargin(0.30) c1.SetRightMargin(0.01) c1.SetBottomMargin(0.15) c1.SetTopMargin(0.01) c1.SetGridx() htot.GetYaxis().SetLabelSize(0.20) htot.GetXaxis().SetLabelSize(0.16) htot.GetXaxis().SetNdivisions(4) htot.SetStats(0) gtot.SetLineColor(2) gtot.SetFillColor(2) gtot.SetMarkerColor(2) gtot.SetTitle('') gsys.SetLineColor(4) gsys.SetFillColor(4) gsys.SetMarkerColor(5) leg.AddEntry(gtot,'Total Error','f') leg.AddEntry(gsys,'Statistical Error','f') htot.Draw() gtot.Draw('E2same') gsys.Draw('E2same') l1.Draw('lsame') l2.Draw('lsame') if plotleg: leg.Draw('same') c1.Print(output) def main(argv): # create a series of parameter lists for the input c1 = TCanvas() lg1 = [] # lg1.append([systematic, rho, delta, xi]) legend = ['#rho^{#club}','#delta^{#club}','#xi^{#diamond}'] lg1.append(['Fringe Field Depol.', 0.0, 0.0, 34]) lg1.append(['Target Depol.', 0.0, 0.0, 14]) lg1.append(['Chamber Response', 2.9, 5.2, 10]) lg1.append(['Energy Scale',2.9, 4.1, 2.0]) lg1.append(['Positron Interactions', 1.6, 0.9, 3.0]) lg1.append(['Resolution', 1.2, 1.4, 0.5]) lg1.append(['Beam Intensity', 0.3, 0.3, 0.2]) lg1.append(['Alignment/length', 0.1, 0.2, 0.7]) lg1.append(['#eta correlations', 1.1, 0.1, 0.1]) # systematics_hist(lg1[2:], 'Systematics.eps', legend, 10.5) # systematics that only effect the asymmetry systematics_hist(lg1[:2], 'Systematics_pol.eps', legend, 35, 1) lg2 = [] lg2.append(['Chamber Response', 1.02, 1.79, 0.88]) lg2.append(['Energy Scale',1.17, 1.21, 1.46]) lg2.append(['Positron Interactions', 1.84, 1.61, 0.65]) lg2.append(['Resolution', 0.56, 0.71, 1.52]) lg2.append(['Beam Intensity', 0.16, 0.0, 0.34]) lg2.append(['Alignment/length', 0.38, 0.36, 0.30]) lg2.append(['#eta correlations', 1.05, 0.12, 1.05]) legend2 = ['#rho (published)','#delta (published)','#xi (published)', \ '#rho (final)','#delta (final)', \ '#xi (final)'] systematics_hist(lg2, 'Systematics.eps', legend, 10.5) systematics_hist2(lg1[2:], lg2, 'Systematics_est.eps', legend2, 10.5) lro = [] # ldk.append([source, central value, total error, sys error]) lro.append(['Derenzo, 1969',0.7518, 0.0026, 0.0026, 0.0026]) lro.append(['Musser, 2005', 0.7508, 0.00105, 0.00032, 0.001]) lro.append(['MacDonald, 2008', 0.75014, 0.000484, 0.00017, 0.00044]) lro.append(['TWIST Final Results', 0.74991, 0.00029, 0.000089, 0.00028]) decay_par_graph(lro, 'rhoparvals.eps', 0.744, 0.756) ldt = [] ldt.append(['Balke, 1988', 0.7486, 0.0038, 0.0026, 0.0028]) ldt.append(['Gaponenko, 2005',0.74964, 0.0013, 0.00066, 0.00112]) ldt.append(['MacDonald, 2008', 0.75067, 0.000734, 0.00030, 0.00067]) ldt.append(['TWIST Final Results', 0.750723, 0.00038, 0.00016, 0.00029]) decay_par_graph(ldt, 'deltaparvals.eps', 0.744, 0.756, 0) lxi = [] lxi.append(['Beltrami, 1987', 1.0027, 0.00845, 0.0079, 0.0030]) lxi.append(['Jamieson, 2006', 1.0003, 0.0038, 0.0006, 0.0038]) lxi.append(['TWIST Final Results', 1.00109, 0.00119, 0.00035, 0.00114]) decay_par_graph(lxi, 'xiparvals.eps', 0.985, 1.015, 0) cplleg = ['PDG 2004', 'MacDonald 2008', 'Current Results'] lcc = [] lcc.append(['|g^{S}_{RR}|', 0.166, 0.062, 0.031]) lcc.append(['|g^{V}_{RR}|', 0.033, 0.031, 0.015]) lcc.append(['|g^{S}_{LR}|', 0.125, 0.074, 0.041]) lcc.append(['|g^{V}_{LR}|', 0.060, 0.025, 0.018]) lcc.append(['|g^{T}_{LR}|', 0.036, 0.021, 0.012]) lcc.append(['|g^{S}_{RL}|', 0.424, 0.412, 0.412]) lcc.append(['|g^{V}_{RL}|', 0.110, 0.104, 0.103]) lcc.append(['|g^{T}_{RL}|', 0.122, 0.104, 0.103]) lcc.append(['|g^{S}_{LL}|', 0.550, 0.550, 0.550]) systematics_hist(lcc, 'coupling_constants.eps', cplleg, 1) if __name__=='__main__': main(sys.argv[1:])