Dear @linjunze,
Your simulation input (isotope beam, RADDECAY card specifying semi-analogue mode) will indeed take into account the emission of cascade gamma-rays. However, as specified in the manual this happens in an uncorrelated way. This means that the different possible gamma rays are sampled from the yield distribution in the nuclear database of FLUKA and are in no way in direct correlation to the preceding decay in the chain. Since the different isotopes are defined in your input file, I assume that you are checking for each of those specific isotopes how much the cascade gamma rays contribute? You might need to sample many primaries to see the desired effect in a statistically converged way. Is this why you are not certain about your method of simulation?
I’m also happy to point you to two related but more detailed posts on this topic, one about printing info on primary and secondary particles produced in radioactive decay, and one about detailed scoring of the specific gamma cascade energy using an mgdraw routine.
Hope this helps, please let us know in case you need more information,
Andreas