(* SumRPMSbALEPHM2m090O10.m: as prg SumRPMSbALEPHM2m090.m, *)
(* but with O10 (D=10) terms included in the Borel-Laplace (global) fits *)
(*-------------------------------------------------------------_*)
(* 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-resummation, or renormalon-motivated sum *)
(*-------------------------------------------------*)
 <</home/gorazd/3danQCD.dir/AdlerFunction4lMiniMOM.m;
<</home/gorazd/3danQCD.dir/MSbarRenMod5A.save;

 (*  alphaMSMZ=0.XXXX; *)

 (* Inputs: *)
mtauphys=1.77686;  (*= 1.77686+-0.00012 *)
mtau=mtauphys;

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.65 =< M2min < 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) *)
 sigmax=2.80; 

Save["O4O6fitMinALEPH1179090.save", {M2min,alphaMSMZ,sigmax}]; 
(*--------------------------------__*)
(* ALEPH experimental data: *)
<</home/gorazd/3danQCD.dir/expdataALEPH.m;
(* NNAL=77 will be now, and not NNAL=79 bins *)
NNAL=77;
sigmax=2.80;
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,-1.,+1.}];
(* ' *)
  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 *)

(* ff[2,x]=bta[x]; ff[3,x]=bta[x] bta'[x]; etc. ' *)
ff[nn_, xx_] := 
  Module[{ffp, nnp, xxp}, ffp[2, xxp_] := bta[xxp]; 
   ffp[nnp_, xxp_] := 
    ffp[nnp, xxp] = bta[xxp] D[ffp[nnp - 1, xxp], xxp]; ffp[nn, xx]];

(* ta[1,x_]:=x; ta[2,x]=((-1)/be0 1!) bta[x];  ta[3,x]=((-1)^2/be0^2 2!) bta[x] bta'[x]; ' *)
taprov[nn_,x_]:=((-1)^(nn-1)/(be0^(nn-1) (nn-1)!)) ff[nn,x];  (* nn=2,3,4, *)
taprov[1,x_]:=x;
ta[nn_,x_,NN_]:=Normal[Series[Expand[taprov[nn,x]], {x, 0, NN}]];
  (* i.e., we will expand ta_n only up to powers of a^NN (not higher), *)
  (* because we will not need higher in practice, and we save time now *)
ta[1,x_,NN_]:=x;
  (*----------------------------------------------*)
<</home/gorazd/3danQCD.dir/DiegoT.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 *)
 BPFOD0fast[M2_,Nn_,aALf_]:=Sum[coeffBPFO[M2,nn] x^nn,{nn,1,Nn}]/.{x->aALf};
(* aALf=alpha_s(kappa*sigmax)/Pi *)
 (* And now for tFO, i.e., in {tilde a}_nn(kappa sigmax)=taprov[nn,aALf]: *)
  coeffBPtFO[M2f_,nn_]:=Sum[tddtk[nn-kk,tkap] be0^kk cKP[kk,M2f]  Binomial[nn,kk],{kk,0,nn}];
BPtFOD0[M2_,Nn_,aALf_]:=Sum[coeffBPtFO[M2,nn-1] taprov[nn,x],{nn,1,Nn}]/.{x->aALf};  (* This is double-pinched D=0 part of the Borel-Laplace in the tFO ({tilde FOPT}) approach, truncated at {tilde a}_{Nn}(kappa sigmax) *)
(* 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 data *)

ReBPFOfast[M2f_,NN_,aALf_,O4_,O6_,O8_,O10_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0fast[M2f,NN,aALf]+ 2 Pi^2 ( (1/(M2f sigmax))(2 O4 + O6/sigmax) + (1/M2f^2) (O4 + 2 (O6/sigmax) + (O8/sigmax^2)) + (1/M2f^3) (1/2)(O6 + 2 O8/sigmax + O10/sigmax^2) + (1/M2f^4) (1/6) (O8 +2 O10/sigmax) + (1/M2f^5) (1/24) O10 )];  (* This is the FOPT approach when the order of truncation is a^NN; with a=a(sigmax); the higher order coefficients are those of the renormalon-motivated model in MSbar; the higher-twist terms are as determied by double-pinched Borel-Laplace  *) 

ReBPtFO[M2f_,Nn_,aALf_,O4_,O6_,O8_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPtFOD0[M2f,Nn,aALf]+ 2 Pi^2 ( (1/(M2f sigmax))(2 O4 + O6/sigmax) + (1/M2f^2) (O4 + 2 (O6/sigmax) + (O8/sigmax^2)) + (1/M2f^3) (1/2)(O6 + 2 O8/sigmax + O10/sigmax^2) + (1/M2f^4) (1/6) (O8 +2 O10/sigmax) + (1/M2f^5) (1/24) O10 )];  (* This is the FOPT approach when the order of truncation is a^Nn; with a=a(sigmax); the higher order coefficients are those of the renormalon-motivated model in MSbar; the higher-twist terms are as determined by double-pinched Borel-Laplace  *)

(* Now search for condensate values with this Borel-Laplace (double) pinched ('P') approach, with FO method for D=0 part: *)
NchP=10; (* 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): *)

