
(* SumRMSbALEPHp3ka1ka2O6ka2.m: summary of the program: *)
(* The experimental data are of ALEPH, for the (V+A)-channel; the number of experimental bins is 77, i.e., the last two experimental bins are not included.  *) 
(*  D=0 part of Adler function, d(Q^2)_{D=0}, is generated by the 'p3ka1ka2' form of Borel of the auxiliary quantity td(Q^2), i.e., the Borel transform B[td](u) has, among other terms (u=2, u=-1 renormalon terms): two u=3 terms, ~1/(3-u)^ka1 and ~1/(3-u)^ka2, with: ka1=0.375; ka2=0.778; this is as suggested by the work of Boito, Hornung and Jamin (BHJ), 1510.03812v2; OPE is taken with zero anom. dimensions, except for D=6 term which has the form ~ a(Q^2)^k2 <O6>/(Q^2)^3 with k2=1-ka2 =0.222; we take d4=(275 +- 63) for the coefficient at a^5=(alpha_s(Q^2)/Pi)^5 in d(Q^2)_{D=0}; the renormalisation scale parameter kap=mu^2/Q^2 varies in the range kap=1^{+1}_{-1/2}; experimentral error of extracted parameters takes into account the correlation matrix (77x77) of the data. *)
(*--------------------------------------------------*)
(* SumRMSbALEPHp3ka1ka2O6ka2.m: as SumRMSbALEPHp3ka2O6ka2.m, only that now the D=0 part is taken to be the 'p3ka1ka2' case, i.e., B[td](u) has two u=3 terms, 1/(3-u)^ka1 and 1/(3-u)^ka2, (i.e., the D=0 part is better behaved than in the 'p3ka2p41' case where B[td](u) has 1/(3-u)^ka2 and 1/(4-u)) *) 
(*--------------------------------------------------*)
(*  SumRMSbALEPHp3ka2O6ka2.m: as  SumRMSbALEPHO6ka12Err.m, but with D=0 Adler function u=3 IR sector having only ka2=0.778(=k3) term (B[td](u) ~ 1/(3-u)^ka2) (ONE term, 'p3ka2'), no ka1=0.375 term in  B[td](u) (i.e., no 1/(3-u)^ka1 term); and in OPE D=6 contribution of V+A channel Adler function we have only ONE condensate ~ a(Q^2)^(1-ka2) <O6>/Q^6 ('O6ka2') *)
(*--------------------------------------------------*)
(* SumRMSbALEPHO6ka12Err.m: as SumRMSbALEPHO6ka12.m, but with more detailed experimental error analysis of the extracted values *)
(*--------------------------------------------------*)
(* SumRMSbALEPHO6ka12.m: here we use NchP=2 [Nd=3*(NchP+1)=9 points in the complex M^2-place] in chi^2 of Borel-Laplace sum rules that we minimise; all collateral test parts are deleted to make the program shorter *) 
(* SumRMSbALEPHO6ka12.m: as SumMSbALEPHM2O61O62.m (d4MSbar=275.+ -63.), but now the IR p=3 part of B[td][u]  has 1/(3-u)^ka2 and 1/(3-u)^ka1, where: ka2= 0.778 =app (0.8024+0.7527)/2, and ka1 = 0.375 =app (0.4173+0.3320)/2 (as suggested by the work of BHJ (Boito, Hornung and Jamin, 1510.03812v2); in SumMSbALEPHM2O61O62.m, we had large-beta0 values ka2=2 and ka1=1 *)
(* SumMSbALEPHM2O61O62.m: as SumRMSbALEPHM2O61O62d4338.m, but for: d4MSbar=275.+ -63.19, and not 338.19+-63.19 *)
(* Note, the string 'OK' indicates that the part there has been converted to the deMSbar=275 case (from d4MSbar=338.19 case) *)
(*-------------------------------------------------*)
(* SumRMSbALEPHM2O61O62d4338.m: as SumRPMSbALEPHM2m090.m, but instead of O6 and O8, we have O6^{(1)} and O6^{(2)} condensates *)
(* Look for: NEW *)
(*-------------------------------------------------*)
(* SumRPMSbALEPHM2m090.m:  with (double) pinched ('P') Borel-Laplace sum rules *)
(* and alpha_s(mtau^2; MSbar) is adjustable (as input, or as result fo global fit) *)
(*-------------------------------------------------*)
(* Borel-Laplace fitted in |M^2|-intervals: (M2min,1.5); M2min=0.65 or 0.90; *)
(* here: M2min=0.90, i.e., YYY=090 *)
(* in pQCD, MSbar scheme *)

(* This prg. calculates BT sum rule, when in Adler function MSbar pQCD coupling *)
(* aMSpert is used, for ALEPH data; either Adler is TPS *)
(* (truncated series in power of aMSpert, up to aMSpert^5 or higher order, *)
 (* either CIPT, or FOPT, or PV Borel-resumation, or renormalon-motivated sum *)
(*-------------------------------------------------*)

 (* PDG2021, average from hadronic tau decays: alpha_s(MZ^2)=0.1187 +- 0.0018; [PDG2021 says that this corresponds to alpha_s(mtau^2)=0.325±0.016 (not quite true!); my RGE running is different and basically agrees with Pich RGE running] *)
 (* PDG2021, world average: alpha_s(MZ^2)=0.1176±0.0011 (without lattice); 0.1182±0.0008 (lattice); 0.1179±0.0010 (unweighted average of "without lattice" and "lattice" results) *)
 (*-------------------------------------------------*)
NchP=2; (* number of intervals in: (M2min < |M^2} < 1.5) *)
M2FP[ii_]:=M2min+ii*(1.5-M2min)/NchP; (* ii=0,1,2,...,NchP *)
(* We will use sum of squares minimization in this interval, for M^2=|M^2| Exp[I Psi] with Psi=0, Pi/6, Pi/4 one after another in the sum (it turns out that it is not good to include Psi=Pi/3, this value appears to be too large, and xi^2 become large): *)
Nm=NchP+1; Nd=3 Nm;  
M2al0[ii_]:=0.90+(ii-1)*0.60/(Nm-1);
Do[M2al[jj]=M2al0[jj],{jj,1,Nm}]; Do[M2al[jj+Nm]=M2al0[jj] Exp[I Pi/6],{jj,1,Nm}]; Do[M2al[jj+2 Nm]=M2al0[jj] Exp[I Pi/4],{jj,1,Nm}];

 (*---------*)
 (* ka1=0.375; *) ka2=0.778;  (* April 2022 *)
 k3=ka2;
 (* --> 1-ka2=0.222, this is the power in the OPE at the condensate O6; see JHEP 2015 (1510.03812v2) of Boito, Hornung and Jamin, eq. (2.31) *)
 
 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275 (ECH-predicted): *)
 <</home/gorazd/SumRpQCD.dir/MSbarRenMod5Ap3ka1ka2d4275.save;

 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275.+63.=338.: *)
 (* <</home/gorazd/SumRpQCD.dir/MSbarRenMod5Ap3ka1ka2d4338.save; when d4barPREDICTED=338. (RenMod-predicted in MiniMOM): *)
 
 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275.-63.=212.: *)
 (* <</home/gorazd/SumRpQCD.dir/MSbarRenMod5Ap3ka1ka2d4212.save; *)

  (*---------*)
 
 aMSp[Q2f_, fif_, alsinf_] := Module[{Q2ff = Q2f, fiff = fif, alsinff = alsinf}, Clear[alsin] ; alsin = alsinff; <</home/gorazd/SumRpQCD.dir/aMSQ2complS_almtauinput.m; aMSpert[Q2ff, fiff]];
(* aMSp[0.8, 0.4, 0.320] is then a(Q^2=0.8 Exp[I 0.4]; MSbar) with initial condition alpha_s(mtau^2;MSbar)=0.320  (all with nf=3) *)

(* The other way: when we want to fix alpha_s(mtau^2;MSbar)=alsin to a value fixed value, say 0.320: alsin=0.320; <</home/gorazd/SumRpQCD.dir/aMSQ2complS_almtauinput.m; Then MSpert[Q2,phi]=a(Q^2=Q2 Exp[I phi]; MSbar; nf=3) such that Pi a(mtau^2;MSbar)=0.320 *)


mtau=1.77686;  (*= 1.77686+-0.00012 *)

M2min=0.90; (* --> YYY=090 *)  (* M2min=0.65; M2min=1.00; *)
 (* M2min is |M^2|_{min} (in GeV^2) which we will use in the fit procedures; *)
 (* 0.90 =< M2 < 1.5; note that |M^2|_{max}=1.5 always, in our approach *)
(* Note: for wolrd average value alpha_s(MZ^2)=0.1179+-0.0010 -> *)
(* alpha_s(mtau^2)=0.31363+0.00801-0.00772 =appr 0.314 +-0.008 *)
 (* ALEPH: *)
 (* sigmax=2.80 (GeV^2) *)

(*--------------------------------__*)
(* ALEPH experimental data: *)
<</home/gorazd/SumRpQCD.dir/expdataALEPH.m;
(* NNAL=77 will be now, and not NNAL=79 bins *)
NNAL=77;
 sigmax=sbin[[NNAL]] + dsbin[[NNAL]]/2;  (* =2.80 GeV^2 *)
ddel=fPi^2 mPi^2;
fPi = 0.13050;
mPi = 0.13957;
(* gives: OmegaVplA[ii],sbin[[ii]  and dsbin[[ii] (ii=1,...,77) *)
(*--------------------------------------------*)
 (* Now for FOPT: *)

nqu=3; be0=(11 -2 nqu/3)/4; 
be1=(102-38 nqu/3)/16; ce1=be1/be0;
bett2MS=(2857/2 - 5033 nqu/18 + 325 nqu^2/54)/64;
bett3MS=(1/4^4)( (149753/6) + 3564 Zeta[3] - nqu ( (1078361/162) + (6508/27) Zeta[3] )  + nqu^2 ( (50065/162) + (6472/81) Zeta[3] )  + nqu^3 (1093/729) );
bett4MS[Nff_]:=(1/4^5) ( (8157455/16) + (621885/2) Zeta[3] - (88209/2) Zeta[4] - 288090 Zeta[5] + Nff ( -(336460813/1944) - (4811164/81) Zeta[3] + (33935/6) Zeta[4] + (1358995/27) Zeta[5] ) + Nff^2 ( (25960913/1944) + (698531/81) Zeta[3] - (10526/9) Zeta[4] - (381760/81) Zeta[5] ) + Nff^3 (-(630559/5832) - (48722/243) Zeta[3] + (1618/27) Zeta[4] + (460/9) Zeta[5] ) + Nff^4 ( (1205/2916) - (152/81) Zeta[3] ) );
 ce2MS=bett2MS/be0;
ce3MS=bett3MS/be0;
ce4MS=bett4MS[nqu]/be0;

betaMSbar[x_]:=-be0 x^2 - be1 x^3 - bett2MS x^4 - bett3MS x^5 - bett4MS[3.] x^6;
bta[x_]:=betaMSbar[x];
(* Note: x =a(mu^2) = alpha_s(mu^2)/Pi *)
(* x=Log[Q2/sigmax];  *) delext=5.; Clear[x];
solQ2pos[aALf_]:= NDSolve[{y'[x]==bta[y[x]], y[0]==aALf}, y, {x,-0.9,+1.5}];
(* ' *)
(*  aMSext[xx_,aAL0f_? NumberQ]:=First[y[Log[xx]]/.solQ2pos[aAL0f]]; *)
 aMSext[xx_,aAL0f_]:=First[y[Log[xx]]/.solQ2pos[aAL0f]];
			     (* aMSext[x,aAL0]=a(Q^2=x Q20), where Q2=x*Q20 and a(Q^2=Q20) = aAL0 *)

(*----------------------------------------------*)
<</home/gorazd/SumRpQCD.dir/Prgs.dir/SumRthMSbar.save;
  (* This gives coeffext coefficients:  coeffext[jj,nn] is the coefficient of a^nn in ta[jj]={tilde a}_{jj} (nn <= 16) *)
  (*-----------------------------------_*)  
  (* (Double-)Pinched ('P') Borel-Laplace in FO approach (PFO): *)

  kappa=1;  (* kappa = mu^2/Q^2, where mu^2 is RScl *)
  tkap=kappa Exp[tKK];
  
(* Now (double) pinched ('P') Re[BPFOfast]: *)
cKP[kp_,M2p_]:=(1/(2 Pi)) (-I)^kp NIntegrate[ th^kp ( ( (1 + Exp[I th])^2 - 2 (M2p/sigmax) (1 + Exp[I th]) + 2  (M2p/sigmax)^2 ) Exp[(sigmax/M2p) Exp[I th]] - 2 (M2p/sigmax)^2 Exp[-(sigmax/M2p)] ),{th, -Pi,Pi}];
  coeffBPFO[M2f_,Nn_]:=Sum[tddtk[nn,tkap] Binomial[nn + kk, kk] be0^kk coeffext[nn + kk + 1, Nn] cKP[kk,M2f], {nn, 0, Nn - 1}, {kk, 0, Nn - 1 - nn}]; 
 (* coeffBPFO[M2,Nn] is the coefficient at x^Nn in the expansion of BPFOexpD0[x,M2,Nn] or of any BPFOexpD0[x,M2,MM,M2] with MM >= Nn *)
 BPFOD0[M2_,Nn_,aALf_]:=Sum[coeffBPFO[M2,nn] x^nn,{nn,1,Nn}]/.{x->aALf};
(* aALf=alpha_s(kappa*sigmax)/Pi *)
  
(* experimental values of double-pinched Borel-Laplace: *)  
BPexpBinAL[M2_]:=(1/M2) (Sum[ (1 - sbin[[ii]]/sigmax)^2 Exp[-sbin[[ii]]/M2]*OmegaVplA[ii]*dsbin[[ii]],{ii,1,NNAL}] + 2 Pi^2 fPi^2  (1 - mPi^2/sigmax)^2 Exp[-mPi^2/M2]);
ReBPexpBinAL[M2_]:=Re[(1/M2) ( Sum[ (1 - sbin[[ii]]/sigmax)^2 Exp[-sbin[[ii]]/M2]*OmegaVplA[ii]*dsbin[[ii]],{ii,1,NNAL}] + 2 Pi^2 fPi^2 (1 - mPi^2/sigmax)^2 Exp[-mPi^2/M2])];
delReBPexpBinAL[M2_]:=(2 Pi^2) Sqrt[Sum[Re[(1 - sbin[[ii]]/sigmax)^2 Exp[-sbin[[ii]]/M2]/M2] Re[(1 - sbin[[jj]]/sigmax)^2 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}]];
(* end of experimental values *)

(* Covariant-matrix values of the Re(BPexp): *) 
CovReBPexpBinAL[M2A_,M2B_]:=(2 Pi^2)^2 (Sum[Re[(1 - sbin[[ii]]/sigmax)^2 Exp[-sbin[[ii]]/M2A]/M2A] Re[(1 - sbin[[jj]]/sigmax)^2 Exp[-sbin[[jj]]/M2B]/M2B] dsbin[[ii]] dsbin[[jj]] (  covVV[ii,jj]+covAA[ii,jj]+covVA[ii,jj]+ covAV[ii,jj]),{ii,1,NNAL},{jj,1,NNAL}]);
(* CovReBPexpBinAL[M2,M2]=delReBPexpBinAL[M2]^2 *)
  (*-----------*) 
(* Finding of Re[B] in FO approach; *)
(* Now the dfi integration in D=6 term with a(sigmax Exp[I fi])^(1 -ka2) and a(sigmax Exp[I fi])^(1 -ka2) 5-loop CI(RGE)-evolved (' Num5;'):  *)
sol1a[aALf_]:= NDSolve[{-I y'[th] == bta[y[th]], y[0] == aALf}, y, {th, -Pi, +Pi}];
 aCI[th_, aALf_] := First[y[th] /. sol1a[aALf]];
(* aALf=a(kappa sigmax); aCI[th,aALf]=a(Q2=kappa sigma Exp[I th]) ' *)
  
(* Now the dfi integration in D=6 term with 1/a(sigmax Exp[I fi]) 5-loop CI(RGE)-evolved; note that aALf = a(kappa sigmax), and is used in the D=0 part; so a(sigmax)=aMSext[(1/kappa), aALf];  in D=6 we use a(sigmax Exp[I fi]) = aCI[fi, a(sigmax)] = aCI[fi, aMSext[(1/kappa), aALf]] *)
ReBPFONumO85l[M2f_,Nn_,aALf_,O4_,O6_,O8_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/6) (1/M2f^4) (1+ 6 (M2f/sigmax) + 6 (M2f/sigmax)^2) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O62 aCI[fi,aMSext[(1/kappa), aALf]]^(1-ka2) + O61 aCI[fi,aMSext[(1/kappa), aALf]]^(1-ka1)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]]; 

(* Now somewhat simplified ('s') form, which is valid only for kappa=1: *)
ReBPFONumO85ls[M2f_,Nn_,aALf_,O4_,O6_,O8_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/6) (1/M2f^4) (1+ 6 (M2f/sigmax) + 6 (M2f/sigmax)^2) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 aCI[fi,aALf]^(1-ka2)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]]; 

(* Now the form in D=6 with 1-loop running in D=6 ('1l'); when kappa not=1, RGE from aALf=a(kappa sigmax) to a(sigmax) is performed by 5-loop RGE: aALf=a(kappa sigmax); a(sigmax)=aMSext[(1/kappa),aALf]; (note: aALf=a(kappa sigmax) in this general case): *)
ReBPFONumO8[M2f_,Nn_,aALf_,O4_,O6_,O8_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/6) (1/M2f^4) (1+ 6 (M2f/sigmax) + 6 (M2f/sigmax)^2) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aMSext[(1/kappa), aALf]/(1 + I be0 aMSext[(1/kappa), aALf] fi))^(1-ka2)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
			     
(* Now the simplified (faster) form in D=6 and 1-loop running in D=6 ('O8s'), which is valid only for kappa=1: (note: aALf=a(sigmax) in this case)*)
ReBPFONumO8s[M2f_,Nn_,aALf_,O4_,O6_,O8_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/6) (1/M2f^4) (1+ 6 (M2f/sigmax) + 6 (M2f/sigmax)^2)  + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aALf/(1 + I aALf be0 fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]]; 

(* Now the case of FO where we go beyond O8, all with zero anom dimension: *)	
 (* Now the case with simple O6, O8 and O10 (all with zero anom. dimensions): *)
 ReBPFONumO10[M2f_,Nn_,aALf_,O4_,O6_,O8_,O10_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/6) (1/M2f^4) (1+ 6 (M2f/sigmax) + 6 (M2f/sigmax)^2) + 2 Pi^2 O10  (1/24) (1/M2f^5) (1+ 8 (M2f/sigmax) + 12 (M2f/sigmax)^2) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aMSext[(1/kappa), aALf]/(1 + I be0 aMSext[(1/kappa), aALf] fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
(* aALf=a(kappa sigmax); aMSext[(1/kappa), aALf]=a(sigmax) *)
 ReBPFONumO12[M2f_,Nn_,aALf_,O4_,O6_,O8_,O10_,O12_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/3!) (1/M2f^4) (1+ 2*3 (M2f/sigmax) + 3*2 (M2f/sigmax)^2) + 2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2)+ 2 Pi^2 O12  (1/5!) (1/M2f^6) (1+ 2*5 (M2f/sigmax) + 5*4 (M2f/sigmax)^2) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aALf/(1 + I aALf be0 fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
 ReBPFONumO14[M2f_,Nn_,aALf_,O4_,O6_,O8_,O10_,O12_,O14_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/3!) (1/M2f^4) (1+ 2*3 (M2f/sigmax) + 3*2 (M2f/sigmax)^2) + 2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2)+ 2 Pi^2 O12  (1/5!) (1/M2f^6) (1+ 2*5 (M2f/sigmax) + 5*4 (M2f/sigmax)^2)+ 2 Pi^2 O14  (1/6!) (1/M2f^7) (1+ 2*6 (M2f/sigmax) + 6*5 (M2f/sigmax)^2) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aALf/(1 + I aALf be0 fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
 ReBPFONumO16[M2f_,Nn_,aALf_,O4_,O6_,O8_,O10_,O12_,O14_,O16_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/3!) (1/M2f^4) (1+ 2*3 (M2f/sigmax) + 3*2 (M2f/sigmax)^2) + 2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2)+ 2 Pi^2 O12  (1/5!) (1/M2f^6) (1+ 2*5 (M2f/sigmax) + 5*4 (M2f/sigmax)^2)+ 2 Pi^2 O14  (1/6!) (1/M2f^7) (1+ 2*6 (M2f/sigmax) + 6*5 (M2f/sigmax)^2)+ 2 Pi^2 O16  (1/7!) (1/M2f^8) (1+ 2*7 (M2f/sigmax) + 7*6 (M2f/sigmax)^2) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aALf/(1 + I aALf be0 fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];

			     
 ReBPO4[M2f_,O4_]:=Re[2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )];
 ReBPO6[M2f_,O6_]:=Re[(6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aALf/(1 + I aALf be0 fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
 ReBPO8[M2f_,O8_]:=Re[2 Pi^2 O8  (1/6) (1/M2f^4) (1+ 6 (M2f/sigmax) + 6 (M2f/sigmax)^2)];
 ReBPO10[M2f_,O10_]:=Re[2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2)];
 ReBPO12[M2f_,O12_]:=Re[ 2 Pi^2 O12  (1/5!) (1/M2f^6) (1+ 2*5 (M2f/sigmax) + 5*4 (M2f/sigmax)^2)];
 ReBPO14[M2f_,O14_]:=Re[2 Pi^2 O14  (1/6!) (1/M2f^7) (1+ 2*6 (M2f/sigmax) + 6*5 (M2f/sigmax)^2)];			     
(*-------------------*)
(* Now search for condensate values with this Borel-Laplace (double) pinched ('P') approach, with FO method for D=0 part: *)
 
Clear[aALf];
 O4init=-0.003; O6init=-0.001; aALfinit=0.105; 
O8init=0.001; O6init = 0.003; O10init=0.0005; O12init=0.0005; O14init=0.0003; O16init=0.0001;
(* FO approach: *)
 alO8fitPFOMinNum5l[Nn_]:= FindMinimum[Sum[(ReBPFONum5l[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}];
 alO6fitPFOMinNum5l[Nn_]:= FindMinimum[Sum[(ReBPFONum5l[M2al[jj], Nn, aALf,O4fitf,O6fitf,0] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init}}];
			     
 (* Now the form ('1l') with 1-loop RGE running in D=6, valid for general kappa: *)
 alO8fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO8[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}]; 
(* Now O8=0: *) 
 alO6fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO8[M2al[jj], Nn, aALf,O4fitf,O6fitf,0] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init}}];
			     
(* Now additional condensates (D>8) with zero anom. dimension ('an0'): *)
alO10fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO10[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init}}];
alO12fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO12[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init}}];
alO14fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO14[M2al[jj],Nn,aALf,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf,O14fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init},{O14fitf, O14init}}];
alO16fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO16[M2al[jj],Nn,aALf,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf,O14fitf,O16fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init},{O14fitf, O14init},{O16fitf, O16init}}];
			     
 (* Now the simplified form ('1ls') with 1-loop RGE running in D=6, which is valid only for kappa=1: *)
