(* 4l3dN.m: as 4l3dA00nf3rt0XXXal0YYYYNsh, but longer version *)
(* In comparison with older versions (prgs. 4l3dA00nf3rt0201al01185.m, etc.), it has several improvements: *)
(* 1. In this program in MSbar we run from Q^2=MZ^2 down to (2 mc)^2 with 5-loop RGE (not any more: 4-loop) and the corresponding 4-loop quark mass matching (not any more 3-loop). *)
(* 2. c2 and c3 of LMM are taken as exact, and not as rounded, values. *)
(* 3. The RSch-change from MSbar to LMM (at Q^2=(2 mc)^2, nf=3) is made with 5-loop MSbar beta function. *)
(* 4. 'N' in the program name stands for the new evaluation of rtau(D=0) - via my renormalon model, prg. Adl4l3dRenMod4A(sh).m *)
(*-------------------------------------------*)
 (* Note, BEFORE calling this program, we must specify the OUTER INPUTS: *)
 (* 1. alphaMSMZ (=alpha_s(MZ^2;MSbar)). *)
 (* 2. s0 *)
 (* 3. s1in (=s1s0=s1); in addition, also the initial approximate values (for iterative seach of the solution) are needed: f1sin=...; f2sin=...; f3sin=...; 
s2in=...; s3in =...;*)
 (*-------------------------------------------------------------------*)
(* Search for strings 'OUTER INPUT' to put in the values of the Input parameters *)
(*-------------------------------------------------------------------*)
Needs["NumericalCalculus`"]

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] ) );
(* see 1606.08659, Baikov, Chetyrkin and Kuhn *)
(* checked, is alright *)

ce2MS=bett2MS/be0;
ce3MS=bett3MS/be0;
ce4MS=bett4MS[nqu]/be0;
betaMS4l[x_]:=- be0 x^2 (1 + ce1 x + ce2MS x^2 + ce3MS x^3);
betaMS5l[x_]:=- be0 x^2 (1 + ce1 x + ce2MS x^2 + ce3MS x^3 + ce4MS x^4);

(* See Eq.(55) of 0007088.ChetyRet.MiniMOM.pdf for be2MMChe[Nff]: *)
(* Their be3 is our 4^3 be2, in MiniMOM which they denote as {tilde MOM} *)
be2MMChe[Nff_]:=(1/4^3)*( ((28965/8)-(3861/8) Zeta[3]) + Nff(-(7715/12)+(175/12) Zeta[3]) + Nff^2 (+(989/54)+(8/9) Zeta[3]));
(* See Eq.(56) of 0007088.ChetyRet.MiniMOM.pdf for be2MMChe[Nff]: *)
(* Their be4 is our 4^4 be3, in MiniMOM which they denote as {tilde MOM} *)
be3MMChe[Nff_]:=(1/4^4)*( ((1380469/8) - (625317/16) Zeta[3] - (772695/32) Zeta[5]) + Nff ( - (970819/24) + (516881/72) Zeta[3] + (1027375/144) Zeta[5]) + Nff^2 (+(736541/324) - (6547/27) Zeta[3] - (9280/27) Zeta[5]) + Nff^3 (-(800/27) + (16/9) Zeta[3]) );

(*  Input c2: *)
(* c2=9.3; *)
 c2=N[be2MMChe[3]/be0];  (* = 9.297; *)

(* Input c3: *)
 c3=N[be3MMChe[3]/be0]; (* 71.4538; *)

 
(* Both c2=be2/be0 and c3=be3/be0 taken from PLB681,336 (Smekal, Maltman, Sternbeck), i.e., the values in MiniMOM scheme for Nf=3 *)
(*-------------------------------------------*)
(* Now the 4-loop MiniMOM beta function, extended to the 4-loop Lambert-MiniMOM scheme: *)
om1=c2/ce1^2; om2=c3/ce1^3; sqom2=+Sqrt[om2];
a0=1+sqom2; a1=om1 +sqom2;

beta[x_] := -be0 x^2  (1 + (a0 ce1) x + (a1 ce1^2) x^2)/(1 + (a0-1)ce1 x + (a1 ce1^2) x^2)/(1 - (a1 ce1^2) x^2);
bt[x_]:=beta[x];
betared[x_]:= x^2  (1 + (a0 ce1) x + (a1 ce1^2) x^2)/(1 + (a0-1)ce1 x + (a1 ce1^2) x^2)/(1 - (a1 ce1^2) x^2);
(*-------------------------------------------*)
c2num=Coefficient[Series[betared[x],{x,0,5}],x^4];
c3num=Coefficient[Series[betared[x],{x,0,5}],x^5];
c4num=Coefficient[Series[betared[x],{x,0,6}],x^6];
c5num=Coefficient[Series[betared[x],{x,0,7}],x^7];
c4=c4num; c5=c5num; 
(*----------------------*)
zu[u_, fi_] := - (1/ce1/E) u^(-be0/ce1) Exp[- I*(be0/ce1) fi];
 (* u=|Q2|/L2 = absQ2/L2, where L2=LambertLambda^2 *)
zup[u_, fi_] := (1/ce1/E) u^(-be0/ce1) Exp[I*(Pi - (be0/ce1) fi)];
zdo[u_, fi_] := (1/ce1/E) u^(-be0/ce1) Exp[I*(-Pi - (be0/ce1) fi)];
(* Note, zup=zdo=zu, for all u and fi, in the evaluation by Mathematica; so I can use zu[u,fi] throughout *)
(*----------------------*)
a2l[u_, fi_] := If[fi >= 0, -(1/ce1) 1/(1 + ProductLog[-1, zu[u, fi]]), -(1/ce1) 1/(1 + ProductLog[1, zu[u, fi]])];

a3l[u_, fi_] := If[fi >= 0, -(1/ce1) 1/(1 -c2/ce1^2 + ProductLog[-1, zu[u, fi]]), -(1/ce1) 1/(1 - c2/ce1^2 + ProductLog[1, zu[u, fi]])];
(* u = |Q^2|/L2 > 0; Q^2=Q2 = L2 Exp[I fi]; Q^2/L^2 = u Exp[I fi] *)
(*-------------------------------------*)
det[u_,fi_]:=If[fi >= 0,(sqom2+1+ProductLog[-1, zup[u, fi]])^2- 4(om1+sqom2),(sqom2+1+ProductLog[1, zdo[u, fi]])^2 - 4(om1+sqom2)];
(* First Sqrt[det[u,fi]] in aLamb[u,fi] when fi >=0: *)
	       sqrtfipos[u_,fi_]:=If[Arg[det[u,fi]]>=0,Sqrt[Abs[det[u,fi]]] Exp[ I (1/2) Arg[det[u,fi]]], Sqrt[Abs[det[u,fi]]] Exp[ I (1/2) (Arg[det[u,fi]]+ 2 Pi)]];
 (* Now, Sqrt[det[u,fi]] in aLamb[u,fi] when fi < 0: *)
	       sqrtfineg[u_,fi_]:=If[Arg[det[u,fi]]<=0,Sqrt[Abs[det[u,fi]]] Exp[ I (1/2) Arg[det[u,fi]]], Sqrt[Abs[det[u,fi]]] Exp[ I (1/2) (Arg[det[u,fi]]- 2 Pi)]];
(* Now 4-loop Lambert-MiniMOM a[u,fi] = a[Q2=u L2 Exp[I fi]]: *) 
 a[u_,fi_]:=If[fi >= 0, (2/ce1) 1/(-(sqom2+1+ProductLog[-1, zup[u, fi]])+sqrtfipos[u,fi]), (2/ce1) 1/(-(sqom2+1+ProductLog[1, zdo[u, fi]])+sqrtfineg[u,fi])];
(* aLamb[u,fi] = aLamb4l(Q^2), where: Q^2 = L2 u Exp[I fi]; |Q^2| = L2 u (in GeV^2) *)

(* a[u,fi] =alpha_s[Q2]/Pi is the solution with beta[x], i.e., an effective "4-loop" solution, in the Lambert-MiniMOM  RSch with a chosen c2 and c3 *)
(*----------------------*) 
(* Beware: ProductLog[1,z] cannot be applied for z real, *)
  (* z has to be in the lower half plane; *)
  (* ProductLog[-1,z] can be applied for z in the upper half AND for z real *)
(*----------------------*)
 (* Let us find uL = QL^2/L^2, i.e. zL: *)
  (* we take into account that c2>0, so for a3l[u,fi] we have no addition Landau pole uL to worry about *)
  
sL3l=ce1^(-ce1/be0); (* sL3l =approx 0.634696; sL such that: tr13l[s]=0 for s<-sL *)
sL=ce1^(-ce1/be0);
(*--------------------*)
tr1[s_]:=If[s > 0, Im[a[s,-Pi]],Im[a[-s,-10.^-8]]];
tr1pl[s_? NumberQ]:=Im[a[s,-Pi]];  (* this is only for s>0 *)

