(* SumRD0ka12O6ka2CombV.m: as SumRMSbCombVp3ka1ka2O6ka2.m, but a more COMPACT, shorter version *)
(* The experimental data based on a combination ('Comb') of ALEPH and OPAL and ee->hadron production data, for the V-channel; cf. PRD103,034028(2012) (see their Table I on page 15; Boito, Golterman, Maltman, Peris); and the correlation matrix obtained directly from these authors by request. *)
(*  D=0 part of Adler function, d(Q^2)_{D=0}, is generated by the 'p3ka1ka2' form of Borel of the auxiliary quantity td(Q^2), i.e., the Borel transform B[td](u) has, among other terms (u=2, u=-1 renormalon terms): two u=3 terms, ~1/(3-u)^ka1 and ~1/(3-u)^ka2, with: ka1=0.375; ka2=0.778; this is as suggested by the work of Boito, Hornung and Jamin (BHJ), 1510.03812v2; OPE is taken with zero anom. dimensions, except for D=6 term which has the form ~ a(Q^2)^k2 <O6>/(Q^2)^3 with k2=1-ka2 =0.222; we take d4=(275 +- 63) for the coefficient at a^5=(alpha_s(Q^2)/Pi)^5 in d(Q^2)_{D=0}; the renormalisation scale parameter kap=mu^2/Q^2 varies in the range kap=1^{+1}_{-1/2}; experimentral error of extracted parameters takes into account the correlation matrix (68x68) of the data. *)
(*--------------------------------------------------*)
NchP=2; (* number of intervals in: (M2min < |M^2} < 1.5) *)
M2min=0.90;
M2FP[ii_]:=M2min+ii*(1.5-M2min)/NchP; (* ii=0,1,2,...,NchP *)
(* We will use sum of squares minimization in this interval, for M^2=|M^2| Exp[I Psi] with Psi=0, Pi/6, Pi/4 one after another in the sum (it turns out that it is not good to include Psi=Pi/3, this value appears to be too large, and xi^2 become large): *)
Nm=NchP+1; Nd=3 Nm;  
M2al0[ii_]:=0.90+(ii-1)*0.60/(Nm-1);
Do[M2al[jj]=M2al0[jj],{jj,1,Nm}]; Do[M2al[jj+Nm]=M2al0[jj] Exp[I Pi/6],{jj,1,Nm}]; Do[M2al[jj+2 Nm]=M2al0[jj] Exp[I Pi/4],{jj,1,Nm}];
 (*---------*)
 (* ka1=0.375; *) ka2=0.778;  (* April 2022 *)
 k3=ka2;
 (* --> 1-ka2=0.222, this is the power in the OPE at the condensate O6; see JHEP 2015 (1510.03812v2) of Boito, Hornung and Jamin, eq. (2.31) *)
 
 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275 (ECH-predicted): *)
<</home/gorazd/SumRpQCD.dir/MSbarRenMod5Ap3ka1ka2d4275.save;
                           
 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275.+63.=338.: *)
 (* <</home/gorazd/SumRpQCD.dir/MSbarRenMod5Ap3ka1ka2d4338.save; when d4barPREDICTED=338. (RenMod-predicted in MiniMOM): *)
 
 (* when d4barPREDICTED(=ddtk[4, Exp[tKK]])=275.-63.=212.: *)
 (*  <</home/gorazd/SumRpQCD.dir/MSbarRenMod5Ap3ka1ka2d4212.save; *)

  (*---------*)
mtau=1.77686;  (*= 1.77686+-0.00012 *)
(*--------------------------------__*)
 (* Combined ('Comb') experimental data: *)
 <</home/gorazd/SumRpQCD.dir/V_CombPRD103;
 (* NNAL=68 is the number of bins *)
NNAL=68;
 sigmax=TabsjComb[[NNAL]];  (* =3.0574 GeV^2 *) (* note, mtau^2=3.15723 GeV^2 *)
 Do[sbin[ii]= TabsjComb[[ii]], {ii, 1, NNAL}];
 Do[OmegaV[ii]= TabOmegaVComb[[ii]], {ii, 1, NNAL}]; 
 sbin[0]=0; OmegaV[0]=0;
 (* These are the original s_i and OmegaV(s_i) *)
 (* Now the (effective) delta(s_j)=dsbin[j] that appear in integrals of Integrate[f(sig) OmegaV(sig), {sig,0, sigmax}] = Sum[f[sbin[k]] OmegaV[sbin[k]] dsbin[k],{k,1,NNAL}]: *) 
 Clear[dsbin];
 Do[dsbin[ii]=(1/2)(sbin[ii+1]-sbin[ii-1]), {ii, 1, NNAL-1}];
 dsbin[NNAL]=(1/2)(sbin[NNAL]-sbin[NNAL-1]); 
 (* Note, Integrate[dsig f(sig) OmegaV(sig),{s,0,sigmax}] -> Sum_{k=1}^{NNAL} (sbin[k]-sbin[k-1]) (1/2)( f(sbin[k]) OmegaV[k] + f(sbin[k-1]) OmegaV[k-1]); according to instruction of Peris, email of 22 Aug. 2022 *)
 (* This gives then simply: =Sum{k=1}^{NNAL} f(sbin[k]) OmegaV[k] dsbin[k] *)
 (* i.e., Integrate[dsig f(sig) OmegaV(sig),{s,0,sigmax}] -> Sum{k=1}^{NNAL} f(sbin[k]) OmegaV[k] dsbin[k] *)
 Do[covVV[ii,jj]=TabCovOmegaVComb[[(ii-1)*NNAL+jj]],{ii,1,NNAL},{jj,1,NNAL}];

 (*--------------------------------------------*)
 (* Now for FOPT: *)

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

betaMSbar[x_]:=-be0 x^2 - be1 x^3 - bett2MS x^4 - bett3MS x^5 - bett4MS[3.] x^6;
bta[x_]:=betaMSbar[x];
(* Note: x =a(mu^2) = alpha_s(mu^2)/Pi *)
(* x=Log[Q2/sigmax];  *) delext=5.; Clear[x];
solQ2pos[aALf_]:= NDSolve[{y'[x]==bta[y[x]], y[0]==aALf}, y, {x,-0.9,+1.5}];
(* ' *)
(*  aMSext[xx_,aAL0f_? NumberQ]:=First[y[Log[xx]]/.solQ2pos[aAL0f]]; *)
 aMSext[xx_,aAL0f_]:=First[y[Log[xx]]/.solQ2pos[aAL0f]];
			     (* aMSext[x,aAL0]=a(Q^2=x Q20), where Q2=x*Q20 and a(Q^2=Q20) = aAL0 *)
(*----------------------------------------------*)
<</home/gorazd/SumRpQCD.dir/Prgs.dir/SumRthMSbar.save;
  (* This gives coeffext coefficients:  coeffext[jj,nn] is the coefficient of a^nn in ta[jj]={tilde a}_{jj} (nn <= 16) *)
  (*-----------------------------------_*)  
  (* (Double-)Pinched ('P') Borel-Laplace in FO approach (PFO): *)

  kappa=1;  (* kappa = mu^2/Q^2, where mu^2 is RScl *)
  tkap=kappa Exp[tKK];
  
(* Now (double) pinched ('P') Re[BPFOfast]: *)
cKP[kp_,M2p_]:=(1/(2 Pi)) (-I)^kp NIntegrate[ th^kp ( ( (1 + Exp[I th])^2 - 2 (M2p/sigmax) (1 + Exp[I th]) + 2  (M2p/sigmax)^2 ) Exp[(sigmax/M2p) Exp[I th]] - 2 (M2p/sigmax)^2 Exp[-(sigmax/M2p)] ),{th, -Pi,Pi}];
  coeffBPFO[M2f_,Nn_]:=Sum[tddtk[nn,tkap] Binomial[nn + kk, kk] be0^kk coeffext[nn + kk + 1, Nn] cKP[kk,M2f], {nn, 0, Nn - 1}, {kk, 0, Nn - 1 - nn}]; 
 (* coeffBPFO[M2,Nn] is the coefficient at x^Nn in the expansion of BPFOexpD0[x,M2,Nn] or of any BPFOexpD0[x,M2,MM,M2] with MM >= Nn *)
 BPFOD0[M2_,Nn_,aALf_]:=Sum[coeffBPFO[M2,nn] x^nn,{nn,1,Nn}]/.{x->aALf};
(* aALf=alpha_s(kappa*sigmax)/Pi *)
  (* -------------------------------------*)	
(* experimental values of double-pinched Borel-Laplace: *)
 (* Note that, in comparison to the V+A case, we replace OmegaVplA[ii] by: 2 OmegaV[ii]; in this was, the D=0 part of OPE is unchanged, the canonical Adler function *)
 BPexpBinComb[M2_]:=(1/M2) Sum[ (1 - sbin[ii]/sigmax)^2 Exp[-sbin[ii]/M2]*2*OmegaV[ii]*dsbin[ii],{ii,1,NNAL}];
 ReBPexpBinComb[M2_]:=Re[(1/M2) Sum[ (1 - sbin[ii]/sigmax)^2 Exp[-sbin[ii]/M2]*2*OmegaV[ii]*dsbin[ii],{ii,1,NNAL}]];
delReBPexpBinComb[M2_]:= Sqrt[4 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],{ii,1,NNAL},{jj,1,NNAL}]]; 
(* end of experimental values *)

(* Covariant-matrix values of the Re(BPexp): *) 
 CovReBPexpBinComb[M2A_,M2B_]:=4 Sum[Re[(1 - sbin[ii]/sigmax)^2 Exp[-sbin[ii]/M2A]/M2A] Re[(1 - sbin[jj]/sigmax)^2 Exp[-sbin[jj]/M2B]/M2B] dsbin[ii] dsbin[jj] covVV[ii,jj],{ii,1,NNAL},{jj,1,NNAL}];
 (* Note, factor 4 because we have 2 OmegaV[ii] in ReBPexpBinComb[M2_], i.e., factor 2 *)
(* CovReBPexpBinComb[M2,M2]=delReBPexpBinComb[M2]^2 *)
  (*-----------*) 
(* Finding of Re[B] in FO approach; *)

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

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

ReBPO4[M2f_,O4_]:=Re[2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )];
 ReBPO6[M2f_,O6_]:=Re[(6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aALf/(1 + I aALf be0 fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
 ReBPO8[M2f_,O8_]:=Re[2 Pi^2 O8  (1/6) (1/M2f^4) (1+ 6 (M2f/sigmax) + 6 (M2f/sigmax)^2)];
 ReBPO10[M2f_,O10_]:=Re[2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2)];
 ReBPO12[M2f_,O12_]:=Re[ 2 Pi^2 O12  (1/5!) (1/M2f^6) (1+ 2*5 (M2f/sigmax) + 5*4 (M2f/sigmax)^2)];
 ReBPO14[M2f_,O14_]:=Re[2 Pi^2 O14  (1/6!) (1/M2f^7) (1+ 2*6 (M2f/sigmax) + 6*5 (M2f/sigmax)^2)];
 ReBPO16[M2f_,O16_]:=Re[2 Pi^2 O16  (1/7!) (1/M2f^8) (1+ 2*7 (M2f/sigmax) + 7*6 (M2f/sigmax)^2)];
ReBPO18[M2f_,O18_]:=Re[2 Pi^2 O18  (1/8!) (1/M2f^9) (1+ 2*8 (M2f/sigmax) + 8*7 (M2f/sigmax)^2)];
(*-------------------*)
(* Now search for condensate values with this Borel-Laplace (double) pinched ('P') approach, with FO method for D=0 part: *)
 
Clear[aALf];
aALfinit=0.105; O4init=-0.003; O6init=-0.001; 
 O8init=0.001; O6init = 0.003; O10init=0.005; O12init=0.004; O14init=0.005; O16init=0.004; O18init=0.005;
(* FO approach: *)
 alO8fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO8[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}]; 
(* Now O8=0: *)
 alO6fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO8[M2al[jj], Nn, aALf,O4fitf,O6fitf,0] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init}}];
			     
(* Now additional condensates (D>8) with zero anom. dimension ('an0'): *)
alO10fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO10[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf,O10fitf] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init}}];
alO12fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO12[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init}}];
 alO14fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO14[M2al[jj],Nn,aALf,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf,O14fitf] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init},{O14fitf, O14init}}];
 alO14fitPFOMinNum1lalt[Nn_]:= FindMinimum[Sum[(ReBPFONumO16[M2al[jj],Nn,aALf,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf,O14fitf,0] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init},{O14fitf, O14init}}]; (* gives the same as alO14fitPFOMinNum[Nn], I checked *)
 alO16fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO16[M2al[jj],Nn,aALf,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf,O14fitf,O16fitf] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init},{O14fitf, O14init},{O16fitf, O16init}}];
 alO18fitPFOMinNum[Nn_]:= FindMinimum[Sum[(ReBPFONumO18[M2al[jj],Nn,aALf,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf,O14fitf,O16fitf,O18fitf] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init},{O14fitf, O14init},{O16fitf, O16init},{O18fitf, O18init}}];
 

 (* Now the simplified form ('1ls') with 1-loop RGE running in D=6, which is valid only for kappa=1: *)
alO8fitPFOMinNum1ls[Nn_]:= FindMinimum[Sum[(ReBPFONum1ls[M2al[jj], Nn, aALf,O4fitf,O6fitf,O8fitf] - ReBPexpBinAL[M2al[jj]])^2/delReBPexpBinAL[M2al[jj]]^2, {jj, 1, Nd}], {{aALf,aALfinit},{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init}}];

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

(* First OPE up to O8: *)
(* Do[alO8fitPFOMinnum[Nn] = alO8fitPFOMinNum[Nn], {Nn, 3, 12}]: *)
 alO8fitPFOMinnum[3]={0.0372992, {aALf -> 0.118814, O4fitf -> -0.0012288, 
  O6fitf -> -0.00474407, O8fitf -> 0.0025538}};
alO8fitPFOMinnum[4]={0.0279993, {aALf -> 0.112934, O4fitf -> -0.00245831, 
  O6fitf -> -0.00324364, O8fitf -> 0.00204874}};
 alO8fitPFOMinnum[5]={0.0217684, {aALf -> 0.110865, O4fitf -> -0.00375331, 
  O6fitf -> -0.00180631, O8fitf -> 0.00161335}};
alO8fitPFOMinnum[6]={0.0177579, {aALf -> 0.110293, O4fitf -> -0.00515511, 
  O6fitf -> -0.000464651, O8fitf -> 0.00123969}};
alO8fitPFOMinnum[7]={0.0151768, {aALf -> 0.11054, O4fitf -> -0.00678522, 
  O6fitf -> 0.00082451, O8fitf -> 0.000910491}};
alO8fitPFOMinnum[8]={0.0138528, {aALf -> 0.111239, O4fitf -> -0.00878507, 
  O6fitf -> 0.00203821, O8fitf -> 0.000636454}};
alO8fitPFOMinnum[9]={0.0136393, {aALf -> 0.112249, O4fitf -> -0.0114423, 
  O6fitf -> 0.0031634, O8fitf -> 0.000428706}};
alO8fitPFOMinnum[10]={0.0150586, {aALf -> 0.113127, O4fitf -> -0.0149524, 
  O6fitf -> 0.0039334, O8fitf -> 0.00036895}};
alO8fitPFOMinnum[11]={0.0171913, {aALf -> 0.114114, O4fitf -> -0.0204556, 
  O6fitf -> 0.00449757, O8fitf -> 0.000414598}};
alO8fitPFOMinnum[12]={0.0213671, {aALf -> 0.113497, O4fitf -> -0.0262468, 
  O6fitf -> 0.00378263, O8fitf -> 0.000780983}};
  Clear[Nn]; Do[Module[{}, xi2FOO8[Nn]=alO8fitPFOMinnum[Nn][[1]]; aALfFOO8[Nn]=aALf/.alO8fitPFOMinnum[Nn][[2]]; alsmtauFOO8[Nn]=aMSext[mtau^2/(kappa sigmax),aALfFOO8[Nn]] Pi;  O4FOO8[Nn]=O4fitf/.alO8fitPFOMinnum[Nn][[2]]; O6FOO8[Nn]=O6fitf/.alO8fitPFOMinnum[Nn][[2]]; O8FOO8[Nn]=O8fitf/.alO8fitPFOMinnum[Nn][[2]]],{Nn,3,12}]; (* OK *)
 (* 20 Sept. 2022 *)
 (*-----*)  
  (* When including O10: *)