alO8fitPFOMinNum1ls[Nn_]:= FindMinimum[Sum[(ReBPFONum1ls[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}];

  
   ms=mtau^2/sigmax;
    alsmt[asigmax_]:=aMSext[ms,asigmax] Pi;
 (*-----------------------__*)
(* BL results, global FO approach: *)
(* NchP=2 (-> Nd=9;) *) 

(* First OPE up to O6: *)
(* Do[alO6fitPFOMinnum[Nn] = alO6fitPFOMinNum[Nn], {Nn, 3, 12}]: *)
(*
alO6fitPFOMinnum[3]={0.00183444, {aALf -> 0.109719, O4fitf -> 0.00108318, 
  O6fitf -> -0.0000666842}};
alO6fitPFOMinnum[4]={0.013117, {aALf -> 0.104272, O4fitf -> 0.000523923, 
  O6fitf -> 0.000317927}};			     
alO6fitPFOMinnum[5]={0.0548481, {aALf -> 0.102166, O4fitf -> -0.0000320614, 
  O6fitf -> 0.000663257}};
alO6fitPFOMinnum[6]={0.104638, {aALf -> 0.101353, O4fitf -> -0.000596561, 
  O6fitf -> 0.000967138}};
alO6fitPFOMinnum[7]={0.151688, {aALf -> 0.101186, O4fitf -> -0.00120153, 
  O6fitf -> 0.00123956}};
alO6fitPFOMinnum[8]={0.187692, {aALf -> 0.101341, O4fitf -> -0.00187671, 
  O6fitf -> 0.00147986}};
alO6fitPFOMinnum[9]={0.20875, {aALf -> 0.101676, O4fitf -> -0.00267471, 
  O6fitf -> 0.00168898}};
alO6fitPFOMinnum[10]={0.208774, {aALf -> 0.101997, O4fitf -> -0.0036408, 
  O6fitf -> 0.00184885}};
alO6fitPFOMinnum[11]={0.197968, {aALf -> 0.102361, O4fitf -> -0.00495312, 
  O6fitf -> 0.00197854}};
alO6fitPFOMinnum[12]={0.166852, {aALf -> 0.102385, O4fitf -> -0.00663458, 
  O6fitf -> 0.00201692}};
 sm=2.8;
Clear[Nn]; Do[Module[{}, xi2FOO6[Nn]=alO6fitPFOMinnum[Nn][[1]]; aALfFOO6[Nn]=aALf/.alO6fitPFOMinnum[Nn][[2]]; alsmtauFOO6[Nn]=aMSext[mtau^2/sm,aALfFOO6[Nn]] Pi;  O4FOO6[Nn]=O4fitf/.alO6fitPFOMinnum[Nn][[2]]; O6FOO6[Nn]=O6fitf/.alO6fitPFOMinnum[Nn][[2]]],{Nn,3,12}]; (* 3  Oct. 2022; OK *)
 (*------*)	   
(* OPE up to O8: *)			     
(* Do[alO8fitPFOMinnum[Nn] = alO8fitPFOMinNum[Nn], {Nn, 3, 12}]: *)
 alO8fitPFOMinnum[3]={0.0003320765408670988`, {aALf -> 0.10947625819470701`, O4fitf -> 0.0011424824605536655`, O6fitf -> -0.00021115694489708654`, O8fitf -> 0.00007126760422125499`}};
alO8fitPFOMinnum[4]={0.0000437671657334451`, {aALf -> 0.10491004047412288`, O4fitf -> 0.00032971779393080997`, O6fitf -> 0.0007703998633271148`, O8fitf -> -0.00021601526934385918`}};
alO8fitPFOMinnum[5]={0.0006056795956205499`, {aALf -> 0.1034084886515888`, O4fitf -> -0.00047394353176959515`, O6fitf -> 0.0016362077198946686`, O8fitf -> -0.0004528400264098766`}};
alO8fitPFOMinnum[6]={0.0013458631511188827`, {aALf -> 0.10305103906030373`, O4fitf -> -0.00128248576620942`, O6fitf -> 0.0023793000827803837`, O8fitf -> -0.0006423457407170692`}};
alO8fitPFOMinnum[7]={0.001966420198896209`, {aALf -> 0.10326279325376696`, O4fitf -> -0.0021426733795407373`, O6fitf -> 0.003024818494579634`, O8fitf -> -0.0007938140996529453`}};
alO8fitPFOMinnum[8]={0.0022829918575089287`, {aALf -> 0.1037264221189971`, O4fitf -> -0.0030897230158157766`, O6fitf -> 0.0035590553739760244`, O8fitf -> -0.0009034138999636468`}};
alO8fitPFOMinnum[9]={0.002266948631210346`, {aALf -> 0.10430084852705807`, O4fitf -> -0.004191499978876304`, O6fitf -> 0.003973551717158389`, O8fitf -> -0.0009694176771084516`}};
alO8fitPFOMinnum[10]={0.0018919440718063092`, {aALf -> 0.10471293449912358`, O4fitf -> -0.0054673208772661125`, O6fitf -> 0.004192309505371942`, O8fitf -> -0.0009736714722481764`}};
alO8fitPFOMinnum[11]={0.0015073584473872487`, {aALf -> 0.10509782201559896`, O4fitf -> -0.00718713212121124`, O6fitf -> 0.00429579956707237`, O8fitf -> -0.0009436177514408772`}};
alO8fitPFOMinnum[12]={0.000999633729238889`, {aALf -> 0.10482784259342191`, O4fitf -> -0.009160146189470817`, O6fitf -> 0.004091296272991851`, O8fitf -> -0.0008410482908200126`}};
  sm=2.8;
  Clear[Nn]; Do[Module[{}, xi2FOO8[Nn]=alO8fitPFOMinnum[Nn][[1]]; aALfFOO8[Nn]=aALf/.alO8fitPFOMinnum[Nn][[2]]; alsmtauFOO8[Nn]=aMSext[mtau^2/sm,aALfFOO8[Nn]] Pi;  O4FOO8[Nn]=O4fitf/.alO8fitPFOMinnum[Nn][[2]]; O6FOO8[Nn]=O6fitf/.alO8fitPFOMinnum[Nn][[2]]; O8FOO8[Nn]=O8fitf/.alO8fitPFOMinnum[Nn][[2]]],{Nn,3,12}]; (* 19 Sept. 2022; OK *)
 (*-----*)  
  (* When including O10: *)
alO10fitPFOMinnum[3]={6.11766569371541`*^-6, {aALf -> 0.10873977272220776`, O4fitf -> 0.0013589500250509796`, O6fitf -> -0.0008828923876971016`, O8fitf -> 0.0005955917442750289`, O10fitf -> -0.0003281795611166152`}};
alO10fitPFOMinnum[4]={3.31314*10^-7, {aALf -> 0.105155, O4fitf -> 0.000240625, O6fitf -> 0.00103396, O8fitf -> -0.000416054, O10fitf -> 0.00012374}};
alO10fitPFOMinnum[5]={7.48823*10^-6, {aALf -> 0.104293, O4fitf -> -0.00084901, O6fitf -> 0.00268136, O8fitf -> -0.00122807, O10fitf -> 0.000473255}};
alO10fitPFOMinnum[6]={0.000016857910673044713`, {aALf -> 0.10437876547460796`, O4fitf -> -0.0019204515004994278`, O6fitf -> 0.004034835903445157`, O8fitf -> -0.0018438441890750155`, O10fitf -> 0.000723842369965663`}};
alO10fitPFOMinnum[7]={0.000023814383921098708`, {aALf -> 0.10491855920446051`, O4fitf -> -0.003034067860368765`, O6fitf -> 0.005143751252951458`, O8fitf -> -0.002297489698624533`, O10fitf -> 0.0008940258827032286`}};
alO10fitPFOMinnum[8]={0.000025670889484749748`, {aALf -> 0.1055909512932032`, O4fitf -> -0.00421526293750236`, O6fitf -> 0.005957039292628349`, O8fitf -> -0.002565953796231432`, O10fitf -> 0.0009767589344250825`}};
alO10fitPFOMinnum[9]={0.000022796191350809326`, {aALf -> 0.10625375384599656`, O4fitf -> -0.005530199678813183`, O6fitf -> 0.006451771734039753`, O8fitf -> -0.0026472468460906054`, O10fitf -> 0.000976554286208314`}};
alO10fitPFOMinnum[10]={0.00001602593796943203`, {aALf -> 0.10654148509571322`, O4fitf -> -0.00692496382585033`, O6fitf -> 0.006468559523793831`, O8fitf -> -0.002486205678118484`, O10fitf -> 0.0008787394478097528`}};
alO10fitPFOMinnum[11]={0.000011359986980929269`, {aALf -> 0.10676493306572749`, O4fitf -> -0.008801443939441496`, O6fitf -> 0.006309868437339813`, O8fitf -> -0.0022611394051049007`, O10fitf -> 0.0007679226612097049`}};
alO10fitPFOMinnum[12]={6.580933702780576`*^-6, {aALf -> 0.10609467180884717`, O4fitf -> -0.01070301031873825`, O6fitf -> 0.0056338982885779645`, O8fitf -> -0.0018573874137424447`, O10fitf -> 0.0006030370265580529`}};
  sm=2.8; Clear[Nn]; Do[Module[{}, xi2FOO10[Nn]=alO10fitPFOMinnum[Nn][[1]]; aALfFOO10[Nn]=aALf/.alO10fitPFOMinnum[Nn][[2]]; alsmtauFOO10[Nn]=aMSext[mtau^2/sm,aALfFOO10[Nn]] Pi;  O4FOO10[Nn]=O4fitf/.alO10fitPFOMinnum[Nn][[2]]; O6FOO10[Nn]=O6fitf/.alO10fitPFOMinnum[Nn][[2]]; O8FOO10[Nn]=O8fitf/.alO10fitPFOMinnum[Nn][[2]]]; O10FOO10[Nn]=O10fitf/.alO10fitPFOMinnum[Nn][[2]],{Nn,3,12}]; (* 19 Sept. 2022, OK *)
  (*-----*)  
 (* When including O12: *)
alO12fitPFOMinnum[3]={5.960320942215042`*^-8, {aALf -> 0.10834244442110494`, O4fitf -> 0.001505742131129495`, O6fitf -> -0.0014675083507159326`, O8fitf -> 0.001237685789067114`, O10fitf -> -0.0011117936583177888`, O12fitf -> 0.0005558112547609793`}};
alO12fitPFOMinnum[4]={7.91855576832514`*^-9, {aALf -> 0.1052392489437132`, O4fitf -> 0.0002026273918172919`, O6fitf -> 0.0011777554572750318`, O8fitf -> -0.000569869565344403`, O10fitf -> 0.0003091533572276277`, O12fitf -> -0.00013050564913669353`}};
alO12fitPFOMinnum[5]={1.009251868258666`*^-7, {aALf -> 0.10468953371771345`, O4fitf -> -0.0010539555369480473`, O6fitf -> 0.003409599606934195`, O8fitf -> -0.0019897983675667215`, O10fitf -> 0.0013793036490508585`, O12fitf -> -0.0006329771389260049`}};
alO12fitPFOMinnum[6]={2.055472978110436`*^-7, {aALf -> 0.10498298081226208`, O4fitf -> -0.002269969485574706`, O6fitf -> 0.0051877836581429625`, O8fitf -> -0.003024193483471328`, O10fitf -> 0.002109773719828882`, O12fitf -> -0.000961766157262958`}};
alO12fitPFOMinnum[7]={2.7082562183851253`*^-7, {aALf -> 0.10566568847989169`, O4fitf -> -0.0035123346566881194`, O6fitf -> 0.0065838122811898185`, O8fitf -> -0.003739667578518527`, O10fitf -> 0.0025665421736791177`, O12fitf -> -0.0011539285921355422`}};
alO12fitPFOMinnum[8]={2.709027887618087`*^-7, {aALf -> 0.10640666986066298`, O4fitf -> -0.0047937389704825325`, O6fitf -> 0.007513177909546774`, O8fitf -> -0.004090640156615137`, O10fitf -> 0.0027272096701284614`, O12fitf -> -0.0012031121097324442`}};
alO12fitPFOMinnum[9]={2.2139585347264137`*^-7, {aALf -> 0.10706341597815373`, O4fitf -> -0.0061744518997120515`, O6fitf -> 0.007954975726548682`, O8fitf -> -0.004090433364888927`, O10fitf -> 0.0026235519019683724`, O12fitf -> -0.0011310882135660222`}};
alO12fitPFOMinnum[10]={1.4171101429339786`*^-7, {aALf -> 0.10723044461611325`, O4fitf -> -0.007554883943158532`, O6fitf -> 0.007716945486272621`, O8fitf -> -0.003671931351807432`, O10fitf -> 0.0022368938628504692`, O12fitf -> -0.0009370818358878611`}};
alO12fitPFOMinnum[11]={1.0005675435782689`*^-7, {aALf -> 0.10735175487354545`, O4fitf -> -0.009445572900498236`, O6fitf -> 0.007340356827400047`, O8fitf -> -0.003234239729927077`, O10fitf -> 0.0018913465298946047`, O12fitf -> -0.0007798882889302377`}};
alO12fitPFOMinnum[12]={6.00365851794509`*^-8, {aALf -> 0.10650224007088809`, O4fitf -> -0.011257984109015011`, O6fitf -> 0.0063703931300022765`, O8fitf -> -0.00256599607431514`, O10fitf -> 0.0014371137768128743`, O12fitf -> -0.0005847366124712279`}};
  sm=2.8;  Clear[Nn]; Do[Module[{}, xi2FOO12[Nn]=alO12fitPFOMinnum[Nn][[1]]; aALfFOO12[Nn]=aALf/.alO12fitPFOMinnum[Nn][[2]]; alsmtauFOO12[Nn]=aMSext[mtau^2/sm,aALfFOO12[Nn]] Pi;  O4FOO12[Nn]=O4fitf/.alO12fitPFOMinnum[Nn][[2]]; O6FOO12[Nn]=O6fitf/.alO12fitPFOMinnum[Nn][[2]]; O8FOO12[Nn]=O8fitf/.alO12fitPFOMinnum[Nn][[2]]]; O10FOO12[Nn]=O10fitf/.alO12fitPFOMinnum[Nn][[2]]; O12FOO12[Nn]=O12fitf/.alO12fitPFOMinnum[Nn][[2]],{Nn,3,12}]; (* 19 Sept. 2022; OK *)
  (*-------*)
  (* When including O14 (Nt=5-8): *)
