! plot_twist_lrs_limit set linthk -6 ! use 68.27%CL (one sigma) !cl = 0.6827 ! use 90%CL cl = 0.9 ! sigma corresponding to this cl ! usual gaussian, 1.64@90% sig_clg = sqrt(2)*aerf(cl) ! or one-sided limit, 1.28@90% sig_clg1 = sqrt(2)*aerf(1 - 2*(1 - cl)) scale_xlo = 200 scale_xhi = 1200 scale_ylo = -0.1 scale_yhi = 0.1 scale scale_xlo scale_xhi 5 scale_ylo scale_yhi 4 x_mass = [scale_xlo;scale_xhi] y_zeta = [scale_ylo;scale_yhi] ! ! rho ! ! Derenzo PR 181 (1969) 1954 rho_der = 0.7518 drho_der = 1.e-4*26. llim_rho_der = rho_der-sig_clg*drho_der ulim_rho_der = rho_der+sig_clg*drho_der ! lrs model, limit on zeta from Derenzo rho; rho>=0.75, so use one-sided limit llim_rho_der = 0.75-sig_clg1*drho_der ! note this doesn't use rho_der zeta_rho_der_1[1] = -sqrt((0.75 - llim_rho_der)*2./3.) zeta_rho_der_1[2] = zeta_rho_der_1[1] zeta_rho_der_2[1] = -zeta_rho_der_1[1] zeta_rho_der_2[2] = zeta_rho_der_2[1] ! musser et al. (2005) 0.7580 pm 0.00032 pm 0.00097 pm 0.00023 rho_twist_muss = 0.75080 drho_twist_muss = 1.e-4*sqrt(3.2**2 + 9.7**2 + 2.3**2) ! NSERC research proposal, 2004, for results projected in 2005 rho_twist_05p = 0.7500 drho_twist_05p = 1.e-4*sqrt(2.5**2 + 3.7**2 + 2.3**2) llim_rho_twist_05p = rho_twist_05p-sig_clg*drho_twist_05p ulim_rho_twist_05p = rho_twist_05p+sig_clg*drho_twist_05p ! lrs model, limit on zeta from final TWIST rho !llim_rho_twist_05p = 0.75-sig_clg1*drho_twist_05p !zeta_rho_twist_1_05p[1] = -sqrt((0.75 - llim_rho_twist_05p)*2./3.) !zeta_rho_twist_1_05p[2] = zeta_rho_twist_1_05p[1] !zeta_rho_twist_2_05p[1] = -zeta_rho_twist_1_05p[1] !zeta_rho_twist_2_05p[2] = zeta_rho_twist_2_05p[1] ! NSERC research proposal, 2004, for projected final results rho_twist_04fp = 0.7500 drho_twist_04fp = 1.e-4*sqrt(1.0**2 + 2.4**2 + 2.3**2) llim_rho_twist_04fp = rho_twist_04fp-sig_clg*drho_twist_04fp ulim_rho_twist_04fp = rho_twist_04fp+sig_clg*drho_twist_04fp ! lrs model, limit on zeta from final TWIST rho !llim_rho_twist_04fp = 0.75-sig_clg1*drho_twist_04fp !zeta_rho_twist_1_04fp[1] = -sqrt((0.75 - llim_rho_twist_04fp)*2./3.) !zeta_rho_twist_1_04fp[2] = zeta_rho_twist_1_04fp[1] !zeta_rho_twist_2_04fp[1] = -zeta_rho_twist_1_04fp[1] !zeta_rho_twist_2_04fp[2] = zeta_rho_twist_2_04fp[1] ! macdonald et al. (2008) 0.75014 pm 0.00017 pm 0.00044 pm 0.00011 rho_twist_macd = 0.75014 drho_twist_macd = 1.e-4*sqrt(1.7**2 + 4.6**2 + 1.1**2) llim_rho_twist_macd = rho_twist_macd-sig_clg*drho_twist_macd ulim_rho_twist_macd = rho_twist_macd+sig_clg*drho_twist_macd ! lrs model, limit on zeta from TWIST rho; rho>=0.75, so use one-sided limit llim_rho_twist_macd = 0.75-sig_clg1*drho_twist_macd ! note this doesn't use rho_twist zeta_rho_twist_macd_1[1] = -sqrt((0.75 - llim_rho_twist_macd)*2./3.) zeta_rho_twist_macd_1[2] = zeta_rho_twist_macd_1[1] zeta_rho_twist_macd_2[1] = -zeta_rho_twist_macd_1[1] zeta_rho_twist_macd_2[2] = zeta_rho_twist_macd_2[1] ! 2008 est, for final results !rho_twist_08fp = 0.7500 !drho_twist_08fp = 1.e-4*sqrt(1.3**2 + 2.4**2 + 1.1**2) rho_twist_08fp = rho_twist_macd drho_twist_08fp = drho_twist_macd llim_rho_twist_08fp = rho_twist_08fp-sig_clg*drho_twist_08fp ulim_rho_twist_08fp = rho_twist_08fp+sig_clg*drho_twist_08fp ! lrs model, limit on zeta from final TWIST rho llim_rho_twist_08fp = 0.75-sig_clg1*drho_twist_08fp zeta_rho_twist_1_08fp[1] = -sqrt((0.75 - llim_rho_twist_08fp)*2./3.) zeta_rho_twist_1_08fp[2] = zeta_rho_twist_1_08fp[1] zeta_rho_twist_2_08fp[1] = -zeta_rho_twist_1_08fp[1] zeta_rho_twist_2_08fp[2] = zeta_rho_twist_2_08fp[1] ! delta ! gaponenko et al. delta_twist_gap = 0.74964 ddelta_twist_gap = 1.e-4*sqrt(6.6**2 + 11.2**2) llim_delta_twist_gap = delta_twist_gap - sig_clg*ddelta_twist_gap ulim_delta_twist_gap = delta_twist_gap + sig_clg*ddelta_twist_gap ! NSERC research proposal, 2004, projected for 2005 result delta_twist_05p = 0.7500 ddelta_twist_05p = 1.e-4*sqrt(5.3**2 + 2.8**2) llim_delta_twist_05p = delta_twist_05p - sig_clg*ddelta_twist_05p ulim_delta_twist_05p = delta_twist_05p + sig_clg*ddelta_twist_05p ! NSERC research proposal, 2004, for final result delta_twist_04f = 0.7500 ddelta_twist_04f = 1.e-4*sqrt(2.2**2 + 3.2**2) llim_delta_twist_04f = delta_twist_04f - sig_clg*ddelta_twist_04f ulim_delta_twist_04f = delta_twist_04f + sig_clg*ddelta_twist_04f ! macdonald et al., (2008) 0.75067 pm 0.00030 pm 0.00067 delta_twist_macd = 0.75067 ddelta_twist_macd = 1.e-4*sqrt(3.0**2 + 6.7**2) llim_delta_twist_macd = delta_twist_macd - sig_clg*ddelta_twist_macd ulim_delta_twist_macd = delta_twist_macd + sig_clg*ddelta_twist_macd ! 2008 est, projected for final result !delta_twist_08fp = 0.7500 !ddelta_twist_08fp = 1.e-4*sqrt(2.3**2 + 2.2**2) delta_twist_08fp = delta_twist_macd ddelta_twist_08fp = ddelta_twist_macd llim_delta_twist_08fp = delta_twist_08fp - sig_clg*ddelta_twist_08fp ulim_delta_twist_08fp = delta_twist_08fp + sig_clg*ddelta_twist_08fp ! pmu*xi ! Beltrami et al., PLB 194, 326 (1987) pmuxi_bel = 1.0027 dpmuxi_bel = 1.e-4*sqrt(79**2 + 30**2) llim_pmuxi_bel = pmuxi_bel - sig_clg*dpmuxi_bel ulim_pmuxi_bel = pmuxi_bel + sig_clg*dpmuxi_bel ! Jamieson et al., 2006 pmuxi_twist_bj = 1.0003 dpmuxi_twist_bj = 1.e-4*sqrt(6**2 + 38**2) llim_pmuxi_twist_bj = pmuxi_twist_bj - sig_clg*dpmuxi_twist_bj ulim_pmuxi_twist_bj = pmuxi_twist_bj + sig_clg*dpmuxi_twist_bj ! estimated final result, 2008 pmuxi_twist_08fp = 1.000 !dpmuxi_twist_08fp = 1.e-4*sqrt(3.**2 + 5.9**2) dpmuxi_twist_08fp = 8.2e-4 llim_pmuxi_twist_08fp = pmuxi_twist_08fp - sig_clg*dpmuxi_twist_08fp ulim_pmuxi_twist_08fp = pmuxi_twist_08fp + sig_clg*dpmuxi_twist_08fp ! pmu*xi*delta/rho ! Jodidio et al., original PRD article pmuxidelrho_jod_o = 0.99863 dpmuxidelrho_jod_o = 1.e-4*sqrt(4.6**2 + 7.5**2) ! Jodidio et al., PRD erratum pmuxidelrho_jod = 0.99790 dpmuxidelrho_jod = 1.e-4*sqrt(4.6**2 + 7.5**2) ! calculated using one-sided limit llim_pmuxidelrho_jod = pmuxidelrho_jod - sig_clg1*dpmuxidelrho_jod ! 0.99677 ! from paper, including corrections to paper using muSR technique (Stoker) llim_pmuxidelrho_jod_stoker = 0.99682 ! Stoker et al., PRL dpmuxidelrho_sto = 1.e-4*sqrt(16**2 + 6**2 + 16**2) ! stat, correction, syst llim_pmuxidelrho_sto = 0.9955 ! work backwards to get 0.9985 as a central value pmuxidelrho_sto_o = llim_pmuxidelrho_sto + sig_clg1*dpmuxidelrho_sto ! but Stoker et al. had same error as Jodidio, so must be modified pmuxidelrho_sto = pmuxidelrho_sto_o*pmuxidelrho_jod/pmuxidelrho_jod_o ! lim_pmuxidelrho = 0.99682 from Jodidio (erratum), including Stoker ! find weighted sum var_jod = dpmuxidelrho_jod**2 var_sto = dpmuxidelrho_sto**2 var_js = 1./((1./var_jod) + (1./var_sto)) pmuxidelrho_js = var_js*(pmuxidelrho_jod/var_jod + pmuxidelrho_sto/var_sto) dpmuxidelrho_js = sqrt(var_js) ! result in paper is 0.99682; this gives 0.996828 (rounding?) lim_pmuxidelrho_js = pmuxidelrho_js - sig_clg1*dpmuxidelrho_js ! calculate pmuxi from this limit and twist rho, delta !pmuxi_twist_js = pmuxidelrho_js*rho_twist/delta_twist !dpmuxi_twist_js = pmuxi_twist_js* - ! sqrt((dpmuxidelrho_js/pmuxidelrho_js)**2 - ! + (drho_twist/rho_twist)**2 - ! + (ddelta_twist/delta_twist)**2 ) ! !llim_pmuxi_twist_js = pmuxi_twist_js - sig_clg*dpmuxi_twist_js !ulim_pmuxi_twist_js = pmuxi_twist_js + sig_clg*dpmuxi_twist_js ! ! calculate xi from Q_R^mu = 0.5*(1 + xi/3 -16*xi*delta/9) >=0 !ulim_xi_qr = 1./((16.*llim_delta_twist/9.) - (1./3.)) ! ! ! too frustrating... instead of calculation, ! put in numbers from Gaponenko explicitly llim_pmuxi_twist_js = 0.9960 ulim_pmuxi_twist_js = 1.0040 ! NSERC research proposal, 2004, for 2005 result !llim_pmuxi_twist_05p = 1.0 - sig_clg*1.e-4*sqrt(7.2**2 + 6.9**2) ! NSERC research proposal, 2004, for final result !llim_pmuxi_twist_04f = 1.0 - sig_clg*1.e-4*sqrt(3.0**2 + 3.0**2) ! return @pmuxi_lim llim_pmuxi_bel x_pmuxi_bel = massvar_nonmanif y_pmuxi_bel = zetavar_nonmanif x_pmuxi_bel_m = massvar_manif y_pmuxi_bel_m = zetavar_manif @pmuxidelrho_lim llim_pmuxidelrho_jod_stoker x_pmuxi_jod_stoker = massvar_nonmanif y_pmuxi_jod_stoker = zetavar_nonmanif x_pmuxi_jod_stoker_m = massvar_manif y_pmuxi_jod_stoker_m = zetavar_manif @pmuxi_lim llim_pmuxi_twist_js x_pmuxi_twist_js = massvar_nonmanif y_pmuxi_twist_js = zetavar_nonmanif x_pmuxi_twist_js_m = massvar_manif y_pmuxi_twist_js_m = zetavar_manif @pmuxi_lim llim_pmuxi_twist_bj x_pmuxi_twist_bj = massvar_nonmanif y_pmuxi_twist_bj = zetavar_nonmanif x_pmuxi_twist_bj_m = massvar_manif y_pmuxi_twist_bj_m = zetavar_manif @pmuxi_lim llim_pmuxi_twist_08fp x_pmuxi_twist_08fp = massvar_nonmanif y_pmuxi_twist_08fp = zetavar_nonmanif x_pmuxi_twist_08fp_m = massvar_manif y_pmuxi_twist_08fp_m = zetavar_manif !!@pmuxi_lim llim_pmuxi_twist_f !x_pmuxi_twist_f = massvar_nonmanif !y_pmuxi_twist_f = zetavar_nonmanif !x_pmuxi_twist_f_m = massvar_manif !y_pmuxi_twist_f_m = zetavar_manif ! set up plot shape, size, etc. set linthk -6 enable border orientation portrait window 21 0 20 100 100 !window 21 0 0 100 100 window 21 set xtica 90 set ytica -90 if (~exist(`letter_size')) then letter_size = 3 set %txthit letter_size set %xlabsz letter_size set %ylabsz letter_size set %xnumsz letter_size set %ynumsz letter_size set %xlaxis 25 set %ylaxis 25 pl_xlab = `(g<_>L<^>/g<_>R<^>)m<_>2<^> (GeV/c<^>2<_>)' pl_ylab = `(g<_>R<^>/g<_>L<^>)' pl_titl = `General LRS model' ! pl_titl = ` ' set histyp 0 label\xaxis pl_xlab label\yaxis pl_ylab clear legend frame 30 77 85 87 legend frame off !!! plot_beltrami: legend_text = - `Beltrami et al. P<_><^> ' !set lintyp 299 !color magenta !gra\noaxes x_pmuxi_bel y_pmuxi_bel !set lintyp 1 !color black !gra\noaxes x_pmuxi_bel y_pmuxi_bel set lintyp 4 color white legend on gra\noaxes legend_text x_pmuxi_bel y_pmuxi_bel legend off !!! plot_rho_derenzo: legend_text = - `Derenzo et al. ' !set lintyp 1 !color black !gra\noaxes x_mass zeta_rho_der_1 !gra\noaxes x_mass zeta_rho_der_2 set lintyp 4 color white legend on !JB gets rid of horizontal white line one !gra\noaxes legend_text x_mass zeta_rho_der_1 legend off !JB gets rid of horizontal white line two !gra\noaxes x_mass zeta_rho_der_2 goto plot_twist !!! plot_jodidio: legend_text = - `Jodidio et al. P<_><^>/ ' !set lintyp 282 !color green !gra\noaxes x_pmuxi_jod_stoker y_pmuxi_jod_stoker !set lintyp 1 !color black !gra\noaxes x_pmuxi_jod_stoker y_pmuxi_jod_stoker set lintyp 6 color green legend on gra\noaxes legend_text x_pmuxi_jod_stoker y_pmuxi_jod_stoker legend off !!! plot_twist: ! plot TWIST current plot_bj: legend frame 30 21 95 37 legend_text = - `TWIST P<_><^> (Jamieson) ' !set lintyp 277 !color blue !gra\noaxes legend_text x_pmuxi_twist_bj y_pmuxi_twist_bj !set lintyp 1 !color black !gra\noaxes x_pmuxi_twist_bj y_pmuxi_twist_bj set lintyp 3 color blue legend on gra\noaxes legend_text x_pmuxi_twist_bj y_pmuxi_twist_bj legend off !!! plot_twist_js: legend_text = - `TWIST , (MacDonald) + Jodidio' !set lintyp 255 !color blue !gra\noaxes x_pmuxi_twist_js y_pmuxi_twist_js !set lintyp 5 !color black !gra\noaxes x_pmuxi_twist_js y_pmuxi_twist_js set lintyp 5 color blue legend on gra\noaxes legend_text x_pmuxi_twist_js y_pmuxi_twist_js legend off !!! plot_rho_twist: legend_text = - `TWIST ' set lintyp 5 color blue !legend on !gra\noaxes legend_text x_mass zeta_rho_twist_macd_1 !JB: following two lines do horizontal blue !gra\noaxes x_mass zeta_rho_twist_macd_1 !gra\noaxes x_mass zeta_rho_twist_macd_2 !legend off !!! plot_08fp: ! plot TWIST expected final values as of 2008 legend_text = - `TWIST , (MacDonald) + P<_><^> this measurement' !`TWIST P<_><^> estimated ' !set lintyp 277 !color white !gra\noaxes legend_text x_pmuxi_twist_bj y_pmuxi_twist_bj !set lintyp 1 !color black !gra\noaxes x_pmuxi_twist_bj y_pmuxi_twist_bj set lintyp 1 color red legend on gra\noaxes legend_text x_pmuxi_twist_08fp y_pmuxi_twist_08fp legend off legend_text = - `TWIST ' set lintyp 1 color red !JB: following two lines do horizontal red !gra\noaxes x_mass zeta_rho_twist_1_08fp !gra\noaxes x_mass zeta_rho_twist_2_08fp !set lintyp 1 !color white !legend on !gra\noaxes legend_text x_mass zeta_rho_twist_1_08fp !gra\noaxes x_mass zeta_rho_twist_1_08fp !legend off !gra\noaxes x_mass zeta_rho_twist_2_08fp ! set >0 if you don't want the non_manifest LRS contours plotted no_non_manifest = 1 if (no_non_manifest > 0) then goto no_non_manifest_1 ! now plot non-manifest profile vector spacevec 1 spacevec = [.8] set hatch 2 spacevec 45 !set lintyp 277 set lintyp 102 legend_text = `TWIST P<_><^> (NMLRS), final ' !gra\noaxes x_pmuxi_twist y_pmuxi_twist legend_text = - `TWIST P<_><^> (V<^>R<_><_>ud<^>=V<^>L<_><_>ud<^>)' color blue gra\noaxes x_pmuxi_twist y_pmuxi_twist set lintyp 1 color black gra\noaxes x_pmuxi_twist y_pmuxi_twist set lintyp 5 color blue gra\noaxes x_pmuxi_twist y_pmuxi_twist no_non_manifest_1: window\noconfirm 99 70 64 95 67 erasewindow window 21 color white set %xloc 90 set %yloc 57 set linthk -10 set cursor -11 text\noconfirm `allowed' graph\axesonly x_mass zeta_rho_twist_2_08fp