Simulation of neutral Higgs pairs production processes in PYTHIA using HPAIR matrix elements

Double Higgs production is currently not fully implemented in PYTHIA, as only the resonant $H\rightarrow hh$ MSSM processes and $f_i\bar{f}_i$ production modes can be generated. Using matrix elements from the HPAIR program we have added in PYTHIA the gluon fusion production modes in Standard Model and MSSM. Those matrix elements include resonance process as well as continuum production, which is relevant for SM and MSSM high $\tan\beta$ scenarios. Comparison between PYTHIA Monte Carlo results and VEGAS integrated cross sections are shown.