alO14fitPFOMinnum[5]={4.416554616332084`*^-11, {aALf -> 0.10505247024628854`, O4fitf -> -0.0012534825518710159`, O6fitf -> 0.004176669019744087`, O8fitf -> -0.0029003727823937514`, O10fitf -> 0.0027603300901433785`, O12fitf -> -0.002263957863509443`, O14fitf -> 0.0011090527133476064`}};
alO14fitPFOMinnum[6]={7.670209762355503`*^-11, {aALf -> 0.10551259050690996`, O4fitf -> -0.0025952413039984877`, O6fitf -> 0.006346060193002115`, O8fitf -> -0.004369654457840858`, O10fitf -> 0.004122919090738741`, O12fitf -> -0.0033222153579891144`, O14fitf -> 0.001597886875541521`}};
alO14fitPFOMinnum[7]={9.236039221995693`*^-11, {aALf -> 0.10630180883165337`, O4fitf -> -0.003943959483179591`, O6fitf -> 0.007982995857956851`, O8fitf -> -0.005329000249234984`, O10fitf -> 0.004914574696272011`, O12fitf -> -0.003890486146477185`, O14fitf -> 0.0018460755980352506`}};
alO14fitPFOMinnum[8]={8.589511358215219`*^-11, {aALf -> 0.10707754283713948`, O4fitf -> -0.005296942127302009`, O6fitf -> 0.0089665953445103`, O8fitf -> -0.005705787867099653`, O10fitf -> 0.005090304310580939`, O12fitf -> -0.003947661723053665`, O14fitf -> 0.0018487004776619588`}};
alO14fitPFOMinnum[9]={6.631259247424694`*^-11, {aALf -> 0.10770102443675497`, O4fitf -> -0.006709540069331502`, O6fitf -> 0.009292863759926582`, O8fitf -> -0.005548592031185111`, O10fitf -> 0.004747113050889547`, O12fitf -> -0.0035979213878794126`, O14fitf -> 0.001663433218844904`}};
 sm=2.8;  Clear[Nn]; Do[Module[{}, xi2FOO14[Nn]=alO14fitPFOMinnum[Nn][[1]]; aALfFOO14[Nn]=aALf/.alO14fitPFOMinnum[Nn][[2]]; alsmtauFOO14[Nn]=aMSext[mtau^2/sm,aALfFOO14[Nn]] Pi;  O4FOO14[Nn]=O4fitf/.alO14fitPFOMinnum[Nn][[2]]; O6FOO14[Nn]=O6fitf/.alO14fitPFOMinnum[Nn][[2]]; O8FOO14[Nn]=O8fitf/.alO14fitPFOMinnum[Nn][[2]]]; O10FOO14[Nn]=O10fitf/.alO14fitPFOMinnum[Nn][[2]]; O12FOO14[Nn]=O12fitf/.alO14fitPFOMinnum[Nn][[2]]; O14FOO14[Nn]=O14fitf/.alO14fitPFOMinnum[Nn][[2]],{Nn,5,9}];
  (* 19 Sept. 2022 *)
 (*-------*)
  (* When including O16 (Nt=5-8): *)
