reconstruction/tracking_efficiency/TrackingEfficiency_processor.cc

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