(* ::Package:: *)
(*
EKO NLO helper skript
=====================

Before running
--------------
- Compile FormGet (if necessary) (some instructions are given below)
- Adjust the path to FormGet below
- Download the file hepres.h and hepres2.h from A. Vogt
- Strip the "L" in this files
- make sure a directory "ad_nlo" exists - here the output will be stored
- run this file by bringing up Mathematica and then issue "<<NLO.m"

What it does
------------
This script brings the next-to-leading order non-singlet [MVV04] and singlet [VMV04] expresisons
into an executable form. A. Vogt provides the FORM expressions of their computations on his
website: https://www.liverpool.ac.uk/~avogt/split.html. From these files the FORM variable
definition "L" has to be stripped in order to be readable for FormGet [HPS16] as it expects an
FORM *output* file.

References
----------
[MVV04] S. Moch, J. A. M. Vermaseren, and A. Vogt. The Three loop splitting functions in QCD: The Nonsinglet case.
        Nucl. Phys., B688:101–134, 2004. arXiv:hep-ph/0403192, doi:10.1016/j.nuclphysb.2004.03.030.
[VMV04] A. Vogt, S. Moch, and J. A. M. Vermaseren. The Three-loop splitting functions in QCD: The Singlet case.
        Nucl. Phys., B691:129–181, 2004. arXiv:hep-ph/0404111, doi:10.1016/j.nuclphysb.2004.04.024.
[HPS16] Thomas Hahn, Sebastian Paßehr, and Christian Schappacher. FormCalc 9 and Extensions.
        PoS, LL2016:068, 2016. arXiv:1604.04611, doi:10.1088/1742-6596/762/1/012065.

*)

(* add the formal harmonic sums definitions to allow to evaluate the raw input  *)
S[{m_},M_]:=Sum[Sign[m]^k/k^Abs@m,{k,1,M}]
S[{m_,m1_},M_]:=Sum[Sign[m]^k/k^Abs@m S[{m1},k],{k,1,M}]

(*
In order to compile the FormGet executable - adjust and uncomment the following two statements

Needs["CCompilerDriver`"];
CreateExecutable[{"/home/felix/Physik/Mathematica/FormGet.tm"}, "FormGet", "ShellOutputFunction" :> Print,  "ShellCommandFunction" :> Print]
*)

(*
Load raw data
- adjust the path to FormGet
- strip the L from the Form files
*)
Install["/home/felix/Physik/Mathematica/FormGet"]
hepres = FormGet["hepres.h"];
hepres2 = FormGet["hepres2.h"];

