Hi,
For my purpose i need to retrieve the true interactions positions of gamma rays interacting with a silicon detector as well as the deposited energies. For this reason i was playing around with MCTrack and MCParticle. The simulation consists in 2 point sources, 9 cm away from the detector surface.
Although i though MCTrack was enough for this purpose i thought of having a look anyway at MCParticle but i found some discrepancies and unexpected behaviors and i am not sure if those arise from errors in the allpix functions or from my understanding of them.
Using a python script i was getting data from MCTrack and MCParticle as
for mc_track in br_mc_track:
#here we record all the interactions the event had with the detector togheter with the tracks generated
data_array[1] = mc_track.getStartPoint().x()
data_array[2] = mc_track.getStartPoint().y()
data_array[3] = mc_track.getStartPoint().z()
data_array[4] = mc_track.getGlobalStartTime()
data_array[5] = mc_track.getKineticEnergyInitial()
data_array[6] = mc_track.getCreationProcessName()
data_array[7] = mc_track.getParticleID()
data_array[8] = hex(ROOT.addressof(mc_track))
data_array[9] = hex(ROOT.addressof(mc_track.getParent()))
data_array[10] = mc_track.getOriginatingVolumeName()
# Create a dictionary from the array and column names
new_row = dict(zip(df.columns, data_array))
print(iev)
# Convert the dictionary to a DataFrame and append it
df = pd.concat([df, pd.DataFrame([new_row])], ignore_index=True)
for mc_particle in br_mc_particle:
data_array[1] = mc_particle.getGlobalStartPoint().x()
data_array[2] = mc_particle.getGlobalStartPoint().y()
data_array[3] = mc_particle.getGlobalStartPoint().z()
data_array[4] = mc_particle.getGlobalTime()
data_array[5] = mc_particle.getKineticEnergyStart()
data_array[6] = f"mc_particle, is it primary? {mc_particle.isPrimary()}. Connected track: {hex(ROOT.addressof(mc_particle.getTrack()))}"
data_array[7] = mc_particle.getParticleID()
data_array[8] = hex(ROOT.addressof(mc_particle))
data_array[9] = hex(ROOT.addressof(mc_particle.getParent()))
data_array[10] = "parent particle: {hex(ROOT.addressof(mc_particle.getPrimary()))}"
new_row = dict(zip(df.columns, data_array))
print(iev)
df = pd.concat([df, pd.DataFrame([new_row])], ignore_index=True)
respectively where iev is the event ID.
An example of results is shown below
We can see that for event IDs 104,113,123 everything looks fine, the first row(s) in each event are from MCTrack, and the initial kinetic energy (init_kin) is 511 keV as set in the .conf file for the source. the same is shown in the MCParticle row(s), with the difference that the initial coordinates are in proximity with the detector surface (which is at z=0) since MCParticle start saving informations of the particle once this enters the detector.
But then when we start to have interactions in the active area of the detector as for event 124, in the MCParticle rows (distinguishable from having mc_particle under the column process_name), we see that the starting kinetic energy of the primary particle is not 511 kev as expected but lower like 482keV.
If we look at the deposited energy from the first compton interaction by summing the initial kinetic energy of the secondary particles related to the interaction we notice that 482keV should rather be the energy of the primary particle after the interaction (so after having deposited some of its energy since 511-28.36=482.64keV).
Similarly happens for all other cases in which there is an interaction in the detector, where the init_kin from MCParticle is lower than 511kev and very similar to the final energy of the primary particle after having interacted.
On a last note if i calculate the deposited energy in the compton interaction looking at MCTrack for event 124, i would sum 28.5keV+0.1keV=28.6keV, while from MCParticle= 28.36keV+0.0keV =28.36keV. Similarly for all other interaction cases the results seem to be slightly different. Could this be because MCParticle looks at the true deposited energy while MCTrack has the true deposited energy + other effects present in the physics model such as doppler broadening? (i am using FTFP_BERT_PEN)
I also looked at getTotalDepositedCharge() for MCParticle and from that divide by 2 and multiply by 3.64eV and i do get similar results as from the other two ways mentioned before to measure the deposited energy although not exactly the same and i am not sure what is the more accurate/best way among those 3.
Bests
Luca