Clear[aAlf];
O4init=0.001; O6init=-0.001; O8init=0.; O10init=0.; aALfinit=0.100;
alO4O6O8O10fitPFOMinN[NN_]:= FindMinimum[Sum[(ReBPFOfast[M2FP[jj] Exp[I Pi/6], NN, aALf,O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfast[M2FP[jj] Exp[I Pi/4], NN, aALf, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfast[M2FP[jj], NN, aALf, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init}}];
(* When O10 not included in the fit (O10=0): *)
(* alO4O6O8fitPFOMinN[NN_]:= FindMinimum[Sum[(ReBPFOfast[M2FP[jj] Exp[I Pi/6], NN, aALf,O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfast[M2FP[jj] Exp[I Pi/4], NN, aALf, O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfast[M2FP[jj], NN, aALf, O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}]; *)
(* aALf=alpha_s(sigmax)/Pi *)
(* tFO ({tilde FOPT}) approach: *)
alO4O6O8O10fitPtFOMinN[Nn_]:= FindMinimum[Sum[(ReBPtFO[M2FP[jj] Exp[I Pi/6], Nn, aALf,O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPtFO[M2FP[jj] Exp[I Pi/4], Nn, aALf, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPtFO[M2FP[jj], Nn, aALf, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init}}];
(* aALf=alpha_s(kappa*sigmax)/Pi *)

 (* Do[alO4O6O8O10fitPFOMinnum[nn]=alO4O6O8O10fitPFOMinN[nn], {nn, 3, 15}]; *)
(* Do[alO4O6O8O10fitPtFOMinnum[nn]=alO4O6O8O10fitPtFOMinN[nn], {nn, 3, 12}]; *)			     
(* Save["alO4O6O8fitMinALEPH090O10.save", {M2min,sigmax}];  *)
 (* Save["alO4O6O8fitMinALEPH090O10.save", {M2min, sigmax, NchP, alO4O6O8O10fitPFOMinnum}]; *)
 (* Save["alO4O6O8fitMinALEPH090O10.save", {M2min, sigmax, NchP, alO4O6O8O10fitPtFOMinnum}]; *)
 
(* FO global, NchP=10; M2min=0.9, (and kappa=1), aALf VARIED: *)
(*
 kappa=1; 
alO4O6O8O10fitPFOMinnum[3]={0.0000181596, {aALf -> 0.108732, O4fitf -> 0.00136836, O6fitf -> -0.000532766, O8fitf -> 0.000545872, O10fitf -> -0.000315657}}; 
alO4O6O8O10fitPFOMinnum[4]={8.8112*10^-7, {aALf -> 0.105174, O4fitf -> 0.000225725, O6fitf -> 0.000626003, O8fitf -> -0.000366152, O10fitf -> 0.00011218}}; 
alO4O6O8O10fitPFOMinnum[5]={0.0000195842, {aALf -> 0.104093, O4fitf -> -0.000848814, O6fitf -> 0.00157982, O8fitf -> -0.00106692, O10fitf -> 0.000427771}}; 
alO4O6O8O10fitPFOMinnum[6]={0.0000405758, {aALf -> 0.10358, O4fitf -> -0.00182466, O6fitf -> 0.00229841, O8fitf -> -0.00154483, O10fitf -> 0.000630505}}; 
alO4O6O8O10fitPFOMinnum[7]={0.0000543979, {aALf -> 0.103123, O4fitf -> -0.00274054, O6fitf -> 0.00283537, O8fitf -> -0.00186052, O10fitf -> 0.000754875}}; 
alO4O6O8O10fitPFOMinnum[8]={0.0000635937, {aALf -> 0.102593, O4fitf -> -0.0036694, O6fitf -> 0.00326539, O8fitf -> -0.00208509, O10fitf -> 0.00083781}}; 
 (* compare with (global fit when when O10=0): *)
 (* alO4O6O8fitPFOMinnum[8]={0.004287816184596373`, {aALf -> 0.10122255774076128`, O4fitf -> -0.0027756129809026844`, O6fitf -> 0.002062451844738742`,  O8fitf -> -0.0007675526757783264`}}; not very different O4 and O6; O8 is different (because includes effects of O8 and O10 of the "full" fit with O10) *)
alO4O6O8O10fitPFOMinnum[9]={0.0000749698, {aALf -> 0.102091, O4fitf -> -0.00476467, O6fitf -> 0.0037006, O8fitf -> -0.00230686, O10fitf -> 0.000920885}}; 
alO4O6O8O10fitPFOMinnum[10]={0.0000908152, {aALf -> 0.101678, O4fitf -> -0.0062089, O6fitf -> 0.00420978, O8fitf -> -0.00256626, O10fitf -> 0.00101999}}; 
alO4O6O8O10fitPFOMinnum[11]={0.000117097, {aALf -> 0.101585, O4fitf -> -0.00844193, O6fitf -> 0.0049338, O8fitf -> -0.00293857, O10fitf -> 0.00116411}}; 
alO4O6O8O10fitPFOMinnum[12]={0.000136613, {aALf -> 0.101583, O4fitf -> -0.0117986, O6fitf -> 0.00576187, O8fitf -> -0.00330336, O10fitf -> 0.00129185}}; 
alO4O6O8O10fitPFOMinnum[13]={0.000171149, {aALf -> 0.102189, O4fitf -> -0.0182629, O6fitf -> 0.00711369, O8fitf -> -0.00387178, O10fitf -> 0.00148846}}; 
alO4O6O8O10fitPFOMinnum[14]={0.000112264, {aALf -> 0.101424, O4fitf -> -0.0251768, O6fitf -> 0.00740164, O8fitf -> -0.00365639, O10fitf -> 0.00133653}}; 
alO4O6O8O10fitPFOMinnum[15]={0.000342205, {aALf -> 0.105098, O4fitf -> -0.0652053, O6fitf -> 0.0131878, O8fitf -> -0.00597386, O10fitf -> 0.0021719}};
Clear[Nn]; Do[Module[{}, xi2FO[Nn]=alO4O6O8O10fitPFOMinnum[Nn][[1]]; aALfFO[Nn]=aALf/.alO4O6O8O10fitPFOMinnum[Nn][[2]]; alsmtauFO[Nn]=aMSext[mtauphys^2/(kappa sigmax),aALfFO[Nn]] Pi;  O4FO[Nn]=O4fitf/.alO4O6O8O10fitPFOMinnum[Nn][[2]]; O6FO[Nn]=O6fitf/.alO4O6O8O10fitPFOMinnum[Nn][[2]]; O8FO[Nn]=O8fitf/.alO4O6O8O10fitPFOMinnum[Nn][[2]]]; O10FO[Nn]=O10fitf/.alO4O6O8O10fitPFOMinnum[Nn][[2]],{Nn,3,15}];
 *)
 (* End of FO Borel-Laplace (PFO) *)
(*---------------------------------------------------*)
(* global fits FO, with O10 varied: The behaviour of a^{(2,1)}(aAL,O6,O8) and a^{(2,0)}(aAL,O4,O6) suggests that the truncation order Nt in the FO global approach (i.e., the approach of minimization with alpha_s variable) is Nt=8; this same conclusion, Nt=8, is also valid in the "canonically truncated OPE case," i.e., when O10=0 is kept;  *)
(* For the latter reason, we will now keep Nt=8 in the FO global fit with O10 varied; this will be done in complete analogy with the "canonically truncated OPE case" (O10=0), see prg SUMRPMSbALEPHM2m090.m *)

(* Thus the central case of FO global fit is: *)			     
alO4O6O8O10fitPFOMinnum[8]={0.0000635937, {aALf -> 0.102593, O4fitf -> -0.0036694, O6fitf -> 0.00326539, O8fitf -> -0.00208509, O10fitf -> 0.00083781}};

(* Now the estimate of experimental uncertainty in this global approach with double-pinched FO Borel-Laplace (with Nt=8): *)
Nt=8;
(* Number of equidistant points Nm in the interval 0.9 =< |M^2| =< 1.5 GeV^2: *)
Nm=2; Nd=3 Nm;  
M2cov1F[ii_]:=0.90+(ii-1)*0.60/(Nm-1); (* ii=1,..,Nm; this gives M2covF[1]=0.90, ... M2covF[Nm]=1.5 *)
Do[M2covF[jj]=M2cov1F[jj],{jj,1,Nm}]; Do[M2covF[jj+Nm]=M2cov1F[jj] Exp[I Pi/6],{jj,1,Nm}]; Do[M2covF[jj+2 Nm]=M2cov1F[jj] Exp[I Pi/4],{jj,1,Nm}]; (* in this way we have 3 Nm=Nd points M^2 in the complex M^2-plane, three with Psi=0, three with Psi=Pi/6, and three with Psi=Pi/4 *)

(* Now vary O4 in UNCORRELATED chi ('var'): *)
(* Note: VBBexpFAux[M2,M2] = delReBPexpBinAL[M2]^2, as it should be *)
Nt0=8;
chi2o2O4[O4f_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALfFO[Nt0],O4f,O6FO[Nt0],O8FO[Nt0],O10FO[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O4FOvarfitMin=FindMinimum[chi2o2O4[O4fitf],{O4fitf,O4FO[Nt0]}];
O4FOvarfit=O4fitf/.O4FOvarfitMin[[2]];
Xi2O4FOvaro77=O4FOvarfitMin[[1]];                       
 delO4FO= delO4f/.FindRoot[chi2o2O4[O4FOvarfit+delO4f]==Xi2O4FOvaro77+1,{delO4f, 3. 10^-5}];
(* --> for Nm=2 (Nd=6) we obtain: delO4FO=0.0000304; note that O4FO[8]=-0.0036694, by global analysis *)
(* I will take uncorrelated uncertainty,  Nm=2 (Nd=6) delO4FO=0.0000304, to avoid possible subestimations (the correlated method gives very small uncertainties) *)

chi2o2O6[O6f_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALfFO[Nt0],O4FO[Nt0],O6f,O8FO[Nt0],O10FO[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O6FOvarfitMin=FindMinimum[chi2o2O6[O6fitf],{O6fitf,O6FO[Nt0]}];
O6FOvarfit=O6fitf/.O6FOvarfitMin[[2]];
Xi2O6FOvaro77=O6FOvarfitMin[[1]];
 delO6FO= delO6f/.FindRoot[chi2o2O6[O6FOvarfit+delO6f]==Xi2O6FOvaro77+1,{delO6f, 3. 10^-5}];
(* Nm=2; Nd=6; delO6FO=0.0000451; *)

chi2o2O8[O8f_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALfFO[Nt0],O4FO[Nt0],O6FO[Nt0],O8f,O10FO[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O8FOvarfitMin=FindMinimum[chi2o2O8[O8fitf],{O8fitf,O8FO[Nt0]}];
O8FOvarfit=O8fitf/.O8FOvarfitMin[[2]];
Xi2O8FOvaro77=O8FOvarfitMin[[1]];
 delO8FO= delO8f/.FindRoot[chi2o2O8[O8FOvarfit+delO8f]==Xi2O8FOvaro77+1,{delO8f, 3. 10^-5}];
(* Nm=2; Nd=6; delO8FO=0.0000809; *)

chi2o2O10[O10f_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALfFO[Nt0],O4FO[Nt0],O6FO[Nt0],O8FO[Nt0],O10f] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O10FOvarfitMin=FindMinimum[chi2o2O10[O10fitf],{O10fitf,O10FO[Nt0]}];
O10FOvarfit=O10fitf/.O10FOvarfitMin[[2]];
Xi2O10FOvaro77=O10FOvarfitMin[[1]];
 delO10FO= delO10f/.FindRoot[chi2o2O10[O10FOvarfit+delO10f]==Xi2O10FOvaro77+1,{delO10f, 3. 10^-5}];
(* Nm=2; Nd=6; delO10FO=0.000189576; *)

chi2o2aAL[aALf_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALf,O4FO[Nt0],O6FO[Nt0],O8FO[Nt0],O10FO[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
aALFOvarfitMin=FindMinimum[chi2o2aAL[aALfitf],{aALfitf,aALfFO[Nt0]}];
aALFOvarfit=aALfitf/.aALFOvarfitMin[[2]];
Xi2aALFOvaro77=aALFOvarfitMin[[1]];
 delaALFO= delaALf/.FindRoot[chi2o2aAL[aALFOvarfit+delaALf]==Xi2aALFOvaro77+1,{delaALf, 3. 10^-5}];
(* Nm=2; Nd=6; Nt=8; delaALFO=0.0000961907 *)
(*-------*)
(* Together, the estimates of EXPERIMENTAL uncertainties in the FO global approach (with Nt=6) are (Nm=2; Nd=6), FO global:  *)
(* delO4FO= 0.000030; delO6FO=0.000045; delO8FO=0.000081; delO10FO=0.000190; delaALFO=0.000096(=delta(a(sigmax))) *)
(* --> alpha_s(mtau^2)=0.311494+{^{+0.000279}_{-0.000279}}_{exp}; *)
(* 10^3 O4=-3.7+{^{+0.1}_{-0.0}}_{exp}; 10^3 O6=+3.3+{{+0.0}_{-0.1}}_{exp}; 10^3 O8=-2.1+{^{+0.1}_{-0.1}}_{exp}; 10^3 O10=+0.8+{^{+0.2}_{-0.2}}_{exp};*)

(* The precise central value result of the FO global fit:
alO4O6O8O10fitPFOMinnum[8]={0.0000635937, {aALf -> 0.102593, O4fitf -> -0.0036694, O6fitf -> 0.00326539, O8fitf -> -0.00208509, O10fitf -> 0.00083781}};
 (*-----*)
(* Variation of kappa, FO global: *)
(* kappa=2.0: 
alO4O6O8O10fitPFOMinnum[8]={0.0000396325, {aALf -> 0.0864191, O4fitf -> -0.00227956, O6fitf -> 0.00246831, O8fitf -> -0.00160916, O10fitf -> 0.00064656}};
  -> alpha_s(mtau^2)=0.314407 *)
(* kappa=1.0 (the central case of the global fit): 
alO4O6O8O10fitPFOMinnum[8]={0.0000635937, {aALf -> 0.102593, O4fitf -> -0.0036694, O6fitf -> 0.00326539, O8fitf -> -0.00208509, O10fitf -> 0.00083781}};
  -> alpha_s(mtau^2)=0.311494 *)
(* kappa=0.5:
alO4O6O8O10fitPFOMinnum[8]={0.0000787165, {aALf -> 0.127161, O4fitf -> -0.00599366, O6fitf -> 0.00415249, O8fitf -> -0.00249493, O10fitf -> 0.000978365}};
 -> alpha_s(mtau^2)=0.307221 *)
(* --> thus the variations: *)
(* alpha_s(mtau^2)=0.311494 + {^{+0.002913}_{-0.004273}}_{kappa} *)
(* 10^3 O4=-3.7+{^{+1.4}_{-2.3}}_{kappa}; 10^3 O6=+3.3+{{-0.8}_{+0.9}}_{kappa}; 10^3 O8=-2.1+{^{+0.5}_{-0.4}}_{kappa}; 10^3 O10=+0.8+{^{-0.2}_{+0.2}}_{kappa};*)
(*-----*)
(* Variation of d4, FO global (note, now: Nt=5, not Nt=8): *)
(* d4=2 338.19 (tddtk[4,tkap]=37.77+338.19): 
alO4O6O8O10fitPFOMinN[5]={0.0000123295, {aALf -> 0.102874, O4fitf -> -0.000658707, O6fitf -> 0.0013847, O8fitf -> -0.000911732, O10fitf -> 0.000354432}};
---> alpha_s(mtau^2)=0.312311; *) 
(* d4=338.19;  (tddtk[4,tkap]=37.77) (the central case of the global fit):
alO4O6O8O10fitPFOMinN[5]={0.0000195842, {aALf -> 0.104093, O4fitf -> -0.000848814, O6fitf -> 0.00157982, O8fitf -> -0.00106692, O10fitf -> 0.000427771}};
---> alpha_s(mtau^2)=0.31585; *) 
(* d4=0 (tddtk[4,tkap]=37.77-338.19): 
alO4O6O8O10fitPFOMinN[5]={0.0000300075, {aALf -> 0.105451, O4fitf -> -0.00106834, O6fitf -> 0.00180484, O8fitf -> -0.00124575, O10fitf -> 0.000512247}}; 
---> alpha_s(mtau^2)=0.319786; *)
(* --> thus the variations: *)
 (* alpha_s(mtau^2)=0.31585+{^{-0.003539}_{+0.003936}}_{d4}; *)
(* 10^3 O4=-0.85+{^{+0.2}_{-0.2}}_{d4}; 10^3 O6= +1.6+{{-0.2}_{+0.2}}_{d4}; 10^3 O8=-1.1+{^{+0.2}_{-0.2}}_{d4}; 10^3 O10=+0.4+{^{-0.1}_{+0.1}}_{d4};*)
(*-----*)
(* Variation of Nt(truncation order, Nt=8+-1), FO global: *)
(* Nt=9: alO4O6O8O10fitPFOMinnum[9]={0.0000749698, {aALf -> 0.102091, O4fitf -> -0.00476467, O6fitf -> 0.0037006, O8fitf -> -0.00230686, O10fitf -> 0.000920885}}; --> alpha_s(mtau^2)=0.310035; *)
(* Nt=8: alO4O6O8O10fitPFOMinnum[8]={0.0000635937, {aALf -> 0.102593, O4fitf -> -0.0036694, O6fitf -> 0.00326539, O8fitf -> -0.00208509, O10fitf -> 0.00083781}}; --> alpha_s(mtau^2)=0.311494; *)
(* Nt=7: alO4O6O8O10fitPFOMinnum[7]={0.0000543979, {aALf -> 0.103123, O4fitf -> -0.00274054, O6fitf -> 0.00283537, O8fitf -> -0.00186052, O10fitf -> 0.000754875}}; --> alpha_s(mtau^2)=0.313034; *)
(* --> thus the variations: *)
(* alpha_s(mtau^2)=0.311494 + {^{-0.001459}_{+0.00154}}_{Nt} *)
(* 10^3 O4=-3.7+{^{-1.1}_{+0.9*}}_{Nt}; 10^3 O6=+3.3+{{+0.4}_{-0.4*}}_{Nt}; 10^3 O8=-2.1+{^{-0.2}_{+0.2}}_{Nt}; 10^3 O10=+0.8+{^{+0.1}_{-0.1}}_{Nt};*) 
(*-----------------------------------------*)
(* Altogether, the FO global fit approach gives: *)
(* Nt=8; *)
(* alpha_s(mtau^2)=0.311494 +{^{+0.000279}_{-0.000279}}_{exp} + {^{+0.002913}_{-0.004273}}_{kappa} + {^{-0.003539}_{+0.003936}}_{d4} + {^{-0.001459}_{+0.00154}}_{Nt}  = 0.311494 +{^{+0.005141}_{-0.005744}} =app 0.3115+{^{+0.0051}_{-0.0057}} =app 0.311 + {^{+0.006}_{-0.005}};
10^3 O4= -3.7 + {^{+0.1}_{-0.0}} {exp} + {^{+1.4}_{-2.3}}_{kappa} + {^{+0.2}_{-0.2}}_{d4} + {^{-1.1}_{+0.9*}}_ {Nt} =  -3.7 + {^{+1.7}_{-2.6 }}
10^3 O6 = +3.3 + {{+0.0}_{-0.1}}_{exp} + {{-0.8}_{+0.9}}_{kappa} + {{-0.2}_{+0.2}}_{d4} + {{+0.4}_{-0.4*}}_{Nt} =  +3.3 + {{+1.0}_{-0.9*}}
10^3 O8 = -2.1 + {^{+0.1}_{-0.1}}_{exp} + {^{+0.5}_{-0.4}}_{kappa} + {^{+0.2}_{-0.2}}_{d4} + {^{-0.2}_{+0.2}}_{Nt} =  -2.1 + {^{+0.6}_{-0.5}}
10^3 O10 = +0.8 + {^{+0.2}_{-0.2}}_{exp} + {^{-0.2}_{+0.2}}_{kappa} + {^{-0.1}_{+0.1}}_{d4} + {^{+0.1}_{-0.1}}_{Nt} = +0.8 + {^{+0.3}_{-0.3}}

 (* Central case: *) 
(* Nt=8: alO4O6O8O10fitPFOMinnum[8]={0.0000635937, {aALf -> 0.102593, O4fitf -> -0.0036694, O6fitf -> 0.00326539, O8fitf -> -0.00208509, O10fitf -> 0.00083781}}; --> alpha_s(mtau^2)=0.311494; *) 

 (* End of global FO fit *)
 (*---------------------------------------------------------------------_*)
(*---------------------------------------------------------------------_*)
 (* 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} *)
(*---*)
(* Now the experimental ALEPH values for (1,n) moments a^{(1,n)} (see also the file rtauALEPH.m): *)
 grs0a[1,nn_][sig_]:=((nn+2)/(nn+1)) (1/sigmax) (1 - (sig/sigmax)^(nn+1));
aexpBin77[1,nn_]:=(Sum[ dsbin[[ii]] grs0a[1,nn][sbin[[ii]]] OmegaVplA[ii],{ii,1,77}] + 2 Pi^2 fPi^2 grs0a[1,nn][mPi^2]) -  1;
delaexpBin77[1,nn_]:=(2 Pi^2) Sqrt[Sum[ grs0a[1,nn][sbin[[ii]]] grs0a[1,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}]];
(*  a^{(1,n)}_{exp} = aexpBin77[1,n] +- delaexpBin77[1,n] *)
(*--------*)
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] *)
aaFOD0fast[2,n_,aALf_][NN_]:=Sum[coeffaaFO[2,n][nn] x^nn,{nn,1,NN}]/.{x->aALf};
(* aALf=a(kappa*sigmax) *)
aaFOD0fastIR2[2,n_,aALf_][NN_]:=Sum[coeffaaFOIR2[2,n][nn] x^nn,{nn,1,NN}]/.{x->aALf};
aaFOD0fastIR3[2,n_,aALf_][NN_]:=Sum[coeffaaFOIR3[2,n][nn] x^nn,{nn,1,NN}]/.{x->aALf};
aaFOD0fastUV1[2,n_,aALf_][NN_]:=Sum[coeffaaFOUV1[2,n][nn] x^nn,{nn,1,NN}]/.{x->aALf}; 
 aaFOfast[2,n_,aALf_][NN_,O2np4_,O2np6_]:=aaFOD0fast[2,n,aALf][NN] + ((n+3)/(n+1)) 2 Pi^2 (-1)^n ( (n+2) O2np4/sigmax^(n+2) + (n+1) O2np6/sigmax^(n+3) );
(* O2np4= <O_{2 n + 4}>_{V+A}; O2np6= <O_{2 n + 6}>_{V+A}; *)
(* Note, on 2 March 2021 I compared with Cesar and found out an error in my prg in aaFOfast[2,n], aaCIfast[2,n], aaPVfast[2,n], aaRenfast[2,n]: the factor in front of the condensates must be (n+3)/(n+1), and NOT (n+3)/(n+2); this is numerically important *)
 (* aALf=alpha_s(kappa*sigmax)/Pi *)
(* moments a^{(2,n}}, tFO: *)
coeffaatFO[2,n_][nn_]:=Sum[tddtk[nn-kk,tkap] be0^kk  cAg[2,n][kk] Binomial[nn,kk],{kk,0,nn}];
coeffaatFOIR2[2,n_][nn_]:=Sum[tddtkIR2[nn-kk,tkap] be0^kk  cAg[2,n][kk] Binomial[nn,kk],{kk,0,nn}];
coeffaatFOIR3[2,n_][nn_]:=Sum[tddtkIR3[nn-kk,tkap] be0^kk  cAg[2,n][kk] Binomial[nn,kk],{kk,0,nn}];
coeffaatFOUV1[2,n_][nn_]:=Sum[tddtkUV1[nn-kk,tkap] be0^kk  cAg[2,n][kk] Binomial[nn,kk],{kk,0,nn}]; 
aatFOD0[2,n_,aALf_][NN_]:=Sum[coeffaatFO[2,n][nn-1] taprov[nn,x],{nn,1,NN}]/.{x->aALf};
(* aALf=a(kappa*sigmax) *)
aatFOD0IR2[2,n_,aALf_][NN_]:=Sum[coeffaatFOIR2[2,n][nn-1] taprov[nn,x],{nn,1,NN}]/.{x->aALf};
aatFOD0IR3[2,n_,aALf_][NN_]:=Sum[coeffaatFOIR3[2,n][nn-1] taprov[nn,x],{nn,1,NN}]/.{x->aALf};
aatFOD0UV1[2,n_,aALf_][NN_]:=Sum[coeffaatFOUV1[2,n][nn-1] taprov[nn,x],{nn,1,NN}]/.{x->aALf}; 
aatFO[2,n_,aALf_][NN_,O2np4_,O2np6_]:=aatFOD0[2,n,aALf][NN] + ((n+3)/(n+1)) 2 Pi^2 (-1)^n ( (n+2) O2np4/sigmax^(n+2) + (n+1) O2np6/sigmax^(n+3) );
 (*----*)
(* moments a^{(1,n}}, FO: *)			     
cAg[1,nn_][kp_]:=(1/(2 Pi)) ((nn+2)/(nn+1)) Re[ (-I)^kp NIntegrate[ th^kp ( (1+Exp[I th]) - (1/(nn+2))(1 - (-Exp[I th])^(nn+2)) ),{th, -Pi,Pi}] ];  
(* kp = 0, 1, 2, ... *)
coeffaaFO[1,n_][NN_]:=Sum[tddtk[nn,tkap] Binomial[nn + kk, kk] be0^kk coeffext[nn + kk + 1, NN] cAg[1,n][kk], {nn, 0, NN - 1}, {kk, 0, NN - 1 - nn}];
(* coeffaaFO[1,n][NN] is the coefficient at x^NN in the expansion of aa[1,n] *)
aaFOD0fast[1,n_,aALf_][NN_]:=Sum[coeffaaFO[1,n][nn] x^nn,{nn,1,NN}]/.{x->aALf};
aaFOfast[1,n_,aALf_][NN_, O2np4_]:=aaFOD0fast[1,n,aALf][NN]+((n+2)/(n+1)) 2 Pi^2 (-1)^n O2np4/sigmax^(n+2);
(* moments a^{(1,n}}, tFO: *)				     
coeffaatFO[1,n_][nn_]:=Sum[tddtk[nn-kk,tkap] be0^kk  cAg[1,n][kk] Binomial[nn,kk],{kk,0,nn}];
aatFOD0[1,n_,aALf_][NN_]:=Sum[coeffaatFO[1,n][nn-1] taprov[nn,x],{nn,1,NN}]/.{x->aALf};			     
aatFO[1,n_,aALf_][NN_, O2np4_]:=aatFOD0[1,n,aALf][NN]+((n+2)/(n+1)) 2 Pi^2 (-1)^n O2np4/sigmax^(n+2);
(*--------------*)
(*----------------------------------------*)
 (* When aALf=a(sigmax) is varied ('avar'), i.e.,using alO4O6O8O10fitPFOMinnum[N] (N=3,4,...15) from global FO: *)			     
(*
 Tab20D0FOavar=Table[{NN,aaFOfast[2,0,aALfFO[NN]][NN,0,0]},{NN,3,15}];
 (* Tab20FOavar=Table[{NN,aaFOfast[2,0,aALfFO[NN]][NN,O4FO[NN],O6FO[NN]]},{NN,3,15}];*)   Tab20FOavar=Table[{NN, Tab20D0FOavar[[NN - 2]][[2]] + ((0 + 3)/(0 + 1)) 2 Pi^2 (-1)^0 ((0 + 2) O4FO[NN]/sigmax^(0 + 2) + (0 + 1) O6FO[NN]/sigmax^(0 + 3))}, {NN, 3, 15}];
 delTab20FOavar= Table[{Nn, Tab20FOavar[[Nn - 2]][[2]] - Tab20FOavar[[Nn - 3]][[2]]}, {Nn, 4, 15}];
 Tab21D0FOavar=Table[{NN,aaFOfast[2,1,aALfFO[NN]][NN,0,0]},{NN,3,15}]; 
 (* Tab21FOavar=Table[{NN,aaFOfast[2,1,aALfFO[NN]][NN,O6FO[NN],O8FO[NN]]},{NN,3,15}];*)
Tab21FOavar=Table[{NN, Tab21D0FOavar[[NN - 2]][[2]] + ((1 + 3)/(1 + 1)) 2 Pi^2 (-1)^1 ((1 + 2) O6FO[NN]/sigmax^(1 + 2) + (1 + 1) O8FO[NN]/sigmax^(1 + 3))}, {NN, 3, 15}];
delTab21FOavar= Table[{Nn, Tab21FOavar[[Nn - 2]][[2]] - Tab21FOavar[[Nn - 3]][[2]]}, {Nn, 4, 15}];
Tab10D0FOavar=Table[{NN,aaFOfast[1, 0, aALfFO[NN]][NN, 0]},{NN,3,15}];
Tab10FOavar=Table[{NN,Tab10D0FOavar[[NN - 2]][[2]] + ((0+2)/(0+1)) 2 Pi^2 (-1)^0 O4FO[NN]/sigmax^(0 + 2)}, {NN, 3, 15}];
delTab10FOavar= Table[{Nn, Tab10FOavar[[Nn - 2]][[2]] - Tab10FOavar[[Nn - 3]][[2]]}, {Nn, 4, 15}];  
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,15}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,15}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,15}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,15}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,15}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,15}]; ff10exp=Table[{NN,aexpBin77[1,0]},{NN,3,15}]; ff10plexp=Table[{NN,aexpBin77[1,0]+delaexpBin77[1,0]},{NN,3,15}]; ff10miexp=Table[{NN,aexpBin77[1,0]-delaexpBin77[1,0]},{NN,3,15}];
*)
 (*-------*)
 (*-------*)   
(*
 (* global FO: a21, a20, a10: *)
 a20FOavar=ListPlot[{Tab20D0FOavar, Tab20FOavar, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 15.2}, {0.205, 0.35}}, PlotStyle -> {Black, Red, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N,t], "FO:" Superscript[(a), "(2,0)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
 a21FOavar=ListPlot[{Tab21D0FOavar, Tab21FOavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 15.2}, {0.19, 0.24}}, PlotStyle -> {Black, Red, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N,t], "FO:" Superscript[a, "(2,1)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
a20FOavarZoom =ListPlot[{Tab20D0FOavar, Tab20FOavar, ff20exp, ff20miexp, ff20plexp},PlotRange -> {{2.8, 12.2}, {0.235, 0.250}}, PlotStyle -> {Black, Red, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N,t], "FO:" Superscript[a, "(2,0)"]},Epilog -> {Inset[Style["(b)", 13], Scaled[{0.20, 0.85}]]}]
a21FOavarZoom=ListPlot[{Tab21D0FOavar, Tab21FOavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.20, 0.215}}, PlotStyle -> {Black, Red, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N,t], "FO:" Superscript[a, "(2,1)"]},Epilog -> {Inset[Style["(b)", 13], Scaled[{0.20, 0.85}]]}]
a10FOavar=ListPlot[{Tab10D0FOavar, Tab10FOavar, ff10exp, ff10miexp, ff10plexp}, PlotRange -> {{2.8, 15.2}, {0.18, 0.30}}, PlotStyle -> {Black, Red, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N,t], "FO:" Superscript[(a), "(1,0)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
 a10FOavarZoom= ListPlot[{Tab10D0FOavar, Tab10FOavar, ff10exp, ff10miexp, ff10plexp}, PlotRange -> {{2.8, 12.2}, {0.195, 0.215}}, PlotStyle -> {Black, Red, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}, {Blue, PointSize[0.002]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N,t], "FO:" Superscript[(a), "(1,0)"]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.20, 0.85}]]}]
 *)

(*---------------------------------------------------------------------_*)
(*---------------------------------------------------------------------_*)
(*---------------------------------------------------------------------_*)
(*---------------------------------------------------------------------_*)
(* Double-Pinched ('P') Borel-Laplace in CI approach (PCI): *)

  kappa=1;  (* kappa = mu^2/Q^2, where mu^2 is RScl *)
  tkap=kappa Exp[tKK];
 
 (* Now (double) pinched ('P') Re[BPCI]: *)

(* Integration of 5-loop MSbar RGE in the circle with |Q^2|=kappa sigmax: *)
(* x=Log[kappa sigmax Exp[I th]]; dx= I dth = I dth;  *)
(* a[kappa sigmax Exp[I th]] = y[yh]; *) 
 (* sol1a[aALf_? NumberQ]:=NDSolve[{-I y'[th]==bta[y[th]], y[0]==aALf}, y, {th,-Pi,+Pi}]; *)
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]) *)
 
  (*  aMSp[Q2f_, fif_, alsinf_] := Module[{Q2ff = Q2f, fiff = fif, alsinff = alsinf}, Clear[alsin] ; alsin = alsinff; <</home/gorazd/3danQCD.dir/aMSQ2complS_almtauinput.m; aMSpert[Q2ff, fiff]]; such that aMSp[Q2,fi,als] = a(Q2 Exp[I*fi]) such that alpha_s(mtau^2)=als *)

     dff[th_, Nn_, aALf_] := Module[{aa}, aa = aCI[th, aALf]; Sum[ddtk[nn - 1,tkap] aa^nn, {nn, 1, Nn}]]; (* = d_{D=0}(Q^2=sigmax Exp[I th])^{[Nn]} (TPS [Nn}) such that a[sigmax]=aALf *)

   BPCID0[M2_?NumberQ, Nn_, aALf_]:= (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)]) dff[th,Nn,aALf], {th, -Pi, Pi}];  
   (* aALf= a(kappa sigmax) *)

   ReBPCI[M2f_,Nn_,aALf_,O4_,O6_,O8_,O10_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPCID0[M2f,Nn,aALf]+ 2 Pi^2 ( (1/(M2f sigmax))(2 O4 + O6/sigmax) + (1/M2f^2) (O4 + 2 (O6/sigmax) + (O8/sigmax^2)) + (1/M2f^3) (1/2)(O6 + 2 O8/sigmax + O10/sigmax^2) + (1/M2f^4) (1/6) (O8 +2 O10/sigmax) + (1/M2f^5) (1/24) O10 )]; (* Re[B[M2,Nn,O4,O6,O8,O10]] for CIPT truncated at a^Nn, RScl kap=1 *)
    (* aALf= a(sigmax) *)
  
(* Now search for condensate values with this Borel-Laplace (double) pinched ('P') approach, with FO method for D=0 part: *)
NchP=10; (* 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): *)

O4init = -0.0040; O6init = 0.00100; O8init = -0.0040;  O10init = 0.0030; aALfinit = 0.1110;
alO4O6O8O10fitPCIMinN[NN_]:= FindMinimum[Sum[(ReBPCI[M2FP[jj] Exp[I Pi/6], NN, aALf,O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPCI[M2FP[jj] Exp[I Pi/4], NN, aALf, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPCI[M2FP[jj], NN, aALf, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init}}];
(* aALf=alpha_s(sigmax)/Pi *)
alO4O6O8fitPCIMinN[NN_]:= FindMinimum[Sum[(ReBPCI[M2FP[jj] Exp[I Pi/6], NN, aALf,O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPCI[M2FP[jj] Exp[I Pi/4], NN, aALf, O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPCI[M2FP[jj], NN, aALf, O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}];
(* aALf=alpha_s(sigmax)/Pi *)
(* Do[O4O6O8O10fitPCIMinnum[nn]=O4O6O8O10fitPCIMinN[nn], {nn, 3, 15}]; *)

(* NchP = 10; kappa=1; aALf=alpha_s(sigmax)/Pi is VARIED, and O10 is varied (nonzero), PCI approach: *)
(*
alO4O6O8O10fitPCIMinnum[3]={5.55398*10^-6, {aALf -> 0.117577, O4fitf -> -0.00203946, O6fitf -> 0.000990239, O8fitf -> -0.000853276, O10fitf -> 0.00017157}};
alO4O6O8O10fitPCIMinnum[4]={8.22072*10^-6, {aALf -> 0.113295, O4fitf -> -0.00237773, O6fitf -> 0.00102901, O8fitf -> -0.000969021, O10fitf -> 0.000194119}};
alO4O6O8O10fitPCIMinnum[5]={0.0000117136, {aALf -> 0.111146, O4fitf -> -0.00274081, O6fitf -> 0.000966968, O8fitf -> -0.0010859, O10fitf -> 0.000168532}};
 (* To be compared with (when O10=0 fixed): *)
(* alO4O6O8fitPCIMinnum[5]={0.00023087384541395543`, {aALf -> 0.11068523409892024`, O4fitf -> -0.0025578308620901797`, O6fitf -> 0.0007549465433464493`, O8fitf -> -0.0008047841334926128`}}; *)
alO4O6O8O10fitPCIMinnum[6]={0.0000187257, {aALf -> 0.109448, O4fitf -> -0.00323489, O6fitf -> 0.000722631, O8fitf -> -0.00125299, O10fitf -> 0.0000435186}};
alO4O6O8O10fitPCIMinnum[7]={0.0000302372, {aALf -> 0.108534, O4fitf -> -0.00362907, O6fitf -> 0.000339252, O8fitf -> -0.0014097, O10fitf -> -0.000199733}};
alO4O6O8O10fitPCIMinnum[8]={0.0000756182, {aALf -> 0.107495, O4fitf -> -0.0041504, O6fitf -> -0.000600016, O8fitf -> -0.00167349, O10fitf -> -0.000942878}}; 
 Clear[Nn]; Do[Module[{}, xi2CI[Nn]=alO4O6O8O10fitPCIMinnum[Nn][[1]]; aALfCI[Nn]=aALf/.alO4O6O8O10fitPCIMinnum[Nn][[2]]; alsmtauCI[Nn]=aMSext[mtauphys^2/(kappa sigmax),aALfCI[Nn]] Pi;  O4CI[Nn]=O4fitf/.alO4O6O8O10fitPCIMinnum[Nn][[2]]; O6CI[Nn]=O6fitf/.alO4O6O8O10fitPCIMinnum[Nn][[2]]; O8CI[Nn]=O8fitf/.alO4O6O8O10fitPCIMinnum[Nn][[2]]]; O10CI[Nn]=O10fitf/.alO4O6O8O10fitPCIMinnum[Nn][[2]], {Nn,3,8}]; 
 *)			     

 (* End of CI Borel-Laplace (PCI) *)
(*---------------------------------------------------*)
 (* global fits CI, with O10 varied: The behaviour of a^{(2,1)}(aAL,O6,O8) and a^{(2,0)}(aAL,O4,O6) suggests that the truncation order Nt in the CI global approach (i.e., the approach of minimization with alpha_s variable) is Nt=4-5; in the "canonically truncated OPE case," i.e., when O10=0 is kept, we get Nt=5 *)
(* For the latter reason, we will now keep Nt=5 in the CI global fit with O10 varied; this will be done in complete analogy with the "canonically truncated OPE case" (O10=0), see prg SUMRPMSbALEPHM2m090.m *)

(* Thus the central case of CI global fit is: *)			     
alO4O6O8O10fitPCIMinnum[5]={0.0000117136, {aALf -> 0.111146, O4fitf -> -0.00274081, O6fitf -> 0.000966968, O8fitf -> -0.0010859, O10fitf -> 0.000168532}};

 Nt0=5;										(*-----------*)
(* Experimental uncertainties: *)
(* Number of equidistant points Nm in the interval 0.9 =< |M^2| =< 1.5 GeV^2: *)
Nm=2; Nd=3 Nm;  
M2cov1F[ii_]:=0.90+(ii-1)*0.60/(Nm-1); (* ii=1,..,Nm; this gives M2covF[1]=0.90, ... M2covF[Nm]=1.5 *)
Do[M2covF[jj]=M2cov1F[jj],{jj,1,Nm}]; Do[M2covF[jj+Nm]=M2cov1F[jj] Exp[I Pi/6],{jj,1,Nm}]; Do[M2covF[jj+2 Nm]=M2cov1F[jj] Exp[I Pi/4],{jj,1,Nm}]; (* in this way we have 3 Nm=Nd points M^2 in the complex M^2-plane, two with Psi=0, two with Psi=Pi/6, and two with Psi=Pi/4 *)
(* Now vary O4 in UNCORRELATED chi ('var'): *)
(* Note: VBBexpFAux[M2,M2] = delReBPexpBinAL[M2]^2, as it should be *)
chi2o2O4[O4f_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALfCI[Nt0],O4f,O6CI[Nt0],O8CI[Nt0],O10CI[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O4CIvarfitMin=FindMinimum[chi2o2O4[O4fitf],{O4fitf,O4CI[Nt0]}];
O4CIvarfit=O4fitf/.O4CIvarfitMin[[2]];
Xi2O4CIvaro77=O4CIvarfitMin[[1]];
O4CIvarfitnum=	O4CIvarfit;
delO4CI= delO4f/.FindRoot[chi2o2O4[O4CIvarfitnum+delO4f]==Xi2O4CIvaro77+1,{delO4f, 3. 10^-5}];
(* --> for Nm=2 (Nd=6) we obtain: delO4CI=0.0000304; note that O4CI[5]=-0.00274081, by global analysis *)

chi2o2O6[O6f_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALfCI[Nt0],O4CI[Nt0],O6f,O8CI[Nt0],O10CI[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O6CIvarfitMin=FindMinimum[chi2o2O6[O6fitf],{O6fitf,O6CI[Nt0]}];
O6CIvarfit=O6fitf/.O6CIvarfitMin[[2]];
Xi2O6CIvaro77=O6CIvarfitMin[[1]];
 delO6CI= delO6f/.FindRoot[chi2o2O6[O6CIvarfit+delO6f]==Xi2O6CIvaro77+1,{delO6f, 3. 10^-5}];
(* Nm=2; Nd=6; delO6CI=0.0000451; *)

chi2o2O8[O8f_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALfCI[Nt0],O4CI[Nt0],O6CI[Nt0],O8f,O10CI[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O8CIvarfitMin=FindMinimum[chi2o2O8[O8fitf],{O8fitf,O8CI[Nt0]}];
O8CIvarfit=O8fitf/.O8CIvarfitMin[[2]];
Xi2O8CIvaro77=O8CIvarfitMin[[1]];
 delO8CI= delO8f/.FindRoot[chi2o2O8[O8CIvarfit+delO8f]==Xi2O8CIvaro77+1,{delO8f, 3. 10^-5}];
(* Nm=2; Nd=6; delO8CI=0.0000809; *)

chi2o2O10[O10f_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALfCI[Nt0],O4CI[Nt0],O6CI[Nt0],O8CI[Nt0],O10f] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O10CIvarfitMin=FindMinimum[chi2o2O10[O10fitf],{O10fitf,O10CI[Nt0]}];
O10CIvarfit=O10fitf/.O10CIvarfitMin[[2]];
Xi2O10CIvaro77=O10CIvarfitMin[[1]];
 delO10CI= delO10f/.FindRoot[chi2o2O10[O10CIvarfit+delO10f]==Xi2O10CIvaro77+1,{delO10f, 3. 10^-5}];
(* Nm=2; Nd=6; delO10CI=0.000190; *)

 chi2o2aAL[aALf_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALf,O4CI[Nt0],O6CI[Nt0],O8CI[Nt0],O10CI[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
aALCIvarfitMin=FindMinimum[chi2o2aAL[aALfitf],{aALfitf,aALfCI[Nt0]}];
aALCIvarfit=aALfitf/.aALCIvarfitMin[[2]];
Xi2aALCIvaro77=aALCIvarfitMin[[1]];
(* delaALCdelaALCI= delaALf/.FindRoot[chi2o2aAL[aALCIvarfit+delaALf]==Xi2aALCIvaro77+1,{delaALf, 10.^-4}]; delaALCdelaALCI DOES NOT work in practice, so I checked by hand: chi2o2aAL[aALCIvarfit]=0.0000215648;  chi2o2aAL[aALCIvarfit +-0.000132] =app 1.00 *)
(* Nm=2; Nd=6; Nt0=5; delaALCI=0.000132 *)

(* The precise central value result of the CI global fit: *)
alO4O6O8O10fitPCIMinnum[5]={0.0000117136, {aALf -> 0.111146, O4fitf -> -0.00274081, O6fitf -> 0.000966968, O8fitf -> -0.0010859, O10fitf -> 0.000168532}};
 (*----*)
(* Together, the estimates of EXPERIMENTAL uncertainties in the CI global fit approach (with Nt=6) are (Nm=2; Nd=6; uncorrelated xi^2):  *)
(* delO4CI= 0.000030; delO6CI=0.000045; delO8CI=0.000081; delO10CI=0.000190; delaALCI=0.000132; (=delta(a(sigmax))) *)
(* --> alpha_s(mtau^2)=0.336229+{^{+0.000380}_{-0.000380}}_{exp}; *)  
(* 10^3 O4=-2.7+{^{+0.0}_{-0.1}}_{exp}; 10^3 O6=+1.0+{^{+0.0}_{-0.1}}_{exp}; 10^3 O8=-1.1+{^{+0.1}_{-0.1}}_{exp}; 10^3 O10= +0.2 + {^{+0.2}_{-0.2}}_{exp};*)
(*---------*)
(* Variation of kappa, CI global: *)
(* kappa=2.0: alO4O6O8O10fitPCIMinnum[5]={5.7213*10^-6, {aALf -> 0.0928694, O4fitf -> -0.00209872, O6fitf -> 0.00121855, O8fitf -> -0.000878172, O10fitf -> 0.000282039}}; --> alpha_s(mtau^2)= 0.343184; *)
(* kappa=1.0: alO4O6O8O10fitPCIMinnum[5]={0.0000117136, {aALf -> 0.111146, O4fitf -> -0.00274081, O6fitf -> 0.000966968, O8fitf -> -0.0010859, O10fitf -> 0.000168532}};  --> alpha_s(mtau^2)= 0.336229; *)
  (* kappa=0.5: alO4O6O8O10fitPCIMinnum[5]={0.00020785, {aALf -> 0.140935, O4fitf -> -0.00296599, O6fitf -> -0.000279614, O8fitf -> -0.000718691, O10fitf -> -0.00147545}}; --> alpha_s(mtau^2)=0.329624; *)
  (* ---> alpha_s(mtau^2) =  0.336229 + {^{+0.006955}_{-0.006605}}_{kappa};
   10^3 O4=-2.7 + {^{+0.6}_{-0.3}}_{kappa}; 10^3 O6=+1.0 + {^{+0.2}_{-1.3}}_{kappa}; 10^3 O8=-1.1 + {^{+0.2}_{+0.4}}_{kappa}; 10^3 O10= +0.2 + {^{+0.1}_{-1.7}}_{kappa}; *) 
(*-----*)  
 (* Variation of d4=338.19 +- 338.19, (kappa=1) CI global: *)
(* d4=ddtk[tkap,4]=676.38: alO4O6O8O10fitPCIMinnum[5]={0.0000144131, {aALf -> 0.109238, O4fitf -> -0.0030412, O6fitf -> 0.000954052, O8fitf -> -0.00118078, O10fitf -> 0.000175928}}; --> alpha_s(mtau^2)= 0.330732; *)				(* d4=ddtk[tkap,4]=338.19: alO4O6O8O10fitPCIMinnum[5]={0.0000117136, {aALf -> 0.111146, O4fitf -> -0.00274081, O6fitf -> 0.000966968, O8fitf -> -0.0010859, O10fitf -> 0.000168532}};  --> alpha_s(mtau^2)= 0.336229; *)
(* d4=ddtk[tkap,4]=0:   alO4O6O8O10fitPCIMinnum[4]={8.22072*10^-6, {aALf -> 0.113295, O4fitf -> -0.00237773, O6fitf -> 0.00102901, O8fitf -> -0.000969021, O10fitf -> 0.000194119}};  --> alpha_s(mtau^2)=0.342405; (note: when d4=0, Nt=5 is equal to Nt=4 case) *)
(*--> alpha_s(mtau^2) =  0.336229 + {^{-0.005497}_{+0.006176}}_{d4};
   10^3 O4=-2.7 + {^{-0.3}_{+0.3}}_{d4}; 10^3 O6=+1.0 + {^{+0.0}_{-0.0}}_{d4}; 10^3 O8=-1.1 + {^{-0.1}_{+0.1}}_{d4}; 10^3 O10= +0.2 + {^{+0.0}_{-0.0}}_{d4}; *)
 (*-----*)
 (* Variation of Nt=5 +-1 (kappa=1, CI global): *)
 (* Nt=6: alO4O6O8O10fitPCIMinnum[6]={0.0000187257, {aALf -> 0.109448, O4fitf -> -0.00323489, O6fitf -> 0.000722631, O8fitf -> -0.00125299, O10fitf -> 0.0000435186}}; --> alpha_s(mtau^2)=0.331338; *)
 (* Nt=5: alO4O6O8O10fitPCIMinnum[5]={0.0000117136, {aALf -> 0.111146, O4fitf -> -0.00274081, O6fitf -> 0.000966968, O8fitf -> -0.0010859, O10fitf -> 0.000168532}};   --> alpha_s(mtau^2)= 0.336229; *)
 (* Nt=4: alO4O6O8O10fitPCIMinnum[4]={8.22072*10^-6, {aALf -> 0.113295, O4fitf -> -0.00237773, O6fitf -> 0.00102901, O8fitf -> -0.000969021, O10fitf -> 0.000194119}};  --> alpha_s(mtau^2)=0.342405; *)
 (*--> alpha_s(mtau^2) =  0.336229 + {^{-0.004891}_{+0.006176}}_{Nt};
   10^3 O4=-2.7 + {^{-0.5}_{+0.3}}_{Nt}; 10^3 O6=+1.0 + {^{-0.3}_{+0.0}}_{Nt}; 10^3 O8=-1.1 + {^{-0.2}_{+0.1}}_{Nt}; 10^3 O10= +0.2 + {^{-0.2}_{+0.0}}_{Nt}; *)
 (*------------------------*)
(* Combining together, the CI global fit approach (with O10 varied nonzero) gives: *)
(*
 alpha_s(mtau^2) =  0.336229 + {^{+0.000380}_{-0.000380}}_{exp} + {^{+0.006955}_{-0.006605}}_{kappa} + {^{-0.005497}_{+0.006176}}_{d4} + ^{-0.004891}_{+0.006176}}_{Nt} = 0.336229 + {^{+0.011172}_{-0.009895}} =app 0.3362 + {^{+0.0112}_{-0.0099}} =app 0.336 + {^{+0.011}_{-0.010}}
 10^3 O4=-2.7 + {^{+0.0}_{-0.0*}}_{exp} + {^{+0.64*}_{-0.2*}}_{kappa} + {^{-0.3}_{+0.4*}}_{d4} + {^{-0.49*}_{+0.36*}}_{Nt} = -2.7 + {^{+0.8*}_{-0.6*}}
 10^3 O6=+1.0 + {^{+0.0}_{-0.0}}_{exp} + {^{+0.2}_{-1.25*}}_{kappa} + {^{+0.0}_{-0.0}}_{d4} + {^{-0.25*}_{+0.0}}_{Nt}  =+1.0 + {^{+0.2}_{-1.3}}
 10^3 O8 = -1.1 + {^{+0.1}_{-0.1}}_{exp} + {^{+0.2}_{+0.36*}}_{kappa} + {^{-0.1}_{+0.1}}_{d4} + {^{-0.2}_{+0.1}}_{Nt} =  -1.1 + {^{+0.4}_{-0.2}}
 10^3 O10= +0.2 + {^{+0.2}_{-0.2}}_{exp} + ^{+0.1}_{-1.64*}}_{kappa} +  {^{+0.0}_{-0.0}}_{d4} + {^{-0.1*}_{+0.0}}_{Nt} =  +0.2 + {^{+0.2}_{-1.7}}
 *)

(* Central case: *)
(* Nt=5: alO4O6O8O10fitPCIMinnum[5]={0.0000117136, {aALf -> 0.111146, O4fitf -> -0.00274081, O6fitf -> 0.000966968, O8fitf -> -0.0010859, O10fitf -> 0.000168532}};   --> alpha_s(mtau^2)= 0.336229; *)

(* End of global CI fit *)
  (*---------------------------------------------------------------------_*)

 (* 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} *)
(*---*)
(* Now the experimental ALEPH values for (1,n) moments a^{(1,n)} (see also the file rtauALEPH.m): *)
 grs0a[1,nn_][sig_]:=((nn+2)/(nn+1)) (1/sigmax) (1 - (sig/sigmax)^(nn+1));
aexpBin77[1,nn_]:=(Sum[ dsbin[[ii]] grs0a[1,nn][sbin[[ii]]] OmegaVplA[ii],{ii,1,77}] + 2 Pi^2 fPi^2 grs0a[1,nn][mPi^2]) -  1;
delaexpBin77[1,nn_]:=(2 Pi^2) Sqrt[Sum[ grs0a[1,nn][sbin[[ii]]] grs0a[1,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}]];
(*  a^{(1,n)}_{exp} = aexpBin77[1,n] +- delaexpBin77[1,n] *)
(*--------*)
(* Now the theoretical FESR (2,n) moments a^{(2,n)} in CI: *)

kappa=1;  (* kappa = mu^2/Q^2, where mu^2 is RScl *)
tkap=kappa Exp[tKK];
 
(* Integration of 5-loop MSbar RGE in the circle with |Q^2|=kappa sigmax: *)
(* x=Log[kappa sigmax Exp[I th]]; dx= I dth = I dth;  *)
(* a[kappa sigmax Exp[I th]] = y[yh]; *) 
(*'*)
 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]) *)
   
     dff[th_, Nn_, aALf_] := Module[{aa}, aa = aCI[th, aALf]; Sum[ddtk[nn - 1,tkap] aa^nn, {nn, 1, Nn}]]; (* = d_{D=0}(Q^2=sigmax Exp[I th])^{[Nn]} (TPS [Nn}) such that a[sigmax]=aALf *)
  (*-----*)  
(* moments a^{(2,n}}, CI: *)			     
   aaCID0fast[2,nn_,aALf_][Nn_]:=(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)) ) dff[th,Nn,aALf], {th, -Pi, Pi}];
   (* aALf= a(kappa sigmax) *)
   (*  aaCID0fast[2,nn,aALf][NN] is a^{(2,nn)} moment calculated with the value of  a(kappa sigmax) equal to aALf, calculated to truncation order Nn (~a^Nn) *)
aaCIfast[2,n_,aALf_][NN_,O2np4_,O2np6_]:=aaCID0fast[2,n,aALf][NN] + ((n+3)/(n+1)) 2 Pi^2 (-1)^n ( (n+2) O2np4/sigmax^(n+2) + (n+1) O2np6/sigmax^(n+3) );
(* O2np4= <O_{2 n + 4}>_{V+A}; O2np6= <O_{2 n + 6}>_{V+A}; *)
  (* Note, on 2 March 2021 I compared with Cesar and found out an error in my prg in aaFOfast[2,n], aaCIfast[2,n], aaPVfast[2,n], aaRenfast[2,n]: the factor in front of the condensates must be (n+3)/(n+1), and NOT (n+3)/(n+2); this is numerically important *)
  (*----*)
(* moments a^{(1,n}}, CI: *)
  aaCID0fast[1,nn_,aALf_][Nn_]:=(1/(2 Pi)) ((nn+2)/(nn+1))  NIntegrate[ ( (1 + Exp[I th]) - (1/(nn+2))(1 - (-Exp[I th])^(nn+2)) ) dff[th,Nn,aALf], {th, -Pi, Pi}];
   (* aALf= a(kappa sigmax) *)
aaCIfast[1,n_,aALf_][Nn_, O2np4_]:=aaCID0fast[1,n,aALf][NN]+((n+2)/(n+1)) 2 Pi^2 (-1)^n O2np4/sigmax^(n+2); 
  (*----------*)
(* When aALf=a(sigmax) is varied ('avar'), i.e.,using alO4O6O8O10fitPFOMinnum[N] (N=3,4,...15): *)			     
(*
 Tab20D0CIavar=Table[{NN,Re[aaCID0fast[2,0,aALfCI[NN]][NN]]},{NN,3,8}];
 (* Tab20CIavar=Table[{NN,aaCIfast[2,0,aALfCI[NN]][NN,O4CI[NN],O6CI[NN]]},{NN,3,8}];*)
 Tab20CIavar=Table[{NN, Tab20D0CIavar[[NN - 2]][[2]] + ((0 + 3)/(0 + 1)) 2 Pi^2 (-1)^0 ((0 + 2) O4CI[NN]/sigmax^(0 + 2) + (0 + 1) O6CI[NN]/sigmax^(0 + 3))}, {NN, 3, 8}];
 delTab20CIavar= Table[{NN, Tab20CIavar[[NN - 2]][[2]] - Tab20CIavar[[NN - 3]][[2]]}, {NN, 4, 8}]; (* --> minimum at Nt=5 (i.e., practically no change from Nt=4 to Nt=5) *)
 Tab21D0CIavar=Table[{NN,Re[aaCID0fast[2,1,aALfCI[NN]][NN]]},{NN,3,8}]; 
 (* Tab21CIavar=Table[{NN,aaCIfast[2,1,aALfCI[NN]][NN,O6CI[NN],O8CI[NN]]},{NN,3,8}];*)
  Tab21CIavar=Table[{NN, Tab21D0CIavar[[NN - 2]][[2]] + ((1 + 3)/(1 + 1)) 2 Pi^2 (-1)^1 ((1 + 2) O6CI[NN]/sigmax^(1 + 2) + (1 + 1) O8CI[NN]/sigmax^(1 + 3))}, {NN, 3, 8}];
 delTab21CIavar= Table[{NN, Tab21CIavar[[NN - 2]][[2]] - Tab21CIavar[[NN - 3]][[2]]}, {NN, 4, 8}];  (*-->  local minimum at about Nt=5 (almost no change when from Nt=4 to Nt=5) *)
 Tab10D0CIavar=Table[{NN,Re[aaCID0fast[1, 0, aALfCI[NN]][NN]]},{NN,3,8}];
  Tab10CIavar=Table[{NN,Tab10D0CIavar[[NN - 2]][[2]] + ((0+2)/(0+1)) 2 Pi^2 (-1)^0 O4CI[NN]/sigmax^(0 + 2)}, {NN, 3, 8}];
delTab10CIavar= Table[{NN, Tab10CIavar[[NN - 2]][[2]] - Tab10CIavar[[NN - 3]][[2]]}, {NN, 4, 8}];  (* --> most stable when going from Nt=4 -> Nt=5 --> so Nt=5 is the best *)
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,15}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,15}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,15}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,15}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,15}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,15}];
ff10exp=Table[{NN,aexpBin77[1,0]},{NN,3,15}]; ff10plexp=Table[{NN,aexpBin77[1,0]+delaexpBin77[1,0]},{NN,3,15}]; ff10miexp=Table[{NN,aexpBin77[1,0]-delaexpBin77[1,0]},{NN,3,15}];
  *)
  (*
 a20CIavar=ListPlot[{Tab20D0CIavar, Tab20CIavar, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 8.2}, {0.205, 0.35}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "CI:" Superscript[(a), "(2,0)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
 a20CIavarZoom =ListPlot[{Tab20D0CIavar, Tab20CIavar, ff20exp, ff20miexp, ff20plexp},PlotRange -> {{2.8, 8.2}, {0.235, 0.245}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "CI:" Superscript[a, "(2,0)"]},Epilog -> {Inset[Style["(b)", 13], Scaled[{0.15, 0.95}]]}]
a21CIavar=ListPlot[{Tab21D0CIavar, Tab21CIavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 8.2}, {0.19, 0.22}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "CI:" Superscript[a, "(2,1)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
a21CIavarZoom=ListPlot[{Tab21D0CIavar, Tab21CIavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 8.2}, {0.20, 0.210}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "CI:" Superscript[a, "(2,1)"]},Epilog -> {Inset[Style["(b)", 13], Scaled[{0.40, 0.95}]]}]
 a10CIavar=ListPlot[{Tab10D0CIavar, Tab10CIavar, ff10exp, ff10miexp, ff10plexp}, PlotRange -> {{2.8, 8.2}, {0.19, 0.24}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "CI:" Superscript[(a), "(1,0)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
a10CIavarZoom =ListPlot[{Tab10D0CIavar, Tab10CIavar, ff10exp, ff10miexp, ff10plexp},PlotRange -> {{2.8, 8.2}, {0.20, 0.21}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "CI:" Superscript[a, "(1,0)"]},Epilog -> {Inset[Style["(b)", 13], Scaled[{0.20, 0.85}]]}]
   *) 
(* a20 and a21 (and a10) show stability at Nt=5 (the smallest change when Nt=4->Nt=5) --> so we will choose Nt=5 in the CI global approach *)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(* 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.127;  *)
  aAL=0.10482;
(* 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*Pi/Nn,  DBorelS[N[jj*Pi/Nn]]}, {jj, 1, Nn - 1}]; 

Clear[deltabDBorelS];
tabdelDBorel=Table[{jj*Pi/Nn,  delDBorel[N[jj*Pi/Nn]]}, {jj, 1, Nn - 1}];

Clear[tabaCI];
Do[tabaCI[nn]=Table[{jj*Pi/Nn,  aCI[N[jj*Pi/Nn],aAL]^nn}, {jj, 1, Nn - 1}],{nn,1,15}];

Clear[tabdelD];
  Do[tabdelD[nn]=Table[{jj*Pi/Nn,  Sum[delddtk[mm-1,tkap] tabaCI[mm][[jj]][[2]],{mm,1,nn}]}, {jj, 1, Nn - 1}],{nn,3,15}];

  Do[tabDBorel[nn]=Table[{jj*Pi/Nn, tabDBorelS[[jj]][[2]]+tabdelD[nn][[jj]][[2]]},{jj,1,Nn-1}],{nn,3,15}]; 
  (* Note: tabDBorel[N]= table of the values {jj*Pi/Nn,DBorel[th=jj*Pi/Nn,N]} (jj=1, ...,Nn-1), 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}];

 
ReBPPV[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 ( (1/(M2f sigmax))(2 O4 + O6/sigmax) + (1/M2f^2) (O4 + 2 (O6/sigmax) + (O8/sigmax^2)) + (1/M2f^3) (1/2)(O6 + 2 O8/sigmax + O10/sigmax^2) + (1/M2f^4) (1/6) (O8 +2 O10/sigmax) + (1/M2f^5) (1/24) O10 )]; (* Re[B[M2,Nn,O4,O6,O8,O10]] for PV with delD truncated at a^Nn, RScl kappa *)
    (* aALf= a(kappa sigmax) *)

 (* Now search for condensate values with this Borel-Laplace (double) pinched ('P') approach, with (Borel-resummed) PV method for D=0 part: *)
NchP=10; (* 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): *)

 O4init = -0.0020; O6init = 0.0030; O8init = -0.0020; O10init = 0.0010; (* aALfinit = aAL; *)
O4O6O8O10fitPPVMinN[Nn_]:= FindMinimum[Sum[(ReBPPV[M2FP[jj] Exp[I Pi/6], Nn, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPPV[M2FP[jj] Exp[I Pi/4], Nn, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPPV[M2FP[jj], Nn, O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init}}];
(* aALf=alpha_s(kappa*sigmax)/Pi *)
O4O6O8fitPPVMinN[Nn_]:= FindMinimum[Sum[(ReBPPV[M2FP[jj] Exp[I Pi/6], Nn, O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPPV[M2FP[jj] Exp[I Pi/4], Nn, O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPPV[M2FP[jj], Nn, O4fitf,O6fitf,O8fitf,0] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}];
  
  (* for kappa=1 and NchP=10, PV(+CI) approach: *)

 (* Nt=3 (and O10 varies; kappa=1): *)
 (*
 aAL=0.1450; O4O6O8O10fitPPVMinnum[3]={0.0464116, {O4fitf -> 0.0330244, O6fitf -> 0.0256436, O8fitf -> 0.0217437, O10fitf -> -0.00896405}}
aAL=0.1430; O4O6O8O10fitPPVMinnum[3]={0.0167762, {O4fitf -> 0.0313563, O6fitf -> 0.0230933, O8fitf -> 0.0197076, O10fitf -> -0.0079595}}
aAL=0.1420; O4O6O8O10fitPPVMinnum[3]={0.0109288, {O4fitf -> 0.0305629, O6fitf -> 0.0218301, O8fitf -> 0.0187645, O10fitf -> -0.00748941}}
aAL=0.1411; O4O6O8O10fitPPVMinnum[3]={0.00747757, {O4fitf -> 0.0298438, O6fitf -> 0.0207034, O8fitf -> 0.018109, O10fitf -> -0.00729544}}
aAL=0.1410; O4O6O8O10fitPPVMinnum[3]={0.00735844, {O4fitf -> 0.0297656, O6fitf -> 0.02058, O8fitf -> 0.0180301, O10fitf -> -0.00726297}}
aAL=0.1409; O4O6O8O10fitPPVMinnum[3]={0.00727494, {O4fitf -> 0.0296876, O6fitf -> 0.0204568, O8fitf -> 0.017952, O10fitf -> -0.00723098}}
aAL=0.1408; O4O6O8O10fitPPVMinnum[3]={0.00722655, {O4fitf -> 0.0296098, O6fitf -> 0.0203339, O8fitf -> 0.0178745, O10fitf -> -0.00719945}}
aAL=0.1407; O4O6O8O10fitPPVMinnum[3]={0.00721277, {O4fitf -> 0.0295321, O6fitf -> 0.020211, O8fitf -> 0.0177976, O10fitf -> -0.00716839}}
aAL=0.1406; O4O6O8O10fitPPVMinnum[3]={0.00723309, {O4fitf -> 0.0294547, O6fitf -> 0.0200884, O8fitf -> 0.0177215, O10fitf -> -0.0071378}}
aAL=0.1405; O4O6O8O10fitPPVMinnum[3]={0.007287, {O4fitf -> 0.0293773, O6fitf -> 0.019966, O8fitf -> 0.017646, O10fitf -> -0.00710766}}
aAL=0.1400; O4O6O8O10fitPPVMinnum[3]={0.00804798, {O4fitf -> 0.0289935, O6fitf -> 0.0193549, O8fitf -> 0.0172816, O10fitf -> -0.00696924}}
aAL=0.1390; O4O6O8O10fitPPVMinnum[3]={0.0118026, {O4fitf -> 0.0282376, O6fitf -> 0.0181511, O8fitf -> 0.0165915, O10fitf -> -0.00671082}}
aAL=0.1350; O4O6O8O10fitPPVMinnum[3]={0.0501494, {O4fitf -> 0.0253856, O6fitf -> 0.0135304, O8fitf -> 0.0144313, O10fitf -> -0.00607523}}   
aAL=0.1300; O4O6O8O10fitPPVMinnum[3]={0.127053, {O4fitf -> 0.0222193, O6fitf -> 0.00819888, O8fitf -> 0.0128949, O10fitf -> -0.00595153}}   

  (* ---> when Nt=3 (and O10 varied), xi^2=0.00721277 minimal for aAL=0.1407 -> alpha_s(mtau^2)=0.41966 *)
  (* aALfPV[3](* =aAL *)=0.1407; O4O6O8O10fitPPVMinnum[3]={0.00721277, {O4fitf -> 0.0295321, O6fitf -> 0.020211, O8fitf -> 0.0177976, O10fitf -> -0.00716839}} 
  
(* Nt=4 (and O10 varies; kappa=1): *)
 (* 
aAL=0.11215; O4O6O8O10fitPPVMinnum[4]={0.000101763, {O4fitf -> 0.00367654, O6fitf -> 0.00538088, O8fitf -> -0.000846921, O10fitf -> 0.0002077}};
aAL=0.11220; O4O6O8O10fitPPVMinnum[4]={0.0000941022, {O4fitf -> 0.00367562, O6fitf -> 0.00541847, O8fitf -> -0.000866493, O10fitf -> 0.000215321}};
  aAL=0.11223; O4O6O8O10fitPPVMinnum[4]={0.0000917686, {O4fitf -> 0.00367429, O6fitf -> 0.00544327, O8fitf -> -0.000884241, O10fitf -> 0.000226591}};
aAL=0.11224; O4O6O8O10fitPPVMinnum[4]={0.000091527, {O4fitf -> 0.00367385, O6fitf -> 0.00545155, O8fitf -> -0.000890156, O10fitf -> 0.000230346}};
aAL=0.11225; O4O6O8O10fitPPVMinnum[4]={0.000091554, {O4fitf -> 0.00367341, O6fitf -> 0.00545982, O8fitf -> -0.00089607, O10fitf -> 0.000234101}};  
aAL=0.11226; O4O6O8O10fitPPVMinnum[4]={0.0000918495, {O4fitf -> 0.00367298, O6fitf -> 0.00546809, O8fitf -> -0.000901983, O10fitf -> 0.000237854}};
aAL=0.11229; O4O6O8O10fitPPVMinnum[4]={0.0000943486, {O4fitf -> 0.00367167, O6fitf -> 0.00549293, O8fitf -> -0.00091972, O10fitf -> 0.000249112}}
aAL=0.11230; O4O6O8O10fitPPVMinnum[4]={0.0000957195, {O4fitf -> 0.00367124, O6fitf -> 0.00550121, O8fitf -> -0.000925632, O10fitf -> 0.000252864}}; 
  aAL=0.11235; O4O6O8O10fitPPVMinnum[4]={0.000106612, {O4fitf -> 0.00366908, O6fitf -> 0.00554264, O8fitf -> -0.000955178, O10fitf -> 0.000271611}};
aAL=0.1124; O4O6O8O10fitPPVMinnum[4]={0.000124243, {O4fitf -> 0.00366696, O6fitf -> 0.00558411, O8fitf -> -0.000984711, O10fitf -> 0.000290343}};
aAL=0.1125; O4O6O8O10fitPPVMinnum[4]={0.000179771, {O4fitf -> 0.00366278, O6fitf -> 0.00566717, O8fitf -> -0.00104373, O10fitf -> 0.000327759}};

  (* ---> when Nt=4 (and O10 varied), xi^2=0.000091527 minimal for aAL=0.11224 -> alpha_s(mtau^2)=0.339375 *)
    (*
aALfPV[4](* =aAL *)=0.11224; O4O6O8O10fitPPVMinnum[4]={0.000091527, {O4fitf -> 0.00367385, O6fitf -> 0.00545155, O8fitf -> -0.000890156, O10fitf -> 0.000230346}};
   *)
  
(* when Nt=5: *)
(*
aAL=0.10630; O4O6O8O10fitPPVMinnum[5]={0.00027358, {O4fitf -> 0.000764494, O6fitf -> 0.00371865, O8fitf -> -0.00129767, O10fitf -> 0.000229267}};
aAL=0.10650; O4O6O8O10fitPPVMinnum[5]={0.0000967106, {O4fitf -> 0.00071263, O6fitf -> 0.00389608, O8fitf -> -0.00148235, O10fitf -> 0.000349864}};
aAL=0.10655; O4O6O8O10fitPPVMinnum[5]={0.000083969, {O4fitf -> 0.000699741, O6fitf -> 0.0039404, O8fitf -> -0.00152823, O10fitf -> 0.000379552}};
aAL=0.10656; O4O6O8O10fitPPVMinnum[5]={0.0000829603, {O4fitf -> 0.000697165, O6fitf -> 0.00394927, O8fitf -> -0.0015374, O10fitf -> 0.00038549}};
aAL=0.10657; O4O6O8O10fitPPVMinnum[5]={0.0000824651, {O4fitf -> 0.000694589, O6fitf -> 0.00395814, O8fitf -> -0.00154658, O10fitf -> 0.000391426}};
aAL=0.10658; O4O6O8O10fitPPVMinnum[5]={0.0000824836, {O4fitf -> 0.000692014, O6fitf -> 0.00396701, O8fitf -> -0.00155576, O10fitf -> 0.000397363}};
aAL=0.10660; O4O6O8O10fitPPVMinnum[5]={0.0000840624, {O4fitf -> 0.000686865, O6fitf -> 0.00398476, O8fitf -> -0.00157411, O10fitf -> 0.000409236}};
 aAL=0.10665; O4O6O8O10fitPPVMinnum[5]={0.0000969624, {O4fitf -> 0.000674009, O6fitf -> 0.00402935, O8fitf -> -0.0016212, O10fitf -> 0.000440614}};
aAL=0.10670; O4O6O8O10fitPPVMinnum[5]={0.000122484, {O4fitf -> 0.000661158, O6fitf -> 0.00407377, O8fitf -> -0.00166711, O10fitf -> 0.000470291}};
 *)

 (* ---> when Nt=5 (and O10 varied), xi^2=0.0000824651 minimal for aAL=0.10657 -> alpha_s(mtau^2)=0.323026 *)
 (*
aALfPV[5](* =aAL *)=0.10657; O4O6O8O10fitPPVMinnum[5]={0.0000824651, {O4fitf -> 0.000694589, O6fitf -> 0.00395814, O8fitf -> -0.00154658, O10fitf -> 0.000391426}};
  *)
  
  (* Nt=6 (truncation number): *)
  (*
aALfPV[6](* =aAL *)=0.10405; O4O6O8O10fitPPVMinnum[6]={0.00107164, {O4fitf -> -0.000186278, O6fitf -> 0.00285717, O8fitf -> -0.00119258, O10fitf -> 0.0000140319}}; 
(* when O10=0: O4O6O8fitPPVMinnum[6]={0.00111994, {O4fitf -> -0.000184268, O6fitf -> 0.00285025, O8fitf -> -0.00117891}} *)
aALfPV[6](* =aAL *)=0.10415; O4O6O8O10fitPPVMinnum[6]={0.000768567, {O4fitf -> -0.000219468, O6fitf -> 0.00294075, O8fitf -> -0.00128592, O10fitf -> 0.0000717951}};
aALfPV[6](* =aAL *)=0.10455; O4O6O8O10fitPPVMinnum[6]={0.000202773, {O4fitf -> -0.000341809, O6fitf -> 0.00323602, O8fitf -> -0.00158191, O10fitf -> 0.000232565}};
aALfPV[6](* =aAL *)=0.10460; O4O6O8O10fitPPVMinnum[6]={0.000158912, {O4fitf -> -0.000357498, O6fitf -> 0.00327473, O8fitf -> -0.00162225, O10fitf -> 0.000255592}}
aALfPV[6](* =aAL *)=0.10470; O4O6O8O10fitPPVMinnum[6]={0.0000981246, {O4fitf -> -0.000388838, O6fitf -> 0.00335217, O8fitf -> -0.00170289, O10fitf -> 0.000301393}};
aALfPV[6](* =aAL *)=0.10475; O4O6O8O10fitPPVMinnum[6]={0.0000784076, {O4fitf -> -0.000404741, O6fitf -> 0.00339158, O8fitf -> -0.0017458, O10fitf -> 0.000324634}};
aALfPV[6](* =aAL *)=0.10480; O4O6O8O10fitPPVMinnum[6]={0.0000736056, {O4fitf -> -0.000420171, O6fitf -> 0.00342976, O8fitf -> -0.00178369, O10fitf -> 0.000347371}};
aALfPV[6](* =aAL *)=0.10481; O4O6O8O10fitPPVMinnum[6]={0.0000731467, {O4fitf -> -0.000423303, O6fitf -> 0.00343752, O8fitf -> -0.00179177, O10fitf -> 0.000351966}};
aALfPV[6](* =aAL *)=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}};
aALfPV[6](* =aAL *)=0.10483; O4O6O8O10fitPPVMinnum[6]={0.0000733176, {O4fitf -> -0.000429565, O6fitf -> 0.00345305, O8fitf -> -0.00180794, O10fitf -> 0.000361155}};
aALfPV[6](* =aAL *)=0.10485; O4O6O8O10fitPPVMinnum[6]={0.0000749411, {O4fitf -> -0.000435827, O6fitf -> 0.00346859, O8fitf -> -0.00182411, O10fitf -> 0.000370343}};
aALfPV[6](* =aAL *)=0.10490; O4O6O8O10fitPPVMinnum[6]={0.0000853616, {O4fitf -> -0.000451476, O6fitf -> 0.00350745, O8fitf -> -0.00186454, O10fitf -> 0.000393305}};
aALfPV[6](* =aAL *)=0.10500; O4O6O8O10fitPPVMinnum[6]={0.000133528, {O4fitf -> -0.000482753, O6fitf -> 0.00358524, O8fitf -> -0.00194546, O10fitf -> 0.00043919}};
aALfPV[6](* =aAL *)=0.10510; O4O6O8O10fitPPVMinnum[6]={0.00021824, {O4fitf -> -0.000514003, O6fitf -> 0.00366314, O8fitf -> -0.00202642, O10fitf -> 0.000485026}};
*)

  (* Thus, when Nt=6; the global fit (with O10 included; kappa=1) gives the minimim xi^2= 0.0000730506 for aALfPV[6](* =aAL *)=0.10482, i.e., alpha_s(mtau^2)=0.317958 *)
  (*
aALfPV[6](* =aAL *)=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}};
   *)   

  (* Nt=7 (truncation number): *)
  (*
aALfPV[7](* =aAL *)=0.10400;O4O6O8O10fitPPVMinnum[7]={0.000104978, {O4fitf -> -0.000930034, O6fitf -> 0.00296518, O8fitf -> -0.00184993, O10fitf -> 0.000190924}}
aALfPV[7](* =aAL *)=0.10405;O4O6O8O10fitPPVMinnum[7]={0.0000879307, {O4fitf -> -0.00094857, O6fitf -> 0.00300604,  O8fitf -> -0.00189842, O10fitf -> 0.000219151}};  ***
aALfPV[7](* =aAL *)=0.10410;O4O6O8O10fitPPVMinnum[7]={0.0000856729, {O4fitf -> -0.000967103, O6fitf -> 0.00304692, O8fitf -> -0.00194693, O10fitf -> 0.000247362}};
aALfPV[7](* =aAL *)=0.10412;O4O6O8O10fitPPVMinnum[7]={0.0000889169, {O4fitf -> -0.000974515, O6fitf -> 0.00306328, O8fitf -> -0.00196634, O10fitf -> 0.000258642}};
aALfPV[7](* =aAL *)=0.10415;O4O6O8O10fitPPVMinnum[7]={0.0000982319, {O4fitf -> -0.000985633, O6fitf -> 0.00308782, O8fitf -> -0.00199547, O10fitf -> 0.000275557}};

   aALfPV[7](* =aAL *)=0.10430;O4O6O8O10fitPPVMinnum[7]={0.000223335, {O4fitf -> -0.00104124, O6fitf -> 0.00321091,O8fitf -> -0.00214177, O10fitf -> 0.000361016}};
aALfPV[7](* =aAL *)=0.10432;O4O6O8O10fitPPVMinnum[7]={0.000250267, {O4fitf -> -0.00104865, O6fitf -> 0.0032273, O8fitf -> -0.00216122, O10fitf -> 0.000372269}};
aALfPV[7](* =aAL *)=0.10433;O4O6O8O10fitPPVMinnum[7]={0.000160503, {O4fitf -> -0.00104599, O6fitf -> 0.0032104, O8fitf -> -0.00212213, O10fitf -> 0.000334446}};
aALfPV[7](* =aAL *)=0.10434;O4O6O8O10fitPPVMinnum[7]={0.000168565, {O4fitf -> -0.00104952, O6fitf -> 0.00321793,O8fitf -> -0.00213056, O10fitf -> 0.000338918}};
aALfPV[7](* =aAL *)=0.10435;O4O6O8O10fitPPVMinnum[7]={0.000177011, {O4fitf -> -0.00105306, O6fitf -> 0.00322546, O8fitf -> -0.00213899, O10fitf -> 0.00034339}};
aALfPV[7](* =aAL *)=0.10438;O4O6O8O10fitPPVMinnum[7]={0.000206537, {O4fitf -> -0.00106357, O6fitf -> 0.00324772, O8fitf -> -0.00216356, O10fitf -> 0.00035569}};
aALfPV[7](* =aAL *)=0.1044;O4O6O8O10fitPPVMinnum[7]={0.000226984, {O4fitf -> -0.00107064, O6fitf -> 0.00326279, O8fitf -> -0.00218043, O10fitf -> 0.000364625}};

   (* Thus, when Nt=7; the global fit (with O10 included; kappa=1) gives the minimim xi^2= 0.000160503 for aALfPV[7](* =aAL *)=0.10410, i.e., alpha_s(mtau^2)=0.316537 *)
(*
 aALfPV[7](* =aAL *)=0.10410; ;O4O6O8O10fitPPVMinnum[7]={0.0000856729, {O4fitf -> -0.000967103, O6fitf -> 0.00304692, O8fitf -> -0.00194693, O10fitf -> 0.000247362}};
*)

 (* Nt=8 (truncation number): *)
  (*
aALfPV[8](* =aAL *)=0.10365; O4O6O8O10fitPPVMinnum[8]={0.000125631, {O4fitf -> -0.0013299, O6fitf -> 0.00252617, O8fitf -> -0.00199521, O10fitf -> -0.0000634271}}
aALfPV[8](* =aAL *)=0.10370; O4O6O8O10fitPPVMinnum[8]={0.0000999185, {O4fitf -> -0.00134988, O6fitf -> 0.00256519, O8fitf -> -0.00204514, O10fitf -> -0.000036911}}
aALfPV[8](* =aAL *)=0.10374; O4O6O8O10fitPPVMinnum[8]={0.0000903695, {O4fitf -> -0.00136586, O6fitf -> 0.00259642, O8fitf -> -0.0020851, O10fitf -> -0.0000157208}}
aALfPV[8](* =aAL *)=0.10375; O4O6O8O10fitPPVMinnum[8]={0.000089515, {O4fitf -> -0.00136986, O6fitf -> 0.00260423, O8fitf -> -0.0020951, O10fitf -> -0.0000104264}}
aALfPV[8](* =aAL *)=0.10376; O4O6O8O10fitPPVMinnum[8]={0.0000892742, {O4fitf -> -0.00137385, O6fitf -> 0.00261204, O8fitf -> -0.00210509, O10fitf -> -5.13327*10^-6}}
aALfPV[8](* =aAL *)=0.10377; O4O6O8O10fitPPVMinnum[8]={0.0000896474, {O4fitf -> -0.00137785, O6fitf -> 0.00261985, O8fitf -> -0.00211509, O10fitf -> 1.58633*10^-7}}
aALfPV[8](* =aAL *)=0.10378; O4O6O8O10fitPPVMinnum[8]={0.0000906347, {O4fitf -> -0.00138185, O6fitf -> 0.00262765, O8fitf -> -0.00212508, O10fitf -> 5.44919*10^-6}}   
aALfPV[8](* =aAL *)=0.10380; O4O6O8O10fitPPVMinnum[8]={0.0000944529, {O4fitf -> -0.00138984,O6fitf -> 0.00264327, O8fitf -> -0.00214508, O10fitf -> 0.0000160264}}  
aALfPV[8](* =aAL *)=0.10390; O4O6O8O10fitPPVMinnum[8]={0.000150114, {O4fitf -> -0.00142983, O6fitf -> 0.00272149,O8fitf -> -0.00224532,O10fitf -> 0.0000691018}}
aALfPV[8](* =aAL *)=0.10400; O4O6O8O10fitPPVMinnum[8]={0.000267748, {O4fitf -> -0.00146978, O6fitf -> 0.00279966, O8fitf -> -0.00234548, O10fitf -> 0.00012178}}
 aALfPV[8](* =aAL *)=0.10405; O4O6O8O10fitPPVMinnum[8]={0.000349809, {O4fitf -> -0.00148976, O6fitf -> 0.00283877, O8fitf -> -0.0023956, O10fitf -> 0.00014807}}
aALfPV[8](* =aAL *)=0.10410; O4O6O8O10fitPPVMinnum[8]={0.00044741, {O4fitf -> -0.00150974, O6fitf -> 0.00287788, O8fitf -> -0.00244574, O10fitf -> 0.000174325}}
aALfPV[8](* =aAL *)=0.10415; O4O6O8O10fitPPVMinnum[8]={0.000560585, {O4fitf -> -0.00152971, O6fitf -> 0.00291701,O8fitf -> -0.00249591, O10fitf -> 0.000200545}}
aALfPV[8](* =aAL *)=0.10420; O4O6O8O10fitPPVMinnum[8]={0.00068937, {O4fitf -> -0.00154969, O6fitf -> 0.00295614, O8fitf -> -0.00254611, O10fitf -> 0.000226731}}
aALfPV[8](* =aAL *)=0.10440; O4O6O8O10fitPPVMinnum[8]={0.00102509, {O4fitf -> -0.00162198, O6fitf -> 0.00308297, O8fitf -> -0.00268911, O10fitf -> 0.00027947}}
aALfPV[8](* =aAL *)=0.10450; O4O6O8O10fitPPVMinnum[8]={0.00132018, {O4fitf -> -0.00166023, O6fitf -> 0.00315469, O8fitf -> -0.00277682, O10fitf -> 0.000319906}}
   
   (* Thus, when Nt=8; the global fit (with O10 included; kappa=1) gives the minimim xi^2=0.0000892742 for aALfPV[8](* =aAL *)=0.10376, i.e., alpha_s(mtau^2)=0.314883 *)
(*
 aALfPV[8](* =aAL *)=0.10376; O4O6O8O10fitPPVMinnum[8]={0.0000892742, {O4fitf -> -0.00137385, O6fitf -> 0.00261204, O8fitf -> -0.00210509, O10fitf -> -5.13327*10^-6}};
 *) 

(* In Summary, kappa=1, aALfPV[Nt]=a(kappa sigmax), at Nt truncation in del(Adler), O10 is also varied: *)
 (*   
 aALfPV[3]=0.1407; O4O6O8O10fitPPVMinnum[3]={0.00721277, {O4fitf -> 0.0295321, O6fitf -> 0.020211, O8fitf -> 0.0177976, O10fitf -> -0.00716839}};
aALfPV[4]=0.11224; O4O6O8O10fitPPVMinnum[4]={0.000091527, {O4fitf -> 0.00367385, O6fitf -> 0.00545155, O8fitf -> -0.000890156, O10fitf -> 0.000230346}};
aALfPV[5]=0.10657; O4O6O8O10fitPPVMinnum[5]={0.0000824651, {O4fitf -> 0.000694589, O6fitf -> 0.00395814, O8fitf -> -0.00154658, O10fitf -> 0.000391426}};
aALfPV[6]=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}};
aALfPV[7]=0.10410; ;O4O6O8O10fitPPVMinnum[7]={0.0000856729, {O4fitf -> -0.000967103, O6fitf -> 0.00304692, O8fitf -> -0.00194693, O10fitf -> 0.000247362}};
aALfPV[8]=0.10376; O4O6O8O10fitPPVMinnum[8]={0.0000892742, {O4fitf -> -0.00137385, O6fitf -> 0.00261204, O8fitf -> -0.00210509, O10fitf -> -5.13327*10^-6}};
  kappa=1; Clear[Nn]; Do[Module[{}, xi2PV[Nn]=O4O6O8O10fitPPVMinnum[Nn][[1]]; alsmtauPV[Nn]=aMSext[mtau^2/(kappa sigmax),aALfPV[Nn]] Pi;  O4PV[Nn]=O4fitf/.O4O6O8O10fitPPVMinnum[Nn][[2]]; O6PV[Nn]=O6fitf/.O4O6O8O10fitPPVMinnum[Nn][[2]]; O8PV[Nn]=O8fitf/.O4O6O8O10fitPPVMinnum[Nn][[2]]]; O10PV[Nn]=O10fitf/.O4O6O8O10fitPPVMinnum[Nn][[2]],{Nn,3,8}];
   *)

   (* End of PV Borel-Laplace (PPV) *) 
(*---------------------------------------------------------------------_*)
 (* global fits PV, with O10 varied: The behaviour of a^{(2,1)}(aAL,O6,O8) and a^{(2,0)}(aAL,O4,O6) suggests that the truncation order Nt in the PV global approach (i.e., the approach of minimization with alpha_s variable) is Nt=5; in the "canonically truncated OPE case," i.e., when O10=0 is kept, we get Nt=6 *)
(* For the latter reason, we will now keep Nt=6 in the PV global fit with O10 varied; this will be done in complete analogy with the "canonically truncated OPE case" (O10=0), see prg SUMRPMSbALEPHM2m090.m *)

(* Thus the central case of PV global fit is: *)			     
aALfPV[6]=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}};

 Nt0=6;	(* aALfPV[6]=0.10482;  *)
 (*-----------*)
   (* Experimental uncertainties (estimate): *)
(* Number of equidistant points Nm in the interval 0.9 =< |M^2| =< 1.5 GeV^2: *)
Nm=2; Nd=3 Nm;  
M2cov1F[ii_]:=0.90+(ii-1)*0.60/(Nm-1); (* ii=1,..,Nm; this gives M2covF[1]=0.90, ... M2covF[Nm]=1.5 *)
Do[M2covF[jj]=M2cov1F[jj],{jj,1,Nm}]; Do[M2covF[jj+Nm]=M2cov1F[jj] Exp[I Pi/6],{jj,1,Nm}]; Do[M2covF[jj+2 Nm]=M2cov1F[jj] Exp[I Pi/4],{jj,1,Nm}]; (* in this way we have 3 Nm=Nd points M^2 in the complex M^2-plane, two with Psi=0, two with Psi=Pi/6, and two with Psi=Pi/4 *)
(* Now vary O4 in UNCORRELATED chi ('var'): *)
(* Note: VBBexpFAux[M2,M2] = delReBPexpBinAL[M2]^2, as it should be *)
chi2o2O4[O4f_]:=Sum[(ReBPPV[M2covF[jj],Nt0,O4f,O6PV[Nt0],O8PV[Nt0],O10PV[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O4PVvarfitMin=FindMinimum[chi2o2O4[O4fitf],{O4fitf,O4PV[Nt0]}];
O4PVvarfit=O4fitf/.O4PVvarfitMin[[2]];
Xi2O4PVvaro77=O4PVvarfitMin[[1]];
O4PVvarfitnum=	O4PVvarfit;
delO4PV= delO4f/.FindRoot[chi2o2O4[O4PVvarfitnum+delO4f]==Xi2O4PVvaro77+1,{delO4f, 3. 10^-5}];
(* --> for Nm=2 (Nd=6) we obtain: delO4PV=0.0000304; note that O4PV[6]=-0.000426434, by global analysis *)

chi2o2O6[O6f_]:=Sum[(ReBPPV[M2covF[jj],Nt0,O4PV[Nt0],O6f,O8PV[Nt0],O10PV[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O6PVvarfitMin=FindMinimum[chi2o2O6[O6fitf],{O6fitf,O6PV[Nt0]}];
O6PVvarfit=O6fitf/.O6PVvarfitMin[[2]];
Xi2O6PVvaro77=O6PVvarfitMin[[1]];
 delO6PV= delO6f/.FindRoot[chi2o2O6[O6PVvarfit+delO6f]==Xi2O6PVvaro77+1,{delO6f, 3. 10^-5}];
(* Nm=2; Nd=6; delO6PV=0.0000451; *)

chi2o2O8[O8f_]:=Sum[(ReBPPV[M2covF[jj],Nt0,O4PV[Nt0],O6PV[Nt0],O8f,O10PV[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O8PVvarfitMin=FindMinimum[chi2o2O8[O8fitf],{O8fitf,O8PV[Nt0]}];
O8PVvarfit=O8fitf/.O8PVvarfitMin[[2]];
Xi2O8PVvaro77=O8PVvarfitMin[[1]];
 delO8PV= delO8f/.FindRoot[chi2o2O8[O8PVvarfit+delO8f]==Xi2O8PVvaro77+1,{delO8f, 3. 10^-5}];
(* Nm=2; Nd=6; delO8PV=0.0000809; *)

chi2o2O10[O10f_]:=Sum[(ReBPPV[M2covF[jj],Nt0,O4PV[Nt0],O6PV[Nt0],O8PV[Nt0],O10f] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O10PVvarfitMin=FindMinimum[chi2o2O10[O10fitf],{O10fitf,O10PV[Nt0]}];
O10PVvarfit=O10fitf/.O10PVvarfitMin[[2]];
Xi2O10PVvaro77=O10PVvarfitMin[[1]];
 delO10PV= delO10f/.FindRoot[chi2o2O10[O10PVvarfit+delO10f]==Xi2O10PVvaro77+1,{delO10f, 3. 10^-5}];
(* Nm=2; Nd=6; delO10PV=0.000190; *)

 chi2o2aAL=Sum[(ReBPPV[M2covF[jj],Nt0,O4PV[Nt0],O6PV[Nt0],O8PV[Nt0],O10PV[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
  (* I have to vary aAL (=a(sigmax)) around aALfPV[Nt0]=0.10482, (without varying O4PV[Nt0], O6PV[Nt0], O8PV[Nt0], O10PV[Nt0]), create thus DBorfast[th], evaluate chi2o2aAL, to see when chi2o2aAL increases by unity *)
   (* aAL=0.10482 -> chi2o2aAL=0.0000332221; *)
      (* aAL=0.10482-0.00043 -> chi2o2aAL=1.0222; *)
   (* aAL=0.10482-0.00042 -> chi2o2aAL=0.975291; *)
   (* aAL=0.10482+0.00042 -> chi2o2aAL=0.999063; *)
   
 (* aAL=0.1045 -> chi2o2aAL=; *)
   (* -->  Nm=2; Nd=6; Nt0=6; --> delaALPV=0.00042 (O10=0:  delaALPV=0.0004) *)
   (*---*)
(* The precise central value result of the PV global fit: *)
aALfPV[6]=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}};
 (*----*)
(* Together, the estimates of EXPERIMENTAL uncertainties in the PV global fit approach (with Nt=6) are (Nm=2; Nd=6; uncorrelated xi^2):  *)
(* delO4PV= 0.000030; delO6PV=0.000045; delO8PV=0.000081; delO10PV=0.000190; delaALPV=0.00042; (=delta(a(sigmax))) *)
(* --> alpha_s(mtau^2)=0.3180+{^{+0.0012}_{-0.0013}}_{exp}; *)  
(* 10^3 O4=-0.4+{^{+0.0}_{-0.1}}_{exp}; 10^3 O6=+3.4+{^{+0.1}_{-0.0}}_{exp}; 10^3 O8=-1.8+{^{+0.1}_{-0.1}}_{exp}; 10^3 O10= +0.4 + {^{+0.2}_{-0.2}}_{exp};*)
(*---------*)
 (*---------------------------------------------------*)
(* Variation of kappa, PV global, Nt0=6: *)
   (* kappa=2.0  (note: now aAL=a(2.0 sigmax) : *)
  (*
   aAL=0.08820: xi2[Nt0]=0.000380731;
   aAL=0.08830: xi2[Nt0]=0.000161404;
   aAL=0.08835: xi2[Nt0]=0.0000931314;
   aAL=0.08840: xi2[Nt0]=0.0000525105
   aAL=0.08842: xi2[Nt0]=0.0000443769;
   aAL=0.08845: xi2[Nt0]=0.0000632203;
   aAL=0.08850: xi2[Nt0]=0.0000529446
   aAL=0.08852: xi2[Nt0]=0.0000534826;
   aAL=0.08855: xi2[Nt0]=0.0000599523;
   aAL=0.08857: xi2[Nt0]=0.0000680461;
   aAL=0.08860: xi2[Nt0]=0.0000858676;
   aAL=0.08870: xi2[Nt0]=0.000194664;
   (* xi^2 oscillates around minimum for aAL=0.0884-0.8855; I take aAL(=a(2 sigmax))=0.08850 --> alpha_s(mtau^2)=0.3236 *)
   (* kappa=1.0 (the central case of the global fit): *)
   (* aAL=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}}; *)
   (* kappa=0.5  (note: now aAL=a(0.5 sigmax), Nt0=6 : *)
  (*
   aAL=0.1290; xi2[Nt0]=0.00229163;
   aAL=0.1300; xi2[Nt0]=0.000754336;
   aAL=0.1301; xi2[Nt0]=0.000709915;
   aAL=0.1302; xi2[Nt0]=0.000670159;
   aAL=0.1303; xi2[Nt0]=0.000646416;
   aAL=0.1304; xi2[Nt0]=0.000629703;
   aAL=0.1305; xi2[Nt0]=0.000626583;
   aAL=0.1306; xi2[Nt0]=0.000633829;
   aAL=0.1307; xi2[Nt0]=0.000654224;
   aAL=0.1306; xi2[Nt0]=0.000682013;
   aAL=0.1310; xi2[Nt0]=0.00105243;
   aAL=0.1315; xi2[Nt0]=0.00193193;
   (* xi^2 minimal for aAL=0.1305(=a(0.5 sigmax))--> alpha_s(mtau^2)=0.3128 *)

   (* Thus: *)
 (* kappa=2.0: aAL=0.0885; O4O6O8O10fitPPVMinN[6]={0.0000529446, {O4fitf -> 0.000483166, O6fitf -> 0.00405471, O8fitf -> -0.0017799, O10fitf -> 0.000677834}}; alpha_s(mtau^2)=0.3236; *)
   (* kappa=1: aAL=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}}; alpha_s(mtau^2)=0.3180; *)
   (* kappa=0.5: aAL=0.1305; O4O6O8O10fitPPVMinN[6]={0.000626583, {O4fitf -> -0.000396605, O6fitf -> 0.00148616, O8fitf -> -0.000863364, O10fitf -> -0.0024436}}; alpha_s(mtau^2)=0.3128; *)
   
  (* ---> alpha_s(mtau^2) =  0.3180 + {^{+0.0056}_{-0.0052}}_{kappa};
   10^3 O4=-0.4 + {^{0.9}_{-0.0}}_{kappa}; 10^3 O6=+3.4 + {^{+0.6}_{-1.9}}_{kappa}; 10^3 O8=-1.8 + {^{+0.0}_{+0.9}}_{kappa}; 10^3 O10= +0.4 + {^{+0.3}_{-2.8}}_{kappa}; *) 
(*-----*)
   (* Variation of d4=338.19 +- 338.19, (kappa=1) PV global (note: Nt=5 kept here, NOT Nt=6): *)
   (* delddtk[4,tkap]=1207.16+338.19 =1545.35 (d4=2 338.19= 676.38), Nt=5: *)
(* 
 aAL=0.1050; xi2[5]=0.000252094;
 aAL=0.1051; xi2[5]=0.000156074;
 aAL=0.1052; xi2[5]=0.0000947366;
 aAL=0.1053; xi2[5]=0.0000682017;
 aAL=0.10535; xi2[5]=0.0000680222;
 aAL=0.1054; xi2[5]=0.0000768167;
 aAL=0.1055; xi2[5]=0.00012052;
 (* --> minimum xi^2 at aAL= 0.10535: O4O6O8O10fitPPVMinN[5]={0.0000680222, {O4fitf -> 0.0000886255, O6fitf -> 0.00377405, O8fitf -> -0.00168207, O10fitf -> 0.000444409}}; alpha_s(mtau^2)=0.31949; *)
 
(* delddtk[4,tkap]=1207.16 (d4=338.19); (the central case of the global fit): *)
 (* aAL=aALfPV[5]=0.10657; O4O6O8O10fitPPVMinnum[5]={0.0000824651, {O4fitf -> 0.000694589, O6fitf -> 0.00395814, O8fitf -> -0.00154658, O10fitf -> 0.000391426}};  --> alpha_s(mtau^2)=0.32303; *)

 (* delddtk[4,tkap]=1207.16-338.19=868.969 (d4=0), Nt=5: *)
   (*
    aAL=0.1082: xi2[5]=0.000103124;
    aAL=0.1081: xi2[5]=0.000076859;
    aAL=0.1080: xi2[5]=0.0000818406;
    aAL=0.1079: xi2[5]=0.00011791;
    (* --> minimum xi^2 at aAL=0.1081: O4O6O8O10fitPPVMinN[5]={0.000076859, {O4fitf -> 0.00135556, O6fitf -> 0.00432698, O8fitf -> -0.00146239, O10fitf -> 0.000364965}}; --> alpha_s(mtau^2)=0.327448; *)
    
(*--> alpha_s(mtau^2) =  0.3230 + {^{-0.0035}_{+0.0045}}_{d4};
   10^3 O4=+0.7 + {^{+0.2}_{+0.7}}_{d4}; 10^3 O6=+4.0 + {^{-0.2}_{+0.3}}_{d4}; 10^3 O8=-1.5 + {^{-0.1}_{+0.1}}_{d4}; 10^3 O10= +0.4 + {^{+0.0}_{-0.0}}_{d4}; *)
 (*--------*)
 (* Variation of Nt=6 +-1 (kappa=1, PV global): *)
    (* Nt=7: aAL=aALfPV[7]=0.10410; ;O4O6O8O10fitPPVMinnum[7]={0.0000856729, {O4fitf -> -0.000967103, O6fitf -> 0.00304692, O8fitf -> -0.00194693, O10fitf -> 0.000247362}};  --> alpha_s(mtau^2)=0.3159 *)
 (* Nt=6: aAL=aALfPV[6]=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}}; *)  --> alpha_s(mtau^2)= 0.3180; *)
 (* Nt=5: aALfPV[5]=0.10657; O4O6O8O10fitPPVMinnum[5]={0.0000824651, {O4fitf -> 0.000694589, O6fitf -> 0.00395814, O8fitf -> -0.00154658, O10fitf -> 0.000391426}};  --> alpha_s(mtau^2)=0.32303; *)
 (*--> alpha_s(mtau^2) = 0.3180 + {^{-0.0021*}_{+0.0050}}_{Nt};
   10^3 O4=-0.4 + {^{-0.6}_{+1.1}}_{Nt}; 10^3 O6=+3.4 + {^{-0.4*}_{+0.6}}_{Nt}; 10^3 O8=-1.8 + {^{-0.1*}_{+0.3}}_{Nt}; 10^3 O10= +0.4 + {^{-0.1*}_{0.0}}_{Nt}; *)
(*-----------------------------------------*)
     (* Adding renormalon ambiguity ('amb') to Adler: in ReBPPV we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn] + BPPVdelD0[M2f]) *)
(* Nt=6; O4O6O8fitPPVMinN[6] gives the following xi^2 =xi2[6]: *)
 (*
 aAL=0.10390: xi2[6]=0.000562824;
 aAL=0.10400: xi2[6]=0.000348192;
 aAL=0.10410: xi2[6]=0.000197597;
 aAL=0.10420: xi2[6]=0.000111264;
 aAL=0.10430: xi2[6]=0.0000888236;
 aAL=0.10435: xi2[6]=0.0000816233;
 aAL=0.10440: xi2[6]=0.0000793033;  ***
 aAL=0.10445: xi2[6]=0.0000860671;
 aAL=0.10450: xi2[6]=0.000102407;
 aAL=0.10455: xi2[6]=0.000130435;
  *)
 (* Minimal xi^2 for aAL=0.10440: O4O6O8O10fitPPVMinN[6]={0.0000793033, {O4fitf -> -0.0013398, O6fitf -> 0.00434042, O8fitf -> -0.00184565, O10fitf -> 0.000631921}}; --> alpha_s(mtau^2)=0.3167; *)
 (*---*)
 (* Central value (no ambiguity BPPVdelD0[M2f] added): *)
 (* aAL=aALfPV[6]=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}}; *)  --> alpha_s(mtau^2)= 0.3180; *)
   (*---*)
(* When  in ReBPPV we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn] - BPPVdelD0[M2f]) *)    
(* Nt=6; O4O6O8fitPPVMinN[6] gives the following xi^2 =xi2[6]: *)
 *
   aAL=0.10510: xi2[6]=0.000118398;
   aAL=0.10520: xi2[6]=0.0000760625;
   aAL=0.10524: xi2[6]=0.000068722;
   aAL=0.10525: xi2[6]=0.0000677449;
   aAL=0.10526: xi2[6]=0.0000671113;
   aAL=0.10527: xi2[6]=0.0000668213; ***
   aAL=0.10528: xi2[6]=0.000066875;
   aAL=0.10530: xi2[6]=0.0000680141;   
   aAL=0.10535: xi2[6]=0.0000768851;
   aAL=0.10540: xi2[6]=0.0000942168;
   aAL=0.10550: xi2[6]=0.00015513;
   (* Minimal xi^2 for aAL=0.10527: O4O6O8O10fitPPVMinN[6]={0.0000668213, {O4fitf -> 0.000532149, O6fitf -> 0.00248666, O8fitf -> -0.00176425, O10fitf -> 0.0000514196}}; alpha_s(mtau^2)= 0.3180 + {{-0.0013}_{+0.0013}}_{amb}; *)
   10^3 O4=-0.4 + {^{-0.9}_{+0.9}}_{amb}; 10^3 O6=+3.4 + {^{+0.9}_{-0.9}}_{amb}; 10^3 O8=-1.8 + {^{-0.0}_{+0.0}}_{amb}; 10^3 O10= +0.4 + {^{+0.2}_{-0.3}}_{amb}; *)  (*------------------------*)
(* Combining together, the PV global fit approach (with O10 varied nonzero) gives: *)
(*
 alpha_s(mtau^2)= 0.3180+{^{+0.0012}_{-0.0013}}_{exp}  + {^{+0.0056}_{-0.0052}}_{kappa} + {^{-0.0035}_{+0.0045}}_{d4} + {^{-0.0021*}_{+0.0050}}_{Nt} + {{-0.0013}_{+0.0013}}_{amb} =0.3180+{^{+0.0089}_{-0.0069*}} =app 0.318 + {^{+0.009}_{-0.007}}
 10^3 O4=-0.4 + {^{+0.0}_{-0.1}}_{exp} +  {^{+0.9*}_{-0.0}}_{kappa} + {^{-0.6*}_{+0.7}}_{d4} + {^{-0.6}_{+1.1}}_{Nt} +  {^{-0.9}_{+0.9}}_{amb} = -0.4 + {^{+1.8*}_{-1.2}}
 10^3 O6=+3.4 + {^{+0.1}_{-0.0}}_{exp} + {^{+0.6}_{-1.9}}_{kappa} + {^{-0.2}_{+0.3}}_{d4} + {^{-0.4*}_{+0.6}}_{Nt} + {^{+0.9}_{-0.9}}_{amb} = +3.4 + {^{+1.3}_{-2.1}}
 10^3 O8 =-1.8 + {^{+0.1}_{-0.1}}_{exp} + {^{+0.0}_{+0.9}}_{kappa} + {^{-0.1}_{+0.1}}_{d4} + {^{-0.1*}_{+0.3}}_{Nt} + {^{-0.0}_{+0.0}}_{amb} = -1.8 + {^{+1.0}_{-0.2*}}
 10^3 O10= +0.4 + {^{+0.2}_{-0.2}}_{exp} + {^{+0.3}_{-2.8*}}_{kappa} +  {^{+0.0}_{-0.0}}_{d4} + {^{-0.1*}_{0.0}}_{Nt} + {^{+0.2}_{-0.3}}_{amb} = +0.4 + {^{+0.4}_{-2.8*}}   (asterisk where typos were later corrected)
 *)

(* Central case: *)
 (* Nt=6: aAL=aALfPV[6]=0.10482; O4O6O8O10fitPPVMinnum[6]={0.0000730506, {O4fitf -> -0.000426434, O6fitf -> 0.00344529, O8fitf -> -0.00179985, O10fitf -> 0.000356561}}; *)  --> alpha_s(mtau^2)= 0.3180; *)
   
 (* End of global PV fit *)
  (*---------------------------------------------------------------------_*)
  (*---------------------------------------------------------------------_*)
  (* 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} *)
(*---*)
(* Now the experimental ALEPH values for (1,n) moments a^{(1,n)} (see also the file rtauALEPH.m): *)
 grs0a[1,nn_][sig_]:=((nn+2)/(nn+1)) (1/sigmax) (1 - (sig/sigmax)^(nn+1));
aexpBin77[1,nn_]:=(Sum[ dsbin[[ii]] grs0a[1,nn][sbin[[ii]]] OmegaVplA[ii],{ii,1,77}] + 2 Pi^2 fPi^2 grs0a[1,nn][mPi^2]) -  1;
delaexpBin77[1,nn_]:=(2 Pi^2) Sqrt[Sum[ grs0a[1,nn][sbin[[ii]]] grs0a[1,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}]];
(*  a^{(1,n)}_{exp} = aexpBin77[1,n] +- delaexpBin77[1,n] *)
(*-----------------------*)
(* Now the theoretical FESR (2,n) moments a^{(2,n)} in CI: *)

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

  Nt=6;  (* the truncation order of delD[th,Nn], at ~a^Nn *)

  aAL=aALfPV[Nt];
 (* aAL=a(kappa*sigmax), FIXED here *)  
  
(* 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]) *)

(* aCI[th,aAL] = a(kappa sigmax Exp[I th]), and aAL=a(kappa sigmax) *)
  
delD[th_,Nt]:= Sum[delddtk[jj-1,tkap] aCI[th,aAL]^jj, {jj, 1, Nt}];

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}]);
DBorelS[th_]:=(DBorel1pl[th]+DBorel1mi[th])/2;

DBorel[th_,Nt]:=DBorelS[th]+delD[th,Nt];(* This is PV *)
  (* DBorel[-th,Nt]=Conjugate[DBorel[th,Nt]]; *)
  (* aALf=a(kappa sigmax)= a(tkap Exp[tKK]) *)
 delDBorel[th_]:=(DBorel1pl[th]-DBorel1mi[th])/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=aALtPV[Nt] (where Nt is a fixed truncation order) : *)

 Nn=60; 
  
Clear[tabDBorelS];  
tabDBorelS=Table[{jj*Pi/Nn,  DBorelS[N[jj*Pi/Nn]]}, {jj, 1, Nn - 1}]; 

Clear[tabdelDBorel];
tabdelDBorel=Table[{jj*Pi/Nn,  delDBorel[N[jj*Pi/Nn]]}, {jj, 1, Nn - 1}];

Clear[tabaCI];
Do[tabaCI[nn]=Table[{jj*Pi/Nn,  aCI[N[jj*Pi/Nn],aAL]^nn}, {jj, 1, Nn - 1}],{nn,1,Nt}];

Clear[tabdelD];
  tabdelD[Nt]=Table[{jj*Pi/Nn,  Sum[delddtk[mm-1,tkap] tabaCI[mm][[jj]][[2]],{mm,1,Nt}]}, {jj, 1, Nn - 1}];

tabDBorel[Nt]=Table[{jj*Pi/Nn, tabDBorelS[[jj]][[2]]+tabdelD[Nt][[jj]][[2]]},{jj,1,Nn-1}];
  (* Note: tabDBorel[N]= table of the values {jj*Pi/Nn,DBorel[th=jj*Pi/Nn,N]} (jj=1, ...,Nn-1), 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 *)

  
DBornn[Nt]=Interpolation[tabDBorel[Nt]];
DBorfast[th_,Nt] :=If[th > 0, DBornn[Nt][th], Conjugate[DBornn[Nt][-th]]];
  (* DBorfast[th,Nt]: This is PV Borel resummation of Adler D(Q^2=sigmax Exp[I th]) (D=0 Adler without LO term 1), FAST evaluation, with the polynomial (correction) part delD[th,Nt] truncated at a[kappa sigmax Exp[I th]]^Nt; 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 *)
 (*---------------------------------------------------*)

   (* 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) *)
aaPVfast[2,n_][Nt,O2np4_,O2np6_]:=aaPVD0fast[2,n][Nt] + ((n+3)/(n+1)) 2 Pi^2 (-1)^n ( (n+2) O2np4/sigmax^(n+2) + (n+1) O2np6/sigmax^(n+3) );
(* O2np4= <O_{2 n + 4}>_{V+A}; O2np6= <O_{2 n + 6}>_{V+A}; *)
   (* aaPVD0Sfast[2,nn_]:=(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)) ) DBorSfast[th], {th, -Pi, Pi}]; *)
   (* Note, on 2 March 2021 I compared with Cesar and found out an error in my prg in aaFOfast[2,n], aaCIfast[2,n], aaPVfast[2,n], aaRenfast[2,n]: the factor in front of the condensates must be (n+3)/(n+1), and NOT (n+3)/(n+2); this is numerically important *)
   (*----*)
(* moments a^{(1,n}}, PV (at Nt truncation order): *)
aaPVD0fast[1,nn_][Nt]:=(1/(2 Pi)) ((nn+2)/(nn+1))  NIntegrate[ ( (1 + Exp[I th]) - (1/(nn+2))(1 - (-Exp[I th])^(nn+2)) ) DBorfast[th,Nt], {th, -Pi, Pi}];
   (* aAL (=a(kappa sigmax))=aALfPV[Nt] *)
aaPVfast[1,n_][Nt, O2np4_]:=aaPVD0fast[1,n][Nt]+((n+2)/(n+1)) 2 Pi^2 (-1)^n O2np4/sigmax^(n+2); 
   (* aaPVD0Sfast[1,nn_]:=(1/(2 Pi)) ((nn+2)/(nn+1))  NIntegrate[ ( (1 + Exp[I th]) - (1/(nn+2))(1 - (-Exp[I th])^(nn+2)) ) DBorSfast[th], {th, -Pi, Pi}]; *)
   (* aAL (=a(kappa sigmax))=aALfPV[Nt] *)
(*------*)
   (* Now I have to repeat, for each truncation order Nt=3,4,5,...8, the calculation of  aaPVD0fast[2,0][Nt],  aaPVD0fast[2,1][Nt] and  aaPVD0fast[1,0][Nt] (and their extensions aaPVfast[2,0][Nt,O4PV[Nt],O6PV[Nt]], aaPVfast[2,1][Nt,O6PV[Nt],O8PV[Nt]] and aaPVfast[1,0][Nt,O4PV[Nt]]): *)
(* When aALf=a(sigmax) is varied ('avar'), i.e.,using alO4O6O8fitPPVMinnum[N] (N=3,4,...15): *)
 (*
Tab20D0PVavar= (*=Table[{Nt,Re[aaPVD0fast[2,0][Nt]]},{Nt,3,8}]*){{3,-0.261275},{4,0.170029},{5,0.21906},{6,0.237395},{7, 0.246632},{8,0.253943}};
Tab21D0PVavar= (*=Table[{Nt,Re[aaPVD0fast[2,1][Nt]]},{Nt,3,8}]*){{3,0.335151},{4,0.233925},{5,0.225017},{6,0.221968},{7, 0.219612},{8,0.217049}};
Tab10D0PVavar= (*=Table[{Nt,Re[aaPVD0fast[1,0][Nt]]},{Nt,3,8}]*){{3,0.0479665},{4,0.184494},{5,0.199494},{6,0.205181},{7, 0.207825},{8,0.209846}};
Tab20PVavar=Table[{NN, Tab20D0PVavar[[NN - 2]][[2]] + ((0 + 3)/(0 + 1)) 2 Pi^2 (-1)^0 ((0 + 2) O4PV[NN]/sigmax^(0 + 2) + (0 + 1) O6PV[NN]/sigmax^(0 + 3))}, {NN, 3, 8}];
  delTab20PVavar=Table[{NN, Tab20PVavar[[NN-2]][[2]]-Tab20PVavar[[NN-3]][[2]]},{NN,4,8}];
Tab21PVavar=Table[{NN, Tab21D0PVavar[[NN - 2]][[2]] + ((1 + 3)/(1 + 1)) 2 Pi^2 (-1)^1 ((1 + 2) O6PV[NN]/sigmax^(1 + 2) + (1 + 1) O8PV[NN]/sigmax^(1 + 3))}, {NN, 3, 8}]; delTab21PVavar= Table[{NN, Tab21PVavar[[NN - 2]][[2]] - Tab21PVavar[[NN - 3]][[2]]}, {NN, 4, 8}];
Tab10PVavar=Table[{NN,Tab10D0PVavar[[NN - 2]][[2]] + ((0+2)/(0+1)) 2 Pi^2 (-1)^0 O4PV[NN]/sigmax^(0 + 2)}, {NN, 3, 8}];
  delTab10PVavar=Table[{NN, Tab10PVavar[[NN-2]][[2]]-Tab10PVavar[[NN-3]][[2]]},{NN,4,8}];
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,15}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,15}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,15}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,15}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,15}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,15}];
ff10exp=Table[{NN,aexpBin77[1,0]},{NN,3,15}]; ff10plexp=Table[{NN,aexpBin77[1,0]+delaexpBin77[1,0]},{NN,3,15}]; ff10miexp=Table[{NN,aexpBin77[1,0]-delaexpBin77[1,0]},{NN,3,15}];
 *)   
  (*
 a20PVavar=ListPlot[{Tab20D0PVavar, Tab20PVavar, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 8.2}, {0.16, 0.30}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "PV:" Superscript[(a), "(2,0)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
 a20PVavarZoom =ListPlot[{Tab20D0PVavar, Tab20PVavar, ff20exp, ff20miexp, ff20plexp},PlotRange -> {{2.8, 8.2}, {0.235, 0.245}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "PV:" Superscript[a, "(2,0)"]},Epilog -> {Inset[Style["(b)", 13], Scaled[{0.15, 0.95}]]}]
a21PVavar=ListPlot[{Tab21D0PVavar, Tab21PVavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 8.2}, {0.18, 0.235}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "PV:" Superscript[a, "(2,1)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
a21PVavarZoom=ListPlot[{Tab21D0PVavar, Tab21PVavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 8.2}, {0.200, 0.210}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "PV:" Superscript[a, "(2,1)"]},Epilog -> {Inset[Style["(b)", 13], Scaled[{0.40, 0.95}]]}]
a10PVavar=ListPlot[{Tab10D0PVavar, Tab10PVavar, ff10exp, ff10miexp, ff10plexp}, PlotRange -> {{2.8, 8.2}, {0.18, 0.23}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "PV:" Superscript[(a), "(1,0)"]},Epilog -> {Inset[Style["(a)", 13], Scaled[{0.20, 0.85}]]}]
a10PVavarZoom =ListPlot[{Tab10D0PVavar, Tab10PVavar, ff10exp, ff10miexp, ff10plexp},PlotRange -> {{2.8, 8.2}, {0.200, 0.210}}, PlotStyle -> {Blue, Red, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}, {Black, PointSize[0.005]}}, Filling -> {4 -> {5}}, AxesLabel -> {N, "PV:" Superscript[a, "(1,0)"]},Epilog -> {Inset[Style["(b)", 13], Scaled[{0.20, 0.85}]]}]
   *) 
   (* at Nt=6 we had local minimum of a20, a21 and a10 when O10=0 (now the local minimum is at Nt=5), so I will choose Nt=6 in the PV(+CI) global approach *)
   
    (*-------------------------------------------------*)  