alO16fitPFOMinnum[5]={8.803371462624767`*^-12, {aALf -> 0.1054045252381648`, O4fitf -> -0.0014627825814064262`, O6fitf -> 0.0050566824803591455`, O8fitf -> -0.004078929482066488`, O10fitf -> 0.00492295191307302`, O12fitf -> -0.005723100754318357`, O14fitf -> 0.005304547100085862`, O16fitf -> -0.00287054653674426`}};
alO16fitPFOMinnum[6]={1.4430444405032987`*^-11, {aALf -> 0.10599313591599337`, O4fitf -> -0.0029128372740897424`, O6fitf -> 0.007581048178949539`, O8fitf -> -0.005988216747253488`, O10fitf -> 0.007053164716447614`, O12fitf -> -0.007976632791047626`, O14fitf -> 0.007219984162352174`, O16fitf -> -0.00383688030130917`}};
alO16fitPFOMinnum[7]={1.6824827466513236`*^-11, {aALf -> 0.10685839641086918`, O4fitf -> -0.004348350976421998`, O6fitf -> 0.009412220240140148`, O8fitf -> -0.007160985941115983`, O10fitf -> 0.008190526499392323`, O12fitf -> -0.009064679172050082`, O14fitf -> 0.008076957904330825`, O16fitf -> -0.004244851295592659`}};
alO16fitPFOMinnum[8]={1.5231998944619104`*^-11, {aALf -> 0.10764592674795546`, O4fitf -> -0.005751106241713302`, O6fitf -> 0.0103955801953494`, O8fitf -> -0.007499260837870933`, O10fitf -> 0.00826957074549516`, O12fitf -> -0.008955100881719089`, O14fitf -> 0.007872630681522098`, O16fitf -> -0.004102379035386993`}};
alO16fitPFOMinnum[9]={1.1531004224272422`*^-11, {aALf -> 0.10822571830152997`, O4fitf -> -0.00717606894087501`, O6fitf -> 0.010564834381426087`, O8fitf -> -0.0071175621645237835`, O10fitf -> 0.007520547484351395`, O12fitf -> -0.007971146277187306`, O14fitf -> 0.006932480573786434`, O16fitf -> -0.0035929196267454533`}}; 
 sm=2.8;  Clear[Nn]; Do[Module[{}, xi2FOO16[Nn]=alO16fitPFOMinnum[Nn][[1]]; aALfFOO16[Nn]=aALf/.alO16fitPFOMinnum[Nn][[2]]; alsmtauFOO16[Nn]=aMSext[mtau^2/sm,aALfFOO16[Nn]] Pi;  O4FOO16[Nn]=O4fitf/.alO16fitPFOMinnum[Nn][[2]]; O6FOO16[Nn]=O6fitf/.alO16fitPFOMinnum[Nn][[2]]; O8FOO16[Nn]=O8fitf/.alO16fitPFOMinnum[Nn][[2]]]; O10FOO16[Nn]=O10fitf/.alO16fitPFOMinnum[Nn][[2]]; O12FOO16[Nn]=O12fitf/.alO16fitPFOMinnum[Nn][[2]]; O14FOO16[Nn]=O14fitf/.alO16fitPFOMinnum[Nn][[2]]; O16FOO16[Nn]=O16fitf/.alO16fitPFOMinnum[Nn][[2]], {Nn,5,9}];
 (* 9 Oct. 2022 *)
 *)
  (* The extracted values of alpha_s do not change much when O12 or O12&O14 are included in the OPE; this is GOOD! *)			     
(*-----*)  
(* The convergence of alpha_s(mtau^2) when we increase the number of OPE terms: *) 			     
(* alsmtauFOO8[8]=0.314786; alsmtauFOO10[8]=0.320192; alsmtauFOO12[8]=0.322553; alsmtauFOO14[8]=0.324493;  alsmtauFOO16[8]=0.326136; convergence is slow starting from O10 truncation on; I will choose O10 truncation, because from O8 to O10 no sign of convergence; O10->O12 and O12->O14 and O14 -> O16 transition show almost equally reasonable convergenvce *)
     
(* end of BL results, global FO approach with 'p3ka1ka2' & 'O6ka2'; 19 Sept. 2022) *)
  (*----------------------------------------------------*)
(*---------------------------------------------------------------------_*)
 (* Now FESR (2,n) moments a^{(2,n)} (i.e., double-pinched) in FO: *)

 (*--------*)
 (* First the experimental ALEPH values of (2,n) moments a^{(2,n)} (see also file rtauALEPH.m): *)
grs0a[2,nn_][sig_]:=((nn+3)/(nn+1)) (1/sigmax) (1 - (nn+2) (sig/sigmax)^(nn+1) + (nn+1) (sig/sigmax)^(nn+2) );
aexpBin77[2,nn_]:=(Sum[ dsbin[[ii]] grs0a[2,nn][sbin[[ii]]] OmegaVplA[ii],{ii,1,77}] + 2 Pi^2 fPi^2 grs0a[2,nn][mPi^2]) -  1;
delaexpBin77[2,nn_]:=(2 Pi^2) Sqrt[Sum[ grs0a[2,nn][sbin[[ii]]] grs0a[2,nn][sbin[[jj]]] dsbin[[ii]] dsbin[[jj]](  covVV[ii,jj]+covAA[ii,jj]+ covVA[ii,jj]+covAV[ii,jj] ),{ii,1,77},{jj,1,77}]];
(*  aexpBin77[2,n] +- delaexpBin77[2,n] = a^{(2,n)}_{exp} *)
(*--------*)
kappa=1.0;
 tkap=kappa*Exp[tKK];
 (*-----*)
(* moments a^{(2,n}}, FO: *)			     
 cAg[2,nn_][kp_]:=(1/(2 Pi)) ((nn+3)/(nn+1)) Re[ (-I)^kp NIntegrate[ th^kp ( Exp[I th] (1 - (-Exp[I th])^(nn+1)) + ((nn+1)/(nn+3))(1 - (-Exp[I th])^(nn+3)) ),{th, -Pi,Pi}] ];  
(* kp = 0, 1, 2, ... *)
 coeffaaFO[2,n_][NN_]:=Sum[tddtk[nn,tkap] Binomial[nn + kk, kk] be0^kk coeffext[nn + kk + 1, NN] cAg[2,n][kk], {nn, 0, NN - 1}, {kk, 0, NN - 1 - nn}];
 coeffaaFOIR2[2,n_][NN_]:=Sum[tddtkIR2[nn,tkap] Binomial[nn + kk, kk] be0^kk coeffext[nn + kk + 1, NN] cAg[2,n][kk], {nn, 0, NN - 1}, {kk, 0, NN - 1 - nn}];
  coeffaaFOIR3[2,n_][NN_]:=Sum[tddtkIR3[nn,tkap] Binomial[nn + kk, kk] be0^kk coeffext[nn + kk + 1, NN] cAg[2,n][kk], {nn, 0, NN - 1}, {kk, 0, NN - 1 - nn}];
  coeffaaFOUV1[2,n_][NN_]:=Sum[tddtkUV1[nn,tkap] Binomial[nn + kk, kk] be0^kk coeffext[nn + kk + 1, NN] cAg[2,n][kk], {nn, 0, NN - 1}, {kk, 0, NN - 1 - nn}];
(* coeffaaFO[2,n][NN] is the coefficient at x^NN in the expansion of aa[2,n] *)
aaFOD0[2,n_,aALf_][NN_]:=Sum[coeffaaFO[2,n][nn] x^nn,{nn,1,NN}]/.{x->aALf};
(* aALf=a(kappa*sigmax) *)
(* aaFOD0[2,n,aALf][NN]: this is the D=0 contribution to a^{(2,n)}(sigmax), as a function of aALf=a(sigmax), in the FO approach with truncation index Nt=NN *)

(* Now including D >= 4 contributions to a^{(2,n)}(sigmax): *)

(* First a^{(2,0)}(sigmax)  (Nn is the truncation index Nt of the D=0 part): *)
 aaFO20[aALf_][Nn_, O4f_,O6f_]:=aaFOD0[2,0,aALf][Nn]+ ((0+3)/(0+1)) 2 Pi^2 (-1)^0 ( (0+2) O4f/sigmax^(0+2) ) + 3 (1/2/Pi) (6 Pi^2/sigmax^3) ( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[I 3 fi]) ) (O6f (aALf/(1 + I aALf be0 fi))^(1-ka2)),{fi,-Pi,Pi}] );
(* a^{(2,1)}(sigmax): *)
aaFO21[aALf_][Nn_, O6f_,O8f_]:=aaFOD0[2,1,aALf][Nn]+ 2 (1/2/Pi) (6 Pi^2/sigmax^3)*( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 - Exp[2 I fi]) + (1/2) (1 - Exp[I 4 fi]) ) (O6f (aALf/(1 + I aALf be0 fi))^(1-ka2)),{fi,-Pi,Pi}] ) + ((1+3)/(1+1)) 2 Pi^2 (-1)^1 ( (1+1) O8f/sigmax^(1+3) );			     
(* a^{(2,2)}(sigmax): *)
aaFO22[aALf_][Nn_, O6f_,O8f_,O10f_]:=aaFOD0[2,2,aALf][Nn]+ (6 Pi^2 O6f/sigmax^3)*(1/2/Pi)*( NIntegrate[ Exp[-I 3 fi] ((aALf/(1 + I aALf be0 fi))^(1-ka2)) ( ((2+3)/(2+1)) Exp[I fi] (1 +(-1)^2 Exp[(2+1) I fi]) +  (1 + (-1)^2 Exp[I (2+3) fi]) ),{fi,-Pi,Pi}] ) + ((2+3)/(2+1)) 2 Pi^2 (-1)^2 ( (2+2) O8f/sigmax^(2+2) + (2+1) O10f/sigmax^(2+3));
(* a^{(2,3)}(sigmax): *)
aaFO23[aALf_][Nn_, O6f_,O10f_,O12f_]:=aaFOD0[2,3,aALf][Nn]+ (6 Pi^2 O6f/sigmax^3)*(1/2/Pi)*( NIntegrate[ Exp[-I 3 fi] ((aALf/(1 + I aALf be0 fi))^(1-ka2)) ( ((3+3)/(3+1)) Exp[I fi] (1 +(-1)^3 Exp[(3+1) I fi]) +  (1 + (-1)^3 Exp[I (3+3) fi]) ),{fi,-Pi,Pi}] ) + ((3+3)/(3+1)) 2 Pi^2 (-1)^3 ( (3+2) O10f/sigmax^(3+2) + (3+1) O12f/sigmax^(3+3));
(* a^{(2,4)}(sigmax): *)
aaFO24[aALf_][Nn_, O6f_,O12f_,O14f_]:=aaFOD0[2,4,aALf][Nn]+ (6 Pi^2 O6f/sigmax^3)*(1/2/Pi)*( NIntegrate[ Exp[-I 3 fi] ((aALf/(1 + I aALf be0 fi))^(1-ka2)) ( ((4+3)/(4+1)) Exp[I fi] (1 +(-1)^4 Exp[(4+1) I fi]) +  (1 + (-1)^4 Exp[I (4+3) fi]) ),{fi,-Pi,Pi}] ) + ((4+3)/(4+1)) 2 Pi^2 (-1)^4 ( (4+2) O12f/sigmax^(4+2) + (4+1) O14f/sigmax^(4+3));
(*----------------------------------------------------*)
(* For the 'O10' case (i.e., our MAIN considered case), finding optimal Nt: *)

 (*   
Tab20D0FOO10avar=Table[{NN,Re[aaFOD0[2,0,aALfFOO10[NN]][NN]]},{NN,3,12}]; 
Tab21D0FOO10avar=Table[{NN,Re[aaFOD0[2,1,aALfFOO10[NN]][NN]]},{NN,3,12}]; 
Tab22D0FOO10avar=Table[{NN,Re[aaFOD0[2,2,aALfFOO10[NN]][NN]]},{NN,3,12}]; *)
 (* Tab20FOO10avar=Table[{NN,Re[aaFO20[aALfFOO10[NN]][NN,O4FOO10[NN],O6FOO10[NN]]]},{NN,3,12}];
Tab21FOO10avar=Table[{NN,Re[aaFO21[aALfFOO10[NN]][NN,O6FOO10[NN],O8FOO10[NN]]]},{NN,3,12}]; 
 Tab22FOO10avar=Table[{NN,Re[aaFO22[aALfFOO10[NN]][NN,O6FOO10[NN],O8FOO10[NN],O10FOO10[NN]]]},{NN,3,12}];
 *)			     

(*
 Tab20FOO10avar(*={{3, 0.240443}, {4, 0.240364}, {5, 0.240307}, {6, 0.240274}, {7, 0.240256}, {8, 0.240253}, {9, 0.240262}, {10, 0.240282}, {11, 0.240298}, {12, 0.240317}} *)={{3, 0.2404433656481368`}, {4, 0.2403641748034315`}, {5, 0.24030725083709822`}, {6, 0.24027390115792013`}, {7, 0.24025567153593705`}, {8, 0.24025278052724663`}, {9, 0.24026187978013136`}, {10, 0.24028249062697601`}, {11, 0.24029789858441927`}, {12, 0.24031718778148875`}}; (*  minimum {8, 0.240253} at Nt=8; experimental value 0.2404+-0.0026 reproduced to high accuracy; see also: delTab20FOO10avar= Table[{Nn, Tab20FOO10avar[[Nn - 2]][[2]] - Tab20FOO10avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; -this shows that the minimal local sensitivity is achieved at Nt=8 *)
 Tab21FOO10avar(* ={{3, 0.206288}, {4, 0.206039}, {5, 0.205865}, {6, 0.20576}, {7, 0.205703}, {8, 0.205695}, {9, 0.205723}, {10, 0.205787}, {11, 0.205835}, {12, 0.205895}} *)={{3, 0.2062878139081978`}, {4, 0.20603880734016033`}, {5, 0.2058648015126896`}, {6, 0.20575961782811167`}, {7, 0.20570327952006257`}, {8, 0.20569463188895185`}, {9, 0.20572310689628914`}, {10, 0.20578714316146599`}, {11, 0.2058347999449503`}, {12, 0.2058945542254531`}}; (*  minimum {8, 0.205695} at Nt=8; experimental value 0.2061+-0.0026 reproduced to high accuracy; see also: delTab21FOO10avar= Table[{Nn, Tab21FOO10avar[[Nn - 2]][[2]] - Tab21FOO10avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; -this shows that the minimal local sensitivity is achieved at Nt=8 *)
 Tab22FOO10avar(*={{3, 0.196504}, {4, 0.195886}, {5, 0.195458}, {6, 0.195199}, {7, 0.195062}, {8, 0.195042}, {9, 0.195113}, {10, 0.19527}, {11, 0.195387}, {12, 0.195533}} *)={{3, 0.19650383805580934`}, {4, 0.19588600132780692`}, {5, 0.19545823138178284`}, {6, 0.19519900929312317`}, {7, 0.19506154529914765`}, {8, 0.19504162914128295`}, {9, 0.19511266151635243`}, {10, 0.19527048631550425`}, {11, 0.1953869777584099`}, {12, 0.19553326145037297`}}; (* minimum {8, 0.195042} at Nt=8; experimental value 0.1960+-0.0028 reproduced reasonably well (minimum at Nt=8 by 0.36 sigma lower); see also: delTab22FOO10avar= Table[{Nn, Tab22FOO10avar[[Nn - 2]][[2]] - Tab22FOO10avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; -this shows that the minimal local sensitivity is achieved at Nt=8 *)
 *)
 (* Note: aexpBin77[2, 0] + -delaexpBin77[2, 0] = 0.24038 + -0.00261;
          aexpBin77[2, 1] + -delaexpBin77[2, 1] = 0.20610 + -0.00260; 
          aexpBin77[2, 2] + -delaexpBin77[2, 2] = 0.19603 + -0.00282; *)

 (* Therefore, in our FO O10 case, we will choose the truncation index N0=Nt=8 (although Nt=7 would also be almost equally acceptable); note: alsmtauFOO10[8]=0.320192; while alsmtauFOO10[7]=0.318243, quite close; 19 Sept. 2022 *)
(*-----------------------------------------------------------*)
 (* For the 'O12' case, the optimal Nt can be Nt=7 or Nt=8; Nt=7 slightly preferred in this case; note: alsmtauFOO12[7]=0.320408 (very close to Nt=8 case of O10 approach; good!); alsmtauFOO12[8]=0.322553; 19 Sept. 2022 *)
(*-----------------------------------------------------------*)
 (*----------------------------------------------------*)
 (* I WILL TAKE O10 OPE APPROACH (with Nt=8); this because the momenta a20 and a21 are reproduced to high accuracy; and also we have a reasonably good convergence of alpha_s when the OPE truncation order increases to beyond D=10: alsmtauFOO8[8]=0.314786; alsmtauFOO10[8]=0.320192; alsmtauFOO12[8]=0.322553 (alsmtauFOO12[7]=0.320408); alsmtauFOO14[8]=0.324493; 19 Sept. 2022 *)
 (*-------------------------------------------*)
  (*-------------------------------------------*)
 (*-------------------------------------------*)
 (* Now the estimate of EXPERIMENTAL uncertainty in this global approach with double-pinched FO Borel-Laplace (with Nt=N0(=8)) approach 'O10': *)
 (* See Appendix  PRD84, 113006 (Boito et al., 2011): *)
 
(* Now the M^2_{al}: all 9 values (three rays in the complex M^2-plane) as one array, as M2al[jj] (jj=1,...,9): *)
 (* Note, we have NchP=2; Nm=NchP+1; Nd=3 Nm (=9); *)

 (* alO10fitPFOMinnum[8]={0.000025670889484749748`, {aALf -> 0.1055909512932032`, O4fitf -> -0.00421526293750236`, O6fitf -> 0.005957039292628349`, O8fitf -> -0.002565953796231432`, O10fitf -> 0.0009767589344250825`}};  alsmtauFOO10[8]=0.32019158726284414`; *)
 N0=8; epp1=0.005; Np=5;
(* First derivatives of ReBPFONNumO10[M2] with respect to p1=aAL, p2=O4, p3=O6, p4=O8, p5=O10: *)
 dReBPFOd[1][M2_]:=(ReBPFONumO10[M2, N0, aALfFOO10[N0](1+epp1),0,O6FOO10[N0],0,0]-ReBPFONumO10[M2, N0, aALfFOO10[N0](1-epp1),0,O6FOO10[N0],0,0])/(2 epp1 aALfFOO10[N0]);  
dReBPFOd[2][M2_]:=2 Pi^2 Re[ (1/M2^2) ( 1 + 2 M2/sigmax )];
dReBPFOd[3][M2_]:=(6 Pi^2/sigmax^3) (1/2/Pi) Re[NIntegrate[ ((aALfFOO10[N0]/(1 + I be0 aALfFOO10[N0] fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2/sigmax) (1 + Exp[I fi]) + 2 (M2/sigmax)^2) Exp[(sigmax/M2) Exp[I fi]] - 2 (M2/sigmax)^2 Exp[- sigmax/M2] ),{fi, -Pi, Pi}]]; 
dReBPFOd[4][M2_]:=2 Pi^2  (1/3!) Re[(1/M2^4) (1+ 6 (M2/sigmax) + 6 (M2/sigmax)^2)]; 
dReBPFOd[5][M2_]:=2 Pi^2 (1/4!) Re[(1/M2^5) (1+ 8 (M2/sigmax) + 12 (M2/sigmax)^2)];
(* Now Tables of derivatives, for all Mal[ii]: *) 
 tabdReBPFOd[1] = Table[{M2al[ii], dReBPFOd[1][M2al[ii]]}, {ii, 1, Nd}]; (* This one takes some time to evaluate *)
 tabdReBPFOd[2] = Table[{M2al[ii], dReBPFOd[2][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[3] = Table[{M2al[ii], dReBPFOd[3][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[4] = Table[{M2al[ii], dReBPFOd[4][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[5] = Table[{M2al[ii], dReBPFOd[5][M2al[ii]]}, {ii, 1, Nd}];

(* The Table tabCovReBPexpBinAL requires longer time, so better evaluate it only once, namely here, and then use it for all the methods: *)
 tabCovReBPexpBinAL = Table[{M2al[ii],M2al[jj],CovReBPexpBinAL[M2al[ii],M2al[jj]]},{ii,1,Nd},{jj,1,Nd}];

 Np=5; (* This is the number of parameters that we use, in our OPE O10 approach (p1=aAL, p2=O4, p3=O6, p4=O8, p5=O10) *)
 (* Now the NpxNp mtrx A: *)
 (* Now calculate the NpxNp mtrx AA(=A), for the FO case, using tables obtained above: *)
 tabAAFO=Table[{ik, jk, Sum[(tabdReBPFOd[ik][[jj]][[2]] tabdReBPFOd[jk][[jj]][[2]])/tabCovReBPexpBinAL[[jj]][[jj]][[3]], {jj, 1, Nd}]}, {ik, 1, Np}, {jk, 1, Np}];
 Do[AAFOa[ii,jj]=tabAAFO[[ii]][[jj]][[3]],{ii,1,Np},{jj,1,Np}];  (* This is NpxNp matrix A (or: AA), calculated more efficiently *)
 AAFOm=Table[AAFOa[ik,jk],{ik,1,Np},{jk,1,Np}];
(* AAFOm[[ik,jk]] = AA[ik,jk]; AAFOm is the NpxNp matrix *)
AAFOmI=Inverse[AAFOm];
 (* Now the NpxNp covariance matrix of the fitted parameters pj: *) 
 (* p1=aAL, p2=<O4>, p3=<O6>, p4=<O8>, p5=<O10>; where O_{2n} means: O_{2n}=2 <O_{2n}>_V  *)
(* PP[i,j]= < (delta p_i) (delta p_j) >: *)
  (* We calulate PP[i,j] covariance matrix of the extracted parameters, using the tables obtained above: *) 
 PPtabFO[jp_,kp_]:=Sum[AAFOmI[[jp,jpp]] AAFOmI[[kp,kpp]] tabdReBPFOd[jpp][[Al]][[2]] tabdReBPFOd[kpp][[Be]][[2]] tabCovReBPexpBinAL[[Al]][[Be]][[3]]/(tabCovReBPexpBinAL[[Al]][[Al]][[3]]  tabCovReBPexpBinAL[[Be]][[Be]][[3]]),{jpp,1,Np},{kpp,1,Np},{Be,1,Nd},{Al,1,Nd}];
 (* epp1=0.005; Table[{Nn, Sqrt[PPtabFO[Nn, Nn]]}, {Nn, 1, Np}]={{1, 0.000975657}, {2, 0.000476129}, {3, 0.000766927}, {4, 0.000429174}, {5, 0.000201231}}; *)
 (* The experimental uncertainty of aAL=a(sigmax) is: del(aAL)=Sqrt[PPtabFO[1,1]]; the experimental uncertainty for O4 is del(O4)=Sqrt[PPtabFO[2,2]]; etc. *)
 (* --> a(sigmax)=aALfFOO10[8] +- Sqrt[PPtabFO[1, 1]]=0.105591+-0.000975657_{exp}; --> alpha_s(mtau^2)=0.320192+-0.002827_{exp};  19 Sept. 2022 *)

 (* When we increase NNAL=77 to NNAL=79 (i.e., include the last two bins) --> alO10fitPFOMinnum1lN79[8]={0.0000936873, {aALf -> 0.102933, O4fitf -> -0.00453871, O6fitf -> 0.00711644, O8fitf -> -0.00340316, O10fitf -> 0.00145154}}; and we get  Sqrt[PPtabFO[1,1]]=0.00138495; alpha_s(mtau^2)=0.323369+-0.004351_{exp}; i.e., del(alpha_s)_{exp} increases (in comparison to the NNAL=77 case) by factor 0.004351/ 0.002827=1.539; *)
 (*-------------------------*)
 (* Summary of the experimental uncertainties, 'O10' approach with Nt=8: *)
 (* alpha_s(mtau^2)=0.320192+-0.002827_{exp} 
 10^3 O4= -4.2 +- 0.5_{exp};
 10^3 O6= +6.0 +- 0.8_{exp};
 10^3 O8= -2.6 +-0.4_{exp};
  10^3 O10= +1.0 +- 0.2_{exp}; *)
 (* 19 Sept. 2022; checked on 26 Sept. 2022 *)
   (*---------------------------------------------*)
 (* The precise central value result of the FO global fit: *)
 (* alO10fitPFOMinnum[8]={0.000025670889484749748`, {aALf -> 0.1055909512932032`, O4fitf -> -0.00421526293750236`, O6fitf -> 0.005957039292628349`, O8fitf -> -0.002565953796231432`, O10fitf -> 0.0009767589344250825`}}; *)
 (* 19 Sept. 2022 *)
 (*---------*)
(* Variation of kappa, FO global (at Nt=8): *)
 (* kappa=2.0:    (N0=8)
  alO10fitPFOMinnum[N0]={0.0000158274, {aALf -> 0.0874284, O4fitf -> -0.00241598, O6fitf -> 0.00431012, O8fitf -> -0.00191064, O10fitf -> 0.000734355}};  aALf(sigmax)=0.105123; alpha_s(mtau^2)=0.318834; *)
(* kappa=1.0 (the central case of the global fit):
 alO10fitPFOMinnum[N0]={0.0000256709, {aALf -> 0.105591, O4fitf -> -0.00421526, O6fitf -> 0.00595704, O8fitf -> -0.00256595, O10fitf -> 0.000976759}};  aALf(sigmax)=0.105591; alpha_s(mtau^2)=0.320192; *)
(* kappa=2/3: 
 alO10fitPFOMinnum[N0]={0.0000260995, {aALf -> 0.122394, O4fitf -> -0.00633303, O6fitf -> 0.00724463, O8fitf -> -0.00295056, O10fitf -> 0.00108021}}; a(sigmax)=0.106853; *)
 (* kappa=0.5:  
    alO10fitPFOMinnum[N0]={0.0000125881, {aALf -> 0.140486, O4fitf -> -0.00937426, O6fitf -> 0.00812192, O8fitf -> -0.00293817, O10fitf -> 0.000972757}}; a(sigmax)=0.108612; alpha_s(mtau^2)=0.328926; *)
   
(* Note that aALf=a(kappa*sigmax) *)
(* --> thus the variations: (kappa=1+{^{+1}_{-1/2}}) *)
(*    alpha_s(mtau^2)=0.320192+{^{-0.001358}_{+0.008734}}_{kap};
      10^3 O4= -4.2+{^{+1.8}_{-5.2}}_{kap};  
      10^3 O6= +6.0+{^{-1.7}_{+2.1}}_{kap};  
      10^3 O8 =-2.6+{^{+0.7}_{-0.3}}_{kap};
      10^3 O10=+1.0+{^{-0.3}_{+0.0}}_{kap}; 
 *)
 (* Checked on 26 Sept. 2022; correction found in del(alpha_s):^{-0.001358} *)
 
 (* if kappa=1+{^{+1}_{-1/3}}: *)
 (* alpha_s(mtau^2)=0.320192+{^{-0.001358}_{+0.003653}}_{kap};
      10^3 O4= -4.2+{^{+1.8}_{-2.1}}_{kap};  
      10^3 O6= +6.0+{^{-1.7}_{+1.2}}_{kap};  
      10^3 O8 =-2.6+{^{+0.7}_{-0.4}}_{kap};
      10^3 O10=+1.0+{^{-0.3}_{+1.0}}_{kap}; 
  *)
 (* 19 Sept. 2022 *)

(*-----------*)
(* Variation of d4, FO global (NOTE, at Nt=8 now) d4=275 +- 63: *)

 (* when ddtk[4, tkap]=275.+63.=338.: 
  alO10fitPFOMinnum[8]= {0.0000228748, {aALf -> 0.102592, O4fitf -> -0.00359833, O6fitf -> 0.00543239, O8fitf -> -0.00234516, O10fitf -> 0.000900422}}; alpha_s(mtau^2)=0.311491; *)
  (* when the central case d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=275.: 
   alO10fitPFOMinnum[8]={0.0000256709, {aALf -> 0.105591, O4fitf -> -0.00421526, O6fitf -> 0.00595704, O8fitf -> -0.00256595, O10fitf -> 0.000976759}}; alpha_s(mtau^2)=0.320192; *)
  (* when ddtk[4, tkap]=ddtk[4,Exp[tKK]]=275.-63.=212.: 
   alO10fitPFOMinnum[8]={0.0000245752, {aALf -> 0.107399, O4fitf -> -0.00451787, O6fitf -> 0.00611531, O8fitf -> -0.00261271, O10fitf -> 0.000982845}};  alpha_s(mtau^2)=0.325422; *)

(* alpha_s(mtau^2)=0.320192+{^{-0.008701}}_{+0.005230}}_{d4}; 
 10^3 O4= -4.2+{^{+0.6}_{-0.3}}_{d4}; 
 10^3 O6= +6.0+{^{-0.6}_{+0.1}}_{d4}; 
 10^3 O8 =-2.6+{^{+0.3}_{-0.0}}_{d4}; 
 10^3 O10=+1.0+{^{-0.1}_{+0.0}}_{d4};
*)
 (* 19 Sept. 2022; checked on 26 Sept. 2022 *)
 (*--------------*)
(* Variation of Nt(truncation order, Nt=8+-1 and Nt=8 +-2, Nt=8-3), FO global: *)
(*
alO10fitPFOMinnum[10]={0.0000160259, {aALf -> 0.106541, O4fitf -> -0.00692496, O6fitf -> 0.00646856, O8fitf -> -0.00248621, O10fitf -> 0.000878739}}; alsmtauFOO10[10]=0.322943;
alO10fitPFOMinnum[9]={0.0000227962, {aALf -> 0.106254, O4fitf -> -0.0055302, O6fitf -> 0.00645177, O8fitf -> -0.00264725, O10fitf -> 0.000976554}}; alsmtauFOO10[9]=0.32211;
The central case:
 alO10fitPFOMinnum[8]={0.0000256709, {aALf -> 0.105591, O4fitf -> -0.00421526, O6fitf -> 0.00595704, O8fitf -> -0.00256595, O10fitf -> 0.000976759}}; alsmtauFOO10[8]=0.320192;
alO10fitPFOMinnum[7]={0.0000238144, {aALf -> 0.104919, O4fitf -> -0.00303407, O6fitf -> 0.00514375, O8fitf -> -0.00229749, O10fitf -> 0.000894026}};  alsmtauFOO10[7]=0.318243;
alO10fitPFOMinnum[6]={0.0000168579, {aALf -> 0.104379, O4fitf -> -0.00192045, O6fitf -> 0.00403484, O8fitf -> -0.00184384, O10fitf -> 0.000723842}}; alsmtauFOO10[6]=0.316678;
alO10fitPFOMinnum[5]={7.48823*10^-6, {aALf -> 0.104293, O4fitf -> -0.00084901, O6fitf -> 0.00268136, O8fitf -> -0.00122807, O10fitf -> 0.000473255}}; alsmtauFOO10[5]=0.31643;
 *)

(*  alpha_s(mtau^2)=0.320192+{^{+0.002751}_{-0.003762}}_{Nt};
 10^3 O4= -4.2+{^{-2.7}_{+3.4}}_{Nt}; 
 10^3 O6= +6.0+{^{+0.5}_{-3.3}}_{Nt}; 
 10^3 O8= -2.6 +{^{+0.1}_{+1.4}}_{Nt}; 
 10^3 O10=+1.0+{^{-0.1}_{-0.5}}_{Nt};
*)
 (* 19 Sept. 2022; checked on 26 Sept. 2022 *)

 (* Note: Nt=8^{+2}_{-3} *)
 (*--------------*)
 (* Variation when going beyond the O10 truncation, including O12 in OPE: *)
 (* N0=8; alO12fitPFOMinNum[N0]=alO12fitPFOMinnum[N0];
  N0=8; *)
 (* alO12fitPFOMinnum[8]={2.70903*10^-7, {aALf -> 0.106407, O4fitf -> -0.00479374, O6fitf -> 0.00751318, O8fitf -> -0.00409064, O10fitf -> 0.00272721, O12fitf -> -0.00120311}};  alsmtauFOO12[N0]=0.322553; *)
 (* The central case:
  alO10fitPFOMinnum[8]={0.0000256709, {aALf -> 0.105591, O4fitf -> -0.00421526, O6fitf -> 0.00595704, O8fitf -> -0.00256595, O10fitf -> 0.000976759}}; alsmtauFOO10[N0]=0.320192; *)
 (* --> alsmtauFOO12[N0]-alsmtauFOO10[N0]=0.322553-0.320192=+0.002361; *)
 (* 10^3 (O4FOO12[N0]-O4FOO10[N0])=-0.57 =app -0.6 *)
 (* 10^3 (O6FOO12[N0]-O6FOO10[N0])= +1.56 =app +1.6 *)
 (* 10^3 (O8FOO12[N0]-O8FOO10[N0])= -1.52 =app -1.5 *)
 (* 10^3 (O10FOO12[N0]-O10FOO10[N0])= +1.75 =app +1.8 *)
 (* -->
  alpha_s(mtau^2)=0.320192+-0.002361_{O12};
  10^3 O4= -4.2 -+0.6_{O12}; 
  10^2 O6= +6.0 +-1.6_{O12};
  10^2 O8= -2.6 -+1.5_{O12};
  10^3 O10=+1.0 +-1.8_{O12};
*)
  (* 27 Sept. 2022 *)
 (*--------------*)
 (* Variation when increasing the number of data bins from  NNAL=77 tyo NNAL=79 (i.e., including the last two bins, which have bad data): *)
 (* N0=8; NNAL=79: alO10fitPFOMinnum1lNNAL79[N0] = alO10fitPFOMinNum[N0] *)
 (* alO10fitPFOMinnum1lNNAL79[8]={0.0000936873, {aALf -> 0.102933, O4fitf -> -0.00453871, O6fitf -> 0.00711644, O8fitf -> -0.00340316, O10fitf -> 0.00145154}}; alpha_s(mtau^2)=0.32337; (note that sigmax=3.15708 when NNAL=79); aALf=a(sigmax) *)
  (* The central case:
  alO10fitPFOMinnum[N0]={0.0000256709, {aALf -> 0.105591, O4fitf -> -0.00421526, O6fitf -> 0.00595704, O8fitf -> -0.00256595, O10fitf -> 0.000976759}}; alsmtauFOO10[N0]=0.320192; *)
 (* --> alpha_s(mtau^2)-alsmtauFOO10[N0]=+0.003178; *)
 (* 10^3 (O4fitf-O4FOO10[N0])= -0.32 =app -0.3 *)
 (* 10^3 (O6fitf-O6FOO10[N0])= +1.16 =app +1.2 *)
 (* 10^3 (O8fitf-O8FOO10[N0])= -0.84 =app -0.8 *)
 (* 10^3 (O10fitf-O10FOO10[N0])= +0.47 =app +0.5 *)
  
 (* -->
  alpha_s(mtau^2)=0.320192+-0.003178_{Nbin}
  10^3 O4= -4.2 -+0.3_{Nbin}; 
  10^2 O6= +6.0 +-1.2_{Nbin};
  10^2 O8= -2.6 -+0.8_{Nbin};
  10^3 O10=+1.0 +-0.5_{Nbin};
*) 
 (* 27 Sept. 2022 *)
  (*-----------------------------------------*)
 (* Altogether, the FO global fit approach with (O6-O10; O6 with nonzero anom.dim) [and a(Q^2) in D=6 part runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ] gives: *)
(* FO
 Nt=8;  (kap=1+{^{+1}_{-1/2}}; d4=275+-63; Nt=8^{+2}_{-3}):
 alpha_s(mtau^2)=0.320192+-0.002827_{exp}+{^{-0.001358}_{+0.008734}}_{kap}+{^{-0.008701}}_{+0.005231}}_{d4}+{^{+0.002751}_{-0.003762}}_{Nt}+-0.002361_{O12}+-0.003178_{Nbin} =0.320192+{^{+0.011614}_{-0.0107412}}
 = 0.3202+{^{+0.0116}_{-0.0107}} =app 0.320+{^{+0.012}_{-0.011}};
 10^3 O4=-4.2 +- 0.5_{exp}+{^{+1.8}_{-5.2}}_{kap}+{^{+0.6}_{-0.3}}_{d4}+{^{-2.7}_{+3.4}}_{Nt}-+0.6_{O12}-+0.3_{Nbin} =-4.2+{^{+4.0}_{-5.9}};
  10^3 O6=+6.0 +- 0.8_{exp}+{^{-1.7}_{+2.1}}_{kap}+{^{-0.6}_{+0.1}}_{d4}+{^{+0.5}_{-3.3}}_{Nt}+-1.6_{O12}+-1.2_{Nbin}  =+6.0+{^{+2.3}_{-3.8}};
 10^3 O8=-2.6 +-0.4_{exp}+{^{+0.7}_{-0.3}}_{kap}+{^{+0.3}_{-0.0}}_{d4}+{^{+0.1}_{+1.4}}_{Nt}-+1.5_{O12}-+0.8_{Nbin} =-2.6+{^{+2.3}_{-1.8}};
10^3 O10=+1.0 +- 0.2_{exp}+{^{-0.3}_{+0.0}}_{kap}+{^{-0.1}_{+0.0}}_{d4}+{^{-0.1}_{-0.5}}_{Nt}+-1.8_{O12}+-0.5_{Nbin} =+1.0+{^{+1.9}_{-2.0}};
 chi^2=0.0000256709;
 *)

 (* 19 Sept. 2022; checked on 26 Sept. 2022; added _{O12} and _{Nbin} on 27 Sept. 2022  *)

 (* End of FO global fit results [with 'p3ka1ka2' (in D=0), O6 with a(Q^2)^(1-ka2) ('O6ka2'), OPE up to O10;  a(Q^2)^(1-ka2) in D=6 OPE term runs by one-loop RGE from a(sigmax)=aALfFO[Nt] to a(sigmax Exp[I fi]) ], OK *)

(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
 (*-----------------------------------------------------------------------*)
 
 (* Double-pinched ('P') Borel-Laplace in PV(+CI) approach (PPV): *)

(* kappa=Exp[-tKK]; *)
  kappa=1;  (* kappa = mu^2/Q^2, where mu^2 is RScl *)
  tkap=kappa Exp[tKK];

(* First for FIXED a(kappa sigmax)=aAL *)
sol2a[aALf_]:= NDSolve[{-I y'[th] == bta[y[th]], y[0] == aALf}, y, {th, -Pi, +Pi}];
(*'*) (* aALf=a[kappa sigmax] *)
   (*  aCI[th_? NumberQ,aALf_? NumberQ]:=First[y[th]/.sol2a[aALf]]; *)
 aCI[th_? NumberQ, aALf_] := First[y[th] /. sol2a[aALf]];
(* aALf=a(kappa sigmax); aCI[th,aALf]=a(Q2=kappa sigma Exp[I th]) *)

(* aAL=a(kappa*sigmax), FIXED here: *)  
  (* aAL=0.10;  *)
 aAL=0.1060905;
(* aCI[th,aAL] = a(kappa sigmax Exp[I th]), and aAL=a(kappa sigmax) *)
  
  delD[th_,Nn_]:= Sum[delddtk[jj-1,tkap] aCI[th,aAL]^jj, {jj, 1, Nn}];
  (* delD[th,N] is the polynomial part of the Borel od Adler D, of order up to a(kappa sigmax Exp[I th])^N; note that kappa=tkap Exp[-tKK] (kappa=mu^2/Q^2) *)
  delDext[th_,Nn_,aALf_]:= Sum[delddtk[jj-1,tkap] aCI[th,aALf]^jj, {jj, 1, Nn}];

  DBorelep[th_, ep_]:=(1/be0)*(NIntegrate[BDsing[u,tkap] Exp[-u/(be0 aCI[th,aAL])], {u, 0., I ep}] + NIntegrate[BDsing[u,tkap] Exp[-u/(be0 aCI[th,aAL])], {u, I ep, I ep +10^3}]);
  
th0=0.5;
DBorel1pl[th_]:=(1/be0*Exp[I th0])*(NIntegrate[BDsing[u Exp[I th0],tkap] Exp[-u Exp[I th0]/(be0 aCI[th,aAL])], {u, 0., 10.^3}]);
DBorel1mi[th_]:=(1/be0*Exp[-I th0])*(NIntegrate[BDsing[u Exp[-I th0],tkap] Exp[-u Exp[-I th0]/(be0 aCI[th,aAL])], {u, 0., 10.^3}]);
  
DBorel1plext[th_,aALf_]:=(1/be0*Exp[I th0])*(NIntegrate[BDsing[u Exp[I th0],tkap] Exp[-u Exp[I th0]/(be0 aCI[th,aALf])], {u, 0., 10.^3}]);
DBorel1miext[th_,aALf_]:=(1/be0*Exp[-I th0])*(NIntegrate[BDsing[u Exp[-I th0],tkap] Exp[-u Exp[-I th0]/(be0 aCI[th,aALf])], {u, 0., 10.^3}]);
(* Note, by construction,  DBorel1mi[-th]=Conjugate[DBorel1pl[th]] *)
  (* aALf=a(kappa sigmax)= a(tkap Exp[tKK]) *)
  
DBorelS[th_]:=(DBorel1pl[th]+DBorel1mi[th])/2;
DBorel[th_,Nm_]:=DBorelS[th]+delD[th,Nm];(* This is PV *)
  (* DBorel[-th,Nm]=Conjugate[DBorel[th,Nm]]; *)
  (* aALf=a(kappa sigmax)= a(tkap Exp[tKK]) *)
  
DBorelSext[th_,aALf_]:=(DBorel1plext[th,aALf]+DBorel1miext[th,aALf])/2;
DBorelext[th_,Nm_,aALf_]:=DBorelSext[th,aALf]+delDext[th,Nm,aALf];(* This is PV *)

  delDBorel[th_]:=(DBorel1pl[th]-DBorel1mi[th])/2/Pi/I; 
(* Note, delD[-th,10]=Conjugate[delD[th,10]];  *)
(* delDBorel[-th]=Conjugate[delDBorel[th]] *) 
(* delDBorel[th]=delDBorelIR2plplIR3mimi[th], i.e., the difference when going above (Im u > 0) for both IRp=2 and IRp=3, and going below (Im u < 0) for both IRp=2 and IRp=3 *)

 delDBorelext[th_,aALf_]:=(DBorel1plext[th,aALf]-DBorel1miext[th,aALf])/2/Pi/I; 
  (*------------*)
  (* Now generation of the tables of Borel-resummed Adler function along the contour of radius kappa*sigmax, for fixed a(kappa sigmax)=aAL: *)

 Nn=60; 
  
Clear[tabDBorelS];  
  tabDBorelS=Table[{(jj-1)*Pi/Nn,  DBorelS[N[(jj-1)*Pi/Nn]]}, {jj, 1, Nn}]; 
  (* tabDBorel[[jj]][[2]]=DBorel[(jj-1)Pi/Nn] *)
  
Clear[deltabDBorelS];
  tabdelDBorel=Table[{(jj-1)*Pi/Nn,  delDBorel[N[(jj-1)*Pi/Nn]]}, {jj, 1, Nn}];
  (* tabdelDBorel[[jj]][[2]]=delDBorel[(jj-1)Pi/Nn] *)
  
Clear[tabaCI];
  Do[tabaCI[nn]=Table[{(jj-1)*Pi/Nn,  aCI[N[(jj-1)*Pi/Nn],aAL]^nn}, {jj, 1, Nn}],{nn,1,15}];
  (* tabaCI[nn][[jj]][[2]] = a(Q^2=sigmax (jj-1)Pi/Nn)^nn, such that aCI(sigmax)=aAL *)
  
Clear[tabdelD];
  Do[tabdelD[nn]=Table[{(jj-1)*Pi/Nn,  Sum[delddtk[mm-1,tkap] tabaCI[mm][[jj]][[2]],{mm,1,nn}]}, {jj, 1, Nn}],{nn,3,15}];
  (* tabdelD[nn][[jj]][[2]]]= Sum[delddtk[mm-1,tkap] a(Q^2=(jj-1)Pi/Nn)^mm,{mm,1,nn} *)
  
   Do[tabDBorel[nn]=Table[{(jj-1)*Pi/Nn, tabDBorelS[[jj]][[2]]+tabdelD[nn][[jj]][[2]]},{jj,1,Nn}],{nn,3,15}];
   (* Note: tabDBorel[nn][[jj]][[2]]= the value DBorel[th=(jj-1)*Pi/Nn,N]} (jj=1, ...,Nn), where DBorel has the polynomial corrective part with power up to (and including) aCI[th,aAL]^N *)

Clear[Nn];
  
DBorSup=Interpolation[tabDBorelS];
DBorSfast[th_] :=If[th >= 0, DBorSup[th], Conjugate[DBorSup[-th]]];
(* DBorSfast[th] is PV  Borel resummation of Adler d(Q^2=sigmax Exp[I th]) , FAST evaluation, the singular (S) part, i.e., with NO polynomial (correction) part *)

  
Do[DBornn[nn]=Interpolation[tabDBorel[nn]],{nn,3,15}];
DBorfast[th_,nn_] :=If[th >= 0, DBornn[nn][th], Conjugate[DBornn[nn][-th]]];
(* DBorfast[th,nn]: This is PV Borel resummation of Adler d(Q^2=sigmax Exp[I th]) , FAST evaluation, with the polynomial (correction) part delD[th,nn] truncated at a[kappa sigmax Exp[I th]]^nn; note that Adler d(Q^2) = a(kappa Q^2) + O(a^2) i.e., without the term 1 (=O(a^0) *)
 (* It should approximate well the ('exact') expression DBorel[th,nn], which is evaluated much more SLOWLY; *)

delDBorSup=Interpolation[tabdelDBorel];
delDBorfast[th_] :=If[th >= 0, delDBorSup[th], Conjugate[delDBorSup[-th]]]; 
(* delDBorfast[th] is the renormalon-uncertainty measure of DBorfast[th] *)
 (* It is equal to delDBorel[th], but much faster evaluation *)
 (*-----------*)
 (* Now double-pinched ('P') Borel-Laplace, in the Borel-resummed ('PV') approach  (-> 'PPV'): *) 
(* Here: M2 = M^2 = |M^2| Exp[I Psi] (complex), where later Psi=Pi/6,Pi/4,0: *)
 (* fixed aAL=a(kappa sigmax) *)
 

BPPVD0[M2_?NumberQ, Nn_]:= (1/(2 Pi)) NIntegrate[(((1 + Exp[I th])^2 - 2 (M2/sigmax) (1 + Exp[I th]) + 2 (M2/sigmax)^2) Exp[(sigmax/M2) Exp[I th]] - 2 (M2/sigmax)^2 Exp[-(sigmax/M2)]) DBorfast[th,Nn], {th, -Pi, Pi}];
 (*  BPPVD0sl[M2_?NumberQ, Nn_,aALf_?NumberQ]:= (1/(2 Pi)) NIntegrate[(((1 + Exp[I th])^2 - 2 (M2/sigmax) (1 + Exp[I th]) + 2 (M2/sigmax)^2) Exp[(sigmax/M2) Exp[I th]] - 2 (M2/sigmax)^2 Exp[-(sigmax/M2)]) DBorelext[th,Nn,aALf], {th, -Pi, Pi}]; *)
 
BPPVdelD0[M2_?NumberQ]:= (1/(2 Pi)) NIntegrate[(((1 + Exp[I th])^2 - 2 (M2/sigmax) (1 + Exp[I th]) + 2 (M2/sigmax)^2) Exp[(sigmax/M2) Exp[I th]] - 2 (M2/sigmax)^2 Exp[-(sigmax/M2)]) delDBorfast[th], {th, -Pi, Pi}];

 
 (* This is NEW (16 Sept. 2021) (approach with O4-O10, with D=6 with a(Q^2)^(1-ka2) O6): *)    
 ReBPPVNumO10[M2f_,Nn_,O4_,O6_,O8_,O10_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPPVD0[M2f,Nn]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/3!) (1/M2f^4) (1+ 2*3 (M2f/sigmax) + 3*2 (M2f/sigmax)^2) + 2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2) + (6 Pi^2 O6/sigmax^3) (1/2/Pi) NIntegrate[ ((aMSext[(1/kappa), aAL]/(1 + I be0 aMSext[(1/kappa), aAL] fi))^(1-ka2))  Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
 ReBPPVNumO12[M2f_,Nn_,O4_,O6_,O8_,O10_,O12_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPPVD0[M2f,Nn]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/3!) (1/M2f^4) (1+ 2*3 (M2f/sigmax) + 3*2 (M2f/sigmax)^2) + 2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2)+ 2 Pi^2 O12  (1/5!) (1/M2f^6) (1+ 2*5 (M2f/sigmax) + 5*4 (M2f/sigmax)^2) + (6 Pi^2 O6/sigmax^3) (1/2/Pi) NIntegrate[ ((aMSext[(1/kappa), aAL]/(1 + I be0 aMSext[(1/kappa), aAL] fi))^(1-ka2))  Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
 
 (* Now search for condensate values with this Borel-Laplace (double) pinched ('P') approach, with (Borel-resummed) PV method for D=0 part: *)
 (* Nm=NchP+1; Nd=3 Nm;  
M2al0[ii_]:=0.90+(ii-1)*0.60/(Nm-1);
  Do[M2al[jj]=M2al0[jj],{jj,1,Nm}]; Do[M2al[jj+Nm]=M2al0[jj] Exp[I Pi/6],{jj,1,Nm}]; Do[M2al[jj+2 Nm]=M2al0[jj] Exp[I Pi/4],{jj,1,Nm}]; *)
  (* We will use sum of squares minimization in this interval, for M^2=|M^2| Exp[I Psi] with Psi=0, Pi/6, Pi/4 one after another in the sum (it turns out that it is not good to include Psi=Pi/3, this value appears to be too large, and xi^2 become large): *)

O4init=-0.003; O6init=-0.001; 
 O8init=0.001; O6init = 0.003; O10init=0.005; O12init=0.002;
 (* aAL=a(kappa*sigmax) is fixed  above *)
 O12fitPPVMinNum[Nn_]:= FindMinimum[Sum[(ReBPPVNumO12[M2al[jj],Nn,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init}}];
 O10fitPPVMinNum[Nn_]:= FindMinimum[Sum[(ReBPPVNumO10[M2al[jj],Nn,O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init}}];
 O8fitPPVMinNum[Nn_]:= FindMinimum[Sum[(ReBPPVNumO10[M2al[jj],Nn,O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}];
 
(* O10fitPPVMinnum[nn]=O10fitPPVMinNum[nn] (nn=3,...,12) *)
 (* 16 Sept. 2022 *)

 (* Do[O10fitPPVMinnum[Nn] = O10fitPPVMinNum[Nn], {Nn, 3, 12}];
  Do[xi2PV[Nn]=O10fitPPVMinnum[Nn][[1]], {Nn, 3, 12}]; *)
 (* This then gives results for O10fitPFOMinnum[3], ..., O10fitPFOMinnum[12]; at a(sigmax)=aAL value *)
 (* end of BL results (O6-O10) approach in PPV *)
 (* 16 Sept. 2022 *)
 
 (*---------------------------------------------------------------------_*)
  (*---------------------------------------------------------------------_*)
  (*---------------------------------------------------------------------_*)

  (*--------------------------*)
(* NEW BL results, PPV approach with (O61,O62): *)
 (* Summary, for kappa=1 and NchP=2 (Nd=9), PV(+CI) approach, case 'p3ka1ka2' (D=0) & 'O6ka1' O10 (O6 with a(Q^2)^(1-ka2), OPE includes up to O10): *)

  (* NchP=2: *)
 (*
 aALfPV[3]=0.107331; O10fitPPVMinnum[3]={0.000012560181137155481`, {O4fitf -> -0.0028228346980238937`, O6fitf -> 0.00564744262708358`, O8fitf -> -0.0018061898714896896`, O10fitf -> 0.0007619409644375851`}};
  aALfPV[4]=0.106206; O10fitPPVMinnum[4]={0.000012046548491781471`, {O4fitf -> -0.0029177004063317747`, O6fitf -> 0.005429924152810957`, O8fitf -> -0.0018060500625300354`, O10fitf -> 0.0007353989048967403`}};
  aALfPV[5]=0.106151; O10fitPPVMinnum[5]={0.0000118900239272905`, {O4fitf -> -0.0029311823066489155`, O6fitf -> 0.005417618551128278`, O8fitf -> -0.0018082414972161534`, O10fitf -> 0.0007329289124612002`}};
  aALfPV[6]=0.1060765; O10fitPPVMinnum[6]={0.000011761439320603452`, {O4fitf -> -0.0029640001358094306`, O6fitf -> 0.005384564342874001`, O8fitf -> -0.001816019355122791`, O10fitf -> 0.0007251522512722274`}};
  aALfPV[7]=0.1060905;O10fitPPVMinnum[7]={0.000011740745332238454`, {O4fitf -> -0.0029547054662461976`, O6fitf -> 0.005398958409024923`, O8fitf -> -0.0018134852463415506`, O10fitf -> 0.0007297809837811238`}};
  aALfPV[8]=0.106073;O10fitPPVMinnum[8]={0.000011929089705661197`, {O4fitf -> -0.0029693501508607154`, O6fitf -> 0.005361127805697869`, O8fitf -> -0.001818293096269184`, O10fitf -> 0.0007133425310656665`}};
  aALfPV[9]=0.106087;O10fitPPVMinnum[9]={0.000011520622985178095`, {O4fitf -> -0.0029599736474723525`, O6fitf -> 0.0054095220030782815`, O8fitf -> -0.00181536054757717`, O10fitf -> 0.0007413539664213614`}};
  aALfPV[10]=0.1060655;O10fitPPVMinnum[10]={0.000012557827754857217`, {O4fitf -> -0.002963918897805264`, O6fitf -> 0.005328057951768351`, O8fitf -> -0.0018125594894732406`, O10fitf -> 0.0006805244433931174`}};
  aALfPV[11]=0.1061055;O10fitPPVMinnum[11]={0.000010344642055837876`, {O4fitf -> -0.002972553658784546`, O6fitf -> 0.005457890098766144`, O8fitf -> -0.0018346398662830942`, O10fitf -> 0.0007994127239800274`}};
  aALfPV[12]=0.106053;O10fitPPVMinnum[12]={0.000013471260329678357`, {O4fitf -> -0.0029512540319802092`, O6fitf -> 0.005317179576628503`, O8fitf -> -0.0017887744202082865`, O10fitf -> 0.0006509648150917149`}}; 

 Clear[Nn]; Do[Module[{}, xi2PVO10[Nn]=O10fitPPVMinnum[Nn][[1]]; aALfPVO10[Nn]=aALfPV[Nn]; alsmtauPVO10[Nn]=aMSext[mtau^2/(kappa sigmax),aALfPVO10[Nn]] Pi;  O4PVO10[Nn]=O4fitf/.O10fitPPVMinnum[Nn][[2]]; O6PVO10[Nn]=O6fitf/.O10fitPPVMinnum[Nn][[2]]; O8PVO10[Nn]=O8fitf/.O10fitPPVMinnum[Nn][[2]]]; O10PVO10[Nn]=O10fitf/.O10fitPPVMinnum[Nn][[2]],{Nn,3,12}];
  *)
 
 (* 19 Sept. 2022, OK *)

 (* end of NEW BL results approach in PPV, for the case 'p3ka1ka2' & 'O6ka2' (D=6 OPE with a(Q^2)^(1-ka2) <O6>/Q^6 ('O6ka2')) and OPE up to O10 *)
  (*---------------------------------------------------------------------_*)
  (*---------------------------------------------------------------------_*)
 (* Now FESR (2,n) moments a^{(2,n)} (i.e., double-pinched) in CI: *)
 (*--------*)
 (* First the experimental ALEPH values of (2,n) moments a^{(2,n)} (see also file rtauALEPH.m): *)
grs0a[2,nn_][sig_]:=((nn+3)/(nn+1)) (1/sigmax) (1 - (nn+2) (sig/sigmax)^(nn+1) + (nn+1) (sig/sigmax)^(nn+2) );
aexpBin77[2,nn_]:=(Sum[ dsbin[[ii]] grs0a[2,nn][sbin[[ii]]] OmegaVplA[ii],{ii,1,77}] + 2 Pi^2 fPi^2 grs0a[2,nn][mPi^2]) -  1;
delaexpBin77[2,nn_]:=(2 Pi^2) Sqrt[Sum[ grs0a[2,nn][sbin[[ii]]] grs0a[2,nn][sbin[[jj]]] dsbin[[ii]] dsbin[[jj]](  covVV[ii,jj]+covAA[ii,jj]+ covVA[ii,jj]+covAV[ii,jj] ),{ii,1,77},{jj,1,77}]];
(*  aexpBin77[2,n] +- delaexpBin77[2,n] = a^{(2,n)}_{exp} *)
(* Note: aexpBin77[2, 0] + -delaexpBin77[2, 0] = 0.24038 + -0.00261;
          aexpBin77[2, 1] + -delaexpBin77[2, 1] = 0.20610 + -0.00260; 
          aexpBin77[2, 2] + -delaexpBin77[2, 2] = 0.19603 + -0.00282;
          aexpBin77[2, 3] + -delaexpBin77[2, 3] = 0.18965 + -0.00313; *)	(* 19 Sept. 2022 *)
 (*---*)
  (*-----------------------*)
 (* moments a^{(2,n}}, PV, at fixed truncation order and for aAL=(a(kappa*sigmax))=aALfPV[Nt]: *)			     

 aaPVD0fast[2,nn_][Nt_]:=(1/(2 Pi)) ((nn+3)/(nn+1))  NIntegrate[ ( Exp[I th] (1 - (-Exp[I th])^(nn+1)) + ((nn+1)/(nn+3))(1 - (-Exp[I th])^(nn+3)) ) DBorfast[th,Nt], {th, -Pi, Pi}];
(* aAL(=a(kappa sigmax)) =aALfPV[Nt], with Nt chosen FIXED truncation order *)
(*  aaPVD0fast[2,nn][Nt] is a^{(2,nn)} moment calculated with the value of  a(kappa sigmax) equal to aALf, calculated to truncation order Nt (~a^Nt) *)
 (*-------------------------*)
(* First a^{(2,0)}(sigmax)  (Nn is the truncation index Nt of the D=0 part): *)
aaPV20[Nn_, O4f_,O6f_]:=aaPVD0fast[2,0][Nn]+ ((0+3)/(0+1)) 2 Pi^2 (-1)^0 ( (0+2) O4f/sigmax^(0+2) ) + 3 (1/2/Pi) (6 Pi^2 O6f/sigmax^3) ( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[I 3 fi]) ) ((aAL/(1 + I aAL be0 fi))^(1-ka2)),{fi,-Pi,Pi}] );
aaPV21[Nn_, O6f_,O8f_]:=aaPVD0fast[2,1][Nn]+ 2 (1/2/Pi) (6 Pi^2 O6f/sigmax^3)*( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 - Exp[2 I fi]) + (1/2) (1 - Exp[I 4 fi]) ) ( (aAL/(1 + I aAL be0 fi))^(1-ka2)),{fi,-Pi,Pi}] ) + ((1+3)/(1+1)) 2 Pi^2 (-1)^1 ( (1+1) O8f/sigmax^(1+3) );			     
(* a^{(2,2)}(sigmax): *)
aaPV22[Nn_, O6f_,O8f_,O10f_]:=aaPVD0fast[2,2][Nn]+ (6 Pi^2 O6f/sigmax^3)*(1/2/Pi)*( NIntegrate[ Exp[-I 3 fi] ((aAL/(1 + I aAL be0 fi))^(1-ka2)) ( ((2+3)/(2+1)) Exp[I fi] (1 +(-1)^2 Exp[(2+1) I fi]) +  (1 + (-1)^2 Exp[I (2+3) fi]) ),{fi,-Pi,Pi}] ) + ((2+3)/(2+1)) 2 Pi^2 (-1)^2 ( (2+2) O8f/sigmax^(2+2) + (2+1) O10f/sigmax^(2+3));
(* 19 Sept. 2022 *)
 (* end of a2n in global PV approach with (O6ka1=O61, O6ka2=O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ] *)
 
 (* Note that aaPV2n are valid only when kappa=1, i.e., when aAL=a(sigmax) *)
  (*------------------------*)
   (* Tables of a20, a21 vs Nt in PV when aALf=a(sigmax) is VARIED ('avar'), i.e.,using alO4O6fitPPVMinnum[N] (N=3,4,...12) from global PVfit (with a(sigmax Exop[I fi]) in D=6 part of OPE running by one-loop RGE): *)			     

 (* Note, have to adjust aAL to aAL=aALfPV[NN] at each NN and calculate at that NN first the DBorfast etc. *)
 
 (*
  Tab20D0PVO10avar(*=Table[{NNt,Re[aaPVD0fast[2, 0][NNt]]},{NNt,3,12}] ={{3,0.274445}, {4,0.276218}, {5,0.276441}, {6,0.276988}, {7,0.276826}, {8,0.277104}, {9,0.27689}, {10,0.277071}, {11,0.27701}, {12,0.276894}} *)={{3,0.27444535682789517`}, {4,0.2762184471671655`}, {5,0.2764414750949396`}, {6,0.27698836534923693`}, {7,0.2768262026535546`}, {8,0.27710404542873784`}, {9,0.2768904788836259`}, {10,0.2770705342011809`}, {11,0.27701012812005876`}, {12,0.27689411746580644`}};
  Tab20PVO10avar(*=Table[{NNt,Re[aaPV20[NNt,O4PVO10[NNt],O6PVO10[NNt]]]},{NNt,3,12}] ={{3, 0.233913}, {4, 0.239662}, {5, 0.239961}, {6, 0.24037}, {7, 0.240293}, {8, 0.240389}, {9, 0.240312}, {10, 0.240429}, {11, 0.240213}, {12, 0.240496}} *)={{3, 0.23391293216089695`}, {4, 0.23966193213334558`}, {5, 0.2399607079020776`}, {6, 0.24036976248077535`}, {7, 0.24029250776096564`}, {8, 0.24038888813730283`}, {9, 0.24031229516121666`}, {10, 0.24042921114043955`}, {11, 0.24021276122216786`}, {12, 0.24049648924616127`}};
  (* Several local minima of a20 when Nt increases; the first local minimum at Nt=7; then at Nt=9, 11; see also: delTab20PVO10avar= Table[{Nn, Tab20PVO10avar[[Nn - 2]][[2]] - Tab20PVO10avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; -from here we see that, with increasing Nt, at Nt=7 we get for the first time the reduced local sensitiivity  *)

  Tab21D0PVO10avar(*=Table[{NNt,Re[aaPVD0fast[2, 1][NNt]]},{NNt,3,12}] ={{3,0.221982}, {4,0.221234}, {5,0.22119}, {6,0.221072}, {7,0.221122}, {8,0.220991}, {9,0.221156}, {10,0.220886}, {11,0.221297}, {12,0.220876}}*)={{3,0.2219816615991318`}, {4,0.22123398556691598`}, {5,0.22119013816511046`}, {6,0.22107163962650178`}, {7,0.22112238981469343`}, {8,0.2209910008208125`}, {9,0.22115588074965425`}, {10,0.22088647721260582`}, {11,0.2212970740896767`}, {12,0.22087625736225389`}};
  Tab21PVO10avar(*=Table[{NNt,Re[aaPV21[NNt,O6PVO10[NNt],O8PVO10[NNt]]]},{NNt,3,12}]={{3, 0.200988}, {4, 0.205348}, {5, 0.205571}, {6, 0.205875}, {7, 0.205818}, {8, 0.205888}, {9, 0.205834}, {10, 0.205914}, {11, 0.205766}, {12, 0.20596}}*)={{3, 0.20098810073840523`}, {4, 0.20534757655767047`}, {5, 0.20557120156323372`}, {6, 0.20587496292587354`}, {7, 0.20581829318941888`}, {8, 0.20588796060800008`}, {9, 0.20583400002962315`}, {10, 0.20591447401046092`}, {11, 0.20576624882317118`}, {12, 0.20596023248559753`}}; (* Various local minima of a21 when Nt increases; the first local minimum at Nt=7; then at Nt=9, 11; see also: delTab21PVO10avar= Table[{Nn, Tab21PVO10avar[[Nn - 2]][[2]] - Tab21PVO10avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; -from here we see that, with increasing Nt, at Nt=7 we get for the first time the reduced local sensitiivity  *)

  Tab22D0PVO10avar(*=Table[{NNt,Re[aaPVD0fast[2, 2][NNt]]},{NNt,3,12}] ={{3,0.197853}, {4,0.197923}, {5,0.197935}, {6,0.197966}, {7,0.197956}, {8,0.197979}, {9,0.197955}, {10,0.197982}, {11,0.197973}, {12,0.197938}}*)={{3,0.19785259481867953`}, {4,0.19792305128994148`}, {5,0.19793458580109674`}, {6,0.1979661834247426`}, {7,0.19795569507476365`}, {8,0.19797906100587784`}, {9,0.19795500186328935`}, {10,0.19798213537926626`}, {11,0.19797270112405713`}, {12,0.19793835146193617`}};
  Tab22PVO10avar(*=Table[{NNt,Re[aaPV22[NNt,O6PVO10[NNt],O8PVO10[NNt],O10PVO10[NNt]]]},{NNt,3,12}]={{3, 0.191686}, {4, 0.194975}, {5, 0.195151}, {6, 0.19539}, {7, 0.195346}, {8, 0.195398}, {9, 0.19536}, {10, 0.195414}, {11, 0.19532}, {12, 0.195441}}*)={{3, 0.19168566382293484`}, {4, 0.19497532076485732`}, {5, 0.19515072383166085`}, {6, 0.19538988737990767`}, {7, 0.19534568540520503`}, {8, 0.19539848439118498`}, {9, 0.19536028971546562`}, {10, 0.19541352106162693`}, {11, 0.1953195795610797`}, {12, 0.19544142091167221`}}; (* Various local minima of a22 when Nt increases; the first local minimum at Nt=7; then at Nt=9, 11; see also: delTab22PVO10avar= Table[{Nn, Tab22PVO10avar[[Nn - 2]][[2]] - Tab22PVO10avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; -from here we see that, with increasing Nt, at Nt=7 we get for the first time the reduced local sensitiivity  *)
 ff20exp=Table[{NNt,aexpBin77[2,0]},{NNt,3,12}]; ff20plexp=Table[{NNt,aexpBin77[2,0]+delaexpBin77[2,0]},{NNt,3,12}]; ff20miexp=Table[{NNt,aexpBin77[2,0]-delaexpBin77[2,0]},{NNt,3,12}];
ff21exp=Table[{NNt,aexpBin77[2,1]},{NNt,3,12}];  ff21plexp=Table[{NNt,aexpBin77[2,1]+delaexpBin77[2,1]},{NNt,3,12}]; ff21miexp=Table[{NNt,aexpBin77[2,1]-delaexpBin77[2,1]},{NNt,3,12}]; 
ff22exp=Table[{NNt,aexpBin77[2,2]},{NNt,3,12}];  ff22plexp=Table[{NNt,aexpBin77[2,2]+delaexpBin77[2,2]},{NNt,3,12}]; ff22miexp=Table[{NNt,aexpBin77[2,2]-delaexpBin77[2,2]},{NNt,3,12}];
  (* 19 Sept. 2022 *)
  *)
 
 (* NNt=Nt=7 will be our preferred truncation index because this is the smallest Nt where stationarity is achieved; 19 Sept. 2022 *)
 (* OK *)

(*  Note: aexpBin77[2, 0] + -delaexpBin77[2, 0] = 0.24038 + -0.00261;
          aexpBin77[2, 1] + -delaexpBin77[2, 1] = 0.20610 + -0.00260; 
          aexpBin77[2, 2] + -delaexpBin77[2, 2] = 0.19603 + -0.00282; 	                  aexpBin77[2, 3] + -delaexpBin77[2, 3] = 0.18965 + -0.00313; *)
(* 19 Sept. 2022 *)
(* at Nt=7: a20={7, 0.240086}; a21={7, 0.205666}; a22={7, 0.195227}; *)
	  (*-------*)   
(*---------------------------------------*)
 (* the central value of the PV global fit result, at Nt=7, is ('p3ka1ka2' & 'O6ka2', and truncation at O10): *)
 (* aAL=aALfPV[7]=0.1060905; O10fitPPVMinnum[7]={0.000011740745332238454`, {O4fitf -> -0.0029547054662461976`, O6fitf -> 0.005398958409024923`, O8fitf -> -0.0018134852463415506`, O10fitf -> 0.0007297809837811238`}}; alpha_s(mtau^2)(=alsmtauPVO10[7])=0.321638; (to be compared with the central FO case: alsmtauFOO10[8]=0.320192; the two values are close to each other) *)
 (* 19 Sept. 2022 *)
 
 (*------------------------*) 
 (* Now the estimate of EXPERIMENTAL uncertainty in this global approach with double-pinched PV Borel-Laplace (with Nt=5): *)
 (* See Appendix  PRD84, 113006 (Boito et al., 2011): *)
 
(* Now the M^2_{al}: all Nd=9 values (three rays in the complex M^2-plane) as one array, as M2al[jj] (jj=1,...,9): *)
 (* Nm=NchP+1; Nd=3 Nm;  
M2al0[ii_]:=0.90+(ii-1)*0.60/(Nm-1);
Do[M2al[jj]=M2al0[jj],{jj,1,Nm}]; Do[M2al[jj+Nm]=M2al0[jj] Exp[I Pi/6],{jj,1,Nm}]; Do[M2al[jj+2 Nm]=M2al0[jj] Exp[I Pi/4],{jj,1,Nm}];
  *)
 
 N0=7; epp1=0.001; (* aALfPV[N0]=0.1060905; *)
 (* First derivatives of ReBPFONum[M2] with respect to p1=aAL, p2=O4, p3=O6, p4=O8, p5=O10  [note, we must interpret in this program: O2n = 2 <O_{2 n}>_V]: *)
 (* For variation of ReBPPVNumO10[M2al[jj], N0, O4PV[N0],O6PV[N0],O8PV[N0],O10PV[N0]] when aAL is varied 0.1% up and one 0.1% down: *)
 (* aAL=aALfPV[N0](1+epp1)=0.1061965905: run the program for PV case with this aAL (and kappa=1)  (between 'Double-pinched' and at least until 'Now search for')  and evaluate: *)
 (* aAL=aALfPV[N0](1+epp1): tabReBPPVUP= Table[{M2al[jj], ReBPPVNumO10[M2al[jj], N0, 0, O6PVO10[N0], 0, 0]}, {jj, 1,Nd}] *) tabReBPPVUP={{0.9`, 0.8880902909859436`}, {1.2`, 0.7091878366336621`}, {1.5`, 0.5959843660618702`}, {0.7794228634059949` + 0.44999999999999996` I, 0.8048992752276432`}, {1.0392304845413265` + 0.5999999999999999` I, 0.6560088465812144`}, {1.299038105676658` + 0.7499999999999999` I, 0.5535250423048319`}, {0.6363961030678928` + 0.6363961030678927` I, 0.7295476843443451`}, {0.848528137423857` + 0.8485281374238569` I, 0.6013877893999541`}, {1.0606601717798214` + 1.0606601717798212` I, 0.5058555051115319`}}; 
 (*  aAL=aALfPV[N0](1-epp1)=0.1059844095:  run the program for PV case with this aAL (and kappa=1)  (between 'Double-pinched' and 'Now search for') and evaluate: *)
 (* aAL=aALfPV[N0](1-epp1): tabReBPPVDO= Table[{M2al[jj], ReBPPVNumO10[M2al[jj], N0, 0, O6PVO10[N0], 0, 0]}, {jj, 1,Nd}] *) tabReBPPVDO={{0.9`, 0.8872290825376388`}, {1.2`, 0.7084518245732581`}, {1.5`, 0.5953592457151683`}, {0.7794228634059949` + 0.44999999999999996` I, 0.8039342710079846`}, {1.0392304845413265` + 0.5999999999999999` I, 0.6552735322280671`}, {1.299038105676658` + 0.7499999999999999` I, 0.5529296236942612`}, {0.6363961030678928` + 0.6363961030678927` I, 0.7285777226505995`}, {0.848528137423857` + 0.8485281374238569` I, 0.6007007117905349`}, {1.0606601717798214` + 1.0606601717798212` I, 0.5053193797171874`}}; 
 (* tabdReBPPVd[1] below is now the table of numerical derivative of  ReBPPVNumO10[M2] with respect to aAL around aAL=aALfPV[N0]=0.1060905: *)
 tabdReBPPVd[1] = Table[{M2al[ii],  (tabReBPPVUP[[ii]][[2]]-tabReBPPVDO[[ii]][[2]])/(2 epp1 aALfPV[N0])}, {ii, 1, Nd}];

(* Now the derivatives with respect to the condensates: *)
 dReBPPVd[2][M2_]:=2 Pi^2 Re[ (1/M2^2) ( 1 + 2 M2/sigmax )];
dReBPPVd[3][M2_]:=(6 Pi^2/sigmax^3) (1/2/Pi) Re[NIntegrate[ ((aALfPVO10[N0]/(1 + I be0 aALfPVO10[N0] fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2/sigmax) (1 + Exp[I fi]) + 2 (M2/sigmax)^2) Exp[(sigmax/M2) Exp[I fi]] - 2 (M2/sigmax)^2 Exp[- sigmax/M2] ),{fi, -Pi, Pi}]]; 
dReBPPVd[4][M2_]:=2 Pi^2  (1/3!) Re[(1/M2^4) (1+ 6 (M2/sigmax) + 6 (M2/sigmax)^2)]; 
dReBPPVd[5][M2_]:=2 Pi^2 (1/4!) Re[(1/M2^5) (1+ 8 (M2/sigmax) + 12 (M2/sigmax)^2)];
 (* Now Tables of derivatives with respect to the condensates, for all Mal[jj]: *)  
 tabdReBPPVd[2] = Table[{M2al[ii], dReBPPVd[2][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPPVd[3] = Table[{M2al[ii], dReBPPVd[3][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPPVd[4] = Table[{M2al[ii], dReBPPVd[4][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPPVd[5] = Table[{M2al[ii], dReBPPVd[5][M2al[ii]]}, {ii, 1, Nd}]; 
 tabCovReBPexpBinAL = Table[{M2al[ii],M2al[jj],CovReBPexpBinAL[M2al[ii],M2al[jj]]},{ii,1,Nd},{jj,1,Nd}];

  Np=5; (* This is the number of parameters that we use, in our OPE O10 approach (p1=aAL, p2=O4, p3=O6, p4=O8, p5=O10) *)
 (* Now the NpxNp mtrx A: *)
 (* Now calculate the NpxNp mtrx AA(=A), for the PV case, using tables obtained above: *)
 tabAAPV=Table[{ik, jk, Sum[(tabdReBPPVd[ik][[jj]][[2]] tabdReBPPVd[jk][[jj]][[2]])/tabCovReBPexpBinAL[[jj]][[jj]][[3]], {jj, 1, Nd}]}, {ik, 1, Np}, {jk, 1, Np}];
 Do[AAPVa[ii,jj]=tabAAPV[[ii]][[jj]][[3]],{ii,1,Np},{jj,1,Np}];  (* This is NpxNp matrix A (or: AA), calculated more efficiently *)
 AAPVm=Table[AAPVa[ik,jk],{ik,1,Np},{jk,1,Np}];
(* AAPVm[[ik,jk]] = AA[ik,jk]; AAPVm is the NpxNp matrix *)
AAPVmI=Inverse[AAPVm];
 (* beware if message:  "badly condition mtrx," in such case  one of Det[AAPVm] and Det[AAPVmI] is extremely close to zero and the other is extremely large *)

 (* Now the NpxNp covariance matrix of the fitted parameters pj: *) 
 (* p1=aAL, p2=<O4>, p3=<O6>, p4=<O8>, p5=<O10>; where O_{2n} means: O_{2n}=2 <O_{2n}>_V  *)
(* PP[i,j]= < (delta p_i) (delta p_j) >: *)
  (* We calulate PP[i,j] covariance matrix of the extracted parameters, using the tables obtained above: *) 
 PPtabPV[jp_,kp_]:=Sum[AAPVmI[[jp,jpp]] AAPVmI[[kp,kpp]] tabdReBPPVd[jpp][[Al]][[2]] tabdReBPPVd[kpp][[Be]][[2]] tabCovReBPexpBinAL[[Al]][[Be]][[3]]/(tabCovReBPexpBinAL[[Al]][[Al]][[3]]  tabCovReBPexpBinAL[[Be]][[Be]][[3]]),{jpp,1,Np},{kpp,1,Np},{Be,1,Nd},{Al,1,Nd}];
 (* epp1=0.001: Table[{Nn, Sqrt[PPtabPV[Nn, Nn]]}, {Nn, 1, Np}]= {{1, 0.000920796}, {2, 0.000310954}, {3, 0.00067329}, {4, 0.000364165}, {5, 0.000188452}}; *)
 (* The experimental uncertainty of aAL=a(sigmax) is: del(aAL)=Sqrt[PPtabPV[1,1]]; the experimental uncertainty for O4 is del(O4)=Sqrt[PPtabPV[2,2]]; etc. *)
 (* --> a(sigmax)=aALfPVO10[N0] +- Sqrt[PPtabPV[1, 1]]=0.106091+-0.000920796_{exp}; --> alpha_s(mtau^2)=0.321638+-0.002667_{exp};  20 Sept. 2022 *)
 (*-------------------------*)
 (* Summary of the experimental uncertainties, 'O10' PV approach with Nt=N0=7*)
 (* 'p3ka1ka2'&'O6ka2' O10 PV approach with Nt=N0(=7): *)
 (* alpha_s(mtau^2)=0.321638+-0.002667_{exp} 
 10^3 O4= -3.0 +- 0.3_{exp};
 10^3 O6= +5.4 +- 0.7_{exp};
 10^3 O8= -1.8 +-0.4_{exp};
  10^3 O10= +0.7 +- 0.2_{exp}; *)
 (* 20 Sept. 2022; checked on 26 Sept. 2022 *)
   (*---------------------------------------------*)
(* Variation of kappa, PV global: *)
(* kappa=2.0: 
 (* Use prg. ScanSumRMSbALEPHPV.m, with Nt0=7; kappa=2 and d4input=275; [note: aAL=a(kappa*sigmax)]: *)
 aAL=0.087991; O10fitPPVMinN[N0]={0.0000126728, {O4fitf -> -0.0030882, O6fitf -> 0.00523896,  O8fitf -> -0.0018595, O10fitf -> 0.000724899}}; alpha_s(mtau^2)=0.317997; a(sigmax)=0.105979;  
kappa=1.0 (the central case of the global fit): 
 aAL=0.1060905; O10fitPPVMinnum[N0]={0.000011740745332238454`, {O4fitf -> -0.0029547054662461976`, O6fitf -> 0.005398958409024923`, O8fitf -> -0.0018134852463415506`, O10fitf -> 0.0007297809837811238`}}; a(sigmax)=0.1060905; alpha_s(mtau^2)(=alsmtauPVO10[7])=0.321638; 
 kappa=0.50:
 aAL=0.136007; O10fitPPVMinnum[N0]={8.50278*10^-6, {O4fitf -> -0.00289374, O6fitf -> 0.00572602, O8fitf -> -0.00183141, O10fitf -> 0.000923249}}; a(sigmax)=0.106161;
 *)
 (* 20 Sept. 2022; checked on 26 Sept. 2022 ---> *)
 (* kap=1^{+1}-{1/2}:*)
 (*
alpha_s(mtau^2)=0.321638+{^{-0.000324}_{+0.000205}}_{kap};
10^3 O4= -3.0+{^{-0.1}_{+0.1}}_{kap};
10^3 O6= +5.4+{^{-0.2}_{+0.3}}_{kap};
10^3 O8= -1.8+{^{-0.1}_{+0.0}}_{kap};
  10^3 O10=+0.7+{^{-0.0}_{+0.2}}_{kap};
  *)  
  (* 20 Sept. 2022; checked on 25 Sept. 2022 *)
  (*----------------*)
 (* Variation of d4, PV global (d4=227 +- 63; kappa=1): 
  When ddtk[4,tkap]=227+63=338:
  aAL=0.102830; O10fitPPVMinnum[7]={0.0000278475, {O4fitf -> -0.00229926, O6fitf -> 0.00823351, O8fitf -> -0.00300158, O10fitf -> 0.0010119}}; alpha_s(mtau^2)=0.312183;
  Central case: ddtk[4,tkap]=275.:
aAL=0.1060905; O10fitPPVMinnum[7]={0.000011740745332238454`, {O4fitf -> -0.0029547054662461976`, O6fitf -> 0.005398958409024923`, O8fitf -> -0.0018134852463415506`, O10fitf -> 0.0007297809837811238`}}; a(sigmax)=0.1060905; alpha_s(mtau^2)(=alsmtauPVO10[7])=0.321638;
  When ddtk[4,tkap]=227-63=212:
  aAL=0.107090; O10fitPPVMinnum[7]={0.0000111401, {O4fitf -> -0.00281572, O6fitf -> 0.00488116, O8fitf -> -0.00176815, O10fitf -> 0.000656549}}; alpha_s(mtau^2)=0.324529;
  (* checked on 26 Sept. 2022 *)  
  *)
 (*
  alpha_s(mtau^2)=0.321638+{^{-0.009455}_{+0.002891}}_{d4};
10^3 O4= -3.0+{^{+0.7}_{+0.2}}_{d4};
10^3 O6= +5.4+{^{+2.8}_{-0.5}}_{d4};
10^3 O8= -1.8+{^{-1.2}_{+0.0}}_{d4};
10^3 O10=+0.7+{^{+0.3}_{+0.0}}_{d4};
  *)
  (* 20 Sept. 2022;  checked on 26 Sept. 2022 *)  
  (*---------------------*)
  (* Variation of Nt (truncation order, Nt=7+-1, 7+-2: *)
   (*
   aALfPV[9]=0.106087; O10fitPPVMinnum[9]={0.0000115206, {O4fitf -> -0.00295997, O6fitf -> 0.00540952, O8fitf -> -0.00181536, O10fitf -> 0.000741354}}; alpha_s(mtau^2)=alsmtauPVO10[9]=0.321628;
   aALfPV[8]=0.106073; O10fitPPVMinnum[8]={0.0000119291, {O4fitf -> -0.00296935, O6fitf -> 0.00536113, O8fitf -> -0.00181829, O10fitf -> 0.000713343}}; alpha_s(mtau^2)=alsmtauPVO10[8]=0.321587;
   (* The central case:  *)
   aALfPV[7]=0.1060905; O10fitPPVMinnum[7]={0.0000117407, {O4fitf -> -0.00295471, O6fitf -> 0.00539896, O8fitf -> -0.00181349, O10fitf -> 0.000729781}}; alpha_s(mtau^2)=alsmtauPVO10[7]=0.321638;
   aALfPV[6]=0.106077; O10fitPPVMinnum[6]={0.0000117614, {O4fitf -> -0.002964, O6fitf -> 0.00538456, O8fitf -> -0.00181602, O10fitf -> 0.000725152}}; alpha_s(mtau^2)=alsmtauPVO10[6]=0.321597;   
   aALfPV[5]=0.106151; O10fitPPVMinnum[5]={0.00001189, {O4fitf -> -0.00293118, O6fitf -> 0.00541762, O8fitf -> -0.00180824, O10fitf -> 0.000732929}}; alpha_s(mtau^2)=alsmtauPVO10[5]=;0.321813
    *) 
  (* Nt=7^{+2}_{-2}:*)
 (*
alpha_s(mtau^2)=0.321638+{^{-0.000051}_{+0.000175}}_{Nt};
10^3 O4= -3.0+{^{-0.0}_{+0.1}}_{Nt};
10^3 O6= +5.4+-0.0_{Nt};
10^3 O8= -1.8+-0.0_{Nt};
10^3 O10=+0.7+-0.0_{Nt};
  *)
  (* 20 Sept. 2022; checked on 26 Sept. 2022 *)
  (*-------------------------*)
  (* Variation due to the renormalon ambiguity ('amb') to Adler: in ReBPPVNumO10 we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn]+BPPVdelD0[M2f]) *)
  (*
   aAL=0.1060705; O10fitPPVMinnum[7]={0.0000117887, {O4fitf -> -0.00382497, O6fitf -> 0.00539665, O8fitf -> -0.00182161, O10fitf -> 0.000756584}}; alpha_s(mtau^2)=0.321580; *)
(*---*)
 (* The central value  (BBPVD0[M2f,Nn]): 
  aAL=aALfPV[7]=0.1060905; O10fitPPVMinnum[7]={0.000011740745332238454`, {O4fitf -> -0.0029547054662461976`, O6fitf -> 0.005398958409024923`, O8fitf -> -0.0018134852463415506`, O10fitf -> 0.0007297809837811238`}}; alpha_s(mtau^2)(=alsmtauPVO10[7])=0.321638; *)
  (*---*)
  (* When  in ReBPPVNumO10 we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn]-BPPVdelD0[M2f]) *)
  (*  aAL=0.1061115; O10fitPPVMinnum[7]={0.0000117157, {O4fitf -> -0.00208282, O6fitf -> 0.00540312, O8fitf -> -0.00180667, O10fitf -> 0.000703613}}; alpha_s(mtau^2)=0.321699; *)
(*
alpha_s(mtau^2)=0.321638+{^{-0.000058}_{+0.000061}}_{amb};
 10^3 O4= -3.0+{^{-0.8}_{+0.9}}_{amb};
 10^3 O6= +5.4+-0.0_{amb};
 10^3 O8= -1.8+-0.0_{amb};
 10^3 O10=+0.7+-0.0_{amb};
 *)
 (* 20 Sept. 2022; checked on 26 Sept. 2022 *)
 (*--------------*)
 (* Variation when going beyond the O10 truncation, including O12 in OPE: *)
 (* N0=7; Using d4input = 275; *)
  (* N0=7; aAL=0.1066415; O12fitPPVMinNum[N0]={1.53425*10^-7, {O4fitf -> -0.00325579, O6fitf -> 0.00643967, O8fitf -> -0.0028478, O10fitf -> 0.00196693, O12fitf -> -0.000857502}}; alpha_s(mtau^2)=0.323232 *)
 (* The central case:
  N0=7; aALfPVO10[N0]=0.106091; O10fitPPVMinnum[N0]={0.0000117407, {O4fitf -> -0.00295471, O6fitf -> 0.00539896, O8fitf -> -0.00181349, O10fitf -> 0.000729781}}; alpha_s(mtau^2)=alsmtauPVO10[N0]=0.321638; *)

 (* --> alpha_s(mtau^2)-alsmtauPVO10[N0]=+0.001594; *)
 (* 10^3 (O4fitf-O4PVO10[N0])=-0.30 =app -0.3 *)
 (* 10^3 (O6fitf-O6PVO10[N0])= +1.04 =app +1.0 *)
 (* 10^3 (O8fitf-O8PVO10[N0])= -1.03 =app -1.0 *)
 (* 10^3 (O10fitf-O10PVO10[N0])= +1.24 =app +1.1 *)
 (* -->
  alpha_s(mtau^2)=0.321638+-0.001594_{O12};
  10^3 O4= -3.0 -+0.3_{O12}; 
  10^2 O6= +5.4 +-1.0_{O12};
  10^2 O8= -1.8 -+1.0_{O12};
  10^3 O10=+0.7 +-1.1_{O12};
*)
  (* 27 Sept. 2022 *)
 (*--------------*)
 (* Variation when increasing the number of data bins from  NNAL=77 to NNAL=79 (i.e., including the last two bins, which have bad data): *)
 (* N0=7; NNAL=79: aAL=0.10340;  (* aAL=a(sigmax); note, aAL determined only up to 5 digits; shallow minimum of chi^2 *); O10fitPPVMinnum[N0]={0.000057653, {O4fitf -> -0.00342623, O6fitf -> 0.00659414, O8fitf -> -0.00254008, O10fitf -> 0.00114513}}; alpha_s(mtau^2)=0.324836; (note that sigmax=3.15708, because NNAL=79) *)
  (* The central case:
  N0=7; aALfPVO10[N0]=0.106091; O10fitPPVMinnum[N0]={0.0000117407, {O4fitf -> -0.00295471, O6fitf -> 0.00539896, O8fitf -> -0.00181349, O10fitf -> 0.000729781}}; alpha_s(mtau^2)=alsmtauPVO10[N0]=0.321638; *)
 (* --> alpha_s(mtau^2)-alsmtauPVO10[N0]=+0.003198; *)
 (* 10^3 (O4fitf-O4PVO10[N0])= -0.47 =app -0.5 *)
 (* 10^3 (O6fitf-O6PVO10[N0])= +1.20 =app +1.2 *)
 (* 10^3 (O8fitf-O8PVO10[N0])= -0.73 =app -0.7 *)
 (* 10^3 (O10fitf-O10PVO10[N0])= +0.42 =app +0.4 *)
  
 (* -->
  alpha_s(mtau^2)=0.321638+-0.003198_{Nbin}
  10^3 O4= -3.0 -+0.5_{Nbin}; 
  10^2 O6= +5.4 +-1.2_{Nbin};
  10^2 O8= -1.8 -+0.7_{Nbin};
  10^3 O10=+0.7 +-0.4_{Nbin};
*) 
 (* 27 Sept. 2022 *)
  (*-----------------------------------------*)
  (* Altogether, the PV(+CI) global fit approach with (O61, O62) [and a(Q^2) at O61 and O62 runs according to the one-loop RGE from a(sigmax)=aALfPV[Nt] to a(Q^2=sigmax Exp[I fi]) ], gives: *)
  (* PV(+CI)
   Nt(=N0)=7;  ( ddtk[4, Exp[tKK]]=275.);
   d4=275.+-63.; kap=1^{+1}_{-1/2};  Nt=7^{+2}_{-2}:

   alpha_s(mtau^2)=0.321638+-0.002667_{exp}+{^{-0.000323}_{+0.000205}}_{kap}+{^{-0.009455}_{+0.002891}}_{d4}+{^{-0.000051}_{+0.000175}}_{Nt}+{^{-0.000058}_{+0.000061}}_{amb}+-0.001594_{O12}+-0.003198_{Nbin} = 0.321638+{^{+0.005321}_{-0.010459}};
 =  0.3216+{^{+0.0054}_{-0.0104}} =app 0.322+{^{+0.005}_{-0.011}}
   10^3 O4= -3.0 +- 0.3_{exp}+{^{-0.1}_{+0.1}}_{kap}+{^{+0.7}_{+0.2}}_{d4}+{^{-0.0}_{+0.1}}_{Nt}+{^{-0.8}_{+0.9}}_{amb}-+0.3_{O12}-+0.5_{Nbin} =  -3.0+{^{+1.3}_{-1.1}};
   10^3 O6= +5.4 +- 0.7_{exp}+{^{-0.2}_{+0.3}}_{kap}+{^{+2.8}_{-0.5}}_{d4}+-0.0_{Nt}+-0.0_{amb}+-1.0_{O12}+-1.2_{Nbin} =  +5.4+{^{+3.3}_{-1.8}};
   10^3 O8= -1.8 +- 0.4_{exp}+{^{-0.1}_{+0.0}}_{kap}+{^{-1.2}_{+0.0}}_{d4}+-0.0_{Nt}+-0.0_{amb}-+1.0_{O12}-+0.7_{Nbin} = -1.8+{^{+1.3}_{-1.8}};
   10^3 O10=+0.7 +- 0.2_{exp}+{^{-0.0}_{+0.2}}_{kap}+{^{+0.3}_{+0.0}}_{d4}+-0.0_{Nt}+-0.0_{amb}+-1.1_{O12}+-0.4_{Nbin} = +0.7+{^{+1.2}_{-1.2}};
   chi^2=0.00001174;
   *)
   
   (* 24 Sept. 2022; alpha_s checked on 26 Sept. 2022; added _{O12} and _{Nbin} uncertainties on 27 Sept. 2022 *)
			  
 (* End of PV(+CI) global fit results with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax)=aALftFO[Nt] to a(Q^2=sigmax Exp[I fi]) ] *)
(*------------------------------------*)
(*------------------------------------*)
(*------------------------------------*)
   (* Some Figures: *)
   a20avar= ListPlot[{Tab20FOO10avar, Tab20PVO10avar, Tab20D0FOO10avar, Tab20D0PVO10avar, ff20exp, ff20miexp, ff20plexp},PlotRange -> {{2.8, 12.2}, {0.22, 0.30}}, PlotStyle -> {{Red, PointSize[0.001]},{Blue, PointSize[0.001]},Red, Blue, {Black, PointSize[0.0005]}, {Black, PointSize[0.0005]}, {Black, PointSize[0.0005]}}, Filling -> {6 -> {7}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True, True, True}, PlotMarkers -> {\[FilledSmallSquare], \[FilledSmallCircle],\[FilledSquare], \[FilledCircle],\[Bullet], \[Bullet],\[Bullet]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.90, 0.95}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.07, 0.95}]], Inset[Style["FO (D=0)", 11], Scaled[{0.47, 0.95}]], Inset[Style["PV (D=0)", 11], Scaled[{0.75, 0.75}]](* , Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.77, 0.60}]]*) },  PlotLegends -> Placed[{"FO", "PV"}, {0.90, 0.45}]];
a20avarZoom= 
 ListPlot[{Tab20FOO10avar, Tab20PVO10avar, ff20exp}, 
  PlotRange -> {{2.8, 12.2}, {0.2399, 0.2405}}, 
  PlotStyle -> {Red, Blue, {Black, PointSize[0.001]}}, 
  AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True}, 
  PlotMarkers -> {\[FilledSquare], \[FilledCircle], \[Bullet], \[Bullet],\[Bullet]}, 
  Epilog -> {Inset[Style["(b)", 13], Scaled[{0.30, 0.95}]],
	   Inset[Style["(V+A)", 13], Scaled[{0.65, 0.95}]],
    Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.07, 0.45}]]
  (* ,Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.40, 0.95}]] *)}, PlotLegends -> Placed[{"FO", "PV"}, {0.90, 0.15}]];
 
 a21avar= ListPlot[{Tab21FOO10avar, Tab21PVO10avar, Tab21D0FOO10avar, Tab21D0PVO10avar, ff21exp, ff21miexp, ff21plexp},PlotRange -> {{2.8, 12.2}, {0.20, 0.23}}, PlotStyle -> {{Red, PointSize[0.001]},{Blue, PointSize[0.001]},Red, Blue, {Black, PointSize[0.0005]}, {Black, PointSize[0.0005]}, {Black, PointSize[0.0005]}}, Filling -> {6 -> {7}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True, True, True}, PlotMarkers -> {\[FilledSmallSquare], \[FilledSmallCircle],\[FilledSquare], \[FilledCircle],\[Bullet], \[Bullet],\[Bullet]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.40, 0.95}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.07, 0.95}]], Inset[Style["FO (D=0)", 11], Scaled[{0.70, 0.88}]], Inset[Style["PV (D=0)", 11], Scaled[{0.15, 0.75}]](* , Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.77, 0.60}]]*) },  PlotLegends -> Placed[{"FO", "PV"}, {0.90, 0.45}]];