(* get expression from hepres *)
getRawNS[name_,coeff_]:=(Coefficient[name/.hepres//.{N->n,delta@_->0,theta@_->1,den@z_:>1/z},coeff])/.Br->Identity
getRawS[name_,coeff_]:=(Coefficient[name/.hepres2//.{N->n,delta@_->0,theta@_->1,den@z_:>1/z},coeff])/.Br->Identity
parseRaw[expr_]:=expr/.{R->List};
(* undo Fortran in Pegasus *)
parsePegasus[expr_]:=Module[{b},
    b=expr;
    b=b/.{SSCHLM->SLC-SLV,
    SSCHLP->SLC+SLV,
    SSTR2M->ZETA@2-DPSI [N1/2,1],
    SSTR3M->1/2*DPSI [N1/2,2]+ZETA@3,
    SSTR2P->ZETA@2-DPSI [N2/2,1],
    SSTR3P->1/2*DPSI [N2/2,2]+ZETA@3
    };
    b=b/.{SLC->-5/8*ZETA@3,
    SLV->-ZETA@2/2*(PSI[N1/2]-PSI[N/2])+S1/NS+SPMOM};
    b=b/.{
    N->n,NS->n^2,NT->n^3,NFO->n^4,NFI->n^5,NSI->n^6,NSE->n^7,NE->n^8,NN->n^9,
    N1->n+1,N1S->(n+1)^2,N1T->(n+1)^3,
    N2->n+2,N2S->(n+2)^2,N2T->(n+2)^3,
    NM->n-1,NMS->(n-1)^2
    };
    b=b/.{ZETA->Zeta,PSI@z_:>PolyGamma[z],DPSI[z_,k_]:>PolyGamma[k,z]};
    b=b/.{S1->S[{1},n],S2->S[{2},n]};
    b
];
(* my transformations inspired by Muselli, though he still has errors *)
transformMe[expr_]:=Module[{me},
    me=expr;
    me=me/.{S[R[1,2],n]->-S[R[2,1],n]+S[R@3,n]+S[R@2,n]S[R@1,n]}//Expand;
    me=me/.{S[R[1,-2],n]->-S[R[-2,1],n]+S[R[-3],n]+S[R[-2],n]S[R@1,n]}//Expand;
    me=me/.{S[R[-2,1],n]->-5/8zeta@3+ zeta@2(S[R[-1],n]-\[Eta]@n/n+Log@2)+\[Eta][n]*(S[R@1,n]/n^2+g3@n)}//Expand;
    me=me/.FunctionExpand@Flatten@Table[S[If[Head@r===List,R@@r,R@r],k+n]->S[If[Head@r===List,r,{r}],k+n]-S[If[Head@r===List,r,{r}],n]+S[If[Head@r===List,R@@r,R@r],n],{r,{1,2,{1,1}}},{k,{2,1,-1}}];
    me=me/.{S[R[-2],n+2]->\[Eta]@n/(n+2)^2+S[R[-2],n+1],S[R[-2],n-1]->-\[Eta]@n/(n)^2+S[R[-2],n]}/.{S[R[-2],n+1]->-\[Eta]@n/(n+1)^2+S[R[-2],n]};
    me=me/.{S[R[1,1],n]->1/2(S[R@1,n]^2+S[R@2,n])};
    me=me/.{EulerGamma->S[R@1,n]-S[{1},n]+EulerGamma};
    me=me/.Table[S[R[-l],n]->2^(1-l)((1+\[Eta]@n)/2S[R[l],n/2]+(1-\[Eta]@n)/2S[R[l],(n-1)/2])-S[R@l,n],{l,3}];
    me
];
(* make my expression executable - preserving eta *)
nMe[expr_]:=Module[{me},
    me=expr;
    me=me/.{S[R[-1],n]->\[Eta]@n*1/2(PolyGamma[0,(1+n+1)/2]-PolyGamma[0,(1+n)/2])-Log@2};
    me=me/.{S[R[-2],n]->-\[Eta]@n*1/4(PolyGamma[1,(1+n+1)/2]-PolyGamma[1,(1+n)/2])-1/2zeta@2};
    me=me/.{S[R[-3],n]->(\[Eta]@n*1/2*1/8(PolyGamma[2,(1+n+1)/2]-PolyGamma[2,(1+n)/2])-3/4zeta@3)};
    me=me/.{R->List,zeta->Zeta};
    me
];
(* write my expression to file *)
exportMe[expr_,fp_]:=Module[{me,out},
    me=expr;
    me=Collect[me,{g3@_,zeta@_,S[R@_,_]},FullSimplify];
    me=me/.{zeta@2->zeta2,zeta@3->zeta3,g3@n->g3n};
    me=me/.{S[R@1,n]->S1,S[R@2,n]->S2,S[R@3,n]->S3};
    me=me/.{S[R@1,(n-1)/2]->Sp1m,S[R@2,(n-1)/2]->Sp2m,S[R@3,(n-1)/2]->Sp3m};
    me=me/.{S[R@1,n/2]->Sp1p,S[R@2,n/2]->Sp2p,S[R@3,n/2]->Sp3p};
    (* write *)
    Print["writing to ",fp];
    out=OpenWrite[fp];
    WriteString[out,ToString@CForm@me];
    Close[out];
];
(* evaluate *)
eval[expr_,nn_]:=Module[{g3n},
    g3n = Integrate[x^(nn-1)PolyLog[2,x]/(1+x),{x,0,1}];
    expr/.{g3@n->g3n,SPMOM->g3n,n->nn}
];

(*
Non-singlet sector
------------------
*)

(* ns+- difference *)
Module[{spmom,rawMCF2,rawMCFCA,rawPCF2,rawPCFCA,MochMCF2,MochMCFCA,MochPCF2,MochPCFCA,PNMA,PNPA,PNSB,PegasusMA,PegasusPA,PegasusB,PegasusMCFCA,PegasusPCFCA,meMCF2,meMCFCA,mePCF2,mePCFCA,nmeMCF2,nmeMCFCA,nmePCF2,nmePCFCA},
    spmom[n_]:= Integrate[x^(n-1)PolyLog[2,x]/(1+x),{x,0,1}];
    (* load raw *)
    rawMCF2=getRawNS[gqq1m,cf^2];
    rawMCFCA=getRawNS[gqq1m,cf ca];
    rawPCF2=getRawNS[gqq1p,cf^2];
    rawPCFCA=getRawNS[gqq1p,cf ca];
    MochMCF2=parseRaw@rawMCF2;
    MochMCFCA=parseRaw@rawMCFCA;
    MochPCF2=parseRaw@rawPCF2;
    MochPCFCA=parseRaw@rawPCFCA;
    (* load Pegasus *)
    PNMA=(16*S1*(2*N+1)/(NS*N1S)+ 
    16*(2*S1-1/(N*N1))*(S2-SSTR2M)+ 
    64*SSCHLM+24*S2-3-8*SSTR3M- 
    8*(3*NT+NS-1)/(NT*N1T)+
    16*(2*NS+2*N+1)/(NT*N1T)
    )*(-1/2);
    PNPA=(16*S1*(2*N+1)/(NS*N1S)+16*(2*S1-1/(N*N1))*(S2-SSTR2P)+64*SSCHLP+24*S2-3-8*SSTR3P-8*(3*NT+NS-1)/(NT*N1T)-16*(2*NS+2*N+1)/(NT*N1T))*(-1/2);
    PNSB=(S1*(536/9+8*(2*N+1)/(NS*N1S))- (16*S1+52/3-8/(N*N1))*S2-43/6-(151*NFO+263*NT+97*NS+3*N+9)*4/(9*NT*N1T))*(-1/2);
    PegasusMA=-parsePegasus@PNMA;
    PegasusPA=-parsePegasus@PNPA;
    PegasusB=-parsePegasus@PNSB;
    PegasusMCFCA=PegasusB-PegasusMA/2;
    PegasusPCFCA=PegasusB-PegasusPA/2;
    (* load me *)
    meMCF2=transformMe@rawMCF2;
    meMCFCA=transformMe@rawMCFCA;
    mePCF2=transformMe@rawPCF2;
    mePCFCA=transformMe@rawPCFCA;
    Print@Collect[meMCF2/.{\[Eta]@_:>-1},{g3@_,zeta@_,S[R@_,_]},FullSimplify];
    Print@Collect[meMCFCA/.{\[Eta]@_:>-1},{g3@_,zeta@_,S[R@_,_]},FullSimplify];
    Print@Collect[mePCF2/.{\[Eta]@_:>1},{g3@_,zeta@_,S[R@_,_]},FullSimplify];
    Print@Collect[mePCFCA/.{\[Eta]@_:>1},{g3@_,zeta@_,S[R@_,_]},FullSimplify];
    exportMe[meMCF2/.{\[Eta]@_:>-1},"ad_nlo/gqq1m_cfcf.c"];
    exportMe[meMCFCA/.{\[Eta]@_:>-1},"ad_nlo/gqq1m_cfca.c"];
    exportMe[mePCF2/.{\[Eta]@_:>1},"ad_nlo/gqq1p_cfcf.c"];
    exportMe[mePCFCA/.{\[Eta]@_:>1},"ad_nlo/gqq1p_cfca.c"];
    nmeMCF2=nMe@meMCF2;
    nmeMCFCA=nMe@meMCFCA;
    nmePCF2=nMe@mePCF2;
    nmePCFCA=nMe@mePCFCA;
    (* check I match to Pegasus *)
    Print[FullSimplify[PegasusMA-nmeMCF2/.{SPMOM->g3@n,\[Eta]@_:>-1}]];
    Print[FullSimplify[PegasusMCFCA-nmeMCFCA/.{SPMOM->g3@n,\[Eta]@_:>-1}]];
    Print[FullSimplify[PegasusPA-nmePCF2/.{SPMOM->g3@n,\[Eta]@_:>1}]];
    Print[FullSimplify[PegasusPCFCA-nmePCFCA/.{SPMOM->g3@n,\[Eta]@_:>1}]];
    (* get some reference numbers *)
    Print@FullSimplify@eval[nmePCF2/.{\[Eta]@_->1},2];
    Print@FullSimplify@eval[nmePCFCA/.{\[Eta]@_->1},2];
    (* compare all 3 expresisons *)
    FullSimplify@Table[With[{g3n=spmom@n},{n,
    PegasusMA-MochMCF2/.{SPMOM->g3n},
    Identity[PegasusMA-nmeMCF2/.{SPMOM->g3n,g3@n->g3n,\[Eta]@_:>-1}],
    Identity[PegasusMCFCA-MochMCFCA/.{SPMOM->g3n}],
    Identity[PegasusMCFCA-nmeMCFCA/.{SPMOM->g3n,g3@n->g3n,\[Eta]@_:>-1}],
    PegasusPA-MochPCF2/.{SPMOM->g3n},
    Identity[PegasusPA-nmePCF2/.{SPMOM->g3n,g3@n->g3n,\[Eta]@_:>1}],
    Identity[PegasusPCFCA-MochPCFCA/.{SPMOM->g3n}],
    Identity[PegasusPCFCA-nmePCFCA/.{SPMOM->g3n,g3@n->g3n,\[Eta]@_:>1}]
    }],{n,1,7}]
]

(* ns+- shared part *)
Module[{rawMCFnf,rawPCFnf,MochMCFnf,MochPCFnf,PNSC,PegasusC},
    rawMCFnf=getRawNS[gqq1m,cf nf];
    rawPCFnf=getRawNS[gqq1p,cf nf];
    exportMe[rawMCFnf,"ad_nlo/gqq1m_cfnf.c"];
    exportMe[rawPCFnf,"ad_nlo/gqq1p_cfnf.c"];
    MochMCFnf=parseRaw@rawMCFnf;
    MochPCFnf=parseRaw@rawPCFnf;
    PNSC=(-160/9*S1+32/3*S2+4/3+16*(11*NS+5*N-3)/(9*NS*N1S))*(-1/2);
    PegasusC=-parsePegasus@PNSC;
    Table[{n,MochMCFnf,MochMCFnf-PegasusC/2,MochPCFnf-PegasusC/2},{n,6}]
]

(*
Singlet sector
--------------
*)

(* pure singlet *)
Module[{rawNFCF,MochNFCF,PPSA,PegasusNFCF,meNFCF,nmeNFCF},
    (* load raw *)
    rawNFCF=getRawS[gqqps1,nf cf];
    MochNFCF=parseRaw@rawNFCF;
    (* load Pegasus *)
    PPSA=(5*NFI+32*NFO+49*NT+38*NS+28*N+8)/(NM*NT*N1T*N2S)*2;
    PegasusNFCF=-1/2*4*parsePegasus@PPSA;
    (* load me *)
    meNFCF=transformMe@rawNFCF;
    nmeNFCF=nMe@meNFCF;
    Print@Collect[meNFCF/.{\[Eta]@_:>1},{S[R@__,_]},FullSimplify];
    exportMe[meNFCF/.{\[Eta]@_:>1},"ad_nlo/gqqps1_nfcf.c"];
    (* check I match to Pegasus *)
    Print[FullSimplify[PegasusNFCF-nmeNFCF/.{\[Eta]@_:>1}]];
    (* get some reference numbers *)
    Print@FullSimplify[nmeNFCF/.{\[Eta]@_->1,n->2}];
    (* compare all 3 expresisons *)
    FullSimplify@Table[With[{},{n,
    PegasusNFCF-MochNFCF,
    Identity[PegasusNFCF-nmeNFCF/.{\[Eta]@_:>1}]
    }],{n,2,7}]
]

(* quark-gluon *)
Module[{rawNFCF,rawNFCA,MochNFCF,MochNFCA,PQGA,PQGB,PegasusNFCF,PegasusNFCA,meNFCF,meNFCA,nmeNFCF,nmeNFCA},
    (* load raw *)
    rawNFCF=getRawS[gqg1,nf cf];
    rawNFCA=getRawS[gqg1,nf ca];
    MochNFCF=parseRaw@rawNFCF;
    MochNFCA=parseRaw@rawNFCA;
    (* load Pegasus *)
    PQGA=(-2*S1*S1+2*S2-2*SSTR2P)*(NS+N+2)/(N*N1*N2)+(8*S1*(2*N+3))/(N1S*N2S)+2*(NN+6*NE+15*NSE+25*NSI+36*NFI+85*NFO+128*NT+104*NS+64*N+16)/(NM*NT*N1T*N2T);
    PQGB=(2*S1*S1-2*S2+5)*(NS+N+2)/(N*N1*N2)-4*S1/NS+(11*NFO+26*NT+15*NS+8*N+4)/(NT*N1T*N2);
    PegasusNFCA=-1/2*4*parsePegasus@PQGA;
    PegasusNFCF=-1/2*4*parsePegasus@PQGB;
    (* load me *)
    meNFCF=transformMe@rawNFCF;
    meNFCA=transformMe@rawNFCA;
    nmeNFCF=nMe@meNFCF;
    nmeNFCA=nMe@meNFCA;
    Print@Collect[meNFCF/.{\[Eta]@_:>1},{S[R@__,_]},FullSimplify];
    Print@Collect[meNFCA/.{\[Eta]@_:>1},{S[R@__,_]},FullSimplify];
    exportMe[meNFCF/.{\[Eta]@_:>1},"ad_nlo/gqg1_nfcf.c"];
    exportMe[meNFCA/.{\[Eta]@_:>1},"ad_nlo/gqg1_nfca.c"];
    (* check I match to Pegasus *)
    Print[FullSimplify[PegasusNFCF-nmeNFCF/.{\[Eta]@_:>1}]];
    Print[FullSimplify[PegasusNFCA-nmeNFCA/.{\[Eta]@_:>1}]];
    (* get some reference numbers *)
    Print@FullSimplify[nmeNFCF/.{\[Eta]@_->1,n->2}];
    Print@FullSimplify[nmeNFCA/.{\[Eta]@_->1,n->2}];
    (* compare all 3 expresisons *)
    FullSimplify@Table[With[{},{n,
    PegasusNFCA-MochNFCA,
    Identity[PegasusNFCA-nmeNFCA/.{\[Eta]@_:>1}],
    PegasusNFCF-MochNFCF,
    Identity[PegasusNFCF-nmeNFCF/.{\[Eta]@_:>1}]
    }],{n,2,7}]
]

(* gluon-quark *)
Module[{rawCFCF,rawCFCA,rawCFNF,MochCFCF,MochCFCA,MochCFNF,PGQA,PGQB,PGQC,PegasusCFCF,PegasusCFCA,PegasusCFNF,meCFCF,meCFCA,meCFNF,nmeCFCF,nmeCFCA,nmeCFNF},
    (* load raw *)
    rawCFCF=getRawS[ggq1,cf^2];
    rawCFCA=getRawS[ggq1,cf ca];
    rawCFNF=getRawS[ggq1,cf nf];
    MochCFCF=parseRaw@rawCFCF;
    MochCFCA=parseRaw@rawCFCA;
    MochCFNF=parseRaw@rawCFNF;
    (* load Pegasus *)
    PGQA=(-S1*S1+5*S1-S2)*(NS+N+2)/(NM*N*N1)-2*S1/N1S-(12*NSI+30*NFI+43*NFO+28*NT-NS-12*N-4)/(2*NM*NT*N1T);
    PGQB=(S1*S1+S2-SSTR2P)*(NS+N+2)/(NM*N*N1)-S1*(17*NFO+41*NS-22*N-12)/(3*NMS*NS*N1)+(109*NN+621*NE+1400*NSE+1678*NSI+695*NFI-1031*NFO-1304*NT-152*NS+432*N+144)/(9*NMS*NT*N1T*N2S);
    PGQC=4/3*((S1-8/3)*(NS+N+2)/(NM*N*N1)+1/N1S);
    PegasusCFCF=-4*parsePegasus@PGQA;
    PegasusCFCA=-4*parsePegasus@PGQB;
    PegasusCFNF=-1/2*4*parsePegasus@PGQC;
    (* load me *)
    meCFCF=transformMe@rawCFCF;
    meCFCA=transformMe@rawCFCA;
    meCFNF=transformMe@rawCFNF;
    nmeCFCF=nMe@meCFCF;
    nmeCFCA=nMe@meCFCA;
    nmeCFNF=nMe@meCFNF;
    Print@Collect[meCFCF/.{\[Eta]@_:>1},{S[R@__,_]},FullSimplify];
    Print@Collect[meCFCA/.{\[Eta]@_:>1},{S[R@__,_]},FullSimplify];
    Print@Collect[meCFNF/.{\[Eta]@_:>1},{S[R@__,_]},FullSimplify];
    exportMe[meCFCF/.{\[Eta]@_:>1},"ad_nlo/ggq1_cfcf.c"];
    exportMe[meCFCA/.{\[Eta]@_:>1},"ad_nlo/ggq1_cfca.c"];
    exportMe[meCFNF/.{\[Eta]@_:>1},"ad_nlo/ggq1_cfnf.c"];
    (* check I match to Pegasus *)
    Print[FullSimplify[PegasusCFCF-nmeCFCF/.{\[Eta]@_:>1}]];
    Print[FullSimplify[PegasusCFCA-nmeCFCA/.{\[Eta]@_:>1}]];
    Print[FullSimplify[PegasusCFNF-nmeCFNF/.{\[Eta]@_:>1}]];
    (* get some reference numbers *)
    Print@FullSimplify[nmeCFCF/.{\[Eta]@_->1,n->2}];
    Print@FullSimplify[nmeCFCA/.{\[Eta]@_->1,n->2}];
    Print@FullSimplify[nmeCFNF/.{\[Eta]@_->1,n->2}];
    (* compare all 3 expresisons *)
    FullSimplify@Table[With[{},{n,
    PegasusCFCF-MochCFCF,
    Identity[PegasusCFCF-nmeCFCF/.{\[Eta]@_:>1}],
    PegasusCFCA-MochCFCA,
    Identity[PegasusCFCA-nmeCFCA/.{\[Eta]@_:>1}],
    PegasusCFNF-MochCFNF,
    Identity[PegasusCFNF-nmeCFNF/.{\[Eta]@_:>1}]
    }],{n,2,7}]
]

(* gluon-gluon *)
Module[{spmom,rawCACA,rawCANF,rawCFNF,MochCACA,MochCANF,MochCFNF,PGGA,PGGB,PGGC,PegasusCACA,PegasusCANF,PegasusCFNF,meCACA,meCANF,meCFNF,nmeCACA,nmeCANF,nmeCFNF},
    spmom[n_]:= Integrate[x^(n-1)PolyLog[2,x]/(1+x),{x,0,1}];
    (* load raw *)
    rawCACA=getRawS[ggg1,ca^2];
    rawCANF=getRawS[ggg1,ca nf];
    rawCFNF=getRawS[ggg1,cf nf];
    MochCACA=parseRaw@rawCACA;
    MochCANF=parseRaw@rawCANF;
    MochCFNF=parseRaw@rawCFNF;
    (* load Pegasus *)
    PGGA=-(2*NFI+5*NFO+8*NT+7*NS-2*N-2)*8*S1/(NMS*NS*N1S*N2S)-67/9*S1+8/3-4*SSTR2P*(NS+N+1)/(NM*N*N1*N2)+2*S1*SSTR2P-4*SSCHLP+1/2*SSTR3P+(457*NN+2742*NE+6040*NSE+6098*NSI+1567*NFI-2344*NFO-1632*NT+560*NS+1488*N+576)/(18*NMS*NT*N1T*N2T);
    PGGB=(38*NFO+76*NT+94*NS+56*N+12)*(-2)/(9*NM*NS*N1S*N2)+20/9*S1-4/3;
    PGGC=(2*NSI+4*NFI+NFO-10*NT-5*NS-4*N -4)*(-2)/(NM*NT*N1T*N2)-1;
    PegasusCACA=-4*parsePegasus@PGGA;
    PegasusCANF=-1/2*4*parsePegasus@PGGB;
    PegasusCFNF=-1/2*4*parsePegasus@PGGC;
    (* load me *)
    meCACA=transformMe@rawCACA;
    meCANF=transformMe@rawCANF;
    meCFNF=transformMe@rawCFNF;
    nmeCACA=nMe@meCACA;
    nmeCANF=nMe@meCANF;
    nmeCFNF=nMe@meCFNF;
    Print@Collect[meCACA/.{\[Eta]@_:>1},{g3@_,zeta@_,S[R@__,_]},FullSimplify];
    Print@Collect[meCANF/.{\[Eta]@_:>1},{g3@_,zeta@_,S[R@__,_]},FullSimplify];
    Print@Collect[meCFNF/.{\[Eta]@_:>1},{g3@_,zeta@_,S[R@__,_]},FullSimplify];
    exportMe[meCACA/.{\[Eta]@_:>1},"ad_nlo/ggg1_caca.c"];
    exportMe[meCANF/.{\[Eta]@_:>1},"ad_nlo/ggg1_canf.c"];
    exportMe[meCFNF/.{\[Eta]@_:>1},"ad_nlo/ggg1_cfnf.c"];
    (* check I match to Pegasus *)
    Print[FullSimplify[PegasusCACA-nmeCACA/.{SPMOM->g3@n,\[Eta]@_:>1}]];
    Print[FullSimplify[PegasusCANF-nmeCANF/.{SPMOM->g3@n,\[Eta]@_:>1}]];
    Print[FullSimplify[PegasusCFNF-nmeCFNF/.{SPMOM->g3@n,\[Eta]@_:>1}]];
    (* get some reference numbers *)
    Print@FullSimplify[nmeCACA/.{g3->spmom}/.{\[Eta]@_->1,n->2}];
    Print@FullSimplify[nmeCANF/.{\[Eta]@_->1,n->2}];
    Print@FullSimplify[nmeCFNF/.{\[Eta]@_->1,n->2}];
    (* compare all 3 expresisons *)
    FullSimplify@Table[With[{g3n=spmom[n]},{n,
    PegasusCACA-MochCACA/.{SPMOM->g3n},
    Identity[PegasusCACA-nmeCACA/.{SPMOM->g3n,g3@n->g3n,\[Eta]@_:>1}],
    PegasusCANF-MochCANF,
    Identity[PegasusCANF-nmeCANF/.{\[Eta]@_:>1}],
    PegasusCFNF-MochCFNF,
    Identity[PegasusCFNF-nmeCFNF/.{\[Eta]@_:>1}]
    }],{n,2,7}]
]