(* SumMSbALEPHM2O61O62.m: *)
(* as SumRMSbALEPHM2O61O62d4338.m, but for: d4MSbar=275.+ -63.19, and not 338.19+-63.19 *)
(* Note, the string 'OK' indicates that the part there has been converted to the deMSbar=275 case (from d4MSbar=338.19 case) *)
(*-------------------------------------------------*)
(* SumRMSbALEPHM2O61O62d4338.m: as SumRPMSbALEPHM2m090.m, but instead of O6 and O8, we have O6^{(1)} and O6^{(2)} condensates *)
(* Look for: NEW *)
(*-------------------------------------------------*)
(* SumRPMSbALEPHM2m090.m:  with (double) pinched ('P') Borel-Laplace sum rules *)
(* and alpha_s(mtau^2; MSbar) is adjustable (as input, or as result fo global fit) *)
(*-------------------------------------------------*)
(* Borel-Laplace fitted in |M^2|-intervals: (M2min,1.5); M2min=0.65 or 0.90; *)
(* here: M2min=0.90, i.e., YYY=090 *)
(* in pQCD, MSbar scheme *)

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

 (* PDG2021, average from hadronic tau decays: alpha_s(MZ^2)=0.1187 +- 0.0018; [PDG2021 says that this corresponds to alpha_s(mtau^2)=0.325±0.016 (not quite true!); my RGE running is different and basically agrees with Pich RGE running] *)
 (* PDG2021, world average: alpha_s(MZ^2)=0.1176±0.0011 (without lattice); 0.1182±0.0008 (lattice); 0.1179±0.0010 (unweighted average of "without lattice" and "lattice" results) *)
 (*-------------------------------------------------*)
 <</home/gorazd/3danQCD.dir/AdlerFunction4lMiniMOMexact.m;
 (*---------*)
 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275 (ECH-predicted): *)
<</home/gorazd/3danQCD.dir/MSbarRenMod5Ad4275.save; 

 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275.+63.19=338.19: *)
 (* <</home/gorazd/3danQCD.dir/MSbarRenMod5Ad4338.save; when d4barPREDICTED=338.19 (RenMod-predicted in MiniMOM): *)
 
 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275.-63.19=211.81 (=approx 212): *)
 (* <</home/gorazd/3danQCD.dir/MSbarRenMod5Ad4212.save; *)

  (*---------*)
 
BtDUV1[u_, tkapp_] := Exp[Log[tkapp]*u]*Pi*(tdUVp12/(1 + u)^2);
BtDIR3[u_, tkapp_] := 
  Exp[Log[tkapp]*u]*Pi*(tdIRp32/(3 - u)^2 + tdIRp31/(3 - u));
BtDIR2[u_, tkapp_] := 
  Exp[Log[tkapp]*u]*Pi*(tdIRp2*(1/(2 - u) + talf*-1*Log[1 - u/2]));
(* BtDUV1[u,tkap]+BtDIR3[u,tkap]+BtDIR2[u,tkap]=BtDgen[u,tkap] *)
tddtkUV1[nn_, tkapp_] := Coefficient[Series[BtDUV1[u, tkapp], {u, 0, 30}], u^nn] nn! be0^nn;
tddtkUV1[0, tkapp_] := BtDUV1[0, tkapp];
tddtkIR3[nn_, tkapp_] := Coefficient[Series[BtDIR3[u, tkapp], {u, 0, 30}], u^nn] nn! be0^nn;
tddtkIR3[0, tkapp_] := BtDIR3[0, tkapp];
tddtkIR2[nn_, tkapp_] := Coefficient[Series[BtDIR2[u, tkapp], {u, 0, 30}], u^nn] nn! be0^nn;
tddtkIR2[0, tkapp_] := BtDIR2[0, tkapp];
(* tddtkUV1[n, tkap] + tddtkIR3[n, tkap] + tddtkIR2[n, tkap] = tddtk[n,tkap]  +{tilde d}_n(kappa=tkap Exp[-tKK]) *)
(* Check for the onset of the UV1-renormalon dominance in {tilde d}_n, for kappa=1: kappa=1; tkap = kappa Exp[tKK]; tabtddUV1vstdd=Table[{nn, tddtkUV1[nn, tkap]/tddtk[nn, tkap]}, {nn, 1, 15}]={{1, 0.0751416}, {2, -0.227154}, {3, 0.257563}, {4, -1.97777}, {5, 0.573134}, {6, 1.55361}, {7, 0.856263}, {8, 1.08529}, {9, 0.964573}, {10, 1.01747}, {11, 0.992026}, {12, 1.00378}, {13, 0.99823}, {14, 1.00084}, {15, 0.999605}}; in a way UV1 dominates starting at nn=5-6 upwards *)
(* Do[ddtkUV1[n - 1, tkap] = Expand[Sum[kk[n - s, s] tddtkUV1[n - 1 - s, tkap], {s, 0, n - 2}]], {n, 2, 31}];  Do[ddtkIR2[n - 1, tkap] = Expand[Sum[kk[n - s, s] tddtkIR2[n - 1 - s, tkap], {s, 0, n - 2}]], {n, 2, 31}]; Do[ddtkIR3[n - 1, tkap] = Expand[Sum[kk[n - s, s] tddtkIR3[n - 1 - s, tkap], {s, 0, n - 2}]], {n, 2, 31}]; tabddUV1vsdd = Table[{nn, ddtkUV1[nn, tkap]/ddtk[nn, tkap]}, {nn, 1, 15}]={{1, 0.0751416}, {2, -0.0888303}, {3, 0.0786098}, {4, -0.100583}, {5, 0.0988345}, {6, -0.170178}, {7, 0.170568}, {8, -0.445023}, {9, 0.325567}, {10, -3.6304}, {11, 0.56694}, {12, 1.81027}, {13, 0.795075}, {14, 1.17141}, {15, 0.924117}} *)
(* All this indicates that {tilde d}_n of Adler d(Q^2)_{D=0} are clearly dominated by UV1 part starting already at n=5-7; while d_n are dominated by UV1 part starting at HIGHER n (n =11-13); this indicates that FOPT evaluations of FESRs and Borel-Laplace are probably even better behaved when using {tilde a}_n(kappa sigmax) [insetad of powers of a(kappa sigmax)], as the UV1 renormalon of {tilde d}_{D=0} is cancelled in the Borel of {tilde d}_{D=0}, i.e., using {tilde FOPT}  *)


 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]];
(* aMSp[0.8, 0.4, 0.320] is then a(Q^2=0.8 Exp[I 0.4]; MSbar) with initial condition alpha_s(mtau^2;MSbar)=0.320  (all with nf=3) *)

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


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

(* <</home/gorazd/3danQCD.dir/aMSpertComplexnf3al01179.save; *)
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["O4O6O8fitMinALEPH090.save", {M2min,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,-0.9,+1.5}];
(* ' *)
  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 values *)

  (*-----------*) 
  JJ[1, AA_] :=E^-AA - AA (-1 + EulerGamma + Gamma[0, AA] + Log[AA]);
JJ[2,AA_]:=1/2 (-1 + AA) E^-AA - 1/4 AA^2 (-3 + 2 EulerGamma - 2 CoshIntegral[AA] + 2 Log[AA] + 2 SinhIntegral[AA]);
 JJ[3,AA_]:=1/6 (2 + (-1 + AA) AA) E^-AA - 1/36 AA^3 (-11 + 6 EulerGamma - 6 CoshIntegral[AA] + 6 Log[AA] + 6 SinhIntegral[AA]);
  (*
   Alternative equivalent form ('a'):
 JJa[1, AA_] :=E^-AA - AA (-1 + EulerGamma + Gamma[0, AA] + Log[AA]);
   JJa[2,AA_]:=1/2 (-1 + AA) E^-AA - 1/4 AA^2 (-3 + 2 EulerGamma +2 Gamma[0,AA] + 2 Log[AA] );  
   JJa[3,AA_]:=1/6 (2 + (-1 + AA) AA) E^-AA - (1/36) AA^3 (-11 + 6 EulerGamma + 6 Gamma[0,AA] + 6 Log[AA] );
   *)
  ReBPFOfast[M2f_,Nn_,aALf_,O4_,O61_,O62_]:=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 + (O61 + O62 (1/aALf - be0 Log[kappa]))/sigmax) + (1/M2f^2) (O4 + 2 ((O61 + O62 (1/aALf - be0 Log[kappa]))/sigmax)) + (1/M2f^3) (1/2)(O61 + O62 (1/aALf - be0 Log[kappa]))  ) + (6 Pi^2/sigmax^3) be0 O62 (JJ[1,sigmax/M2f] + JJ[2,sigmax/M2f] 2 (1 - M2f/sigmax) + JJ[3,sigmax/M2f] (1 - 2 M2f/sigmax + 2 (M2f/sigmax)^2) - (1/3) 2 (M2f/sigmax)^2 Exp[-sigmax/M2f])];  (* 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; a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE with the initial value a(sigmax)=aALf (the same initial value as in D=0 part, but in the D=0 part a(sigmax Exp[I fi]) runs by FOPT 5-loop)  *)
  (* Note that, instaad of simply O62/aALf, we have O62 (1/aALf - be0 Log[kappa]); this is because aALf=a(kappa sigmax), but at O62 we have to have 1/a(sigmax); we use the one-loop RGE running from Q^2=kappa sigmax) to Q^2=sigmax --> 1/a(sigmax)=(1/aALf - be0 Log[kappa]) *)

  (* Now the dfi fi ... integration numerically: *) 
  ReBPFOfastNum[M2f_,Nn_,aALf_,O4_,O61_,O62_]:=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 + (O61 + O62 (1/aALf - be0 Log[kappa]))/sigmax) + (1/M2f^2) (O4 + 2 ((O61 + O62 (1/aALf - be0 Log[kappa]))/sigmax)) + (1/M2f^3) (1/2)(O61 + O62 (1/aALf - be0 Log[kappa]))  ) + (6 Pi^2/sigmax^3) be0 O62 (I/2/Pi) NIntegrate[fi Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
  (* Note that, instead of simply O62/aALf, we have O62 (1/aALf - be0 Log[kappa]); this is because aALf=a(kappa sigmax), but at O62 we have to have 1/a(sigmax); we use the one-loop RGE running from Q^2=kappa sigmax) to Q^2=sigmax --> 1/a(sigmax)=(1/aALf - be0 Log[kappa]) *)

  (* Now the dfi integration in D=6 term with 1/a(sigmax Exp[I fi]) 5-loop FO expanded (in powers of a(sigmax)): 1/a(sigm Exp[I fi]) = (1/a(sigm)) (1 + ... + O(c4 be0 I fi a(sigm)^5)):  *)
  ReBPFOfastNum5l[M2f_,Nn_,aALf_,O4_,O61_,O62_]:=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 + (O61 + O62/aALf)/sigmax) + (1/M2f^2) (O4 + 2 ((O61 + O62/aALf)/sigmax)) + (1/M2f^3) (1/2)(O61 + O62/aALf)  ) + (6 Pi^2/sigmax^3) O62 (1/2/Pi) NIntegrate[ ( I be0 fi + aALf ce1 I be0 fi  + aALf^2 (ce2MS I be0 fi + (1/2) ce1 be0^2 fi^2) + aALf^3 (ce3MS I be0 fi + (1/2) (2 ce2MS + ce1^2) be0^2 fi^2 - (1/3) ce1 I be0^3 fi^3) + aALf^4 (ce4MS I be0 fi + (3/2) (ce3MS + ce1 ce2MS) be0^2 fi^2 - (ce2MS + (5/6) ce1^2) I be0^3 fi^3 -(1/4) ce1 be0^4 fi^4) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
  (* aFO5linv[fi_,aALf_]:=(1/aALf) + ( I be0 fi + aALf ce1 I be0 fi  + aALf^2 (ce2MS I be0 fi + (1/2) ce1 be0^2 fi^2) + aALf^3 (ce3MS I be0 fi + (1/2) (2 ce2MS + ce1^2) be0^2 fi^2 - (1/3) ce1 I be0^3 fi^3) + aALf^4 (ce4MS I be0 fi + (3/2) (ce3MS + ce1 ce2MS) be0^2 fi^2 - (ce2MS + (5/6) ce1^2) I be0^3 fi^3 -(1/4) ce1 be0^4 fi^4) );   (aFO5linv[fi,aALf] = 1/(a(sigmax Exp[I fi])) in 5-loop FO expansion, where a(sigmax)=aALf)  *)

  (* Now the dfi integration in D=6 term with 1/a(sigmax Exp[I fi]) 5-loop CI(RGE)-evolved:  *)
sol1a[aALf_]:= NDSolve[{-I y'[th] == bta[y[th]], y[0] == aALf}, y, {th, -Pi, +Pi}];
 aCI[th_, aALf_] := First[y[th] /. sol1a[aALf]];
(* aALf=a(kappa sigmax); aCI[th,aALf]=a(Q2=kappa sigma Exp[I th]) ' *)
  
  (* Now the dfi integration in D=6 term with 1/a(sigmax Exp[I fi]) 5-loop CI(RGE)-evolved: *)
  ReBPFOfastNum5lCI[M2f_,Nn_,aALf_,O4_,O61_,O62_]:=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 + O61/sigmax) + (1/M2f^2) (O4 + 2 O61/sigmax) + (1/M2f^3) (1/2) O61  ) + (6 Pi^2/sigmax^3) O62 (1/2/Pi) NIntegrate[ (1/aCI[fi,aALf]) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]]; 

  (*-------------------*)
  
ReBPtFO[M2f_,Nn_,aALf_,O4_,O61_,O62_]:=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 + (O61 + O62 (1/aALf - be0 Log[kappa]))/sigmax) + (1/M2f^2) (O4 + 2 ((O61 + O62 (1/aALf - be0 Log[kappa]))/sigmax)) + (1/M2f^3) (1/2)(O61 + O62 (1/aALf - be0 Log[kappa]))  ) + (6 Pi^2/sigmax^3) be0 O62 (JJ[1,sigmax/M2f] + JJ[2,sigmax/M2f] 2 (1 - M2f/sigmax) + JJ[3,sigmax/M2f] (1 - 2 M2f/sigmax + 2 (M2f/sigmax)^2) - (1/3) 2 (M2f/sigmax)^2 Exp[-sigmax/M2f])];  (* This is the tFOPT 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; a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE with the initial value a(sigmax)=aALf (the same initial value as in D=0 part, but in the D=0 part a(sigmax Exp[I fi]) and the log-derivatives run by tFOPT 5-loop)  *) 

  (* Note that, instead of simply O62/aALf, we have O62 (1/aALf - be0 Log[kappa]); this is because aALf=a(kappa sigmax), but at O62 we have to have 1/a(sigmax); we use the one-loop RGE running from Q^2=kappa sigmax) to Q^2=sigmax --> 1/a(sigmax)=(1/aALf - be0 Log[kappa]) *)
 
(* 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; O61init=-0.001; O62init=0.005; aALfinit=0.100;
(* FO approach: *)
alO4O612fitPFOMinN[Nn_]:= FindMinimum[Sum[(ReBPFOfast[M2FP[jj] Exp[I Pi/6], Nn, aALf,O4fitf,O61fitf,O62fitf] - 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,O61fitf,O62fitf] - 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,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O61fitf, O61init},{O62fitf, O62init}}];
(* aALf=alpha_s(kappa*sigmax)/Pi *)

 (* Now the dfi fi ... integration numerically, : *)
