Charge creation type conversion introduces peak shift to lower energies

I am doing simulations with X-rays in the range of 1~30keV. I noticed that there is a consistent peak shift to lower energies, even at full absorption of the gamma ray. I found out that this is a result of how the charge creation is performed. After sampling a Gaussian based on the Fano factor, the result is truncated to form the number of charge carriers. This means that if 10.99 get sampled from the Gaussian distribution, 10 charges will be created. When this happens multiple times over the course of a single event, this effect accumulates and introduces a bias towards lower energies.

Changing this line to instead round down to the nearest integer, i.e. changing line:

To:

auto charge = static_cast<unsigned int>(std::max(charge_fluctuation(random_generator_), 0.) + 0.5);

Will fix this behaviour. The two are compared in the figures below, with the version including the change above labeled as “Improved”.

A comparison of the deposition spectra comparing the result with the modification shown above is shown below. Here we have an incident gamma of 10keV, default charge creation energy (3.64eV) and Fano factor of 0 is shown below. We would expect a mean close to 10,000/3.64 = 2747.25 charges at full absorption. As you can see, there is a bias towards lower energies, as a result of the truncation of the number of charges after the gaussian sampling. This will result in a peak shift of the main energy peak.
num charges

As a test to verify this, a range of simulations were done varying the energy of the incident gamma ray. The gamma rays were shot at a cube of silicon 10x10x10cm. The total number of events for each simulation of events was 300k. The energy depositions were extracted, the main peaks fitted with a Gaussian, and the peak shift with respect to the input energy was measured. They can be reproduced using the configurations below:
model.conf (285 Bytes)
detector.conf (74 Bytes)
setup.conf (577 Bytes)
Below you can see that the peak shifts more and more as we increase the energy. This is caused by more energy depositions happening in the active volume, thus increasing accumulation of truncation errors. By instead rounding the number of charges to the nearest integer, the main peak will instead average to the nearest whole number after division. For example a 10keV incident gamma ray, 10,000/3.64 = 2747.25 will average to 2747.

Hi @sblokhuizen

wow, thank you very much for reporting this - this is a nice investigation and very good conclusion. :slight_smile:

We’d like to fix this int he next release - in order to give you the credit due, would you like to make the change yourself so your contribution is recorded in git? If yes, please let us know - we gladly help you through the steps necessary!

Best regards,
Simon

Hi Simon,

Thank you! I tried to create a merge request through git before writing this post, but I was having trouble with the CERN authentication. I now created a pull request on github instead following the CONTRIBUTING.md document. This is the first time I have done this, so please let me know if I did anything wrong :slight_smile:

1 Like