a21avarZoom= 
 ListPlot[{Tab21FOO10avar, Tab21PVO10avar, ff21exp}, 
  PlotRange -> {{2.8, 12.2}, {0.2055, 0.2063}}, 
  PlotStyle -> {Red, Blue, {Black, PointSize[0.001]}}, 
  AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True}, 
  PlotMarkers -> {\[FilledSquare], \[FilledCircle], \[Bullet], \[Bullet],\[Bullet]},  Epilog -> {Inset[Style["(b)", 13], Scaled[{0.30, 0.95}]],
	     Inset[Style["(V+A)", 13], Scaled[{0.65, 0.95}]],
    Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.07, 0.45}]]
  (* ,Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.40, 0.95}]] *)},  PlotLegends -> Placed[{"FO", "PV"}, {0.90, 0.15}]];

a22avar= ListPlot[{Tab22FOO10avar, Tab22PVO10avar, Tab22D0FOO10avar, Tab22D0PVO10avar, ff22exp, ff22miexp, ff22plexp},PlotRange -> {{2.8, 12.2}, {0.19, 0.20}}, PlotStyle -> {{Red, PointSize[0.001]},{Blue, PointSize[0.001]},Red, Blue, {Black, PointSize[0.0005]}, {Black, PointSize[0.0005]}, {Black, PointSize[0.0005]}}, Filling -> {6 -> {7}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True, True, True}, PlotMarkers -> {\[FilledSmallSquare], \[FilledSmallCircle],\[FilledSquare], \[FilledCircle],\[Bullet], \[Bullet],\[Bullet]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.40, 0.95}]], Inset[Style[Superscript[(a), "(2,2)"], 13], Scaled[{0.07, 0.95}]], Inset[Style["FO (D=0)", 11], Scaled[{0.80, 0.94}]], Inset[Style["PV (D=0)", 11], Scaled[{0.65, 0.75}]](* , Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.77, 0.60}]]*) },  PlotLegends -> Placed[{"FO", "PV"}, {0.90, 0.15}]];
a22avarZoom= 
 ListPlot[{Tab22FOO10avar, Tab22PVO10avar, ff22exp}, 
  PlotRange -> {{2.8, 12.2}, {0.1950, 0.1962}}, 
  PlotStyle -> {Red, Blue, {Black, PointSize[0.001]}}, 
  AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True}, 
  PlotMarkers -> {\[FilledSquare], \[FilledCircle], \[Bullet], \[Bullet],\[Bullet]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.30, 0.95}]],
	     Inset[Style["(V+A)", 13], Scaled[{0.65, 0.95}]], 
    Inset[Style[Superscript[(a), "(2,2)"], 13], Scaled[{0.07, 0.45}]]
