David,
The answer given by jrellin to a question on a recent post does a good job of describing the phenomena which contribute to the variance of a full energy peak Simulated count variance much smaller than experimental — how to model realistic detector fluctuations? . I would say it is a worthwhile read for you. If you are trying to estimate a detector response for a detector you already have, then it would be much easier to simulate the energy deposition of a gamma ray in your system, (which you would statistically characterize) and then statistically broaden the energy deposition tally.
Back to the current question, I think the reason you are still observing a tail is because the optical photon emission spectrum you have defined. However, I would like to reiterate, a bit more strongly, that you should not use the energy deposited in your photomultiplier (PM) as a metric to estimate your full energy peak. While optical photon wavelength impacts the probability of reaching and interacting with and producing a photoelectron in a PM, it does not broaden the full energy peak the way that is shown in your graphs.
As a though experiment, let us assume that both a 3 and 6 eV photon have a 100% probability to reach the PM, likewise the PM sensitivity at these energies is also 100%. Furthermore we have magically altered our scintillator to only yield 3 and 6 eV photons. The question to ask then is, would energy deposition in the PM by photons be representative of our full energy peak. The answer would be no, this is true even if the variance introduced by the PM and readout electronics is negligible. This is because both a 3 and 6 eV photon would produce a single photoelectron in the PM. In short by using the energy deposition in the PM from optical photons as your metric, you are incorrectly introducing a broadening term.
If you set the scintillation parameters so that you only emit monoenergetic optical photons, then I think that high energy tail would go away.