I despair of the definition of primary particles, which Allpix^2 seems to use.
My simulation consists of several successive layers of pixel sensors (ibl-planar-model) and a proton source (200MeV).
Since I want to display the track of the particles, I read the hits from MCPixelHit and get the corresponding particles with hit->getMCParticles() . Afterwards I get the type of the parent particle with particle->getTrack()->getParent()->getParticleID(). Additionally, I get the volume in which the particle was created with particle->getTrack()->getOriginatingVolumeName(). I have inserted the whole function below:
void test() {
// create ofstream and open outputfile
ofstream outputFile("output/test.csv");
// write header of the file
outputFile << "Originating_Volume" << "," << "ParentID"<<std::endl;
//-------------------------------------------------------------------
//open File with TTree
TFile *file = new TFile("output/test.root");
// Load PixelHit
TTree* PixelHit_tree = static_cast<TTree*>(file->Get("PixelHit"));
if(!PixelHit_tree) {
std::cout << "Could not read tree PixelHit, cannot continue." << std::endl;
}
// Load MCParticle
TTree* MCParticle_tree = static_cast<TTree*>(file->Get("MCParticle"));
if(!MCParticle_tree) {
std::cout << "Could not read tree MCParticle" << std::endl;
}
// Load MCTrack
TTree* MCTrack_tree = static_cast<TTree*>(file->Get("MCTrack"));
if(!MCTrack_tree) {
std::cout << "Could not read tree MCTrack" << std::endl;
}
//------------------------------------------------------------------------
// number of detectors
int detector_number = PixelHit_tree->GetListOfBranches()->GetEntries();
// number of events
int event_number = PixelHit_tree->GetEntries();
//-----------------------------------------------------------------------
//Loop over all detectors
for(int det=0; det!=detector_number ;det++){
const char* detector =PixelHit_tree->GetListOfBranches()->At(det)->GetName();
//------------------------------------------
// Get branches
// PixelHit
TBranch* PixelHit_branch = PixelHit_tree->FindBranch(detector);
//MCParticle
TBranch* MCParticle_branch = MCParticle_tree->FindBranch(detector);
//MCTrack
TBranch* MCTrack_branch = MCTrack_tree->FindBranch("global");
//-------------------------------------------
//PixelHit
std::vector<allpix::PixelHit*> Hits;
PixelHit_branch->SetObject(&Hits);
//MCParticle
std::vector<allpix::MCParticle*> Particles;
MCParticle_branch->SetObject(&Particles);
//MCTrack
std::vector<allpix::MCTrack*> Tracks;
MCTrack_branch->SetObject(&Tracks);
//-----------------------------------
//Loop over all events
for(int event=0; event<event_number; event++){
// Get Entry of event
PixelHit_tree->GetEntry(event);
MCParticle_tree->GetEntry(event);
MCTrack_tree->GetEntry(event);
//-----------------------------------------
//Loop over all hits
for(auto& hit : Hits){
// get vector of related particles
std::vector<const allpix::MCParticle*> hit_particles = hit->getMCParticles();
//------------------------------------
//Loop over all related particles
for(auto& particle : hit_particles){
std::string volume = particle->getTrack()->getOriginatingVolumeName();
int parent = 0;
if (particle->getTrack()->getParent() !=nullptr){
parent = particle->getTrack()->getParent()->getParticleID();
}
// write this into the outputfile
outputFile << volume << "," << parent<< std::endl;
}
}
}
PixelHit_tree->ResetBranchAddresses();
MCParticle_tree->ResetBranchAddresses();
MCTrack_tree->ResetBranchAddresses();
}
// close the output file
outputFile.close();
}
Particles which were started as primary particles in the simulation are created in the “world” volume and have no parent particle. Particles that are not from the “world” volume are in my opinion secondary particles and should be assigned to a parent particle. This is also the case in most cases. However, there are a small number of particles that are not produced in the “world” volume and do not have a parent particle. Do I have a mistake? What can be the reason for this or how can these particles be interpreted physically?
sorry for your desparation, I hope I can shed some light onto this issue. In APSQ the definition of “particle” deviates slightly from the one in Geant4, simply because we are usually interested in what happens inside the sensor.
Hence, we define MCParticles per sensor and call all those primary that enter the sensor from the outside. All particles generated by interactions within the sensor, like delta rays, are consequently called secondary. This I would argue is a quite common nomenclature.
I guess this becomes a bit more complicated as particles propagate further. Any particle not absorbed in the first sensor also penetrate the second detector. There they are counted as primary particle because they originate outside this sensor’s volume - but their originating volume is not the world volume because they were created in the first sensor.
This approach allows to unanimously distinguish for each sensor whether a certain particle was created inside or outside its volume - which can be quite handy. Now in order to be able to track particles between sensors, we have the MCTrack object which follows the physical particle throughout its lifetime.
Tl;dr: depending on what your goal of filtering is, you either have to look at MCParticle (per sensor) or MCTrack (for the full setup). Does this help?
Thank you very much for your answer. Unfortunately, it doesn’t help me quite yet. I want to track the particles through the whole setup, so I use the getParent() function of MCTrack. This should really mean the parent particle, right?
However, as described above, there are a few particles that do not have a “parent track” and still do not come from the world volume but from the detector volume. For me this point is not understandable.
I would take MCParticles from each detector, filter them for primaries (empty parentID) and then check their assigned MCTrack. If the MCTrack for two MCParticles in different detectors is the same, they have been created by the same physical particle - and this is the information you are after, right?
for MCTracks I do not really understand where objects without parent should come from, if not from the world. If you could provide me a simplified version of your configuration files which reproduce this issue, I can have a look.
Also, if you have suggestions how we could make this whole MC truth information more accessible, let me know! Maybe adding some more detailed description to the user manual would already help?
[detector1]
type = "ibl_planar"
position = 0 0 0
orientation = -180deg 0 0
[detector2]
type = "ibl_planar"
position = 0 0 1mm
orientation = -180deg 0 0
And with this example script I find particles without parent track:
//open File with TTree
TFile *file = new TFile("output/test.root");
// Load MCParticle
TTree* MCParticle_tree = static_cast<TTree*>(file->Get("MCParticle"));
// Load MCTrack
TTree* MCTrack_tree = static_cast<TTree*>(file->Get("MCTrack"));
//------------------------------------------------------------------------
// number of detectors
int detector_number = MCParticle_tree->GetListOfBranches()->GetEntries();
// number of events
int event_number = MCParticle_tree->GetEntries();
//-----------------------------------------------------------------------
//Loop over all detectors
for(int det=0; det!=detector_number ;det++){
const char* detector =MCParticle_tree->GetListOfBranches()->At(det)->GetName();
//------------------------------------------
// Get branches
//MCParticle
TBranch* MCParticle_branch = MCParticle_tree->FindBranch(detector);
//MCTrack
TBranch* MCTrack_branch = MCTrack_tree->FindBranch("global");
//-------------------------------------------
//MCParticle
std::vector<allpix::MCParticle*> Particles;
MCParticle_branch->SetObject(&Particles);
//MCTrack
std::vector<allpix::MCTrack*> Tracks;
MCTrack_branch->SetObject(&Tracks);
//-----------------------------------
//Loop over all events
for(int event=0; event<event_number; event++){
// Get Entry of event
MCParticle_tree->GetEntry(event);
MCTrack_tree->GetEntry(event);
//------------------------------------
//Loop over all particles
for(auto& particle : Particles){
std::string volume = particle->getTrack()->getOriginatingVolumeName();
int parent = 0;
if (particle->getTrack()->getParent() !=nullptr){ //check if particle has parent track
parent = particle->getTrack()->getParent()->GetUniqueID();
}
else if (volume != "World"){ // if particle has no parent track, check if it originates from world
std::cout << particle->getParticleID() << ", " << volume << ", " << particle->getTrack()->getCreationProcessName()<< std::endl;
}
}
}
MCParticle_tree->ResetBranchAddresses();
MCTrack_tree->ResetBranchAddresses();
}
I think the manual has a small paragraph in the FAQ about how to find primaries. In my opinion, it is definitely a good idea to go into this in more detail, maybe similar to what you described above.
@pschutze and me took a deep dive into the code and believe we found the issue you are facing. The way we currently follow track creation and MCParticles relies on them interacting in the sensitive volume of your detector in order to be registered. This greatly simplifies the way we can treat event information and reduces the amount of data we need to temporarily cache and sift through. The problem you discovered lies directly here:
If particle A interacts with detector 1, then creates a secondary, particle B, in some support material (e.g. ASIC), which in turn creates particle C, and we only see particle C in detector 2 we have a problem: we will only see particle A and particle C within the sensors but particle B is lost. Therefore also the relation C → B → A is lost and C suddenly has no parent despite not originating from the world volume.
To some extend, this issue is addressed in the manual:
The MCParticle additionally stores a parent MCParticle object, if available. The lack of a parent doesn’t guarantee that this MCParticle originates from a primary particle, but only means that no parent on the given detector exists. Also, the MCParticle stores a reference to the MCTrack it is associated with.
We agree that from a physical perspective this would definitely make sense and we’ll continue to look into a solution but can’t promise anything right now…