(* ,Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.40, 0.95}]] *)},  PlotLegends -> Placed[{"FO", "PV"}, {0.90, 0.65}]];

(*   Export["a20avarVAa.pdf", a20avar]; Export["a20avarVAb.pdf",a20avarZoom]; Export["a21avarVAa.pdf", a21avar]; Export["a21avarVAb.pdf",a21avarZoom]; Export["a22avarVAa.pdf", a22avar]; Export["a22avarVAb.pdf",a22avarZoom]; *) 

aAL=aALfFO[8](=0.104654), [ alpha_s(mtau^2)=alsmtauFO[8]=0.317477,i.e., the coupling alpha_s taken FIXED at value as obtained from global BL FO fits at Nt=8: *)(*------------------------------------*)
  (* When alpha_s(mtau^2) is FIXED ('afxd'), Tables of a20(D=0), a21(D=0) vs Nt:
   i.e., alpha_s(mtau^2) is the average of all four approaches (FO and PV for (V+A)-channel and for V-channel [for that channel, see prg. SumRMSbCombVp3ka1ka2O6ka2.m]): alpha_s=0.3169:

 ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}]; (* OK *)
ff22exp=Table[{NNt,aexpBin77[2,2]},{NNt,3,12}];  ff22plexp=Table[{NNt,aexpBin77[2,2]+delaexpBin77[2,2]},{NNt,3,12}]; ff22miexp=Table[{NNt,aexpBin77[2,2]-delaexpBin77[2,2]},{NNt,3,12}];
   
   alsmtC = 0.3169;
  aAL = aMSext[sigmax/mtau^2, alsmtC/Pi];
 Tab20D0FOO10afxd= Table[{NN, Re[aaFOD0[2, 0, aAL][NN]]}, {NN, 3, 12}];
 Tab21D0FOO10afxd= Table[{NN, Re[aaFOD0[2, 1, aAL][NN]]}, {NN, 3, 12}];
 Tab22D0FOO10afxd= Table[{NN, Re[aaFOD0[2, 2, aAL][NN]]}, {NN, 3, 12}]; 

   kappa=1;
   alsmtC = 0.3169; aAL=aMSext[sigmax/mtau^2, alsmtC/Pi];
   d4input=275;
   <<SumRMSbALEPHp3ka1ka2O6ka2PVsh.m;
 aaPVD0fast[2,nn_][Nt_]:=(1/(2 Pi)) ((nn+3)/(nn+1))  NIntegrate[ ( Exp[I th] (1 - (-Exp[I th])^(nn+1)) + ((nn+1)/(nn+3))(1 - (-Exp[I th])^(nn+3)) ) DBorfast[th,Nt], {th, -Pi, Pi}];