alO4O612fitPFOMinNNum[Nn_]:= FindMinimum[Sum[(ReBPFOfastNum[M2FP[jj] Exp[I Pi/6], Nn, aALf,O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfastNum[M2FP[jj] Exp[I Pi/4], Nn, aALf, O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfastNum[M2FP[jj], Nn, aALf, O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O61fitf, O61init},{O62fitf, O62init}}];
(* gives the same result as alO4O6O8fitPFOMinN, only about 5% longer time needed *)			     

  (* Now the dfi integration in D=6 term with 1/a(sigmax Exp[I fi]) 5-loop (in powers of a(sigmax)): *)
alO4O612fitPFOMinNNum5l[Nn_]:= FindMinimum[Sum[(ReBPFOfastNum5l[M2FP[jj] Exp[I Pi/6], Nn, aALf,O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfastNum5l[M2FP[jj] Exp[I Pi/4], Nn, aALf, O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfastNum5l[M2FP[jj], Nn, aALf, O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O61fitf, O61init},{O62fitf, O62init}}];


 (* Now the dfi integration in D=6 term with 1/a(sigmax Exp[I fi]) 5-loop CI(RGE)-evolved: *)
 alO4O612fitPFOMinNNum5lCI[Nn_]:= FindMinimum[Sum[(ReBPFOfastNum5lCI[M2FP[jj] Exp[I Pi/6], Nn, aALf,O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/6]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/6]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfastNum5lCI[M2FP[jj] Exp[I Pi/4], Nn, aALf, O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj] Exp[I Pi/4]])^2/delReBPexpBinAL[M2FP[jj] Exp[I Pi/4]]^2, {jj, 0, NchP}]+Sum[(ReBPFOfastNum5lCI[M2FP[jj], Nn, aALf, O4fitf,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O61fitf, O61init},{O62fitf, O62init}}];

  (*-----------------------__*)
  (* BL results, global FO approach with (O61, O62): *)
  (* Note: the results of 'NNum5l' [where 1/a(sigmax Exp[I fi]) ] at <O_6^{(2)}> runs by 5-loop FO expansion (in powers of a(sigmax)) is NOT GOOD; 5-loop CI running of a(sigmax Exp[I fi}) [in  <O_6^{(2)}>/a(sigmax Exp[I fi]) ] ('NNum5lCI') gives good results, similar to the case of one-loop running [in  <O_6^{(2)}>/a(sigmax Exp[I fi]) ] ('N', 'NNum') *)
(*
alO4O612fitPFOMinnum[3]={0.000563168, {aALf -> 0.109216, O4fitf -> 0.00126284, O61fitf -> 0.000261931, O62fitf -> -0.0000559142}};
alO4O612fitPFOMinnum[4]={1.2064*10^-7, {aALf -> 0.10544, O4fitf -> 0.0000153059, O61fitf -> -0.000642673, O62fitf -> 0.000148444}};
alO4O612fitPFOMinnum[5]={0.000304264, {aALf -> 0.10453, O4fitf -> -0.00121531, 
  O61fitf -> -0.00143785, O62fitf -> 0.000326175}};
alO4O612fitPFOMinnum[6]={0.000807655, {aALf -> 0.104682, O4fitf -> -0.00245071,  O61fitf -> -0.00205447, O62fitf -> 0.000476713}};
alO4O612fitPFOMinnum[7]={0.00120741, {aALf -> 0.105423, O4fitf -> -0.00377787, 
  O61fitf -> -0.0025026, O62fitf -> 0.000606378}};
alO4O612fitPFOMinnum[8]={0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, 
  O61fitf -> -0.00275591, O62fitf -> 0.000710728}};
alO4O612fitPFOMinnum[9]={0.00108708, {aALf -> 0.107909, O4fitf -> -0.00710709, 
  O61fitf -> -0.00278058, O62fitf -> 0.00078443}};
alO4O612fitPFOMinnum[10]={0.000509472, {aALf -> 0.109315, O4fitf -> -0.00934643,  O61fitf -> -0.00244755, O62fitf -> 0.00078981}};
alO4O612fitPFOMinnum[11]={0.0000454886, {aALf -> 0.111281, O4fitf -> -0.0129596, O61fitf -> -0.00169077, O62fitf -> 0.000721737}};
alO4O612fitPFOMinnum[12]={0.000263186, {aALf -> 0.111324, O4fitf -> -0.0159239, O61fitf -> -0.000390038, O62fitf -> 0.000441971}}; (* OK *)
   
 Clear[Nn]; Do[Module[{}, xi2FO[Nn]=alO4O612fitPFOMinnum[Nn][[1]]; aALfFO[Nn]=aALf/.alO4O612fitPFOMinnum[Nn][[2]]; alsmtauFO[Nn]=aMSext[mtauphys^2/(kappa sigmax),aALfFO[Nn]] Pi;  O4FO[Nn]=O4fitf/.alO4O612fitPFOMinnum[Nn][[2]]; O61FO[Nn]=O61fitf/.alO4O612fitPFOMinnum[Nn][[2]]; O62FO[Nn]=O62fitf/.alO4O612fitPFOMinnum[Nn][[2]]],{Nn,3,12}]; (* OK *)
   *)
  (* end of BL results, global FO approach with (O61,O62) *)
  (*----------------------------------------------------*)

  (* when d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=275.: alO4O612fitPFOMinN[8]={0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, O61fitf -> -0.00275591, O62fitf -> 0.000710728}}; alpha_s(mtau^2)=alsmtauFO[8]-0.322843;
   aaFO21[0.106507][8, -0.00275591, 0.000710728]=0.204956;
   aaFO20[0.106507][8, -0.00527267, -0.00275591, 0.000710728]=0.239921; OK  *)

  (* when d4=4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=401.38: alO4O612fitPFOMinN[8]={0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, O61fitf -> -0.00275591, O62fitf -> 0.000710728}}; alpha_s(mtau^2)=0.303406;
   aaFO21[0.0998138][8, -0.00258002, 0.000561191]=0.205148;
   aaFO20[0.0998138][8, -0.00367209, -0.00258002, 0.000561191]=0.240002; 
   *)
  
  (*----------------------------------------------------*)			        
(* tFO ({tilde FOPT}) approach: *)
(* with (O61, O62): in O62/a(sigmax Exp[I fi]) the running is one-loop, but with a(sigmax) the same as in the D=0 part *)
Clear[aALf];
O4init=0.001; O61init=-0.001; O62init=0.005; aALfinit=0.100;
alO4O612fitPtFOMinN[Nn_]:= FindMinimum[Sum[(ReBPtFO[M2FP[jj] Exp[I Pi/6], Nn, aALf,O4fitf,O61fitf,O62fitf] - 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,O61fitf,O62fitf] - 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,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O61fitf, O61init},{O62fitf, O62init}}];
(* aALf=alpha_s(kappa*sigmax)/Pi *)

  (* Do[alO4O612fitPtFOMinnum[Nn]=alO4O612fitPtFOMinN[Nn],{Nn,3,12}]; *)
  (*-----------------------__*)
  (* BL results, global tFO approach with (O61, O62): in O62/a(sigmax Exp[I fi]) the running is one-loop, but with a(sigmax) the same as in the D=0 part *)

  (*
alO4O612fitPtFOMinnum[3]={0.0000550704, {aALf -> 0.103353, O4fitf -> 0.000594932, O61fitf -> -0.000289906, O62fitf -> 0.0000639169}};
alO4O612fitPtFOMinnum[4]={0.000894446, {aALf -> 0.102959, O4fitf -> -0.00162577, O61fitf -> -0.00191124, O62fitf -> 0.000410859}};
alO4O612fitPtFOMinnum[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, 
  O61fitf -> -0.00324893, O62fitf -> 0.000772663}};
alO4O612fitPtFOMinnum[6]={0.00174365, {aALf -> 0.108768, O4fitf -> -0.00790883, 
  O61fitf -> -0.00348822, O62fitf -> 0.000975418}};
alO4O612fitPtFOMinnum[7]={0.0216534, {aALf -> 0.116036, O4fitf -> -0.0173627, 
  O61fitf -> -0.00198169, O62fitf -> 0.00111877}};
alO4O612fitPtFOMinnum[8]={0.00185682, {aALf -> 0.109561, O4fitf -> -0.0130742, 
  O61fitf -> -0.0000822622, O62fitf -> 0.000402686}};
alO4O612fitPtFOMinnum[9]={3.0397, {aALf -> 0.101047, O4fitf -> -0.00395687, 
  O61fitf -> 0.00796974, O62fitf -> -0.00110049}};
alO4O612fitPtFOMinnum[10]={0.00355735, {aALf -> 0.097154, O4fitf -> -0.00545422, O61fitf -> 0.00176413, O62fitf -> -0.000224147}};
alO4O612fitPtFOMinnum[11]={11.2221, {aALf -> 0.0848126, O4fitf -> 0.00638193, 
  O61fitf -> 0.0205593, O62fitf -> -0.00249062}};
   alO4O612fitPtFOMinnum[12]={0.00865425, {aALf -> 0.0857143, O4fitf -> -0.000304392, O61fitf -> 0.00312514, O62fitf -> -0.000428014}}; (* OK *)

 Clear[Nn]; Do[Module[{}, xi2tFO[Nn]=alO4O612fitPtFOMinnum[Nn][[1]]; aALftFO[Nn]=aALf/.alO4O612fitPtFOMinnum[Nn][[2]]; alsmtautFO[Nn]=aMSext[mtauphys^2/(kappa sigmax),aALftFO[Nn]] Pi;  O4tFO[Nn]=O4fitf/.alO4O612fitPtFOMinnum[Nn][[2]]; O61tFO[Nn]=O61fitf/.alO4O612fitPtFOMinnum[Nn][[2]]; O62tFO[Nn]=O62fitf/.alO4O612fitPtFOMinnum[Nn][[2]]],{Nn,3,12}];  (* OK *)
  *)
  (* end of BL results, global tFO approach with (O61, O62) *)
  (*---------------------------------------------------------------------_*)
  (* when d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=275. (tddtk[4, tkap]=-25.4181): alO4O612fitPtFOMinN[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, O61fitf -> -0.00324893, O62fitf -> 0.000772663}}; alpha_s(mtau^2)=alsmtautFO[5]=0.320945;
   aatFO21[0.105851][5, -0.00324893, 0.000772663]=0.204368;
   aatFO20[0.105851][5, -0.00457748, -0.00324893, 0.000772663]=0.23969; OK *)

  (* when d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=401.38 (tddtk[4, tkap]=100.962): alO4O612fitPtFOMinN[5]={0.00160019, {aALf -> 0.103349, O4fitf -> -0.00363644, O61fitf -> -0.00280385, O62fitf -> 0.000635493}}; alpha_s(mtau^2)=0.31369;
   aatFO21[0.10334889326990178][5, -0.0028038536467282884, 0.0006354932519756344]=0.20482;  note:  aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
   aatFO20[0.10334889326990178][5, -0.00363643742791981, -0.0028038536467282884, 0.0006354932519756344]=0.23987; note: aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026
   *)
   (*---------------------------------------------------------------------_*)
  

  alsmtaufFO[8]=aMSext[mtau^2/sigmax,aALfFO[8]] Pi;

Plot0FOo77= 
Plot[{ReBPexpBinAL[M2] + delReBPexpBinAL[M2], 
  ReBPexpBinAL[M2] - delReBPexpBinAL[M2], 
  ReBPexpBinAL[M2],
  ReBPFOfast[M2, 8, aALfFO[8], O4FO[8], O61FO[8], O62FO[8]]},
  {M2, M2min, 1.5}, PlotRange -> {{M2min, 1.5}, {0.50, 0.75}}, 
 AxesOrigin -> {M2min, 1.09}, 
 PlotStyle -> {{Black, Thin}, {Black, Thin}, {Black, 
    Thickness[0.003]}, {Red, Dashed, 
    Thickness[0.003]}},(*{Black,Dotted,Thickness[0.003]},{Black,DotDashed,\
Thickness[0.003]}},*)Filling -> {1 -> {2}}, Frame -> True, 
 FrameLabel -> {Superscript[Abs[M], 2][Superscript[GeV, 2]], 
   Re B[Superscript[M, 2]]}, 
 PlotLegends -> Placed[{None, None, Style["ALEPH exp", 10], 
    Style["FO (Nt=8)", 10]}, {0.25, 0.25}], 
 Epilog -> {Inset[Style["(a)", 11], Scaled[{0.50, 0.88}]], 
   Inset[Style[\[Psi] == 0, 11], Scaled[{0.70, 0.88}]], 
   Inset[Style[Subscript[\[Alpha], s], 11][
      Subsuperscript[m, \[Tau], 2]] == 
     Style[Round[alsmtaufFO[8] 10^4]/10.^4, 11], 
    Scaled[{0.77, 0.74}]], 
   Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], 
	 Scaled[{0.77, 0.60}]]}];

PlotPi6FOo77= 
Plot[{ReBPexpBinAL[M2 Exp[I Pi/6]] + delReBPexpBinAL[M2 Exp[I Pi/6]], 
      ReBPexpBinAL[M2 Exp[I Pi/6]] - delReBPexpBinAL[M2 Exp[I Pi/6]],
      ReBPexpBinAL[M2 Exp[I Pi/6]],
    ReBPFOfast[M2 Exp[I Pi/6], 8, aALfFO[8], O4FO[8], O61FO[8], O62FO[8]]},
  {M2, M2min, 1.5}, PlotRange -> {{M2min, 1.5}, {0.50, 0.75}}, 
 AxesOrigin -> {M2min, 1.09}, 
 PlotStyle -> {{Black, Thin}, {Black, Thin}, {Black, 
    Thickness[0.003]}, {Red, Dashed, 
    Thickness[0.003]}},(*{Black,Dotted,Thickness[0.003]},{Black,DotDashed,\
Thickness[0.003]}},*)Filling -> {1 -> {2}}, Frame -> True, 
 FrameLabel -> {Superscript[Abs[M], 2][Superscript[GeV, 2]], 
   Re B[Superscript[M, 2]]}, 
 PlotLegends -> Placed[{None, None, Style["ALEPH exp", 10], 
    Style["FO (Nt=8)", 10]}, {0.25, 0.25}], 
 Epilog -> {Inset[Style["(b)", 11], Scaled[{0.50, 0.88}]], 
   Inset[Style[\[Psi] == \[Pi]/6 , 11], Scaled[{0.70, 0.88}]], 
   Inset[Style[Subscript[\[Alpha], s], 11] [
      Subsuperscript[m, \[Tau], 2]] == 
     Style[Round[alsmtaufFO[8] 10^4]/10.^4, 11], Scaled[{0.77, 0.74}]], 
   Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], 
	 Scaled[{0.77, 0.60}]]}];

PlotPi4FOo77= 
 Plot[{ReBPexpBinAL[M2 Exp[I Pi/4]] + delReBPexpBinAL[M2 Exp[I Pi/4]],
    ReBPexpBinAL[M2 Exp[I Pi/4]] - delReBPexpBinAL[M2 Exp[I Pi/4]], 
   ReBPexpBinAL[M2 Exp[I Pi/4]], 
   ReBPFOfast[M2 Exp[I Pi/4], 8, aALfFO[8], O4FO[8], O61FO[8], 
    O62FO[8]]}, {M2, M2min, 1.5}, 
  PlotRange -> {{M2min, 1.5}, {0.50, 0.75}}, 
  AxesOrigin -> {M2min, 1.09}, 
  PlotStyle -> {{Black, Thin}, {Black, Thin}, {Black, 
     Thickness[0.003]}, {Red, Dashed, Thickness[0.003]}},(*{Black,
  Dotted,Thickness[0.003]},{Black,DotDashed,Thickness[0.003]}},*)
  Filling -> {1 -> {2}}, Frame -> True, 
  FrameLabel -> {Superscript[Abs[M], 2][Superscript[GeV, 2]], 
    Re B[Superscript[M, 2]]}, 
  PlotLegends -> 
   Placed[{None, None, Style["ALEPH exp", 10], 
     Style["FO (Nt=8)", 10]}, {0.25, 0.25}], 
  Epilog -> {Inset[Style["(c)", 11], Scaled[{0.50, 0.88}]], 
    Inset[Style[\[Psi] == \[Pi]/4, 11], Scaled[{0.70, 0.88}]], 
    Inset[Style[Subscript[\[Alpha], s], 11][
       Subsuperscript[m, \[Tau], 2]] == 
      Style[Round[alsmtaufFO[8] 10^4]/10.^4, 11], Scaled[{0.77, 0.74}]], 
   Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], 
	 Scaled[{0.77, 0.60}]]}];

(* Export["FigPsi0PFOo77O612d4275.pdf",Plot0FOo77];  Export["FigPsiPi6PFOo77O612d4275.pdf",PlotPi6FOo77];  Export["FigPsiPi4PFOo77O612d4275.pdf", PlotPi4FOo77]; *)
(*---------------------------------------------------------------------_*)
 (* 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] *)
(*---*)
(* Now the experimental ALEPH values for (0,0) moment a^{(0,0)} (see also the file rtauALEPH.m): *)
grs0a[0,0][sig_]:=(1/sigmax);
aexpBin77[0,0]=(Sum[ dsbin[[ii]] grs0a[0,0][sbin[[ii]]] OmegaVplA[ii],{ii,1,77}] + 2 Pi^2 fPi^2 grs0a[0,0][mPi^2]) -  1;
delaexpBin77[0,0]=(2 Pi^2) Sqrt[Sum[ grs0a[0,0][sbin[[ii]]] grs0a[0,0][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^{(0,0)}_{exp} = aexpBin77[0,0] +- delaexpBin77[0,0] *)
(*--------*)

 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) *)
(*--------------*)
(*  This is the NEW (30 Aug.2021) in FO approach with (O61, O62): *)
(* a21 and a20, this is the case when a(sigmax Exp[I fi]) in [O_6^{(2)}/a(sigmax Exp[I fi])] runs by one-loop RGE: *)
aaFO21[aALf_][NN_,O61_,O62_]:=aaFOD0fast[2,1,aALf][NN]+ (6 Pi^2/sigmax^3)*( (-2)(O61+O62/aALf) + O62 be0 (1/3) );  (* this is a^{(2,1)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE  *)
aaFO20[aALf_][NN_,O4_,O61_,O62_]:=aaFOD0fast[2,0,aALf][NN]+ (6 Pi^2) ( 2 O4/sigmax^2 + (1/sigmax^3)( (O61+O62/aALf) + O62 be0 (11/6) ) );   (* this is a^{(2,0)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE  *)

(* Note that the expressions aaFO21[aALf_][NN_,O61_,O62_] and aaFO20[aALf_][NN_,O4_,O61_,O62_] are valid only when kappa=1, i.e., when aALf=a(sigmax) *)

(* This is the case when a(sigmax Exp[I fi]) in [O_6^{(2)}/a(sigmax Exp[I fi])] runs by 5-loop RGE (i.e., 5-l CI(RGE): '5lCI'): *)
Clear[sol1]; Clear[aCI];
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]) ' *)
aaFO215lCI[aALf_][NN_,O61_,O62_]:=aaFOD0fast[2,1,aALf][NN]+ (6 Pi^2/sigmax^3)*( (-2) O61 + O62 2 (1/2/Pi) NIntegrate[(1/aCI[fi,aALf])*( (Exp[- I 2 fi] -1) + (1/2) Exp[I fi] (Exp[-I 4 fi] -1) ),{fi,-Pi,Pi}] );  (* this is a^{(2,1)} in the approach with O61 and O62 *)
aaFO205lCI[aALf_][NN_,O4_,O61_,O62_]:=aaFOD0fast[2,0,aALf][NN]+ (6 Pi^2) ( 2 O4/sigmax^2 + (1/sigmax^3) (O61 + O62 3 (1/2/Pi) NIntegrate[(1/aCI[fi,aALf])*( Exp[-I 2 fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[-I 3 fi]) ),{fi,-Pi,Pi}] ) );
  (*-----------------------------_*)
  (* Using the results at Nt=8: *)
  (*  Timing[alO4O612fitPFOMinN[8]]={22.2212, {0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, O61fitf -> -0.00275591, O62fitf -> 0.000710728}}} (alpha_s(mtau^2)=0.322843) OK *)

   (* ---> *)
   (* aaFO21[aALfFO[8]][8, O61FO[8], O62FO[8]]=0.204956;
    aaFO20[aALfFO[8]][8, O4FO[8], O61FO[8], O62FO[8]]=0.239921; OK *)

   VERY GOOD, because:
   aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
   aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026
   *)
  (*------------------------*)
  (*  This is for global FO approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ]: *)
  (* Tables of a20, a21 vs Nt in FO when aALf=a(sigmax) is VARIED ('avar'), i.e.,using alO4O612fitPFOMinnum[N] (N=3,4,...12) from global FO: *)			     
