
<<VplusA_ALEPH;

mttau=1.77686; (* This is PDG2016 value; I did this change on Dec.3, 2016 *) 
(* mttau=1.77682; tau mass, in GeV PDG2014 value, used here until Dec.2, 2016 *)
smaxtau=3.07854 + 0.157089/2;  (* = mttau^2 *)

BBe=0.17827;
VudCKM=0.97425;
SEW=1.0201;
mPi=0.13957018;  (* mass of charged pion, in GeV *)
fPismall=0.09221;  (* (small) decay constant of pion; *)
fPi=Sqrt[2.] fPismall; (* (large) decay constant of pion *)

(* On 6 April 2017, I change the data to: *)
mPi=0.13957;
fPi=0.13050;
(* In this way the same data as for OPAL; from PRD 2016 *)
(* The value of fPi=0.13050; (and not:0.130405) is from  [http://pdg.lbl.gov/2016/reviews/rpp2016-rev-pseudoscalar-meson-decay-cons.pdf], their Eq.(6), and PDG 2016 value Vud=0.97417 +- 0.00021, giving fPi=(0.12713/0.97417)=0.130501 GeV, with uncertainties +-0.00014 GeV *) 

wt[sig_]:=(1 - sig/mttau^2)^2 (1 + 2 sig/mttau^2);

Do[dwt[j]=NIntegrate[wt[sig],{sig,sbin[[j]]-dsbin[[j]]/2,sbin[[j]]+dsbin[[j]]/2}],{j,1,79}];
(* Note, the s=sigma inteval j is: (sbin[[j]]-dsbin[[j]]/2,sbin[[j]]+dsbin[[j]]/2), j=1,2,...,79 *)
(* sigma_max=sbin[[79]]+dsbin[[79]]/2 =3.1570845  = mttau^2; where mttau=1.77682; *)
   
ffactor=100 12 Pi^2 VudCKM^2 SEW BBe/mttau^2;
resc=0.9987;
(* resc is the additional rescaling factor for spectral functions, needed to update the ALEPH results -- see Santi peris email of Nov.8, 2016, and their paper (with Boito et al.) PRD91,034003 (2015) [arXiv:1410.3528] Sec. III.A *) 
(* the covariance matrix elements are multiplied by the correction factor resc^2 *)

Do[rhoVplA[j]=resc sfm2[[j]]/ffactor/dwt[j],{j,1,79}]; 
Do[delrhoVplA[j]=resc derr[[j]]/ffactor/dwt[j],{j,1,79}]; 

(* covAL[ii,jj] will be the ALEPH 160x160 covariance mtrx *)
Do[covAL[ii,jj]=covj[jj][[ii]],{ii,1,160},{jj,1,160}];
Do[covALVV[ii,jj]=covAL[ii,jj],{ii,1,79},{jj,1,79}];
Do[covALAA[ii,jj]=covAL[80+ii,80+jj],{ii,1,79},{jj,1,79}];
Do[covALVA[ii,jj]=covAL[ii,80+jj],{ii,1,79},{jj,1,79}];
Do[covALAV[ii,jj]=covAL[80+ii,jj],{ii,1,79},{jj,1,79}];

(* Note: derr[[jj]]*0.01=Sqrt[covALVV[jj,jj]+covALAA[jj,jj] + 2 covALAV[jj,jj]]; for any jj=1,2,...,79 *)

Do[OmegaVplA[j]= 2 Pi^2 rhoVplA[j],{j,1,79}];
Do[delOmegaVplA[j]= 2 Pi^2 delrhoVplA[j],{j,1,79}]; 

(* the covariance matrix elements are multiplied by the correction factor resc^2: *)
Do[covVV[ii,jj]=resc^2 10^4 covALVV[ii,jj]/ffactor^2/dwt[jj]/dwt[ii],{ii,1,79},{jj,1,79}];
Do[covAA[ii,jj]=resc^2 10^4 covALAA[ii,jj]/ffactor^2/dwt[jj]/dwt[ii],{ii,1,79},{jj,1,79}];
Do[covVA[ii,jj]=resc^2 10^4 covALVA[ii,jj]/ffactor^2/dwt[jj]/dwt[ii],{ii,1,79},{jj,1,79}];
Do[covAV[ii,jj]=resc^2 10^4 covALAV[ii,jj]/ffactor^2/dwt[jj]/dwt[ii],{ii,1,79},{jj,1,79}];
(* Check this *)
(* covVV, covAA, covVA, covAV are the VV, AA, AV, and VA covariance matrices, for ALEPH case *)

(* Note: covVA[i,j]=covAV[j,i], as it should be *)
 
delOmegaALCovBin[ii_]:= 2 Pi^2 Sqrt[(  covVV[ii,ii]+covAA[ii,ii]+covVA[ii,ii] + covAV[ii,ii] )];
(* Note: delOmegaALCovBin[j]=delOmegaVplA[j], j=1,2,...79; (to nonzero 5 digits coincide) as it should be *)

OmegaVplAint=Table[{sbin[[j]],OmegaVplA[j]},{j,1,79}];
OmegaVplAf[sig_]:=Interpolation[OmegaVplAint][sig];
delOmegaVplAint=Table[{sbin[[j]],delOmegaVplA[j]},{j,1,79}];
delOmegaVplAf[sig_]:=Interpolation[delOmegaVplAint][sig];

(* plOmegaVplANoPiALEPH=Plot[{OmegaVplAf[sig], OmegaVplAf[sig] + delOmegaVplAf[sig], OmegaVplAf[sig] -  delOmegaVplAf[sig]}, {sig, 0, 3.12}, PlotRange -> {{0, 3.16}, {0., 2.8}}, PlotStyle -> {Black, Dashed, Dashed}, Frame -> True, FrameLabel -> {\[Sigma][Superscript[GeV, 2]], Subscript[\[Omega],V+A]}]; *)

(* NOw, for more concise plot, as ListPlot in Bins, we create tables: *)
OmegaTabAL=Table[{sbin[[jj]],OmegaVplA[jj]},{jj,1,79}];
delOmegaTabAL=Table[{sbin[[jj]],delOmegaALCovBin[jj]},{jj,1,79}];
OmegaTabALPlus= Table[{sbin[[jj]],OmegaVplA[jj] +delOmegaALCovBin[jj]},{jj,1,79}];
OmegaTabALMinus=Table[{sbin[[jj]],OmegaVplA[jj] -delOmegaALCovBin[jj]},{jj,1,79}];

liplOmegaVplANoPiALEPH = 
 ListPlot[{OmegaTabAL, OmegaTabALPlus, OmegaTabALMinus}, 
  PlotRange -> {{0, 3.12}, {0., 2.8}}, Filling -> {3 -> {2}}, 
  PlotStyle -> {{Black, PointSize[0.01]}, {Red, 
     PointSize[0.006]}, {Red, PointSize[0.006]}}, 
  FillingStyle -> {Black, Thickness[0.01]}, Frame -> True, 
	  FrameLabel -> {\[Sigma][Superscript[GeV, 2]], 
	  Subscript[\[Omega], V + A]},
	  Epilog -> {Inset[
			   Style[ALEPH, 11], Scaled[{0.70, 0.85}]]}];

(* Export["OmegaVANoPiALEPH.pdf", liplOmegaVplANoPiALEPH] *)


NNAL=79; (* This is the total number of bins; the 80th bin is empty, and entirely at sig>mtau^2 (which is weird that they included) *)
(*----------------------------------------_*)

BexpBinAL[M2_? NumberQ]:=(1/M2) (Sum[ Exp[-sbin[[ii]]/M2]*OmegaVplA[ii]*dsbin[[ii]],{ii,1,NNAL}] + 2 Pi^2 fPi^2 Exp[-mPi^2/M2]);
ReBexpBinAL[M2_]:=Re[(1/M2) (Sum[ Exp[-sbin[[ii]]/M2]*OmegaVplA[ii]*dsbin[[ii]],{ii,1,NNAL}] + 2 Pi^2 fPi^2 Exp[-mPi^2/M2])];
(* ReBexpBin[M2] is the integral of R when the intergation is performed as a sum over the OPAL bins; *)
 
   (* BexpBinALb[M2_? NumberQ]:=(1/M2) (Sum[ NIntegrate[Exp[-sbin[[ii]]/M2]*OmegaVA[ii],{sig, sbin[[ii]] - dsbin[[ii]]/2, sbin[[ii]] + dsbin[[ii]]/2}],{ii,1,NNAL}] + 2 Pi^2 fPi^2 Exp[-mPi^2/M2]);  -- gives the same result, but ALb is in principle 'better' *)

delReBexpBinAL[M2_]:=(2 Pi^2) Sqrt[Sum[Re[Exp[-sbin[[ii]]/M2]/M2] Re[Exp[-sbin[[jj]]/M2]/M2] dsbin[[ii]] dsbin[[jj]] (  covVV[ii,jj]+covAA[ii,jj]+covVA[ii,jj]+ covAV[ii,jj]),{ii,1,NNAL},{jj,1,NNAL}]];
(* NOte: delReBexpBinAL[M2] is the experimental uncertainty of ReBexpBinAL[M2]; *)
(* BTs based on ALEPH data: Borel Transforms generated by prg. expdataALEPH.m: (Re)BexpBin[M2] and delReBexpBin[M2] as function of (complex) Borel parameter M2 *)
(* In delReBexpBinAL[M2] (2 Pi^2) factor included *)
(* Note that our "observable" is Re[BT[M2]], and not BT[M2]; so, sigma of Re[BT[M2]] is then: delReBexpBin *)

(*----------------------------------------_*)
(*   plBexpBinALPsi0= Plot[BexpBinAL[M2], {M2, 0.65, 1.5}, PlotRange -> {{0.65, 1.5}, {1.05, 1.30}}, AxesOrigin -> {0.65,1.0}, PlotStyle -> {Black}, Frame -> True, FrameLabel -> {Superscript[Abs[M], 2] [Superscript[GeV, 2]], Re Superscript[Subscript[B, exp],(bin)][Superscript[M, 2]]}, Epilog -> {Inset[Style[\[Psi]==0, 12], Scaled[{0.85, 0.85}]]}]; *)

(* plBexpBinALPsiPi4= Plot[Re[BexpBinAL[M2 Exp[I Pi/4]]], {M2, 0.65, 1.5}, PlotRange -> {{0.65, 1.5}, {1.15, 1.32}}, AxesOrigin -> {0.65,1.0}, PlotStyle -> {Black}, Frame -> True, FrameLabel -> {Superscript[Abs[M], 2] [Superscript[GeV, 2]], Re  Superscript[Subscript[B, exp],(bin)][Superscript[M, 2]]}, Epilog -> {Inset[Style[\[Psi]==\[Pi]/4, 12], Scaled[{0.85, 0.85}]]}]; *)

(* plBexpBinALPsiPi6= Plot[Re[BexpBinAL[M2 Exp[I Pi/6]]], {M2, 0.65, 1.5}, PlotRange -> {{0.65, 1.5}, {1.09, 1.30}}, AxesOrigin -> {0.65,1.0}, PlotStyle -> {Black}, Frame -> True, FrameLabel -> {Superscript[Abs[M], 2] [Superscript[GeV, 2]], Re Subscript[B, exp][Superscript[M, 2]]}, Epilog -> {Inset[Style[\[Psi]==\[Pi]/6, 12], Scaled[{0.85, 0.85}]]}]; *)

 (* Export["BexpBinALPsi0.pdf",plBexpBinALPsi0]; Export["BexpBinALPsiPi4.pdf",plBexpBinALPsiPi4]; Export["BexpBinALPsiPi6.pdf",plBexpBinALPsiPi6]; *)
    
(* Save["BexpBinAL.save",{BexpBinAL}] *)
 (*-----------------------------------------------------*)
  (*-----------------------------------------------------*)
  (*-----------------------------------------------------*)
(* Similary as ReBexpBin[M2] and delReBexpBin[M2] is calculated, we can calculare rtauexp (by changing: ReExp[-sigf[ii]/M2]/M2] -> (2/sigmax)(1 - sigf[ii]/mtau2)^2 (1 + 2 sigf[ii]/mtau2): *)
(* sigmax=2.80; *)
 (* rtauexp =appr (2/sigmax)  NIntegrate[(1 - sig/sigmax)^2 (1 + 2 sig/sigmax) OmegaVplAf[sig], {sig, 0, sigmax}] +  (2/sigmax) 2 Pi^2 fPi^2 (1 - mPi^2/sigmax)^2 (1 + 2 mPi^2/sigmax) - 1; *)
 (*---------------------------------_*)  
 

