In fusion energy systems (FES), high energy neutrons are emitted from the plasma source - due to the D-T fusion reaction - enabling them to penetrate deep in the materials surrounding the core. Energy is then deposited along the path of the neutrons due to interactions with nuclides, resulting in - besides nuclear heating - two main effects; radiation damage and transmutation. Radiation damage causes changes in the macroscopic properties of the materials due to microscopic changes that result from interactions of high energy neutrons with nuclides. Transmutation is caused by the absorption of neutrons by nuclides in the medium and almost always results in the production of radioactive nuclides. Such radioactive nuclides are of importance to FES design and operation as they persist after the shutdown of the facility due to their long half lives. Efforts are directed to quantify the shutdown dose rate (SDR) that results from gamma emitting nuclides produced by transmutation. Monte Carlo (MC) methods are favored over deterministic methods for the simulation of particles transport in FES due to complexity of the models and to reduce the uncertainties/errors of the predicted particle flux distributions due to approximations. The rigorous 2-step method (R2S) relies on dedicated activation calculations to predict the photon emission density distribution, and is widely used for SDR quantification. It involves a neutron transport step, activation analysis to obtain the photon emission density, and a photon transport step to calculate the SDR. It is often the case that neutrons suffer attenuation in traversing the medium from the plasma source - due to interactions with nuclides - and that results in a steep gradient in the neutron flux. Variance reduction (VR) tools have been developed with the primary goal of pushing neutrons - simulated particles - to regions of the phase-space that are of importance for the quantities under consideration in order to reduce the uncertainty in the MC results. The recently developed Group-wise Transmutation - Consistent Adjoint Driven Importance Sampling (GT-CADIS) method provides a capability to obtain the photon emission density distribution as a function of the energy dependent group-wise neutron flux distribution via linearization of the transmutation operator. Using the photon emission density it is possible to overcome previous difficulties of the error propagation in the R2S workflow. One primary concern with the R2S workflow is that only the contribution of the photon transport step is considered as a measure of the uncertainty of the calculated SDR, while the contribution from the neutron transport step remains undefined. Previous methods have tried to tackle this issue but there was always difficulty in obtaining the correlation of the neutron fluxes and that resulted in implementing either impractical approximations or just calculating the upper and lower bounds of the uncertainty of the SDR. In this document, the R2S workflow has been investigated. First, issues related to the neutron transport step and the uncertainty of the photon emission density have been addressed. Second, a scheme was developed to propagate the statistical uncertainty of the neutron transport step to the SDR. Starting with the neutron transport step, a variation of the main R2S that aimed at increasing the resolution while reducing the computational expenses was found to introduce systematic errors that might undermine the gain in the computational cost. One of the difficulties in propagating the neutron flux uncertainty to the photon emission density is obtaining the correlation values between the neutron fluxes in different energy groups and mesh voxels. By utilizing the GT method, an approximation to the calculation of the correlation coefficients has been investigated building on the fact that using group-wise transmutation the correlation terms needed were greatly reduced. It was discovered that the correlation between the neutron fluxes in different energy groups is a function of the material composition. That facilitated obtaining the needed correlation matrix and quantifying the uncertainty of the photon emission density. A method to propagate the photon source uncertainty to the SDR by random sampling was developed and was demonstrated to be efficient on various types of numerical experiments as well as a production level problem.