(*
 Tab20D0FOavar=Table[{NN,aaFOD0fast[2,0,aALfFO[NN]][NN]},{NN,3,12}]; (* Tab20D0FOavar={{3, 0.222904}, {4, 0.236422}, {5, 0.250341}, {6, 0.26499}, {7, 0.281496}, {8, 0.301097}, {9, 0.326499}, {10, 0.359622}, {11, 0.415128}, {12, 0.466577}} *)
 Tab20FOavar=Table[{NN,aaFO20[aALfFO[NN]][NN,O4FO[NN],O61FO[NN],O62FO[NN]]},{NN,3,12}];(* Tab20FOavar={{3, 0.240684}, {4, 0.24037}, {5, 0.240151}, {6, 0.240016}, {7, 0.239939}, {8, 0.239921}, {9, 0.239973}, {10, 0.240106}, {11, 0.240319}, {12, 0.240598}} *)
 delTab20FOavar= Table[{Nn, Tab20FOavar[[Nn - 2]][[2]] - Tab20FOavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];  (* Nn=7 -> 8 has clearly the smallest change; therefore, we take Nt=8 in FO case *)
 Tab21D0FOavar=Table[{NN,aaFOD0fast[2,1,aALfFO[NN]][NN]},{NN,3,12}]; (* Tab21D0FOavar={{3, 0.205618}, {4, 0.209891}, {5, 0.213932}, {6, 0.217695}, {7, 0.221291}, {8, 0.224652}, {9, 0.227714}, {10, 0.229605}, {11, 0.230339}, {12, 0.22504}} *)
 Tab21FOavar=Table[{NN,aaFO21[aALfFO[NN]][NN,O61FO[NN],O62FO[NN]]},{NN,3,12}]; (* Tab21FOavar={{3, 0.206854}, {4, 0.206064}, {5, 0.205514}, {6, 0.205175}, {7, 0.204988}, {8, 0.204956}, {9, 0.205083}, {10, 0.205427}, {11, 0.205929}, {12, 0.206619}} *)
delTab21FOavar= Table[{Nn, Tab21FOavar[[Nn - 2]][[2]] - Tab21FOavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];  (* Nn=7 -> 8 has clearly the samllest change; therefore, we take Nt=8 in FO case *)
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}]; 
*)

  (* Nt=8 is the best Nt in FO case *)
 (*-------*)   
(*
(*  This is for the global FO approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ]: *)
 (* Figures for global FO: a21, a20, vs Nt [aALf=a(sigmax) is VARIED ('avar')], with the (O61, O62) approach : *)
 a20FOavar= ListPlot[{Tab20D0FOavar, Tab20FOavar, ff20exp, ff20miexp, ff20plexp},PlotRange -> {{2.8, 12.2}, {0.205, 0.35}}, PlotStyle -> {Blue, Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.40, 0.85}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.77, 0.60}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"FO (D=0)", "FO"}, {0.17, 0.58}]];
a20FOavarZoom= ListPlot[{Tab20FOavar, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.237, 0.244}}, PlotStyle -> {Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.60, 0.95}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.40, 0.95}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"FO"}, {0.90, 0.95}]];

a21FOavar= ListPlot[{Tab21D0FOavar, Tab21FOavar, ff21exp, ff21miexp, ff21plexp},PlotRange -> {{2.8, 12.2}, {0.20, 0.24}}, PlotStyle -> {Blue, Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.40, 0.85}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.77, 0.60}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"FO (D=0)", "FO"}, {0.17, 0.58}]];
a21FOavarZoom= ListPlot[{Tab21FOavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.202, 0.210}}, PlotStyle -> {Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.60, 0.95}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.40, 0.95}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"FO"}, {0.90, 0.95}]];

 *)
    (* Export["a20FOavarO612d4275.pdf",a20FOavar]; Export["a21FOavarO612d4275.pdf",a21FOavar];  Export["a20FOavarO612Zd4275.pdf",a20FOavarZoom]; Export["a21FOavarO612Zd4275.pdf",a21FOavarZoom];  *)
    (* global FO: a21 and a20 suggest that Nt=8 is the most reasonable truncation (in this global FO fit approach), because the difference from Nt=7 to B=Nt=8 has the smallest absolute value, and a20 and a21 there achieve local minimum as a function of Nt; *)
    (*--------------------------*)
    (* Tables of a20(D=0), a21(D=0) vs Nt in FO when aALf=a(sigmax) is FIXED ('afxd'), i.e., aAL=aALfFO[8](=0.106507), [ alpha_s(mtau^2)=alsmtauFO[8]=0.322843 ],i.e., the coupling alpha_s taken FIXED at value as obtaiuned from global BL FO fits at Nt=8: *)
(*
 Tab20D0FOafxd=Table[{NN,aaFOD0fast[2,0,aALfFO[8]][NN]},{NN,3,12}]; (* Tab20D0FOafxd={{3, 0.213565}, {4, 0.241121}, {5, 0.260746}, {6, 0.276123}, {7, 0.289057}, {8, 0.301097}, {9, 0.313517}, {10, 0.32849}, {11, 0.347835}, {12, 0.377691}}; *)
 Tab21D0FOafxd=Table[{NN,aaFOD0fast[2,1,aALfFO[8]][NN]},{NN,3,12}]; (* Tab21D0FOafxd={{3, 0.197376}, {4, 0.213728}, {5, 0.221716}, {6, 0.225235}, {7, 0.225838}, {8, 0.224652}, {9, 0.222153}, {10, 0.219239}, {11, 0.215153}, {12, 0.211958}}; *)
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}];
 *)
(*
 (* Figures for global FO: a21(D=0), a20(D=0), vs Nt, WHERE aAL=a(sigmax) is FIXED ('aFXD'); a(sigmax)=aALfFO[8](=0.106507), [ alpha_s(mtau^2)=alsmtauFO[8]=0.322843 ]: *)
a20FOafxd= ListPlot[{Tab20D0FOafxd, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.20, 0.40}}, PlotStyle -> {Blue, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.60, 0.85}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.08, 0.95}]],Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]]== Round[10000 alsmtauFO[Nt = 8]]/10000.,13], Scaled[{0.23, 0.65}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.15, 0.50}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"FO (D=0)"}, {0.20, 0.78}]];
a21FOafxd= ListPlot[{Tab21D0FOafxd, ff21exp, ff21miexp, ff21plexp},PlotRange -> {{2.8, 12.2}, {0.18, 0.25}}, PlotStyle -> {Blue, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.70, 0.85}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.10, 0.95}]],Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]]== Round[10000 alsmtauFO[Nt = 8]]/10000.,13], Scaled[{0.23, 0.65}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.15, 0.50}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"FO (D=0)"}, {0.20, 0.78}]];

(* end of the global FO approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ] *)  

 (*---------------------------------------*)

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) );
  (*--------------*)
(*  This is the global tFO approach with (O61, O62) : *)
  (* a21 and a20, this is the case when a(sigmax Exp[I fi]) in [O_6^{(2)}/a(sigmax Exp[I fi])] runs by one-loop RGE: *)
  (* This is NEW: *)
aatFO21[aALf_][NN_,O61_,O62_]:=aatFOD0[2,1,aALf][NN]+ (6 Pi^2/sigmax^3)*( (-2)(O61+O62/aALf) + O62 be0 (1/3) );  (* this is a^{(2,1)} in the tFO approach with O61 and O62, and a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE  *)
aatFO20[aALf_][NN_,O4_,O61_,O62_]:=aatFOD0[2,0,aALf][NN]+ (6 Pi^2) ( 2 O4/sigmax^2 + (1/sigmax^3)( (O61+O62/aALf) + O62 be0 (11/6) ) );   (* this is a^{(2,0)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE  *)
(* end of a21 and a20 in global tFO approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ] *)

(* Note that the expressions aatFO21[aALf_][NN_,O61_,O62_] and aatFO20[aALf_][NN_,O4_,O61_,O62_] are valid only when kappa=1, i.e., when aALf=a(sigmax) *)
(*-----------------------------_*)
(* Using the results at Nt=5 (and d4bar(=ddtk[4, Exp[tKK]])=275.): *)
  (*  alO4O612fitPtFOMinnum[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, O61fitf -> -0.00324893, O62fitf -> 0.000772663}} OK *)

  (* ---> *)
(* aatFO21[aALftFO[5]][5, O61tFO[5], O62tFO[5]]=0.204368;
 aatFO20[aALftFO[5]][5, O4tFO[5], O61tFO[5], O62tFO[5]]=0.23969; OK
 
   GOOD, because:
   aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
   aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026
 *)
  (*------------------------*)
(*  This is for the global tFO approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ]: *)
  (* Tables of a20, a21 vs Nt in tFO when aALf=a(sigmax) is VARIED ('avar'), i.e.,using alO4O612fitPtFOMinnum[N] (N=3,4,...12) from global tFO: *)			     
(*
 Tab20D0tFOavar=Table[{NN,aatFOD0[2,0,aALftFO[NN]][NN]},{NN,3,12}]; (* Tab20D0tFOavar={{3, 0.229886}, {4, 0.254364}, {5, 0.289315}, {6, 0.333717}, {7, 0.469039}, {8, 0.424269}, {9, 0.313244}, {10, 0.327514}, {11, 0.181898}, {12, 0.256028}} *)
 Tab20tFOavar=Table[{NN,aatFO20[aALftFO[NN]][NN,O4tFO[NN],O61tFO[NN],O62tFO[NN]]},{NN,3,12}]; (*  Tab20tFOavar={{3, 0.240471}, {4, 0.239985}, {5, 0.23969}, {6, 0.239878}, {7, 0.239861}, {8, 0.240937}, {9, 0.233344}, {10, 0.241161}, {11, 0.226835}, {12, 0.241626}} *)
 delTab20tFOavar= Table[{Nn, Tab20tFOavar[[Nn - 2]][[2]] - Tab20tFOavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];
 Tab21D0tFOavar=Table[{NN,aatFOD0[2,1,aALftFO[NN]][NN]},{NN,3,12}];  (* Tab21D0tFOavar={{3, 0.207954}, {4, 0.215482}, {5, 0.224658}, {6, 0.232449}, {7, 0.244095}, {8, 0.226065}, {9, 0.177595}, {10, 0.205586}, {11, 0.134793}, {12, 0.200042}} *)
 Tab21tFOavar=Table[{NN,aatFO21[aALftFO[NN]][NN,O61tFO[NN],O62tFO[NN]]},{NN,3,12}]; (* Tab21tFOavar={{3, 0.206311}, {4, 0.205095}, {5, 0.204368}, {6, 0.204858}, {7, 0.205032}, {8, 0.207494}, {9, 0.191128}, {10, 0.208062}, {11, 0.177269}, {12, 0.209256}} *) 
delTab21tFOavar= Table[{Nn, Tab21tFOavar[[Nn - 2]][[2]] - Tab21tFOavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}]; 
*)

 (*-------*)   
(*
(*  This is for the global tFO approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ]: *)
 (* Figures for global tFO: a21, a20, vs Nt [aALf=a(sigmax) is VARIED ('avar')], with the global tFO approach with (O61, O62): *)
a20tFOavar= ListPlot[{Tab20D0tFOavar, Tab20tFOavar, ff20exp, ff20miexp, ff20plexp},PlotRange -> {{2.8, 12.2}, {0.205, 0.35}}, PlotStyle -> {Blue, Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.40, 0.85}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.77, 0.60}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"tFO (D=0)", "tFO"}, {0.17, 0.58}]];
a20tFOavarZoom= ListPlot[{Tab20tFOavar, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.232, 0.244}}, PlotStyle -> {Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.60, 0.95}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.40, 0.95}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"tFO"}, {0.90, 0.95}]];

a21tFOavar= ListPlot[{Tab21D0tFOavar, Tab21tFOavar, ff21exp, ff21miexp, ff21plexp},PlotRange -> {{2.8, 12.2}, {0.20, 0.24}}, PlotStyle -> {Blue, Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.40, 0.85}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.40, 0.95}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"tFO (D=0)", "tFO"}, {0.17, 0.58}]];
a21tFOavarZoom= ListPlot[{Tab21tFOavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.198, 0.210}}, PlotStyle -> {Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.60, 0.95}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.08, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.40, 0.95}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"tFO"}, {0.90, 0.95}]];

 *)
    (* Export["a20tFOavarO612d4275.pdf",a20tFOavar]; Export["a21tFOavarO612d4275.pdf",a21tFOavar];  Export["a20tFOavarO612Zd4275.pdf",a20tFOavarZoom]; Export["a21tFOavarO612Zd4275.pdf",a21tFOavarZoom];  *)
    (* global tFO: a21 and a20 suggest that Nt=5 is the most reasonable truncation (in this global tFO fit approach), because there we get local minimum *)
    (*-------------------*)
    (* Tables of a20(D=0), a21(D=0) vs Nt in tFO when aALf=a(sigmax) is FIXED ('afxd'), i.e., aAL=aALftFO[5](=0.105851), [ alpha_s(mtau^2)=alsmtautFO[5]=0.320945],i.e., the coupling alpha_s taken FIXED at value as obtaiuned from global BL tFO fits at Nt=5: *)
(*
Tab20D0tFOafxd=Table[{NN,aatFOD0[2,0,aALftFO[5]][NN]},{NN,3,12}]; (* Tab20D0tFOafxd={{3, 0.241427}, {4, 0.271948}, {5, 0.289315}, {6, 0.3062}, {7, 0.324714}, {8, 0.360754}, {9, 0.388832}, {10, 0.61617}, {11, 0.372224}, {12, 3.96521}}; *)
Tab21D0tFOafxd=Table[{NN,aatFOD0[2,1,aALftFO[5]][NN]},{NN,3,12}]; (* Tab21D0tFOafxd={{3, 0.2177}, {4, 0.228198}, {5, 0.224658}, {6, 0.220728}, {7, 0.21373}, {8, 0.214066}, {9, 0.174309}, {10, 0.298048}, {11, -0.347865}, {12, 3.03195}}; *)
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}];
 *)
(*
 (* Figures for global tFO: a21(D=0), a20(D=0), vs Nt, WHERE aAL=a(sigmax) is FIXED ('aFXD'); a(sigmax)=aALftFO[5](=0.104501), [ alpha_s(mtau^2)=alsmtautFO[5]= ]: *)
a20tFOafxd= ListPlot[{Tab20D0tFOafxd, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.20, 0.42}}, PlotStyle -> {Blue, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.70, 0.90}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.08, 0.95}]], Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]] == Round[10000 alsmtautFO[5]]/10000., 13], Scaled[{0.75, 0.45}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.80, 0.35}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"tFO (D=0)"}, {0.77, 0.58}]];
a21tFOafxd= ListPlot[{Tab21D0tFOafxd, ff21exp, ff21miexp, ff21plexp},PlotRange -> {{2.8, 12.2}, {0.18, 0.30}}, PlotStyle -> {Blue, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.50, 0.90}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.08, 0.95}]],Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]]== Round[10000 alsmtautFO[5]]/10000.,13], Scaled[{0.20, 0.65}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.12, 0.53}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"tFO (D=0)"}, {0.17, 0.78}]];

    (*------*)

  (* end of the global tFO approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ] *)  

 (*---------------------------------------*)
(*---------------------------------------------*)
(*---------------------------------------------*)
(* moment a^{(0,0}}, FO: *)			     
cAg[0,0][kp_]:=(1/(2 Pi)) Re[ (-I)^kp NIntegrate[ th^kp ( (1+Exp[I th]) ),{th, -Pi,Pi}] ];  
(* kp = 0, 1, 2, ... *)
coeffaaFO[0,0][NN_]:=Sum[tddtk[nn,tkap] Binomial[nn + kk, kk] be0^kk coeffext[nn + kk + 1, NN] cAg[0,0][kk], {nn, 0, NN - 1}, {kk, 0, NN - 1 - nn}];
(* coeffaaFO[0,0][NN] is the coefficient at x^NN in the expansion of aa[0,0] *)
aaFOD0fast[0,0,aALf_][NN_]:=Sum[coeffaaFO[0,0][nn] x^nn,{nn,1,NN}]/.{x->aALf};
aaFOfast[0,0,aALf_][NN_]:=aaFOD0fast[0,0,aALf][NN];
(* moments a^{(0,0}}, tFO: *)				     
coeffaatFO[0,0][nn_]:=Sum[tddtk[nn-kk,tkap] be0^kk  cAg[0,0][kk] Binomial[nn,kk],{kk,0,nn}];
aatFOD0[0,0,aALf_][NN_]:=Sum[coeffaatFO[0,0][nn-1] taprov[nn,x],{nn,1,NN}]/.{x->aALf};			     
 (*---------------------------------------*)
 (*global FO with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ]: *)
 (* The above behaviour of a^{(2,1)}(aAL,O61,O62) and a^{(2,0)}(aAL,O4,O61,O62) suggests that the truncation order Nt in the FO global approach (i.e., the approach of minimization with alpha_s variable) is Nt=8. *)

 (*------*) 
 (* Now the estimate of EXPERIMENTAL uncertainty in this global approach with double-pinched FO Borel-Laplace (with Nt=8): *)

(* VBBexp[Mal2,Mbe2]= < (ReB[Mal2]-<ReB[Mbe2]>) (ReB[Mal2]-<ReB[Mbe2]>) > *)
 VBBexpFAux[Mal2_,Mbe2_]:=(2 Pi^2)^2 Sum[(1-sbin[[ii]]/sigmax)^2 (1-sbin[[jj]]/sigmax)^2 Re[Exp[-sbin[[ii]]/Mal2]/(Mal2)] Re[Exp[-sbin[[jj]]/Mbe2]/(Mbe2)]  dsbin[[ii]] dsbin[[jj]] (  covVV[ii,jj]+covAA[ii,jj]+covVA[ii,jj] + covAV[ii,jj]),{ii,1,NNAL},{jj,1,NNAL}];
(* Note: VBBexpFAux[M2,M2] = delReBPexpBinAL[M2]^2, as it should be *)
(* Further, automatically we have: VBBexpFAux[Mal2,Mbe2]=VBBexpFAux[Mbe2,Mal2] *)
(* VBBexpF[Mal2_,Mbe2_]:=If[Mal2<=Mbe2,VBBexpFAux[Psi,Mal2,Mbe2],VBBexpFAux[Psi,Mbe2,Mal2]]; this is probably not necessary; anyway, mal2, Mbe2 will be complex in general *)

