  We consider power corrections due to a finite top quark mass M_t to the
production of a Higgs boson pair within the Standard Model at next-to-leading
order (NLO) in QCD. Previous calculations for this process and at this
precision were done in the limit of an inifinitely heavy top quark. Our results
for the inclusive production cross section at NLO include terms up to
O(1/M_t^12). We present the Mathematica package TopoID which for arbitrary
processes aims to perform the necessary steps from Feynman diagrams to
unrenormalized results expressed in terms of master integrals. We employ it for
advancing in this process towards next-to-next-to-leading order (NNLO) where
further automatization is needed.
