We have developed a new diagnostic tool for the study of gamma-ray emission lines from radioactive isotopes (such as 26Al and 60Fe) in conjunction with other multi-wavelength observables of Galactic clusters, associations, and alike objects. Our evolutionary synthesis models are based on the code of \cite{CMH94}, which has been updated to include recent stellar evolution tracks, new stellar atmospheres for OB and WR stars, and nucleosynthetic yields from massive stars during hydrostatic burning phases and explosive SN II and SN Ib events. The temporal evolution of 26Al and 60Fe production, the equivalent yield of 26Al per ionising O7 V star (Y26O7 V), and other observables are predicted for a coeval population. The main results are: The emission of the 26Al 1.809 MeV line is characterised by four phases: stellar wind dominated phase (la 3 Myr), SN Ib dominated phase ( ~ 3-7 Myr), SN II dominated phase ( ~ 7-37 Myr), and exponential decay phase (ga 37 Myr). The equivalent yield Y26O7 V is an extremely sensitive age indicator for the stellar population which can be used to discriminate between Wolf-Rayet star and SN II 26Al nucleosynthesis in the association. The ratio of the 60Fe/26Al emissivity is also an age indicator that constrains the contribution of explosive nucleosynthesis to the total 26Al production. We also employed our model to estimate the steady state nucleosynthesis of a population of solar metallicity. In agreement with other works, we predict the following relative contributions to the 26Al production: ~ 9 % from stars before the WR phase, ~ 33 % from WR stars, ~ 14 % from SN Ib, and ~ 44 % from SN II. For 60Fe we estimate that ~ 39 % are produced by SN Ib while ~ 61 % come from SN II. Normalising on the total ionising flux of the Galaxy, we predict total production rates of 1.5 Msun Myr-1 and 0.8 Msun Myr-1 for 26Al and 60Fe, respectively. This corresponds to 1.5 Msun of 26Al and 1.7 Msun of 60Fe in the present interstellar medium. To allow for a fully quantitative analysis of existing and future multi-wavelength observations, we propose a Bayesian approach that allows the inclusion of IMF richness effects and observational uncertainties in the analysis. In particular, a Monte Carlo technique is adopted to estimate probability distributions for all observables of interest. We outline the procedure of exploiting these distributions by applying our model to a fictive massive star association. Applications to existing observations of the Cygnus and Vela regions will be discussed in companion papers.