(* 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 *)

(* For delO4 perhaps only Psi=Pi/6 points: Nm=2; Nd=Nm;  M2cov1F[ii_]:=0.90+(ii-1)*0.60/(Nm-1); Do[M2covF[jj]=M2cov1F[jj] Exp[I Pi/6],{jj,1,Nm}];  *)
(* Or Psi=0 & Pi/6: Nm=2; Nd=4; M2cov1F[ii_]:=0.90+(ii-1)*0.60/(Nm-1); Do[M2covF[jj]=M2cov1F[jj],{jj,1,Nm}]; Do[M2covF[jj+Nm]=M2cov1F[jj] Exp[I Pi/6],{jj,1,Nm}]; *)
(* or all Psi=0, Pi/6, Pi/4, but only with one (middle) point for each Psi: Nd=3; M2covF[1]=(0.9+1.5)/2;  M2covF[2]=M2covF[1] Exp[I Pi/6]; M2covF[3]=M2covF[1] Exp[I Pi/4]; *)

 
VBBexpMatr=Table[If[ii<=jj,VBBexpFAux[M2covF[ii],M2covF[jj]],VBBexpFAux[M2covF[jj],M2covF[ii]]],{ii,1,Nd},{jj,1,Nd}];
(* With such construction we ensure that VBBexpMatrPsi is exactly symmetric, as it should be; this allows then for a more stable investion of the matrix *)
VBBexpMatrRes=10^14 VBBexpMatr;
(* NOte, for Nd=9: 10^6 VBBexpMatr has all elements ~1, but Det[10^6 VBBexpMatr]=
 -1.55766*10^-62, practically zero, a degenerate system; *)

(* VBBexpMatr[[i,j]] = VBBexpF[M2covF[ii],M2covF[jj]], this is the Nd x Nd covariance matrix for Re BTs *)

IVBBexpMatrRes= LinearSolve[VBBexpMatrRes, IdentityMatrix[Nd]];
IVBBexpMatr=10^14 IVBBexpMatrRes;
(* IVBBexpMatr[[i,j]] is inverse of the ReBT-covariance matrix, NdxNd matrix again, for BTs *)

(* Instabilities in the Nd=9 case: IVBBexpMatr.VBBexpMatr not=IdentityMatrix[Nd] (some elements are ~10^-5-10^-4) and Det[VBBexpMatr]=app 3 10^-60 very small; however, when Nd=3, then IVBBexpMatr.VBBexpMatr=IdentityMatrix[Nd] to ~10^12 precision;  *)  
(*---*)
 (* Now estimation of delO4 with the (NdxNd) covariance matrix  ('varcov')
  (and using FO Nt=8(=Nt0) global analysis for theoretical values): *)
Nt0=8;
chi2o3O4[O4f_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALfFO[Nt0],O4f,O61FO[Nt0],O62FO[Nt0]] - ReBPexpBinAL[M2covF[jj]]) IVBBexpMatr[[jj,kk]] (ReBPFOfast[M2covF[kk],Nt0,aALfFO[Nt0],O4f,O61FO[Nt0],O62FO[Nt0]] - ReBPexpBinAL[M2covF[kk]]),{jj,1,Nd},{kk,1,Nd}];
O4FOvarcovfitMin=FindMinimum[chi2o3O4[O4fitf],{O4fitf,O4FO[Nt0]+0.0001}];
O4FOvarcovfit=O4fitf/.O4FOvarcovfitMin[[2]];
Xi2O4FOvarcovo77=O4FOvarcovfitMin[[1]]; 
 delO4FOcovpl= delO4f/.FindRoot[chi2o3O4[O4FOvarcovfit+delO4f]==Xi2O4FOvarcovo77+1,{delO4f, 1.6 10^-6}];
 delO4FOcovmi= delO4f/.FindRoot[chi2o3O4[O4FOvarcovfit+delO4f]==Xi2O4FOvarcovo77+1,{delO4f, 1.6 10^-6}];
delO4FOcov=Max[Abs[delO4FOcovpl],Abs[delO4FOcovmi]];
(* for Nm=2 (Nd=6) we get delO4FOcov =app 1.64*10^-6; for Nm=3 (Nd=6), but possible instabilities *)
(* For Nd=3  (or all Psi=0, Pi/6, Pi/4, but only with one (middle) point for each Psi: Nd=3; M2covF[1]=(0.9+1.5)/2;  M2covF[2]=M2covF[1] Exp[I Pi/6]; M2covF[3]=M2covF[1] Exp[I Pi/4]): delO4FOcov =app 6.2 10^-6; so, rather unstable (increased) result when Nd=9 -> Nd=3 *)
(*----------------*)
(*----------------*)
(* Now we will vary O4, O61, O62, alpha_s(sigmax) in UNCORRELATED chi ('var'), to get the EXPERIMENTAL uncertainties (in the global FO approach with (O61,O62)): *)
(* Note: VBBexpFAux[M2,M2] = delReBPexpBinAL[M2]^2, as it should be *)

(* 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 *)


