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