evaluation/source/evaluation.cc

0001
0002 #include <stdio.h>
0003 #include <stdlib.h>
0004
0005 // ROOT
0006 #include "TTree.h"
0007 #include "TFile.h"
0008
0009 //#define _USE_RECONSTRUCTED_TRACKS_
0010
0011 // NPdet
0012 #include "dd4pod/Geant4ParticleCollection.h"
0013 #include "edm4eic/ReconstructedParticleCollection.h"
0014 #include "edm4eic/CherenkovParticleIDCollection.h"
0015 #include "edm4eic/CherenkovPdgHypothesis.h"
0016
0017 int main(int argc, char** argv)
0018 {
0019   // Command line "parser";
0020   if (argc != 2) {
0021     printf("usage: %s <root-data-file>\n", argv[0]);
0022     exit(0);
0023   } //if
0024
0025   // .root file with event tree;
0026   auto fdata = new TFile(argv[1]);
0027   if (!fdata) {
0028     printf("input file '%s' does not exist\n", argv[1]);
0029     exit(0);
0030   } //if
0031   TTree *t = dynamic_cast<TTree*>(fdata->Get("events"));
0032   if (!t) {
0033     printf("input file '%s' does not have \"events\" tree\n", argv[1]);
0034     exit(0);
0035   } //if
0036
0037   // Use MC truth particles for a "main" loop;
0038   auto mctracks   = new std::vector<dd4pod::Geant4ParticleData>();
0039   auto rctracks   = new std::vector<edm4eic::ReconstructedParticleData>();
0040   auto cherenkov  = new std::vector<edm4eic::CherenkovParticleIDData>();
0041   t->SetBranchAddress("mcparticles", &mctracks);
0042
0043   // FIXME: or whatever the branches are called;
0044 #ifdef _USE_RECONSTRUCTED_TRACKS_
0045   t->SetBranchAddress("rcparticles", &rctracks);
0046 #endif
0047   t->SetBranchAddress("PFRICHPID",   &cherenkov);
0048   auto options = new std::vector<edm4eic::CherenkovPdgHypothesis>();
0049   t->SetBranchAddress("PFRICHPID_0", &options);
0050
0051   // Loop through all events;
0052   unsigned false_assignment_stat = 0;
0053   for(int ev=0; ev<t->GetEntries(); ev++) {
0054     t->GetEntry(ev);
0055
0056 #ifdef _USE_RECONSTRUCTED_TRACKS_
0057     // First populate the reconstructed-to-simulated particle mapping table;
0058     std::map<eic::Index, const edm4eic::ReconstructedParticleData*> mc2rc;
0059     for(const auto &rctrack: *rctracks)
0060       mc2rc[rctrack.mcID] = &rctrack;
0061 #endif
0062
0063     // Then the Cherenkov-to-reconstructed mapping; FIXME: may want to use Cherenkov-to-simulated
0064     // mapping to start with, for the debugging purposes;
0065     std::map<eic::Index, const edm4eic::CherenkovParticleIDData*> rc2cherenkov;
0066     for(const auto &pid: *cherenkov)
0067       rc2cherenkov[pid.recID] = &pid;
0068
0069     // Loop through all MC tracks;
0070     for(auto mctrack: *mctracks) {
0071       // FIXME: consider only primaries for now?;
0072       if (mctrack.g4Parent) continue;
0073
0074 #ifdef _USE_RECONSTRUCTED_TRACKS_
0075       // Find a matching reconstructed track;
0076       auto rctrack = mc2rc.find(mctrack.ID) == mc2rc.end() ? 0 : mc2rc[mctrack.ID];
0077       if (!rctrack) continue;
0078
0079       // Find a matching Cherenkov PID record;
0080       auto cherenkov = rc2cherenkov.find(rctrack.ID) == rc2cherenkov.end() ? 0 : rc2cherenkov[rctrack.ID];
0081 #else
0082       auto cherenkov = rc2cherenkov.find(mctrack.ID) == rc2cherenkov.end() ? 0 : rc2cherenkov[mctrack.ID];
0083 #endif
0084       if (!cherenkov) continue;
0085
0086       // Loop through all of the mass hypotheses available for this reconstructed track;
0087       {
0088     const edm4eic::CherenkovPdgHypothesis *best = 0;
0089
0090     //printf("%d %d\n", cherenkov->options_begin, cherenkov->options_end);
0091     for(unsigned iq=cherenkov->options_begin; iq<cherenkov->options_end; iq++) {
0092       const auto &option = (*options)[iq];
0093
0094       // Skip electron hypothesis; of no interest here;
0095       if (abs(option.pdg) == 11) continue;
0096
0097       if (!best || option.weight > best->weight) best = &option;
0098       printf("radiator %3d (pdg %5d): weight %7.2f, npe %7.2f\n",
0099          option.radiator, option.pdg, option.weight, option.npe);
0100     } //for
0101     printf("\n");
0102
0103     // Check whether the true PDG got a highest score;
0104     if (!best || best->pdg != mctrack.pdgID) false_assignment_stat++;
0105       }
0106     } //for track
0107   } //for ev
0108
0109   printf("%d false out of %lld\n", false_assignment_stat, t->GetEntries());
0110
0111   return 0;
0112 } // main()