Nt0=8;
(* Note: alO4O612fitPFOMinnum[8]={0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, O61fitf -> -0.00275591, O62fitf -> 0.000710728}}; alsmtauFO[8]=0.322843; *)
chi2o2O4[O4f_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALfFO[Nt0],O4f,O61FO[Nt0],O62FO[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; OK; note that O4FO[8]=-0.00527267 by global analysis (and O4FOvarfit=-0.00527235, very close to O4FO[8] as it should be)  *)
(* I will take uncorrelated uncertainty,  Nm=2 (Nd=6) delO4FO=0.0000304, to avoid possible subestimations (the correlated method gives very small uncertainties) *)

chi2o2O61[O61f_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALfFO[Nt0],O4FO[Nt0],O61f,O62FO[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O61FOvarfitMin=FindMinimum[chi2o2O61[O61fitf],{O61fitf,O61FO[Nt0]}];
O61FOvarfit=O61fitf/.O61FOvarfitMin[[2]];
Xi2O61FOvaro77=O61FOvarfitMin[[1]];
 delO61FO= delO61f/.FindRoot[chi2o2O61[O61FOvarfit+delO61f]==Xi2O61FOvaro77+1,{delO61f, 3. 10^-5}];
(* Nm=2; Nd=6; delO61FO=0.0000451; OK *)

chi2o2O62[O62f_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALfFO[Nt0],O4FO[Nt0],O61FO[Nt0],O62f] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O62FOvarfitMin=FindMinimum[chi2o2O62[O62fitf],{O62fitf,O62FO[Nt0]}];
O62FOvarfit=O62fitf/.O62FOvarfitMin[[2]];
Xi2O62FOvaro77=O62FOvarfitMin[[1]];
 delO62FO= delO62f/.FindRoot[chi2o2O62[O62FOvarfit+delO62f]==Xi2O62FOvaro77+1,{delO62f, 3. 10^-6}];
(* Nm=2; Nd=6; delO62FO=0.0000050; OK *)

chi2o2aAL[aALf_]:=Sum[(ReBPFOfast[M2covF[jj],Nt0,aALf,O4FO[Nt0],O61FO[Nt0],O62FO[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.000102; OK *)
(*-------*)
(* Together, the estimates of EXPERIMENTAL uncertainties in the FO global approach (with Nt=8) are (Nm=2; Nd=6), FO global:  *)
(* delO4FO= 0.0000304; delO61FO=0.000045; delO62FO=0.000005; delaALFO=0.000102(=delta(a(sigmax))) *)
(* --> *)
(*
 alpha_s(mtau^2)=0.322843+{^{+0.00030}_{-0.00030}}_{exp}; OK 
10^3 O4=-5.27+{^{+0.03}_{-0.03}}_{exp} =app -5.3+-0.0_{exp}; OK
10^3 O61=-2.76+{{+0.04}_{-0.04}}_{exp} =app -2.8+-0.0_{exp}; OK
10^4 O62=+7.11+{^{+0.05}_{-0.05}}_{exp} =app +7.1 +- 0.1_{exp}; OK *)
 
(* The precise central value result of the FO global fit:
alO4O612fitPFOMinnum[8]={0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, O61fitf -> -0.00275591, O62fitf -> 0.000710728}}; OK *)
(*---------*)
(* Variation of kappa, FO global (at Nt=8): *)
(* kappa=2.0: *)
(* alO4O612fitPFOMinN[8]={0.000746938, {aALf -> 0.0877519, O4fitf -> -0.00304484, O61fitf -> -0.00229023, O62fitf -> 0.000511572}}; alsmtauFO[8]=0.320259; OK *)

(* kappa=1.0 (the central case of the global fit): *)
alO4O612fitPFOMinnum[8]={0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, O61fitf -> -0.00275591, O62fitf -> 0.000710728}}; alsmtauFO[8]=0.322843; OK *)

(* kappa=0.5: *)
(* alO4O612fitPFOMinnum[8]={0.0528115, {aALf -> 0.151719, O4fitf -> -0.0163689, 61fitf -> -0.00131501, O62fitf -> 0.00104236}}; alsmtauFO[8]=0.345762; OK, this is STRONG variation of alsmtauFO[8]=alpha_s(mtau^2)  *)
(* Therefore, I will do an alternative: kappa=2/3 (instead of 0.5) *)
(* kappa=2/3: *)
(* alO4O612fitPFOMinnum[8]={0.00135843, {aALf -> 0.125285, O4fitf -> -0.00849201, O61fitf -> -0.00286953, O62fitf -> 0.000915006}}; alsmtauFO[8]=0.329881; OK *)

(* Note that aALf=a(kappa*sigmax) *)
(* --> thus the variations: (kappa=1+{^{+1}_{-1/3}}) *)
(*    alpha_s(mtau^2)=0.322843+{^{-0.002584}_{+0.007038}}_{kap}; OK 
      10^3 O4= -5.3+{^{+2.3}_{-3.2}}_{kap};  OK
      10^3 O61=-2.8+{^{+0.5}_{-0.1}}_{kap};  OK
      10^4 O62=+7.1+{^{-2.0}_{+2.1}}_{kap}; OK
 *)
(*-----------*)
(* Variation of d4, FO global (NOTE, at Nt=8 now) d4=338.19 +- 63.19: *)

(* when ddtk[4, tkap]=275.+63.19=338.19: 
alO4O612fitPFOMinnum[8]={0.00108794, {aALf -> 0.102946, O4fitf -> -0.00437326`, O61fitf -> -0.00267426, O62fitf -> 0.000628783}}; alpha_s(mtau^2)=0.31252; OK *)

(* when the central case d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=275. (tddtk[4, tkap]=-25.4181):
alO4O612fitPFOMinN[8]={0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, O61fitf -> -0.00275591, O62fitf -> 0.000710728}}; alpha_s(mtau^2)=0.322843; OK *)

(* when ddtk[4, tkap]=ddtk[4,Exp[tKK]]=275.-63.19=211.81: 
 alO4O612fitPFOMinN[8]={0.00134764, {aALf -> 0.109242, O4fitf -> -0.00595179, O61fitf -> -0.00273685, O62fitf -> 0.000762533}}; ; alpha_s(mtau^2)=0.330744; OK*)

(* alpha_s(mtau^2)=0.322843+{^{-0.010323}}_{+0.007901}}_{d4}; OK 
 10^3 O4= -5.3+{^{+0.9}_{-0.7}}_{d4}; OK
 10^3 O61=-2.8+{^{+0.1}_{0.0}}_{d4}; OK
 10^4 O62=+7.1+{^{-0.8}_{+0.5}}_{d4}; OK *) 

(*--------------*)
(* Variation of Nt(truncation order, Nt=8+-1 and Nt=8 +-2, Nt=8-3), FO global: *)
(*
 alO4O612fitPFOMinnum[10]={0.000509472, {aALf -> 0.109315, O4fitf -> -0.00934643, O61fitf -> -0.00244755, O62fitf -> 0.00078981}}; alsmtauFO[10]=0.330954; OK
 alO4O612fitPFOMinnum[9]={0.00108708, {aALf -> 0.107909, O4fitf -> -0.00710709, O61fitf -> -0.00278058, O62fitf -> 0.00078443}}; alsmtauFO[9]=0.326896; OK
 central case:
 alO4O612fitPFOMinnum[8]={0.00131593, {aALf -> 0.106507, O4fitf -> -0.00527267, O61fitf -> -0.00275591, O62fitf -> 0.000710728}}; alsmtauFO[8]=0.322843; OK
 alO4O612fitPFOMinnum[7]={0.00120741, {aALf -> 0.105423, O4fitf -> -0.00377787, O61fitf -> -0.0025026, O62fitf -> 0.000606378}}; alsmtauFO[7]=0.319705; OK
 alO4O612fitPFOMinnum[6]={0.000807655, {aALf -> 0.104682, O4fitf -> -0.00245071, O61fitf -> -0.00205447, O62fitf -> 0.000476713}}; alsmtauFO[6]=0.317558; OK
 alO4O612fitPFOMinnum[5]={0.000304264, {aALf -> 0.10453, O4fitf -> -0.00121531, O61fitf -> -0.00143785, O62fitf -> 0.000326175}}; alsmtauFO[5]=0.317117; OK
 *)

(*  alpha_s(mtau^2)=0.322843+{^{+0.008111}_{-0.005726}}_{Nt}; OK 
 10^3 O4= -5.3+{^{-4.0}_{+4.1}}_{Nt}; OK
 10^3 O61=-2.8+{^{+0.4}_{+0.7}}_{Nt}; OK
 10^3 O61=-2.8+{^{+0.4}_{+1.4}}_{Nt}; OK
 10^4 O62=+7.1+{^{+0.8}_{-3.8}}_{Nt}; OK
 (* Note: Nt=8^{+2}_{-3} *)

(*-----------------------------------------*)
 (* Altogether, the FO global fit approach with (O61,O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ] gives: *)
(* FO
 Nt=8;  (kap=1+{^{+1}_{-1/3}}; d4=275+-63.19; Nt=8^{+2}_{-3}:
 alpha_s(mtau^2)=0.322843 + {^{+0.00030}_{-0.00030}}_{exp} + {^{-0.002584}_{+0.007038}}_{kap} + {^{-0.010323}}_{+0.007901}}_{d4} + {^{+0.008111}_{-0.005726}}_{Nt} = 0.32284 +{^{+0.013336}_{-0.01209}} =app 0.3228 +{^{+0.0134}_{-0.0121}} =app 0.323+{^{+0.013}_{-0.012}}; OK
 10^3 O4=-5.3 +- 0.0_{exp} + {^{+2.3}_{-3.2}}_{kap} + {^{+0.9}_{-0.7}}_{d4} + {^{-4.0}_{+4.1}}_{Nt}  = -5.3 + {^{+4.8}_{-5.2}}; OK
 10^3 O61=-2.8 +- 0.0_{exp} + +{^{+0.5}_{-0.1}}_{kap} + {^{+0.1}_{-0.0}}_{d4} + {^{+0.4}_{+1.4}}_{Nt} =-2.8 + {^{+1.5}_{-0.1}}; OK
 10^4 O62=+7.1 +- 0.1_{exp} + {^{-2.0}_{+2.1}}_{kap} + {^{-0.8}_{+0.5}}_{d4} + {^{+0.8}_{-3.8}}_{Nt} =+7.1 + {^{+2.3}_{-4.4}}; OK
 OK *) 

  (* Checked twice, OK *) 

  (* End of FO global fit results [with O61 and O62, and a(*Q^2) in O62/a(Q^2) runs by one-loop RGE from a(sigmax)=aALfFO[Nt] to a(sigmax Exp[I fi]) ], OK *)
(*------------------------------------*)
(*------------------------------------*) 
(*---------------------------------------*)
 (* global tFO: The above behaviour of a^{(2,1)}(aAL,O61,O62) and a^{(2,0)}(aAL,O61,O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax)=aALftFO[Nt] to a(Q^2=sigmax Exp[I fi]) ] suggests that the truncation order Nt=5 in the tFO global approach (i.e., the approach of minimization with alpha_s variable) is Nt=5. *)

 (* the central value of the tFO global fit result, at Nt=5, is: *)
 (* alO4O612fitPtFOMinnum[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, O61fitf -> -0.00324893, O62fitf -> 0.000772663}}; alsmtautFO[5]=0.320945 (ddtk[4,Exp[tKK]]=275.) OK *)
 
  (* Now the estimate of the EXPERIMENTAL uncertainty in this global tFO approach with double-pinched tFO Borel-Laplace (with Nt=5): *)

(* 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 estimation of delO4 with the (NdxNd) covariance matrix  ('varcov')
  (and using tFO Nt=5 global analysis for theoretical values): *)
Nt0=5;
(* Note: alO4O612fitPtFOMinnum[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, O61fitf -> -0.00324893, O62fitf -> 0.000772663}}; alsmtautFO[5]=0.320945 (ddtk[4,Exp[tKK]]=275.) OK *)

(* We vary O4 in UNCORRELATED chi ('var'): *)
chi2o2O4[O4f_]:=Sum[(ReBPtFO[M2covF[jj],Nt0,aALftFO[Nt0],O4f,O61tFO[Nt0],O62tFO[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O4tFOvarfitMin=FindMinimum[chi2o2O4[O4fitf],{O4fitf,O4tFO[Nt0]}];
O4tFOvarfit=O4fitf/.O4tFOvarfitMin[[2]];
Xi2O4tFOvaro77=O4tFOvarfitMin[[1]];                       
 delO4tFO= delO4f/.FindRoot[chi2o2O4[O4tFOvarfit+delO4f]==Xi2O4tFOvaro77+1,{delO4f, 3. 10^-5}];
(* --> for Nm=2 (Nd=6) we obtain: delO4tFO=0.0000304; OK; note that O4tFO[5]=-0.00457748, by global tFO analysis, very close to O4tFOvarfit=-0.00457721 as it should be *)

chi2o2O61[O61f_]:=Sum[(ReBPtFO[M2covF[jj],Nt0,aALftFO[Nt0],O4tFO[Nt0],O61f,O62tFO[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O61tFOvarfitMin=FindMinimum[chi2o2O61[O61fitf],{O61fitf,O61tFO[Nt0]}];
O61tFOvarfit=O61fitf/.O61tFOvarfitMin[[2]];
Xi2O61tFOvaro77=O61tFOvarfitMin[[1]];
 delO61tFO= delO61f/.FindRoot[chi2o2O61[O61tFOvarfit+delO61f]==Xi2O61tFOvaro77+1,{delO61f, 3. 10^-5}];
(* Nm=2; Nd=6; delO61tFO=0.0000451; OK *)

chi2o2O62[O62f_]:=Sum[(ReBPtFO[M2covF[jj],Nt0,aALftFO[Nt0],O4tFO[Nt0],O61tFO[Nt0],O62f] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O62tFOvarfitMin=FindMinimum[chi2o2O62[O62fitf],{O62fitf,O62tFO[Nt0]}];
O62tFOvarfit=O62fitf/.O62tFOvarfitMin[[2]];
Xi2O62tFOvaro77=O62tFOvarfitMin[[1]];
delO62tFO= delO62f/.FindRoot[chi2o2O62[O62tFOvarfit+delO62f]==Xi2O62tFOvaro77+1,{delO62f, 3. 10^-5}];  (* O62tFOvarfit=0.000772658; O62tFO[5]=0.000772663; very close to each other, as it should be *)
(* Nm=2; Nd=6; delO8tFO=4.923*10^-6; OK *)

chi2o2aAL[aALf_]:=Sum[(ReBPtFO[M2covF[jj],Nt0,aALf,O4tFO[Nt0],O61tFO[Nt0],O62tFO[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
aALtFOvarfitMin=FindMinimum[chi2o2aAL[aALfitf],{aALfitf,aALftFO[Nt0]}];
aALtFOvarfit=aALfitf/.aALtFOvarfitMin[[2]];
Xi2aALtFOvaro77=aALtFOvarfitMin[[1]];
(* aALtFOvarfit=0.105853; aALftFO[5]=0.105851; very close to each other, as it should be *)
 delaALtFO= delaALf/.FindRoot[chi2o2aAL[aALtFOvarfit+delaALf]==Xi2aALtFOvaro77+1,{delaALf, 3. 10^-5}];
(* Nm=2; Nd=6; Nt0=6; delaALtFO=0.000118027; OK *)
(*-------*)
(* Together, the estimates of EXPERIMENTAL uncertainties in the tFO global approach (with Nt=5) are (Nm=2; Nd=6), in tFO global:  *)
(* delO4tFO= 0.0000304; delO61tFO=0.0000451; delO62tFO=4.92 10^-6; delaALtFO=0.000118027; (=delta(a(sigmax)))  OK  --> *)

(* alpha_s(mtau^2)=aMSext[mtau^2/sigmax, aALftFO[5]] Pi  = alsmtautFO[5]=0.320945, the central value, OK *)


(*
alpha_s(mtau^2)=0.320945+{^{+0.000342}_{-0.000342}}_{exp}; OK 
10^3 O4= -4.58+{^{+0.03}_{-0.03}}_{exp} =app -4.6+{^{+0.0}_{-0.0}}_{exp}; OK
10^3 O61=-3.25+{^{+0.05}_{-0.05}}_{exp} =app -3.2+{^{+0.1}_{-0.1}}_{exp}; OK
10^4 O62=+7.73+{^{+0.05}_{-0.05}}_{exp} =app +7.7+{^{+0.1}_{-0.1}}_{exp}; OK *)

(* The precise central value result of the tFO global fit:
 alO4O612fitPtFOMinnum[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, O61fitf -> -0.00324893, O62fitf -> 0.000772663}}; alsmtautFO[5]=0.320945 (ddtk[4,Exp[tKK]]=275.) OK *)
(*----------*)
(* Variation of kappa, tFO global: *)
(* kappa=2.0: 
 alO4O612fitPtFOMinN[5]={0.000553921, {aALf -> 0.0870384, O4fitf -> -0.00168137,O61fitf -> -0.00186779, O62fitf -> 0.000392021}}; alpha_s(mtau^2)=0.317121; OK
 kappa=1.0 (the central case of the global fit): 
 alO4O612fitPtFOMinnum[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, O61fitf -> -0.00324893, O62fitf -> 0.000772663}}; alsmtautFO[5]=0.320945; OK
 kappa=0.5: 
 alO4O612fitPtFOMinN[5]={2.62648, {aALf -> 0.130615, O4fitf -> -0.00315874, O61fitf -> 0.0044223, O62fitf -> -0.000471432}}; alpha_s(mtau^2)=0.31303; OK *)
(* Due to too large variation of alpha_s in FU when kappa=1 -> 0.5, we will use kappa_min=2/3 (instead of kappa_min=0.5): *)
kappa=2/3:
  alO4O612fitPtFOMinN[5]={0.237291, {aALf -> 0.130744, O4fitf -> -0.0113959, O61fitf -> -0.00319109, O62fitf -> 0.00125953}}; alpha_s(mtau^2)=0.34106; very large (?) OK *)

(* --> thus the variations  (kap=1.+{^{+1}_{-1/3}}: *)
alpha_s(mtau^2)=0.320945+{^{-0.003824}_{+0.020115}}_{kap}; OK
10^3 O4= -4.6+{^{+2.9}_{-6.8}}_{kap}; OK
10^3 O61=-3.2+{^{+1.3}_{-0.0}}_{kap}; OK
10^4 O62=+7.7+{^{-3.8}_{+4.9}}_{kap}; OK
*)
 (*--------*)
(* Variation of d4, tFO global (note, since out Nt is Nt=5, this mean that the corresponding variations of d5, d6, etc. do not have any effect here): *)
 (* d4=338.19 +- 63.19 (td4=37.772+-63.10, not a small variation for td4(!)): *)
(*
When ddtk[4, tkap]=275+63.19=338.19 (tddtk[4, tkap]=37.772): 
 alO4O612fitPtFOMinnum[5]={0.00213488, {aALf -> 0.104501, O4fitf -> -0.0040562, O61fitf -> -0.00300781, O62fitf -> 0.000696829}}; alpha_s(mtau^2)=0.317033; OK
When central d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=275. (tddtk[4, tkap]=-25.4181):
 alO4O612fitPtFOMinnum[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, O61fitf -> -0.00324893, O62fitf -> 0.000772663}};  alsmtautFO[5]=0.320945; OK
 When ddtk[4, tkap]=275-63.19=211.809 (tddtk[4, tkap]=-88.6087):
  alO4O612fitPtFOMinnum[5]=0.0040513, {aALf -> 0.107492, O4fitf -> -0.00525721, O61fitf -> -0.00354533, O62fitf -> 0.000871015}};  alpha_s(mtau^2)=0.325691; OK
  *)
 (* Note that here in this approach we see that td4 does not just have minor variation, it varies td4=-25.4181+- 63.10, which is roughly 250% variation around its value *)
 (* ---> *)
 (*
alpha_s(mtau^2)=0.320945+{^{-0.003912}_{+0.004746}}_{d4}; OK
10^3 O4= -4.6+{^{+0.5}_{-0.7}}_{d4}; OK
10^3 O61=-3.2+{^{+0.2}_{-0.3}}_{d4}; OK
10^4 O62=+7.7+{^{-0.7}_{+1.0}}_{d4}; OK *)  
 (*----*)
 (* Variation of Nt (truncation order, Nt=Nt0+-2=5+-2 (and 5+-1), i.e., variation Nt=7,6,5,4; tFO global with (O61, O62)) *)
(*
 alO4O612fitPtFOMinnum[7]={0.0216534, {aALf -> 0.116036, O4fitf -> -0.0173627, O61fitf -> -0.00198169, O62fitf -> 0.00111877}}; alpha_s(mtau^2)=0.350259; OK
 alO4O612fitPtFOMinnum[6]={0.00174365, {aALf -> 0.108768, O4fitf -> -0.00790883, O61fitf -> -0.00348822, O62fitf -> 0.000975418}}; alpha_s(mtau^2)=0.329376; OK
alO4O612fitPtFOMinnum[5]={0.00289741, {aALf -> 0.105851, O4fitf -> -0.00457748, 
O61fitf -> -0.00324893, O62fitf -> 0.000772663}}; alpha_s(mtau^2)=0.320945; OK
 alO4O612fitPtFOMinnum[4]={0.000894446, {aALf -> 0.102959, O4fitf -> -0.00162577, O61fitf -> -0.00191124, O62fitf -> 0.000410859}}; alpha_s(mtau^2)=0.312558; OK
 alO4O612fitPtFOMinnum[3]={0.0000550704, {aALf -> 0.103353, O4fitf ->0.000594932, O61fitf -> -0.000289906, O62fitf -> 0.0000639169}}; alpha_s(mtau^2)=0.313702; INCLUDED; OK *)
  -->
  (* Nt = 5+1; 5 +2; Nt=5-1; Nt=5-2 (Note, when d4=275., we do include Nt=5-2 in the variation delta Nt for tFO): *)
alpha_s(mtau^2)=0.320945+{^{+0.029314}_{-0.008382}}_{Nt}; OK
10^3 O4= -4.6+{^{-12.8}_{+5.2}}_{Nt}; OK
10^3 O61=-3.2+{^{+1.2}_{+3.5}}_{Nt}; OK
10^4 O62=+7.7+{^{+3.5}_{-7.1}}_{Nt}; OK *) 
 (*-----------------------------------------*)
 (* Altogether, the tFO global fit approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax)=aALftFO[Nt] to a(Q^2=sigmax Exp[I fi]) ], gives: *)
(* tFO
 Nt=5; ddth[4,tkap=275.; kap=1.+{^{+1}_{-1/3}}; d4=275+-63.19; Nt=5+-2;
 alpha_s(mtau^2)=0.320945+{^{+0.000342}_{-0.000342}}_{exp}+{^{-0.003824}_{+0.020115}}_{kap}+{^{-0.003912}_{+0.004746}}_{d4}+{^{-0.029314}_{-0.008382}}_{Nt} =0.320945+{^{+0.020670}_{-0.029822}} = 0.3209+{^{+0.0207}_{-0.0298}} =app 0.321+{^{+0.021}_{-0.030}}; OK
 10^3 O4= -4.6+{^{+0.0}_{-0.0}}_{exp}+{^{+2.9}_{-6.8}}_{kap}+{^{+0.5}_{-0.7}}_{d4}+{^{-12.8}_{+5.2}}_{Nt} = -4.6+{^{+6.0}_{-14.5}}; OK 
 10^3 O61=-3.2+{^{+0.1}_{-0.1}}_{exp}+{^{+1.3}_{-0.0}}_{kap}+{^{+0.2}_{-0.3}}_{d4}+{^{+1.2}_{+3.5}}_{Nt} = -3.2+{^{+3.7}_{-0.3}}; OK
10^4 O62=+7.7+{^{+0.1}_{-0.1}}_{exp}+{^{-3.8}_{+4.9}}_{kap}+{^{-0.7}_{+1.0}}_{d4}+{^{+3.5}_{-7.1}}_{Nt}  = +7.7+{^{+6.1}_{-8.1}}; OK
OK *)
 (* End of tFO global fit results with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax)=aALftFO[Nt] to a(Q^2=sigmax Exp[I fi]) ] OK *)
(*------------------------------------*)
(*------------------------------------*)
(*------------------------------------*)
 (*------------------------------------------------------------------*)
 (*------------------------------------------------------------------*)
  (* End of  FO and tFO approaches *)
  (*------------------------------------------------------------------*)
 (*------------------------------------------------------------------*)
   (*------------------------------------------------------------------*)
 (*------------------------------------------------------------------*) 
    (* (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) *)

   BPCID0int[M2_?NumberQ, Nn_, aALf_]:=Module[{aa, Nint,tabACI,jj, ACIintup, ACIint,dffint},aa=aALf; Nint=40; tabACI = Table[{jj*Pi/Nint, aCI[jj*Pi/Nint, aa]}, {jj, 0, Nint - 1}]; ACIintup=Interpolation[tabACI]; ACIint[th_]:=If[th > 0, ACIintup[th], Conjugate[ACIintup[-th]]]; dffint[th_]:= Sum[ddtk[nn - 1,tkap] ACIint[th]^nn, {nn, 1, 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)]) dffint[th], {th, -Pi, Pi}]];
   (* Gives the same result of contour integratioon as BPCID), but needs about 50% more time; how come? I wanted to have a faster CI-evaluation of the contour integral *)
   
  JJ[1, AA_] :=E^-AA - AA (-1 + EulerGamma + Gamma[0, AA] + Log[AA]);
JJ[2,AA_]:=1/2 (-1 + AA) E^-AA - 1/4 AA^2 (-3 + 2 EulerGamma - 2 CoshIntegral[AA] + 2 Log[AA] + 2 SinhIntegral[AA]);
 JJ[3,AA_]:=1/6 (2 + (-1 + AA) AA) E^-AA - 1/36 AA^3 (-11 + 6 EulerGamma - 6 CoshIntegral[AA] + 6 Log[AA] + 6 SinhIntegral[AA]);

   ReBPCI[M2f_,Nn_,aALf_,O4_,O61_,O62_]:=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 + (O61 + O62 (1/aALf - be0 Log[kappa]))/sigmax) + (1/M2f^2) (O4 + 2 ((O61 + O62 (1/aALf - be0 Log[kappa]))/sigmax)) + (1/M2f^3) (1/2)(O61 + O62 (1/aALf - be0 Log[kappa]))  ) + (6 Pi^2/sigmax^3) be0 O62 (JJ[1,sigmax/M2f] + JJ[2,sigmax/M2f] 2 (1 - M2f/sigmax) + JJ[3,sigmax/M2f] (1 - 2 M2f/sigmax + 2 (M2f/sigmax)^2) - (1/3) 2 (M2f/sigmax)^2 Exp[-sigmax/M2f])];
   (* This is the CIPT approach for Borel-Laplace Re[B[M2,Nn,O4,O61,O62]] 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; a(\sigmax Exp[I fi]) in O61/a(sigmax Exp[I fi]) runs by one-loop RGE, and has the initial value a(sigmax)=aALf (the same initial value as in D=0 part, but there it runs by 5-loop RGE) *)
   
(* 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.0030; O61init = -0.001; O62init = 0.001; aALfinit = 0.1100;
alO4O612fitPCIMinN[Nn_]:= FindMinimum[Sum[(ReBPCI[M2FP[jj] Exp[I Pi/6], Nn, aALf,O4fitf,O61fitf,O62fitf] - 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,O61fitf,O62fitf] - 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,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{aALf,aALfinit},{O4fitf, O4init},{O61fitf, O61init},{O62fitf, O62init}}];
(* aALf=alpha_s(kappa*sigmax)/Pi *)
  (*-----------------------__*)
(* NEW BL results, CI approach with (O61, O62): *)
(*

 (* Timing[alO4O612fitPCIMinN[10]]={2746.79, {5.33992, {aALf -> 0.093652, O4fitf -> 0.00150002, O61fitf -> 0.00304009, O62fitf -> -0.000729873}}};
Timing[alO4O612fitPCIMinN[11]]={2698.06, {5.83549, {aALf -> 0.0926053, O4fitf -> 0.00198102, O61fitf -> 0.00403113, O62fitf -> -0.000849245}}};
  Timing[alO4O612fitPCIMinN[12]]={3127.19, {10.2577, {aALf -> 0.0846949, O4fitf -> 0.00542309, O61fitf -> 0.0122664, O62fitf -> -0.00172118}}}; *)

 alO4O612fitPCIMinnum[3]={0.000333408, {aALf -> 0.119688, O4fitf -> -0.00307405, O61fitf -> -0.00244564, O62fitf -> 0.000498155}}; (* OK *)
 alO4O612fitPCIMinnum[4]={0.000497056, {aALf -> 0.115532, O4fitf -> -0.00362321, O61fitf -> -0.00315366, O62fitf -> 0.000581503}}; (* OK *)
alO4O612fitPCIMinnum[5]= {0.00151674, {aALf -> 0.114577, O4fitf -> -0.00444899, O61fitf -> -0.00441847, O62fitf -> 0.000755956}}; (* OK *)
 (* when d4=338.19: alO4O612fitPCIMinnum[5]={0.00175176, {aALf -> 0.114323, O4fitf -> -0.00460931, O61fitf -> -0.00466055, O62fitf -> 0.000788052}}; *)
 alO4O612fitPCIMinnum[6]={0.0280379, {aALf -> 0.116287, O4fitf -> -0.00693134, O61fitf -> -0.00892589, O62fitf -> 0.00141576}}; (* OK *)
alO4O612fitPCIMinnum[7]={0.230799, {aALf -> 0.113335, O4fitf -> -0.00609996, O61fitf -> -0.00851868, O62fitf -> 0.0012252}}; (* OK *)
alO4O612fitPCIMinnum[8]={1.3501, {aALf -> 0.106106, O4fitf -> -0.00320983, O61fitf -> -0.00506373, O62fitf -> 0.000504002}}; (* OK *)
alO4O612fitPCIMinnum[9]={1.54723, {aALf -> 0.105263, O4fitf -> -0.00283077, O61fitf -> -0.00452995, O62fitf -> 0.000411221}}; (* OK *)
alO4O612fitPCIMinnum[10]={4.7654, {aALf -> 0.0956961, O4fitf -> 0.0013899, O61fitf -> 0.00245993, O62fitf -> -0.000636686}}; (* OK *)
alO4O612fitPCIMinnum[11]={0.00726296, {aALf -> 0.108183, O4fitf -> -0.00437088, O61fitf -> -0.00475524, O62fitf -> 0.000460737}}; (* OK *)
 alO4O612fitPCIMinnum[12]={10.0229, {aALf -> 0.0856433, O4fitf -> 0.00564068, O61fitf -> 0.0122339, O62fitf -> -0.00171901}}; (* OK *)

 Clear[Nn]; Do[Module[{}, xi2CI[Nn]=alO4O612fitPCIMinnum[Nn][[1]]; aALfCI[Nn]=aALf/.alO4O612fitPCIMinnum[Nn][[2]]; alsmtauCI[Nn]=aMSext[mtauphys^2/(kappa sigmax),aALfCI[Nn]] Pi;  O4CI[Nn]=O4fitf/.alO4O612fitPCIMinnum[Nn][[2]]; O61CI[Nn]=O61fitf/.alO4O612fitPCIMinnum[Nn][[2]]; O62CI[Nn]=O62fitf/.alO4O612fitPCIMinnum[Nn][[2]]],{Nn,3,12}];

*)
(* end of NEW BL results, CI approach with (O61,O62)) *)
   (*------------------------------------*)
 (*---------------------------------------------------------------------_*)
 (* 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 experimental ALEPH values for (0,0) moment a^{(0,0)} (see also the file rtauALEPH.m): *)
grs0a[0,0][sig_]:=(1/sigmax);
aexpBin77[0,0]=(Sum[ dsbin[[ii]] grs0a[0,0][sbin[[ii]]] OmegaVplA[ii],{ii,1,77}] + 2 Pi^2 fPi^2 grs0a[0,0][mPi^2]) -  1;
delaexpBin77[0,0]=(2 Pi^2) Sqrt[Sum[ grs0a[0,0][sbin[[ii]]] grs0a[0,0][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^{(0,0)}_{exp} = aexpBin77[0,0] +- delaexpBin77[0,0] *)
(*--------*)
(* 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); 
(*---------*)
(* moments a^{(0,0}}, CI: *)				     
aaCID0fast[0,0,aALf_][Nn_]:=(1/(2 Pi)) NIntegrate[ (1 + Exp[I th])  dff[th,Nn,aALf], {th, -Pi, Pi}];
(*-------------------------*)
(*  This is the NEW (7 Sept.2021) (O61, O62) approach: *)
(* This is the case when a(sigmax Exp[I fi]) in [O_6^{(2)}/a(sigmax Exp[I fi])] runs by one-loop RGE: *)
aaCI21[aALf_][NN_,O61_,O62_]:=aaCID0fast[2,1,aALf][NN]+ (6 Pi^2/sigmax^3)*( (-2)(O61+O62/aALf) + O62 be0 (1/3) );  (* this is a^{(2,1)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE  *)
aaCI20[aALf_][NN_,O4_,O61_,O62_]:=aaCID0fast[2,0,aALf][NN]+ (6 Pi^2) ( 2 O4/sigmax^2 + (1/sigmax^3)( (O61+O62/aALf) + O62 be0 (11/6) ) );   (* this is a^{(2,0)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE  *)
(*----------*)
  (*-----------------------------_*)
  (* Using the results at Nt=5: 
   alO4O612fitPCIMinnum[5]={0.00151674, {aALf -> 0.114577, O4fitf -> -0.00444899,O61fitf -> -0.00441847, O62fitf -> 0.000755956}}; (* OK *)
     (* ---> *)
   aaCI21[aALfCI[5]][5, O61CI[5], O62CI[5]]=0.207000;  (* OK *)
   aaCI20[aALfCI[5]][5, O4CI[5], O61CI[5], O62CI[5]]=0.240764;  (* OK *)

 Using the results at Nt=4: 
    alO4O612fitPCIMinnum[4]={0.000497056, {aALf -> 0.115532, O4fitf -> -0.00362321, O61fitf -> -0.00315366, O62fitf -> 0.000581503}}; (* OK *)
     (* --> *)
   aaCI21[aALfCI[4]][4, O61CI[4], O62CI[4]] =0.206608; (* OK *)
   aaCI20[aALfCI[4]][4, O4CI[4], O61CI[4], O62CI[4]]=0.240598; (* OK *)
      
    VERY GOOD, because:
   aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
   aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026
   *)
  (*------------------------*)
(*  This is the NEW (7 Sept.2021) (O61, O62) approach: *)
   (* Tables of a20, a21 vs Nt in CI when aALf=a(sigmax) is VARIED ('avar'), i.e.,using alO4O612fitPCIMinnum[N] (N=3,4,...12) from global CI (with a(sigmax Exop[I fi]) in O62/a(sigmax Exp[I fi]) running by one-loop RGE): *)			     
(*
 Tab20D0CIavar=Table[{NN,Re[aaCID0fast[2,0,aALfCI[NN]][NN]]},{NN,3,12}]; (* Tab20D0CIavar={{3, 0.276829}, {4, 0.283791}, {5, 0.293682}, {6, 0.321435}, {7, 
  0.312362}, {8, 0.280885}, {9, 0.276687}, {10, 0.230251}, {11, 0.3018}, {12, 0.183764}} (* OK *) *)
 Tab20CIavar=Table[{NN,Re[aaCI20[aALfCI[NN]][NN,O4CI[NN],O61CI[NN],O62CI[NN]]]},{NN,3,12}]; (* Tab20CIavar={{3, 0.240564}, {4, 0.240598}, {5, 0.240764}, {6, 0.241245}, {7, 0.240029}, {8, 0.237158}, {9, 0.236818}, {10, 0.232851}, {11, 0.239559}, {12, 0.228704}}}; (* OK *) *)
 delTab20CIavar= Table[{Nn, Tab20CIavar[[Nn - 2]][[2]] - Tab20CIavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];
 Tab21D0CIavar=Table[{NN,Re[aaCID0fast[2,1,aALfCI[NN]][NN]]},{NN,3,12}]; (*  Tab21D0CIavar={{3, 0.214777}, {4, 0.215572}, {5, 0.217228}, {6, 0.222925}, {7, 0.215609}, {8, 0.196912}, {9, 0.194695}, {10, 0.1691}, {11, 0.200591}, {12, 0.142778}}; (* OK *) *)
 Tab21CIavar=Table[{NN,Re[aaCI21[aALfCI[NN]][NN,O61CI[NN],O62CI[NN]]]},{NN,3,12}]; (* Tab21CIavar={{3, 0.206524}, {4, 0.206608}, {5, 0.207}, {6, 0.208261}, {7, 
  0.205724}, {8, 0.199625}, {9, 0.19889}, {10, 0.190435}, {11, 0.204201}, {12, 0.181586}}; (* OK *) *)
delTab21CIavar= Table[{Nn, Tab21CIavar[[Nn - 2]][[2]] - Tab21CIavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}]; 
(* OK *) *)

 (*-------*)   
(*
(*  This is the NEW (30 Aug.2021) (O61, O62) approach: *)
 (* Figures for global CI: a21, a20, a10 vs Nt [aALf=a(sigmax) is VARIED ('avar')], with the NEW (30 Aug.2021) (O61, O62) approach : *)
a20CIavar= ListPlot[{Tab20D0CIavar, Tab20CIavar, ff20exp, ff20miexp, ff20plexp},PlotRange -> {{2.8, 12.2}, {0.18, 0.33}}, PlotStyle -> {Blue, Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.9, 0.95}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.08, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.70, 0.95}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"CI (D=0)", "CI"}, {0.2, 0.15}]];
a20CIavarZoom= ListPlot[{Tab20CIavar, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.228, 0.244}}, PlotStyle -> {Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.40, 0.40}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.10, 0.45}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.15, 0.10}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"CI"}, {0.40, 0.10}]];
a21CIavar= ListPlot[{Tab21D0CIavar, Tab21CIavar, ff21exp, ff21miexp, ff21plexp},PlotRange -> {{2.8, 12.2}, {0.14, 0.225}}, PlotStyle -> {Blue, Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.50, 0.40}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.08, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.15, 0.4}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"CI (D=0)", "CI"}, {0.17, 0.15}]];
a21CIavarZoom= ListPlot[{Tab21CIavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.18, 0.210}}, PlotStyle -> {Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.50, 0.40}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.08, 0.55}]],Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 11], Scaled[{0.15, 0.30}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"CI"}, {0.15, 0.10}]];
(* OK *)
 *)
    (* Export["a20CIavarO612d4275.pdf",a20CIavar]; Export["a21CIavarO612d4275.pdf",a21CIavar];  Export["a20CIavarO612Zd4275.pdf",a20CIavarZoom]; Export["a21CIavarO612Zd4275.pdf",a21CIavarZoom];  (* OK *) *)
    (* global CI: a21 and a20 suggest that Nt=4 is the most reasonable truncation (in this global CI fit approach), because the difference from Nt toNt+1 has the smallest absolute value there; *)
    (*------)*

  (* end of gobal CI approach *)  
  (*-------------------*)
    (* Tables of a20(D=0), a21(D=0) vs Nt in CI when aALf=a(sigmax) is FIXED ('afxd'), i.e., aAL=aALfCI[4](=0.115532), [ alpha_s(mtau^2)=alsmtauCI[4]=0.348817],i.e., the coupling alpha_s taken FIXED at value as obtained from global BL CI fits at Nt=4: *)
(*
Tab20D0CIafxd = 
 Table[{NN, Re[aaCID0fast[2, 0, aALfCI[4]][NN]]}, {NN, 3, 12}]; (* Tab20D0CIafxd={{3, 0.260487}, {4, 0.283791}, {5, 0.298644}, {6, 0.316883}, {7, 0.326313}, {8, 0.345514}, {9, 0.347478}, {10, 0.379583}, {11, 0.356853}, {12, 0.447219}}; *) (* OK *)
 Tab21D0CIafxd=Table[{NN,Re[aaCID0fast[2,1,aALfCI[4]][NN]]},{NN,3,12}]; (* Tab21D0CIafxd={{3, 0.205763}, {4, 0.215572}, {5, 0.21951}, {6, 0.221229}, {7, 0.220326}, {8, 0.214399}, {9, 0.213319}, {10, 0.186833}, {11, 0.212778}, {12, 0.0764723}}; *) (* OK *)
ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}];
 *)
(*
 (* Figures for global CI: a21(D=0), a20(D=0), vs Nt, WHERE aAL=a(sigmax) is FIXED ('aFXD'); a(sigmax)=aALfCI[5](=0.104501), [ alpha_s(mtau^2)=alsmtauCI[5]= ]: *)
a20CIafxd= ListPlot[{Tab20D0CIafxd, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.20, 0.40}}, PlotStyle -> {Blue, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.70, 0.90}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]] == Round[10000 alsmtauCI[4]]/10000., 13], Scaled[{0.75, 0.45}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.75, 0.35}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"CI (D=0)"}, {0.77, 0.58}]];
a21CIafxd= ListPlot[{Tab21D0CIafxd, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.18, 0.23}}, PlotStyle -> {Blue, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.50, 0.90}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]] == Round[10000 alsmtauCI[4]]/10000., 13], Scaled[{0.22, 0.22}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.15, 0.10}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"CI (D=0)"}, {0.17, 0.35}]];
(* OK *)
    (*------*)
 (*---------------------------------------*)
(*-------------------------------------------------*)									
    (* Now continuing the global fit CI approach (with (O61,O62)): Nt=4 is the preferred order of truncation, according to our analysis of a^{(2,1)}, a^{(2,0). *)

Nt0=4;									
(* The precise central value result of the CI global fit: *)
(* alO4O612fitPCIMinnum[4]={0.000497056, {aALf -> 0.115532, O4fitf -> -0.00362321, O61fitf -> -0.00315366, O62fitf -> 0.000581503}}; *)
									
(* First the estimate fo EXPERIMENTAL uncertainties of alAL, O4, O61, O62 in the Nt=4 CI global Borel-Laplace fit; due to instabilities and possible underestimations, we do not use the correlated (with variances) chi^2, but rather the uncorrelated chi^2:  *)									
(* 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 *)

chi2o2O4[O4f_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALfCI[Nt0],O4f,O61CI[Nt0],O62CI[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[4]=-0.00362321, by global analysis *)

chi2o2O61[O61f_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALfCI[Nt0],O4CI[Nt0],O61f,O62CI[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O61CIvarfitMin=FindMinimum[chi2o2O61[O61fitf],{O61fitf,O61CI[Nt0]}];
O61CIvarfit=O61fitf/.O61CIvarfitMin[[2]];
Xi2O61CIvaro77=O61CIvarfitMin[[1]];
 delO61CI= delO61f/.FindRoot[chi2o2O61[O61CIvarfit+delO61f]==Xi2O61CIvaro77+1,{delO61f, 3. 10^-5}];
(* Nm=2; Nd=6; delO61CI=0.0000451; *)

chi2o2O62[O62f_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALfCI[Nt0],O4CI[Nt0],O61CI[Nt0],O62f] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O62CIvarfitMin=FindMinimum[chi2o2O62[O62fitf],{O62fitf,O62CI[Nt0]}];
O62CIvarfit=O62fitf/.O62CIvarfitMin[[2]];
Xi2O62CIvaro77=O62CIvarfitMin[[1]];
 delO62CI= delO62f/.FindRoot[chi2o2O62[O62CIvarfit+delO62f]==Xi2O62CIvaro77+1,{delO62f, 3. 10^-5}];
(* Nm=2; Nd=6; delO62CI=5.38325*10^-6; *)

chi2o2aAL[aALf_]:=Sum[(ReBPCI[M2covF[jj],Nt0,aALf,O4CI[Nt0],O61CI[Nt0],O62CI[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.000229774;  chi2o2aAL[aALCIvarfit +-0.000168] =app 1.00 *)
(* Nm=2; Nd=6; Nt0=5; delaALCI=0.000168 *)

(* The precise central value result of the CI global fit: 
 alO4O612fitPCIMinnum[4]={0.000497056, {aALf -> 0.115532, O4fitf -> -0.00362321, O61fitf -> -0.00315366, O62fitf -> 0.000581503}}; *)
(*----*)
(* Together, the estimates of EXPERIMENTAL uncertainties in the CI global fit approach (with Nt=4) are (Nm=2; Nd=6; uncorrelated xi^2):  *)
(* delO4CI= 0.0000303853; delO61CI=0.0000450575; delO62CI=5.38325*10^-; delaALCI=0.000168 (=delta(a(sigmax))) --> *)
(* 
alpha_s(mtau^2)=0.348817+{^{+0.000481}_{-0.000481}}_{exp}; 
10^3 O4= -3.6+{^{+0.0}_{-0.0}}_{exp};
10^3 O61=-3.2+{^{+0.0}_{-0.1}}_{exp};
10^4 O62=+5.81+{^{+0.05}_{-0.05}}_{exp} =app +5.8+{^{+0.1}_{-0.1}}_{exp};
OK
 *)

(*---------*)

(* Variation of kappa, CI global fit, Nt=4  (note that aALf=a(kappa sigmax) : *)
(*
kappa=2.0: 
 alO4O612fitPCIMinnum[4]={1.73859*10^-6, {aALf -> 0.0957803, O4fitf -> -0.00221221, O61fitf -> -0.0013134, O62fitf -> 0.000306937}}; alpha_s(mtau^2)=0.356569;
kappa=1.0:
 alO4O612fitPCIMinnum[4]={0.000497056, {aALf -> 0.115532, O4fitf -> -0.00362321, O61fitf -> -0.00315366, O62fitf -> 0.000581503}}; alpha_s(mtau^2)=0.348817;
 kappa=2/3:
 alO4O612fitPCIMinnum[4]={0.168056, {aALf -> 0.134815, O4fitf -> -0.00544054, 
  O61fitf -> -0.00689505, O62fitf -> 0.00114122}}; alpha_s(mtau^2)=0.349204; (* OK *)
 kappa=0.5:
alO4O612fitPCIMinnum[4]={1.28048, {aALf -> 0.137865, O4fitf -> -0.0023611, O61fitf -> -0.00320336, O62fitf -> 0.00035881}}; alpha_s(mtau^2)=0.324807; 
 (* I will use kap_min=2/3, NOT kap_min=0.5, OK *)
 
 *)
alpha_s(mtau^2)=0.348817+{^{+0.007752}_{+0.000387}}_{kap}; 
10^3 O4= -3.6+{^{+1.4}_{-1.8}}_{kap};
10^3 O61=-3.2+{^{+1.9}_{-3.7}}_{kap};
10^4 O62=+5.8+{^{-2.7}_{+5.6}}_{kap};
(* I use kap_min=2/3, NOT kap_min=0.5, OK *)
OK  *)

(* Variation of d4, CI global fit (with (O61,O62)): *)
(* Since Nt=4 in CI global, variation of d4 does NOT affect the results; OK *)

(* Variation of Nt(truncation order, Nt=4+-2), CI global: Nt=6,5,4,3 (we do not go down to Nt=2, an extreme) *)
(* 
alO4O612fitPCIMinnum[6]={0.0280379, {aALf -> 0.116287, O4fitf -> -0.00693134, O61fitf -> -0.00892589, O62fitf -> 0.00141576}}; alsmtauCI[6]=0.350977; (* OK *) 
alO4O612fitPCIMinnum[5]={0.00151674, {aALf -> 0.114577, O4fitf -> -0.00444899, O61fitf -> -0.00441847, O62fitf -> 0.000755956}}; alsmtauCI[5]=0.346082;(* OK *) 
alO4O612fitPCIMinnum[4]={0.000497056, {aALf -> 0.115532, O4fitf -> -0.00362321, O61fitf -> -0.00315366, O62fitf -> 0.000581503}}; alsmtauCI[4]=0.348817;(* OK *)alO4O612fitPCIMinnum[3]={0.000333408, {aALf -> 0.119688, O4fitf -> -0.00307405, O61fitf -> -0.00244564, O62fitf -> 0.000498155}}; alsmtauCI[3]=0.360682;(* OK *) 
(* --> *)
(*
 alpha_s(mtau^2)=0.348817+{^{-0.00273519}_{+0.011865}}_{Nt}; 
10^3 O4= -3.6+{^{-3.3}_{+0.5}}_{Nt};
10^3 O61=-3.2+{^{-5.7}_{+0.8}}_{Nt};
10^4 O62=+5.8+{^{+8.4}_{-0.8}}_{Nt};
OK
 *)
(*----------------*)
(* Altogether, the CI global fit approach with (O61,O62) [ where in O62/a(Q^2) the coupling a(Q^2) runs according to the one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]), with a(sigmax)=aALfCI[Nt] ] gives: *)
(* CI
 Nt=4; d4=275.; kap=1.^{+1}_{-1/3}; Nt=4^{+2}_{-1};
alpha_s(mtau^2)=0.348817+{^{+0.000481}_{-0.000481}}_{exp}+{^{+0.007752}_{+0.000387}}_{kap}+{^{+0.0}_{-0.0}}_{d4}+{^{-0.00273519}_{+0.011865}}_{Nt} = 0.348817+{^{+0.014181}_{-0.002777}}  =app 0.3488+{^{+0.0142}_{-0.0028}} =app 0.349+{^{+0.014}_{-0.003}}; (* OK *)
 10^3 O4= -3.6+{^{+0.0}_{-0.0}}_{exp}+{^{+1.4}_{-1.8}}_{kap}+{^{+0.0}_{-0.0}}_{d4}+{^{-3.3}_{+0.5}}_{Nt}  = -3.6+{^{+1.5}_{-3.8}}; (* OK *)
10^3 O61=-3.2+{^{+0.0}_{-0.1}}_{exp}+{^{+1.9}_{-3.7}}_{kap}_{kap}+{^{+0.0}_{-0.0}}_{d4}+{^{-5.7}_{+0.8}}_{Nt} = -3.2+{^{+2.1}_{-6.8}}; (* OK *)
 10^4 O62=+5.8+{^{+0.1}_{-0.1}}_{exp}+{^{-2.7}_{+5.6}}_{kap}+{^{+0.0}_{-0.0}}_{d4}+{^{+8.4}_{-0.8}}_{Nt} = +5.8+{^{+10.1}_{-2.8}};  (* OK *)
 (* OK *)
 
(* End of CI global fit results approach with (O61,O62) [ where in O62/a(Q^2) the coupling a(Q^2) runs according to the one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi])*)
(*-------------------------------------------------*)
(*-------------------------------------------------*)

 (* End of CI 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.10;  *)
  aAL=0.1079;
(* aCI[th,aAL] = a(kappa sigmax Exp[I th]), and aAL=a(kappa sigmax) *)
  
  delD[th_,Nn_]:= Sum[delddtk[jj-1,tkap] aCI[th,aAL]^jj, {jj, 1, Nn}];
  (* delD[th,N] is the polynomial part of the Borel od Adler D, of order up to a(kappa sigmax Exp[I th])^N; note that kappa=tkap Exp[-tKK] (kappa=mu^2/Q^2) *)
  delDext[th_,Nn_,aALf_]:= Sum[delddtk[jj-1,tkap] aCI[th,aALf]^jj, {jj, 1, Nn}];

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

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

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

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

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

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

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

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

 
 (* This is NEW (7 Sept. 2021) (O61, O62 approach): *)    
ReBPPV[M2f_,Nn_,O4_,O61_,O62_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPPVD0[M2f,Nn]+ 2 Pi^2 ( (1/(M2f sigmax))(2 O4 + (O61 + O62 (1/aAL - be0 Log[kappa]))/sigmax) + (1/M2f^2) (O4 + 2 ((O61 + O62 (1/aAL - be0 Log[kappa]))/sigmax)) + (1/M2f^3) (1/2)(O61 + O62 (1/aAL - be0 Log[kappa]))  ) + (6 Pi^2/sigmax^3) be0 O62 (JJ[1,sigmax/M2f] + JJ[2,sigmax/M2f] 2 (1 - M2f/sigmax) + JJ[3,sigmax/M2f] (1 - 2 M2f/sigmax + 2 (M2f/sigmax)^2) - (1/3) 2 (M2f/sigmax)^2 Exp[-sigmax/M2f])]; (* Re[B[M2,Nn,O4,O61,O62]] for PV with delD truncated at a^Nn, RScl kappa; a(\sigmax Exp[I fi]) in O61/a(sigmax Exp[I fi]) runs by one-loop RGE, and has the initial value a(sigmax)=aAL (the same initial value as in D=0 part, but there it runs by 5-loop RGE) *) 
    (* aALf= a(kappa sigmax) *)

  (* Note that, instead of simply O62/aAL, we have O62 (1/aAL - be0 Log[kappa]); this is because aAL=a(kappa sigmax), but at O62 we have to have 1/a(sigmax); we use the one-loop RGE running from Q^2=kappa sigmax) to Q^2=sigmax --> 1/a(sigmax)=(1/aAL - be0 Log[kappa]) *)
 
 (* 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): *)

(* This is NEW (7 Sept. 2021) ((O61, O62) approach): *)
  O4init = -0.0030; O61init = -0.003; O62init = 0.001; (* aALfinit = 0.1100; *)
O4O612fitPPVMinN[Nn_]:= FindMinimum[Sum[(ReBPPV[M2FP[jj] Exp[I Pi/6], Nn, O4fitf,O61fitf,O62fitf] - 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,O61fitf,O62fitf] - 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,O61fitf,O62fitf] - ReBPexpBinAL[M2FP[jj]])^2/delReBPexpBinAL[M2FP[jj]]^2, {jj, 0, NchP}], {{O4fitf, O4init},{O61fitf, O61init},{O62fitf, O62init}}];
(* aALf=alpha_s(kappa*sigmax)/Pi *)

(* Do[O4O612fitPPVMinnum[nn]=O4O612fitPPVMinN[nn], {nn, 4, 10}]; *)

(* NEW BL results, PPV approach with (O61,O62): *)
  (* Summary, for kappa=1 and NchP=10, PV(+CI) approach: *)

(*
 aALfPV[3]=0.093; O4O612fitPPVMinnum[3]={0.00701044, {O4fitf -> -0.00463517, O61fitf -> -0.00409335, O62fitf -> 0.000571855}}; (* OK *)
 aALfPV[4]=0.1046; O4O612fitPPVMinnum[4]={5.96207*10^-6, {O4fitf -> -0.00684832, O61fitf -> -0.00931311, O62fitf -> 0.0012368}}; (* OK *)
 aALfPV[5]=0.1082; O4O612fitPPVMinnum[5]={0.000464085, {O4fitf -> -0.00616911, O61fitf -> -0.00968949, O62fitf -> 0.00127478}}; (* OK *)
 aALfPV[6]=0.1083; O4O612fitPPVMinnum[6]={0.00195454, {O4fitf -> -0.00516294, O61fitf -> -0.00800203, O62fitf -> 0.00105598}}; (* OK *)
 aALfPV[7]=0.1077; O4O612fitPPVMinnum[7]={0.00356788, {O4fitf -> -0.00437873, O61fitf -> -0.00631874, O62fitf -> 0.00083764}}; (* OK *)
 aALfPV[8]=0.1079; O4O612fitPPVMinnum[8]={0.00310712, {O4fitf -> -0.00455454, O61fitf -> -0.00673939, O62fitf -> 0.000892787}}; (* OK *)
 aALfPV[9]=0.1080; O4O612fitPPVMinnum[9]={0.00287242, {O4fitf -> -0.00463592, O61fitf -> -0.00695691, O62fitf -> 0.000921253}}; (* OK *)
 aALfPV[10]=0.1111; O4O612fitPPVMinnum[10]={0.390203, {O4fitf -> -0.00741201, O61fitf -> -0.0157418, O62fitf -> 0.00208312}}; (* OK *)
 aALfPV[11]=0.1085; O4O612fitPPVMinnum[11]={0.536637, {O4fitf -> -0.00567423, O61fitf -> -0.0117142, O62fitf -> 0.00148205}}; (* OK *)
 aALfPV[12]=0.095; O4O612fitPPVMinnum[12]={4.71276, {O4fitf -> 0.00162842, O61fitf -> 0.0028044, O62fitf -> -0.000607198}}; (* OK *)

 Clear[Nn]; Do[Module[{}, xi2PV[Nn]=O4O612fitPPVMinnum[Nn][[1]]; alsmtauPV[Nn]=aMSext[mtauphys^2/(kappa sigmax),aALfPV[Nn]] Pi;  O4PV[Nn]=O4fitf/.O4O612fitPPVMinnum[Nn][[2]]; O61PV[Nn]=O61fitf/.O4O612fitPPVMinnum[Nn][[2]]; O62PV[Nn]=O62fitf/.O4O612fitPPVMinnum[Nn][[2]]],{Nn,3,12}];
*)

(* end of NEW BL results (O61,O62) approach in PPV *)
  (*---------------------------------------------------------------------_*)
  (*---------------------------------------------------------------------_*)
 (* 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 experimental ALEPH values for (0,0) moment a^{(0,0)} (see also the file rtauALEPH.m): *)
grs0a[0,0][sig_]:=(1/sigmax);
aexpBin77[0,0]=(Sum[ dsbin[[ii]] grs0a[0,0][sbin[[ii]]] OmegaVplA[ii],{ii,1,77}] + 2 Pi^2 fPi^2 grs0a[0,0][mPi^2]) -  1;
delaexpBin77[0,0]=(2 Pi^2) Sqrt[Sum[ grs0a[0,0][sbin[[ii]]] grs0a[0,0][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^{(0,0)}_{exp} = aexpBin77[0,0] +- delaexpBin77[0,0] *)
  (*-----------------------*)
 (* 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) *)
   (* 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 *)

   (* (* For FIXED aAL=aALfPV[6]:*) Clear[Nt]; Clear[aaPVD0fast]; Clear[aaPVfast]; Do[Module[{Nti}, Nti = Ntt; aaPVD0fast[2, nn_][Nti] := (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, Nti], {th, -Pi, Pi}];aaPVfast[2, n_][Nti, O2np4_, O2np6_] := aaPVD0fast[2, n][Nti] + ((n + 3)/(n + 1)) 2 Pi^2 (-1)^n ((n + 2) O2np4/sigmax^(n + 2) + (n + 1) O2np6/sigmax^(n + 3))], {Ntt, 3, 15}]; *)
   (*-----*)
(* 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] *)
   (* 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] *)
   (*---------*)
   (* moments a^{(0,0}}, PV (at Nt truncation order): *)
 aaPVD0fast[0,0][Nt_]:=(1/(2 Pi)) NIntegrate[ (1 + Exp[I th])   DBorfast[th,Nt], {th, -Pi, Pi}];
(*-------------------------*)
(*  This is the NEW (9 Sept.2021) (O61, O62) approach: *)
(* This is the case when a(sigmax Exp[I fi]) in [O_6^{(2)}/a(sigmax Exp[I fi])] runs by one-loop RGE: *)
aaPV21[Nt_,O61_,O62_]:=aaPVD0fast[2,1][Nt]+ (6 Pi^2/sigmax^3)*( (-2)(O61+O62/aAL) + O62 be0 (1/3) );  (* this is a^{(2,1)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE  *)
aaPV20[Nt_,O4_,O61_,O62_]:=aaPVD0fast[2,0][Nt]+ (6 Pi^2) ( 2 O4/sigmax^2 + (1/sigmax^3)( (O61+O62/aAL) + O62 be0 (11/6) ) );   (* this is a^{(2,0)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O62/a(sigmax Exp[I fi]) runs by one-loop RGE  *)

 (* Note that aaPV21 and aaPV20 are valid only when kappa=1, i.e., when aAL=a(sigmax) *)
 (*----------*)
  (*-----------------------------*)
  (* Using the results at Nt=6: *)
 (*  aALfPV[6]=0.1083; O4O612fitPPVMinnum[6]={0.00195454, {O4fitf -> -0.00516294, O61fitf -> -0.00800203, O62fitf -> 0.00105598}}; OK 
        (* ---> *)
   (* aaPV21[6, O61PV[6], O62PV[6]]=0.20507; OK
    aaPV20[6, O4PV[6], O61PV[6], O62PV[6]]=0.239949; OK *)


    GOOD, because:
   aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
   aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026
   *)
  (*------------------------*)
(*  This is the NEW (9 Sept.2021) (O61, O62) approach: *)
   (* Tables of a20, a21 vs Nt in CI when aALf=a(sigmax) is VARIED ('avar'), i.e.,using alO4O612fitPCIMinnum[N] (N=3,4,...12) from global CI (with a(sigmax Exop[I fi]) in O62/a(sigmax Exp[I fi]) running by one-loop RGE): *)			     
(*
 (* Note, have to adjust aAL to aAL=aALfPV[NN] at each NN and calculate at that NN first the DBorfast etc. *)

 Tab20D0PVavar(*=Table[{NN,Re[aaPVD0fast[2,0][NN]]},{NN,3,12}]*)={{3,0.298009},{4,0.323306},{5,0.313525},{6,0.301476},{7,0.292626},{8,0.294517},{9,0.295342},{10,0.320053},{11,0.302579},{12,0.224631}}; (* OK *)
 Tab20PVavar(*=Table[{NN,Re[aaPV20[NN,O4PV[NN],O61PV[NN],O62PV[NN]]]},{NN,3,12}]*)={{3,0.239896},{4,0.240388},{5,0.240161},{6,0.239949},{7,0.239734},{8,0.239789},{9,0.239805},{10,0.239378},{11,0.2386},{12,0.232798}};  (* OK *)
 delTab20PVavar= Table[{Nn, Tab20PVavar[[Nn - 2]][[2]] - Tab20PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];

 Tab21D0PVavar(*=Table[{NN,Re[aaPVD0fast[2,1][NN]]},{NN,3,12}]; *) ={{3,0.214934},{4,0.217151},{5,0.214291},{6,0.212367},{7,0.210746},{8,0.211171},{9,0.211363},{10,0.216357},{11,0.210102},{12,0.172166}};  (* OK *)
 Tab21PVavar(*=Table[{NN,Re[aaPV21[NN,O61PV[NN],O62PV[NN]]]},{NN,3,12}];*) ={{3,0.205},{4,0.206106},{5,0.205582},{6,0.20507},{7,0.204571},{8,0.204697},{9,0.204739},{10,0.204342},{11,0.202605},{12,0.190291}};  (* OK *)
delTab21PVavar= Table[{Nn, Tab21PVavar[[Nn - 2]][[2]] - Tab21PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];

 ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}]; 
 (* OK *)
 *)
 (*-------*)   
(*
(*  This is the NEW (9 Sept.2021) (O61, O62) approach: *)
 (* Figures for global PV: a21, a20, a10 vs Nt [aALf=a(sigmax) is VARIED ('avar')], with the NEW (30 Aug.2021) (O61, O62) approach : *)
 a20PVavar= ListPlot[{Tab20D0PVavar, Tab20PVavar, ff20exp, ff20miexp, ff20plexp},PlotRange -> {{2.8, 12.2}, {0.22, 0.35}}, PlotStyle -> {Blue, Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.80, 0.85}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.10, 0.95}]],Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], 
 Scaled[{0.50, 0.85}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"PV (D=0)", "PV"}, {0.17, 0.4}]]; (* OK *)
 a20PVavarZoom= ListPlot[{Tab20PVavar, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.237, 0.244}}, PlotStyle -> {Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.60, 0.95}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.10, 0.95}]],Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.40, 0.95}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"PV"}, {0.90, 0.95}]]; (* OK *)

a21PVavar= ListPlot[{Tab21D0PVavar, Tab21PVavar, ff21exp, ff21miexp, ff21plexp},PlotRange -> {{2.8, 12.2}, {0.20, 0.22}}, PlotStyle -> {Blue, Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {4 -> {5}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.70, 0.85}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.08, 0.95}]],Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], 
 Scaled[{0.90, 0.85}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"PV (D=0)", "PV"}, {0.45, 0.80}]];  (* OK *)
a21PVavarZoom= ListPlot[{Tab21PVavar, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.202, 0.210}}, PlotStyle -> {Red, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True, True}, PlotMarkers -> {\[FilledSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.60, 0.95}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.10, 0.95}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.35, 0.95}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"PV"}, {0.90, 0.95}]]; (* OK *)

 *)

 (* At Nt=8 (aAL=0.1079) local stability (also Nt=9 could be; aAL=0.1080) see a20PVavarZoom and  a21PVavarZoom; so choose Nt=8 *)
 (* Export["a20PVavarO612d4275.pdf",a20PVavar]; Export["a21PVavarO612d4275.pdf",a21PVavar];  Export["a20PVavarO612Zd4275.pdf",a20PVavarZoom]; Export["a21PVavarO612Zd4275.pdf",a21PVavarZoom];  *)
 (*-------------------*)
    (* Tables of a20(D=0), a21(D=0) vs Nt in PV when aALf=a(sigmax) is FIXED ('afxd'), i.e., aAL=aALfPV[8](=0.1079), [ alpha_s(mtau^2)=alsmtauPV[8]=0.32687],i.e., the coupling alpha_s taken FIXED at value as obtained from global BL PV fits at Nt=8: *)
 (*
  (* Note, now we do not adjust for each NN to aAL=aALfPV[NN]; now aAL=0.1079 (=aALfPV[6]) and we use throughout DBorfast for that aAL, at each NN: *)
  Tab20D0PVafxd(*=Table[{NN,Re[aaPVD0fast[2,0][NN]]},{NN,3,12}]*)={{3, 0.462316}, {4, 0.3517}, {5, 0.311372}, {6, 0.298846}, {7, 0.293882}, {8, 0.294517}, {9, 0.294708}, {10, 0.298383}, {11, 0.298549}, {12, 0.305392}}; (* OK *)
Tab21D0PVafxd(*=Table[{NN, Re[aaPVD0fast[2, 1][NN]]}, {NN, 3, 12}]*)={{3, 0.279237}, {4, 0.226856}, {5, 0.213504}, {6, 0.211307}, {7, 0.211291}, {8, 0.211171}, {9, 0.211094}, {10, 0.208703}, {11, 0.208549}, {12, 0.199968}}; (* OK *)
 ff20exp=Table[{NN,aexpBin77[2,0]},{NN,3,12}]; ff20plexp=Table[{NN,aexpBin77[2,0]+delaexpBin77[2,0]},{NN,3,12}]; ff20miexp=Table[{NN,aexpBin77[2,0]-delaexpBin77[2,0]},{NN,3,12}];
ff21exp=Table[{NN,aexpBin77[2,1]},{NN,3,12}];  ff21plexp=Table[{NN,aexpBin77[2,1]+delaexpBin77[2,1]},{NN,3,12}]; ff21miexp=Table[{NN,aexpBin77[2,1]-delaexpBin77[2,1]},{NN,3,12}]; (* OK *)
 *)
(*
 (* Figures for global PV: a21(D=0), a20(D=0), vs Nt, WHERE aAL=a(sigmax) is FIXED ('aFXD'); a(sigmax)=aALfPV[8](=0.1079), [ alpha_s(mtau^2)=alsmtauPV[8]==0.32687]: *)
a20PVafxd= ListPlot[{Tab20D0PVafxd, ff20exp, ff20miexp, ff20plexp}, PlotRange -> {{2.8, 12.2}, {0.23, 0.35}}, PlotStyle -> {Blue, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(a)", 13], Scaled[{0.70, 0.90}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.08, 0.85}]], Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]] == Round[10000 alsmtauPV[8]]/10000., 13], Scaled[{0.75, 0.28}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.15, 0.40}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"PV (D=0)"}, {0.77, 0.40}]]; (* OK *)

a21PVafxd= ListPlot[{Tab21D0PVafxd, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.19, 0.26}}, PlotStyle -> {Blue, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}, {Blue, PointSize[0.001]}}, Filling -> {3 -> {4}}, AxesLabel -> {Subscript[N, t]}, Joined -> {True, True, True, True}, PlotMarkers -> {\[FilledCircle], \[FilledVerySmallSquare], \[FilledVerySmallSquare], \[FilledVerySmallSquare]}, Epilog -> {Inset[Style["(b)", 13], Scaled[{0.50, 0.90}]], Inset[Style[Superscript[(a), "(2,1)"], 13], Scaled[{0.18, 0.95}]], Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]] == Round[10000 alsmtauPV[8]]/10000., 13], Scaled[{0.80, 0.78}]], Inset[Style[Subscript[d, 4] == ddtk[4, Exp[tKK]], 13], Scaled[{0.85, 0.50}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"PV (D=0)"}, {0.85, 0.90}]]; (* OK *)

    (*------*)

 (*------*)

  (* end of the NEW (O61, O62) PV approach *)
(*------------------------------------*) 
(*---------------------------------------*)
(* global PV: The above behaviour of a^{(2,1)}(aAL,O61,O62) and a^{(2,0)}(aAL,O61,O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax)=aALfPV[Nt] to a(Q^2=sigmax Exp[I fi]) ] suggests that the truncation order Nt=8 in the PV(+CI) global approach (i.e., the approach of minimization with alpha_s variable) is Nt=8. (note, we have here the case: ddtk[4,Exp[tKK]=275.) OK *)

 (* the central value of the PV global fit result, at Nt=8, is:
aAL=aALfPV[8]=0.1079; O4O612fitPPVMinnum[8]={0.00310712, {O4fitf -> -0.00455454, O61fitf -> -0.00673939, O62fitf -> 0.000892787}};  alpha_s(mtau^2)(=alsmtauPV[8])=0.32687; OK *)
 
  (* Now the estimate of the EXPERIMENTAL uncertainty in this global PV approach with double-pinched tFO Borel-Laplace (with Nt=6): *)

(* 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 estimation of delO4 with the (NdxNd) covariance matrix  ('varcov')
  (and using PV Nt=8 global analysis for theoretical values): *)
Nt0=8;
 (* Use the Borel transform functions, etc., generated with aAL=aALfPV[Nt0] *)
 (* We vary O4 in UNCORRELATED chi ('var'): *)
chi2o2O4[O4f_]:=Sum[(ReBPPV[M2covF[jj],Nt0,O4f,O61PV[Nt0],O62PV[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]];
 (* O4PVvarfit=-0.00455414; very close to ,O4PV[Nt0]=-0.00455454, as it should be, OK *)
 delO4PV= delO4f/.FindRoot[chi2o2O4[O4PVvarfit+delO4f]==Xi2O4PVvaro77+1,{delO4f, 3. 10^-5}];
 (* --> for Nm=2 (Nd=6) we obtain: delO4PV=0.0000304; OK *)

chi2o2O61[O61f_]:=Sum[(ReBPPV[M2covF[jj],Nt0,O4PV[Nt0],O61f,O62PV[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O61PVvarfitMin=FindMinimum[chi2o2O61[O61fitf],{O61fitf,O61PV[Nt0]}];
O61PVvarfit=O61fitf/.O61PVvarfitMin[[2]];
Xi2O61PVvaro77=O61PVvarfitMin[[1]];
 (* Note that O61PVvarfit=-0.0067395, very close to O61PV[Nt0]=-0.00673939, as it should be, OK *)
 delO61PV= delO61f/.FindRoot[chi2o2O61[O61PVvarfit+delO61f]==Xi2O61PVvaro77+1,{delO61f, 3. 10^-5}];
(* Nm=2; Nd=6; delO61PV=0.0000451; OK *)

chi2o2O62[O62f_]:=Sum[(ReBPPV[M2covF[jj],Nt0,O4PV[Nt0],O61PV[Nt0],O62f] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}];
O62PVvarfitMin=FindMinimum[chi2o2O62[O62fitf],{O62fitf,O62PV[Nt0]}];
O62PVvarfit=O62fitf/.O62PVvarfitMin[[2]];
Xi2O62PVvaro77=O62PVvarfitMin[[1]];
 (* Note that O62PVvarfit=0.000892796, very close to O62PV[Nt0]=0.000892787, as it should be, OK *)
 delO62PV= delO62f/.FindRoot[chi2o2O62[O62PVvarfit+delO62f]==Xi2O62PVvaro77+1,{delO62f, 3. 10^-5}];
 (* Nm=2; Nd=6; delO8PV=5.02047*10^-6; OK *)

 aAL=0.1079; (*= aALfPV[Nt0], Nt0=8, OK *)
 chi2o2aAL=Sum[(ReBPPV[M2covF[jj],Nt0,O4PV[Nt0],O61PV[Nt0],O62PV[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}]; (* -> chi2o2aAL=0.00158713; OK *)
 aAL=0.1079 + 0.000112: chi2o2aAL=1.013;  (* OK *)
aAL=0.1079 - 0.000109: chi2o2aAL=1.002;   (* OK *)
(* --> delaALPV=0.00020; *)
(* Nm=2; Nd=6; Nt0=8; delaALPV=+0.000112; -0.000109; OK *)
(*-------*)
(* Together, the estimates of EXPERIMENTAL uncertainties in the PV global approach (with Nt=8) are (Nm=2; Nd=6), in PV global:  *)
(* delO4PV= 0.0000303853; delO61PV=0.0000450575; delO62PV=5.02047*10^-6; delaALPV=0.00011; (=delta(a(sigmax)))  OK --> *)

(*
alpha_s(mtau^2)=0.32687+{^{+0.000318}_{-0.000318}_{exp};
10^3 O4= -4.55+{^{+0.03}_{-0.03}}_{exp} =app -4.6+{^{+0.0}_{-0.0}}_{exp};
10^3 O61=-6.74+{^{+0.05}_{-0.05}}_{exp} =app -6.7+{^{+0.0}_{-0.0}}_{exp};
10^4 O62=+8.93+{^{+0.05}_{-0.05}}_{exp} =app +8.9+{^{+0.1}_{-0.1}}_{exp};
OK *)

(* The precise central value result of the PV global fit:
 (aAL=)aALfPV[8]=0.1079; O4O612fitPPVMinnum[8]={0.00310712, {O4fitf -> -0.00455454, O61fitf -> -0.00673939, O62fitf -> 0.000892787}}; OK *)
 
(*--------------*)
(* Variation of kappa, PV global: *)
(* kappa=2.0: 
 aAL=0.0894;  O4O612fitPPVMinnum[8]={0.00117666, {O4fitf -> -0.00491979, O61fitf -> -0.00514513, O62fitf -> 0.000758806}}; alpha_s(mtau^2)=0.327562; (* OK *)
kappa=1.0 (the central case of the global fit): 
aAL=0.1079; O4O612fitPPVMinnum[8]={0.00310712, {O4fitf -> -0.00455454, O61fitf -> -0.00673939, O62fitf -> 0.000892787}}; alpha_s(mtau^2)=0.32687; (* OK *)
kappa=2/3:
 aAL=0.1288;O4O612fitPPVMinN[8]={0.496186, {O4fitf -> -0.00770929, O61fitf -> -0.0186842, O62fitf -> 0.00248721}}; alpha_s(mtau^2)=0.337111; (* OK *) 
 ---> *)

 (*
alpha_s(mtau^2)=0.32687+{^{+0.000692}_{+0.010241}}_{kap};
10^3 O4=  -4.6+{^{-0.3}_{-3.1}}_{kap};
10^3 O61= -6.7+{^{+1.6}_{-12.0}}_{kap};
10^4 O62=+8.9+{^{-1.3}_{+16.0}}_{kap};
 OK *)
(*----------------*)
(* Variation of d4, PV global (d4=338.19+-63.19): 
When ddtk[4,tkap]=227+63.19=338.19 (delddtk[4, tkap]=1207.16):
 aAL=0.1057; O4O612fitPPVMinnum[8]={0.530249, {O4fitf -> -0.00333074, O61fitf -> -0.00846099, O62fitf -> 0.00135968}}; alpha_s(mtau^2)=0.320507; (* OK *)
 Central case: ddtk[4,tkap]=275. (delddtk[4, tkap]=-1016.32; kappa=1):
 aAL=aALfPV[8]=0.1079; O4O612fitPPVMinnum[8]={0.00310712, {O4fitf -> -0.00455454, O61fitf -> -0.00673939, O62fitf -> 0.000892787}}; alpha_s(mtau^2)=0.32687; (* OK *)
 When ddtk[4,tkap]=227-63.19=211.81 (delddtk[4, tkap]=-103.10):
 aAL=0.1132; O4O612fitPPVMinnum[8]={0.00244176, {O4fitf -> -0.00732388, O61fitf -> -0.00395464, O62fitf -> 0.00117182}}; alpha_s(mtau^2)=0.342132; (* OK *)
 (*
alpha_s(mtau^2)=0.32687+{^{-0.006363}_{+0.015262}}_{d4};
10^3 O4=  -4.6+{^{+1.3}_{-2.7}}_{d4};
10^3 O61= -6.7+{^{-1.8}_{+2.7}}_{d4};
10^4 O62=+8.9+{^{+4.7}_{+2.8}}_{d4};
  OK *) 
  (*---------------------*)
 (* Variation of Nt (truncation order, Nt=Nt0+2-3=8+2-3, i.e., variation Nt=10,9, 8, 7,6,5; PV global with (O61, O62); *)
(*
aALfPV[10]=0.1111; O4O612fitPPVMinnum[10]={0.390203, {O4fitf -> -0.00741201, O61fitf -> -0.0157418, O62fitf -> 0.00208312}}; alpha_s(mtau^2)=0.336096; (* OK *)
aALfPV[9]=0.1080; O4O612fitPPVMinnum[9]={0.00287242, {O4fitf -> -0.00463592, O61fitf -> -0.00695691, O62fitf -> 0.000921253}}; alpha_s(mtau^2)=0.327159;(* OK *)
aALfPV[8]=0.1079; O4O612fitPPVMinnum[8]={0.00310712, {O4fitf -> -0.00455454, O61fitf -> -0.00673939, O62fitf -> 0.000892787}}; alpha_s(mtau^2)=0.32687; (* OK *)
aALfPV[7]=0.1077; O4O612fitPPVMinnum[7]= {0.00356788, {O4fitf -> -0.00437873, O61fitf -> -0.00631874, O62fitf -> 0.00083764}}; alpha_s(mtau^2)=0.326292;(* OK *)
aALfPV[6]=0.1083; O4O612fitPPVMinnum[6]={0.00195454, {O4fitf -> -0.00516294, O61fitf -> -0.00800203, O62fitf -> 0.00105598}}; alpha_s(mtau^2)=0.328025; (* OK *)
aALfPV[5]=0.1082; O4O612fitPPVMinnum[5]={0.000464085, {O4fitf -> -0.00616911, O61fitf -> -0.00968949, O62fitf -> 0.00127478}}; alpha_s(mtau^2)=0.327737; (* OK *)
*) 

 (*
alpha_s(mtau^2)=0.32687+{^{+0.009226}_{+0.001155}}_{Nt};
10^3 O4=  -4.6+{^{-2.8}_{-1.5}}_{Nt};
10^3 O61= -6.7+{^{-9.0}_{-3.0}}_{Nt};
10^4 O62= +8.9+{^{+11.9}_{+3.8}}_{Nt};
 Note: Nt=8^{+2}_{-3}
 OK  *)   
(*-------------------------*)
     (* Variation due to the renormalon ambiguity ('amb') to Adler: in ReBPPV we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn] + BPPVdelD0[M2f]) *)
(* Nt=8; O4O6O8fitPPVMinN[8] gives the following xi^2 =xi2[8]: *)
(*
 aAL=0.1140: xi2[8]=0.00129068;
 aAL=0.1138: xi2[8]=0.00112753;
 aAL=0.1137: xi2[8]=0.00080744;  ***
 aAL=0.1136: xi2[8]=0.000985077;
 aAL=0.1135: xi2[6]=0.00127044;
  *)
 (* Therefore: aAL=0.1037 where xi2[8] is minimal: *)
 (* BBPVD0[M2f,Nn] -> (BPPVD0[M2f, Nn] + BPPVdelD0[M2f]); Nt=8:
  aAL=0.1137;  O4O6O8fitPPVMinN[8]={0.00080744, {O4fitf -> -0.00947019, O61fitf -> -0.0182213, O62fitf -> 0.00256145}}; alpha_s(mtau^2)=0.343567; OK *)
(*---*)
 (* The central value  (BBPVD0[M2f,Nn]): 
 aAL=0.1079; O4O612fitPPVMinnum[8]={0.00310712, {O4fitf -> -0.00455454, O61fitf -> -0.00673939, O62fitf -> 0.000892787}}; alpha_s(mtau^2)=0.32687; (* OK *)
  (*---*)
(* When  in ReBPPV we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn] - BPPVdelD0[M2f]) *)
(* Nt=8; O4O6O8fitPPVMinN[8] gives the following xi^2 =xi2[8]: *)
(*
 aAL=0.1050: xi2[8]=0.0194491;
 aAL=0.1054: xi2[8]=0.00943972;
 aAL=0.1055: xi2[8]=0.00857864;
 aAL=0.1056: xi2[8]=0.00839214; ***
 aAL=0.1057: xi2[8]=0.00889092;
 aAL=0.1058: xi2[8]=0.010152;
 aAL=0.1060: xi2[8]=0.0173995;
  *)
 (* Therefore: aAL=0.1056 where xi2[8] is minimal: *)
(* BBPVD0[M2f,Nn] -> (BPPVD0[M2f, Nn] - BPPVdelD0[M2f]); Nt=6:
 aAL=0.1056;  O4O6O8fitPPVMinN[8]={0.00839214, {O4fitf -> -0.00255204, O61fitf -> -0.00242706, O62fitf -> 0.000301289}}; alpha_s(mtau^2)=0.320218; OK *)
  (*---*)
   (*------*)
(* --> aAL=0.1079 + {^{+0.0058}_{-0.0033}}_{amb} OK  --> *)     
(*
alpha_s(mtau^2)=0.32687+{^{+0.016697}_{-0.006652}}_{amb};
10^3 O4=  -4.6+{^{-4.9}_{+2.0}}_{amb};
10^3 O61= -6.7+{^{-11.5}_{+4.3}_{amb};
10^4 O62= +8.9+{^{+16.7}_{-5.9}}_{amb};
OK
*)
 (*-----------------------------------------*)
 (* Altogether, the PV(+CI) global fit approach with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax)=aALftFO[Nt] to a(Q^2=sigmax Exp[I fi]) ], gives: *)
  (* PV(+CI)
   Nt=8;  ( ddtk[4, Exp[tKK]]=275.);
   d4=275.+-63.19; kap=1^{+1}_{-1/3};  Nt=8^{+2}-{-3};
   alpha_s(mtau^2)=0.32687+{^{+0.000318}_{-0.000318}_{exp}+{^{+0.000692}_{+0.010241}}_{kap}+{^{-0.006363}_{+0.015262}}_{d4}+{^{+0.009226}_{+0.001155}}_{Nt}+{^{+0.016697}_{-0.006652}}_{amb} = 0.32687+{^{+0.026492}_{-0.009211}} =app 0.3269+{^{+0.0265}_{-0.0092}} =app 0.327+{^{+0.026}_{-0.009}} OK
10^3 O4=-4.6+{^{+0.0}_{-0.0}}_{exp}+{^{-0.3}_{-3.1}}_{kap}+{^{+1.3}_{-2.7}}_{d4}+{^{-2.8}_{-1.5}}_{Nt}+{^{-4.9}_{+2.0}}_{amb} = -4.6+{^{+2.4}_{-7.0}}  OK 
10^3 O61=-6.7+{^{+0.0}_{-0.0}}_{exp}+{^{+1.6}_{-12.0}}_{kap}+{^{-1.8}_{+2.7}}_{d4}+{^{-9.0}_{-3.0}}_{Nt}+{^{-11.5}_{+4.3}_{amb}  = -6.7+{^{+5.3}_{-19.0}} OK
10^4 O62=+8.9+{^{+0.1}_{-0.1}}_{exp}+{^{-1.3}_{+16.0}}_{kap}+{^{+4.7}_{+2.8}}_{d4}+{^{+11.9}_{+3.8}}_{Nt}+{^{+16.7}_{-5.9}}_{amb} =+8.9+{^{+26.4}_{-6.0}} OK
OK  *)
 (* End of PV(+CI) global fit results with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax)=aALftFO[Nt] to a(Q^2=sigmax Exp[I fi]) ] *)
(*------------------------------------*)
(*------------------------------------*)
(*------------------------------------*)
 (* To relate alpha_s(MZ^2) with alpha_s(mtau^2), via 5-loop MSbar RGE: *)
 (*
mtau=1.77686; mc=1.27; mb=4.20; kappa=2.; Q=mtau;
alphaMSMZ=0.11...; <</home/gorazd/3danQCD.dir/RGErundown5l.m;  aaMSQ*Pi
 Then: aaMSQ*Pi is the value of alpha_s(mtau^2) *)
