(* SumRMSbALEPHO6ka12.m: *)
(* SumRMSbALEPHO6ka12.m: as SumMSbALEPHM2O61O62.m (d4MSbar=275.+ -63.), but now the IR p=3 part of B[td][u]  has 1/(3-u)^ka2 and 1/(3-u)^ka1, where: ka2= 0.778 =app (0.8024+0.7527)/2, and ka1 = 0.375 =app (0.4173+0.3320)/2 (as suggested by the work of BHJ (Boito, Hornung and Jamin, 1510.03812v2); in SumMSbALEPHM2O61O62.m, we had large-beta0 values ka2=2 and ka1=1 *)
(* SumMSbALEPHM2O61O62.m: as SumRMSbALEPHM2O61O62d4338.m, but for: d4MSbar=275.+ -63.19, and not 338.19+-63.19 *)
(* Note, the string 'OK' indicates that the part there has been converted to the deMSbar=275 case (from d4MSbar=338.19 case) *)
(*-------------------------------------------------*)
(* SumRMSbALEPHM2O61O62d4338.m: as SumRPMSbALEPHM2m090.m, but instead of O6 and O8, we have O6^{(1)} and O6^{(2)} condensates *)
(* Look for: NEW *)
(*-------------------------------------------------*)
(* SumRPMSbALEPHM2m090.m:  with (double) pinched ('P') Borel-Laplace sum rules *)
(* and alpha_s(mtau^2; MSbar) is adjustable (as input, or as result fo global fit) *)
(*-------------------------------------------------*)
(* Borel-Laplace fitted in |M^2|-intervals: (M2min,1.5); M2min=0.65 or 0.90; *)
(* here: M2min=0.90, i.e., YYY=090 *)
(* in pQCD, MSbar scheme *)

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

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

 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275 (ECH-predicted): *)
<</home/gorazd/3danQCD.dir/MSbarRenMod5Ad4275Up.save; 

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

  (*---------*)
 
BtDUV1[u_, tkapp_] := Exp[Log[tkapp]*u]*Pi*(tdUVp12/(1 + u)^2);
BtDIR3[u_, tkapp_] := 
  Exp[Log[tkapp]*u]*Pi*(tdIRp3ka2/(3 - u)^ka2 + tdIRp3ka1/(3 - u)^ka1);
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;

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 *)
 BPFOD0[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 *)

  (*-----------*) 
(* Finding of Re[B] in FO approach; *)
(* Now the dfi integration in D=6 term with a(sigmax Exp[I fi])^(1 -ka2) and a(sigmax Exp[I fi])^(1 -ka1) 5-loop CI(RGE)-evolved (' Num5;'):  *)
sol1a[aALf_]:= NDSolve[{-I y'[th] == bta[y[th]], y[0] == aALf}, y, {th, -Pi, +Pi}];
 aCI[th_, aALf_] := First[y[th] /. sol1a[aALf]];
(* aALf=a(kappa sigmax); aCI[th,aALf]=a(Q2=kappa sigma Exp[I th]) ' *)
  
(* Now the dfi integration in D=6 term with 1/a(sigmax Exp[I fi]) 5-loop CI(RGE)-evolved; note that aALf = a(kappa sigmax), and is used in the D=0 part; so a(sigmax)=aMSext[(1/kappa), aALf];  in D=6 we use a(sigmax Exp[I fi]) = aCI[fi, a(sigmax)] = aCI[fi, aMSext[(1/kappa), aALf]] *)
ReBPFONum5l[M2f_,Nn_,aALf_,O4_,O61_,O62_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPFOD0[M2f,Nn,aALf]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax ) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O62 aCI[fi,aMSext[(1/kappa), aALf]]^(1-ka2) + O61 aCI[fi,aMSext[(1/kappa), aALf]]^(1-ka1)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]]; 

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

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


(*-------------------*)
(* Now the same, but with tFO in the D=0 part: *)
ReBPtFONum5l[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 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax ) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O62 aCI[fi,aMSext[(1/kappa), aALf]]^(1-ka2) + O61 aCI[fi,aMSext[(1/kappa), aALf]]^(1-ka1)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];

(* Now the simplified ('s') form, which is valid only for kappa=1: *)
ReBPtFONum5ls[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 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax ) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O62 aCI[fi,aALf]^(1-ka2) + O61 aCI[fi,aALf]^(1-ka1)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];			     

