The current concordance lCDM cosmological model describes a universe where cold dark matter seeds structure formation and a cosmological constant drives its accelerated expansion. Precise measurements of various astronomical observables allow us to test this model and any deviations, if found, may lead to an improved cosmological theory. Ongoing and planned large scale surveys of the skies have the power to study the lCDM model. However the data sets they generate will be dominated by complex systematic uncertainties. One probe of cosmological parameters, the evolution of clusters of galaxies, has the power to differentiate simple models of dark energy, like the cosmological constant, from possible modifications of the theory of gravity since in can be used to trace the expansion as well as growth history of the universe. The formation of galaxy clusters at early times contains information about their gravitational environment and the late time halting of cluster mergers probes the accelerated expansion of the universe. However galaxy clusters are very complex systems. The underlying dark matter halos hosting them cannot be described analytically and require numerical simulations to understand. In addition the infalling baryonic matter is heated and capable of radiating energy leading to star formation and various feedback systems. Untangling these astrophysical complications can be accomplished by observing galaxy clusters in multiple wavelengths that reveal their different aspects as well as studying them via numerical simulations. These different observables must be combined appropriately, accounting for their different selection functions, correlations as well as observational effects. I will discuss one such study performed using optical, X-Ray and Sunyaev-Zel’dovich (SZ) simulations in order to explain a discrepancy seen by the Planck satellite when stacking its SZ observations. Careful accounting of systematic effects reduces this discrepancy beyond statistical significance and demonstrates the importance of careful modeling of the various selection functions as well as accounting for the uncertainties in all of the observables.