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.