Dear colleagues,
I’m simulation the response of the Timepix3 detector to gammas. So, I have a very simple setup: beam source (1mm size) and detector (from default example) placed in 2 cm from the source. Then I would like to calculate the efficiency of gammas registration for different energies, separately for events produced by photoelectric absorption and by Compton scattering.
To do this for each pixel with charge above threshold I call MCparticles associated with this pixel in the following way:
for (auto& pixel_charge : input_charges) { /// Over all pixel charges in event
int x = pixel_charge->getIndex().x();
int y = pixel_charge->getIndex().y();
int charge = pixel_charge->getCharge();
for (auto& mc_part : pixel_charge->getMCParticles()) {
auto track = mc_part->getTrack();
auto particleId = track->getParticleID();
auto processName = track->getCreationProcessName();
auto processId = track->getCreationProcessType();
}
}
But, sometimes with one pixel several MC particles are associated:
ev: 51
x y particleId processName processId
130 132 22 none -1
130 132 22 phot 2
130 132 11 phot 2
130 132 11 phot 2
It’s clear how it’s possible from physics. But, should I select only the last entry (particleId = photoelectron 11 + processName = phot ) for the proper assignment “pixel – interaction process”?
For Compton events the situation is a little bit more complicated:
ev: 9892
127 128 22 none -1
127 128 11 compt 2
131 122 22 none -1
131 122 11 phot 2
131 122 11 phot 2
131 123 11 phot 2
131 123 11 phot 2
Again, here for pixel (127,128) should I select only the last entry to determine that pixel was fired by Compton scattering?
Is it right that entries in MCparticles object are written in chronological order of the interaction process?
Hope my question is more or less clear.
Thank you in advance,
Petr
1 Like
Dear @psmolyan
I am not 100% sure I fully understand your goal, but I nevertheless have a few comments or maybe directions you might want to look into:
-
There is no guarantee that interactions or MCParticles are ordered by time, but since this is usually the order in which Geant4 produces them and the order is preserved, that is correct.
-
If I understand correctly, you are interested in the first interaction the photon from your source undergoes in the detector, irrespectively of what happens afterwards (i.e. if you have a compton scattering and subsequent absorption you still want the compton scattering). In that case you could restrict your MCParticles to those being primary particles by checking:
if(particle->getParent() == nullptr) {
// this is a primary particle, meaning it came "from the outside" of the sensor
}
-
You refer to thresholds in your text but are operating on PixelCharge objects. Have you considered adding a simple digitizer to the end of your simulation chain that would create PixelHit objects? In that case you would also directl have the method getPrimaryMCParticles() which would only give you the above mentioned subset of particles.
-
Whether a pixel fired “because of the compton scattering” is a more complicated question, because
- The Energy resulting from the Compton scattering might be shared between multiple pixels
- The energy in a single pixel might not be sufficient for crossing the threshold, but may contribute significantly
Here of course things can easily become complicated, so it depends on how precise your answer has to be in order ti judge how deep you would like to dig.
I hope this helped you a bit - please add some more explanation if I misunderstood your question or my interpretation is wrong or the answer isn’t sufficient.
All the best,
Simon
Dear Simon,
thank you very much for your quick answer!
Sure, I apply the threshold after getting the charge.
My main goal is to know how many pixels (actually clusters) were activated by photoelectric absorption and by Compton scattering.
If I understand correctly, you are interested in the first interaction the photon from your source undergoes in the detector, irrespectively of what happens afterwards (i.e. if you have a compton scattering and subsequent absorption you still want the compton scattering).
Exactly.
Initially, I was confused by the fact that one pixel can be associated with several MC particles with different particle ID and the same process name (like for event 51 above). But, typically the last entry for such MC particle vector was particleId = 11 with processName = phot. So, I decided to assign this last entry with pixel.
Now I found good example for Compton event:
ev: 979
195 183 11 phot 2 202.669
195 184 11 phot 2 202.669
196 183 11 phot 2 202.669
196 183 11 compt 2 70.4692
196 184 22 none -1 300
196 184 11 phot 2 202.669
196 184 11 compt 2 70.4692
196 185 11 phot 2 202.669
196 185 11 compt 2 70.4692
197 183 22 phot 2 23.1526
197 183 11 phot 2 18.2384
197 183 11 phot 2 202.669
197 183 11 compt 2 70.4692
197 184 22 none -1 300
197 184 22 phot 2 23.1526
197 184 11 phot 2 18.2384
197 184 11 phot 2 202.669
197 184 11 compt 2 70.4692
198 183 11 phot 2 18.2384
198 183 11 compt 2 70.4692
198 184 22 phot 2 23.1526
198 184 11 phot 2 18.2384
198 184 11 compt 2 70.4692
Now the last column is deposited_energy_mc_keV = (track->getKineticEnergyInitial() - track->getKineticEnergyFinal() )*1000
One can see that e.g. pixel (197;183) is connected with photoelectric absorption and Compton scattering. Now my approach with the last entry gives an ambiguous result.
Best regards,
Petr
Hi @psmolyan
let’s have a closer look at the even you are showing - could you also print the IDs of the MCParticles as well as their parents?
std::cout << "ID: " << mc_part << " Parent: " << mc_part->getParent();
This should give is the memory addresses which we can use as IDs. We could then see the inheritance tree - some of the particles you’re seeing might be results from the primary interaction.
Best,
Simon
Hi @psmolyan
so, interestingly your primary gamma is 0x556e4bde7240 and all others are created as daughter particles of it. However, the only two track segments appearing for this are the ones stating none as their interaction.
If we look at
196 183 11 compt 2 70.4692 0x556e4c04e350 0x556e4bde7240
196 184 11 compt 2 70.4692 0x556e4c04e350 0x556e4bde7240
we see that one interaction of a single photon via Compton effect contributed to the total charge in two adjacent pixels.
The issue with the energy difference you are adding is that this is the energy initially left in the sensor by the interaction but does not tell you how much was distributed to each of the pixels. If you are really interested in looking at the processes in this detail, one option would be to also keep the intermediate PropagatedCharge objects. For each of your PixelCharge you could then do:
for(auto prop_charge : pixel_charge->getPropagatedCharges()) {
// let's see how much charge this sub-contirbution has:
auto charge = prop_charge->getCharge();
// let's see from which process it came:
auto process_id = prop_charge->getMCParticle()->getTrack->getCreationProcessType();
// now we can build a detailed table for each pixel of how much energy came from which process
}
Does that make sense?
Simon
Hi @simonspa
Great! I do understand. Thank you very much for this suggestion.
Unfortunately, for this run I excluded saving of the Propagated charge to save the disk space. I will rerun it and will try to dig into this issue.
Best regards,
Petr
Now I guess my issue is solved. I summed propagated charge for each pixel for two processes separately:
if (event_num == 980)
{
for (auto prop_charge : pixel_charge->getPropagatedCharges()) {
auto charge = prop_charge->getCharge();
auto process_id = prop_charge->getMCParticle()->getTrack()->getCreationProcessName();
if (process_id == “compt”) sumcharge_compt += charge;
if (process_id == “phot”) sumcharge_phot += charge;
}
cout << “prop charge: " << x << " " << y << " compt: " << sumcharge_compt * Epair << " phot: " << sumcharge_phot * Epair<<”\n" << endl;
}
And final output is:
mc_part: 94 39 22 phot 2 23.1526 0x557fe2b8aad0 0x557fe2a02220
mc_part: 94 39 11 phot 2 19.1585 0x557fe2a011e0 0x557fe2b8aad0
prop charge: 94 39 compt: 0 phot: 19.2218
mc_part: 97 44 11 compt 2 55.3611 0x557fe2ac3850 0x557fe2a02220
prop charge: 97 44 compt: 22.9607 phot: 0
mc_part: 97 45 22 none -1 300 0x557fe2a02220 0
mc_part: 97 45 11 compt 2 55.3611 0x557fe2ac3850 0x557fe2a02220
prop charge: 97 45 compt: 28.2412 phot: 0
mc_part: 99 41 22 none -1 300 0x557fe2a02220 0
mc_part: 99 41 11 phot 2 68.7847 0x557fe2a012a0 0x557fe2a02220
mc_part: 99 41 11 compt 2 79.8813 0x557fe2a01420 0x557fe2a02220
prop charge: 99 41 compt: 0.03101 phot: 66.1355
mc_part: 100 40 22 none -1 300 0x557fe2a02220 0
mc_part: 100 40 11 phot 2 68.7847 0x557fe2a012a0 0x557fe2a02220
mc_part: 100 40 11 compt 2 79.8813 0x557fe2a01420 0x557fe2a02220
prop charge: 100 40 compt: 38.7448 phot: 0.11075
mc_part: 100 41 11 phot 2 68.7847 0x557fe2a012a0 0x557fe2a02220
mc_part: 100 41 11 compt 2 79.8813 0x557fe2a01420 0x557fe2a02220
prop charge: 100 41 compt: 38.1201 phot: 0.77525
mc_part: 105 33 22 none -1 300 0x557fe2a02220 0
mc_part: 105 33 11 compt 2 64.8215 0x557fe2ac39d0 0x557fe2a02220
prop charge: 105 33 compt: 36.3526 phot: 0
mc_part: 105 34 11 compt 2 64.8215 0x557fe2ac39d0 0x557fe2a02220
prop charge: 105 34 compt: 12.9622 phot: 0
mc_part: 106 33 11 compt 2 64.8215 0x557fe2ac39d0 0x557fe2a02220
prop charge: 106 33 compt: 6.09125 phot: 0
Now I can simply select the highest summed propagated charge per pixel and assign corresponding process to this pixel.
Thank you, Simon, for your help!
Best regards,
Petr
You are most welcome, I am happy we could solve this! 