examples/track_matching/reco_particles_track_matching.cc

0001 // Very minimal and verbose example of how to use association IDs for MCParticles and ReconstructedParticles
0002 // In general there are two ways associations may work: through PODIO references and IDs
0003 // This example illustrates how to use IDs
0004 //
0005 // Created by Dmitry Romanov
0006 // Subject to the terms in the LICENSE file found in the top-level directory.
0007 //
0008
0009 #include <TFile.h>
0010 #include <TTree.h>
0011 #include <TTreeReader.h>
0012 #include <TTreeReaderArray.h>
0013 #include <fmt/core.h>
0014 #include <string>
0015
0016 void reco_particles_track_matching(const std::string &file_name) {
0017     auto *file = new TFile(file_name.c_str());
0018     auto *tree = (TTree *) file->Get("events");
0019     TTreeReader tree_reader(tree);       // !the tree reader
0020     tree->Print();
0021
0022     // Reconstructed particles pz array for each reconstructed particle
0023     TTreeReaderArray<float> reco_pz_array = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
0024
0025     // MC particle pz array for each MC particle
0026     TTreeReaderArray<float> mc_pz_array = {tree_reader, "MCParticles.momentum.z"};
0027
0028     // Next arrays correspond to particle associations
0029     // Each association has 2 ids - indexes in corresponding reco and MC arrays
0030     TTreeReaderArray<unsigned int> rec_id = {tree_reader, "ReconstructedChargedParticleAssociations.recID"};
0031     TTreeReaderArray<unsigned int> sim_id = {tree_reader, "ReconstructedChargedParticleAssociations.simID"};
0032
0033
0034     // Read 100 events
0035     tree_reader.SetEntriesRange(0, 100);
0036     while (tree_reader.Next()) {
0037
0038         // Number of mc particles, reco particles and associations may differ
0039         fmt::print("New event. N reco particles: {}, N mc particles: {}, N assoc: {}\n",
0040                    reco_pz_array.GetSize(), mc_pz_array.GetSize(), rec_id.GetSize());
0041
0042         // Iterate over associations
0043         for(unsigned int i=0; i < rec_id.GetSize(); i++) {
0044
0045             // For each association pull index of reco and MC array
0046             auto reco_array_index = rec_id[i];
0047             auto mc_array_index = sim_id[i];
0048
0049             float reco_pz = reco_pz_array[reco_array_index];
0050             float mc_pz = mc_pz_array[mc_array_index];
0051             fmt::print("   reco={:>10.4f} mc={:>10.4f}\n", reco_pz, mc_pz);
0052         }
0053     }
0054 }
0055
0056
0057 int main() {
0058     reco_particles_track_matching("input.edm4eic.root");
0059     return 0;
0060 }