Tab20D0PVO10afxd=Table[{NNt,Re[aaPVD0fast[2, 0][NNt]]},{NNt,3,12}];
Tab21D0PVO10afxd=Table[{NNt,Re[aaPVD0fast[2, 1][NNt]]},{NNt,3,12}];
Tab22D0PVO10afxd=Table[{NNt,Re[aaPVD0fast[2, 1][NNt]]},{NNt,3,12}];

a20FOafxd= 
 ListPlot[{Tab20D0FOO10afxd, Tab20D0PVO10afxd, ff20exp, ff20miexp, 
   ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.23, 0.28}}, 
  PlotStyle -> {Red, 
    Blue, {Black, PointSize[0.0005]}, {Black, 
     PointSize[0.0005]}, {Black, PointSize[0.0005]}}, 
  Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, 
  Joined -> {True, True, True, True, True}, 
  PlotMarkers -> {\[FilledSquare], \[FilledCircle], \[Bullet], \
\[Bullet], \[Bullet]}, 
  Epilog -> {Inset[Style["(a)", 13], Scaled[{0.85, 0.95}]], 
    Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.07, 0.95}]],
     Inset[Style[
      Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]] == 
       Round[10000 alsmtC]/10000., 12], Scaled[{0.46, 0.40}]]}, 
  PlotLegends -> Placed[{"FO(D=0)", "PV(D=0)"}, {0.85, 0.45}]];