alO10fitPFOMinnum[3]={0.000468416, {aALf -> 0.108554, O4fitf -> 0.0023653, O6fitf -> -0.0165085, O8fitf -> 0.0120746, O10fitf -> -0.00619667}};
alO10fitPFOMinnum[4]={0.000363131, {aALf -> 0.104901, O4fitf -> 0.00108062, O6fitf -> -0.0142531, O8fitf -> 0.010715, O10fitf -> -0.0055669}};
alO10fitPFOMinnum[5]={0.000296743, {aALf -> 0.103993, O4fitf -> -0.000174432, O6fitf -> -0.0122475, O8fitf -> 0.00963993, O10fitf -> -0.00508076}};
alO10fitPFOMinnum[6]={0.000258065, {aALf -> 0.104045, O4fitf -> -0.00141047, O6fitf -> -0.0105707, O8fitf -> 0.00882783, O10fitf -> -0.00473109}};
alO10fitPFOMinnum[7]={0.000237254, {aALf -> 0.104557, O4fitf -> -0.00269582, O6fitf -> -0.00917748, O8fitf -> 0.00822993, O10fitf -> -0.00449205}};
alO10fitPFOMinnum[8]={0.000232008, {aALf -> 0.10521, O4fitf -> -0.00406041, O6fitf -> -0.00813787, O8fitf -> 0.00787249, O10fitf -> -0.00437269}};
alO10fitPFOMinnum[9]={0.000239035, {aALf -> 0.105862, O4fitf -> -0.00558015, O6fitf -> -0.00748477, O8fitf -> 0.00775786, O10fitf -> -0.0043674}};
alO10fitPFOMinnum[10]={0.000258906, {aALf -> 0.10616, O4fitf -> -0.00719903, O6fitf -> -0.00742232, O8fitf -> 0.0079557, O10fitf -> -0.00449688}};
alO10fitPFOMinnum[11]={0.000276271, {aALf -> 0.106391, O4fitf -> -0.00936452, O6fitf -> -0.00758875, O8fitf -> 0.0082454, O10fitf -> -0.00464896}};
alO10fitPFOMinnum[12]={0.000299255, {aALf -> 0.105764, O4fitf -> -0.0115836, O6fitf -> -0.00841774, O8fitf -> 0.00876595, O10fitf -> -0.00487739}};
  
 Clear[Nn]; Do[Module[{}, xi2FOO10[Nn]=alO10fitPFOMinnum[Nn][[1]]; aALfFOO10[Nn]=aALf/.alO10fitPFOMinnum[Nn][[2]]; alsmtauFOO10[Nn]=aMSext[mtau^2/(kappa sigmax),aALfFOO10[Nn]] Pi;  O4FOO10[Nn]=O4fitf/.alO10fitPFOMinnum[Nn][[2]]; O6FOO10[Nn]=O6fitf/.alO10fitPFOMinnum[Nn][[2]]; O8FOO10[Nn]=O8fitf/.alO10fitPFOMinnum[Nn][[2]]]; O10FOO10[Nn]=O10fitf/.alO10fitPFOMinnum[Nn][[2]],{Nn,3,12}]; (* OK *)
 (* 20 Sept. 2022 *)
 (*------*)
(* When including O12: *)
alO12fitPFOMinnum[3]={3.0098373736638357`*^-6, {aALf -> 0.10425427189742154`, 
  O4fitf -> 0.00418040954845473`, O6fitf -> -0.02438604094557204`, 
  O8fitf -> 0.02100184930748344`, O10fitf -> -0.01773835385294055`, 
  O12fitf -> 0.00856269983874853`}};
alO12fitPFOMinnum[4]={2.1529443700413555`*^-6, {aALf -> 0.10141412359999684`, 
  O4fitf -> 0.0028654921382851957`, O6fitf -> -0.02161623218020555`, 
  O8fitf -> 0.018866698883443524`, O10fitf -> -0.01597455051865921`, 
  O12fitf -> 0.007663684779880596`}};
alO12fitPFOMinnum[5]={1.6647123145421099`*^-6, {aALf -> 0.10089593943075899`, 
  O4fitf -> 0.0016373835534777526`, O6fitf -> -0.019279603147665973`, 
  O8fitf -> 0.017269329206383365`, O10fitf -> -0.014696254514302594`, 
  O12fitf -> 0.007029606071590681`}};
alO12fitPFOMinnum[6]={1.4151084101495323`*^-6, {aALf -> 0.10111807041751376`, 
  O4fitf -> 0.0004948314894119508`, O6fitf -> -0.01746205859028248`, 
  O8fitf -> 0.01615964576306968`, O10fitf -> -0.013859596995574226`, 
  O12fitf -> 0.006631556579731124`}};
alO12fitPFOMinnum[7]={1.3063525829312999`*^-6, {aALf -> 0.10165139303784408`, 
  O4fitf -> -0.0006161592345017561`, O6fitf -> -0.0160892463953796`, 
  O8fitf -> 0.015437587730085741`, O10fitf -> -0.013365431872314955`, 
  O12fitf -> 0.006411991600158035`}};
alO12fitPFOMinnum[8]={1.3090077454738742`*^-6, {aALf -> 0.1022092839761624`, 
  O4fitf -> -0.0017008605723491004`, O6fitf -> -0.015217865999921027`,
   O8fitf -> 0.015111574818264691`, O10fitf -> -0.01320644737408753`, 
  O12fitf -> 0.006361776637205605`}};
alO12fitPFOMinnum[9]={1.3857888138244514`*^-6, {aALf -> 0.10268010677140371`, 
  O4fitf -> -0.002804517427768928`, O6fitf -> -0.014824256141643266`, 
  O8fitf -> 0.015128481321918448`, O10fitf -> -0.013317283514439957`, 
  O12fitf -> 0.006440925470563028`}};
alO12fitPFOMinnum[10]={1.524123013580342`*^-6, {aALf -> 0.10280567816532336`, 
  O4fitf -> -0.003884116189909939`, O6fitf -> -0.015004509865655458`, 
  O8fitf -> 0.015495398028431638`, O10fitf -> -0.0136792063694428`, 
  O12fitf -> 0.006633365443849604`}};
alO12fitPFOMinnum[11]={1.6140656443903257`*^-6, {aALf -> 0.1028913471133843`, 
  O4fitf -> -0.005293536749927028`, O6fitf -> -0.015299417624083994`, 
  O8fitf -> 0.01587032350050324`, O10fitf -> -0.013995345341799562`, 
  O12fitf -> 0.0067858771626714805`}};
alO12fitPFOMinnum[12]={1.7229522135570993`*^-6, {aALf -> 0.10237297808734626`,  O4fitf -> -0.006770224121266429`, O6fitf -> -0.016047711762919917`, 
  O8fitf -> 0.01641800422160808`, O10fitf -> -0.014402672471641806`, 
  O12fitf -> 0.006972443135782378`}};

 Clear[Nn]; Do[Module[{}, xi2FOO12[Nn]=alO12fitPFOMinnum[Nn][[1]]; aALfFOO12[Nn]=aALf/.alO12fitPFOMinnum[Nn][[2]]; alsmtauFOO12[Nn]=aMSext[mtau^2/(kappa sigmax),aALfFOO12[Nn]] Pi;  O4FOO12[Nn]=O4fitf/.alO12fitPFOMinnum[Nn][[2]]; O6FOO12[Nn]=O6fitf/.alO12fitPFOMinnum[Nn][[2]]; O8FOO12[Nn]=O8fitf/.alO12fitPFOMinnum[Nn][[2]]]; O10FOO12[Nn]=O10fitf/.alO12fitPFOMinnum[Nn][[2]]; O12FOO12[Nn]=O12fitf/.alO12fitPFOMinnum[Nn][[2]],{Nn,3,12}]; (* OK *)
 (* 20 Sept. 2022 *)
 (*-------*)
(* When including O14: *)
alO14fitPFOMinnum[3]={1.1383111947888302`*^-10, {aALf -> 0.10166154974500308`,
  O4fitf -> 0.005315827138057051`, O6fitf -> -0.029798030291659058`, 
  O8fitf -> 0.027898302944348984`, O10fitf -> -0.02919156228567578`, 
  O12fitf -> 0.023016217095238267`, 
  O14fitf -> -0.010363229861687997`}};
 alO14fitPFOMinnum[4]={4.432950338674559`*^-11, {aALf -> 0.09936857087090613`,
  O4fitf -> 0.003939695798263747`, O6fitf -> -0.026482035301300205`, 
  O8fitf -> 0.024938104748950134`, O10fitf -> -0.025926946566216615`, 
  O12fitf -> 0.020123394724284573`, 
  O14fitf -> -0.008884644245567648`}};
 alO14fitPFOMinnum[5]={2.4821421978499674`*^-11, {aALf -> 0.09911353885417123`, O4fitf -> 0.0026965972962134177`, O6fitf -> -0.023800645452445746`, 
  O8fitf -> 0.022806873723527094`, O10fitf -> -0.023657928734582893`, 
  O12fitf -> 0.018166728957687985`, 
  O14fitf -> -0.007902845351734255`}};
 alO14fitPFOMinnum[6]={2.162145517677897`*^-11, {aALf -> 0.09944498261132279`,  O4fitf -> 0.0015896048791685647`, O6fitf -> -0.02183316122649106`, 
  O8fitf -> 0.021418817457182576`, O10fitf -> -0.022273934027325634`, 
  O12fitf -> 0.01702502218873528`, O14fitf -> -0.007347029504348342`}};
 alO14fitPFOMinnum[7]={2.1390737001023034`*^-11, {aALf -> 0.0999866925852624`,  O4fitf -> 0.0005668811803354681`, O6fitf -> -0.020460466045113167`, 
  O8fitf -> 0.020604470638858785`, O10fitf -> -0.021551303699834778`, 
  O12fitf -> 0.016476080094497646`, 
  O14fitf -> -0.0070950969627611715`}};
alO14fitPFOMinnum[8]={2.106191760880208`*^-11, {aALf -> 0.10047836592504042`, 
  O4fitf -> -0.0003696719554787936`, O6fitf -> -0.019713130748522524`,
   O8fitf -> 0.020339051491874722`, O10fitf -> -0.02143311489658673`, 
  O12fitf -> 0.016452202725026902`, 
  O14fitf -> -0.007106774667222815`}};
alO14fitPFOMinnum[9]={2.1362987802889103`*^-11, {aALf -> 0.10083953524458814`,
  O4fitf -> -0.0012640718790727563`, O6fitf -> -0.019500418408007567`,
   O8fitf -> 0.0204942590228039`, O10fitf -> -0.021740501021064906`, 
  O12fitf -> 0.016776360728737415`, 
  O14fitf -> -0.007285814596871546`}};
alO14fitPFOMinnum[10]={2.339321134711136`*^-11, {aALf -> 0.10086973054763466`,
  O4fitf -> -0.0020919134022933375`, O6fitf -> -0.019840822753056823`,
   O8fitf -> 0.02101513134594308`, O10fitf -> -0.02238015984059938`, 
  O12fitf -> 0.017352733208600812`, 
  O14fitf -> -0.007579977188341899`}};
alO14fitPFOMinnum[11]={2.5054901433463565`*^-11, {aALf -> 0.1008937219993657`,
  O4fitf -> -0.0031831619830700168`, O6fitf -> -0.020194212642430515`,
   O8fitf -> 0.021447019739176686`, O10fitf -> -0.02283905210354748`, 
  O12fitf -> 0.017726992008524806`, 
  O14fitf -> -0.007757828299180296`}};
alO14fitPFOMinnum[12]={2.7366434548673548`*^-11, {aALf -> 0.10042781094585672`,   O4fitf -> -0.00433785808755618`, O6fitf -> -0.020908347192280356`, 
  O8fitf -> 0.02202141736250955`, O10fitf -> -0.023392030516719977`, 
  O12fitf -> 0.01816164450305739`, O14fitf -> -0.007961826004702167`}};
 
 Clear[Nn]; Do[Module[{}, xi2FOO14[Nn]=alO14fitPFOMinnum[Nn][[1]]; aALfFOO14[Nn]=aALf/.alO14fitPFOMinnum[Nn][[2]]; alsmtauFOO14[Nn]=aMSext[mtau^2/(kappa sigmax),aALfFOO14[Nn]] Pi;  O4FOO14[Nn]=O4fitf/.alO14fitPFOMinnum[Nn][[2]]; O6FOO14[Nn]=O6fitf/.alO14fitPFOMinnum[Nn][[2]]; O8FOO14[Nn]=O8fitf/.alO14fitPFOMinnum[Nn][[2]]]; O10FOO14[Nn]=O10fitf/.alO14fitPFOMinnum[Nn][[2]]; O12FOO14[Nn]=O12fitf/.alO14fitPFOMinnum[Nn][[2]]; O14FOO14[Nn]=O14fitf/.alO14fitPFOMinnum[Nn][[2]],{Nn,3,12}]; (* OK *)
 (* 20 Sept. 2022 *)
 (*-------*)
 (* When including O16: *)
 alO16fitPFOMinnum[5]={1.1754514010318373`*^-11, {aALf -> 0.09887716595232399`, O4fitf -> 0.0028469300399070743`, O6fitf -> -0.024513152839930732`, O8fitf -> 0.023801881802943497`, O10fitf -> -0.02563088002028852`, O12fitf -> 0.02152772618228135`, O14fitf -> -0.012201740188257619`, O16fitf -> 0.003072335227133809`}}; alsmtauFOO16[5]=alsmt[0.09887716595232399`]; (* =0.307895 *)
 alO16fitPFOMinnum[6]={1.6548567345953157`*^-11, {aALf -> 0.09929450678561229`, O4fitf -> 0.0016939731634630782`, O6fitf -> -0.022296849519078257`, O8fitf -> 0.022055516038874075`, O10fitf -> -0.023523143067072087`, O12fitf -> 0.019141588537602453`, O14fitf -> -0.010045664714618537`, O16fitf -> 0.00192489379197372`}};  alsmtauFOO16[6]=alsmt[0.09929450678561229`]; (*=0.30918; *)
 alO16fitPFOMinnum[7]={1.8271396866499047`*^-11, {aALf -> 0.09986437802022469`,  O4fitf -> 0.0006580994291502563`, O6fitf -> -0.020837022556830127`, O8fitf -> 0.021113218686053682`, O10fitf -> -0.022540954354466312`, O12fitf -> 0.018146363550343853`, O14fitf -> -0.009220324071451852`, O16fitf -> 0.0015140970971112798`}}; alsmtauFOO16[7]=alsmt[0.09986437802022469`]; (* =0.310934; *)
  alO16fitPFOMinnum[8]={1.7124969483566423`*^-11, {aALf -> 0.1003360120389139`, O4fitf -> -0.0002560190232202324`, O6fitf -> -0.020144996753400472`, O8fitf -> 0.02091445229853505`, O10fitf -> -0.02254674751257154`, O12fitf -> 0.018328861713913822`, O14fitf -> -0.009493436863190489`, O16fitf -> 0.0017001440849350555`}}; alsmtauFOO16[8]=alsmt[0.1003360120389139`]; (* =0.312386) *)
alO16fitPFOMinnum[9]={1.470836093557136`*^-11, {aALf -> 0.10064938755617471`, O4fitf -> -0.0011008472998680057`, O6fitf -> -0.020064390393143363`, O8fitf -> 0.021238115481375074`, O10fitf -> -0.023178882749123286`,O12fitf -> 0.019202660666026578`, O14fitf -> -0.0103747745374971`, O16fitf -> 0.0022022958626741326`}};  alsmtauFOO16[9]=alsmt[0.10064938755617471`]; (*=0.31335; *)

 Clear[Nn]; Do[Module[{}, xi2FOO16[Nn]=alO16fitPFOMinnum[Nn][[1]]; aALfFOO16[Nn]=aALf/.alO16fitPFOMinnum[Nn][[2]]; alsmtauFOO16[Nn]=aMSext[mtau^2/(kappa sigmax),aALfFOO16[Nn]] Pi;  O4FOO16[Nn]=O4fitf/.alO16fitPFOMinnum[Nn][[2]]; O6FOO16[Nn]=O6fitf/.alO16fitPFOMinnum[Nn][[2]]; O8FOO16[Nn]=O8fitf/.alO16fitPFOMinnum[Nn][[2]]]; O10FOO16[Nn]=O10fitf/.alO16fitPFOMinnum[Nn][[2]]; O12FOO16[Nn]=O12fitf/.alO16fitPFOMinnum[Nn][[2]]; O14FOO16[Nn]=O14fitf/.alO16fitPFOMinnum[Nn][[2]]; O16FOO16[Nn]=O16fitf/.alO16fitPFOMinnum[Nn][[2]],{Nn,5,9}]; (* OK *)
  (* 3 Oct. 2022 *)
