Owing to its flexibility in modeling real-world spatial configurations of users and base stations (BSs), the Poisson cluster process (PCP) has recently emerged as an appealing way to model and analyze heterogeneous cellular networks (HetNets). Despite its undisputed relevance to HetNets -- corroborated by the models used in industry -- the PCP's use in performance analysis has been limited. This is primarily because of the lack of analytical tools to characterize performance metrics such as the coverage probability of a user connected to the strongest BS. In this paper, we develop an analytical framework for the evaluation of the coverage probability, or equivalently the complementary cumulative density function (CCDF) of signal-to-interference-and-noise-ratio (SINR), of a typical user in a K-tier HetNet under a max power-based association strategy, where the BS locations of each tier follow either a Poisson point process (PPP) or a PCP. The key enabling step involves conditioning on the parent PPPs of all the PCPs which allows us to express the coverage probability as a product of sum-product and probability generating functionals (PGFLs) of the parent PPPs. In addition to several useful insights, our analysis provides a rigorous way to study the impact of the cluster size on the SINR distribution, which was not possible using existing PPP-based models.