#!/bin/env python # sw_plot_glrs.py # Begun: 20100204 # A Python/ROOT script to plot an exclusion plot for # LRS model parameters # PYTHON IMPORTS import os, sys # for exit etc import math from array import array from ROOT import * # from optparse import OptionParser mwl = 80.398 # mass of W_L, GeV/c^2 cabb = 0.97 ############################################################################# ## PLOT THE GRAPHS ############################################################################# def initialize(canvas_name="",wcanvas=800,hcanvas=600): print "Initializing ROOT canvas... ", gROOT.Reset() c1 = TCanvas(canvas_name, canvas_name, wcanvas, hcanvas) c1.SetWindowSize(wcanvas + (wcanvas - c1.GetWw()), hcanvas + (hcanvas - c1.GetWh())) c1.SetGridx() c1.SetGridy() #print "Done." return c1 def makegstyle(cStyle): # ROOT documentation TStyle.html, TStyle.h.html cStyle.SetCanvasBorderMode(0) cStyle.SetCanvasColor(10) ## Background color (white) cStyle.SetCanvasDefH(800) cStyle.SetCanvasDefW(600) cStyle.SetPadColor(10) ## Background color (white) cStyle.SetPadBorderSize(1) cStyle.SetPadBorderMode(0) #cStyle.SetTitleFillColor(10) ## Background color (white) cStyle.SetTitleBorderSize(0) ## cStyle.SetOptStat(0) ## Hide stats #cStyle.SetOptTitle(0) ## Do not display histogram titles. #cStyle.SetStatColor(10) ## Background color (white) #cStyle.SetPalette(1) ## Nicer colour scale for 2D histograms. cStyle.SetPaperSize(20,24) ## same as kUSLetter, cm with some border #cStyle.SetPaperSize(20,20) ## same as kUSLetter, cm with some border #cStyle.SetStripDecimals(False) # for publishing cStyle.SetFrameLineWidth(4) cStyle.SetLineWidth(2) return cStyle ## for pmu*xi (half-)limit; note coefficient of zeta^2 is 2 ## (1_r[0]^2/cos^2(theta_cabb)*t^2 + 2*t*zeta*(r[0]/theta_cabb)*r[1] ## + 2*zeta^2 = hlim ## r[0] = cos(theta_R) ## r[1] = cos(alpha + omega) (CP violation) ## cos(theta_cabb) = 0.97, Cabbibbo angle ## t = ((gL*m1)/gR*m2))^2 def swlrs_pmuxi_max(x,r,hlim): rc = r[0]/cabb a = 2. b = 2.*x*rc*r[1] c = (1. + rc*rc)*x*x - hlim det = (b*b) - (4.*a*c) if (det < 0.): return float("nan") elif (det == 0.): return -1.*b/(2.*a) else: return (-1.*b + sqrt(det))/(2.*a) def swlrs_pmuxi_min(x,r,hlim): rc = r[0]/cabb a = 2. b = 2.*x*rc*r[1] c = (1. + rc*rc)*x*x - hlim det = (b*b) - (4.*a*c) if (det < 0.): return float("nan") elif (det == 0.): return -1.*b/(2.*a) else: return (-1.*b - sqrt(det))/(2.*a) ## for pmu*xi*delta/rho (half-)limit; note coefficient of zeta^2 is 1 ## (1+r[0]^2/cos^2(theta_cabb))*t^2 + 2*t*zeta*(r[0]/theta_cabb)*r[1] ## + zeta^2 = hlim ## r[0] = cos(theta_R) ## r[1] = cos(alpha + omega) (CP violation) ## cos(theta_cabb) = 0.97, Cabbibbo angle ## t = ((gL*m1)/gR*m2))^2 def swlrs_pmxdr_max(x,r,hlim): rc = r[0]/cabb a = 1. b = 2.*x*rc*r[1] c = (1. + rc*rc)*x*x - hlim det = (b*b) - (4.*a*c) if (det < 0.): return float("nan") elif (det == 0.): return -1.*b/(2.*a) else: return (-1.*b + sqrt(det))/(2.*a) def swlrs_pmxdr_min(x,r,hlim): rc = r[0]/cabb a = 1. b = 2.*x*rc*r[1] c = (1. + rc*rc)*x*x - hlim det = (b*b) - (4.*a*c) if (det < 0.): return float("nan") elif (det == 0.): return -1.*b/(2.*a) else: return (-1.*b - sqrt(det))/(2.*a) if __name__ == '__main__': llcomp = [] ## pre-TWIST, Beltrami pmx_bel_ll = 0.98765 llcomp_bel = 1. - pmx_bel_ll llcomp.append(llcomp_bel) ## pre-TWIST, Jodidio pmxdr_jod_ll = 0.99682 llcomp_jod = 1. - pmxdr_jod_ll llcomp.append(llcomp_jod) ## this is Carl's estimate for TWIST 2010 (xi90) pmx_tw_ll = 0.99922 llcomp_tw = 1. - pmx_tw_ll ## this is Carl's estimate for TWIST 2010 (xdr90) #pmxdr_tw_ll = 0.999564 ## this is Carl's estimate after considering rho correlations ## 20100208, p. 4 pmxdr_tw_ll = 0.99962 #llcomp_tw = 1. - pmxdr_tw_ll llcomp.append(llcomp_tw) r = [] r.append(1.) r.append(1.) ntries = 1200 xlrs = [] xlrslo = [] xlrshi = [] ylrs = [] ylrslo = [] ylrshi = [] npts = [] ncurves = 3 for i in range(0,ncurves): xlrs.append(array('d')) xlrslo.append(array('d')) xlrshi.append(array('d')) ylrs.append(array('d')) ylrslo.append(array('d')) ylrshi.append(array('d')) npts.append(0.) ######## load data for plots ic = 0 print "<<< graph ", ic, 1.-llcomp[ic] hlim = llcomp[ic]/2. print ">>>" mwr_min = 10000. zeta_min = 0. zeta_max = 0. npts[ic] = 0 for i in range(ntries): mass = 1.*(i+1) #mass = 635. + 1.*(i+1) #r[0] = cabb ##MLRS, t_theta=t since cos_theta_C(R)=cos_theta_C(L) #r[1] = 1. ##MLRS, cos(alpha+omega) = 1 since no CP violation x = (mwl/mass)**2 zeta_hi_max = -0.1 zeta_lo_min = 0.1 mass_for_zeta = -1. for ir0 in range(-10,11): r[0] = 0.1*ir0/cabb for ir1 in range(-10,11): # find min and max for 21 values of each cosine r[1] = 0.1*ir1 zeta_hi = swlrs_pmuxi_max(x,r,hlim) if (str(zeta_hi) == 'nan'): continue ## this should work in general zeta_hi_max = max(zeta_hi_max,zeta_hi) zeta_lo = swlrs_pmuxi_min(x,r,hlim) if (str(zeta_lo) == 'nan'): continue zeta_lo_min = min(zeta_lo_min,zeta_lo) mass_for_zeta = mass #print i, mass, mass_for_zeta, x, r[0], r[1], zeta_hi, zeta_lo if (mass_for_zeta < 0.): continue mwr_min = min(mass_for_zeta,mwr_min) zeta_min = min(zeta_lo_min,zeta_min) ## store min and max for all masses zeta_max = max(zeta_hi_max,zeta_max) #load values into plot arrays xlrslo[ic].append(mass_for_zeta) xlrshi[ic].append(mass_for_zeta) ylrslo[ic].append(zeta_lo_min) ylrshi[ic].append(zeta_hi_max) npts[ic] = npts[ic] + 1 #print "===", npts[ic], mass_for_zeta, mwr_min, zeta_lo_min, zeta_hi_max print ic, npts[ic], zeta_min, zeta_max, mwr_min xlrslo[ic].reverse() xlrs[ic] = xlrslo[ic] + xlrshi[ic] ylrslo[ic].reverse() ylrs[ic] = ylrslo[ic] + ylrshi[ic] ### ic = 1 print "<<< graph ", ic, 1.-llcomp[ic] hlim = llcomp[ic]/2. print ">>>" mwr_min = 10000. zeta_min = 0. zeta_max = 0. npts[ic] = 0 for i in range(ntries): mass = 1.*(i+1) #mass = 635. + 1.*(i+1) #r[0] = cabb ##MLRS, t_theta=t since cos_theta_C(R)=cos_theta_C(L) #r[1] = 0. ##MLRS, cos(alpha+omega) = 0 since no CP violation x = (mwl/mass)**2 zeta_hi_max = -0.1 zeta_lo_min = 0.1 mass_for_zeta = -1. for ir0 in range(-10,11): r[0] = 0.1*ir0/cabb for ir1 in range(-10,11): # find min and max for 21 values of each cosine r[1] = 0.1*ir1 zeta_hi = swlrs_pmxdr_max(x,r,hlim) if (str(zeta_hi) == 'nan'): continue ## this should work in general zeta_hi_max = max(zeta_hi_max,zeta_hi) zeta_lo = swlrs_pmxdr_min(x,r,hlim) if (str(zeta_lo) == 'nan'): continue zeta_lo_min = min(zeta_lo_min,zeta_lo) mass_for_zeta = mass #print i, mass, mass_for_zeta, x, r[0], r[1], zeta_hi, zeta_lo if (mass_for_zeta < 0.): continue mwr_min = min(mass_for_zeta,mwr_min) zeta_min = min(zeta_lo_min,zeta_min) ## store min and max for all masses zeta_max = max(zeta_hi_max,zeta_max) #load values into plot arrays xlrslo[ic].append(mass_for_zeta) xlrshi[ic].append(mass_for_zeta) ylrslo[ic].append(zeta_lo_min) ylrshi[ic].append(zeta_hi_max) npts[ic] = npts[ic] + 1 #print "===", npts[ic], mass_for_zeta, mwr_min, zeta_lo_min, zeta_hi_max print ic, npts[ic], zeta_min, zeta_max, mwr_min xlrslo[ic].reverse() xlrs[ic] = xlrslo[ic] + xlrshi[ic] ylrslo[ic].reverse() ylrs[ic] = ylrslo[ic] + ylrshi[ic] ### ic = 2 print "<<< graph ", ic, 1.-llcomp[ic] hlim = llcomp[ic]/2. print ">>>" mwr_min = 10000. zeta_min = 0. zeta_max = 0. npts[ic] = 0 for i in range(ntries): mass = 1.*(i+1) #mass = 635. + 1.*(i+1) #r[0] = cabb ##MLRS, t_theta=t since cos_theta_C(R)=cos_theta_C(L) #r[1] = 1. ##MLRS, cos(alpha+omega) = 1 since no CP violation x = (mwl/mass)**2 zeta_hi_max = -0.1 zeta_lo_min = 0.1 mass_for_zeta = -1. for ir0 in range(-10,11): r[0] = 0.1*ir0/cabb for ir1 in range(-10,11): # find min and max for 21 values of each cosine r[1] = 0.1*ir1 zeta_hi = swlrs_pmuxi_max(x,r,hlim) if (str(zeta_hi) == 'nan'): continue ## this should work in general zeta_hi_max = max(zeta_hi_max,zeta_hi) zeta_lo = swlrs_pmuxi_min(x,r,hlim) if (str(zeta_lo) == 'nan'): continue zeta_lo_min = min(zeta_lo_min,zeta_lo) mass_for_zeta = mass #print i, mass, mass_for_zeta, x, r[0], r[1], zeta_hi, zeta_lo if (mass_for_zeta < 0.): continue mwr_min = min(mass_for_zeta,mwr_min) zeta_min = min(zeta_lo_min,zeta_min) ## store min and max for all masses zeta_max = max(zeta_hi_max,zeta_max) #load values into plot arrays xlrslo[ic].append(mass_for_zeta) xlrshi[ic].append(mass_for_zeta) ylrslo[ic].append(zeta_lo_min) ylrshi[ic].append(zeta_hi_max) npts[ic] = npts[ic] + 1 #print "===", npts[ic], mass_for_zeta, mwr_min, zeta_lo_min, zeta_hi_max print ic, npts[ic], zeta_min, zeta_max, mwr_min xlrslo[ic].reverse() xlrs[ic] = xlrslo[ic] + xlrshi[ic] ylrslo[ic].reverse() ylrs[ic] = ylrslo[ic] + ylrshi[ic] ############# gStyle = makegstyle(gStyle) gStyle.SetOptTitle(0) tmarg = 0.1 # fractional margin sizes rmarg = 0.1 bmarg = 0.20 lmarg = 0.20 clrs = initialize("clrs",600,600) clrs.SetTopMargin(tmarg) clrs.SetBottomMargin(bmarg) clrs.SetLeftMargin(lmarg) clrs.SetRightMargin(rmarg) glrs0 = TGraph(2*npts[0],xlrs[0],ylrs[0]) glrs0.SetLineColor(1) #glrs0.SetFillColor(1) gStyle.SetLineStyleString(11,"60 20") glrs0.SetLineStyle(11) glrs0.SetLineWidth(int(4)) #glrs0.SetFillStyle(3005) glrs0.GetXaxis().SetLimits(0.,1200.) #glrs0.GetXaxis().SetLimits(635.,645.) glrs0.GetXaxis().SetTitle("(g_{L}/g_{R})m_{2} (GeV/c^{2})") #glrs0.GetXaxis().SetNdivisions(int(404),False) glrs0.GetXaxis().SetTitleSize(0.06) glrs0.GetXaxis().SetTitleOffset(1.1) glrs0.GetXaxis().SetLabelSize(0.04) glrs0.GetXaxis().CenterTitle(True) glrs0.GetYaxis().SetTitleSize(0.06) glrs0.GetYaxis().SetTitleOffset(1.4) glrs0.GetYaxis().SetLabelSize(0.04) glrs0.GetYaxis().CenterTitle(True) glrs0.GetYaxis().SetTitle("Mixing angle (g_{R}/g_{L})#zeta") glrs0.GetYaxis().SetRangeUser(-0.1,0.1) glrs1 = TGraph(2*npts[1],xlrs[1],ylrs[1]) glrs1.SetLineColor(4) glrs1.SetLineStyle(7) #glrs1.SetFillColor(4) glrs1.SetLineWidth(int(4)) #glrs1.SetFillStyle(3004) glrs2 = TGraph(2*npts[2],xlrs[2],ylrs[2]) glrs2.SetLineColor(2) glrs2.SetLineStyle(1) #glrs2.SetFillColor(2) glrs2.SetLineWidth(int(4)) #glrs2.SetFillStyle(3005) lrsleg = TLegend(lmarg+0.01,bmarg+0.01,lmarg+0.4,bmarg+0.1) lrsleg.SetBorderSize(0) lrsleg.SetFillColor(10) lrsleg.AddEntry(glrs0,"Beltrami 1987 P_{#mu}#xi (LRS)","l") lrsleg.AddEntry(glrs1,"Jodidio 1988 P_{#mu}#xi#delta/#rho (LRS)","l") lrsleg.AddEntry(glrs2,"TWIST 2010 P_{#mu}#xi (LRS)","l") lrs = TMultiGraph() #lrs.SetTitle("LRS") lrs.Add(glrs0,"AC") lrs.Add(glrs1,"C") lrs.Add(glrs2,"C") lrs.Draw() lrsleg.Draw() ask_user = raw_input('press ENTER to quit >')