(*-------*)
 (* When including O18: *)
alO18fitPFOMinnum[5]={3.805791318809418`*^-27, {aALf -> 0.09912056457476118`, 
  O4fitf -> 0.002684453816507215`, O6fitf -> -0.023696925272726297`, 
  O8fitf -> 0.022563550784237114`, O10fitf -> -0.022839983645706965`, 
  O12fitf -> 0.015716995196559228`, O14fitf -> -0.0018271615834885232`,
  O16fitf -> -0.011042597731904264`, O18fitf -> 0.010927811901948392`}}; alsmtauFOO18[5]=alsmt[0.09912056457476118`]; (*= 0.308644 *)
  alO18fitPFOMinnum[6]={6.438161320366478`*^-28, {aALf -> 0.0995910765713232`, O4fitf -> 0.0014788730760735692`, O6fitf -> -0.021284163912297285`, O8fitf -> 0.02054503292118087`, O10fitf -> -0.02015539961400899`, O12fitf -> 0.012170300169985684`, O14fitf -> 0.0023556337177273526`, O16fitf -> -0.014904844111793972`, O18fitf -> 0.013005371580834767`}}; alsmtauFOO18[6]=alsmt[0.0995910765713232`]; (*=0.310093; *)
alO18fitPFOMinnum[7]={5.9010761965196274`*^-27, {aALf -> 0.100188215330497`, 
  O4fitf -> 0.0004064784953273417`, O6fitf -> -0.019735905188814715`, 
  O8fitf -> 0.01949684206758617`, O10fitf -> -0.018967230794953303`, 
  O12fitf -> 0.01077740043423484`, O14fitf -> 0.0038598238216708226`, 
  O16fitf -> -0.016212341128144653`, O18fitf -> 0.013685124759682683`}}; alsmtauFOO18[7]=alsmt[0.100188215330497`]; (* =0.311931 *)
alO18fitPFOMinnum[8]={6.287488558255601`*^-24, {aALf -> 0.10066052143228256`, 
  O4fitf -> -0.000524953642526554`, O6fitf -> -0.019059994224412687`, 
  O8fitf -> 0.01934242864069993`, O10fitf -> -0.019086247090727666`, 
  O12fitf -> 0.01120092950477101`, O14fitf -> 0.003157302814609542`, 
  O16fitf -> -0.015447281536640352`, O18fitf -> 0.013241392301950984`}}; alsmtauFOO18[8]=alsmt[0.10066052143228256`]; (* =0.313384 *)
  alO18fitPFOMinnum[9]={2.2213338641925968`*^-27, {aALf -> 0.1009573233169813`, O4fitf -> -0.0013741984134904076`, O6fitf -> -0.019058494206441963`, O8fitf -> 0.019793258626425638`, O10fitf -> -0.019997808853817786`, O12fitf -> 0.012638977937759536`, O14fitf -> 0.0012944756238937275`, O16fitf -> -0.013637935481476253`, O18fitf -> 0.012246763548894854`}}; alsmtauFOO18[9]=alsmt[0.1009573233169813`]; (*=0.314298 *)
(* 3 Oct. 2022 *)
(*-------*)
 
  (* The convergence of extracted values of alpha_s(mtau^2), in FO at Nt=8, when the number of OPE terms is increased: *)

(* at Nt=8: *)
(* alsmtauFOO8[8]=0.345884; alsmtauFOO10[8]=0.327375; alsmtauFOO12[8]=0.318149; alsmtauFOO14[8]=0.312824; alsmtauFOO16[8]=0.312386; alsmtauFOO18[8]=0.313384; so we see that stabilisation occurs at O14 (inclusion of O16 does not change alpha_s significantly) *)
(* 20 Sept. 2022 *)
(* at Nt=7: *)
(* alsmtauFOO8[7]=0.34374; alsmtauFOO10[7]=0.325368; alsmtauFOO12[7]=0.316433; alsmtauFOO14[7]=0.311311; alsmtauFOO16[7]=0.310934; alsmtauFOO18[7]=0.311931; so we see that stabilisation occurs at O14 (inclusion of O16 does not change alpha_s significantly) *)
(* 20 Sept. 2022 *)
  
(* Now the same at Nt=5: *)
(* alsmtauFOO8[5]=0.344737; alsmtauFOO10[5]=0.323634; alsmtauFOO12[5]=0.314109; alsmtauFOO14[5]=0.308623; alsmtauFOO16[5]=0.307895; alsmtauFOO18[5]=0.308644; so we see, again, that stabilisation occurs at O14 (inclusion of O16 does not change alpha_s significantly) *)

  (* So WE WILL USE the approach with O14, i.e., V-channel with OPE up to (and including) O14; O6 term has nonzero anom. dimension a^(1-ka2)=a^0.222; other OPE contributions (D=4,8,10,12,14) are taken with anom. dimensions zero; Note that D=0 Adler function has IR renormalons u=2,3; u=3 in B[td](u) has the anom. dimension 0.222, and 0.625 i.e., the term in B[td][u] is ~1/(3-u)^ka2 = 1/(3-u)^0.778, which corresponds to a(Q^2)^(1-ka2) <O6>_V/(Q^2)^3 term in OPE needed for the renormalon cancellation; in addition B[td][u] has the term ~1/(3-u)^ka1 = 1/(3-u)^0.375; this term is not reflected in our OPE because FindMinimum in sauch cases leads to numerical instabilities (if we go in OPE up to D=8 or higher) *)
 (* end of BL results, global FO approach with O6) *)
  (*----------------------------------------------------*)
  (*----------------------------------------------------*)
    (*----------------------------------------------------*)
 (* Now FESR (2,n) moments a^{(2,n)} (i.e., double-pinched) in FO: *)

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

(* Now including D >= 4 contributions to a^{(2,n)}(sigmax): *)
(* First a^{(2,0)}(sigmax)  (Nn is the truncation index Nt of the D=0 part): *)
 aaFO20[aALf_][Nn_, O4f_,O6f_]:=aaFOD0[2,0,aALf][Nn]+ ((0+3)/(0+1)) 2 Pi^2 (-1)^0 ( (0+2) O4f/sigmax^(0+2) ) + 3 (1/2/Pi) (6 Pi^2/sigmax^3) ( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 + Exp[I fi]) + (1/3) (1 + Exp[I 3 fi]) ) (O6f (aALf/(1 + I aALf be0 fi))^(1-ka2)),{fi,-Pi,Pi}] );
(* a^{(2,1)}(sigmax): *)
aaFO21[aALf_][Nn_, O6f_,O8f_]:=aaFOD0[2,1,aALf][Nn]+ 2 (1/2/Pi) (6 Pi^2/sigmax^3)*( NIntegrate[ Exp[-I 3 fi] ( Exp[I fi] (1 - Exp[2 I fi]) + (1/2) (1 - Exp[I 4 fi]) ) (O6f (aALf/(1 + I aALf be0 fi))^(1-ka2)),{fi,-Pi,Pi}] ) + ((1+3)/(1+1)) 2 Pi^2 (-1)^1 ( (1+1) O8f/sigmax^(1+3) );			     
(* a^{(2,2)}(sigmax): *)
aaFO22[aALf_][Nn_, O6f_,O8f_,O10f_]:=aaFOD0[2,2,aALf][Nn]+ (6 Pi^2 O6f/sigmax^3)*(1/2/Pi)*( NIntegrate[ Exp[-I 3 fi] ((aALf/(1 + I aALf be0 fi))^(1-ka2)) ( ((2+3)/(2+1)) Exp[I fi] (1 +(-1)^2 Exp[(2+1) I fi]) +  (1 + (-1)^2 Exp[I (2+3) fi]) ),{fi,-Pi,Pi}] ) + ((2+3)/(2+1)) 2 Pi^2 (-1)^2 ( (2+2) O8f/sigmax^(2+2) + (2+1) O10f/sigmax^(2+3));
(* a^{(2,3)}(sigmax): *)
aaFO23[aALf_][Nn_, O6f_,O10f_,O12f_]:=aaFOD0[2,3,aALf][Nn]+ (6 Pi^2 O6f/sigmax^3)*(1/2/Pi)*( NIntegrate[ Exp[-I 3 fi] ((aALf/(1 + I aALf be0 fi))^(1-ka2)) ( ((3+3)/(3+1)) Exp[I fi] (1 +(-1)^3 Exp[(3+1) I fi]) +  (1 + (-1)^3 Exp[I (3+3) fi]) ),{fi,-Pi,Pi}] ) + ((3+3)/(3+1)) 2 Pi^2 (-1)^3 ( (3+2) O10f/sigmax^(3+2) + (3+1) O12f/sigmax^(3+3));
(* a^{(2,4)}(sigmax): *)
aaFO24[aALf_][Nn_, O6f_,O12f_,O14f_]:=aaFOD0[2,4,aALf][Nn]+ (6 Pi^2 O6f/sigmax^3)*(1/2/Pi)*( NIntegrate[ Exp[-I 3 fi] ((aALf/(1 + I aALf be0 fi))^(1-ka2)) ( ((4+3)/(4+1)) Exp[I fi] (1 +(-1)^4 Exp[(4+1) I fi]) +  (1 + (-1)^4 Exp[I (4+3) fi]) ),{fi,-Pi,Pi}] ) + ((4+3)/(4+1)) 2 Pi^2 (-1)^4 ( (4+2) O12f/sigmax^(4+2) + (4+1) O14f/sigmax^(4+3));
(* 10 Sept. 2022 *)
(*----------------------------------------------------*)
(* For the 'O14' case (i.e., our MAIN considered case), finding optimal Nt: *)
(* Tab20FOO14avar=Table[{NN,Re[aaFO20[aALfFOO14[NN]][NN,O4FOO14[NN],O6FOO14[NN]]]},{NN,3,12}];
Tab21FOO14avar=Table[{NN,Re[aaFO21[aALfFOO14[NN]][NN,O6FOO14[NN],O8FOO14[NN]]]},{NN,3,12}]; 
 Tab22FOO14avar=Table[{NN,Re[aaFO22[aALfFOO14[NN]][NN,O6FOO14[NN],O8FOO14[NN],O10FOO14[NN]]]},{NN,3,12}];
 Tab23FOO14avar=Table[{NN,Re[aaFO23[aALfFOO14[NN]][NN,O6FOO14[NN],O10FOO14[NN],O12FOO14[NN]]]},{NN,3,12}];
  Tab24FOO14avar=Table[{NN,Re[aaFO23[aALfFOO14[NN]][NN,O6FOO14[NN],O12FOO14[NN],O14FOO14[NN]]]},{NN,3,12}]; 
 *)	

(* 
 Tab20FOO14avar(*={{3, 0.23069}, {4, 0.230686}, {5, 0.230683}, {6, 0.230682}, {7, 0.230681}, {8, 0.230681}, {9, 0.230682}, {10, 0.230683}, {11, 0.230683}, {12, 0.230684}}*)={{3, 0.23069025186101244`}, {4, 0.23068573432130154`}, {5, 0.230682936249899`}, {6, 0.23068152385736435`}, {7, 0.230681024871472`}, {8, 0.2306812616442459`}, {9, 0.23068188929082137`}, {10, 0.23068273576898882`}, {11, 0.23068311230110386`}, {12, 0.23068357238847348`}};  (* minimum at Nt=7 {7, 0.2306810}, at Nt=8 slightly higher {8, 0.2306812}; but local minimal sensitivity at Nt=8; see also: delTab20FOO14avar= Table[{Nn, Tab20FOO14avar[[Nn - 2]][[2]] - Tab20FOO14avar[[Nn - 3]][[2]]}, {Nn, 4, 12}];  Local minimum of a20 at Nt=7, minimal sensitivity at Nt=7 and above *)
 Tab21FOO14avar(*={{3, 0.232013}, {4, 0.231995}, {5, 0.231984}, {6, 0.231978}, {7, 0.231976}, {8, 0.231977}, {9, 0.23198}, {10, 0.231983}, {11, 0.231984}, {12, 0.231986}}*)={{3, 0.23201261244846627`}, {4, 0.23199478530924272`}, {5, 0.2319837562612334`}, {6, 0.23197819981241696`}, {7, 0.23197624589850063`}, {8, 0.23197719050900453`}, {9, 0.2319796701120207`}, {10, 0.2319830059774161`}, {11, 0.23198448145958328`}, {12, 0.23198628887860745`}}; (* minimum at Nt=7 {7, 0.2319762}, at Nt=8 slightly higher {8, 0.2319771}; but local minimal sensitivity at Nt=8; see also: delTab21FOO14avar= Table[{Nn, Tab21FOO14avar[[Nn - 2]][[2]] - Tab21FOO14avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; *)
 Tab22FOO14avar(*={{3, 0.2037}, {4, 0.203643}, {5, 0.203608}, {6, 0.203591}, {7, 0.203585}, {8, 0.203588}, {9, 0.203595}, {10, 0.203606}, {11, 0.203611}, {12, 0.203616}}*)={{3, 0.20369963535735408`}, {4, 0.2036431259945516`}, {5, 0.20360821969402496`}, {6, 0.2035906795119538`}, {7, 0.20358455039194143`}, {8, 0.2035875880972619`}, {9, 0.20359545951763244`}, {10, 0.2036060140654741`}, {11, 0.20361064637893617`}, {12, 0.203616340835449`}}; (* minimum at Nt=7 {7, 0.203585}; but local minimal sensitivity at Nt=8; see also: delTab22FOO14avar= Table[{Nn, Tab22FOO14avar[[Nn - 2]][[2]] - Tab22FOO14avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; *)
 Tab23FOO14avar={{3, 0.183049}, {4, 0.182899}, {5, 0.182806}, {6, 0.182759}, {7, 0.182743}, {8, 0.182751}, {9, 0.182772}, {10, 0.1828}, {11, 0.182813}, {12, 0.182828}}; (* minimum at Nt=7  {7, 0.182743}; but local minimal sensitivity at Nt=8; see also: delTab23FOO14avar= Table[{Nn, Tab23FOO14avar[[Nn - 2]][[2]] - Tab23FOO14avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; *)
 Tab24FOO14avar={{3, 0.158959}, {4, 0.161586}, {5, 0.163409}, {6, 0.164516}, {7, 0.165088}, {8, 0.165173}, {9, 0.164917}, {10, 0.164398}, {11, 0.164028}, {12, 0.163589}}; (* local maximum at Nt=8 {8, 0.165173} (at Nt=7 a close value {7, 0.165088}); but experimental medium values are not well reproduced at any Nt (!); so O14 OPE contribution (in the O14 OPE truncated appoach) appears to be a collective effect of all terms with D >= 14, and not really D=14 contribution *)
 (* but local minimal sensitivity at Nt=8; see also: delTab24FOO14avar= Table[{Nn, Tab24FOO14avar[[Nn - 2]][[2]] - Tab24FOO14avar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; *)
 *)
 (* aexpBinComb[2, 0]+-delaexpBinComb[2, 0]=0.230681 +-0.088331; 
  aexpBinComb[2, 1]+-delaexpBinComb[2, 1]=0.231976 +-0.085489; 
  aexpBinComb[2, 2]+-delaexpBinComb[2, 2]=0.203590 +- 0.086436;
  aexpBinComb[2, 3]+-delaexpBinComb[2, 3]=0.182778 +- 0.090541;
  aexpBinComb[2, 4]+-delaexpBinComb[2, 4]=0.171290 +- 0.096610; *) 

 (* Therefore, the stability under variation of Nt is achieved at Nt=8 (FO, OPE at O14) *) (* So we take Nt=8, in O14 FO approach ('p3ka1ka2 & O6ka2 O14' FO approach) *)
 (* FO until here on 24 Sept. 2022 *)
(*-----------------------------------------------*)
(*---------------------------------------------*)
(*---------------------------------------------*)

(* Now the estimate of EXPERIMENTAL uncertainty in this global approach with double-pinched FO Borel-Laplace (with Nt=8): *)
 (* See Appendix  PRD84, 113006 (Boito et al., 2011): *)
 
(* Now the M^2_{al}: all 9 values (three rays in the complex M^2-plane) as one array, as M2al[jj] (jj=1,...,9): *)
 (* Note, we have NchP=2; Nm=NchP+1; Nd=3 Nm (=9); *)

(* The central considered FO case: O14 (OPE includes O14) and Nt=8: *) 
(* alO14fitPFOMinnum[8]={2.106191760880208`*^-11, {aALf -> 0.10047836592504042`, O4fitf -> -0.0003696719554787936`, O6fitf -> -0.019713130748522524`, O8fitf -> 0.020339051491874722`, O10fitf -> -0.02143311489658673`, O12fitf -> 0.016452202725026902`, O14fitf -> -0.007106774667222815`}};  alsmtauFOO14[8]=0.312824; *)

N0=8; (* N0=Nt *)  Np=7; (* Np is the number of parameters: a(sigmax) and (6) condensates *)
(* First derivatives of ReBPFONNumO14[M2] with respect to p1=aAL, p2=O4, p3=O6, p4=O8, p5=O10, p6= O12, p7=O14 (Note: O2n=2*<O_{2 n}>_V) : *)
epp1=0.00002;
dReBPFOd[1][M2_]:=(ReBPFONumO14[M2, N0, aALfFOO14[N0](1+epp1),0,O6FOO14[N0],0,0,0,0]-ReBPFONumO14[M2, N0, aALfFOO14[N0](1-epp1),0,O6FOO14[N0],0,0,0,0])/(2 epp1 aALfFOO14[N0]);
dReBPFOd[2][M2_]:= Re[ 2 Pi^2 (1/M2^2) ( 1 + 2 M2/sigmax ) ];
dReBPFOd[3][M2_]:=Re[(6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ ( (aALfFOO14[N0]/(1 + I be0 aALfFOO14[N0] fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2/sigmax) (1 + Exp[I fi]) + 2 (M2/sigmax)^2) Exp[(sigmax/M2) Exp[I fi]] - 2 (M2/sigmax)^2 Exp[- sigmax/M2] ),{fi, -Pi, Pi}]];
dReBPFOd[4][M2_]:=Re[2 Pi^2  (1/3!) (1/M2^4) (1+ 2*3 (M2/sigmax) + 3*2 (M2/sigmax)^2)];
dReBPFOd[5][M2_]:=Re[2 Pi^2  (1/4!) (1/M2^5) (1+ 2*4 (M2/sigmax) + 4*3 (M2/sigmax)^2)];
dReBPFOd[6][M2_]:=Re[ 2 Pi^2  (1/5!) (1/M2^6) (1+ 2*5 (M2/sigmax) + 5*4 (M2/sigmax)^2)];
dReBPFOd[7][M2_]:=Re[2 Pi^2 (1/6!) (1/M2^7) (1+ 2*6 (M2/sigmax) + 6*5 (M2/sigmax)^2)];
(* Now Tables of derivatives, for all Mal[ii]: *)
 tabdReBPFOd[1] = Table[{M2al[ii], dReBPFOd[1][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[2] = Table[{M2al[ii], dReBPFOd[2][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[3] = Table[{M2al[ii], dReBPFOd[3][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[4] = Table[{M2al[ii], dReBPFOd[4][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[5] = Table[{M2al[ii], dReBPFOd[5][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[6] = Table[{M2al[ii], dReBPFOd[6][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[7] = Table[{M2al[ii], dReBPFOd[7][M2al[ii]]}, {ii, 1, Nd}];
 (* The Table tabCovReBPexpBinComb requires longer time, so better evaluate it only once, namely here, and then use it for all the methods: *)
 tabCovReBPexpBinComb = Table[{M2al[ii],M2al[jj],CovReBPexpBinComb[M2al[ii],M2al[jj]]},{ii,1,Nd},{jj,1,Nd}];
mCovReBPexpBinComb=Table[tabCovReBPexpBinComb[[i, j]][[3]],{i,1,Nd},{j,1,Nd}];
(* This is the NdxNd matrix form of the covariances <Re[B(M^2_{al})] Re[B(M^2_{be})]; it shows very strong correlations of various Re[B(M^2_{al})] due to experimental uncertainties of data Omega_V(sig); all terms of mCovReBPexpBinComb are (4-9)x 10^(-6) *)
(* Now calculate the NpxNp mtrx AA(=A), for the FO case, using tables obtained above: *)
 tabAAFO=Table[{ik, jk, Sum[(tabdReBPFOd[ik][[jj]][[2]] tabdReBPFOd[jk][[jj]][[2]])/tabCovReBPexpBinComb[[jj]][[jj]][[3]], {jj, 1, Nd}]}, {ik, 1, Np}, {jk, 1, Np}];
 Do[AAFOa[ii,jj]=tabAAFO[[ii]][[jj]][[3]],{ii,1,Np},{jj,1,Np}];  (* This is NpxNp matrix A (or: AA), calculated more efficiently *)
AAFOm=Table[AAFOa[ik,jk],{ik,1,Np},{jk,1,Np}]; (* all these terms are =app 10^4-10^8 *)
AAFOmI=Inverse[AAFOm];   (* AAFOmI = (AAFOm)^{-1}, inverse NpxNp matrix of AAm *)
(* WARNING APPREARS: Result for Inverse of badly conditioned matrix; then the results are completely unreliable and numerically unstable *)

(* Now the NpxNp covariance matrix of the fitted parameters pj: *) 
(* p1=aAL, p2=<O4>, p3=<O6>, p4=<O8>, p5=<O10>, p6=<O12>, p7=<O14>; where <O_{2n}> means: 2 <O_{2n}>_V  *)
(* PP[i,j]= < (delta p_i) (delta p_j) >: *)
  (* We calulate PP[i,j] covariance matrix of the extracted parameters, using the tables obtained above: *) 
 PPtabFO[jp_,kp_]:=Sum[AAFOmI[[jp,jpp]] AAFOmI[[kp,kpp]] tabdReBPFOd[jpp][[Al]][[2]] tabdReBPFOd[kpp][[Be]][[2]] tabCovReBPexpBinComb[[Al]][[Be]][[3]]/(tabCovReBPexpBinComb[[Al]][[Al]][[3]]  tabCovReBPexpBinComb[[Be]][[Be]][[3]]),{jpp,1,Np},{kpp,1,Np},{Be,1,Nd},{Al,1,Nd}];

(* I then get, for epp1=0.001: 
 Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0. + 0.00061228 I}, {2, 0.00193715}, {3, 0.00158199}, {4, 0. + 0.00239748 I}, {5, 0.00560637}, {6, 0. + 0.00508427 I}, {7, 0. + 0.00813153 I}}; several nonreal (!) values *)
(* I then get, for epp1=0.0005: 
Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0. + 0.00114117 I}, {2, 0.00138266}, {3, 0.00495152}, {4, 0.00410643}, {5, 0. + 0.00901797 I}, {6, 0.00739782}, {7, 0. + 0.00878564 I}}; several nonreal (!) values *)
(* For epp1=0.0001: 
Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0.00325555}, {2, 0. + 0.00120972 I}, {3, 0.00262261}, {4, 0. + 0.00383353 I}, {5, 0.008008}, {6, 0.00757648}, {7, 0.0101994}}; several nonreal (!) values *)
(* For epp1=0.00005: 
 Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0.00258873}, {2, 0.00249193}, {3, 0.00911064}, {4, 0.00862283}, {5, 0.00758644}, {6, 0.00291814}, {7, 0.00691348}}; *)
(* For epp1=0.00002:
 Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0.00216859}, {2, 0.00195746}, {3, 0.00639061}, {4, 0.00277211}, {5, 0.0101939}, {6, 0.0138473}, {7, 0. + 0.00999415 I}}; *)
(* UNSTABLE; but the above numbers suggest that: del(a(sigmax))_{exp}=Sqrt[PPtabFO[1, 1]] is probably =app (1-5)*10^-3; but UNRELIABLE *)

(* When we use Nd=12 (instead of Nd=9) --> alsmtauFOO14[8]=0.312894; O4FOO14[8]=-0.000386797; etc. (i.e., slightly different from the Nd=9 'canonical' case), then we get, for epp1=0.00002, again instability warning  for inversion of AAFOm, and: Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0.00191212}, {2, 0.00114922}, {3, 0.0046047}, {4, 0. + 0.00421972 I}, {5, 0.00620421}, {6, 0.00628008}, {7, 0. + 0.00438692 I}}; i.e., Sqrt[PPtabFO[1, 1]]=0.00191212 close to the case of Nd=9 (Sqrt[PPtabFO[1, 1]]=0.00216859), and even closer to the (stable) Np=6 case (Sqrt[PPtabFO[1, 1]]=0.00171157; where O14 was kept fixed for numerical stability of inversion of AAFOm) than when Nd=9 (and Np=7) *)
(*------------------------------------*)
(* Nonetheless, if we keep O14 FIXED, at the central extracted value O14FOO14[N0](=-0.0070951), then the above formalism leads to 6x6 matrix (nP=6; i.e., variations of parameters up to p6=O12 then I get): *)
 (* EXPERIMENTAL UNCERTAINTIES, KEEPING O14FOO14[N0](=-0.0070951) FIXED, i.e., Np=6: *)
 (* Now calculate the NpxNp mtrx AA(=A), for the FO case, using tables obtained above (for epp1=:0.00002) *)
 Np=6;
 tabAAFO=Table[{ik, jk, Sum[(tabdReBPFOd[ik][[jj]][[2]] tabdReBPFOd[jk][[jj]][[2]])/tabCovReBPexpBinComb[[jj]][[jj]][[3]], {jj, 1, Nd}]}, {ik, 1, Np}, {jk, 1, Np}];
 Do[AAFOa[ii,jj]=tabAAFO[[ii]][[jj]][[3]],{ii,1,Np},{jj,1,Np}];  (* This is NpxNp matrix A (or: AA), calculated more efficiently *)
AAFOm=Table[AAFOa[ik,jk],{ik,1,Np},{jk,1,Np}]; 
AAFOmI=Inverse[AAFOm]; (* AAFOmI = (AAFOm)^-1 Np=6; *)
 tabAAFO=Table[{ik, jk, Sum[(tabdReBPFOd[ik][[jj]][[2]] tabdReBPFOd[jk][[jj]][[2]])/tabCovReBPexpBinComb[[jj]][[jj]][[3]], {jj, 1, Nd}]}, {ik, 1, Np}, {jk, 1, Np}];
 Do[AAFOa[ii,jj]=tabAAFO[[ii]][[jj]][[3]],{ii,1,Np},{jj,1,Np}];  (* This is NpxNp matrix A (or: AA), calculated more efficiently *)
AAFOm=Table[AAFOa[ik,jk],{ik,1,Np},{jk,1,Np}]; 
AAFOmI=Inverse[AAFOm]; (* AAFOmI = (AAFOm)^{-1}, inverse NpxNp matrix of AAm *)
 (* Now: NO WARNING 'Result for Inverse of badly conditioned matrix'; GOOD *)
 (* p1=aAL, p2=O4, p3=O6, p4=O8, p5=O10, p6=O12; where O_{2n} means: 2 <O_{2n}>_V  *)
(* PP[i,j]= < (delta p_i) (delta p_j) >: *)
 PPtabFO[jp_,kp_]:=Sum[AAFOmI[[jp,jpp]] AAFOmI[[kp,kpp]] tabdReBPFOd[jpp][[Al]][[2]] tabdReBPFOd[kpp][[Be]][[2]] tabCovReBPexpBinComb[[Al]][[Be]][[3]]/(tabCovReBPexpBinComb[[Al]][[Al]][[3]]  tabCovReBPexpBinComb[[Be]][[Be]][[3]]),{jpp,1,Np},{kpp,1,Np},{Be,1,Nd},{Al,1,Nd}];
 (* For epp1=0.0001: 
 Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0.00171248}, {2, 0.000825929}, {3, 0.00182014}, {4, 0.00135427}, {5, 0.00125194}, {6, 0.000686739}}; *)
 (* For epp1=0.00002: (I will use this epp1) ******
 Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0.00171157}, {2, 0.000825789}, {3, 0.00181902}, {4, 0.00135496}, {5, 0.00124828}, {6, 0.000690433}}; *)
  ***********   *)
 (* GOOD stability; Sqrt[PPtabFO[1, 1]]=0.00165, really in the above estimated range of (1-3)*10^-3. *)
 (* Further, if repeating the analysis with O12 OPE results (and the correponding Np=6): Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, 6}]={{1, 0.00162681}, {2, 0.000724776}, {3, 0.00174728}, {4, 0.00134608}, {5, 0.00125155}, {6, 0.000690514}} (see below); this is quite similar to the above; so we will use the above for del(a(sigmax)); del(O4), del(O6), del(O8), del(O10); and delP[j] (j=6,7) for del(O12), del(O14) *)

 (* Further, if taking Np=5, i.e., keeping O14 and O12 fixed, we get (for epp1=0.00002): Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0.00152144}, {2, 0.000634067}, {3, 0.00112672}, {4, 0.000606913}, {5, 0.000308437}}; <p1, p1> similar to the above *)
 (*----------------------------------------*)
 (* Another approximate implementation should be to take eq.(15) of my notes [eq. (A4) of Boito et al App. PRD84,113006] with (delta d_{alpha} = delReBPexpBinAL[M2_{alpha}]  (=sigma_{alpha})); *)
 Np=7; (* calculate for Np=7 the mtrx AAFOmI *)
 delP[jp_]:=Sum[AAFOmI[[jp,jk]] tabdReBPFOd[jk][[Al]][[2]]/(Sqrt[CovReBPexpBinComb[M2al[Al],M2al[Al]]]),{jk,1,7},{Al,1,Nd}];
 (* Although AAFOmI is unstable (for Np=7), the above expressions delP[jj] are NOT; this probably has to do with the fact that delP[jj] give the same uncertainties that we obtain when we assume complete correlations of the ReBPexp *)
 (* --> delP[1] (= +-del(a(sigmax))_{exp}) =0.004011; delP[2](=+-del(O4)_{exp})=0.003472; delP[3](=+-del(O6)_{exp})=-0.011835; delP[4](=+-del(O8)_{exp})=0.012947; delP[5](=+-del(O10)_{exp})=-0.018173; delP[6](=+-del(O12)_{exp})=0.018394; delP[7](=+-del(O14)_{exp})=-0.009125; these results suggest that the above  Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}] (Np=6) estimates are probably some underestimation for O4, O6, O8, O10, O12; and here delP[j] are some overestimations *)
 (* The results are stable (independent) of epp1 (the same for epp1=0.0001); this is interesting, because the results of Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}] (for Np=7) are unstable and depend on epp1; it seems that although AAFOmI and tabdReBPFOd[jpp] are unstable under variation epp1 (for Np=7), the combination appearing in delP[jj] is stable; this has to do probably with the fact that del[jj] results are also obtained if doing FindMinimum to the upper and lower bound curve of the experimental band of Re[B{M^2]], which is a stable procedure (even with Np=7):
  alO14fitPFOMinNum1lUP[Nn_]:= FindMinimum[Sum[(ReBPFONumO14[M2al[jj], Nn, aALf, O4fitf, O6fitf, O8fitf, O10fitf, O12fitf, O14fitf] - (ReBPexpBinComb[M2al[jj]]+delReBPexpBinComb[M2al[jj]]))^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf, aALfinit}, {O4fitf, O4init}, {O6fitf, O6init}, {O8fitf, O8init}, {O10fitf, O10init}, {O12fitf, O12init}, {O14fitf, O14init}}];
 alO14fitPFOMinNum1lDO[Nn_]:= FindMinimum[Sum[(ReBPFONumO14[M2al[jj], Nn, aALf, O4fitf, O6fitf, O8fitf, O10fitf, O12fitf, O14fitf] - (ReBPexpBinComb[M2al[jj]]-delReBPexpBinComb[M2al[jj]]))^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{aALf, aALfinit}, {O4fitf, O4init}, {O6fitf, O6init}, {O8fitf, O8init}, {O10fitf, O10init}, {O12fitf, O12init}, {O14fitf, O14init}}];
  and then comparing the results of  alO14fitPFOMinNum1lUP[7] (and alO14fitPFOMinNum1lDO[7]) with those of alO14fitPFOMinNum[7] --> we also get del(a(sigmax)) =app -+0.0040, ... del(O14) =app  -+ 0.0091 
  *)
 (*-----------------------------------------------------------*)
 (* If we do the experimental error analysis in the O12 approach (and the corresponding Np=6), we get similar results (as in the O14 approach with O14 fixed): *)
 N0=8; (* N0=Nt *)  Np=6; (* Np is the number of parameters: a(sigmax) and (6) condensates *)
epp1=0.00002;
dReBPFOd[1][M2_]:=(ReBPFONumO12[M2, N0, aALfFOO12[N0](1+epp1),0,O6FOO12[N0],0,0,0]-ReBPFONumO12[M2, N0, aALfFOO12[N0](1-epp1),0,O6FOO12[N0],0,0,0])/(2 epp1 aALfFOO12[N0]);
dReBPFOd[2][M2_]:= Re[ 2 Pi^2 (1/M2^2) ( 1 + 2 M2/sigmax ) ];
dReBPFOd[3][M2_]:=Re[(6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ ( (aALfFOO12[N0]/(1 + I be0 aALfFOO12[N0] fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2/sigmax) (1 + Exp[I fi]) + 2 (M2/sigmax)^2) Exp[(sigmax/M2) Exp[I fi]] - 2 (M2/sigmax)^2 Exp[- sigmax/M2] ),{fi, -Pi, Pi}]];
dReBPFOd[4][M2_]:=Re[2 Pi^2  (1/3!) (1/M2^4) (1+ 2*3 (M2/sigmax) + 3*2 (M2/sigmax)^2)];
dReBPFOd[5][M2_]:=Re[2 Pi^2  (1/4!) (1/M2^5) (1+ 2*4 (M2/sigmax) + 4*3 (M2/sigmax)^2)];
dReBPFOd[6][M2_]:=Re[ 2 Pi^2  (1/5!) (1/M2^6) (1+ 2*5 (M2/sigmax) + 5*4 (M2/sigmax)^2)];
 tabdReBPFOd[1] = Table[{M2al[ii], dReBPFOd[1][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[2] = Table[{M2al[ii], dReBPFOd[2][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[3] = Table[{M2al[ii], dReBPFOd[3][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[4] = Table[{M2al[ii], dReBPFOd[4][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[5] = Table[{M2al[ii], dReBPFOd[5][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPFOd[6] = Table[{M2al[ii], dReBPFOd[6][M2al[ii]]}, {ii, 1, Nd}];
 tabCovReBPexpBinComb = Table[{M2al[ii],M2al[jj],CovReBPexpBinComb[M2al[ii],M2al[jj]]},{ii,1,Nd},{jj,1,Nd}];
tabAAFO=Table[{ik, jk, Sum[(tabdReBPFOd[ik][[jj]][[2]] tabdReBPFOd[jk][[jj]][[2]])/tabCovReBPexpBinComb[[jj]][[jj]][[3]], {jj, 1, Nd}]}, {ik, 1, Np}, {jk, 1, Np}]; Do[AAFOa[ii,jj]=tabAAFO[[ii]][[jj]][[3]],{ii,1,Np},{jj,1,Np}];  
AAFOm=Table[AAFOa[ik,jk],{ik,1,Np},{jk,1,Np}]; 
AAFOmI=Inverse[AAFOm];   
 (* p1=aAL, p2=<O4>, p3=<O6>, p4=<O8>, p5=<O10>, p6=<O12>: *)
 PPtabFO[jp_,kp_]:=Sum[AAFOmI[[jp,jpp]] AAFOmI[[kp,kpp]] tabdReBPFOd[jpp][[Al]][[2]] tabdReBPFOd[kpp][[Be]][[2]] tabCovReBPexpBinComb[[Al]][[Be]][[3]]/(tabCovReBPexpBinComb[[Al]][[Al]][[3]]  tabCovReBPexpBinComb[[Be]][[Be]][[3]]),{jpp,1,Np},{kpp,1,Np},{Be,1,Nd},{Al,1,Nd}];

(* I get, for epp1=0.00002:
 Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]={{1, 0.0016926}, {2, 0.000864476}, {3, 0.0018469}, {4, 0.0013872}, {5, 0.00126695}, {6, 0.000694515}}; when epp1=0.0001, less than 1% variation; thus STABLE under variation of epp1 *)
 (*-------------------------------------------------------------*)
 (* Thus, summary of experimental uncertainties: *)
 (* for del(a(sigmax)), del(O4), ..., del(O12) we use the matrix correlated fulla analysis with O14 fixed at it central value O14FOO14[8]; and for del(O14) the results of the fully correlated (complete correlation assumption) approach (del(O14))=delP[7]): *)
 (* For epp1=0.00002:    *****
Table[{nn, Sqrt[PPtabFO[nn, nn]]}, {nn, 1, Np}]=1, 0.00171157}, {2, 0.000825789}, {3, 0.00181902}, {4, 0.00135496}, {5, 0.00124828}, {6, 0.000690433}};
delP[7](=+-del(O14)_{exp})=delP7=-0.009125;
  *****   *)
(* N0=8; a(sigmax)=aALfFOO14[N0]+-Sqrt[PPtabFO[1, 1]]=0.100478+-0.001712; O4=O4FOO14[N0]+-Sqrt[PPtabFO[2, 2]]; ...; O12=O12FOO14[N0]+-Sqrt[PPtabFO[6, 6]]; O14=O14FOO14[N0]+-delP[7]: *)
(* Summary of the experimental uncertainties, 'O14' approach with Nt=7: *)
(* 'p3ka1ka2 & O6ka2 & O14' FO Nt=N0(=8) approach: *)
(*
 alpha(mtau^2)=0.312824+-0.005269{exp};
 10^3 O4=  -0.37 +-0.83_{exp};
 10^3 O6= -19.71 +-1.82_{exp};
 10^3 O8= +20.34 +-1.35_{exp};
 10^3 O10=-21.43 +-1.25_{exp};
 10^3 O12=+16.45 +-0.69_{exp};
 10^3 O14= -7.11 +-9.12_{exp};
 *)
(* 24 Sept. 2022; checked 26 Sept. 2022 *)
 (*----------------------------------------------------------_*)
   (* Variation of kappa, FO global (at Nt=N0(=8)): *)
   (* kappa=2.0: (Note, aALf=a(sigmax*kappa))  ****  *)
(* alO14fitPFOMinnum[N0]={2.18365*10^-11, {aALf -> 0.0839948, O4fitf -> 0.00115624, O6fitf -> -0.021589, O8fitf -> 0.0214373, O10fitf -> -0.0223973, O12fitf -> 0.0171835, O14fitf -> -0.00744276}}; a(sigmax)=0.0999768; alpha_s(mtau^2)=0.31128; *)
(* kappa=1.0 (the central case of the global fit): ****  *)
(* alO14fitPFOMinnum[N0]={2.10619*10^-11, {aALf -> 0.100478, O4fitf -> -0.000369672, O6fitf -> -0.0197131, O8fitf -> 0.0203391, O10fitf -> -0.0214331, O12fitf -> 0.0164522, O14fitf -> -0.00710677}}; a(sigmax)=0.100478; alpha_s(mtau^2)=0.312824; *)
(* kappa=0.5: ****  *)
(*  alO14fitPFOMinnum[N0]={2.76171*10^-11, {aALf -> 0.128197, O4fitf -> -0.00298092, O6fitf -> -0.0190505, O8fitf -> 0.0209217, O10fitf -> -0.0225893, O12fitf -> 0.0176882, O14fitf -> -0.0077967}}; a(sigmax)=0.101727; alpha_s(mtau^2)=0.316667; *)

(* --> thus the variations: (kappa=1+{^{+1}_{-1/2}}) *)
(*    alpha_s(mtau^2)=0.312824+{^{-0.001544}_{+0.003843}}_{kap}; 
      10^3 O4=  -0.37+{^{+1.53}_{-2.61}}_{kap};  
      10^3 O6= -19.71+{^{-1.88}_{+0.66}}_{kap};  
      10^3 O8= +20.34+{^{+1.10}_{+0.58}}_{kap}; 
      10^3 O10=-21.43+{^{-0.97}_{-1.16}}_{kap}; 
      10^3 O12=+16.45+{^{+0.73}_{+1.24}}_{kap}; 
      10^3 O14= -7.11+{^{-0.33}_{-0.69}}_{kap};
 *)
(* 24 Sept. 2022; checked 26 Sept. 2022 *)
(*-----------*)
 (* Variation of d4, FO global (NOTE, at Nt=N0(=7)) d4=275 +- 63: *)

 (* when ddtk[4, tkap]=275.+63.=338.: 
  alO14fitPFOMinNum[N0]={2.12759*10^-11, {aALf -> 0.0981615, O4fitf -> 0.0000553346, O6fitf -> -0.0202435, O8fitf -> 0.0205175, O10fitf -> -0.0215588, O12fitf -> 0.0165188, O14fitf -> -0.00712651}}; alpha_s(mtau^2)=0.305691; *)
 
(* when the central case d4=ddtk[4,tkap]=ddtk[4,Exp[tKK]]=275.: 
 alO14fitPFOMinnum[N0]={2.10619*10^-11, {aALf -> 0.100478, O4fitf -> -0.000369672, O6fitf -> -0.0197131, O8fitf -> 0.0203391, O10fitf -> -0.0214331, O12fitf -> 0.0164522, O14fitf -> -0.00710677}};  alpha_s(mtau^2)=0.312824; *)

(* when ddtk[4, tkap]=ddtk[4,Exp[tKK]]=275.-63.=212: 
 alO14fitPFOMinNum[N0]={2.10598*10^-11, {aALf -> 0.101765, O4fitf -> -0.000506141, O6fitf -> -0.019695, O8fitf -> 0.0204765, O10fitf -> -0.0216289, O12fitf -> 0.0166516, O14fitf -> -0.00721686}};  alpha_s(mtau^2)=0.316782; *)

(* alpha_s(mtau^2)=0.312824+{^{-0.00713246}}_{+0.00395858}}_{d4}; 
   10^3 O4=  -0.37+{^{+0.43}_{-0.14}}_{d4};  
   10^3 O6= -19.71+{^{-0.53}_{+0.01}}_{d4};  
   10^3 O8= +20.34+{^{+0.18}_{+0.14}}_{d4}; 
   10^3 O10=-21.43+{^{-0.13}_{-0.20}}_{d4}; 
   10^3 O12=+16.45+{^{+0.07}_{+0.20}}_{d4}; 
   10^3 O14= -7.11+{^{-0.02}_{-0.11}}_{d4}; 
    *)
    (* 24 Sept. 2022; checked 26 Sept. 2022 *)
 (*--------------*)
    (* Variation of Nt(truncation order, Nt=8+-1 and Nt=8+-2, 8-3) --> 8+2-3, FO global: *)
(*
 alO14fitPFOMinnum[10]={2.33932*10^-11, {aALf -> 0.10087, O4fitf -> -0.00209191, O6fitf -> -0.0198408, O8fitf -> 0.0210151, O10fitf -> -0.0223802, O12fitf -> 0.0173527, O14fitf -> -0.00757998}}; alsmtauFOO14[10]=0.314028;
  alO14fitPFOMinnum[9]={2.1363*10^-11, {aALf -> 0.10084, O4fitf -> -0.00126407,  O6fitf -> -0.0195004, O8fitf -> 0.0204943, O10fitf -> -0.0217405, 
  O12fitf -> 0.0167764, O14fitf -> -0.00728581}}; alsmtauFOO14[9]=0.313935;
The central case:
  alO14fitPFOMinnum[8]={2.10619*10^-11, {aALf -> 0.100478, O4fitf -> -0.000369672, O6fitf -> -0.0197131, O8fitf -> 0.0203391, O10fitf -> -0.0214331, 
  O12fitf -> 0.0164522, O14fitf -> -0.00710677}}; alsmtauFOO14[8]=0.312824;
  alO14fitPFOMinnum[7]={2.13907*10^-11, {aALf -> 0.0999867, O4fitf -> 0.000566881, O6fitf -> -0.0204605, O8fitf -> 0.0206045, O10fitf -> -0.0215513, 
  O12fitf -> 0.0164761, O14fitf -> -0.0070951}}; alsmtauFOO14[7]=0.311311;
 alO14fitPFOMinnum[6]={2.16215*10^-11, {aALf -> 0.099445, O4fitf -> 0.0015896, O6fitf -> -0.0218332, O8fitf -> 0.0214188, O10fitf -> -0.0222739, 
  O12fitf -> 0.017025, O14fitf -> -0.00734703}}; alsmtauFOO14[6]=0.309643;
  alO14fitPFOMinnum[5]={2.48214*10^-11, {aALf -> 0.0991135, O4fitf -> 0.0026966, O6fitf -> -0.0238006, O8fitf -> 0.0228069, O10fitf -> -0.0236579, 
  O12fitf -> 0.0181667, O14fitf -> -0.00790285}}; alsmtauFOO14[5]=0.308623;
*)
(*   alpha_s(mtau^2)=0.312824+{^{+0.001204}}_{-0.004201}}_{Nt}; 
      10^3 O4=  -0.37+{^{-1.72}_{+3.07}}_{Nt};  
      10^3 O6= -19.71+{^{-0.13}_{-4.09}}_{Nt};  
      10^3 O8= +20.34+{^{+0.68}_{-2.47}}_{Nt}; 
      10^3 O10=-21.43+{^{-0.95}_{-2.23}}_{Nt}; 
      10^3 O12=+16.45+{^{+0.90}_{+1.72}}_{Nt}; 
      10^3 O14= -7.11+{^{-0.47}_{-0.79}}_{Nt}; 
    *)
    (* 24 Sept. 2022; checked 26 Sept. 2022 *)
    
 (* Note: Nt=8^{+2}_{-3} *)
 (*--------------*)
    (* Variation due to chnage of truncation (from O14 to O16 in OPE):
alO16fitPFOMinnum[8]={1.7125*10^-11, {aALf -> 0.100336, O4fitf -> -0.000256019, O6fitf -> -0.020145, O8fitf -> 0.0209145, O10fitf -> -0.0225467, O12fitf -> 0.0183289, O14fitf -> -0.00949344, O16fitf -> 0.00170014}}; alsmtauFOO16[8]=0.312386; *)
(* The central value:
alO14fitPFOMinnum[8]={2.10619*10^-11, {aALf -> 0.100478, O4fitf -> -0.000369672, O6fitf -> -0.0197131, O8fitf -> 0.0203391, O10fitf -> -0.0214331, O12fitf -> 0.0164522, O14fitf -> -0.00710677}}; alsmtauFOO14[8]=0.312824; *)

(*-->  alsmtauFOO16[8]-alsmtauFOO14[8]=-0.000438;
 10^3 (O4fitf-O4FOO14[8])=+0.11 =app +0.1;
 10^3 (O6fitf-O6FOO14[8])=-0.43 =app -0.4;
 10^3 (O8fitf-O8FOO14[8])=+0.58 =app +0.6;
 10^3 (O10fitf-O10FOO14[8])=-1.11 =app -1.1;
 10^3 (O12fitf-O12FOO14[8])=+1.88 =app +1.9;
 10^3 (O14fitf-O14FOO14[8])=-2.39 =app -2.4; *) 

    (*   alpha_s(mtau^2)=0.312824-+0.000438_{O16};
      10^3 O4=  -0.37 +-0.11_{O16};  
      10^3 O6= -19.71 -+0.43_{O16};  
      10^3 O8= +20.34 +-0.58_{O16};
      10^3 O10=-21.43 -+1.11_{O16};
      10^3 O12=+16.45 +-1.88_{O16};
      10^3 O14= -7.11 -+2.39_{O16};
    *)
    (* 28 Sept. 2022 *)
    
 (*-----------------------------------------*)
    (* Altogether, the FO global fit approach with (p3ka1ka2&O6ka2&O14; FO; Nt=N0 (=7)); in D=6 term a(Q^2)^(1-ka2) runs according to one-loop RGE from a(sigmax) to a(Q^2=sigmax Exp[I fi]) ]: *)
(* FO
 Nt(=N0)=8;  (kap=1+{^{+1}_{-1/2}}; d4=275+-63; Nt=8^{+2}_{-3}):
 
 alpha_s(mtau^2)=0.312824+-0.005269{exp}+{^{-0.001544}_{+0.003843}}_{kap}+{^{-0.007132}}_{+0.003959}}_{d4}+{^{+0.001204}}_{-0.004201}}_{Nt}-+0.000438_{O16}=0.312824+{^{+0.007736}}_{-0.009943}}
= 0.3128+{^{+0.0078}}_{-0.0099}} =app 0.313+{^{+0.008}}_{-0.010}}
10^3 O4=  -0.37 +-0.83_{exp}+{^{+1.53}_{-2.61}}_{kap} +{^{+0.43}_{-0.14}}_{d4}+{^{-1.72}_{+3.07}}_{Nt}+-0.11_{O16}
=  -0.37+{^{+3.56}_{-3.24}} =app  -0.4+{^{+3.6}_{-3.2}}
10^3 O6= -19.71 +-1.82_{exp}+{^{-1.88}_{+0.66}}_{kap}+{^{-0.53}_{+0.01}}_{d4}+{^{-0.13}_{-4.09}}_{Nt}-+0.43_{O16}
= -19.71 +{^{+1.99}_{-4.90}}  =app -19.7 +{^{+2.0}_{-4.9}} 
10^3 O8= +20.34 +-1.35_{exp}+{^{+1.10}_{+0.58}}_{kap}+{^{+0.18}_{+0.14}}_{d4}+{^{+0.68}_{-2.47}}_{Nt}+-0.58_{O16}
= +20.34+{^{+1.97}_{-2.87}} =app  +20.3+{^{+2.0}_{-2.9}}
10^3 O10=-21.43 +-1.25_{exp}+{^{-0.97}_{-1.16}}_{kap}+{^{-0.13}_{-0.20}}_{d4}+{^{-0.95}_{-2.23}}_{Nt}-+1.11_{O16}
= -21.43+{^{+1.67}_{-3.02}} =app -21.4+{^{+1.7}_{-3.0}}
10^3 O12=+16.45 +-0.69_{exp}+{^{+0.73}_{+1.24}}_{kap}+{^{+0.07}_{+0.20}}_{d4}+{^{+0.90}_{+1.72}}_{Nt}+-1.88_{O16}
= +16.45+{^{+3.38}_{+2.00}} =app  +16.5+{^{+3.4}_{+2.0}} 
10^3 O14= -7.11 +-9.12_{exp}+{^{-0.33}_{-0.69}}_{kap}+{^{-0.02}_{-0.11}}_{d4}+{^{-0.47}_{-0.79}}_{Nt}-+2.39_{O16}
= -7.11+{^{+9.43}_{-9.49}} =app -7.1+{^{+9.4}_{-9.5}}
chi^2=2.10619*10^-11;
*)

(* 24 Sept. 2022; checked alpha_s(mtau^2) on 26 Sept. 2022; added _{O16} uncertainties (in quadrature) in the analysis on 28 Sept. 2022 *)

(* End of FO global fit results, the case 'p3ka1ka2 & O6ka2 & O14', a(Q^2) runs in D=6 OPE term by one-loop RGE from a(sigmax)=aALfFO[Nt] to a(sigmax Exp[I fi])  *)
 (*------------------------------------------------------------------*)
   (*------------------------------------------------------------------*)
  (*------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(*---------------------------------------------------------------------*)
(* 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.1005565;
(* 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=kappa sigmax Exp[I (jj-1)Pi/Nn])^nn, such that aCI(kappa 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=kappa sigmax Exp[I (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 (13 Sept. 2022) (O4-O14 approach; O6 with nonzero anom. dim): *)    
 ReBPPVNumO14[M2f_,Nn_,O4_,O6_,O8_,O10_,O12_,O14_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPPVD0[M2f,Nn]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/3!) (1/M2f^4) (1+ 2*3 (M2f/sigmax) + 3*2 (M2f/sigmax)^2) + 2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2)+ 2 Pi^2 O12  (1/5!) (1/M2f^6) (1+ 2*5 (M2f/sigmax) + 5*4 (M2f/sigmax)^2)+ 2 Pi^2 O14  (1/6!) (1/M2f^7) (1+ 2*6 (M2f/sigmax) + 6*5 (M2f/sigmax)^2)+ (6 Pi^2 O6/sigmax^3) (1/2/Pi) NIntegrate[ ((aMSext[(1/kappa), aAL]/(1 + I be0 aMSext[(1/kappa), aAL] fi))^(1-ka2))  Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];
 (* Note: aAL=a(kappa*sigmax); -> aMSext[(1/kappa), aAL] = a(sigmax) *)
 (* BPPVD0[M2f,Nn] has  the polynomial (correction) part delD[th,nn] truncated at a[kappa sigmax Exp[I th]]^nn *)
 ReBPPVNumO16[M2f_,Nn_,O4_,O6_,O8_,O10_,O12_,O14_,O16_]:=Re[((1- 2 M2f/sigmax) + 2 (M2f/sigmax)^2 (1 - Exp[-sigmax/M2f]))+BPPVD0[M2f,Nn]+ 2 Pi^2 O4 (1/M2f^2) ( 1 + 2 M2f/sigmax )  + 2 Pi^2 O8  (1/6) (1/M2f^4) (1+ 6 (M2f/sigmax) + 6 (M2f/sigmax)^2) + 2 Pi^2 O10  (1/4!) (1/M2f^5) (1+ 2*4 (M2f/sigmax) + 4*3 (M2f/sigmax)^2)+ 2 Pi^2 O12  (1/5!) (1/M2f^6) (1+ 2*5 (M2f/sigmax) + 5*4 (M2f/sigmax)^2)+ 2 Pi^2 O14  (1/6!) (1/M2f^7) (1+ 2*6 (M2f/sigmax) + 6*5 (M2f/sigmax)^2)+ 2 Pi^2 O16  (1/7!) (1/M2f^8) (1+ 2*7 (M2f/sigmax) + 7*6 (M2f/sigmax)^2)+ (6 Pi^2/sigmax^3) (1/2/Pi) NIntegrate[ (O6 (aAL/(1 + I aAL be0 fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2f/sigmax) (1 + Exp[I fi]) + 2 (M2f/sigmax)^2) Exp[(sigmax/M2f) Exp[I fi]] - 2 (M2f/sigmax)^2 Exp[- sigmax/M2f] ),{fi, -Pi, Pi}]];

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

 O4init=-0.003; O6init=-0.001; 
 O8init=0.001; O6init = 0.003; O10init=0.005; O12init=0.004; O14init=0.005; O16init=0.004; O18init=0.005;  (* aAL=a(kappa*sigmax) is fixed  above *)
O14fitPPVMinNum[Nn_]:= FindMinimum[Sum[(ReBPPVNumO14[M2al[jj],Nn,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf,O14fitf] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init},{O14fitf, O14init}}];
O16fitPPVMinNum[Nn_]:= FindMinimum[Sum[(ReBPPVNumO16[M2al[jj],Nn,O4fitf,O6fitf,O8fitf,O10fitf,O12fitf,O14fitf,O16fitf] - ReBPexpBinComb[M2al[jj]])^2/delReBPexpBinComb[M2al[jj]]^2, {jj, 1, Nd}], {{O4fitf, O4init},{O6fitf, O6init},{O8fitf, O8init},{O10fitf, O10init},{O12fitf, O12init},{O14fitf, O14init},{O16fitf, O16init}}];

 (* O14fitPPVMinnum[nn]=O14fitPPVMinNum[nn] (nn=3,...,12) *)
 (* 13 Sept. 2022 *)
 
 (*--------------------------*)
(* NEW BL results, PPV approach with (O61,O62): *)
  (* Summary, for kappa=1 and NchP=10, PV(+CI) approach: *)

 (* NchP=2 (Nd=9): *)
 (*
  aALfPVO14[3]= 0.10167; O14fitPPVMinnum[3]={2.77718*10^-11, {O4fitf -> 0.000653185, O6fitf -> -0.0205304, O8fitf -> 0.0217305, O10fitf -> -0.0227294, O12fitf -> 0.0175581, O14fitf -> -0.00764371}};
  aALfPVO14[4]= 0.100662; O14fitPPVMinnum[4]={2.3352569145551875`*^-11, {O4fitf -> 0.0005814344245981375`, O6fitf -> -0.020754671067052`, O8fitf -> 0.021720909822340432`, O10fitf -> -0.022748736707124986`, O12fitf -> 0.017563080013372416`, O14fitf -> -0.007647005573088679`}};
  aALfPVO14[5]= 0.100611;  O14fitPPVMinnum[5]={2.3421814204710233`*^-11, {O4fitf -> 0.0005709890111093192`, O6fitf -> -0.020766581873468117`, O8fitf -> 0.02171746881298613`, O10fitf -> -0.02274863874543183`, O12fitf -> 0.017562795360977743`, O14fitf -> -0.007647652433972282`}};
  aALfPVO14[6]= 0.100545;  O14fitPPVMinnum[6]={2.3325367089664714`*^-11, {O4fitf -> 0.0005430106481995105`, O6fitf -> -0.020788350521124533`, O8fitf -> 0.02170232473307494`, O10fitf -> -0.022740718108200865`, O12fitf -> 0.017550974427852143`, O14fitf -> -0.007641893307011147`}};
  aALfPVO14[7]= 0.1005565; O14fitPPVMinnum[7]={2.330885502093159`*^-11, {O4fitf -> 0.0005512976262463319`, O6fitf -> -0.020779455039514425`, O8fitf -> 0.021707013617120914`, O10fitf -> -0.022741362567020793`, O12fitf -> 0.017553783045483943`, O14fitf -> -0.007642811311534027`}};
  aALfPVO14[8]= 0.100544; O14fitPPVMinnum[8]={2.3313634863924233`*^-11, {O4fitf -> 0.0005375962813667126`, O6fitf -> -0.020803317008168145`, O8fitf -> 0.0216988239456385`, O10fitf -> -0.022747253999690315`, O12fitf -> 0.017552943637182964`, O14fitf -> -0.007645520685722281`}};
  aALfPVO14[9]= 0.100550; O14fitPPVMinnum[9]={2.3317286039498198`*^-11, {O4fitf -> 0.0005485851664904719`, O6fitf -> -0.02077657617747881`, O8fitf -> 0.021709524573795356`, O10fitf -> -0.022737786581531413`, O12fitf -> 0.017555035045615695`, O14fitf -> -0.007641229786146724`}};
  aALfPVO14[10]=0.100547; O14fitPPVMinnum[10]={2.3333928719758416`*^-11, {O4fitf -> 0.0005372284177393538`, O6fitf -> -0.020813203525140444`, O8fitf -> 0.02168982070833387`, O10fitf -> -0.02275486198140519`, O12fitf -> 0.017547423399214713`, O14fitf -> -0.0076495516244590766`}};
  aALfPVO14[11]=0.100549; O14fitPPVMinnum[11]={2.3413287810995115`*^-11, {O4fitf -> 0.0005467261674244906`, O6fitf -> -0.02076561589055176`, O8fitf -> 0.021716905460583424`, O10fitf -> -0.022718917817195212`, O12fitf -> 0.017555432640357253`, O14fitf -> -0.007625818062543652`}};
  aALfPVO14[12]=0.100548; O14fitPPVMinnum[12]={2.342057335122187`*^-11, {O4fitf -> 0.0005417181638602849`, O6fitf -> -0.0208033699580834`, O8fitf -> 0.021691228388978134`, O10fitf -> -0.022760370097479533`, O12fitf -> 0.017545836013425253`, O14fitf -> -0.007659567267569674`}};
  Clear[Nn]; Do[Module[{}, xi2PVO14[Nn]=O14fitPPVMinnum[Nn][[1]]; alsmtauPVO14[Nn]=aMSext[mtau^2/(kappa sigmax),aALfPVO14[Nn]] Pi;  O4PVO14[Nn]=O4fitf/.O14fitPPVMinnum[Nn][[2]]; O6PVO14[Nn]=O6fitf/.O14fitPPVMinnum[Nn][[2]]; O6PVO14[Nn]=O6fitf/.O14fitPPVMinnum[Nn][[2]];O8PVO14[Nn]=O8fitf/.O14fitPPVMinnum[Nn][[2]];O10PVO14[Nn]=O10fitf/.O14fitPPVMinnum[Nn][[2]];O12PVO14[Nn]=O12fitf/.O14fitPPVMinnum[Nn][[2]];O14PVO14[Nn]=O14fitf/.O14fitPPVMinnum[Nn][[2]]],{Nn,3,12}];
  *)
  (* 20 Sept. 2022 *)