(* Now the form in D=6 with 1-loop running in D=6 ('1l'); when kappa not=1, RGE from aALf=a(kappa sigmax) to a(sigmax) is performed by 5-loop RGE: aALf=a(kappa sigmax); a(sigmax)=aMSext[(1/kappa),aALf]; (note: aALf=a(kappa sigmax) in this general case): *)
ReBPtFONum1l[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 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax ) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O62 (aMSext[(1/kappa), aALf]/(1 + I be0 aMSext[(1/kappa), aALf] fi))^(1-ka2) + O61 (aMSext[(1/kappa), aALf]/(1 + I be0 aMSext[(1/kappa), aALf] fi))^(1-ka1)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];			     
(* Now the simplified (faster) form in D=6 and 1-loop running in D=6 ('1ls'), which is valid only for kappa=1: (note: aALf=a(sigmax) in this case)*)
ReBPtFONum1ls[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 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax ) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O62 (aALf/(1 + I aALf be0 fi))^(1-ka2) + O61 (aALf/(1 + I aALf be0 fi))^(1-ka1)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]]; 			     
(*-----------*) 

 
(* Now 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.003; O61init=-0.03; O62init=0.03; aALfinit=0.105;
(* FO approach: *)
 alO4O612fitPFOMinNum5l[Nn_]:= FindMinimum[Sum[(ReBPFONum5l[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[(ReBPFONum5l[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[(ReBPFONum5l[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 simplified form ('s'), which is valid only for kappa=1: *)
 alO4O612fitPFOMinNum5ls[Nn_]:= FindMinimum[Sum[(ReBPFONum5ls[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[(ReBPFONum5ls[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[(ReBPFONum5ls[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 form ('1l') with 1-loop RGE running in D=6, valiud for general kappa: *)
alO4O612fitPFOMinNum1l[Nn_]:= FindMinimum[Sum[(ReBPFONum1l[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[(ReBPFONum1l[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[(ReBPFONum1l[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 simplified form ('1ls') with 1-loop RGE running in D=6, which is valid only for kappa=1: *)
alO4O612fitPFOMinNum1ls[Nn_]:= FindMinimum[Sum[(ReBPFONum1ls[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[(ReBPFONum1ls[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[(ReBPFONum1ls[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' 5-loop (CI) along the contour in D=6 part; 'nnum1l' runs 1-loop (CI) along the contour in D=6 part; *)
 (* in D=0, 5-loop FO running *)
(*
 alO4O612fitPFOMinnum1l[3]={0.000670097, {aALf -> 0.109353, O4fitf -> 0.00118977, O61fitf -> 0.00273347, O62fitf -> -0.00136823}};
 alO4O612fitPFOMinnum1l[4]={0.0000302671, {aALf -> 0.105136, O4fitf -> 0.000213148, O61fitf -> -0.00802502, O62fitf -> 0.00406987}};
 alO4O612fitPFOMinnum1l[5]={0.000784024, {aALf -> 0.103896, O4fitf -> -0.000754103, O61fitf -> -0.0177607, O62fitf -> 0.00894318}};
 alO4O612fitPFOMinnum1l[6]={0.0018433, {aALf -> 0.103745, O4fitf -> -0.00172053, O61fitf -> -0.0256395, O62fitf -> 0.0129543}};
 alO4O612fitPFOMinnum1l[7]={0.00272373, {aALf -> 0.104134, O4fitf -> -0.00274006, O61fitf -> -0.0317815, O62fitf -> 0.016188}};
 alO4O612fitPFOMinnum1l[8]={0.0031155, {aALf -> 0.104745, O4fitf -> -0.00384881, O61fitf -> -0.0360145, O62fitf -> 0.0185475}};
 alO4O612fitPFOMinnum1l[9]={0.00297135, {aALf -> 0.105435, O4fitf -> -0.00512041, O61fitf -> -0.0382952, O62fitf -> 0.0199889}};
 alO4O612fitPFOMinnum1l[10]={0.00227335, {aALf -> 0.105879, O4fitf -> -0.00654659, O61fitf -> -0.0380327, O62fitf -> 0.0201438}};
 alO4O612fitPFOMinnum1l[11]={0.00164024, {aALf -> 0.106256, O4fitf -> -0.00844986, O61fitf -> -0.036347, O62fitf -> 0.0195712}}; 
 alO4O612fitPFOMinnum1l[12]={0.000926816, {aALf -> 0.105807, O4fitf -> -0.0104797, O61fitf -> -0.032101, O62fitf -> 0.0175266}};
 (*----*)
 
 alO4O612fitPFOMinnum5l[3]={0.00075164, {aALf -> 0.109409, O4fitf -> 0.00116851, O61fitf -> 0.0017225, O62fitf -> -0.000912482}};
 alO4O612fitPFOMinnum5l[4]={0.0000822773, {aALf -> 0.105017, O4fitf -> 0.000269104, O61fitf -> -0.00526165, O62fitf -> 0.00284666}};
 alO4O612fitPFOMinnum5l[5]={0.00123414, {aALf -> 0.103656, O4fitf -> -0.000624974, O61fitf -> -0.0116732, O62fitf -> 0.00625624}};
 alO4O612fitPFOMinnum5l[6]={0.00278356, {aALf -> 0.103406, O4fitf -> -0.00151825, O61fitf -> -0.0168272, O62fitf -> 0.00905427}};
  alO4O612fitPFOMinnum5l[7]={0.0041179, {aALf -> 0.103703, O4fitf -> -0.00245914, O61fitf -> -0.0207782, O62fitf -> 0.0112905}};
  alO4O612fitPFOMinnum5l[8]={0.00482343, {aALf -> 0.104231, O4fitf -> -0.00348281, O61fitf -> -0.0234481, O62fitf -> 0.0129149}};
  alO4O612fitPFOMinnum5l[9]={0.00480678, {aALf -> 0.10485, O4fitf -> -0.00466003, O61fitf -> -0.0248577, O62fitf -> 0.0139205}};
  alO4O612fitPFOMinnum5l[10]={0.00396414, {aALf -> 0.10527, O4fitf -> -0.00600092, O61fitf -> -0.0247294, O62fitf -> 0.0141043}};
 alO4O612fitPFOMinnum5l[11]={0.00310521, {aALf -> 0.105646, O4fitf -> -0.00780071, O61fitf -> -0.023748, O62fitf -> 0.0138261}};
 alO4O612fitPFOMinnum5l[12]={0.00196059, {aALf -> 0.105301, O4fitf -> -0.00980948, O61fitf -> -0.0212736, O62fitf -> 0.0126097}};    
 
Clear[Nn]; Do[Module[{}, xi2FO[Nn]=alO4O612fitPFOMinnum1l[Nn][[1]]; aALfFO[Nn]=aALf/.alO4O612fitPFOMinnum1l[Nn][[2]]; alsmtauFO[Nn]=aMSext[mtauphys^2/(kappa sigmax),aALfFO[Nn]] Pi;  O4FO[Nn]=O4fitf/.alO4O612fitPFOMinnum1l[Nn][[2]]; O61FO[Nn]=O61fitf/.alO4O612fitPFOMinnum1l[Nn][[2]]; O62FO[Nn]=O62fitf/.alO4O612fitPFOMinnum1l[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.: alO4O612fitPFOMinnum1l[8]={0.0031155, {aALf -> 0.104745, O4fitf -> -0.00384881, O61fitf -> -0.0360145, O62fitf -> 0.0185475}}; alpha_s(mtau^2)=alsmtauFO[8]=0.31774;

  aaFO21[0.106507][8, -0.00275591, 0.000710728]=0.204403;
 (* aexpBin77[2,1]+-delaexpBin77[2,1]=0.206095+-0.00260488 = 0.2061+-0.0026 *)
  aaFO20[aALfFO[8]][8, O4FO[8], O61FO[8], O62FO[8]]=0.239693;
 (* aexpBin77[2,0]+-delaexpBin77[2,0]=0.240381+-0.00261244= 0.2404+-0.0026 *)
  (* 6 Apr. 2022 *)
  *)
  (*----------------------------------------------------*)			        
(* 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.003; O61init=-0.03; O62init=0.02; aALfinit=0.105;
alO4O612fitPtFOMinNum5l[Nn_]:= FindMinimum[Sum[(ReBPtFONum5l[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[(ReBPtFONum5l[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[(ReBPtFONum5l[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 simplified form ('s'), which is valid only for kappa=1: *)
alO4O612fitPtFOMinNum5ls[Nn_]:= FindMinimum[Sum[(ReBPtFONum5ls[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[(ReBPtFONum5ls[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[(ReBPtFONum5ls[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}}];
 
 (* Do[alO4O612fitPtFOMinnum[Nn]=alO4O612fitPtFOMinN[Nn],{Nn,3,12}]; *)

 (* Now the form ('1l') with 1-loop RGE running in D=6, valiud for general kappa: *)
alO4O612fitPtFOMinNum1l[Nn_]:= FindMinimum[Sum[(ReBPtFONum1l[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[(ReBPtFONum1l[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[(ReBPtFONum1l[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 simplified form ('1ls') with 1-loop RGE running in D=6, which is valid only for kappa=1: *)
alO4O612fitPtFOMinNum1ls[Nn_]:= FindMinimum[Sum[(ReBPtFONum1ls[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[(ReBPtFONum1ls[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[(ReBPtFONum1ls[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 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 *)

 (*
  alO4O612fitPtFOMinnum1l[3]={0.0000269251, {aALf -> 0.103235, O4fitf -> 0.000676531, O61fitf -> -0.00361756, O62fitf -> 0.00181203}};
  alO4O612fitPtFOMinnum1l[4]={0.00178082, {aALf -> 0.10223, O4fitf -> -0.00103472, O61fitf -> -0.0228097, O62fitf -> 0.0113147}};
  alO4O612fitPtFOMinnum1l[5]={0.00542654, {aALf -> 0.104284, O4fitf -> -0.00318952, O61fitf -> -0.039268, O62fitf -> 0.0199033}};
  alO4O612fitPtFOMinnum1l[6]={0.00435203, {aALf -> 0.106112, O4fitf -> -0.00541216, O61fitf -> -0.0451519, O62fitf -> 0.0235571}};
 alO4O612fitPtFOMinnum1l[7]={0.00187271, {aALf -> 0.108184, O4fitf -> -0.00857779, O61fitf -> -0.0439567, O62fitf -> 0.0239095}};
  alO4O612fitPtFOMinnum1l[8]={0.0000168526, {aALf -> 0.104803, O4fitf -> -0.00834969, O61fitf -> -0.026698, O62fitf -> 0.0148481}};
  alO4O612fitPtFOMinnum1l[9]={2.25587, {aALf -> 0.104716, O4fitf -> -0.0108699, 
  O61fitf -> 0.0307139, O62fitf -> -0.0116543}};
  alO4O612fitPtFOMinnum1l[10]={0.00127956, {aALf -> 0.0964448, O4fitf -> -0.00651692, O61fitf -> -0.00275088, O62fitf -> 0.00219543}};
  alO4O612fitPtFOMinnum1l[11]={16.5689, {aALf -> 0.084939, O4fitf -> 0.00446325, O61fitf -> 0.24232, O62fitf -> -0.10045}};
 alO4O612fitPtFOMinnum1l[12]={0.0103765, {aALf -> 0.085722, O4fitf -> -0.00172608, O61fitf -> 0.028309, O62fitf -> -0.0120207}};
  

  alO4O612fitPtFOMinnum5l[3]={0.0000136725, {aALf -> 0.103189, O4fitf -> 0.000699367, O61fitf -> -0.00239949, O62fitf -> 0.00127906}};
alO4O612fitPtFOMinnum5l[4]={0.00254057, {aALf -> 0.101961, O4fitf -> -0.000871933, O61fitf -> -0.015049, O62fitf -> 0.00791375}};
  alO4O612fitPtFOMinnum5l[5]={0.00764821, {aALf -> 0.10374, O4fitf -> -0.00281584, O61fitf -> -0.0254286, O62fitf -> 0.0137147}};
  alO4O612fitPtFOMinnum5l[6]={0.00688282, {aALf -> 0.105362, O4fitf -> -0.00481765, O61fitf -> -0.0289106, O62fitf -> 0.0161532}};
  alO4O612fitPtFOMinnum5l[7]={0.00393308, {aALf -> 0.107223, O4fitf -> -0.00763533, O61fitf -> -0.0280621, O62fitf -> 0.0164848}};
 alO4O612fitPtFOMinnum5l[8]={0.000262874, {aALf -> 0.104442, O4fitf -> -0.00791187, O61fitf -> -0.0180994, O62fitf -> 0.0109788}};
  alO4O612fitPtFOMinnum5l[9]={2.59921, {aALf -> 0.104832, O4fitf -> -0.0111193, 
  O61fitf -> 0.0209544, O62fitf -> -0.00747601}};
  alO4O612fitPtFOMinnum5l[10]={0.00122723, {aALf -> 0.096436, O4fitf -> -0.00650071, O61fitf -> -0.00227198, O62fitf -> 0.00199077}};  
  alO4O612fitPtFOMinnum5l[11]={19.3775, {aALf -> 0.084728, O4fitf -> 0.00428554, O61fitf -> 0.190217, O62fitf -> -0.0803515}};
 alO4O612fitPtFOMinnum5l[12]={0.0116753, {aALf -> 0.0857676, O4fitf -> -0.0018336, O61fitf -> 0.0205285, O62fitf -> -0.00895125}};


 Clear[Nn]; Do[Module[{}, xi2tFO[Nn]=alO4O612fitPtFOMinnum1l[Nn][[1]]; aALftFO[Nn]=aALf/.alO4O612fitPtFOMinnum1l[Nn][[2]]; alsmtautFO[Nn]=aMSext[mtauphys^2/(kappa sigmax),aALftFO[Nn]] Pi;  O4tFO[Nn]=O4fitf/.alO4O612fitPtFOMinnum1l[Nn][[2]]; O61tFO[Nn]=O61fitf/.alO4O612fitPtFOMinnum1l[Nn][[2]]; O62tFO[Nn]=O62fitf/.alO4O612fitPtFOMinnum1l[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): alO4O612fitPtFOMinnum1l[5]={0.00542654, {aALf -> 0.104284, O4fitf -> -0.00318952, O61fitf -> -0.039268, O62fitf -> 0.0199033}}; alpha_s(mtau^2)=alsmtautFO[5]=0.316404;
   aatFO21[aALftFO[NNt]][NNt, O61tFO[NNt], O62tFO[NNt]]= 0.204745(N=4); 0.203809(N=5);  0.204161(N=6); 0.204912(N=7); 0.206072(N=8); 
  (* aexpBin77[2,1]+-delaexpBin77[2,1]=0.206095+-0.00260488 = 0.2061+-0.0026 *)
  aatFO20[aALftFO[NNt]][NNt, O4tFO[NNt], O61tFO[NNt], O62tFO[NNt]]=0.23984(N=4); 0.239456(N=5); 0.239588(N=6); 0.239893(N=7); 0.24037(N=8);
   (* aexpBin77[2,0]+-delaexpBin77[2,0]=0.240381+-0.00261244= 0.2404+-0.0026 *)
 (* 7 Apr. 2022 *)
 
   *)
   (*---------------------------------------------------------------------_*)
  

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

Plot0FOo77= 
Plot[{ReBPexpBinAL[M2] + delReBPexpBinAL[M2], 
  ReBPexpBinAL[M2] - delReBPexpBinAL[M2], 
  ReBPexpBinAL[M2],
  ReBPFONum1l[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, Thickness[0.001]}, {Black, 
   Thickness[0.001]}, {Black, Thickness[0.001]}, {Red, Dashed, 
   Thickness[0.001]}},(*{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]],
    ReBPFONum1l[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, Thickness[0.001]}, {Black, 
   Thickness[0.001]}, {Black, Thickness[0.001]}, {Red, Dashed, 
   Thickness[0.001]}},(*{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]], 
   ReBPFONum1l[M2 Exp[I Pi/4], 8, aALfFO[8], O4FO[8], O61FO[8], 
    O62FO[8]]}, {M2, M2min, 1.5}, 
  PlotRange -> {{M2min, 1.5}, {0.49, 0.72}}, 
  AxesOrigin -> {M2min, 1.09}, 
   PlotStyle -> {{Black, Thickness[0.001]}, {Black, 
   Thickness[0.001]}, {Black, Thickness[0.001]}, {Red, Dashed, 
   Thickness[0.001]}},(*{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["FigPsi0PFOo77O6ka1ka2.pdf",Plot0FOo77];  Export["FigPsiPi6PFOo77O6ka1ka2.pdf",PlotPi6FOo77];  Export["FigPsiPi4PFOo77O6ka1ka2.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] *)
aaFOD0[2,n_,aALf_][NN_]:=Sum[coeffaaFO[2,n][nn] x^nn,{nn,1,NN}]/.{x->aALf};
(* aALf=a(kappa*sigmax) *)
(*--------------*)
(*  This is the NEW (6 Apr.2022) in FO approach with (O6ka1=O61, O6ka2=O62): *)
(* a21 and a20, this is the case when a(sigmax Exp[I fi]) in [O_6ka a(sigmax Exp[I fi])^(1-ka)] runs by one-loop RGE: *)
aaFO21[aALf_][NN_,O61_,O62_]:=aaFOD0[2,1,aALf][NN]+ 2 (1/2/Pi) (6 Pi^2/sigmax^3)*( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 - Exp[2 I fi]) + (1/2) (1 - Exp[I 4 fi]) ) (O62 (aALf/(1 + I aALf be0 fi))^(1-ka2) + O61 (aALf/(1 + I aALf be0 fi))^(1-ka1)),{fi,-Pi,Pi}] );  (* this is a^{(2,1)} in the approach with O61=O6ka1 and O62=O6ka2, and a(sigmax Exp[I fi]) in O6ka (a(sigmax Exp[I fi]))^(1-ka) runs by one-loop RGE; note that aALf=a(kappa sigmax); but the higher-twist (D=6) terms are evaluated with the reasonable a(sigmax Exp[I fi])^(1-ka) (and not: a(kappa sigmax Exp[I fi])^(1-ka)) only if kappa=1.   *)
aaFO20[aALf_][NN_,O4_,O61_,O62_]:=aaFOD0[2,0,aALf][NN]+ (6 Pi^2) ( 2 O4/sigmax^2) + 3 (1/2/Pi) (6 Pi^2/sigmax^3) ( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[I 3 fi]) ) (O62 (aALf/(1 + I aALf be0 fi))^(1-ka2) + O61 (aALf/(1 + I aALf be0 fi))^(1-ka1)),{fi,-Pi,Pi}] ); (* this is a^{(2,0)} in the approach with O61=O6ka1 and O62=O6ka2, and a(sigmax Exp[I fi]) in O6ka (a(sigmax Exp[I fi]))^(1-ka) runs by one-loop RGE; note that aALf=a(kappa sigmax); but the higher-twist (D=6) terms are evaluated with the reasonable a(sigmax Exp[I fi])^(1-ka) (and not: a(kappa sigmax Exp[I fi])^(1-ka)) only if kappa=1.   *)

(* 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_]:=aaFOD0[2,1,aALf][NN]+ 2 (1/2/Pi) (6 Pi^2/sigmax^3)*( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 - Exp[2 I fi]) + (1/2) (1 - Exp[I 4 fi]) ) (O62 aCI[fi,aALf]^(1-ka2) + O61 aCI[fi,aALf]^(1-ka1)),{fi,-Pi,Pi}] );
aaFO205lCI[aALf_][NN_,O4_,O61_,O62_]:=aaFOD0[2,0,aALf][NN]+ (6 Pi^2) ( 2 O4/sigmax^2) + 3 (1/2/Pi) (6 Pi^2/sigmax^3) ( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[I 3 fi]) )  (O62 aCI[fi,aALf]^(1-ka2) + O61 aCI[fi,aALf]^(1-ka1)),{fi,-Pi,Pi}] );
  (*-----------------------------_*)
  (* Using the results at Nt=8: *)
  (*  Timing[alO4O612fitPFOMinNum1ls[8]]={81.1516, {0.0031155, {aALf -> 0.104745, O4fitf -> -0.00384881, 61fitf -> -0.0360145, O62fitf -> 0.0185475}}} (alpha_s(mtau^2)=0.317741) OK *)

  (* ---> *)
  (* aaFO21[aALfFO[8]][8, O61FO[8], O62FO[8]]=0.204403;
    aaFO20[aALfFO[8]][8, O4FO[8], O61FO[8], O62FO[8]]=0.239693; OK *)
   VERY GOOD, because:
   aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026;
     (7 Apr. 2022)  
   *)
  (*------------------------*)
  (*  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,aaFOD0[2,0,aALfFO[NN]][NN]},{NN,3,12}]; (* Tab20D0FOavar={{3, 0.223382}, {4, 0.235097}, {5, 0.247083}, {6, 0.259506}, {7, 0.273097}, {8, 0.288454}, {9, 0.30672}, {10, 0.328}, {11, 0.356851}, {12, 0.388502}} *)
 Tab20FOavar=Table[{NN,Re[aaFO20[aALfFO[NN]][NN,O4FO[NN],O61FO[NN],O62FO[NN]]]},{NN,3,12}];(* Tab20FOavar={{3, 0.240714}, {4, 0.24031}, {5, 0.240026}, {6, 0.239842}, {7, 0.239732}, {8, 0.239693}, {9, 0.239721}, {10, 0.239809}, {11, 0.239903}, {12, 0.240015}} *)
 delTab20FOavar= Table[{Nn, Tab20FOavar[[Nn - 2]][[2]] - Tab20FOavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];  (* at Nn=8 is the local minimum; therefore, we take Nt=8 in FO case *)
 Tab21D0FOavar=Table[{NN,aaFOD0[2,1,aALfFO[NN]][NN]},{NN,3,12}]; (* Tab21D0FOavar={{3, 0.20604}, {4, 0.208808}, {5, 0.211483}, {6, 0.213968}, {7, 0.216284}, {8, 0.218346}, {9, 0.220068}, {10, 0.221078}, {11, 0.221569}, {12, 0.220757}} *)
 Tab21FOavar=Table[{NN,Re[aaFO21[aALfFO[NN]][NN,O61FO[NN],O62FO[NN]]]},{NN,3,12}]; (* Tab21FOavar={{3, 0.206923}, {4, 0.205918}, {5, 0.205213}, {6, 0.20476}, {7, 0.204492}, {8, 0.204403}, {9, 0.204468}, {10, 0.204694}, {11, 0.204913}, {12, 0.205207}} *)
delTab21FOavar= Table[{Nn, Tab21FOavar[[Nn - 2]][[2]] - Tab21FOavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];  (* at Nn=8 is the local minimum; 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["a20FOavarO6ka1ka2.pdf",a20FOavar]; Export["a21FOavarO6ka1ka2.pdf",a21FOavar];  Export["a20FOavarO6ka1ka2Z.pdf",a20FOavarZoom]; Export["a21FOavarO6ka1ka2Z.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.104745), [ alpha_s(mtau^2)=alsmtauFO[8]=0.31774 ],i.e., the coupling alpha_s taken FIXED at value as obtaiuned from global BL FO fits at Nt=8: *)
(*
 Tab20D0FOafxd=Table[{NN,aaFOD0[2,0,aALfFO[8]][NN]},{NN,3,12}]; (* Tab20D0FOafxd={{3, 0.207623}, {4, 0.2334}, {5, 0.251455}, {6, 0.265429}, {7, 0.277205}, {8, 0.288454}, {9, 0.300652}, {10, 0.316181}, {11, 0.337358}, {12, 0.370678}}; *)
 Tab21D0FOafxd=Table[{NN,aaFOD0[2,1,aALfFO[8]][NN]},{NN,3,12}]; (* Tab21D0FOafxd={{3, 0.192122}, {4, 0.207419}, {5, 0.214767}, {6, 0.218014}, {7, 0.218795}, {8, 0.218346}, {9, 0.21729}, {10, 0.216578}, {11, 0.215706}, {12, 0.216539}}; *)
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}; 
  (*--------------*)
(*  This is the global tFO approach with (O61, O62) : *)
  (* a21 and a20, this is the case when a(sigmax Exp[I fi]) in [O_6ka (a(sigmax Exp[I fi]))^(1-ka)] runs by one-loop RGE: *)
  (* This is NEW: *)
  aatFO21[aALf_][NN_,O61_,O62_]:=aatFOD0[2,1,aALf][NN]+ 2 (1/2/Pi) (6 Pi^2/sigmax^3)*( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 - Exp[2 I fi]) + (1/2) (1 - Exp[I 4 fi]) ) (O62 (aALf/(1 + I aALf be0 fi))^(1-ka2) + O61 (aALf/(1 + I aALf be0 fi))^(1-ka1)),{fi,-Pi,Pi}] );  (* this is a^{(2,1)} in the approach with O61=O6ka1 and O62=O6ka2, and a(sigmax Exp[I fi]) in O6ka (a(sigmax Exp[I fi]))^(1-ka) runs by one-loop RGE; but the higher-twist (D=6) terms are evaluated with the reasonable a(sigmax Exp[I fi])^(1-ka) (and not: a(kappa sigmax Exp[I fi])^(1-ka)) only if kappa=1.  *)
  aatFO20[aALf_][NN_,O4_,O61_,O62_]:=aatFOD0[2,0,aALf][NN]+ (6 Pi^2) ( 2 O4/sigmax^2) + 3 (1/2/Pi) (6 Pi^2/sigmax^3) ( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[I 3 fi]) ) (O62 (aALf/(1 + I aALf be0 fi))^(1-ka2) + O61 (aALf/(1 + I aALf be0 fi))^(1-ka1)),{fi,-Pi,Pi}] ); (* this is a^{(2,0)} in the approach with O6ka1=O61 and O6ka2-O62, and a(sigmax Exp[I fi]) in O6ka (a(sigmax Exp[I fi]))^(1-ka) runs by one-loop RGE; note that aALf=a(kappa sigmax); but the higher-twist (D=6) terms are evaluated with the reasonable a(sigmax Exp[I fi])^(1-ka) (and not: a(kappa sigmax Exp[I fi])^(1-ka)) only if kappa=1.   *)
  (* end of a21 and a20 in global tFO approach with (O6ka1=O61, O6ka2=O62) [and a(Q^2) in O6ka (a(Q^2))^(1-ka) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ]; note that aALf=a(kappa sigmax)  *)


(* 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.203809;
 aatFO20[aALftFO[5]][5, O4tFO[5], O61tFO[5], O62tFO[5]]=0.239456; OK
    GOOD, because:
   aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
   aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026
 (7 Apr. 2022)
 *)
  (*------------------------*)
(*  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.229353}, {4, 0.250124}, {5, 0.277873}, {6, 0.309182}, {7, 0.356239}, {8, 0.357906}, {9, 0.400065}, {10, 0.337423}, {11, 0.189485}, {12, 0.271974}}; *)
 Tab20tFOavar=Table[{NN,Re[aatFO20[aALftFO[NN]][NN,O4tFO[NN],O61tFO[NN],O62tFO[NN]]]},{NN,3,12}]; (*  Tab20tFOavar={{3, 0.240444}, {4, 0.23984}, {5, 0.239456}, {6, 0.239588}, {7, 0.239893}, {8, 0.24037}, {9, 0.234093}, {10, 0.240865}, {11, 0.224646}, {12, 0.241743}}; *)
 delTab20tFOavar= Table[{Nn, Tab20tFOavar[[Nn - 2]][[2]] - Tab20tFOavar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; (* at about Nn=5 local minimum *)
 Tab21D0tFOavar=Table[{NN,aatFOD0[2,1,aALftFO[NN]][NN]},{NN,3,12}];  (* Tab21D0tFOavar={{3, 0.207503}, {4, 0.212387}, {5, 0.218038}, {6, 0.222407}, {7, 0.225547}, {8, 0.219982}, {9, 0.194356}, {10, 0.210964}, {11, 0.136943}, {12, 0.204527}}; *)
 Tab21tFOavar=Table[{NN,Re[aatFO21[aALftFO[NN]][NN,O61tFO[NN],O62tFO[NN]]]},{NN,3,12}]; (* Tab21tFOavar={{3, 0.206248}, {4, 0.204745}, {5, 0.203809}, {6, 0.204161}, {7, 0.204912}, {8, 0.206072}, {9, 0.19269}, {10, 0.20732}, {11, 0.173095}, {12, 0.209539}}; *)
delTab21tFOavar= Table[{Nn, Tab21tFOavar[[Nn - 2]][[2]] - Tab21tFOavar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; (* at about Nn=5 local minimum *)
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}]; 

 (* We will take Nt=5 in tFO; 17 Apr. 2022 *)
 *)

 (*-------*)   
(*
(*  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.40}}, 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.13, 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}]];

 *)
 (* We will take Nt=5 in tFO; 17 Apr. 2022 *)

 (* Export["a20tFOavarO6ka1ka2.pdf",a20tFOavar]; Export["a21tFOavarO6ka1ka2.pdf",a21tFOavar];  Export["a20tFOavarO6ka1ka2Z.pdf",a20tFOavarZoom]; Export["a21tFOavarO6ka1ka2Z.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.104284), [ alpha_s(mtau^2)=alsmtautFO[5]=0.316404],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.234128}, {4, 0.262265}, {5, 0.277873}, {6, 0.293208}, {7, 0.311759}, {8, 0.350039}, {9, 0.391058}, {10, 0.611905}, {11, 0.51409}, {12, 3.50569}}; *)
 Tab21D0tFOafxd=Table[{NN,aatFOD0[2,1,aALftFO[5]][NN]},{NN,3,12}]; (* Tab21D0tFOafxd={{3, 0.211541}, {4, 0.221219}, {5, 0.218038}, {6, 0.215151}, {7, 0.211583}, {8, 0.217742}, {9, 0.193932}, {10, 0.303434}, {11, -0.187359}, {12, 2.41643}}; *)
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.40, 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=O6ka1, O62=O6ka2, 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: alO4O612fitPFOMinnum1l[8]={0.0031155, {aALf -> 0.104745, O4fitf -> -0.00384881, O61fitf -> -0.0360145, O62fitf -> 0.0185475}}
chi2o2O4[O4f_]:=Sum[(ReBPFONum1l[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.00384881 by global analysis (and O4FOvarfit=-0.00384839, 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[(ReBPFONum1l[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]];  (* =-0.0360153; O61FO[8]=-0.0360145, very close as it should be *)
Xi2O61FOvaro77=O61FOvarfitMin[[1]];
 delO61FO= delO61f/.FindRoot[chi2o2O61[O61FOvarfit+delO61f]==Xi2O61FOvaro77+1,{delO61f, 3. 10^-5}];
(* Nm=2; Nd=6; delO61FO=0.000181; OK *)

chi2o2O62[O62f_]:=Sum[(ReBPFONum1l[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]]; (*=0.0185473; O62FO[8]=0.0185475, very close, OK *)
Xi2O62FOvaro77=O62FOvarfitMin[[1]];
 delO62FO= delO62f/.FindRoot[chi2o2O62[O62FOvarfit+delO62f]==Xi2O62FOvaro77+1,{delO62f, 3. 10^-6}];
(* Nm=2; Nd=6; delO62FO=0.000073; OK, 11 Apr. 2022 *)

chi2o2aAL[aALf_]:=Sum[(ReBPFONum1l[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]]; (* =0.104747; aALfFO[8]=0.104745, very close, OK *)
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; 11 Apr. 2022 *)
(*-------*)
(* 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.000181; delO62FO=0.000074; delaALFO=0.000102(=delta(a(sigmax))) *)
(* --> *)
(*
 (* From here on do: 11 Apr. 2022 *)
 alpha_s(mtau^2)=0.317740+{^{+0.00030}_{-0.00030}}_{exp}; OK 
10^3 O4=-3.849+{^{+0.03}_{-0.03}}_{exp} =app -3.8+-0.0_{exp}; OK
10^2 O61=-3.60+{{+0.02}_{-0.02}}_{exp} =app -3.6+-0.0_{exp}; OK
10^2 O62=+1.85+{^{+0.01}_{-0.01}}_{exp} =app +1.9 +- 0.0_{exp}; OK
10^3 O6eff=2.45+-0.06_{exp} =app 2.4 +-0.1_{exp}; OK *)

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61FO[8], O62=O62FO[8], and a(sigmax)=aALfFO[8]=0.104745 *)
(* O6eff= O61FO[8] aALfFO[8]^(1-ka1) + O62FO[8] aALfFO[8]^(1-ka2) *)
(* delO6effaAL=O61FO[8] ((aALfFO[8] + delaALFO)^(1 - ka1) - aALfFO[8]^(1 - ka1)) + O62FO[8] ((aALfFO[8] + delaALFO)^(1 - ka2) - aALfFO[8]^(1 - ka2)); (* this is variation of O6eff due to variation of aALfFO[8] *) *)
(* delO6eff=Sqrt[ (delO61FO aALfFO[8]^(1-ka1))^2 + (delO62FO aALfFO[8]^(1-ka2))^2 + delO6effaAL^2] (i.e., summing in quadrature) *)

(* 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: *)
(* alO4O612fitPFOMinNum1l[8]={0.00184412, {aALf -> 0.0870077, O4fitf -> -0.00221699, O61fitf -> -0.0270431, O62fitf -> 0.0138048}}; alsmtau=0.316986; aALf(sigmax)=0.104485 *)

(* kappa=1.0 (the central case of the global fit): *)
alO4O612fitPFOMinnum1l[8]={0.0031155, {aALf -> 0.104745, O4fitf -> -0.00384881,  O61fitf -> -0.0360145, O62fitf -> 0.0185475}}; alsmtau=0.31774; aALf(sigmax)=0.104745; *)

(* kappa=2/3: *)
(* alO4O612fitPFOMinNum1l[8]={0.00354399, {aALf -> 0.121097, O4fitf -> -0.00579784, O61fitf -> -0.0422295, O62fitf -> 0.0221366}}; alsmtau=aMSext[mtau^2/((2/3) sigmax), 0.121097] Pi=0.321108; a(sigmax)=aMSext[(3/2.), 0.121097] = 0.105908; *)

(* Note that aALf=a(kappa*sigmax) *)
(* --> thus the variations: (kappa=1+{^{+1}_{-1/3}}) *)
(*    alpha_s(mtau^2)=0.317740+{^{-0.000754}_{+0.003368}}_{kap}; OK 
      10^3 O4= -3.8+{^{+2.6}_{-2.0}}_{kap};  OK
      10^2 O61=-3.6+{^{+0.9}_{-0.6}}_{kap};  OK
      10^2 O62=+1.9+{^{-0.5}_{+0.3}}_{kap}; OK
 10^3 O6eff=2.45+{^{-0.68}_{+0.62}}_{kap} =app 2.4+{^{-0.6}_{+0.7}}_{kap}; OK
*)

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61FO[8], O62=O62FO[8], and a(sigmax)=aALfFO[8]=0.104745 *)
(* O6eff= O61FO[8] aALfFO[8]^(1-ka1) + O62FO[8] aALfFO[8]^(1-ka2) *)
(* O6eff is O6eff when central value kappa=1 is taken *)
(* O6effpl= -0.0270431 0.104485^(1-ka1) + 0.0138048 0.104485^(1-ka2) *)
(* O6effpl = O6eff when kappa=2 taken *)
(* O6effmi= -0.0422295 0.105908^(1-ka1) + 0.0221366 0.105908^(1-ka2) *)
(* O6effmi = O6eff when kappa=2/3 taken *)
 (* 11 Apr. 2022 *)

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

(* when ddtk[4, tkap]=275.+63.=338. (tddtk[4, tkap]=37.582): 
 alO4O612fitPFOMinnum1l[8]={0.0026605, {aALf -> 0.101936, O4fitf -> -0.00331932, O61fitf -> -0.0349174, O62fitf -> 0.0176313}}; alpha_s(mtau^2)=0.309584; OK *)

(* when the central case d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=275. (tddtk[4, tkap]=-25.4181):
alO4O612fitPFOMinnum1l[8]={0.0031155, {aALf -> 0.104745, O4fitf -> -0.00384881, O61fitf -> -0.0360145, O62fitf -> 0.0185475}}; alpha_s(mtau^2)=0.31774; OK *)

(* when ddtk[4, tkap]=ddtk[4,Exp[tKK]]=275.-63.=212.; tddtk[4, tkap]=-88.418:
alO4O612fitPFOMinnum1l[8]={0.00312669, {aALf -> 0.106468, O4fitf -> -0.00412029, O61fitf -> -0.0358183, O62fitf -> 0.0186822}};  alpha_s(mtau^2)=0.32273; *)

(* alpha_s(mtau^2)=0.31774+{^{-0.008156}}_{+0.00499}}_{d4}; OK 
 10^3 O4= -3.8+{^{+0.5}_{-0.3}}_{d4}; OK
 10^2 O61=-3.6+{^{+0.1}_{0.0}}_{d4}; OK
 10^2 O62=+1.9+{^{-0.1}_{+0.0}}_{d4}; OK *) 
 10^3 O6eff=2.45+{^{-0.21}_{+0.08}}_{d4} =app 2.4+{^{-0.2}_{+0.1}}_{d4}; OK
*)

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61FO[8], O62=O62FO[8], and a(sigmax)=aALfFO[8]=0.104745 *)
(* O6eff= O61FO[8] aALfFO[8]^(1-ka1) + O62FO[8] aALfFO[8]^(1-ka2) *)
(* O6eff is O6eff when central value kappa=1 is taken *)
(* O6effpl= -0.0270431 0.104485^(1-ka1) + 0.0138048 0.104485^(1-ka2) *)
(* O6effpl = O6eff when kappa=2 taken *)
(* O6effmi= -0.0422295 0.105908^(1-ka1) + 0.0221366 0.105908^(1-ka2) *)
(* O6effmi = O6eff when kappa=2/3 taken *)
 (* 12 Apr. 2022 *)
 
 (*--------------*)
(* Variation of Nt(truncation order, Nt=8+-1 and Nt=8 +-2, Nt=8-3), FO global: *)
(*
 alO4O612fitPFOMinnum1l[10]={0.00227335, {aALf -> 0.105879, O4fitf -> -0.00654659, O61fitf -> -0.0380327, O62fitf -> 0.0201438}}; alsmtauFO[10]=0.321026; OK 
alO4O612fitPFOMinnum1l[9]={0.00297135, {aALf -> 0.105435, O4fitf -> -0.00512041, O61fitf -> -0.0382952, O62fitf -> 0.0199889}}; alsmtauFO[9]=0.31974; OK
 central case:
alO4O612fitPFOMinnum1l[8]={0.0031155, {aALf -> 0.104745, O4fitf -> -0.00384881, O61fitf -> -0.0360145, O62fitf -> 0.0185475}}; alsmtauFO[8]=0.31774; OK
alO4O612fitPFOMinnum1l[7]={0.00272373, {aALf -> 0.104134, O4fitf -> -0.00274006, O61fitf -> -0.0317815, O62fitf -> 0.016188}}; alsmtauFO[7]=0.315969; OK
alO4O612fitPFOMinnum1l[6]={0.0018433, {aALf -> 0.103745, O4fitf -> -0.00172053, O61fitf -> -0.0256395, O62fitf -> 0.0129543}}; alsmtauFO[6]=0.31484; OK
alO4O612fitPFOMinnum1l[5]={0.000784024, {aALf -> 0.103896, O4fitf -> -0.000754103, O61fitf -> -0.0177607, O62fitf -> 0.00894318}}; alsmtauFO[5]=0.315278; OK
 *)

(*  alpha_s(mtau^2)=0.31774+{^{+0.003285}_{-0.002901}}_{Nt}; OK 
 10^3 O4= -3.8+{^{-2.7}_{+3.0}}_{Nt}; OK
 10^2 O61=-3.6+{^{-0.2}_{+1.8}}_{Nt}; OK
 10^2 O62=+1.9+{^{+0.1}_{-1.0}}_{Nt}; OK
 10^3 O6eff=2.45+{^{+0.44}_{-1.35}}_{Nt} =app 2.4+{^{+0.5}_{-1.3}}_{Nt}; OK
*)

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61FO[8], O62=O62FO[8], and a(sigmax)=aALfFO[8]=0.104745 *)
(* O6eff= O61FO[8] aALfFO[8]^(1-ka1) + O62FO[8] aALfFO[8]^(1-ka2) *)
(* O6eff is O6eff when central value kappa=1 is taken *)
(* O6effpl= -0.0270431 0.104485^(1-ka1) + 0.0138048 0.104485^(1-ka2) *)
(* O6effpl = O6eff when kappa=2 taken *)
(* O6effmi= -0.0422295 0.105908^(1-ka1) + 0.0221366 0.105908^(1-ka2) *)
(* O6effmi = O6eff when kappa=2/3 taken *)
 (* 12 Apr. 2022 *)

 (* 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; Nt=8^{+2}_{-3}):
 alpha_s(mtau^2)=0.317740 + {^{+0.00030}_{-0.00030}}_{exp} + {^{-0.000754}_{+0.003368}}_{kap} + {^{-0.008156}}_{+0.00499}}_{d4} + {^{+0.003285}_{-0.002901}}_{Nt} = 0.317740 +{^{+0.006865}_{-0.008695}} =app 0.3177 +{^{+0.0069}_{-0.0087}} =app 0.318+{^{+0.007}_{-0.009}}; OK
 10^3 O4=-3.8 +- 0.0_{exp} + {^{+2.6}_{-2.0}}_{kap} + {^{+0.5}_{-0.3}}_{d4} + {^{-2.7}_{+3.0}}_{Nt}  = -3.8 + {^{+4.0}_{-3.4}}; OK
 10^2 O61=-3.6 +- 0.0_{exp} + +{^{+0.9}_{-0.6}}_{kap} + {^{+0.1}_{-0.0}}_{d4} + {^{-0.2}_{+1.8}}_{Nt} =-3.6 + {^{+2.0}_{-0.6}}; OK
 10^2 O62=+1.9 +- 0.0_{exp} + {^{-0.5}_{+0.3}}_{kap} + {^{-0.1}_{+0.0}}_{d4} + {^{+0.1}_{-1.0}}_{Nt} =+1.9 + {^{+0.3}_{-1.1}}; OK
10^3 O6eff=+2.4 +- 0.1_{exp} + {^{-0.6}_{+0.7}}_{kap} + {^{-0.2}_{+0.1}}_{d4} + {^{+0.5}_{-1.3}}_{Nt} =+2.4 + {^{+0.9}_{-1.4}}; OK
 OK *) 

(* O6eff= O61FO[Nt] aALfFO[Nt]^(1-ka1) + O62FO[Nt] aALfFO[Nt]^(1-ka2) *)

(* 12 Apr. 2022 *) 

  (* 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: alO4O612fitPtFOMinnum1l[5]={0.00542654, {aALf -> 0.104284, O4fitf -> -0.00318952, O61fitf -> -0.039268, O62fitf -> 0.0199033}}; alsmtautFO[5]=0.316404 (ddtk[4,Exp[tKK]]=275.) OK *)

(* We vary O4 in UNCORRELATED chi ('var'): *)
chi2o2O4[O4f_]:=Sum[(ReBPtFONum1l[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]]; (*=-0.00318902; very close to  O4tFO[5]=-0.00318952; OK, as it should be *)
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; 17 Apr. 2022 *)

chi2o2O61[O61f_]:=Sum[(ReBPtFONum1l[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]]; (* =-0.0392693; very close to O61tFO[5]=-0.039268; OK, as it should be *)
Xi2O61tFOvaro77=O61tFOvarfitMin[[1]];
 delO61tFO= delO61f/.FindRoot[chi2o2O61[O61tFOvarfit+delO61f]==Xi2O61tFOvaro77+1,{delO61f, 3. 10^-5}];
(* Nm=2; Nd=6; delO61tFO=0.0001814; OK *)

chi2o2O62[O62f_]:=Sum[(ReBPtFONum1l[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]]; (* =0.0199029; very close to O62tFO[5]=0.0199033; OK, very close as it should be *)
Xi2O62tFOvaro77=O62tFOvarfitMin[[1]];
delO62tFO= delO62f/.FindRoot[chi2o2O62[O62tFOvarfit+delO62f]==Xi2O62tFOvaro77+1,{delO62f, 3. 10^-5}]; 
(* Nm=2; Nd=6; delO62tFO=0.0000738; OK; 17 Apr. 2022 *)

chi2o2aAL[aALf_]:=Sum[(ReBPtFONum1l[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]]; (*=0.104287; very close to aALftFO[5]=0.104284; OK, as it should be *)
Xi2aALtFOvaro77=aALtFOvarfitMin[[1]];
 delaALtFO= delaALf/.FindRoot[chi2o2aAL[aALtFOvarfit+delaALf]==Xi2aALtFOvaro77+1,{delaALf, 3. 10^-5}];
(* Nm=2; Nd=6; Nt0=6; delaALtFO=0.000117; 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.0001814; delO62tFO=0.000074; delaALtFO=0.000117; (=delta(a(sigmax)))  OK; 17 Apr. 2022;   --> *)

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


(* tFO, Nt=5, exp:
alpha_s(mtau^2)=0.316404 -+0.000339_{exp}; OK 
10^3 O4= -3.19 +-0.03_{exp} =app -3.2 +- 0.0_{exp}; OK
10^2 O61=-3.93+{^{+0.02}_{-0.01}}_{exp} =app -3.9 +-0.0_{exp}; OK
10^2 O62=+1.99+-0.01_{exp} =app +2.0+-0.0; OK
10^3 O6eff=+2.49 +-0.06_{exp} =app +2.5 +-0.1_{exp}; OK
 (* 17 Apr. 2022 *)
 
 (* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[5], O62=O62tFO[5]; *)
 (* O6eff= O61tFO[5] aALftFO[5]^(1-ka1) + O62tFO[5] aALftFO[5]^(1-ka2) *)
 (* delO6effaAL=O61tFO[5] ((aALftFO[5] + delaALtFO)^(1 - ka1) - aALftFO[5]^(1 - ka1)) + O62tFO[5] ((aALftFO[5] + delaALtFO)^(1 - ka2) - aALftFO[5]^(1 - ka2)); (* this is variation of O6eff due to variation of aALftFO[5] *) *)
 (* delO6eff=Sqrt[ (delO61tFO aALftFO[5]^(1-ka1))^2 + (delO62tFO aALftFO[5]^(1-ka2))^2 + delO6effaAL^2] (i.e., summing in quadrature) *)

(* The precise central value result of the tFO global fit:
 alO4O612fitPtFOMinnum[5]={0.00542654, {aALf -> 0.104284, O4fitf -> -0.00318952, O61fitf -> -0.039268, O62fitf -> 0.0199033}}; alsmtautFO[5]=0.316404 (ddtk[4,Exp[tKK]]=275.) OK *)
 (* 17 Apr. 2022 *)
 (*----------*)
 (* Variation of kappa, tFO global (Nt0=5): *)
(* kappa=2.0: 
 alO4O612fitPtFOMinNum1l[5]={0.00128704, {aALf -> 0.0865354, O4fitf -> -0.00111201, O61fitf -> -0.0212528, O62fitf -> 0.0106972}}; alpha_s(mtau^2)=0.314916; a(sigmax)=0.103771; OK
 kappa=1.0 (the central case of the global fit): 
 alO4O612fitPtFOMinnum1l[5]={0.00542654, {aALf -> 0.104284, O4fitf -> -0.00318952, O61fitf -> -0.039268, O62fitf -> 0.0199033}}; alsmtautFO[5]=0.316404; OK
kappa=2/3:
 alO4O612fitPtFOMinNum1l[5]={0.0200679, {aALf -> 0.128504, O4fitf -> -0.00919209, O61fitf -> -0.072382, O62fitf -> 0.0387196}}; alpha_s(mtau^2)=0.336506; a(sigmax)=0.111243; OK *)

 (* --> tFO, Nt=5; thus the variations  (kap=1.+{^{+1}_{-1/3}}: *)
alpha_s(mtau^2)=0.316404+{^{-0.001488}_{+0.020102}}_{kap}; OK
10^3 O4= -3.2+{^{+2.1}_{-6.0}}_{kap}; OK
10^2 O61=-3.9+{^{+1.8}_{-3.3}}_{kap}; OK
10^2 O62=+2.0+{^{-0.9}_{+1.9}}_{kap}; OK
10^3 O6eff=+2.49+{^{-1.18}_{+2.94}}_{kap}  =app +2.5+{^{-1.2}_{+2.9}}_{kap}; OK

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[5], O62=O62tFO[5], and a(sigmax)=aALftFO[5]=0.104745 *)
(* O6eff= O61tFO[5] aALftFO[5]^(1-ka1) + O62tFO[5] aALftFO[5]^(1-ka2) *)
(* O6eff is O6eff when central value kappa=1 is taken *)
 (* O6effpl= -0.0212528 0.103771^(1 - ka1) + 0.0106972  0.103771^(1 - ka2) =0.00131124 *)
(* O6effpl = O6eff when kappa=2 taken *)
 (* O6effmi=-0.072382 0.111243^(1 - ka1) + 0.0387196 0.111243^(1 - ka2) = 0.00543313 *)
(* O6effmi = O6eff when kappa=2/3 taken *)
 (* 17 Apr. 2022 *)
 *)
 (*--------*)
(* 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.=338. (tddtk[4, tkap]=37.582):
alO4O612fitPtFOMinNum1l[5]={0.00423518, {aALf -> 0.103192, O4fitf -> -0.00286698, O61fitf -> -0.0367594, O62fitf -> 0.018518}}; alpha_s(mtau^2)=0.313235; OK
When central d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=275. (tddtk[4, tkap]=-25.4181):
alO4O612fitPtFOMinnum1l[5]={0.00542654, {aALf -> 0.104284, O4fitf -> -0.00318952, O61fitf -> -0.039268, O62fitf -> 0.0199033}}; alsmtautFO[5]=0.316404; OK
 When ddtk[4, tkap]=275-63.=212. (tddtk[4, tkap]=-88.418):
alO4O612fitPtFOMinNum1l[5]={0.00708347, {aALf -> 0.10555, O4fitf -> -0.00358391, O61fitf -> -0.0422263, O62fitf -> 0.0215561}}; alpha_s(mtau^2)=0.320072; OK
  *)
 (* Note that here in this approach we see that td4 does not just have minor variation, it varies td4=-25.4+- 63.0, which is roughly 250% variation around its value *)
 (* ---> *)
 (* tFO, Nt=5, d4=275+-63:
alpha_s(mtau^2)=0.316404+{^{-0.003169}_{+0.003668}}_{d4}; OK
10^3 O4= -3.2+{^{+0.3}_{-0.4}}_{d4}; OK
10^2 O61=-3.9+{^{+0.2}_{-0.3}}_{d4}; OK
10^2 O62=+2.0+{^{-0.1}_{+0.2}}_{d4}; OK *)  
10^3 O6eff=+2.49+{^{-0.20}_{+0.24}}_{d4}  =app +2.5 +{^{-0.2}_{+0.2}}_{d4}
 
(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[5], O62=O62tFO[5], and a(sigmax)=aALftFO[5]=0.104745 *)
(* O6eff= O61tFO[5] aALftFO[5]^(1-ka1) + O62tFO[5] aALftFO[5]^(1-ka2) *)
(* O6eff is O6eff when central value kappa=1 is taken *)
 (* O6effpl = -0.0367594 0.103192^(1 - ka1) + 0.018518 0.103192^(1 - ka2) =0.00229477 *)
(* O6effpl = O6eff when kappa=2 taken *)
 (* O6effmi = -0.0422263 0.10555^(1 - ka1) + 0.0215561 0.10555^(1 - ka2) =0.00272788 *)
(* O6effmi = O6eff when kappa=2/3 taken *)
 (* 17 Apr. 2022 *)
  (*----*)
 (* 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)) *)
(*
 alO4O612fitPtFOMinnum1l[7]={0.00187271, {aALf -> 0.108184, O4fitf -> -0.00857779, O61fitf -> -0.0439567, O62fitf -> 0.0239095}}; alpha_s(mtau^2)=0.32769; OK
 alO4O612fitPtFOMinnum1l[6]={0.00435203, {aALf -> 0.106112, O4fitf -> -0.00541216, O61fitf -> -0.0451519, O62fitf -> 0.0235571}}; alpha_s(mtau^2)=0.3217; OK
alO4O612fitPtFOMinnum1l[5]={0.00542654, {aALf -> 0.104284, O4fitf -> -0.00318952, O61fitf -> -0.039268, O62fitf -> 0.0199033}}; alsmtautFO[5]=0.316404; OK
alO4O612fitPtFOMinnum1l[4]={0.00178082, {aALf -> 0.10223, O4fitf -> -0.00103472, O61fitf -> -0.0228097, O62fitf -> 0.0113147}}; alpha_s(mtau^2)=0.310439; OK
alO4O612fitPtFOMinnum[3]={0.0000269251, {aALf -> 0.103235, O4fitf -> 0.000676531, O61fitf -> -0.00361756, O62fitf -> 0.00181203}}; alpha_s(mtau^2)=0.313359; OK *)
 
  -->
 (* tFO, 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.316404+{^{+0.011287}_{-0.005964}}_{Nt}; OK
10^3 O4= -3.2+{^{-5.4}_{+3.9}}_{Nt}; OK
10^2 O61=-3.9+{^{-0.6}_{+3.5}}_{Nt}; OK
10^2 O62=+2.0+{^{+0.4}_{-1.8}}_{Nt}; OK *) 
10^3 O6eff=+2.49+{^{+1.15}_{-2.27}}_{Nt}  =app +2.5 +{^{+1.1}_{-2.3}}_{Nt}

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[5], O62=O62tFO[5], and a(sigmax)=aALftFO[5]=0.104745 *)
(* O6effN[NNt_] := O61tFO[NNt] aALftFO[NNt]^(1 - ka1) +  O62tFO[NNt] aALftFO[NNt]^(1 - ka2); *)
(* O6eff is O6eff when central value kappa=1 is taken *)
(*-----------------------------------------*)
 (* 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.; Nt=5+-2;
 alpha_s(mtau^2)=0.316404+-0.000339_{exp}+{^{-0.001488}_{+0.020102}}_{kap}+{^{-0.003169}_{+0.003668}}_{d4}+{^{+0.011287}_{-0.005964}}_{Nt} =0.316404+{^{+0.023346}_{-0.006924}}  =0.3164+{^{+0.0233}_{-0.0069}} =app 0.316+{^{+0.024}_{-0.007}};

 10^3 O4=  -3.2+-0.0_{exp}+{^{+2.1}_{-6.0}}_{kap} +{^{+0.3}_{-0.4}}_{d4}+{^{-5.4}_{+3.9}}_{Nt} = -3.2+{^{+4.4}_{-8.1}}; OK 
 10^2 O61= -3.9+-0.0_{exp}+{^{+1.8}_{-3.3}}_{kap}+{^{+0.2}_{-0.3}}_{d4}+{^{-0.6}_{+3.5}}_{Nt} = -3.9+{^{+3.9}_{-3.4}}; OK
 10^2 O62= +2.0+-0.0_{exp}+{^{-0.9}_{+1.9}}_{kap}+{^{-0.1}_{+0.2}}_{d4}+{^{+0.4}_{-1.8}}_{Nt}  = +2.0+{^{+2.0}_{-2.0}}; OK
10^3 O6eff=+2.5+-0.1_{exp}+{^{-1.2}_{+2.9}}_{kap}+{^{-0.2}_{+0.2}}_{d4}+{^{+1.1}_{-2.3}}_{Nt} =+2.5 +{^{+3.1}_{-2.6}}
   (* 17 Apr. 2022 *)
	    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 *)
   

   (* Now the form in D=6 with 1-loop running in D=6 ('1l'); when kappa not=1, RGE from aALf=a(kappa sigmax) to a(sigmax) is performed by 5-loop RGE: aALf=a(kappa sigmax); a(sigmax)=aMSext[(1/kappa),aALf]; (note: aALf=a(kappa sigmax) in this general case): *)
ReBPCINum1l[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 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax ) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O62 (aMSext[(1/kappa), aALf]/(1 + I be0 aMSext[(1/kappa), aALf] fi))^(1-ka2) + O61 (aMSext[(1/kappa), aALf]/(1 + I be0 aMSext[(1/kappa), aALf] fi))^(1-ka1)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
   (* 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.03; O62init= 0.02; aALfinit= 0.107;
alO4O612fitPCIMinNum1l[Nn_]:= FindMinimum[Sum[(ReBPCINum1l[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[(ReBPCINum1l[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[(ReBPCINum1l[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): *)
(*
alO4O612fitPCIMinnum1l[3]=   {0.000020963, {aALf -> 0.1178, O4fitf -> -0.00218144, O61fitf -> -0.0181919, O62fitf -> 0.00895282}};
alO4O612fitPCIMinnum1l[4]={0.0000175616, {aALf -> 0.113585, O4fitf -> -0.00257029, O61fitf -> -0.0223433, O62fitf -> 0.0106357}};
alO4O612fitPCIMinnum1l[5]={9.40354*10^-6, {aALf -> 0.112084, O4fitf -> -0.00300731, O61fitf -> -0.0277844, O62fitf -> 0.0128366}};
alO4O612fitPCIMinnum1l[6]={0.000700066, {aALf -> 0.111551, O4fitf -> -0.00393563, O61fitf -> -0.0425986, O62fitf -> 0.0188973}};
alO4O612fitPCIMinnum1l[7]= {0.00508007, {aALf -> 0.112546, O4fitf -> -0.00500932,  O61fitf -> -0.0639685, O62fitf -> 0.0279025}};
alO4O612fitPCIMinnum1l[8]=   {0.640571, {aALf -> 0.110626, O4fitf -> -0.00535267, O61fitf -> -0.0853226, O62fitf -> 0.0349643}};
alO4O612fitPCIMinnum1l[9]=   {0.987878, {aALf -> 0.109166, O4fitf -> -0.00477, O61fitf -> -0.0764436, O62fitf -> 0.0305832}};
alO4O612fitPCIMinnum1l[10]={5.56455, {aALf -> 0.0985141, O4fitf -> -0.000528123, O61fitf -> 0.00576458, O62fitf -> -0.00747795}};
alO4O612fitPCIMinnum1l[11]= {0.0114756, {aALf -> 0.115049, O4fitf -> -0.00869463, O61fitf -> -0.144495, O62fitf -> 0.0585079}};
alO4O612fitPCIMinnum1l[12]=   {13.6605, {aALf -> 0.0873606, O4fitf -> 0.00382456, O61fitf -> 0.133152, O62fitf -> -0.0601897}};

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

*)
(* end of NEW BL results, CI approach with (O6ka1=O61,O6ka2=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: *)			     
   aaCID0[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) *)
   (*  aaCID0[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) *)
  (*----*)
(* moments a^{(1,n}}, CI: *)
  aaCID0[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) *)
(*---------*)
(* moments a^{(0,0}}, CI: *)				     
aaCID0[0,0,aALf_][Nn_]:=(1/(2 Pi)) NIntegrate[ (1 + Exp[I th])  dff[th,Nn,aALf], {th, -Pi, Pi}];
(*-------------------------*)
(*  This is the NEW (19 Apr.2022) in FO approach with (O6ka1=O61, O6ka2=O62): *)
(* a21 and a20, this is the case when a(sigmax Exp[I fi]) in [O_6^{(ka)} a(sigmax Exp[I fi])^(1-ka)] runs by one-loop RGE: *)
aaCI21[aALf_][NN_,O61_,O62_]:=aaCID0[2,1,aALf][NN]+ 2 (1/2/Pi) (6 Pi^2/sigmax^3)*( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 - Exp[2 I fi]) + (1/2) (1 - Exp[I 4 fi]) ) (O62 (aALf/(1 + I aALf be0 fi))^(1-ka2) + O61 (aALf/(1 + I aALf be0 fi))^(1-ka1)),{fi,-Pi,Pi}] );  (* this is a^{(2,1)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O6ka (a(sigmax Exp[I fi]))^(1-ka)) runs by one-loop RGE; note that aALf=a(kappa sigmax); but the higher-twist (D=6) terms are evaluated with the reasonable a(sigmax Exp[I fi])^(1-ka) (and not: a(kappa sigmax Exp[I fi])^(1-ka)) only if kappa=1.  *)
aaCI20[aALf_][NN_,O4_,O61_,O62_]:=aaCID0[2,0,aALf][NN]+ (6 Pi^2) ( 2 O4/sigmax^2) + 3 (1/2/Pi) (6 Pi^2/sigmax^3) ( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[I 3 fi]) ) (O62 (aALf/(1 + I aALf be0 fi))^(1-ka2) + O61 (aALf/(1 + I aALf be0 fi))^(1-ka1)),{fi,-Pi,Pi}] );    (* this is a^{(2,0)} in the approach with O61=O6ka1 and O62=O6ka2, and a(sigmax Exp[I fi]) in O6ka (a(sigmax Exp[I fi]))^(1-ka) runs by one-loop RGE; note that aALf=a(kappa sigmax); but the higher-twist (D=6) terms are evaluated with the reasonable a(sigmax Exp[I fi])^(1-ka) (and not: a(kappa sigmax Exp[I fi])^(1-ka)) only if kappa=1.   *)

  (*-----------------------------_*)
  (* Using the results at Nt=5: 
   alO4O612fitPCIMinnum1l[5]={9.40354*10^-6, {aALf -> 0.112084, O4fitf -> -0.00300731, O61fitf -> -0.0277844, O62fitf -> 0.0128366}};
     (* ---> *)
   aaCI21[aALfCI[5]][5, O61CI[5], O62CI[5]]=0.206036;  (* OK *)
   aaCI20[aALfCI[5]][5, O4CI[5], O61CI[5], O62CI[5]]=0.240365;  (* OK *)

 Using the results at Nt=4: 
    alO4O612fitPCIMinnum1l[4]={0.0000175616, {aALf -> 0.113585, O4fitf -> -0.00257029, O61fitf -> -0.0223433, O62fitf -> 0.0106357}};
     (* --> *)
   aaCI21[aALfCI[4]][4, O61CI[4], O62CI[4]] =0.205908;
   aaCI20[aALfCI[4]][4, O4CI[4], O61CI[4], O62CI[4]]=0.240312;
      
    VERY GOOD, because:
   aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
   aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026

   19 Apr. 2022
   *)
  (*------------------------*)
(*  This is the NEW (19 Apr. 2022) (O6ka1=O61, O6ka2=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[aaCID0[2,0,aALfCI[NN]][NN]]},{NN,3,12}]; (* Tab20D0CIavar={{3, 0.269357}, {4, 0.274805}, {5, 0.280956}, {6, 0.293803}, {7, 0.308081}, {8, 0.311891}, {9, 0.304181}, {10, 0.249193}, {11, 0.360613}, {12, 0.194329}}; *)
 Tab20CIavar=Table[{NN,Re[aaCI20[aALfCI[NN]][NN,O4CI[NN],O61CI[NN],O62CI[NN]]]},{NN,3,12}]; (* Tab20CIavar={{3, 0.240312}, {4, 0.240312}, {5, 0.240365}, {6, 0.240605}, {7, 0.241023}, {8, 0.23928}, {9, 0.238511}, {10, 0.232988}, {11, 0.241522}, {12, 0.227502}}; *)
 delTab20CIavar= Table[{Nn, Tab20CIavar[[Nn - 2]][[2]] - Tab20CIavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];
 Tab21D0CIavar=Table[{NN,Re[aaCID0[2,1,aALfCI[NN]][NN]]},{NN,3,12}]; (*  Tab21D0CIavar={{3, 0.210709}, {4, 0.210955}, {5, 0.211186}, {6, 0.211998}, {7, 0.213871}, {8, 0.205918}, {9, 0.202175}, {10, 0.174238}, {11, 0.206588}, {12, 0.145512}}; *)
 Tab21CIavar=Table[{NN,Re[aaCI21[aALfCI[NN]][NN,O61CI[NN],O62CI[NN]]]},{NN,3,12}]; (* Tab21CIavar={{3, 0.205912}, {4, 0.205908}, {5, 0.206036}, {6, 0.2066}, {7, 0.207591}, {8, 0.204294}, {9, 0.202681}, {10, 0.191105}, {11, 0.208882}, {12, 0.179533}}; *)
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 *) *)
(* The smalles variation of a20 and a21 is when Nt=3 -> Nt=4; so I will take Nt=4 *)
(* 19 Apr. 2022 *)
 (*-------*)   
(*
 (* 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["a20CIavarO6ka1ka2.pdf",a20CIavar]; Export["a21CIavarO6ka1ka2.pdf",a21CIavar];  Export["a20CIavarO6ka1ka2Z.pdf",a20CIavarZoom]; Export["a21CIavarO6ka1ka2Z.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.113585), [ alpha_s(mtau^2)=alsmtauCI[4]=0.343237],i.e., the coupling alpha_s taken FIXED at value as obtained from global BL CI fits at Nt=4: *)
(*
Tab20D0CIafxd = 
 Table[{NN, Re[aaCID0[2, 0, aALfCI[4]][NN]]}, {NN, 3, 12}]; (* Tab20D0CIafxd={{3, 0.25297}, {4, 0.274805}, {5, 0.288579}, {6, 0.305586}, {7, 0.314585}, {8, 0.332854}, {9, 0.33575}, {10, 0.366548}, {11, 0.348531}, {12, 0.433047}}; *) (* OK *)
 Tab21D0CIafxd=Table[{NN,Re[aaCID0[2,1,aALfCI[4]][NN]]},{NN,3,12}]; (* Tab21D0CIafxd={{3, 0.201471}, {4, 0.210955}, {5, 0.214838}, {6, 0.216791}, {7, 0.216157}, {8, 0.211074}, {9, 0.209588}, {10, 0.18555}, {11, 0.205147}, {12, 0.0831589}}; *) (* 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[(ReBPCINum1l[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]]; (* =-0.0025703; very close to O4CI[4]=-0.00257029; OK, as it should be *)
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.0000303853; 20 Apr. 2022 *)

chi2o2O61[O61f_]:=Sum[(ReBPCINum1l[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]]; (* =-0.0223435; very close to O61CI[4]=-0.0223433; OK, as it should be *)
Xi2O61CIvaro77=O61CIvarfitMin[[1]];
 delO61CI= delO61f/.FindRoot[chi2o2O61[O61CIvarfit+delO61f]==Xi2O61CIvaro77+1,{delO61f, 3. 10^-5}];
(* Nm=2; Nd=6; delO61CI=0.000171764; 20 Apr. 2022 *)

chi2o2O62[O62f_]:=Sum[(ReBPCINum1l[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]]; (* = 0.0106356; very close to O62CI[4]=0.0106357; OK, as it should be *)
Xi2O62CIvaro77=O62CIvarfitMin[[1]];
 delO62CI= delO62f/.FindRoot[chi2o2O62[O62CIvarfit+delO62f]==Xi2O62CIvaro77+1,{delO62f, 3. 10^-5}];
(* Nm=2; Nd=6; delO62CI=0.0000724793; 20 Apr. 2022 *)

chi2o2aAL[aALf_]:=Sum[(ReBPCINum1l[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]]; (*=0.113585; very close to aALfCI[4]=0.113585; OK, as it should be *)
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]=6.81716*10^-6;  chi2o2aAL[aALCIvarfit +-0.000161] =app 1.00 *)
(* Nm=2; Nd=6; Nt0=5; delaALCI=0.000161; 20 Apr. 2022 *)

(* The precise central value result of the CI global fit: 
 alO4O612fitPCIMinnum1l[4]={0.0000175616, {aALf -> 0.113585, O4fitf -> -0.00257029, O61fitf -> -0.0223433, O62fitf -> 0.0106357}}; *)
(*----*)
(* 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.000171764; delO62CI=0.0000724793; delaALCI=0.000161 (=delta(a(sigmax))) --> *)
(* 
alpha_s(mtau^2)=0.343237+{^{+0.000462}_{-0.000462}}_{exp}; 
 10^3 O4= -2.6+-0.0_{exp};
 10^2 O61=-2.2+-0.0_{exp};
 10^2 O62=+1.1+-0.0_{exp};
10^3 O6eff=+0.8+-0.1_{exp};
 OK; 20 Apr. 2022;
 *)
 (* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61CI[4], O62=O62CI[4], and a(sigmax)=aALfCI[4]=0.113585 *)
(* delO6effaAL = O61CI[4] ((aALfCI[4] + delaALCI)^(1 - ka1) - aALfCI[4]^(1 - ka1)) + O62CI[4] ((aALfCI[4] + delaALCI)^(1 - ka2) - aALfCI[4]^(1 - ka2)); *)
 (* delO6eff=Sqrt[(delO61CI aALfCI[4]^(1 - ka1))^2 + (delO62CI aALfCI[4]^(1 - ka2))^2+delO6effaAL^2]; *)
 
(*---------*)

(* Variation of kappa, CI global fit, Nt=4  (note that aALf=a(kappa sigmax) : *)
(*
kappa=2.0: 
 alO4O612fitPCIMinnum1l[4]={0.000184369, {aALf -> 0.0950995, O4fitf -> -0.00170015, O61fitf -> -0.0127195, O62fitf -> 0.00670731}}; alpha_s(mtau^2)=0.353415; a(sigmax)=0.11714; (* 20 Apr.2022 *)
kappa=1.0:
 alO4O612fitPCIMinnum1l[4]={0.0000175616, {aALf -> 0.113585, O4fitf -> -0.00257029, O61fitf -> -0.0223433, O62fitf -> 0.0106357}}; alpha_s(mtau^2)=0.343237; a(sigmax)=0.113585; (* 20 Apr.2022 *)
 kappa=2/3:
 alO4O612fitPCIMinnum1l[4]={0.00429065, {aALf -> 0.132707, O4fitf -> -0.00410865, O61fitf -> -0.0525222, O62fitf -> 0.0236743}}; alpha_s(mtau^2)=0.345007; a(sigmax)=0.114202; (* OK, 20 Apr. 2022 *)
 *)
alpha_s(mtau^2)=0.343237+{^{+0.010178}_{+0.001770}}_{kap}; 
10^3 O4= -2.6+{^{+0.9}_{-1.5}}_{kap};
10^2 O61=-2.2+{^{+0.9}_{-3.1}}_{kap};
10^2 O62=+1.1+{^{-0.4}_{+1.3}}_{kap};
10^3 O6eff=+0.82+{^{+0.02}_{+0.27}}_{kap}  =app +0.8+{^{+0.0}_{+0.3}}_{kap}; OK

 (* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61CI[4], O62=O62CI[4], and a(sigmax)=aALfCI[4]=0.113585 *)
(* O6eff= O61CI[4] aALfCI[4]^(1-ka1) + O62CI[4] aALfCI[4]^(1-ka2) =0.000824665*)
(* O6eff is O6eff when central value kappa=1 is taken *)
 (* O6effpl = -0.0127195 0.11714^(1 - ka1) + 0.00670731 0.11714^(1 - ka2)=0.00083704 *)
(* O6effpl = O6eff when kappa=2 taken *)
 (* O6effmi = -0.0525222 0.114202^(1 - ka1) + 0.0236743 0.114202^(1 - ka2)=0.00109165 *)
(* O6effmi = O6eff when kappa=2/3 taken *)
 (* 20 Apr. 2022 *)
 
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+-1, 4+2), CI global: Nt=6,5,4,3 (we do not go down to Nt=2, an extreme) *)
(* 
alO4O612fitPCIMinnum1l[6]={0.000700066, {aALf -> 0.111551, O4fitf -> -0.00393563, O61fitf -> -0.0425986, O62fitf -> 0.0188973}}; alsmtauCI[6]=0.337394;
alO4O612fitPCIMinnum1l[5]={9.40354*10^-6, {aALf -> 0.112084, O4fitf -> -0.00300731, O61fitf -> -0.0277844, O62fitf -> 0.0128366}}; alsmtauCI[5]=0.338927;
alO4O612fitPCIMinnum1l[4]={0.0000175616, {aALf -> 0.113585, O4fitf -> -0.00257029, O61fitf -> -0.0223433, O62fitf -> 0.0106357}}; alsmtauCI[4]=0.343237;
alO4O612fitPCIMinnum1l[3]={0.000020963, {aALf -> 0.1178, O4fitf -> -0.00218144, O61fitf -> -0.0181919, O62fitf -> 0.00895282}}; alsmtauCI[3]=0.3553;
(* --> *)
(*
 alpha_s(mtau^2)=0.343237+{^{-0.00584}_{+0.012062}}_{Nt}; 
10^3 O4= -2.6+{^{-1.3}_{+0.4}}_{Nt};
10^2 O61=-2.2+{^{-2.1}_{+0.4}}_{Nt};
10^2 O62=+1.1+{^{+0.8}_{-0.2}}_{Nt};
10^3 O6eff=+0.82+{^{-0.02}_{-0.03}}_{Nt}  =app +0.8+{^{-0.0}_{+0.0}}_{Nt}; OK

 (* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61CI[4], O62=O62CI[4], and a(sigmax)=aALfCI[4]=0.113585 *)
  (* O6effN[NNt_]:= O61CI[NNt] aALfCI[NNt]^(1-ka1) + O62CI[NNt] aALfCI[NNt]^(1-ka2) *)
 (* 20 Apr. 2022 *)
 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+-1, 4+2;
 alpha_s(mtau^2)=0.343237+{^{+0.000462}_{-0.000462}}_{exp}+{^{+0.010178}_{+0.001770}}_{kap}+-0.0_{d4}+{^{-0.00584}_{+0.012062}}_{Nt} =0.343237+{^{+0.015789}_{-0.005858}} =app 0.3432+{^{+0.0158}_{-0.0058}} =app 0.343+{^{+0.016}_{-0.006}};
 10^3 O4= -2.6 +-0.0_{exp}+{^{+0.9}_{-1.5}}_{kap}+{^{-1.3}_{+0.4}}_{Nt} =-2.6+{^{+1.0}_{-2.0}}; 
 10^2 O61=-2.2 +-0.0_{exp}+{^{+0.9}_{-3.1}}_{kap}+-0.0_{d4}+{^{-2.1}_{+0.4}}_{Nt} =-2.2+{^{+1.0}_{-3.7}};
10^2 O62=+1.1 +-0.0_{exp}+{^{-0.4}_{+1.3}}_{kap}+-0.0_{d4}+{^{+0.8}_{-0.2}}_{Nt} =+1.1+{^{+1.5}_{-0.4}};  
10^3 O6eff=+0.8 +-0.1_{exp}+{^{+0.0}_{+0.3}}_{kap}+-0.0_{d4}+-0.0_{Nt} =+0.8+{^{+0.3}_{-0.1}};  
 (* OK *)
 (* 20 Apr. 2022 *)
 
(* 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.10536;
(* 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): *)    
ReBPPVNum1l[M2f_,Nn_,O4_,O61_,O62_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPPVD0[M2f,Nn]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax ) + (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O62 (aMSext[(1/kappa), aAL]/(1 + I be0 aMSext[(1/kappa), aAL] fi))^(1-ka2) + O61 (aMSext[(1/kappa), aAL]/(1 + I be0 aMSext[(1/kappa), aAL] fi))^(1-ka1)) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
 (* Re[B[M2,Nn,O4,O61,O62]] for PV with delD truncated at a^Nn, RScl kappa; a(\sigmax Exp[I fi]) in D=6 part runs by one-loop RGE, and has the initial value a(sigmax)=aMSext[(1/kappa), aAL] [note that aAL=a(kappa sigmax)] *)
    (* aAL= a(kappa sigmax) *)

 (* Now search for condensate values with this Borel-Laplace (double) pinched ('P') approach, with (Borel-resummed) PV method for D=0 part: *)
NchP=10; (* number of intervals in: (M2min < |M^2} < 1.5) *)
M2FP[ii_]:=M2min+ii*(1.5-M2min)/NchP; (* ii=0,1,2,...,NchP *)
(* We will use sum of squares minimization in this interval, for M^2=|M^2| Exp[I Psi] with Psi=0, Pi/6, Pi/4 one after another in the sum (it turns out that it is not good to include Psi=Pi/3, this value appears to be too large, and xi^2 become large): *)

(* 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[(ReBPPVNum1l[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[(ReBPPVNum1l[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[(ReBPPVNum1l[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: *)

 (* (* 7 Apr. 2022 *):
  (* see TabPVO6ka1O67ka2.save for finding to 3 digits aAL for the min(xi^2)=O4O612fitPPVMinN[Nn][[1]]: *)

  (*  aALfPV[3]=0.1061(xi2=0.00708); 0.1062(xi2=0.00458); 0.1063(xi2=0.00300); 0.1064(xi2=0.002286); 0.10641(xi2=0.002269); 0.10642(xi2=0.002258);  0.10643(xi2=0.002256)***; 0.10644(xi2=0.002264); 0.10645(xi2=0.002281); 0.1065(xi2=0.00250); 0.1066(xi2=0.00363); --> *)
  aALfPV[3]=0.10643;  O4O612fitPPVMinnum[3]={0.00225644, {O4fitf -> -0.00252459, O61fitf -> -0.0212409, O62fitf -> 0.0126795}};

 (*  aALfPV[4]=0.1050(xi2=0.00826); 0.1052(xi2=0.00369); 0.1053(xi2=0.00248); 0.1054(xi2=0.002022); 0.10541(xi2=0.002017)***; 0.10542(xi2=0.002019); 0.10543(xi2=0.002029); 0.10544(xi2=0.002047); 0.1055(xi2=0.00231); 0.1056(xi2=0.00335); --> *)
  aALfPV[4]=0.10541; O4O612fitPPVMinnum[4]={0.00201664, {O4fitf -> -0.00265903, O61fitf -> -0.0230276, O62fitf -> 0.0132392}}; 
   
  (*  aALfPV[5]=0.1051(xi2=0.00452); 0.1052(xi2=0.00298); 0.1053(xi2=0.00214); 0.10535(xi2=0.002002); 0.10536(xi2=0.001997)***; 0.10537(xi2=0.001999); 0.10538(xi2=0.002009); 0.1054(xi2=0.002052); 0.1055(xi2=0.00272); 0.1056(xi2=0.00413); --> *) 
  aALfPV[5]=0.10536; O4O612fitPPVMinnum[5]={0.00199687, {O4fitf -> -0.00267484, O61fitf -> -0.023215, O62fitf -> 0.0133054}}; 
    
   (*  aALfPV[6]=0.1051(xi2=0.00357); 0.1052(xi2=0.00239); 0.1053(xi2=0.001929); 0.10531(xi2=0.001924)***; 0.10532(xi2=0.001926); 0.10533(xi2=0.001936); 0.10534(xi2=0.001953); 0.1054(xi2=0.00221); 0.1055(xi2=0.00325); --> *)
  aALfPV[6]=0.10531; O4O612fitPPVMinnum[6]={0.00192396, {O4fitf -> -0.00271981, O61fitf -> -0.0238468, O62fitf -> 0.0135461}};
  
  (*  aALfPV[7]=0.1051(xi2=0.00360); 0.1052(xi2=0.00242); 0.1053(xi2=0.001969); 0.10531(xi2=0.001965)***; 0.10532(xi2=0.001968); 0.10533(xi2=0.001978); 0.10534(xi2=0.001997); 0.1054(xi2=0.00226); 0.1055(xi2=0.00331); --> *)
  aALfPV[7]=0.10531;  O4O612fitPPVMinnum[7]={0.00196458, {O4fitf -> -0.00270326, O61fitf -> -0.0235453, O62fitf -> 0.0134267}};
  
  (* aALfPV[8]=0.1045(xi2=0.02801); 0.1050(xi2=0.00621); 0.1051(xi2=0.00405);  0.1052(xi2=0.00263); 0.1053(xi2=0.001908); 0.10533(xi2=0.001837); 0.10534(xi2=0.001828); 0.10535(xi2=0.001826)****; 0.10536(xi2=0.001832); 0.1054(xi2=0.001929); 0.1055(xi2=0.00269); 0.1056(xi2=0.00419); --> *)
  aALfPV[8]=0.10535;  O4O612fitPPVMinnum[8]={0.00182587, {O4fitf -> -0.00274744, O61fitf -> -0.0245586, O62fitf -> 0.0138481}};

  (* aALfPV[9]=0.1050(xi2=0.00483); 0.1051(xi2=0.00317); 0.1052(xi2=0.00226); 0.10525(xi2=0.002077); 0.10526(xi2=0.002063); 0.10527(xi2=0.002057)***; 0.10528(xi2=0.002058); 0.10529(xi2=0.002067); 0.1053(xi2=0.002084); 0.1054(xi2=0.00267); --> *)
 aALfPV[9]=0.10527; O4O612fitPPVMinnum[9]={0.00205668, {O4fitf -> -0.00268958, O61fitf -> -0.0230107, O62fitf -> 0.0131881}};   

  (* aALfPV[10]= 0.1048(xi2=0.0160); 0.1050(xi2=0.00837); 0.1052(xi2=0.00361); 0.1053(xi2=0.00226); 0.1054(xi2=0.001629); 0.10542(xi2=0.001589); 0.10543(xi2=0.001580);  0.10544(xi2=0.001578)***;  0.10545(xi2=0.001583); 0.10546(xi2=0.001596); 0.10548(xi2=0.001642); 0.1055(xi2=0.00172); 0.1056(xi2=0.00252); --> *)
  aALfPV[10]=0.10544; O4O612fitPPVMinnum[10]={0.00157761, {O4fitf -> -0.00279294, O61fitf -> -0.0260485, O62fitf -> 0.0145066}};   

  (* aALfPV[11]= 0.1051(xi2=0.002574); 0.10512(xi2=0.002536); 0.10513(xi2=0.002529); 0.10514(xi2=0.002514)***; 0.10515(xi2=0.002521); 0.10516(xi2=0.002536); 0.1052(xi2=0.002675); 0.1053(xi2=0.00358); 0.1054(xi2=0.00527); --> *)
  aALfPV[11]=0.10514; O4O612fitPPVMinnum[11]={0.00251446, {O4fitf -> -0.00262481, O61fitf -> -0.0208682, O62fitf -> 0.0122296}};  

   (* aALfPV[12]= 0.1052(xi2=0.00451); 0.1053(xi2=0.00277); 0.1054(xi2=0.00173); 0.10548(xi2=0.001402); 0.10549(xi2=0.001393); 0.10550(xi2=0.001391)***; 0.10551(xi2=0.001396); 0.10552(xi2=0.001408); 0.1056(xi2=.001756); --> *)
   aALfPV[12]=0.10550; O4O612fitPPVMinnum[12]={0.00139098, {O4fitf -> -0.00281688, O61fitf -> -0.0267973, O62fitf -> 0.0148722}};

 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}];
  (* 7 Apr. 2022 *)
  *)      

(* end of NEW BL results (O61,O62) approach in PPV *)
  (*---------------------------------------------------------------------_*)
  (*---------------------------------------------------------------------_*)
 (* Note, for Nt=4, ...12, the extracted values of aALfPV[Nt] are quite stable, almost independent of Nt (and the condensate values also); this is so probably because Table[{jj, delddtk[jj, tkap]/ddtk[jj, tkap]}, {jj, 1, 12}] gives small values for jj >= 4 (< 0.1, usually smaller than 0.05), which means that delddtk[jj,tkap] are well behaved and have also in practice no renormalon increase as jj increases (theoretically it shoulld be so at least at large jj) *)
 (*---------------------------------------------------------------------_*)
  (*---------------------------------------------------------------------_*)
 (* 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 (8 Apr.2022) (O6ka1, O6ka2) 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]+ 2 (1/2/Pi) (6 Pi^2/sigmax^3)*( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 - Exp[2 I fi]) + (1/2) (1 - Exp[I 4 fi]) ) (O62 (aAL/(1 + I aAL be0 fi))^(1-ka2) + O61 (aAL/(1 + I aAL be0 fi))^(1-ka1)),{fi,-Pi,Pi}] );   (* this is a^{(2,1)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O6ka (a(sigmax Exp[I fi]))^(1-ka)) runs by one-loop RGE; note that aALf=a(kappa sigmax); but the higher-twist (D=6) terms are evaluated with the reasonable a(sigmax Exp[I fi])^(1-ka) (and not: a(kappa sigmax Exp[I fi])^(1-ka)) only if kappa=1.  *)
 aaPV20[Nt_,O4_,O61_,O62_]:=aaPVD0fast[2,0][Nt]+ (6 Pi^2) ( 2 O4/sigmax^2) + 3 (1/2/Pi) (6 Pi^2/sigmax^3) ( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[I 3 fi]) ) (O62 (aAL/(1 + I aAL be0 fi))^(1-ka2) + O61 (aAL/(1 + I aAL be0 fi))^(1-ka1)),{fi,-Pi,Pi}] );  (* this is a^{(2,0)} in the approach with O61 and O62, and a(sigmax Exp[I fi]) in O6ka (a(sigmax Exp[I fi]))^(1-ka)) runs by one-loop RGE; note that aALf=a(kappa sigmax); but the higher-twist (D=6) terms are evaluated with the reasonable a(sigmax Exp[I fi])^(1-ka) (and not: a(kappa sigmax Exp[I fi])^(1-ka)) only if kappa=1.  *)

 (* end of a21 and a20 in global P:V approach with (O6ka1=O61, O6ka2=O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ] *)
 
 (* Note that aaPV21 and aaPV20 are valid only when kappa=1, i.e., when aAL=a(sigmax) *)
 (*----------*)
  (*-----------------------------*)
  (* Using the results at Nt=6: *)
 (*  aALfPV[6]=0.10531; O4O612fitPPVMinnum[6]={0.00192396, {O4fitf -> -0.00271981, O61fitf -> -0.0238468, O62fitf -> 0.0135461}}; 
        (* ---> *)
   (* aaPV21[6, O61PV[6], O62PV[6]]=0.204782; OK
    aaPV20[6, O4PV[6], O61PV[6], O62PV[6]]=0.239847; OK *)
    GOOD, because:
   aexpBin77[2, 1]+-delaexpBin77[2, 1]=0.2061+-0.0026;
   aexpBin77[2, 0]+-delaexpBin77[2, 0]= 0.2404+-0.0026
  (* 8 Apr. 2022 *)
  *)
  (*------------------------*)
(*  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[{NNt,Re[aaPVD0fast[2,0][NNt]]},{NNt,3,12}]*)={{3,0.269833},{4,0.271951},{5,0.272188},{6,0.272838},{7,0.272608},{8,0.273183},{9,0.272473},{10,0.273678},{11,0.271788},{12,0.273901}}; (* OK *)
 Tab20PVavar(*=Table[{NNt,Re[aaPV20[NNt,O4PV[NNt],O61PV[NNt],O62PV[NNt]]]},{NNt,3,12}]*)={{3,0.23981},{4,0.239837},{5,0.239839},{6,0.239847},{7,0.239843},{8,0.239862},{9,0.239831},{10,0.239893},{11,0.239787},{12,0.239917}};  (* OK *)
 delTab20PVavar= Table[{Nn, Tab20PVavar[[Nn - 2]][[2]] - Tab20PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; 

 Tab21D0PVavar(*=Table[{NNt,Re[aaPVD0fast[2,1][NNt]]},{NNt,3,12}]; *)={{3,0.218473},{4,0.218117},{5,0.218099},{6,0.21809},{7,0.218082},{8,0.218185},{9,0.21797},{10,0.218467},{11,0.217521},{12,0.21874}};  (* OK *)
 Tab21PVavar(*=Table[{NNt,Re[aaPV21[NNt,O61PV[NNt],O62PV[NNt]]]},{NNt,3,12}];*)={{3,0.204692},{4,0.204757},{5,0.204762},{6,0.204782},{7,0.204772},{8,0.204817},{9,0.204744},{10,0.204888},{11,0.20464},{12,0.204944}};  (* OK *)
delTab21PVavar= Table[{Nn, Tab21PVavar[[Nn - 2]][[2]] - Tab21PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];

 ff20exp=Table[{NNt,aexpBin77[2,0]},{NNt,3,12}]; ff20plexp=Table[{NNt,aexpBin77[2,0]+delaexpBin77[2,0]},{NNt,3,12}]; ff20miexp=Table[{NNt,aexpBin77[2,0]-delaexpBin77[2,0]},{NNt,3,12}];
ff21exp=Table[{NNt,aexpBin77[2,1]},{NNt,3,12}];  ff21plexp=Table[{NNt,aexpBin77[2,1]+delaexpBin77[2,1]},{NNt,3,12}]; ff21miexp=Table[{NNt,aexpBin77[2,1]-delaexpBin77[2,1]},{NNt,3,12}]; 

  (* NNt=Nt = 5 seems to be the truncation index value of best stationarity of a20 and a21 *) (* 8 Apr. 2022 *)
 (* OK *)
 *)
 (*-------*)   
(*
 (* 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.28}}, 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.75}]], 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.65}]]}, 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.65}]], 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.65}]]}, PlotMarkers -> Automatic, PlotLegends -> Placed[{"PV (D=0)", "PV"}, {0.20, 0.65}]];  (* 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=5 (aAL=aALfPV[5]=0.10536) local stability 
 (* Export["a20PVavarO6ka1ka2.pdf",a20PVavar]; Export["a21PVavarO6ka1ka2.pdf",a21PVavar];  Export["a20PVavarO6ka1ka2Z.pdf",a20PVavarZoom]; Export["a21PVavarO6ka1ka2Z.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[5](=0.10536), [ alpha_s(mtau^2)=alsmtauPV[5]=0.319523],i.e., the coupling alpha_s taken FIXED at value as obtained from global BL PV fits at Nt=5: *)
 (*
  (* Note, now we do not adjust for each NN to aAL=aALfPV[NN]; now aAL=0.10536 (=aALfPV[5]) and we use throughout DBorfast for that aAL, at each NN: *)
  Tab20D0PVafxd(*=Table[{NN,Re[aaPVD0fast[2,0][NN]]},{NN,3,12}]*)={{3, 0.26438}, {4, 0.271684}, {5, 0.272188}, {6, 0.273108}, {7, 0.272879}, {8, 0.273237}, {9, 0.272959}, {10, 0.273244}, {11, 0.272977}, {12, 0.273141}}; (* OK *)
Tab21D0PVafxd(*=Table[{NN, Re[aaPVD0fast[2, 1][NN]]}, {NN, 3, 12}]*)={{3, 0.214335}, {4, 0.217921}, {5, 0.218099}, {6, 0.218286}, {7, 0.218278}, {8, 0.218224}, {9, 0.218323}, {10, 0.218155}, {11, 0.218384}, {12, 0.218193}}; (* 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.28}}, 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.40, 0.60}]], Inset[Style[Superscript[(a), "(2,0)"], 13], Scaled[{0.08, 0.95}]], Inset[Style[Subscript[\[Alpha], s][Subsuperscript[m, \[Tau], 2]] == Round[10000 alsmtauPV[8]]/10000., 13], Scaled[{0.75, 0.40}]], 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.60}]]; (* OK *)

a21PVafxd= ListPlot[{Tab21D0PVafxd, ff21exp, ff21miexp, ff21plexp}, PlotRange -> {{2.8, 12.2}, {0.20, 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.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=5. (note, we have here the case: ddtk[4,Exp[tKK]=275.) OK *)

 (* the central value of the PV global fit result, at Nt=5, is:
  aAL=aALfPV[5]=0.10536; O4O612fitPPVMinnum[5]={0.00199687, {O4fitf -> -0.00267484, O61fitf -> -0.023215, O62fitf -> 0.0133054}};  alpha_s(mtau^2)(=alsmtauPV[5])=0.319523; 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=5;
 (* Use the Borel transform functions, etc., generated with aAL=aALfPV[Nt0] *)
 (* We vary O4 in UNCORRELATED chi ('var'): *)
chi2o2O4[O4f_]:=Sum[(ReBPPVNum1l[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.00267456; very close to ,O4PV[Nt0]=-0.00267484; OK, as it should be *)
 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[(ReBPPVNum1l[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.0232158, very close to O61PV[Nt0]=-0.023215, as it should be, OK *)
 delO61PV= delO61f/.FindRoot[chi2o2O61[O61PVvarfit+delO61f]==Xi2O61PVvaro77+1,{delO61f, 3. 10^-5}];
(* Nm=2; Nd=6; delO61PV=0.0001802; OK *)

chi2o2O62[O62f_]:=Sum[(ReBPPVNum1l[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.0133052, very close to O62PV[Nt0]=0.0133054, as it should be, OK *)
 delO62PV= delO62f/.FindRoot[chi2o2O62[O62PVvarfit+delO62f]==Xi2O62PVvaro77+1,{delO62f, 3. 10^-5}];
 (* Nm=2; Nd=6; delO62PV=0.0000737; OK *)

 aAL=0.10536; (*= aALfPV[Nt0], Nt0=5, OK *)
 chi2o2aAL=Sum[(ReBPPVNum1l[M2covF[jj],Nt0,O4PV[Nt0],O61PV[Nt0],O62PV[Nt0]] - ReBPexpBinAL[M2covF[jj]])^2/ delReBPexpBinAL[M2covF[jj]]^2,{jj,1,Nd}]; (* -> chi2o2aAL=0.00101406; OK *)
 aAL=0.10536 + 0.000182: (* go back and run from the beginning of the PV part: chi2o2aAL=1.00278; *)
aAL=0.10536 - 0.000175: chi2o2aAL=0.999833;   (* OK *)
 (* --> delaALPV= ^{+0.000182}_{-0.000175}; *)
(* Nm=2; Nd=6; Nt0=5; delaALPV=+0.000182; -0.000175; OK *)
(*-------*)
(* Together, the estimates of EXPERIMENTAL uncertainties in the PV global approach (with Nt=5) are (Nm=2; Nd=6), in PV global:  *)
 (* delO4PV= 0.0000303853; delO61PV=0.0001802; delO62PV=0.0000737; delaALPV=0.000182, -0.0001751; (=delta(a(sigmax)))  OK --> *)

(*
 alpha_s(mtau^2)=0.319523+{^{+0.000522}_{-0.000507}}_{exp};
10^3 O4= -2.7+-0.0_{exp};
10^2 O61=-2.3 +-0.0_{exp};
10^2 O62=+1.3 +-0.0_{exp};
10^3 O6eff=+2.4 +-0.1{exp};
OK *)
 (* 17 Apr. 2022 *)
 
 (* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[Nt0], O62=O62tFO[Nt0]; *)
 (* O6eff= O61PV[Nt0] aALfPV[Nt0]^(1-ka1) + O62PV[Nt0] aALfPV[Nt0]^(1-ka2) *)
 (* delO6effaAL=O61PV[Nt0] ((aALfPV[Nt0] + delaALPV)^(1 - ka1) - aALfPV[Nt0]^(1 - ka1)) + O62PV[Nt0] ((aALfPV[Nt0] + delaALPV)^(1 - ka2) - aALfPV[Nt0]^(1 - ka2)); (* this is variation of O6eff due to variation of aALfPV[Nt0] *) *)
 (* delO6eff=Sqrt[ (delO61PV aALfPV[Nt0]^(1-ka1))^2 + (delO62PV aALfPV[Nt0(1-ka2))^2 + delO6effaAL^2] (i.e., summing in quadrature) *)
 
 (* 17 Apr. 2022 *)
(* The precise central value result of the PV global fit:
 Nt0=5; (aAL=)aALfPV[Nt0]=0.10536; O4O612fitPPVMinnum[Nt0]={0.00199687, {O4fitf -> -0.00267484, O61fitf -> -0.023215, O62fitf -> 0.0133054}}; OK *){0.00310712, {O4fitf -> -0.00455454, O61fitf -> -0.00673939, O62fitf -> 0.000892787}}; OK *)
 
(*--------------*)
(* Variation of kappa, PV global: *)
(* kappa=2.0: 
 aAL=0.087285;  O4O612fitPPVMinN[5]={0.00175094, {O4fitf -> -0.00298668, O61fitf -> -0.0263085, O62fitf -> 0.0143661}}; alpha_s(mtau^2)=0.318204; a(sigmax)=0.104905; OK
kappa=1.0 (the central case of the global fit): 
aAL=0.10536; O4O612fitPPVMinnum[5]={0.00199687, {O4fitf -> -0.00267484, O61fitf -> -0.023215, O62fitf -> 0.0133054}}; alpha_s(mtau^2)=0.319523; OK
kappa=2/3:
 aAL=0.12024; O4O612fitPPVMinN[5]={0.00243557, {O4fitf -> -0.00251762, O61fitf -> -0.0203363, O62fitf -> 0.0121764}}; alpha_s(mtau^2)=0.319292; a(sigmax)=0.105281; OK
---> *)

 (*
alpha_s(mtau^2)=0.319523+{^{-0.001319}_{-0.000231}}_{kap};
10^3 O4=  -2.7+{^{-0.3}_{+0.2}_{kap};
10^2 O61= -2.3+{^{-0.3}_{+0.3}}_{kap};
10^2 O62= +1.3+{^{+0.1}_{-0.1}}_{kap};
10^3 O6eff=+2.4+{^{-0.1}_{+0.0}}_{kap};

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[Nt0], O62=O62tFO[Nt0]; *)
 (* O6eff= O61PV[Nt0] aALfPV[Nt0]^(1-ka1) + O62PV[Nt0] aALfPV[Nt0]^(1-ka2) *)
(* O6effpl= -0.0263085 0.104905^(1 - ka1) + 0.0143661 0.104905^(1 - ka2) =0.00228048; *)
(* O6effmi= -0.0203363 0.105281^(1 - ka1) + 0.0121764 0.105281^(1 - ka2)=0.00240708; *)
(* 17 Apr. 2022 *)		
(*----------------*)
 (* Variation of d4, PV global (d4=227 +- 63; kappa=1): 
  When ddtk[4,tkap]=227+63=338 (delddtk[4, tkap]=24.3209):
  aAL=0.10238; O4O612fitPPVMinnum[5]={0.0025934, {O4fitf -> -0.00233239, O61fitf -> -0.0530774, O62fitf -> 0.0275458}}; alpha_s(mtau^2)=0.310875; (* OK *) 
  Central case: ddtk[4,tkap]=275. (delddtk[4, tkap]=14.1972):
 aAL=aALfPV[5]=0.10536; O4O612fitPPVMinnum[5]={0.00199687, {O4fitf -> -0.00267484, O61fitf -> -0.023215, O62fitf -> 0.0133054}}; alpha_s(mtau^2)=0.319523; (* OK *)
 When ddtk[4,tkap]=227-63=212 (delddtk[4, tkap]=18.4374):
  aAL=0.10658; O4O612fitPPVMinnum[5]={0.00154571, {O4fitf -> -0.00258177, O61fitf -> -0.024215, O62fitf -> 0.0134492}}; alpha_s(mtau^2)=0.323054;  (* OK *)
 (*
alpha_s(mtau^2)=0.319523+{^{-0.008648}_{+0.003531}}_{d4};
10^3 O4= -2.7+{^{+0.4}_{+0.1}}_{d4};
10^2 O61=-2.3+{^{-3.0}_{-0.1}}_{d4};
10^2 O62=+1.3+{^{+1.5}_{+0.0}}_{d4};
10^3 O6eff=+2.4+{^{+1.4}_{-0.2}}_{d4};

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[Nt0], O62=O62tFO[Nt0]; *)
 (* O6eff= O61PV[Nt0] aALfPV[Nt0]^(1-ka1) + O62PV[Nt0] aALfPV[Nt0]^(1-ka2) *)
  (* O6effpl= -0.0530774 0.10238^(1 - ka1) + 0.0275458 0.10238^(1 - ka2)=0.00383519; *)
  (* O6effmi= -0.024215 0.10658^(1 - ka1) + 0.0134492 0.10658^(1 - ka2)=0.00220604; *)
  (* 17 Apr. 2022 *)
  (* Here I finished on 17 Apr. 2022 *)
  (*---------------------*)
  (* Variation of Nt (truncation order, Nt=Nt0+2-2, i.e., variation Nt=7,6,5,4,3; PV global with (O61, O62); *)
   aALfPV[7]=0.10531;  O4O612fitPPVMinnum[7]={0.00196458, {O4fitf -> -0.00270326, O61fitf -> -0.0235453, O62fitf -> 0.0134267}}; alpha_s(mtau^2)=0.319378;
   aALfPV[6]=0.10531; O4O612fitPPVMinnum[6]={0.00192396, {O4fitf -> -0.00271981, O61fitf -> -0.0238468, O62fitf -> 0.0135461}}; alpha_s(mtau^2)=0.319378; 
  Central:
   aALfPV[5]=0.10536; O4O612fitPPVMinnum[5]={0.00199687, {O4fitf -> -0.00267484, O61fitf -> -0.023215, O62fitf -> 0.0133054}}; alpha_s(mtau^2)=0.319523;
   aALfPV[4]=0.10541; O4O612fitPPVMinnum[4]={0.00201664, {O4fitf -> -0.00265903, O61fitf -> -0.0230276, O62fitf -> 0.0132392}}; alpha_s(mtau^2)=0.319667;   
   aALfPV[3]=0.10643;  O4O612fitPPVMinnum[3]={0.00225644, {O4fitf -> -0.00252459, O61fitf -> -0.0212409, O62fitf -> 0.0126795}}; alpha_s(mtau^2)=0.32262;
*) 

 (*
alpha_s(mtau^2)=0.319523+{^{-0.000145}_{+0.003098}}_{Nt};
10^3 O4=  -2.7+{^{-0.0}_{+0.1}}_{Nt};
10^2 O61= -2.3+{^{-0.1}_{+0.2}}_{Nt};
10^2 O62= +1.3+{^{+0.1}_{-0.1}}_{Nt};
10^3 O6eff=+2.4+{^{+0.0}_{+0.1}}_{Nt};

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[Nt0], O62=O62tFO[Nt0]; *)
  (* O6effN[NNt_]:= O61PV[NNt] aALfPV[NNt]^(1-ka1) + O62PV[NNt] aALfPV[NNt]^(1-ka2); *)
 (* O6effpl= -0.0530774 0.10238^(1 - ka1) + 0.0275458 0.10238^(1 - ka2)=0.00383519; *)
 (* O6effmi= -0.024215 0.10658^(1 - ka1) + 0.0134492 0.10658^(1 - ka2)=0.00220604; *)

  Note: Nt=5^{+2}_{-2}, 5^{+1}_{-1}
  OK  *)
  (* 17 Apr. 2022 *)
  (* I cam to here on 17 April 2022; bthe parts below still to be done *)
  (*-------------------------*)


  (* Variation due to the renormalon ambiguity ('amb') to Adler: in ReBPPVNum1l we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn] + BPPVdelD0[M2f]) *)
(* Nt=5; O4O6O8fitPPVMinN[5] gives the following xi^2 =xi2[5]: *)
(*
 aAL=0.1060:  xi2[5]=0.0211749;
 aAL=0.10536: xi2[5]=0.00251943;
 aAL=0.1053:  xi2[5]=0.00222612;
 aAL=0.10529: xi2[5]=0.00220413;
 aAL=0.10528: xi2[5]=0.0021898; 
 aAL=0.10527: xi2[5]=0.00218314;  *** 
 aAL=0.10526: xi2[5]=0.00218415;
 aAL=0.10525: xi2[5]=0.00219281;
 aAL=0.10524: xi2[5]=0.00220913;
 aAL=0.10523: xi2[5]=0.00223311;
 aAL=0.1052:  xi2[5]=0.0023509;
 aAL=0.1050:  xi2[5]=0.00485702;
...
*)
 (* Therefore: aAL=0.10527 where xi2[5] is minimal: *)
 (* BBPVD0[M2f,Nn] -> (BPPVD0[M2f, Nn] + BPPVdelD0[M2f]); Nt=5:
  aAL=0.10527; O4O612fitPPVMinN[5]={0.00218314, {O4fitf -> -0.00346542, O61fitf -> -0.0221627, O62fitf -> 0.0128103}}; alpha_s(mtau^2)=0.319262; OK; 18 Apr. 2022  *)
(*---*)
 (* The central value  (BBPVD0[M2f,Nn]): 
  aAL=0.10536;  O4O6O8fitPPVMinN[5]={0.00199687, {O4fitf -> -0.00267484, O61fitf -> -0.023215, O62fitf -> 0.0133054}}; alpha_s(mtau^2)=0.319523; OK; 18 Apr. 2022 *)
  (*---*)
(* When  in ReBPPVNum1l we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn] - BPPVdelD0[M2f]) *)
(* Nt=5; O4O6O8fitPPVMinN[5] gives the following xi^2 =xi2[5]: *)
(*
 aAL=0.10551: xi2[5]=0.0019088;
 aAL=0.10550: xi2[5]=0.00187526;
 aAL=0.10549: xi2[5]=0.00184909;
 aAL=0.10548: xi2[5]=0.00183028;
 aAL=0.10547: xi2[5]=0.00181883;
 aAL=0.10546: xi2[5]=0.00181473; *** 
 aAL=0.10545: xi2[5]=0.00181798;
 aAL=0.1054:  xi2[5]=0.00194421;
 aAL=0.10536: xi2[5]=0.00217742;
  *)
 (* Therefore: aAL=0.10546 where xi2[5] is minimal: *)
(* BBPVD0[M2f,Nn] -> (BPPVD0[M2f, Nn] - BPPVdelD0[M2f]); Nt=5:
 aAL=0.10546;  O4O6O8fitPPVMinN[5]={0.00181473, {O4fitf -> -0.0018788, O61fitf -> -0.0243616, O62fitf -> 0.0138468}}; alpha_s(mtau^2)=0.319812; OK; 18 Apr. 2022  *)
  (*---*)
   (*------*)
(* --> aAL=0.10536 + {^{-0.00009}_{+0.0001}}_{amb} OK  --> *)     
(*
alpha_s(mtau^2)=0.319523+{^{-0.000261}_{+0.000290}}_{amb};
10^3 O4=  -2.7+{^{-0.8}_{+0.8}}_{amb};
 10^2 O61= -2.3+{^{+0.1}_{-0.1}}_{amb};
10^2 O62= +1.3+{^{-0.0}_{+0.1}}_{amb};
10^3 O6eff=+2.4+{^{-0.1}_{+0.0}}_{amb};

(* O6eff = O61 a(sigmax)^(1-ka1) + O62 a(sigmax)^(1-ka2); where O61=O61tFO[Nt0], O62=O62tFO[Nt0]; *)
  (* O6effN[NNt_]:= O61PV[NNt] aALfPV[NNt]^(1-ka1) + O62PV[NNt] aALfPV[NNt]^(1-ka2); *)
(* O6effpl=-0.0221627 0.10527^(1 - ka1) + 0.0128103 0.10527^(1 - ka2)=0.00234457; *)
(* O6effmi=-0.0243616 0.10546^(1 - ka1) + 0.0138468 0.10546^(1 - ka2)= 0.00243157; *)
(* 18 Apr. 2022 *)
 (*-----------------------------------------*)
 (* 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)=aALfPV[Nt] to a(Q^2=sigmax Exp[I fi]) ], gives: *)
  (* PV(+CI)
   Nt=5;  ( ddtk[4, Exp[tKK]]=275.);
   d4=275.+-63.19; kap=1^{+1}_{-1/3};  Nt=5+-1, 5+-2:;
   alpha_s(mtau^2)=0.319523+{^{+0.000522}_{-0.000507}}_{exp}+{^{-0.001319}_{-0.000231}}_{kap}+{^{-0.008648}_{+0.003531}}_{d4}+{^{-0.000145}_{+0.003098}}_{Nt}+{^{-0.000261}_{+0.000290}}_{amb}  =0.319523+{^{+0.004735}_{-0.008768}} =0.3195+{^{+0.0048}_{-0.0087}} = 0.320+{^{+0.004}_{-0.009}};
10^3 O4=    -2.7+-0.0_{exp}+{^{-0.3}_{+0.2}_{kap}+{^{+0.4}_{+0.1}}_{d4}+{^{-0.0}_{+0.1}}_{Nt}+{^{-0.8}_{+0.8}}_{amb}  = -2.7+{^{+0.9}_{-0.9};
10^2 O61=   -2.3+-0.0_{exp}+{^{-0.3}_{+0.3}}_{kap}+{^{-3.0}_{-0.1}}_{d4}+{^{-0.1}_{+0.2}}_{Nt}+{^{+0.1}_{-0.1}}_{amb}=  -2.3+{^{+0.4}_{-3.0}};
10^2 O62=   +1.3+-0.0_{exp}+-0.1_{kap}+{^{+1.5}_{+0.0}}_{d4}+-0.1_{Nt}+{^{-0.0}_{+0.1}}_{amb} =+1.3+{^{+1.5}_{-0.1}};
10^3 O6eff =+2.4+-0.1_{exp}+{^{-0.1}_{+0.0}}_{kap}+{^{+1.4}_{-0.2}}_{d4}+{^{+0.0}_{+0.1}}_{Nt}+{^{-0.1}_{+0.0}}_{amb} =+2.4+{^{+1.4}_{-0.3}};

(* OK; 18 Apr. 2022 *)
			  
 (* End of PV(+CI) global fit results with (O61, O62) [and a(Q^2) in O62/a(Q^2) runs according to one-loop RGE from a(sigmax)=aALftFO[Nt] to a(Q^2=sigmax Exp[I fi]) ] *)
(*------------------------------------*)
(*------------------------------------*)
(*------------------------------------*)
 (* 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) *)
