reconstruction/tracking_efficiency/TrackingEfficiency_processor.cc

0001 #include "TrackingEfficiency_processor.h"
0002
0003 #include <Acts/Definitions/TrackParametrization.hpp>
0004 #include <Acts/EventData/TrackProxy.hpp>
0005 #include <Acts/EventData/VectorMultiTrajectory.hpp>
0006 #include <Acts/EventData/VectorTrackContainer.hpp>
0007 #include <ActsExamples/EventData/Track.hpp>
0008 #include <JANA/JApplication.h>
0009 #include <JANA/JApplicationFwd.h>
0010 #include <JANA/JEvent.h>
0011 #include <JANA/Services/JGlobalRootLock.h>
0012 #include <Math/GenVector/Cartesian3D.h>
0013 #include <Math/GenVector/PxPyPzM4D.h>
0014 #include <edm4eic/ReconstructedParticleCollection.h>
0015 #include <edm4hep/MCParticleCollection.h>
0016 #include <edm4hep/Vector3d.h>
0017 #include <edm4hep/Vector3f.h>
0018 #include <fmt/format.h>
0019 #include <spdlog/logger.h>
0020 #include <cassert>
0021 #include <cmath>
0022 #include <cstddef>
0023 #include <map>
0024 #include <string>
0025 #include <vector>
0026
0027 #include "services/log/Log_service.h"
0028 #include "services/rootfile/RootFile_service.h"
0029
0030 //------------------
0031 // Init
0032 //------------------
0033 void TrackingEfficiency_processor::Init() {
0034   std::string plugin_name = ("tracking_efficiency");
0035
0036   // Get JANA application
0037   auto* app = GetApplication();
0038
0039   // Ask service locator a file to write histograms to
0040   auto root_file_service = app->GetService<RootFile_service>();
0041
0042   // Get TDirectory for histograms root file
0043   auto globalRootLock = app->GetService<JGlobalRootLock>();
0044   globalRootLock->acquire_write_lock();
0045   auto* file = root_file_service->GetHistFile();
0046   globalRootLock->release_lock();
0047
0048   // Create a directory for this plugin. And subdirectories for series of histograms
0049   m_dir_main = file->mkdir(plugin_name.c_str());
0050
0051   // Get logger
0052   m_log = app->GetService<Log_service>()->logger(plugin_name);
0053 }
0054
0055 //------------------
0056 // Process
0057 //------------------
0058 void TrackingEfficiency_processor::Process(const std::shared_ptr<const JEvent>& event) {
0059   using namespace ROOT;
0060
0061   // EXAMPLE I
0062   // This is access to for final result of the calculation/data transformation of central detector CFKTracking:
0063   const auto* reco_particles =
0064       event->GetCollection<edm4eic::ReconstructedParticle>("ReconstructedChargedParticles");
0065
0066   m_log->debug("Tracking reconstructed particles N={}: ", reco_particles->size());
0067   m_log->debug("   {:<5} {:>8} {:>8} {:>8} {:>8} {:>8}", "[i]", "[px]", "[py]", "[pz]", "[P]",
0068                "[P*3]");
0069
0070   for (std::size_t i = 0; i < reco_particles->size(); i++) {
0071     const auto& particle = (*reco_particles)[i];
0072
0073     double px = particle.getMomentum().x;
0074     double py = particle.getMomentum().y;
0075     double pz = particle.getMomentum().z;
0076     // ROOT::Math::PxPyPzM4D p4v(px, py, pz, particle.getMass());
0077     ROOT::Math::Cartesian3D p(px, py, pz);
0078     m_log->debug("   {:<5} {:>8.2f} {:>8.2f} {:>8.2f} {:>8.2f} {:>8.2f}", i, px, py, pz, p.R(),
0079                  p.R() * 3);
0080   }
0081
0082   // EXAMPLE II
0083   // This gets access to more direct ACTS results from CKFTracking
0084   auto acts_track_states =
0085       event->Get<Acts::ConstVectorMultiTrajectory>("CentralCKFActsTrackStates");
0086   auto acts_tracks = event->Get<Acts::ConstVectorTrackContainer>("CentralCKFActsTracks");
0087   m_log->debug("ACTS Tracks( track states size: {}, tracks size: {} )", acts_track_states.size(),
0088                acts_tracks.size());
0089   m_log->debug("{:>10} {:>10}  {:>10} {:>10} {:>10} {:>10} {:>12} {:>12} {:>12} {:>8} {:>8}",
0090                "[loc 0]", "[loc 1]", "[phi]", "[theta]", "[q/p]", "[p]", "[err phi]", "[err th]",
0091                "[err q/p]", "[chi2]", "[ndf]");
0092
0093   // Loop over the tracks
0094   if (!acts_track_states.empty() && !acts_tracks.empty()) {
0095     assert(acts_track_states.front() != nullptr &&
0096            "ConstVectorMultiTrajectory pointer should not be null");
0097     assert(acts_tracks.front() != nullptr &&
0098            "ConstVectorTrackContainer pointer should not be null");
0099
0100     // Construct ConstTrackContainer from underlying containers
0101     auto trackStateContainer =
0102         std::make_shared<Acts::ConstVectorMultiTrajectory>(*acts_track_states.front());
0103     auto trackContainer = std::make_shared<Acts::ConstVectorTrackContainer>(*acts_tracks.front());
0104     ActsExamples::ConstTrackContainer track_container(trackContainer, trackStateContainer);
0105
0106     for (const auto& track : track_container) {
0107       // Get the track parameters
0108       const auto& parameter  = track.parameters();
0109       const auto& covariance = track.covariance();
0110       auto chi2              = track.chi2();
0111       auto ndf               = track.nDoF();
0112
0113       m_log->debug("{:>10.2f} {:>10.2f}  {:>10.2f} {:>10.3f} {:>10.4f} {:>10.3f} {:>12.4e} "
0114                    "{:>12.4e} {:>12.4e} {:>8.2f} {:<6}",
0115                    parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1],
0116                    parameter[Acts::eBoundPhi], parameter[Acts::eBoundTheta],
0117                    parameter[Acts::eBoundQOverP], 1.0 / parameter[Acts::eBoundQOverP],
0118                    sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)),
0119                    sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
0120                    sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)), chi2, ndf);
0121     }
0122   }
0123
0124   // EXAMPLE III
0125   // Loop over MC particles
0126   const auto* mc_particles = event->GetCollection<edm4hep::MCParticle>("MCParticles");
0127   m_log->debug("MC particles N={}: ", mc_particles->size());
0128   m_log->debug("   {:<5} {:<6} {:<7} {:>8} {:>8} {:>8} {:>8}", "[i]", "status", "[PDG]", "[px]",
0129                "[py]", "[pz]", "[P]");
0130   for (std::size_t i = 0; i < mc_particles->size(); i++) {
0131     const auto& particle = (*mc_particles)[i];
0132
0133     // GeneratorStatus() == 1 - stable particles from MC generator. 0 - might be added by Geant4
0134     if (particle.getGeneratorStatus() != 1) {
0135       continue;
0136     }
0137
0138     double px = particle.getMomentum().x;
0139     double py = particle.getMomentum().y;
0140     double pz = particle.getMomentum().z;
0141     ROOT::Math::PxPyPzM4D p4v(px, py, pz, particle.getMass());
0142     ROOT::Math::Cartesian3D p(px, py, pz);
0143     if (p.R() < 1) {
0144       continue;
0145     }
0146
0147     m_log->debug("   {:<5} {:<6} {:<7} {:>8.2f} {:>8.2f} {:>8.2f} {:>8.2f}", i,
0148                  particle.getGeneratorStatus(), particle.getPDG(), px, py, pz, p.R());
0149   }
0150 }
0151
0152 //------------------
0153 // Finish
0154 //------------------
0155 void TrackingEfficiency_processor::Finish() {
0156   fmt::print("OccupancyAnalysis::Finish() called\n");
0157
0158   // Next we want to create several pretty canvases (with histograms drawn on "same")
0159   // But we don't want those canvases to pop up. So we set root to batch mode
0160   // We will restore the mode afterwards
0161   //bool save_is_batch = gROOT->IsBatch();
0162   //gROOT->SetBatch(true);
0163
0164   // 3D hits distribution
0165   //      auto th3_by_det_canvas = new TCanvas("th3_by_det_cnv", "Occupancy of detectors");
0166   //      dir_main->Append(th3_by_det_canvas);
0167   //      for (auto& kv : th3_by_detector->GetMap()) {
0168   //              auto th3_hist = kv.second;
0169   //              th3_hist->Draw("same");
0170   //      }
0171   //      th3_by_det_canvas->GetPad(0)->BuildLegend();
0172   //
0173   //      // Hits Z by detector
0174   //
0175   //      // Create pretty canvases
0176   //      auto z_by_det_canvas = new TCanvas("z_by_det_cnv", "Hit Z distribution by detector");
0177   //      dir_main->Append(z_by_det_canvas);
0178   //      th1_hits_z->Draw("PLC PFC");
0179   //
0180   //      for (auto& kv : th1_z_by_detector->GetMap()) {
0181   //              auto hist = kv.second;
0182   //              hist->Draw("SAME PLC PFC");
0183   //              hist->SetFillStyle(3001);
0184   //      }
0185   //      z_by_det_canvas->GetPad(0)->BuildLegend();
0186   //
0187   //      gROOT->SetBatch(save_is_batch);
0188 }