(* end of BL results (O6-O14) approach in PPV *) 
(*---------------------------------------------------------------------_*)
  (*---------------------------------------------------------------------_*)
 (* moments a^{(2,n}}, PV, at fixed truncation order and for aAL=(a(kappa*sigmax))=aALfPV[Nt]: *)			     

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

 (* Note, have to adjust aAL to aAL=aALfPV[NN] at each NN and calculate at that NN first the DBorfast etc. *)

 (*
  Tab20PVavar(*=Table[{NNt,Re[aaPV20[NNt,O4PVO14[NNt],O6PVO14[NNt]]]},{NNt,3,12}]={{3, 0.230683}, {4, 0.230683}, {5, 0.230683}, {6, 0.230683}, {7, 0.230683}, {8, 0.230683}, {9, 0.230683}, {10, 0.230683}, {11, 0.230683}, {12, 0.230683}}*)={{3,0.2306828972734235`},{4,0.2306826214758988`},{5,0.23068256943829674`},{6,0.23068262039024673`},{7,0.23068260842583804`},{8,0.23068257529919403`},{9,0.23068258640347233`},{10,0.23068262986090218`},{11,0.2306826337326977`},{12,0.2306826490285464`}}; (* Use: NumberForm[Tab20PVavar, 8] to better see the difference;  see also: delTab20PVavar= Table[{Nn, Tab20PVavar[[Nn - 2]][[2]] - Tab20PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; this suggests that the lowest Nt where local minimal sensitivity occurs is Nt=7; *)

  Tab21PVavar(*=Table[{NNt,Re[aaPV21[NNt,O6PVO14[NNt],O8PVO14[NNt]]]},{NNt,3,12}]={{3, 0.231983}, {4, 0.231983}, {5, 0.231982}, {6, 0.231983}, {7, 0.231982}, {8, 0.231982}, {9, 0.231982}, {10, 0.231983}, {11, 0.231983}, {12, 0.231983}}*)={{3,0.23198329934802017`},{4,0.23198253407636477`},{5,0.23198238082626363`},{6,0.23198251166709774`},{7,0.23198248073249841`},{8,0.23198238356464182`},{9,0.2319824175416385`},{10,0.2319825357182307`},{11,0.23198255809795854`},{12,0.23198258232329302`}}; (*  Use: NumberForm[Tab21PVavar, 8] to better see the difference; see also: delTab21PVavar= Table[{Nn, Tab21PVavar[[Nn - 2]][[2]] - Tab21PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; this suggests that the lowest Nt where local minimal sensitivity occurs is Nt=7; *)
 
  Tab22PVavar(*=Table[{NNt,Re[aaPV22[NNt,O6PVO14[NNt],O8PVO14[NNt],O10PVO14[NNt]]]},{NNt,3,12}]={{3, 0.203606}, {4, 0.203604}, {5, 0.203604}, {6, 0.203604}, {7, 0.203604}, {8, 0.203604}, {9, 0.203604}, {10, 0.203604}, {11, 0.203604}, {12, 0.203604}}*)={{3,0.20360618820664636`},{4,0.20360444603020778`},{5,0.20360408146611003`},{6,0.20360433275239287`},{7,0.20360427256113423`},{8,0.20360404279095476`},{9,0.203604127263543`},{10,0.20360437870342254`},{11,0.2036044638291934`},{12,0.2036044602774386`}};(*  Use: NumberForm[Tab22PVavar, 8] to better see the difference; see also: delTab22PVavar= Table[{Nn, Tab22PVavar[[Nn - 2]][[2]] - Tab22PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];  this suggests that the lowest Nt where local minimal sensitivity occurs is Nt=7; *)

  Tab23PVavar(*=Table[{NNt,Re[aaPV23[NNt,O6PVO14[NNt],O10PVO14[NNt],O12PVO14[NNt]]]},{NNt,3,12}]={{3, 0.182799}, {4, 0.182796}, {5, 0.182795}, {6, 0.182796}, {7, 0.182795}, {8, 0.182795}, {9, 0.182795}, {10, 0.182796}, {11, 0.182796}, {12, 0.182796}}*)={{3,0.18279928786747707`},{4,0.18279593540866654`},{5,0.18279519373356107`},{6,0.18279555254004573`},{7,0.18279546404726413`},{8,0.18279500071682928`},{9,0.18279518059734784`},{10,0.18279562006626537`},{11,0.1827958654868584`},{12,0.18279571686904056`}}; (*  Use: NumberForm[Tab23PVavar, 8] to better see the difference; see also:  delTab23PVavar= Table[{Nn, Tab23PVavar[[Nn - 2]][[2]] - Tab23PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}]; e[{Nn, Tab22PVavar[[Nn - 2]][[2]] - Tab22PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];  this suggests that the Nt where local minimal sensitivity occurs is Nt=7; *)

  Tab24PVavar(*=Table[{NNt,Re[aaPV24[NNt,O6PVO14[NNt],O12PVO14[NNt],O14PVO14[NNt]]]},{NNt,3,12}]={{3, 0.171264}, {4, 0.171259}, {5, 0.171258}, {6, 0.171258}, {7, 0.171258}, {8, 0.171257}, {9, 0.171257}, {10, 0.171258}, {11, 0.171259}, {12, 0.171258}}*)={{3,0.17126436338692344`},{4,0.17125897545417862`},{5,0.17125767692321792`},{6,0.17125797645944607`},{7,0.1712578944801864`},{8,0.17125709401168734`},{9,0.17125742249673845`},{10,0.1712580468446355`},{11,0.17125860546496494`},{12,0.17125808336226075`}}; (*  Use: NumberForm[Tab24PVavar, 8] to better see the difference; see also: delTab24PVavar= Table[{Nn, Tab24PVavar[[Nn - 2]][[2]] - Tab24PVavar[[Nn - 3]][[2]]}, {Nn, 4, 12}];  this suggests that the Nt where local minimal sensitivity occurs is Nt=7; *)
  *)
  (* 24 Sept. 2022 *)


  (* Note, a24 PV method reproduces very well the experimental medium value of a24 (0.171290); this is in contrast to FO approach where a24 experimental medium value is NOT well reproduced *) 
(* Therefore, I will take Nt=7 in the PV approach (in the considered V-channel, where OPE is truncated at O14) *)	
 (*
 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}]; 
 (* 30 July 2022 *)
  *) 

 (* aexpBinComb[2, 0]+-delaexpBinComb[2, 0]=0.230681 +-0.088331; *)
 (* aexpBinComb[2, 1]+-delaexpBinComb[2, 1]=0.231976 +-0.085489; *)
 (* aexpBinComb[2, 2]+-delaexpBinComb[2, 2]=0.203590 +- 0.086436; *)
 (* aexpBinComb[2, 3]+-delaexpBinComb[2, 3]=0.182778 +- 0.090541; *)
 (* aexpBinComb[2, 4]+-delaexpBinComb[2, 4]=0.171290 +- 0.096610; *)

  (* For PV approach I came down to here on 20 Sept. 2022 *)
 (* Nt=7 is the best Nt *)
  (*-------------------*)   

(*---------------------------------------*)
 (* the central value of the PV global fit result, at Nt=5, is:
  aAL=aALfPVO14[7]=0.1005565; O14fitPPVMinnum[7]={2.330885502093159`*^-11, {O4fitf -> 0.0005512976262463319`, O6fitf -> -0.020779455039514425`, O8fitf -> 0.021707013617120914`, O10fitf -> -0.022741362567020793`, O12fitf -> 0.017553783045483943`, O14fitf -> -0.007642811311534027`}}; alpha_s(mtau^2)(=alsmtauPV[7])=0.313064; *)
  (*------------------------*) 
 (* Now the estimate of EXPERIMENTAL uncertainty in this global approach with double-pinched PV Borel-Laplace (with Nt=5): *)
 (* See Appendix  PRD84, 113006 (Boito et al., 2011): *)
 
(* Now the M^2_{al}: all Nd=9 values (three rays in the complex M^2-plane) as one array, as M2al[jj] (jj=1,...,9): *)
 (* Nm=NchP+1; Nd=3 Nm;  
M2al0[ii_]:=0.90+(ii-1)*0.60/(Nm-1);
Do[M2al[jj]=M2al0[jj],{jj,1,Nm}]; Do[M2al[jj+Nm]=M2al0[jj] Exp[I Pi/6],{jj,1,Nm}]; Do[M2al[jj+2 Nm]=M2al0[jj] Exp[I Pi/4],{jj,1,Nm}];
  *)

 N0=7; (* N0=Nt *)  Np=7; (* Np is the number of parameters: a(sigmax) and (6) condensates *)
epp1=0.0001; (* aALfPVO14[N0]=0.1005565; *)
 (* First derivatives of ReBPFONum[M2] with respect to p1=aAL, p2=O4, p3=O6, p4=O8, p5=O10, p6=O12, p7=O14  [note, we must interpret in this program: O2n = 2 <O_{2 n}>_V]: *)
 (* For variation of ReBPPVNumO14[M2al[jj], N0, O4PV[N0],O6PV[N0],O8PV[N0],O10PV[N0],O12PV[N0],O14PV[N0]] when aAL is varied 0.1% up and one 0.1% down: *)
 (* aAL=aALfPVO14[N0](1+epp1)=0.10056655565: run the program for PV case with this aAL (and kappa=1)  (between 'Double-pinched' and at least until 'Now search for')  and evaluate: *)
 (*  aAL=aALfPVO14[N0](1+epp1)=0.10056655565: tabReBPPVUP= Table[{M2al[jj], ReBPPVNumO14[M2al[jj], N0, 0, O6PVO14[N0], 0, 0,0,0]}, {jj, 1,Nd}] *) tabReBPPVUP={{0.9`, 0.388481836895965`}, {1.2`, 0.46258523829558984`}, {1.5`, 0.4530688361680313`}, {0.7794228634059949` + 0.44999999999999996` I, 0.6674543169318474`}, {1.0392304845413265` + 0.5999999999999999` I, 0.5827046091874306`}, {1.299038105676658` + 0.7499999999999999` I, 0.5103250774876955`}, {0.6363961030678928` + 0.6363961030678927` I, 0.8877222705721356`}, {0.848528137423857` + 0.8485281374238569` I, 0.6746681296676096`}, {1.0606601717798214` + 1.0606601717798212` I, 0.5492643820136833`}};
  (*  aAL=aALfPVO14[N0](1-epp1)=0.10054644435:  run the program for PV case with this aAL (and kappa=1)  (between 'Double-pinched' and 'Now search for') and evaluate: *)
 (*  aAL=aALfPVO14[N0](1-epp1)=0.10054644435: tabReBPPVDO= Table[{M2al[jj], ReBPPVNumO14[M2al[jj], N0, 0, O6PVO14[N0], 0, 0,0,0]}, {jj, 1,Nd}] *) tabReBPPVDO={{0.9`, 0.3884118798733953`}, {1.2`, 0.46251994281396663`}, {1.5`, 0.4530115371434435`}, {0.7794228634059949` + 0.44999999999999996` I, 0.6673618514169609`}, {1.0392304845413265` + 0.5999999999999999` I, 0.5826339746927316`}, {1.299038105676658` + 0.7499999999999999` I, 0.510267836712585`}, {0.6363961030678928` + 0.6363961030678927` I, 0.8876199541982691`}, {0.848528137423857` + 0.8485281374238569` I, 0.6745981807225716`}, {1.0606601717798214` + 1.0606601717798212` I, 0.5492108255431194`}};
 (* tabdReBPPVd[1] below is now the table of numerical derivative of  ReBPPVNumO10[M2] with respect to aAL around aAL=aALfPVO14[N0]=0.1005565: *)
 tabdReBPPVd[1] = Table[{M2al[ii],  (tabReBPPVUP[[ii]][[2]]-tabReBPPVDO[[ii]][[2]])/(2 epp1 aALfPVO14[N0])}, {ii, 1, Nd}];

(* Now the derivatives with respect to the condensates: *)
 dReBPPVd[2][M2_]:=2 Pi^2 Re[ (1/M2^2) ( 1 + 2 M2/sigmax )];
dReBPPVd[3][M2_]:=(6 Pi^2/sigmax^3) (1/2/Pi) Re[NIntegrate[ ((aALfPVO14[N0]/(1 + I be0 aALfPVO14[N0] fi))^(1-ka2) ) Exp[-I 3 fi] ( ((1 + Exp[I fi])^2 - 2 (M2/sigmax) (1 + Exp[I fi]) + 2 (M2/sigmax)^2) Exp[(sigmax/M2) Exp[I fi]] - 2 (M2/sigmax)^2 Exp[- sigmax/M2] ),{fi, -Pi, Pi}]]; 
dReBPPVd[4][M2_]:=2 Pi^2 (1/3!) Re[(1/M2^4) (1+ 2*3 (M2/sigmax) + 3*2 (M2/sigmax)^2)]; 
dReBPPVd[5][M2_]:=2 Pi^2 (1/4!) Re[(1/M2^5) (1+ 2*4 (M2/sigmax) + 4*3 (M2/sigmax)^2)];
dReBPPVd[6][M2_]:=2 Pi^2 (1/5!) Re[(1/M2^6) (1+ 2*5 (M2/sigmax) + 5*4 (M2/sigmax)^2)];
dReBPPVd[7][M2_]:=2 Pi^2 (1/6!) Re[(1/M2^7) (1+ 2*6 (M2/sigmax) + 6*5 (M2/sigmax)^2)];
 (* Now Tables of derivatives with respect to the condensates, for all Mal[jj]: *)  
 tabdReBPPVd[2] = Table[{M2al[ii], dReBPPVd[2][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPPVd[3] = Table[{M2al[ii], dReBPPVd[3][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPPVd[4] = Table[{M2al[ii], dReBPPVd[4][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPPVd[5] = Table[{M2al[ii], dReBPPVd[5][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPPVd[6] = Table[{M2al[ii], dReBPPVd[6][M2al[ii]]}, {ii, 1, Nd}];
 tabdReBPPVd[7] = Table[{M2al[ii], dReBPPVd[7][M2al[ii]]}, {ii, 1, Nd}];  
 tabCovReBPexpBinComb = Table[{M2al[ii],M2al[jj],CovReBPexpBinComb[M2al[ii],M2al[jj]]},{ii,1,Nd},{jj,1,Nd}];
mCovReBPexpBinComb=Table[tabCovReBPexpBinComb[[i, j]][[3]],{i,1,Nd},{j,1,Nd}];
(* This is the NdxNd matrix form of the covariances <Re[B(M^2_{al})] Re[B(M^2_{be})]; it shows very strong correlations of various Re[B(M^2_{al})] due to experimental uncertainties of data Omega_V(sig); all terms of mCovReBPexpBinComb are (4-9)x 10^(-6) *)
Np=7;
 (* Now calculate the NpxNp mtrx AA(=A), for the PV case, using tables obtained above: *)
 tabAAPV=Table[{ik, jk, Sum[(tabdReBPPVd[ik][[jj]][[2]] tabdReBPPVd[jk][[jj]][[2]])/tabCovReBPexpBinComb[[jj]][[jj]][[3]], {jj, 1, Nd}]}, {ik, 1, Np}, {jk, 1, Np}];
 Do[AAPVa[ii,jj]=tabAAPV[[ii]][[jj]][[3]],{ii,1,Np},{jj,1,Np}];  (* This is NpxNp matrix A (or: AA), calculated more efficiently *)
AAPVm=Table[AAPVa[ik,jk],{ik,1,Np},{jk,1,Np}]; (* all these terms are =app 10^4-10^8 *)
AAPVmI=Inverse[AAPVm];   (* AAPVmI = (AAPVm)^{-1}, inverse NpxNp matrix of AAm *)
(* WARNING APPEARS: Result for Inverse of badly conditioned matrix; then the results are completely unreliable and numerically unstable *)

(* Now the NpxNp covariance matrix of the fitted parameters pj: *) 
(* p1=aAL, p2=<O4>, p3=<O6>, p4=<O8>, p5=<O10>, p6=<O12>, p7=<O14>; where <O_{2n}> means: 2 <O_{2n}>_V  *)
(* PP[i,j]= < (delta p_i) (delta p_j) >: *)
  (* We calulate PP[i,j] covariance matrix of the extracted parameters, using the tables obtained above: *) 
 PPtabPV[jp_,kp_]:=Sum[AAPVmI[[jp,jpp]] AAPVmI[[kp,kpp]] tabdReBPPVd[jpp][[Al]][[2]] tabdReBPPVd[kpp][[Be]][[2]] tabCovReBPexpBinComb[[Al]][[Be]][[3]]/(tabCovReBPexpBinComb[[Al]][[Al]][[3]]  tabCovReBPexpBinComb[[Be]][[Be]][[3]]),{jpp,1,Np},{kpp,1,Np},{Be,1,Nd},{Al,1,Nd}];

(* I then get, for epp1=0.001: 
 Table[{nn, Sqrt[PPtabPV[nn, nn]]}, {nn, 1, 7}]={{1, 0.00243024}, {2, 0. + 0.000813379 I}, {3, 0.000778722}, {4, 0. + 0.00255528 I}, {5, 0.00114604}, {6, 0.0079865}, {7, 0.00406901}}; *)
(* I then get, for epp1=0.0001: 
 Table[{nn, Sqrt[PPtabPV[nn, nn]]}, {nn, 1, 7}]={{1, 0.00238761}, {2, 0.000778572}, {3, 0.00574747}, {4, 0. + 0.00388314 I}, {5, 0.00616531}, {6, 0.0129453}, {7, 0. + 0.00274455 I}}; *)
 (* Np=7: numerical instabilities; some Sqrt[PPtabPV[nn, nn]] are complex numbers *)
 (*----------------------------------------*)
 (* Therefore, we will FIX O14 (i.e., NOT vary) to its central value (O14PVO14[N0]=-0.00764281) and NOT vary it, in order to reduce the matix to 6x6 and thus avoid numerical instabilities: *)
 (* Using the above results, with epp1=0.0001: *)
 Np=6;  
 tabAAPV=Table[{ik, jk, Sum[(tabdReBPPVd[ik][[jj]][[2]] tabdReBPPVd[jk][[jj]][[2]])/tabCovReBPexpBinComb[[jj]][[jj]][[3]], {jj, 1, Nd}]}, {ik, 1, Np}, {jk, 1, Np}];
 Do[AAPVa[ii,jj]=tabAAPV[[ii]][[jj]][[3]],{ii,1,Np},{jj,1,Np}];  (* This is NpxNp matrix A (or: AA), calculated more efficiently *)
AAPVm=Table[AAPVa[ik,jk],{ik,1,Np},{jk,1,Np}]; (* all these terms are =app 10^4-10^8 *)
AAPVmI=Inverse[AAPVm];   (* AAPVmI = (AAPVm)^{-1}, inverse NpxNp matrix of AAm *)
 (* Now: NO WARNING APPEARS: Result for Inverse of badly conditioned matrix; *)
(* p1=aAL, p2=<O4>, p3=<O6>, p4=<O8>, p5=<O10>, p6=<O12>, where <O_{2n}> means: 2 <O_{2n}>_V  *)
(* PP[i,j]= < (delta p_i) (delta p_j) >: *)
  (* We calulate PP[i,j] covariance matrix of the extracted parameters, using the tables obtained above: *) 
 PPtabPV[jp_,kp_]:=Sum[AAPVmI[[jp,jpp]] AAPVmI[[kp,kpp]] tabdReBPPVd[jpp][[Al]][[2]] tabdReBPPVd[kpp][[Be]][[2]] tabCovReBPexpBinComb[[Al]][[Be]][[3]]/(tabCovReBPexpBinComb[[Al]][[Al]][[3]]  tabCovReBPexpBinComb[[Be]][[Be]][[3]]),{jpp,1,Np},{kpp,1,Np},{Be,1,Nd},{Al,1,Nd}];
 (* with epp1=0.0001: *)
 (* Table[{nn, Sqrt[PPtabPV[nn, nn]]}, {nn, 1, Np}]={{1, 0.00173226}, {2, 0.000635052}, {3, 0.00171494}, {4, 0.00123742}, {5, 0.00118289}, {6, 0.00065679}}; *)
 (* So I will use this; del(a(sigmax))_{exp}=Sqrt[PPtabPV[1,1]]; del(O4)_{exp}=Sqrt[PPtabPV[2,2]]; etc.; del(O12)_{exp}=Sqrt[PPtabPV[6,6]];

 (*------------------------------------------------------------------*)
  (* For del(O14)_{exp} we take now: +-delP[7] (probably an overtestimate): *)
   Np=7; (* calculate for Np=7 the mtrx AAFOmI, with epp1=0.0001; *)
 delP[jp_]:=Sum[AAPVmI[[jp,jk]] tabdReBPPVd[jk][[Al]][[2]]/(Sqrt[CovReBPexpBinComb[M2al[Al],M2al[Al]]]),{jk,1,Np},{Al,1,Nd}];
  delP[7]=0.00893699;  (* note: delP[1]=-0.0040112 *)
  (*--------------------------------------------------------------------*)
  (* Thus: (N0=7; Np=6:) Table[{nn, Sqrt[PPtabPV[nn, nn]]}, {nn, 1, Np}]={{1, 0.00173226}, {2, 0.000635052}, {3, 0.00171494}, {4, 0.00123742}, {5, 0.00118289}, {6, 0.00065679}}; delP7]=0.00893699; --> *)
 (* a(sigmax)=aALfPVO14[N0]+-Sqrt[PPtabPV[1,1]]=0.100557+-0.00173226_{exp}; O4=O4PVO14[N0]+-Sqrt[PPtabPV[2,2]]; ...; O12=O12PVO14[N0]+-Sqrt[PPtabPV[6,6]]; O14=O14PVO14[N0]+-delP[7] --> *)
(*
 alpha_s(mtau^2)=0.313064+-0.005332_{exp};
 10^3 O4=  +0.55+- 0.64_{exp};
 10^3 O6= -20.78 +-1.71_{exp};
 10^3 O8= +21.71 +-1.24_{exp};
 10^3 O10=-22.74 +-1.18_{exp};
 10^3 O12=+17.554 +-0.66_{exp};
 10^3 O14= -7.64 +-8.94_{exp};
 *)
 (* 23 Sept. 2022; checked 26 Sept. 2022 *)
(*--------------------------------------------------*)
 (* Variation of kappa, PV global: *)
 (* kappa=2.0: ddtk[4, tkap]=936.536; N0=7: 
 aAL=0.084328; O14fitPPVMinN[N0]={2.43293*10^-11, {O4fitf -> 0.000435668, O6fitf -> -0.020889,  O8fitf -> 0.0216552, O10fitf -> -0.0227037, O12fitf -> 0.0174843, O14fitf -> -0.00758868}}; a(sigmax)=0.10047; alpha_s(mtau^2)=0.312799; *)
 (* kappa=1.0; ddtk[4,tkap]=275.; N0=7 (the central case of the global fit): 
  aAL=0.1005565; O14fitPPVMinnum[N0]={2.33089*10^-11, {O4fitf -> 0.000551298, O6fitf -> -0.0207795, O8fitf -> 0.021707, O10fitf -> -0.0227414, O12fitf -> 0.0175538, O14fitf -> -0.00764281}}; a(sigmax)=aALfPVO14[7]=0.1005565; alpha_s(mtau^2)(=alsmtauPVO14[7])=0.313064; *)
 (* kappa=0.5: ddtk[4, tkap]=-1.16848; N0=7:
  aAL=0.126354; (xi2=2.54993*10^-11); O14fitPPVMinN[N0]={2.54993*10^-11, {O4fitf -> 0.000579519, O6fitf -> -0.0204123, O8fitf -> 0.0214909, O10fitf -> -0.0221568, O12fitf -> 0.0168503, O14fitf -> -0.00710454}}; a(sigmax)=0.100651; alpha_s(mtau^2)=0.313356; *)
 (* checked on 26 Sept. 2022 *)
 (* kap=1^{+1}-{1/2} (Nt=N0(=7)):*) 
 (*
  alpha_s(mtau^2)=0.313064+{^{-0.000265}_{+0.000292}}_{kap};
  10^3 O4=  +0.55+{^{-0.11}_{+0.03}}_{kap};
  10^3 O6= -20.78+{^{-0.11}_{+0.37}}_{kap};
  10^3 O8= +21.71+{^{-0.05}_{-0.22}}_{kap};
  10^3 O10=-22.74+{^{+0.04}_{+0.58}}_{kap};
  10^3 O12=+17.55+{^{-0.07}_{-0.70}}_{kap};
  10^3 O14= -7.64+{^{+0.05}_{+0.54}}_{kap};
  *)		  		  
  (* 23 Sept. 2022; checked on 26 Sept. 2022 *) 
(*----------------*)
 (* Variation of d4, PV global (d4=227 +- 63; kappa=1): 
  ddtk[4,tkap]=227+63=338 (N0=7):
  aAL=0.098353; O14fitPPVMinnum[N0]={2.19131*10^-11, {O4fitf -> 0.000745299, O6fitf -> -0.017357, O8fitf -> 0.0196056, O10fitf -> -0.02105, O12fitf -> 0.0161647, O14fitf -> -0.00697184}}; alpha_s(mtau^2)=0.306281;
  Central case: ddtk[4,tkap]=275. (N0=7):
 aAL=0.1005565; O14fitPPVMinnum[N0]={2.33089*10^-11, {O4fitf -> 0.000551298, O6fitf -> -0.0207795, O8fitf -> 0.021707, O10fitf -> -0.0227414, O12fitf -> 0.0175538, O14fitf -> -0.00764281}};  alpha_s(mtau^2)(=alsmtauPVO14[7])=0.313064; 
  ddtk[4,tkap]=227-63=212 (N0=7):
  aAL=0.101320; O14fitPPVMinnum[N0]={2.54618*10^-11, {O4fitf -> 0.000726708, O6fitf -> -0.0213375, O8fitf -> 0.0219576, O10fitf -> -0.0230582, O12fitf -> 0.0178139, O14fitf -> -0.00777985}}; alpha_s(mtau^2)=0.315414;
  *)
  (* checked on 26 Sept. 2022 *) 
 (* d4=275 +- 63: *) 
(*
 alpha_s(mtau^2)=0.313064+{^{-0.006783}_{+0.002349}}_{d4};
 10^3 O4=  +0.55+{^{+0.20}_{+0.18}}_{d4};
 10^3 O6= -20.78+{^{+3.42}_{-0.56}}_{d4};
 10^3 O8= +21.71+{^{-2.10}_{+0.25}}_{d4};
 10^3 O10=-22.74+{^{+1.69}_{-0.32}}_{d4};
 10^3 O12=+17.55+{^{-1.39}_{+0.26}}_{d4};
 10^3 O14= -7.64+{^{+0.67}_{-0.14}}_{d4};
 *)
  (* 23 Sept. 2022; checked on 26 Sept. 2022 *)
   (*---------------------*)
  (* Variation of Nt (truncation order, Nt=Nt0+2-1, i.e., variation Nt=7,6,5,4; PV global with (O61, O62); 
   aALfPVO14[9]=0.10055; O14fitPPVMinnum[9]={2.33173*10^-11, {O4fitf -> 0.000548585, O6fitf -> -0.0207766, O8fitf -> 0.0217095, O10fitf -> -0.0227378, O12fitf -> 0.017555, O14fitf -> -0.00764123}}; alsmtauPVO14[9]=0.313044;
  aALfPVO14[8]=0.100544; O14fitPPVMinnum[8]={2.33136*10^-11, {O4fitf -> 0.000537596, O6fitf -> -0.0208033, O8fitf -> 0.0216988, O10fitf -> -0.0227473, O12fitf -> 0.0175529, O14fitf -> -0.00764552}}; alsmtauPVO14[8]=0.313026;
  aALfPVO14[7]=0.1005565; O14fitPPVMinnum[7]={2.33089*10^-11, {O4fitf -> 0.000551298, O6fitf -> -0.0207795, O8fitf -> 0.021707, O10fitf -> -0.0227414, O12fitf -> 0.0175538, O14fitf -> -0.00764281}}; alsmtauPVO14[7]=0.313064;
  aALfPVO14[6]=0.100545; O14fitPPVMinnum[6]={2.33254*10^-11, {O4fitf -> 0.000543011, O6fitf -> -0.0207884, O8fitf -> 0.0217023, O10fitf -> -0.0227407, O12fitf -> 0.017551, O14fitf -> -0.00764189}}; alsmtauPVO14[6]=0.313029;
  aALfPVO14[5]=0.100611; O14fitPPVMinnum[5]={2.34218*10^-11, {O4fitf -> 0.000570989, O6fitf -> -0.0207666, O8fitf -> 0.0217175, O10fitf -> -0.0227486, O12fitf -> 0.0175628, O14fitf -> -0.00764765}}; alsmtauPVO14[5]=0.313232;
   *)
 (* Nt=7^{+1}_{-1}, 7^{+2}_{-2}:*)
(* d4=275 +- 63: *) 
(*
 alpha_s(mtau^2)=0.313064+{^{-0.000038}_{+0.000168}}_{Nt};
 10^3 O4=  +0.55+{^{-0.00}_{+0.02}}_{Nt};
 10^3 O6= -20.78+{^{-0.00}_{+0.01}}_{Nt};
 10^3 O8= +21.71+{^{-0.00}_{+0.01}}_{Nt};
 10^3 O10=-22.74+{^{-0.00}_{-0.01}}_{Nt};
 10^3 O12=+17.55+{^{+0.00}_{+0.01}}_{Nt};
 10^3 O14= -7.64+{^{-0.00}_{-0.01}}_{Nt};
 *)
  (* 23 Sept. 2022; checked on 26 Sept. 2022 *)
   (*-------------------------*)
 (* Variation due to the renormalon ambiguity ('amb') to Adler: in ReBPPVNumO14we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn]+BPPVdelD0[M2f])  (N0=7) *)
 (*  aAL=0.100539; O14fitPPVMinNum[N0]={2.31433*10^-11, {O4fitf -> -0.000168643, O6fitf -> -0.0207792, O8fitf -> 0.0216949, O10fitf -> -0.0227116, O12fitf -> 0.0175373, O14fitf -> -0.00763357}}; alpha_s(mtau^2)=0.313010; *)
(*---*)
 (* The central value  (BBPVD0[M2f,Nn]) (N0=7) : 
  aALfPVO14[N0]=0.1005565; O14fitPPVMinnum[N0]={2.33089*10^-11, {O4fitf -> 0.000551298, O6fitf -> -0.0207795, O8fitf -> 0.021707, O10fitf -> -0.0227414, O12fitf -> 0.0175538, O14fitf -> -0.00764281}}; alsmtauPVO14[7]=0.313064; *)
  (*---*)
 (* When  in ReBPPVNumO14 we replace BBPVD0[M2f,Nn] by: (BPPVD0[M2f, Nn]-BPPVdelD0[M2f]) (N0=7) *)
 (*  aAL=0.100574; O14fitPPVMinNum[N0]={2.35378*10^-11, {O4fitf -> 0.00127302, O6fitf -> -0.0207798, O8fitf -> 0.0217193, O10fitf -> -0.0227715, O12fitf -> 0.0175706, O14fitf -> -0.00765228}}; alpha_s(mtau^2)=0.313118; *)
  
(*
 alpha_s(mtau^2)=0.313064+{^{-0.000054}_{+0.000054}}_{amb};
 10^3 O4=  +0.55+{^{-0.72}_{+0.72}}_{amb};
 10^3 O6= -20.78+{^{+0.00}_{-0.00}}_{amb};
 10^3 O8= +21.71+{^{-0.02}_{+0.01}}_{amb};
 10^3 O10=-22.74+{^{+0.03}_{-0.03}}_{amb};
 10^3 O12=+17.55+{^{-0.01}_{+0.02}}_{amb};
 10^3 O14= -7.64+{^{+0.01}_{-0.01}}_{amb};
 *)
  (* 24 Sept. 2022; checked on 26 Sept. 2022 *)
 (*--------------*)
    (* Variation due to chnage of truncation (from O14 to O16 in OPE):
aAL=0.100323; O16fitPPVMinNum[7]={1.26198*10^-11, {O4fitf -> 0.000707192, O6fitf -> -0.0214627, O8fitf -> 0.0226195, O10fitf -> -0.0245396, O12fitf -> 0.0206073, O14fitf -> -0.0115441, O16fitf -> 0.0027866}}; alsmtauPVO16[7]=alsmt[aAL]=0.312346; *)
(* The central value:
aAL=aALfPVO14[7]=0.1005565; O14fitPPVMinnum[7]={2.33089*10^-11, {O4fitf -> 0.000551298, O6fitf -> -0.0207795, O8fitf -> 0.021707, O10fitf -> -0.0227414, O12fitf -> 0.0175538, O14fitf -> -0.00764281}}; alsmtauPVO14[7]=0.313064; *)

(*-->  alsmtauPVO16[7]-alsmtauPVO14[7]=-0.000719;
 10^3 (O4fitf-O4PVO14[7])=+0.16 =app +0.2;
 10^3 (O6fitf-O6PVO14[7])=-0.68 =app -0.7;
 10^3 (O8fitf-O8PVO14[7])=+0.91 =app +0.9;
 10^3 (O10fitf-O10PVO14[7])=-1.80 =app -1.8;
 10^3 (O12fitf-O12PVO14[7])=+3.05 =app +3.1;
 10^3 (O14fitf-O14PVO14[7])=-3.90 =app -3.9; *) 

    (*   alpha_s(mtau^2)=0.313064-+0.000719_{O16};
      10^3 O4=  +0.55 +-0.16_{O16};  
      10^3 O6= -20.78 -+0.68_{O16};  
      10^3 O8= +21.71 +-0.91_{O16};
      10^3 O10=-22.74 -+1.80_{O16};
      10^3 O12=+17.55 +-3.05_{O16};
      10^3 O14= -7.64 -+3.90_{O16};
    *)
    (* 28 Sept. 2022 *)

  (*-----------------------------------------*)
(* Altogether, the PV(+CI) global fit approach with ('p3ka1ka2 & O6ka2 O14') V-channel (Peris et al data), Nt=N0(=7), and a(Q^2) at D=6 OPE term runs according to the one-loop RGE from a(sigmax)=aALfPV[Nt] to a(Q^2=sigmax Exp[I fi]): *)
  (* PV(+CI)
   Nt=N0(=7);  ( ddtk[4, Exp[tKK]]=275.);
   d4=275.+-63.; kap=1^{+1}_{-1/2};  Nt=7^{+1}_{-1} and &7^{+2}_{-2}:
   
 alpha_s(mtau^2)=0.313064+-0.005332_{exp}+{^{-0.000265}_{+0.000292}}_{kap}+{^{-0.006783}_{+0.002349}}_{d4}+{^{-0.000038}_{+0.000168}}_{Nt}+{^{-0.000054}_{+0.000054}}_{amb}-+0.000719_{O16} =0.313064+{^{+0.005880}_{-0.008662}}
   =0.3131+{^{+0.0058}_{-0.0087}}  =app 0.313+{^{+0.006}_{-0.009}};
 10^3 O4=   +0.55 +-0.64_{exp}+{^{-0.11}_{+0.03}}_{kap}+{^{+0.20}_{+0.18}}_{d4}+{^{-0.00}_{+0.02}}_{Nt}+{^{-0.72}_{+0.72}}_{amb}+-0.16_{O16}
   =  +0.55+{^{+0.99}_{-0.98}} =app  +0.6+-1.0;
   10^3 O6=  -20.78 +-1.71_{exp}+{^{-0.11}_{+0.37}}_{kap}+{^{+3.42}_{-0.56}}_{d4}+{^{-0.00}_{+0.01}}_{Nt}+{^{+0.00}_{-0.00}}_{amb}-+0.68_{O16}
   = -20.78+{^{+3.90}_{-1.92}}  =app -20.8+{^{+3.9}_{-1.9}};
   10^3 O8=  +21.71 +-1.24_{exp}+{^{-0.05}_{-0.22}}_{kap}+{^{-2.10}_{+0.25}}_{d4}+{^{-0.00}_{+0.01}}_{Nt}+{^{-0.02}_{+0.01}}_{amb}+-0.91_{O16}
   = +21.71+{^{+3.38}_{-3.98}}  =app +21.7+{^{+3.4}_{-4.0}};
   10^3 O10= -22.74 +-1.18_{exp}+{^{+0.04}_{+0.57}}_{kap}+{^{+1.69}_{-0.32}}_{d4}+{^{-0.00}_{-0.01}}_{Nt}+{^{+0.03}_{-0.03}}_{amb}-+1.80_{O16}
   = -22.74+{^{+2.80}_{-2.17}}  =app -22.7+{^{+2.8}_{-2.2}};
   10^3 O12=+17.554 +-0.66_{exp}+{^{-0.07}_{-0.70}}_{kap}+{^{-1.39}_{+0.26}}_{d4}+{^{-0.00}_{+0.01}}_{Nt}+{^{-0.01}_{+0.02}}_{amb}+-3.05_{O16}
   = +17.55+{^{+3.13}_{-3.49}} =app +17.6+{^{+3.1}_{-3.5}};
   10^3 O14=  -7.64 +-8.94_{exp}+{^{+0.05}_{+0.54}}_{kap}+{^{+0.67}_{-0.14}}_{d4}+{^{-0.00}_{-0.01}}_{Nt}+{^{+0.01}_{-0.01}}_amb}-+3.90_{O16}
=  -7.64+{^{+9.79}_{-9.75}}  =app  -7.6+{^{+9.8}_{-9.8}};
chi^2=2.33089*10^-11;
   *)
 
(* 24 Sept. 2022; checked on 26 Sept. 2022; added _{O16} uncertainties (in quadrature) in the analysis, on 28 Sept. 2022 *)
			  
(* End of PV(+CI) global fit results with ('p3ka1ka2' & 'O6ka2' & 'O14') PV approach, and a(Q^2) in D=6 OPE term 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/SumRpQCD.dir/RGErundown5l.m;  aaMSQ*Pi
 Then: aaMSQ*Pi is the value of alpha_s(mtau^2) *)