tr13l[s_]:= If[s!=0,If[s > 0, Im[a3l[s, -Pi]], Im[(1/ce1)  1/(1 - c2/ce1^2 + ProductLog[-1, zu[-s, 0]])]],0];
(* tr13l[s]=rho_1(sigma = s L^2,3l) = Im[a3l[s, -Pi] = Im[a3l[Q^2=-s L^2 - I eps]], when s=sigma/L^2 > 0; *)
    (* s = u Exp[- I Pi] if s > 0; *)
    (* tr13l[s] is applicable even when s<0; s= -u Exp[I 0] if s < 0 *)
    (* Note that tr13l[-sL]=0; tr13l[0]=0 *)
    (* tr13l[s]=0 for all s<-sL exactly *)

    zuw[w_] := - (1/ce1/E) Exp[-w be0/ce1] Exp[ I*(be0/ce1) (Pi-10.^-15)];
    (* zuw[w] = zu[Exp[w],-Pi+10.^-15] *)
    tr13lex[w_]:=-(1/ce1) Im[1/(1 - c2/ce1^2 + ProductLog[1, zuw[w]])];
        (* tr13lexp[w] = tr13l[Exp[w]] for s=Exp[w]>0  *)
    (*-----------------------------*)

 (* Now the dispersively calculated a[u,fi] = a[Q^2=(u*L2)Exp[I fi]]  (u>0) *)
    (*-----*)
   (* adis[u,fi] is dispersively calculated a[u,fi]: *)
adis[u_,fi_]:=(1/Pi) (NIntegrate[tr1[s]/(s + u Exp[I fi]),{s,-100.,0}, MaxRecursion -> 24] + NIntegrate[tr1[s]/(s + u Exp[I fi]),{s,0,100.}] + NIntegrate[tr1[Exp[w]]/(1 + u Exp[I fi] Exp[-w]),{w,Log[100.],10.^6}]); 

 (* adis[u,fi]=a[u,fi]; except on the cut, i.e., for u < sL *)
     (*-------------------*)
 (* Until now it works fine *)
 (* One weird thing: tr1[s] is nonzero down to s =appr -100, while tr13l[s] is nonzero only down to s =appr - 0.63; this indicates that we will have: LLamb >> LLamb3l *)

    (*------------------------------------------------*)
    (* Input: kappa: *)
  (* kappa=2 will be now the only possibility here: *)
(* Q_{thr} = kappa mq (q=b,c) in RGErun.m *)
   kappa=2;
    (* kappa=2; *)
 (*------------------------------------------------*)
 (* Input: alphaMSMZ = 0.1184 +-0.0007, (= alpha_s(MZ^2;bMS), i.e., the PDG2010 world average *)
 (* Input: alphaMSMZ = 0.1185 +-0.0006, (= alpha_s(MZ^2;bMS), i.e., the PDG2015 world average *)
  (* Input: alphaMSMZ = 0.1181 +- 0.0011, PDG2016 (Chinese Phys. C40, 100001 (2016)) *)

  (* OUTER INPUT: alphaMSMZ: *)
  (*  alphaMSMZ=0.1181; *)
  (* alphaMSMZ=0.1185;  *)
 (*------------------------------------------------*)
  mb=4.20; mc=1.27; 
(* mb=mbbar(mbbar); mc=mcbar(mcbar); thresholds at kappa*mb and kappa*mc *)

  <<RGErundown5l.m;  (* Note, 5-loop MSbar RGE running *)
  ainMS=aaMSQ;
  (* This is: ainMS= a(MSbar,Q=kappa*mc-0, with nf=3); obtained bt RGE running down with 5-loop MSbar beta function, and 4-loop threshold matching at Qthr=kappa*mq (mq=mb, mc) *)

  (* Note, our preferred input values are: alphaMSMZ=0.1181 (the world average); mb = mb(mb)=4.20; mc=mc(mc)=1.27; kappa=2; -> ainMS=0.0838279; *)
  
   (* Note that now (RGErundown.m) the thresholds (CSK) are implemented with running DOWNWARD, that is why somewhat different results than in the case when thresholds are implemented in the RGE-running upwards *)

  ainMSnf3=ainMS;
  ainMSnf4=aaQ4plus;
  (* ainMSnf4=a(MSbar,Q=kappa*mc+0, with nf=4) *)

(* Now we calculate ain=al(ourRSch;kappa mc - 0)/Pi: *)
  (* Note: we know LMM beta function correctly only to 4-loop (c3); however, c4 is also there as a Pade extrapolation; therefore, we will use for MSbar beta in the RSch-change also 5-loop beta-function: *)   
  ain0=ainMS;
  (* ain0: this is the first (pert.) approx.to a(kappa mc;nqu=3;c2,c3,...) *)
 fr[x_]:=(betaMS5l[x] + be0 x^2 (1 + ce1 x))/(x^2 (1+ce1 x) betaMS5l[x]); 
rhs1=1/ainMS+ ce1 Log[ce1 ainMS/(1 + ce1 ainMS)] + NIntegrate[fr[x],{x,0,ainMS}];
(* fl1[x_]:=approx (c2 + c3 x)/((1+ce1 x) (1 + ce1 x + c2 x^2 + c3 x^3)); *)
fl1[x_]:=(bt[x] + be0 x^2 (1 + ce1 x))/(x^2 (1+ce1 x) bt[x]); 
 eepss=10.^(-7);
aux1[yyy_]:=1/yyy + ce1 Log[ce1 yyy/(1 + ce1 yyy)];
(* aux1[ain] + NIntegrate[fl1[xx], {xx, eepss, ain}] - rhs1 = 0 has to be *)
 LHS[yyy_? NumberQ]:=aux1[yyy] + NIntegrate[fl1[xxx], {xxx, eepss, yyy}];
sol1= FindRoot[LHS[yyy]-rhs1,{yyy,0.089}];
ain=yyy/.sol1;
   (* ain: This is ain=a(Q^2=(kappa mc)^2-0; nqu=3; c2,c3,c4,...) pQCD, in RSch LMM of our bt[a]; we will use the same LLamb=Lambda(nqu=3) that is obtained from this a(Q^2;c2,...) also in our A1(Q^2;c2,...) *)
  (*-------------------------------------------*)
 (* Input: mu2in = mu_{in}^2 *)
mu2in=(kappa*mc)^2; 
  (* mu2in=(kappa*mc)^2;  *)
  (*-------------------------------------------*)
(* Now calculating ub=(mu2in/L2) value such that a[ub,0]=ain: *)
  ub0=200.;
sol0=FindRoot[a[uub,0]==ain,{uub,ub0}];
ub=uub/.sol0; 
    (* ub is ub=(mu2in/L2) exact, because A1 uses the same LLambda (=LambertLamb=Sqrt[L2]) as the underlying pQCD coupling a[u,0] = a[Q2/L2,0] *)
   L2=mu2in/ub; LLamb=Sqrt[L2]; (* i.e., LambertLambda(nqu=3) for a{Q^2;c2,...) *)
    (*--------------------*)
  (* Now solving equations for A1 - a ~ (L2/Q^2)^5, considering s0 as input: *)
     (* We use notation: int_{-sL}^s0  ds s^j tr1[s] = Pj (j=0,1,2,3) *)

     (* We use notation: int_{-sL}^s0  ds s^j tr1[s] = Pj (j=0,1,2,3) *)


  (* OUTER INPUT: s0: *)
  (* s0=620.0; *)

  (* s0=3.29; *)   (* s0 has to be such that 1s0>0 and s2s0>0 *)
  (* s0 will be adjusted so that rtau=0.203 is obtained *)
    lns0=Log[s0];
  (*----------------------------------*)
   (* Note: correction in the formulas below: sL^k -> (-sL)^k; for k=0,1,2,3; Correction made on March 16, 2012 *)

P0s0=(1/Pi) (NIntegrate[tr1[s], {s, -101., -sL - 0.001}, MaxRecursion -> 24] + NIntegrate[tr1[s], {s, -sL - 0.001, -sL + 0.001}, MaxRecursion -> 24] + NIntegrate[tr1[s], {s, -sL + 0.001, 0}] + NIntegrate[tr1[s], {s, 0, s0}]);
 (* For faster (and slightly less precise) evaluations: P0s0fast=(1/Pi) (NIntegrate[tr1[s], {s, -100., s0}]); etc. *)
P1s0=(1/Pi) (NIntegrate[s tr1[s], {s, -101., -sL - 0.001}, MaxRecursion -> 24] + NIntegrate[s tr1[s], {s, -sL - 0.001, -sL + 0.001}, MaxRecursion -> 24] + NIntegrate[s tr1[s], {s, -sL + 0.001, 0}] + NIntegrate[s tr1[s], {s, 0, s0}]);
P2s0=(1/Pi) (NIntegrate[s^2 tr1[s], {s, -101., -sL - 0.001}, MaxRecursion -> 24] + NIntegrate[s^2 tr1[s], {s, -sL - 0.001, -sL + 0.001}, MaxRecursion -> 24] + NIntegrate[s^2 tr1[s], {s, -sL + 0.001, 0}] + NIntegrate[s^2 tr1[s], {s, 0, s0}]);
P3s0=(1/Pi) (NIntegrate[s^3 tr1[s], {s, -101., -sL - 0.001}, MaxRecursion -> 24] + NIntegrate[s^3 tr1[s], {s, -sL - 0.001, -sL + 0.001}, MaxRecursion -> 24] + NIntegrate[s^3 tr1[s], {s, -sL + 0.001, 0}] + NIntegrate[s^3 tr1[s], {s, 0, s0}]);
  (* P4s0=(1/Pi) (NIntegrate[s^4 tr1[s], {s, -101., -sL - 0.001}, MaxRecursion -> 12] + NIntegrate[s^4 tr1[s], {s, -sL - 0.001, -sL + 0.001}, MaxRecursion -> 12] + NIntegrate[s^4 tr1[s], {s, -sL + 0.001, 0}] + NIntegrate[s^4 tr1[s], {s, 0, s0}]); *)

    sm=10.^40; lnsm=Log[sm];
  (* J0=(1/Pi) NIntegrate[tr1pl[Exp[w]],{w,lnsm,Infinity}]; *)
  J0=(1/Pi) NIntegrate[tr1pl[Exp[w]],{w,lnsm,10.^10},MaxRecursion->20, WorkingPrecision -> 16];
    Pm1s0= -(1/Pi) NIntegrate[tr1pl[Exp[w]],{w,lns0,lnsm}, MaxRecursion -> 24] - J0;

  (* The relative numerical precision of calculation of Pjs0 is (del Pjs0)/Pjs0 ~10^-10; *)
(* This is IMPORTANT to keep in mind *)

   P0s0=Re[P0s0]; P1s0=Re[P1s0]; P2s0=Re[P2s0]; P3s0=Re[P3s0];
  (* P4s0=Re[P4s0]; *)
  Pm1s0=Re[Pm1s0];

  (*--------------------------------------------------------------*)
  (*--------------------------------------------------------------*)
  (* Now demanding A[u]-a[u] ~ 1/u^5, i.e., s1 fixed (varied) by hand: *) 


  (* OUTER INPUT: s1in(=s1s0=s_1)*)
  (* s1in=5.0; *)

  (* NOte that, in the first step, we should adjust also f1sin, f2sin, f3sin, s2in, s3in --> so that aaa procedure (below) would give us a physically acceptable solutions *)

  (* OUTER INPUT: f1sin; f2sin; f3sin; s2in; s3in: *)  
  (* f1sin=-2.0; f2sin=10.5; f3sin=5.57; s2in=22.2; s3in =510.9; *)
   (* f1sin=-0.5; f2sin=10.5; f3sin=6.67; s2in=71.2; s3in =740.9; *)

s1=s1in;
 
  aaa= FindRoot[{f1s s1^-1 + f2s s2^-1 + f3s s3^-1  == Pm1s0 && f1s + f2s + f3s == P0s0 && f1s s1 + f2s s2 + f3s s3  == P1s0 && f1s s1^2 + f2s s2^2 + f3s s3^2  == P2s0 && f1s s1^3 + f2s s2^3 + f3s s3^3 == P3s0}, {{f1s,f1sin},{f2s,f2sin},{f3s,f3sin},{s2,s2in},{s3,s3in}}];

s1s0=s1in;

f1ss0=f1s/.aaa;
f2ss0=f2s/.aaa; 
f3ss0=f3s/.aaa; 
s2s0=s2/.aaa; 
s3s0=s3/.aaa; 

  (* Checking: f1ss0 s1s0^-1 + f2ss0 s2s0^-1  + f3ss0 s3s0^-1 - Pm1s0 (=0); f1ss0 + f2ss0 + f3ss0 - P0s0 (=0); f1ss0 s1s0 + f2ss0 s2s0 + f3ss0 s3s0 - P1s0 (=0); f1ss0 s1s0^2 + f2ss0 s2s0^2 + f3ss0 s3s0^2 - P2s0 (=0); f1ss0 s1s0^3 + f2ss0 s2s0^3  + f3ss0 s3s0^3 - P3s0 (=0); *)
  (* Checking: *)
  checkm1=f1ss0 s1s0^-1 + f2ss0 s2s0^-1 +  + f3ss0 s3s0^-1  - Pm1s0; 
  check0=f1ss0 + f2ss0 + f3ss0 - P0s0;
  check1=f1ss0 s1s0 + f2ss0 s2s0 + f3ss0 s3s0  - P1s0;
  check2=f1ss0 s1s0^2 + f2ss0 s2s0^2  + f3ss0 s3s0^2  - P2s0;
  check3=f1ss0 s1s0^3 + f2ss0 s2s0^3 + f3ss0 s3s0^3  - P3s0;
  check=Sqrt[checkm1^2/Pm1s0^2+check0^2/P0s0^2+check1^2/P1s0^2+check2^2/P2s0^2+check3^2/P3s0^2];
  (*--------------------------------------------------------------*)

   
(* u=Q2/L2,  Q2 can be complex, i.e., u can be complex *)
  (* LogLinearPlot[Aprel[Q2/L2], {Q2, 10.^-3, 10.^2}] *)
  (* The is maximum is at Q2=appr 0.03 GeV^2; we will need max at 0.09 GeV^2 *)

  sm=10.^40; lnsm=Log[sm];
  (*    J0=(1/Pi) NIntegrate[tr1pl[Exp[w]],{w,lnsm,Infinity}]; *)
    J1=(1/Pi) NIntegrate[tr1pl[Exp[w]] Exp[-w],{w,lnsm,Infinity}];
    J2=(1/Pi) NIntegrate[tr1pl[Exp[w]] Exp[-2 w],{w,lnsm,Infinity}];
    J3=(1/Pi) NIntegrate[tr1pl[Exp[w]] Exp[-3 w],{w,lnsm,Infinity}];
  
  A1an[u_? NumberQ,fi_? NumberQ]:= ( (1/Pi) NIntegrate[tr1pl[Exp[w]]/(1 + u Exp[I fi] Exp[-w]),{w,lns0,lnsm},MaxRecursion->30, WorkingPrecision -> 16] + f1ss0/(u Exp[I fi] +s1s0) + f2ss0/(u Exp[I fi] +s2s0) + f3ss0/(u Exp[I fi] +s3s0) + J0 - u Exp[I fi] J1 + u^2 Exp[I 2 fi] J2 );
  (* s= Exp[w]; tr1pl[Exp[w]]=tr1ex[w]: *)

  (* Calculation of our {tilde A2}^{analytic}[u,fi], in our model, *)
    (* where u=|Q^2|/LLamb^2 and fi=Arg[Q^2]: *)
tA2an[u_? NumberQ,fi_? NumberQ]:=(1/be0/Pi) (u Exp[I fi]) ( NIntegrate[ Exp[-w] tr1pl[Exp[w]]/(1 + u Exp[I fi] Exp[-w])^2,{w,lns0,lnsm},MaxRecursion->30, WorkingPrecision -> 16] + (f1ss0*Pi)/(u Exp[I fi]+s1s0)^2 + (f2ss0*Pi)/(u Exp[I fi]+s2s0)^2+ (f3ss0*Pi)/(u Exp[I fi]+s3s0)^2  + J1 - 2  u Exp[I fi] J2 + 3 u^2 Exp[I 2 fi] J3 );

   (* Calculation of our {tilde A3}^{analytic}[u,fi], in our model, *)
    (* where u=|Q^2|/LLamb^2 and fi=Arg[Q^2]: *)
  tA3an[u_? NumberQ,fi_? NumberQ]:=(1/2/be0^2/Pi)(u Exp[I fi]) ( NIntegrate[ Exp[-w] tr1pl[Exp[w]]/(1 + u Exp[I fi] Exp[-w])^2,{w,lns0,lnsm},MaxRecursion->30, WorkingPrecision -> 16] -2  NIntegrate[ Exp[- w] tr1pl[Exp[w]]/(1 + u Exp[I fi] Exp[-w])^3,{w,lns0,lnsm},MaxRecursion->30, WorkingPrecision -> 16] + (f1ss0*Pi) (-s1s0 + u Exp[I fi])/(s1s0 + u Exp[I fi])^3  + (f2ss0*Pi) (-s2s0 + u Exp[I fi])/(s2s0 + u Exp[I fi])^3  + (f3ss0*Pi) (-s3s0 + u Exp[I fi])/(s3s0 + u Exp[I fi])^3  );

   (* Calculation of our {tilde A4}^{analytic}[u,fi], in our model, *)
    (* where u=|Q^2|/LLamb^2 and fi=Arg[Q^2]: *)
  tA4an[u_? NumberQ,fi_? NumberQ]:=(1/6/be0^3/Pi)(u Exp[I fi]) ( NIntegrate[ Exp[-w] (1 - 4 Exp[-w] u Exp[I fi] + Exp[-2 w] u^2 Exp[2 I fi]) tr1pl[Exp[w]]/(1 + u Exp[I fi] Exp[-w])^4,{w,lns0,lnsm},MaxRecursion->30, WorkingPrecision -> 16] + (f1ss0*Pi) (s1s0^2 - 4 s1s0 u Exp[I fi] + u^2 Exp[2 I fi] )/(s1s0 + u Exp[I fi])^4  + (f2ss0*Pi) (s2s0^2 - 4 s2s0 u Exp[I fi] + u^2 Exp[2 I fi] )/(s2s0 + u Exp[I fi])^4  + (f3ss0*Pi) (s3s0^2 - 4 s3s0 u Exp[I fi] + u^2 Exp[2 I fi] )/(s3s0 + u Exp[I fi])^4  );
   (* tAjan[u,fi] = {tilde A}_j(Q^2) with Q^2 = L^2 (u Exp[I fi])^2 *) 
 

  (* Plot[{Re[tA2an[(mtau2/L2), fi]], Im[tA2an[(mtau2/L2), fi]]}, {fi, 0, Pi - eps}] *)
  (*----------------------------------------------------_*)
M0=Sqrt[s0 L2];

(* Aour=LogLinearPlot[A1an[Q2/L2, 0], {Q2, 10.^-3, 10.^1}] *)

  fm=FindMaximum[{A1an[Q2/L2,0],Q2>0}, {Q2, 0.1}];
  (*  fm=FindMaximum[A1an[Q2/L2,0], {Q2, 0.14}, AccuracyGoal -> 18, PrecisionGoal -> 20, WorkingPrecision -> 40]; *)
		Q2max=Q2/.fm[[2]];		
(* alphaLatt=LogLinearPlot[Pi A1an[Q2latt/1.82^2/L2, 0], {Q2latt, 10.^-3, 10.^2},PlotRange->{{10.^-3,10.^2},{0.,1.5}},Frame->True] *)

(*----------------------------------------------------_*)
(*----------------------------------------------------_*)
(*-----------------------------------------------------*)
		 (* Now the "exact" calculation of rtau(D=0) with 4A Renormalon approach: *)
  A1[Q2_? NumberQ,fi_? NumberQ]:=A1an[Q2/L2,fi];
 tA2[Q2_? NumberQ,fi_? NumberQ]:=tA2an[Q2/L2,fi];
 tA3[Q2_? NumberQ,fi_? NumberQ]:=tA3an[Q2/L2,fi];
 tA4[Q2_? NumberQ,fi_? NumberQ]:=tA4an[Q2/L2,fi];
		 
 <<Adl4l3dRenMod4Ash.m;
  (* write: rtau  --> this gives rtau(D=0) in this approach *)		 

  rtauRenM=rtau;


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

apert[Q2_? NumberQ,fi_? NumberQ]:=a[Q2/L2,fi];
  
 diff[Q2_,fi_]:=A1[Q2,fi]-apert[Q2,fi];

  rdiff[Q2_,fi_]:=1 - apert[Q2,fi]/A1[Q2,fi];
  (* |rdiff| cannot be smaller than ~`0^-9, because of the limited relative numerical precision of the Pjs0 qunatities *)

  coeff[Q2_,fi_]:= (Q2 Exp[I fi])^5 diff[Q2,fi]; (* coeff should be roughly constant at Q2>>L2 *)
  
   (* Apparently works for |Q2/L2| < 10.^38 or so; and any fi *)

 (* LogLinearPlot[A1[Q2, 0], {Q2, 10.^-3, 10.^1}] *)
  
  (* Plot[{Re[tA2[mtau2, fi]], Im[tA2[mtau2, fi]]}, {fi, 0, Pi - eps}] *)
  (*----------------------------------------------------_*)
M0=Sqrt[s0 L2];

(* Aour=LogLinearPlot[A1[Q2, 0], {Q2, 10.^-3, 10.^1}] *)
  (*-----------------------------------------------------*)
(*-----------------------------------------------------*)
(* Now calculating rtau using [2/2] diagonal Baker-Gammel (dBG) for Adler: *)

  (* First the diagonal Baker-Gammel (dBG) of DAdlan in 4danQCD: *)
   dvmptser[xx_]:=xx + td1 xx^2 + td2 xx^3 + td3 xx^4;
		(* Adler is D[Q2]= A1[Q2]+ td1 tA2[Q2] + td2 tA3[Q2] + td3 tA4[Q2]; note, here and from now on we assume that eCC=1, i.e., mu^2=Q^2 *)

   (* This dvmptser[xx] is the AUXILIARY power series -> leads via dPA to dBGA: *)
 Needs["NumericalCalculus`"];

dPAdvmpt[nn_,xx_]:=PadeApproximant[SetPrecision[dvmptser[xxp],30],{xxp,0,nnp}]/.{nnp->nn,xxp->xx};
(*-------------*)
   (* Now calculation of alt[j] and ut[j] (j=1,2), where: *)
   (* dPAdvmpt[n,x] = Sum_{j=1}^n alt[n,j] x/(1 + x ut[n,j] be0) = Sum_{j=1}^n alt[n,j] a(1-l.)(Q[n,j]^2), where Q[n,j]^2 = Q^2 Exp[ut[n,j]]; note, we have mu^2=Q^2 in our dv[Q^2]=dvtmpt[Q^2] *)

  (* We use: Sum_{j=1}^2 alt[n,j]/(1 + x ut[n,j] be0)= Apart[dPAdvmpt[n,x]/x] in order to identify alt[n,j] and ut[n,j]'s;  for n=1; and n=2; *)
(* ' *)

(* Apart[(A + B x)/((x - k1) (x - k2))]; Note: k1 and k2 are complex (conjugate) not real; use the formula obtained from: Apart[(A + B x)/((x - k1) (x - k2))] *)

AA=Numerator[Simplify[dPAdvmpt[2, x]/x]]/.{x->0};
BB=Coefficient[Numerator[Simplify[dPAdvmpt[2, x]/x]],x];
deno=NSolve[Denominator[Simplify[dPAdvmpt[2, x]/x]] == 0, x];		
kk1=x/.deno[[1]];
kk2=x/.deno[[2]];
alt[2,1]= (AA + BB kk1)/(kk1-kk2)/(-kk1);
alt[2,2]= (-AA - BB kk2)/(kk1-kk2)/(-kk2);
ut[2,1]=-1/(kk1)/be0;
ut[2,2]=-1/(kk2)/be0;

 (* this is for [2/2] dBG of Adler function *)

(* Cross-checking: *)
Appt[nn_,xx_] := Sum[alt[nn,jj]/(1 + xx ut[nn,jj] be0), {jj, 1, nn}];
(* Then check: Together[Simplify[N[Appt[n,x] - Apart[dPAdvmpt[n,x]/x]]]] is practically zero, for n=1,2  *)
 (* Yes, I checked, high precision  ~ 10^-16 *)

  Do[Do[kapps[nn,jj]=Exp[ut[nn,jj]],{jj,1,nn}],{nn,1,2}];

(* D[Q2]=A1(Q2) + td1 {tilde A}_2(Q2) + td2 {tilde A}_3(Q2) + td3 {tilde A}_4(Q2) + O(A_5) = alt[2,1] A1(kapps[2,1]Q2) + alt[2,2] A1([kapps[2,2]Q2) + O(A5), where the last expression is [2/2] dBG *) 
 (*---------*)
 mtau2= 1.77686^2;
 a4dk11complex[ffi_]:=A1[kapps[1,1] mtau2,ffi];
(* a4dk11complex[fi] = A1[Q^2=kapps[1,1] mtau2 Exp[i fi]]; -Pi+eeps < fi < Pi - eeps *)
 (*---------*)
   a4dk22=A1[kapps[2,2] mtau2,0];
 a4dk22complex[ffi_]:=A1[kapps[2,2] mtau2,ffi];
 (* aEEk22complex[fi] = aEE[Q^2=kapps[2,2] mtau2 Exp[i fi]]; -Pi+eeps < fi < Pi - eeps *)

  a4dk21=A1[kapps[2,1] mtau2,0];
a4dk21complex[ffi_]:=A1[kapps[2,1] mtau2,ffi];
  (* a4dk21complex[fi] = A1[Q^2=kapps[2,1] mtau2 Exp[i fi]]; -Pi+eeps < fi < Pi - eeps *)
  (*----------------------------*)
	(* Beware: only correct if eCC=1: *)
	       mxpts=1000;
     (* theta_max=Pi-eeps taken *)
		 eeps=10.^-4;

tauk22= (1/2/Pi) NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) a4dk22complex[fi],{fi,-Pi+eeps,+Pi-eeps},MaxPoints->mxpts];
tauk21= (1/2/Pi) NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) a4dk21complex[fi],{fi,-Pi+eeps,+Pi-eeps},MaxPoints->mxpts];

rtaudBG22=alt[2,1] tauk21 +  alt[2,2] tauk22;
rtaudBG22

(* rtaudBG22 = rtau(D=0) with dBG[2/2] used for Adler D(Q^2=mtau^2 Exp[I fi]) *)
 (*----------------------------*)
(*-----------------------------------------------------*)
(* Now an alternative determination of the rtau: Adler function is calculated as (m)TPS in tAn[Q2]'s: ' *)

	       mxpts=1000;
     (* theta_max=Pi-eeps taken *)
		 eeps=10.^-4;
 CC=0;
 eCC=Exp[CC];

		 
tau1=2 (1/2/Pi) Re[ NIntegrate[(1 + Exp[I fi])^3 (1 - Exp[I fi]) (A1an[eCC (mtau2/L2),fi]),{fi,0.,+Pi-eeps},MaxPoints->mxpts]];
tau2v3 = 2 (1/2/Pi) NIntegrate[Re[(1 + Exp[I fi])^3 (1 - Exp[I fi]) (tA2an[eCC (mtau2/L2),fi])],{fi,0.,+Pi-eeps},MaxPoints->mxpts];
tau3v3 = 2 (1/2/Pi) NIntegrate[Re[(1 + Exp[I fi])^3 (1 - Exp[I fi]) (tA3an[eCC (mtau2/L2),fi])],{fi,0.,+Pi-eeps},MaxPoints->mxpts];
tau4v3 = 2 (1/2/Pi) NIntegrate[Re[(1 + Exp[I fi])^3 (1 - Exp[I fi]) (tA4an[eCC (mtau2/L2),fi])],{fi,0.,+Pi-eeps},MaxPoints->mxpts];
  
rtau0=tau1;
rtau1=tau2v3 td1;
rtau2=tau3v3 td2;
rtau3=tau4v3 td3;

rtauTPS=rtau0+rtau1+rtau2+rtau3;
rtau0
rtau1
rtau2
rtau3
rtauTPS
(*-------------------------------------------------------*)
(*-------------------------------------------------------*)
(*-------------------------------------------------------*)
(*-------------------------------------------------------*)
(*---------------------------------------------------------*)
Clear[Q];
		 
(* shortcut from here, call: << 4l3danQCDcouplrt0XXXal0YYYYN.m; *)
(* and then continue below *)

 ff[x_]:=0.8;
 
AourPi = LogLinearPlot[Pi A1[Q2, 0], {Q2, 10^-3, 10^1}, 
  PlotRange -> {{10^(-3), 10^1}, {0, 1.2}}, 
  AxesOrigin -> {10^(-3), 0.}, Frame -> True, 
  FrameLabel -> {Superscript[Q, 2][Superscript[GeV, 2]], \[Pi] A == 
     Subscript[\[Alpha], s]}, PlotStyle -> Black, 
  Epilog -> {Inset[
     Style[Superscript[Subscript[\[Alpha], s], Style[OverBar[MS], 6]][
        Subscript[M, Z]] == 0.1181, 10], Scaled[{0.81, 0.93}]], 
    Inset[Style[Subsuperscript[r, \[Tau], "(D=0)"] == 0.201, 10], 
     Scaled[{0.84, 0.82}]]}];
  AourPi4l3drt0XXXal0YYYYN=AourPi;
(* Save["AourPi4l3drt0XXXal0YYYYN.save",AourPi4l3drt0XXXal0YYYYN]; *)
 
AourPib=LogLinearPlot[Pi A1[Q2, 0], {Q2, 10^-3, 10^1}, 
 PlotRange -> {{10^(-3), 10^1}, {0, 1.2}}, 
 AxesOrigin -> {10^(-3), 0.}, Frame -> True, 
 FrameLabel -> {Superscript[Q, 2][Superscript[GeV, 2]], \[Pi] A == 
    Subscript[\[Alpha], s]}, PlotStyle -> Black, 
 Epilog -> {Inset[Style["(b)", 14], Scaled[{0.9, 0.75}]], 
   Inset[Style[
     Superscript[Subscript[\[Alpha], s], Style[OverBar[MS], 6]][
       Subscript[M, Z]] == 0.1181 \[PlusMinus] 0.004, 9], 
			      Scaled[{0.74, 0.93}]]}];
  AourPib4l3drt0XXXal0YYYYN=AourPib;
(* Save["AourPib4l3drt0XXXal0YYYYN.save",AourPib4l3drt0XXXal0YYYYN]; *)

 << FigAlphaLatt.save;

 FigComb=Show[AourPi, FigAlphaLattvsQ2];
(* Note: FigAlphaLattvsQ is lattice data, generated in prg. IlgenfritzData.m *)
  (* Export["FigComb4l3drt0XXXal0YYYYN.pdf",FigComb] *)
If[alphaMSMZ==0.1181,aMSnf3plot[Q2_? NumberQ]:=If[Q2>0.405,aMSnf3[Q2],10.]];
(* aMSnf3plot[Q2] = a(Q2) in nf=3 MSbar scheme, for Q2>0 *)
 
If[alphaMSMZ==0.1181,aplot[Q2_? NumberQ]:=If[Q2>1.273,apert[Q2,0],10.]];
(* aplot=aLambplot, for a(Q2) in Lambert MiniMOM scheme, with nf=3, and Q2>0 *)

Clear[Q];
AaaMSPi = LogLinearPlot[{A1[Q2, 0], aplot[Q2], aMSnf3plot[Q2]}, {Q2, 10^-3, 10^2}, 
  PlotRange -> {{10^(-3), 10^2}, {0, 1.2}}, 
  AxesOrigin -> {10^(-3), 0.}, Frame -> True, 
  FrameLabel -> {Superscript[Q, 2][Superscript[GeV, 2]]},PlotStyle->{Black,Dashed,Dotted},Epilog -> {Inset[Style[A, 12], Scaled[{0.38, 0.6}]], Inset[Style[a, 12], Scaled[{0.68, 0.55}]],Inset[Style[OverBar[a], 12], Scaled[{0.55, 0.9}]],Inset[Style[
  Superscript[Subscript[\[Alpha], s], Style[OverBar[MS], 6]][
    Subscript[M, Z]] == alphaMSMZ, 9], Scaled[{0.22, 0.93}]], Inset[
 Style[Subscript[Subscript[r, \[Tau]]["D=0"], "th"] == 
   Round[Re[rtaudBG22]*10^4]/10.^4, 9], Scaled[{0.24, 0.80}]]}];

  (* Export["AaaMSPi4lrt0XXXal0YYYYN.pdf",AaaMSPi] *)

AaaMSPiint = 
    LogLinearPlot[{A1[Q2,0], A3o[2,Q2,0], A1[Q2,0]^2, 
   aplot[Q2], aMSnf3plot[Q2]}, {Q2, 10^-3, 10^2}, 
  PlotRange -> {{10^(-3), 10^2}, {-0.1, 0.4}}, 
  AxesOrigin -> {10^(-3), 0.}, Frame -> True, 
  FrameLabel -> {Superscript[Q, 2][Superscript[GeV, 2]]}, 
   PlotStyle -> {Black, {Black,Dashed}, DotDashed, Dashed, Dotted}, 
  Epilog -> {Inset[Style[A, 12], Scaled[{0.37, 0.7}]], 
    Inset[Style[Subscript[A, 2], 12], Scaled[{0.40, 0.28}]], 
    Inset[Style[Superscript[A, 2], 12], Scaled[{0.35, 0.38}]], 
    Inset[Style[a, 12], Scaled[{0.65, 0.75}]], 
    Inset[Style[OverBar[a], 12], Scaled[{0.54, 0.9}]], 
    Inset[Style[
      Superscript[Subscript[\[Alpha], s], Style[OverBar[MS], 6]][
        Subscript[M, Z]] == alphaMSMZ, 9], Scaled[{0.18, 0.93}]], Inset[Style[
    Subscript[Subsuperscript[r, \[Tau], "(D=0)"], "th"] == 0.201, 9], 
      Scaled[{0.17, 0.80}]]}];
 

  (* Export["AaaMSPi4lrt0201al01181Nint.pdf",AaaMSPiint]; *)
(*---------------------------------------------------------------------------_*)
(*---------------------------------------------------------------------------_*)
(*---------------------------------------------------------------------------_*)
(*---------------------------------------------------------------------------_*)
(*---------------------------------------------------------------------------_*)
(*---------------------------------------------------------------------------_*)
(*---------------------------------------------------------------------------_*)
(*---------------------------------------------------------------------------_*)
(*---------------------------------------------------------------------------_*)

  (* Now calculation of the log-deriv analogs and power analogs, for noninteger indices in general (see, e.g., prg. 2danQCDfor Mathematica_ext.m as an example): *)

   (* A1[Q2,fi] = A(Q^2 = Q2 Exp[I fi]);  Q2=|Q^2|, in GeV^2; fi=Arg(Q^2), -Pi =< fi < Pi *)
  (* apert[Q2,fi] = a(Q^2=Q2 Exp[I fi]), i.e., the (underlying) perturbative a=alpha_s/Pi in the Lambert-MiniMOM renormalization scheme *)

  (* s = sigma/L2; u = |Q^2|/L2; Arg[Q^2]=fi  *)
  tr1MA1l[s_]:=(Pi/be0)(1/(Log[s]^2+Pi^2));

  (* In  tAMA1l[nu,u,fi] we take into account that: *)
  (* Limit[1/Log[u]^(nu + 1) - PolyLog[-nu, 1/u]/Gamma[nu + 1], u -> 1] = -Zeta[-nu]/Gamma[nu+1], which holds for all real nu, it seems *)

Delt[nu_]:=10.^-2 + HeavisideTheta[nu-3] (nu-3) 2 10.^-2; 

 Limt1l[u_,nu_,fi_]:=-(Zeta[1 - nu]/Gamma[nu]) + ((u Exp[I fi]-1) Zeta[-nu])/Gamma[nu] - ((u Exp[I fi]-1)^2 (Zeta[-1 - nu] + Zeta[-nu]))/(2 Gamma[nu]);
  tAMA1l[nu_,u_,fi_]:=If[Abs[u Exp[I fi]-1]>Delt[nu],(1/be0^nu)(1/(Log[u] + I fi)^nu - PolyLog[-nu+1,(1/u) Exp[-I fi]]/Gamma[nu]),(1/be0^nu) Limt1l[u,nu,fi]];
(* tAMA1l[nu,u,fi]={tilde A}_{nu}^{MA,1-l.}(Q^2= u L2 Exp[I fi]) *)

 (* NOte, at u =approx +1, we have singularity of the type infty-infty, so around there we replace tAMA1l[nu,u,fi] with Limit[tAMA1l[nu,u,fi], u Exp[I fi] ->1] *)

  (* tA[nu,Q2,fi] = {tilde A}_{nu}(Q^2=Q2 Exp[I f]), with Q2=|Q^2|; Q^2 = Q2 Exp[I fi] (u=Q2/L2) *)
 (* ssig Exp[-I fi]/u = sigma/Q^2 *)

 (* The above integration cannot be applied for nupp=1 (indeterminate expressions) --> there we have to choose, e.g., nupp=1.000000001 *)

  (* Below: alternative tA (tAalt), where nupp can be any integer without danger -- using PolyLog[nu,x] instead of GenLia[-nu,-x]; (however, tAalt is may not good for making derivatives of tAalt[nu] with respect to nu - check this) *)
  (* It seems to be OK, we will declare now tAal as: tA: *)

f1pss0=0; f1ppss0=0;
wm[nupp_]:=If[nupp>-0.45,10.^6,10.^8];
  tA[nupp_? NumberQ,Q2_? NumberQ,fi_? NumberQ]:= tAMA1l[nupp,Q2/L2,fi] + ((-1)/(Pi*(be0^(nupp-1))*Gamma[nupp]))*((Pi f1ss0/s1s0) PolyLog[-nupp+1,-s1s0 Exp[-I fi] L2/Q2] + (Pi f2ss0/s2s0) PolyLog[-nupp+1,-s2s0 Exp[-I fi] L2/Q2] + (Pi f3ss0/s3s0) PolyLog[-nupp+1,-s3s0 Exp[-I fi] L2/Q2] - NIntegrate[(tr1MA1l[ssig]/ssig) PolyLog[-nupp+1,-ssig Exp[-I fi] L2/Q2],{ssig,0,s0},MaxRecursion->30,WorkingPrecision->16] + NIntegrate[((tr1pl[ssig]-tr1MA1l[ssig])/ssig) PolyLog[-nupp+1,-ssig Exp[-I fi] L2/Q2],{ssig,s0,10 s0},MaxRecursion->30,WorkingPrecision->16] + NIntegrate[(tr1pl[Exp[w]]-tr1MA1l[Exp[w]]) PolyLog[-nupp+1,-Exp[w] Exp[-I fi] L2/Q2],{w,Log[10 s0],wm[nupp]}, MaxRecursion -> 30,WorkingPrecision->16]);
  (* tA[nu,Q2,fi]= {tilde A}_{nu}(Q^2=Q2 Exp[I fi]); Q2=|Q^2| *)
  (* Valid for all -1 < nupp; even for integer nupp *)
  (* Perhaps 10^6 is too high value of smax in the integration? *)
  
  (* I have to check, this still does not work well; Aug.16, 2016 *)
(* Now as the above tA, but with the upper integration limit (wm) adjustable: *)
  (* usually wm >= 10.^5 *)

  (* For example: tA[4,Q2,fi] = tA4[Q2,fi] (coincide to high precision *)
  (* For example: tA[3,Q2,fi] = tA3[Q2,fi] (coincide to high precision *)
  (* For example: tA[2,Q2,fi] = tA2[Q2,fi] (coincide to high precision *)
  (* For example: tA[1,Q2,fi] = A1[Q2,fi] (coincide to high precision *)

   (* However, for Q^2 near the negative semiaxis: *)
   (* However, if |Pi-fi|< 10^-2 -> tA[4,Q2,fi] become unreliable (but tA4[Q2,fi] remains reliable *)
    (* If |Pi-fi| < 0.4 10^-2 ->  tA[3,Q2,fi] become unreliable (but tA3[Q2,fi] remains reliable *)
  (* If |Pi-fi| < 0.3 10^-2 ->  tA[2,Q2,fi] become unreliable (but tA2[Q2,fi] remains reliable *)
  (* If |Pi-fi| < 0.1 10^-2 ->  tA[1,Q2,fi] become unreliable (but A1[Q2,fi] remains reliable *)

   (* tA[0,Q2,fi]=1, as it should be *)
   
   <<Li__nu.m;
 wm[nupp_]:=If[nupp>-0.45,10.^6,10.^8];
 tAalt[nupp_,Q2_,fi_]:= tAMA1l[nupp,Q2/L2,fi] + ((-1)/(Pi*(be0^(nupp-1))*Gamma[nupp]))*((Pi f1ss0/s1s0) polylog[-nupp+1,-s1s0 Exp[-I fi] L2/Q2] + (Pi f2ss0/s2s0) polylog[-nupp+1,-s2s0 Exp[-I fi] L2/Q2]+ (Pi f3ss0/s3s0) polylog[-nupp+1,-s3s0 Exp[-I fi] L2/Q2]- (Pi f1pss0/s1s0^2) (polylog[-nupp,-s1s0 Exp[-I fi] L2 /Q2] - polylog[-nupp+1,-s1s0 Exp[-I fi] L2 /Q2]) + (Pi f1ppss0/s1s0^3) (polylog[-nupp-1,-s1s0 Exp[-I fi] L2 /Q2] - 3 polylog[-nupp,-s1s0 Exp[-I fi] L2 /Q2] + 2 polylog[-nupp+1,-s1s0 Exp[-I fi] L2 /Q2]) - NIntegrate[(tr1MA1l[ssig]/ssig) polylog[-nupp+1,-ssig Exp[-I fi] L2/Q2],{ssig,0,s0}] + NIntegrate[((tr1pl[ssig]-tr1MA1l[ssig])/ssig) polylog[-nupp+1,-ssig Exp[-I fi] L2/Q2],{ssig,s0,10 s0}] + NIntegrate[(tr1pl[Exp[w]]-tr1MA1l[Exp[w]]) polylog[-nupp+1,-Exp[w] Exp[-I fi] L2/Q2],{w,Log[10 s0],wm[nupp]}, MaxRecursion -> 20]); 

   (* This gives the same result as tA[nupp,Q2,fi] *)

   (*--------------------------------------------------------------------*)
(* Now calculation of tkNf[s] = {tilde k}_N(s): *)

	       (* For this we first need: the coefficients k1[m],k2[m],k3[m],k4[m] which allow us to relate {tilde A_m} with (A_1)^m, for m not necessarily integer: *)	       
(* Note: m=1,2,3,...; n can be real noninteger: n=s: *)
(* ES[m,n]=S_m(n) = Sum_{k=1^n}(1/k^m) *)
ES[m_,n_]:= ((-1)^(m-1)/(m-1)!)*( PolyGamma[m-1,n+1] - PolyGamma[m-1,1] );
(* m=1,2,3,... ; Note: PolyGamma[0,x]=PolyGamma[x] *)
 
(* Now functions Z[m,n]=Zm(n) = (1/n!) (d^m/dx^m) (Gamma[n+1+x]/Gamma[1+x])|_{x=0} : *)
Z1[n_]:=ES[1,n];
Z2[n_]:= ES[1,n]^2 - ES[2,n];
Z3[n_]:=ES[1,n]^3 - 3 ES[1,n] ES[2,n] + 2 ES[3,n];
Z4[n_]:=ES[1,n]^4 - 6 ES[1,n]^2 ES[2,n] + 3 ES[2,n]^2 + 8 ES[1,n] ES[3,n] - 6 ES[4,n];
(* The above Z-ES relations I obtained by my own calculations *) 
(*--------------------------------------------------------*)
 k1f[s_]:=ce1 s (-1 + HarmonicNumber[s]);

k2f[s_]:= ( ((s-1) s/2) c2 + ce1^2( (s(s+1)/2)( (PolyGamma[s+1]-PolyGamma[1])^2 + (PolyGamma[1,s+1]-PolyGamma[1,1]) + 1) - s^2(PolyGamma[s+1]-PolyGamma[1])) );

k3f[s_]:= ( (s(s-1)(s+4)/12) c3 + (s(s+1)(s+2)/2) ce1 c2 ( (s/(s+2)) ES[1,s+1] + (2/(s+1)) + (1/(s+2)) - (11/6) ) + (s(s+1)(s+2)/6) ce1^3 ( Z3[s+2] - 3 Z2[s+1] + 3 ES[1,s] - 1 ) );


tk1f[s_]:=-k1f[s];
tk2f[s_]:=k1f[s] k1f[s+1]-k2f[s];
tk3f[ s_]:=-k1f[s] k1f[s+1] k1f[s+2] + k1f[s] k2f[s+1] + k1f[s+2] k2f[s] - k3f[s];

(*----------------------------------------------------------------------*)
   (* Now calculation of A(nu,Q^2)=A_{nu}(Q^2) which are the analogs of the underlying pQCD powers (a(Q^2))^nu = (alpha_s(Q^2)/Pi)^nu: *) 
  (* --------------------------------------------------_*)
   (* The 0th order: A0o(nu,Q^2)=A_{nu}(Q^2}^{(0)}={tilde A}_{nu}(Q^2) *)
   (* The 1st order:  A1o(nu,Q^2)=A_{nu}(Q^2}^{(1)}={tilde A}_{nu}(Q^2) + {tilde k}_1(nu) {tilde A}_{nu+1}(Q^2) *)
   (* The 2nd order:  A2o(nu,Q^2)=A_{nu}(Q^2}^{(2)}={tilde A}_{nu}(Q^2) + {tilde k}_1(nu) {tilde A}_{nu+1}(Q^2) + {tilde k}_2(nu) {tilde A}_{nu+2}(Q^2)*)
   (* The 3rd order:  A3o(nu,Q^2)=A_{nu}(Q^2}^{(3)}={tilde A}_{nu}(Q^2) + Sum{s=1,2,3} {tilde k}_s(nu) {tilde A}_{nu+s}(Q^2) *)

A0o[nupp_,Q2_,fi_]:=tA[nupp,Q2,fi];
A1o[nupp_,Q2_,fi_]:=tA[nupp,Q2,fi] + tk1f[nupp] tA[1+nupp,Q2,fi];
A2o[nupp_,Q2_,fi_]:=tA[nupp,Q2,fi] + tk1f[nupp] tA[1+nupp,Q2,fi] + tk2f[nupp] tA[2+nupp,Q2,fi];
A3o[nupp_,Q2_,fi_]:=tA[nupp,Q2,fi] + tk1f[nupp] tA[1+nupp,Q2,fi] + tk2f[nupp]  tA[2+nupp,Q2,fi]+ tk3f[nupp] tA[3+nupp,Q2,fi];
  (*-----------------------------------------------------------------*)
   Nf=nqu;
   (* Save["4l3dAQCDcouplrt0XXXal0YYYYN.m", {s0, s1s0, s2s0, s3s0, f1ss0, f2ss0, f3ss0, M0, LLamb, L2, Q2max, rtauRenM, be0, ce1, c2, c3,c4, c5, ce2MS,ce3MS, bt, betaMS5l, beta, Nf, mu2in,  alphaMSMZ, ainMS, ain, apert, tr1, tr1pl, A1, tA2, tA3, tA4, tA, A0o, A1o, A2o, A3o, aMSnf3}] *)
  (* Save["AdlerFunction4lMiniMOM.m", {nf, d1bar,d2bar,d3bar,td1bar,td2bar,td3bar,td1,td2,td3,d1,d2,d3,alt,kapps}] *)
  (*-----------------------------------------------------------------*)
    (*-----------------------------------------------------------------*)
    (*-----------------------------------------------------------------*)
    (*-----------------------------------------------------------------*)
    (*-----------------------------------------------------------------*)
    (*-----------------------------------------------------------------*)
    (*-----------------------------------------------------------------*)
    (*-----------------------------------------------------------------*)
    (*-----------------------------------------------------------------*)
  (* Now some optional plots: *)
  (* The following taken from prg. 2del3lPadbetc2.m and adjusted *)

   (* First I have to simulate the delta functions, with a trick: ver large values are not displayed and appear as values zero: *)
(* This simulates a peak (pole) of r(s) at s=-uL: fuL[s] = 10.^2 (invisible) for s<-uL; =0 for -uL < s  *)
  fdels1s2s3[s_]:=If[s>s3s0,0,If[s>s2s0,10,If[s>s1s0,0,10]]]; (* This simulates a peaks at s=s3s0, s2s0,s1s0; note: s1 < s2 < s3 *)
  
  
 rho1model[sig_] := If[sig/L2 > s0, tr1pl[sig/L2],fdels1s2s3[sig/L2]];
  rho1aux[sig_]:= If[sig/L2 > s0, tr1pl[sig/L2],0];

 (* rho1model[sig_] := If[sig > s0 L2, tr1[sig/L2], (Pi f1ss0)/(Sqrt[4 Pi tteps1]) Exp[-(sig/L2 - s1s0)^2/(4 tteps1)] + (Pi f2ss0)/(Sqrt[4 Pi tteps2]) Exp[-(sig/L2 - s2s0)^2/(4 tteps2)]]; *)
(*rho1model[sigma] is the rho_ 1(sigma=s L2) of our anQCD model;tteps->0,i.e.,(Pi f1o2s0)/(Sqrt[4 Pi tteps]) Exp[-(s-s1s0)^2/(4 tteps)] is an approximation for Pi f1o2[s0] delta (s-s1[s0])*)


   (* Note: r1perttot[sig] = tr1[sig/L2] + Pi Res delta(sig/L2 + uL), i.e., the true perturbative rho_1(sig), including the singularity at sig=-uL L2: *)
(*   r1perttot[sig_]:=tr1[sig/L2] + (Pi Res)/(Sqrt[4 Pi ttepsuL]) Exp[-(sig/L2 + uL)^2/(4 ttepsuL)]; *)
tr1[s_] := If[s > 0, Im[a[s, -Pi]], Im[a[-s, -10.^-8]]];
   r1perttot[sig_]:=tr1[sig/L2];

rho1MMFiga = 
 rho1modFiga = 
  Plot[{r1perttot[sig]}, {sig, -1.5, 2.0}, 
   PlotRange -> {{-1.5, 2.0}, {0., 0.23}}, AxesOrigin -> {0., 0.}, 
   Frame -> True, 
   FrameLabel -> {\[Sigma] "  "[Superscript[GeV, 2]], 
     Subscript[\[Rho], a]}, PlotStyle -> {Black}, 
   Epilog -> {Inset[Style["(a)", 15], Scaled[{0.8, 0.94}]],(*Inset[
     Style[Subscript[c,2]\[Equal]c2,11],Scaled[{0.19,0.95}]],*)
     Inset[Style[Subscript[\[Rho], a], 10], 
      Scaled[{0.7, 0.35}]], 
     Inset[Style[Subscript[\[Rho], a], 10], 
      Scaled[{0.2, 0.6}]], 
     Inset[Style[
       Superscript[Subscript[\[Alpha], s], Style[OverBar[MS], 6]][
         Subscript[M, Z]] == alphaMSMZ, 11], Scaled[{0.69, 0.80}]], 
     Inset[Style[
       Subscript[N, f] == Nf ", " Subscript[\[CapitalLambda], L] == 
        NumberForm[LLamb, {5, 4}] GeV, 11], Scaled[{0.75, 0.70}]], 
   Inset[Style["MiniMOM", 11], Scaled[{0.62, 0.60}]]}];
  
line1 = Line[{{Log[M0^2], 0}, {Log[M0^2], r1perttot[M0^2]}}];
rho1MMFigb = 
 LogLinearPlot[{rho1model[sig], r1perttot[sig]}, {sig, 0.01, 100000}, 
  PlotRange -> {{0.01, 10000.}, {0., 0.12}}, AxesOrigin -> {0.01, 0.},
   Frame -> True, 
  FrameLabel -> {\[Sigma] "  "[Superscript[GeV, 2]], 
    Subscript[\[Rho], A]}, PlotStyle -> {Black, {Black, Dotted}}, 
    Epilog -> {{Directive[Black,Thick], line1}, Inset[Style["(b)", 15], Scaled[{0.8, 0.94}]],
	    Inset[Superscript[Subscript[\[Alpha], s], Style[OverBar[MS], 6]][
       Subscript[M, Z]] == alphaMSMZ, Scaled[{0.75, 0.80}]], 
    Inset[Subscript[N, f] == 
      3 ", " Subscript[\[CapitalLambda], L] == 
      NumberForm[LLamb, {5, 4}] GeV, Scaled[{0.77, 0.70}]], 
    Inset["MiniMOM", Scaled[{0.67, 0.62}]], 
    Inset[Subscript[s, 0] == s0 ", " Subscript[s, 1] == s1s0, 
     Scaled[{0.75, 0.56}]], 
    Inset[Style[Subsuperscript[M, 0, 2], 13], Scaled[{0.525, 0.06}]], 
    Inset[Style[Subsuperscript[M, 3, 2], 13], Scaled[{0.43, 0.06}]], 
    Inset[Style[Subsuperscript[M, 2, 2], 13], Scaled[{0.27, 0.06}]], 
    Inset[Style[Subsuperscript[M, 1, 2], 13], Scaled[{0.085, 0.06}]], 
    Inset[Style[Subscript[\[Rho], a], 10], Scaled[{0.35, 0.61}]], 
    Inset[Style[Subscript[\[Rho], a], 10], Scaled[{0.7, 0.20}]], 
    Inset[Style["(=" Subscript[ \[Rho], A], 10], 
     Scaled[{0.78, 0.20}]], 
    Inset[Style[")", 10], Scaled[{0.84, 0.20}]], 
    Inset[Style[Subscript[\[Rho], A], 10], Scaled[{0.10, 0.93}]], 
    Inset[Style[Subscript[\[Rho], A], 10], Scaled[{0.26, 0.9}]], 
    Inset[Style[Subscript[\[Rho], A], 10], Scaled[{0.45, 0.9}]]}];


rho1MMFigbaux = 
 LogLinearPlot[rho1aux[sig], {sig, 0.01, 100000.}, 
  PlotRange -> {{0.01, 100000.}, {0., 0.12}}, 
  AxesOrigin -> {0.01, 0.}, Frame -> True, 
  FrameLabel -> {\[Sigma] "  "[Superscript[GeV, 2]], 
    Subscript[\[Rho], A]}, PlotStyle -> Black];

rho1MMFigbComb=Show[rho1MMFigb,rho1MMFigbaux];


(* Export["rho1MMFigart0201al01181N.pdf",rho1MMFiga]; Export["rho1MMFigbrt0201al01181N.pdf",rho1MMFigbComb]; *)
(*-------------------------------------------------------*)
(* Calculation of rtau in MSbar scheme: *)
(* A shortcut for this part: <<4l3danQCDcouplrt0XXXal0YYYY.m; << AdlerFunction4lMiniMOM.m; and then execute the parts below: *)

mtau2= 1.77686^2;
(* Note: for ADLER case, use: mtau2=2.80 *)
mxpts=1000;

 aMSs0=aMSnf3[mtau2];  (* Pi aMSnf3[mtau2]=0.319, which agrees with the worls average, I checked, also Boito et al get this <-> Pi aMSnf5[MZ^2]=0.1185 *)

		   (* Now we integrate MSbar coupling along the circle of radius |Q^2|=s0=mtau2 in the complex Q^2-plane: aMScomplexs0[fi] = aMS[Q^2=s0 Exp[I fi]] (x=|Q^2| in GeV^2; fi angle in radians) *)

eeps=10.^-6;
aabeMScomplex=NDSolve[{yy'[fi]==I betaMS4l[yy[fi]],yy[0]==aMSs0},yy,{fi,0,Pi-eeps}];  
(*'*)
aMScomppluss0[fi_]:=yy[fi]/.aabeMScomplex[[1]]; 

aMScomplexs0[fi_]:=If[fi >=0,aMScomppluss0[fi],Conjugate[aMScomppluss0[-fi]]];

(* aMScomplexs0[fi] = aMS[Q^2=s0 Exp[i fi]]; -Pi+eeps < fi < Pi - eeps *)
(* s0=mtau^2 here *)
(* This are the barMS coupling values on the circle Q^2 = s0 Exp[I fi] *)

tauMSp1 = 2 (1/2/Pi) NIntegrate[
    Re[(1 + Exp[I fi])^3 (1 - Exp[I fi]) (aMScomplexs0[fi])], {fi, 
     0., +Pi - eeps}, MaxPoints -> mxpts];
tauMSp2 = 2 (1/2/Pi) NIntegrate[
    Re[(1 + Exp[I fi])^3 (1 - Exp[I fi])  (aMScomplexs0[fi])^2], {fi, 
     0., +Pi - eeps}, MaxPoints -> mxpts];
tauMSp3 = 2 (1/2/Pi) NIntegrate[
    Re[(1 + Exp[I fi])^3 (1 - Exp[I fi])  (aMScomplexs0[fi])^3], {fi, 
     0., +Pi - eeps}, MaxPoints -> mxpts];
tauMSp4 = 2 (1/2/Pi) NIntegrate[
    Re[(1 + Exp[I fi])^3 (1 - Exp[I fi])  (aMScomplexs0[fi])^4], {fi, 
     0., +Pi - eeps}, MaxPoints -> mxpts];

rtauMSp1 = tauMSp1;
rtauMSp2 = tauMSp2 d1bar;
rtauMSp3 = tauMSp3 d2bar;
rtauMSp4 = tauMSp4 d3bar;

rtauMSpTPS = rtauMSp1 + rtauMSp2 + rtauMSp3 + rtauMSp4;
rtauMSp1
rtauMSp2
rtauMSp3
rtauMSp4
rtauMSpTPS
  (*--------------------------------------------------*)
 (*--------------------------------------------------*)
 (*--------------------------------------------------*)
(* V-channel Adler function (Nf=3) -- theory vs experiment: *)
  (* a shortcut: start from here, call: *)
  (* <<4l3danQCDcouplrt0XXXal0YYYY.m; <<AdlerFunction4lMiniMOM.m; *)
  (* and then continue below: *)

DdBG22[Q2_,ffi_]:=1 + alt[2,1] A1[kapps[2,1] Q2,ffi] + 	 alt[2,2] A1[kapps[2,2] Q2,ffi];
(* This is D[Q^2=Q2 Exp[i ffi]] = 1 + d[Q^2=Q2 Exp[i ffi]], the full nf=3 massless Adler function in the complex plane -- the dBG[2/2] resummed version, in Lambert MimiMOM scheme  *)

DTPS[Q2_,ffi_]:=1 + A1[Q2,ffi] + td1 tA2[Q2,ffi]+ td2 tA3[Q2,ffi]+ td3 tA4[Q2,ffi];
(* This is D[Q^2=Q2 Exp[i ffi]] = 1 + d[Q^2=Q2 Exp[i ffi]], the full nf=3 massless Adler function in the complex plane -- the TPS version, with RScl mu^2=Q^2,  in Lambert MimiMOM scheme*)

DMSbarTPS[Q2_]:=1 + aMSnf3[Q2] + d1bar aMSnf3[Q2]^2  + d2bar aMSnf3[Q2]^3 + d3bar aMSnf3[Q2]^4;
(* This is D[Q^2=Q2] = 1 + d[Q^2], the full nf=3 massless Adler function for positive Q^2=Q2, in MSbar scheme, TPS *) 
(*-----------*)
(* Now higher twist (HT)terms, AQCD (in Lambert MiniMOM) + OPE: *)
fPi=0.13050; (* This value UPDATED on March 22, 2017; it is from PDG, Rosner et al, Leptonic decays of cherged pseudoscalar mesons, revised March 2016, [http://pdg.lbl.gov/2016/reviews/rpp2016-rev-pseudoscalar-meson-decay-cons.pdf], their Eq.(6), and PDG 2016 value Vud=0.97417 +- 0.00021, giving fPi=(0.12713/0.97417)=0.130501 GeV, with uncertainties +-0.00014 GeV *) 
mPi=0.13957; (* This value UPDATED March 22, 2017; PDG2014 and PDG2016 *)
(*-----------*)
(*--------------------------*)(* HT extracted from OPAL data, sigmax=3.136 GeV^2;  the case 4l3drt00201al01181: *)
(* First V+A channel: *)
(* AQCD: <aGG>OPAL =-0.0010349, from OPAL, V+A channel *)
aGGAQCDOPAL=-0.00103498;
O4AQCDOPAL=-0.0005042; (* =(1/6)aGGAQCDOPAL-fPi^2 mPi^2; *)
O6AQCDOPAL=0.00112; (* OPAL extraction fit, V+A channel *)
(* Note, the above values are for V+A channel *)
(*---*)
(* Now going from V+A to 2 V (because we want to represent D(Q^2)_{2 V}, not D(Q^2)_{V+A}  (2 V, not V, because we want to keep the same normalization: D(Q^2)_{2 V} = 1 + A + O(A2)  as is D(Q^2)_{V+A} *)
O4AQCDOPAL2V=O4AQCDOPAL; (* we know that O4_{V} = O4_{A} -> O4_{2V} = O4_{V+A} *)
O6AQCDOPAL2V=-(7/2) O6AQCDOPAL; (* this follows from factorization hypothesis, see Ioffe his PPNP56,232(2006) review *)	      
(*---*)			      
HTOPAL2V[Q2_,ffi_]:= 2 Pi^2 ( 2 O4AQCDOPAL2V/(Q2 Exp[I ffi])^2 +3 O6AQCDOPAL2V/(Q2 Exp[I ffi])^3 );
(*-----------*)
(* Now for MSbar pQCD approach (+OPE) (the case rt0201al01181): *)
aGGMSbarOPAL=+0.00694; (* OPAL extraction, in MSbar pQCD approach, for V+A channel *)
O4MSbarOPAL=+0.000825; (* (1/6)aGGMSbarOPAL-fPi^2 mPi^2; *)
O6MSbarOPAL=-0.00165; (* OPAL extraction fit, for V+A channel *)
O4MSbarOPAL2V=O4MSbarOPAL;
O6MSbarOPAL2V=(-7/2) O6MSbarOPAL;		      
HTMSbarOPAL2V[Q2_]:= 2 Pi^2 ( 2 O4MSbarOPAL2V/Q2^2 +3 O6MSbarOPAL2V/Q2^3 );
(*----------------------------------------------*)	
(* The leading-twist part, for the considered case (4l3drt0201al01181): *)
(* tabDdBG22o4l3drt0201al01181=Table[{0.5 + 0.1 jj,Re[DdBG22[(0.5+0.1 jj)^2,0]]},{jj,0,25}]; *)
(* Save["tabDdBG22o4l3drt0201al01181.save",tabDdBG22o4l3drt0201al01181]; *)
<<tabDdBG22o4l3drt0201al01181.save;
DdBG22int[Q_]:=Interpolation[tabDdBG22o4l3drt0201al01181][Q];
(*----------------------------------------------------*)
(* Now the experimental curve for D_v(Q^2) -- upper and lower bounds, as functions of x=Q (in GeV}: *)
(* From Fig.1.7 of Nesterenko s book *)

 (* Corrected points from Nesterenkos book, April 12, 207: *)
FigDVexp=Import["scanDV2.pdf"]; (* point0o0={67.81,141.0}; point0o0={64.61,141.2}; point2o0={227.8,141.0}; point0o1={67.81,350.0}; *)
point0o0={64.61,140.4}; point2o0={274.2,140.4};point0o1={64.61,350.};
(* DVUB={{67.81,141.0},{72.81,144.},{75.81,147},{82.81,166.},{92.81,199.},{104.8,249.},{114.8,283.},{122.8,306.},{131.8,327.},{144.8,350.},{156.8,362.9},{172.8,372.9},{189.8,378.9},{214.8,382.9},{239.8,383.9},{270.,383.9},{310.,383.9}}; *)
DVUB={{64.61,140.4},{84.6,179.6},{90.2,202.},{99.8,238.8},{115.,290.8},{134.2,337.2},{142.2,350.},{157.4,365.1},{171.,373.1},{189.4,379.5},{205.4,381.9},{236.6,383.5},{259.8,383.5},{313.4,382.7},{379.8,381.9},{456.6,380.3},{481.,380.3}};
(* DVLB={{67.81,141.0},{72.81,144.},{75.81,147},{84.81,168.},{94.81,200.},{110.8,252.},{120.8,280.},{130.8,302.},{148.8,328.},{166.8,342.},{187.8,350.},{208.8,354.},{255.8,355.},{270.,355.},{310.,355.}}; *)
DVLB={{64.61,140.4},{84.6,174.8},{90.2,193.2},{100.6,227.6},{116.6,274.8},{138.2,318.},{148.6,330.},{163.8,341.2},{173.4,346.},{190.2,350.8},{213.4,354.},{238.2,354.8},{265.4,354.8},{313,354.},{400.,353.2},{444.,353.2},{481.,353.2}};
																			     
(* Now linear trasformation to the actual (true) values: *)
(* xlinf=x = x(true); xp=x(picture); ylinf =y = y(true); yp=y(picture): *)
(* x = kx xp + bx; y = ky yp + by; *)
kx=(2-0)/(274.2-64.61);
ky=(1 - 0)/(350.0-140.4);
bx=-kx 64.61;
by=-ky 140.4;

xlinf[xp_]:=kx xp + bx;
ylinf[yp_]:=ky yp + by;
DVUBLin=Table[{xlinf[DVUB[[jj]][[1]]],ylinf[DVUB[[jj]][[2]]]},{jj,1,17}];
DVLBLin=Table[{xlinf[DVLB[[jj]][[1]]],ylinf[DVLB[[jj]][[2]]]},{jj,1,17}];

DVUBLinf=Interpolation[DVUBLin];
DVLBLinf=Interpolation[DVLBLin];
(* End of experimental curve part *)
(*---------------------------*)
figAdlerDVvsQOPALandTheorOPE4l3drt0201al01181 = 
Plot[{DdBG22int[Q] + HTOPAL2V[Q^2, 0], 
    DMSbarTPS[Q^2] + HTMSbarOPAL2V[Q^2], DVUBLinf[Q], 
   DVLBLinf[Q]}, {Q, 0.9, 2.}, PlotRange -> {{0.9, 2.}, {0., 2.}}, 
  AxesOrigin -> {0.9, 0.}, Filling -> {3 -> {4}}, Frame -> True, 
  FrameLabel -> {Q[GeV], Subscript[D, V][Superscript[Q, 2]]}, 
  PlotStyle -> { {Black, Thickness[0.003]}, 
    DotDashed, {Brown, Thin}, {Brown, Thin}}];

(* Save["figAdlerDVvsQOPALandTheorOPE4l3drt0201al01181.save",figAdlerDVvsQOPALandTheorOPE4l3drt0201al01181]; *)

 
