Problem About Calculation of double Higgs production Cross Section (for gluon fusion) with FeynArts+FormCalc

Question

I want to mention about my problem about FeynArts/FormCalc.I have tried to calculate LHC cross section for "gluon gluon -> H H" process at Standard Model. However, my result XS is quite different from the ones in the literature. I have no idea about this difference or where is the problem. I'm using FeynArts 3.9, FormCalc 9.3 and LoopTools 2.13. I attached the my FeynArts mathematica file.

In literature

XS(gluon gluon -> H H at 14 TeV at LHC)=~0,033 pb

but it is ~0,95 pb in my result.

Could you give me some suggestions? Please, check my mathematica file. Is there any mistakes or missing in my codes.

Thank you very much

<< FeynArts.m

<< FormCalc.m

InitialDirectory;name="glglHH−SM"process=V[5],V[5]−>S[1],>S[1];SetOptions[InsertFields,Model−>"SMQCD",InsertionLevel−>>Particles,ExcludeFieldPoints−>FieldPoint[0][V[5],V[5],>V[5]]];

PaintSE = MkDir[name <> ".diagrams"]; DoPaint[diags_, file_, opt___] := Paint[diags, opt, DisplayFunction -> (Export[ ToFileName[$PaintSE, file <> ".ps"], #] &)]

tops = CreateTopologies[0, 2 -> 2]; ins = InsertFields[tops, process]; DoPaint[ins, "born"]; born = CalcFeynAmp[CreateFeynAmp[ins]]

Print["Counter terms"]

tops = CreateCTTopologies[1, 2 -> 2, ExcludeTopologies -> {TadpoleCTs, WFCorrectionCTs}]; ins = InsertFields[tops, process]; DoPaint[ins, "counter"]; counter = CalcFeynAmp[CreateFeynAmp[ins]]

Print["Self energies"]

tops = CreateTopologies[1, 2 -> 2, SelfEnergiesOnly]; ins = InsertFields[tops, process]; DoPaint[ins, "self"]; self = CalcFeynAmp[CreateFeynAmp[ins]]

Print["Vertices"]

tops = CreateTopologies[1, 2 -> 2, TrianglesOnly]; ins = InsertFields[tops, process, ExcludeParticles -> {F[3 | 4, {1 | 2}]}]; DoPaint[ins, "vert"]; vert = CalcFeynAmp[CreateFeynAmp[ins]];

Paint[ins]

Print["Boxes"]

tops = CreateTopologies[1, 2 -> 2, BoxesOnly]; ins = InsertFields[tops, process, ExcludeParticles -> {F[3 | 4, {1 | 2}]}]; DoPaint[ins, "box"]; box = CalcFeynAmp[CreateFeynAmp[ins]];

Paint[ins]

col = ColourME[All, born]

amps = {born, self, vert, box};

{born, self, vert, box} = Abbreviate[amps, 6, Preprocess -> OnSize[100, Simplify, 500, DenCollect]];

abbr = OptimizeAbbr[Abbr[]];

subexpr = OptimizeAbbr[Subexpr[]];

dir = SetupCodeDir[name <> ".fortran"]

WriteSquaredME[born, {self, vert, box}, col, abbr, subexpr,
LoopSquare -> True, dir]

WriteRenConst[amps, dir]