a21FOafxd= 
 ListPlot[{Tab21D0FOO10afxd, Tab21D0PVO10afxd, ff21exp, ff21miexp, 
   ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.200, 0.220}}, 
  PlotStyle -> {Red, 
    Blue, {Black, PointSize[0.0005]}, {Black, 
     PointSize[0.0005]}, {Black, PointSize[0.0005]}}, 
  Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, 
  Joined -> {True, True, True, True, True}, 
  PlotMarkers -> {\[FilledSquare], \[FilledCircle], \[Bullet], \
\[Bullet], \[Bullet]}, 
  Epilog -> {Inset[Style["(b)", 13], Scaled[{0.85, 0.95}]], 
    Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.07, 0.95}]],
     Inset[Style[
      Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]] == 
       Round[10000 alsmtC]/10000., 12], Scaled[{0.40, 0.50}]]}, 
  PlotLegends -> Placed[{"FO(D=0)", "PV(D=0)"}, {0.85, 0.55}]];



  (* Export["a20afixedO6ka2.pdf", a20FOafxd]; Export["a21afixedO6ka2.pdf", a21FOafxd]; *)
  
(* 6 Oct. 2022 *)
(*------------------------------------*)
(*------------------------------------*)
 (* To relate alpha_s(MZ^2) with alpha_s(mtau^2), via 5-loop MSbar RGE: *)
 (*
mtau=1.77686; mc=1.27; mb=4.20; kappa=2.; Q=mtau;
alphaMSMZ=0.11...; <</home/gorazd/SumRpQCD.dir/RGErundown5l.m;  aaMSQ*Pi
 Then: aaMSQ*Pi is the value of alpha_s(mtau^2) *)
