reconstruction/tracking_occupancy/HitReconstructionAnalysis.cc

0001 // Created by Dmitry Romanov
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003 //
0004
0005 #include <JANA/JEvent.h>
0006 #include <TDirectory.h>
0007 #include <edm4eic/TrackerHitCollection.h>
0008 #include <edm4hep/Vector3f.h>
0009 #include <cmath>
0010 #include <cstddef>
0011 #include <exception>
0012 #include <map>
0013 #include <memory>
0014
0015 #include "HitReconstructionAnalysis.h"
0016
0017 void HitReconstructionAnalysis::init(JApplication* /* app */, TDirectory* plugin_tdir) {
0018
0019   auto* dir = plugin_tdir->mkdir("RecOccupancies"); // TODO create directory for this analysis
0020
0021   auto z_limit_min = -2000;
0022   auto z_limit_max = 2000;
0023   auto r_limit_min = 0;
0024   auto r_limit_max = 1200;
0025
0026   m_total_occup_th2 = new TH2F("total_occup", "Occupancy plot for all readouts", 200, z_limit_min,
0027                                +z_limit_max, 100, r_limit_min, r_limit_max);
0028   m_total_occup_th2->SetDirectory(dir);
0029
0030   for (auto& name : m_data_names) {
0031     auto count_hist = std::make_shared<TH1F>(("count_" + name).c_str(),
0032                                              ("Count hits for " + name).c_str(), 100, 0, 30);
0033     count_hist->SetDirectory(dir);
0034     m_hits_count_hists.push_back(count_hist);
0035
0036     auto occup_hist =
0037         std::make_shared<TH2F>(("occup_" + name).c_str(), ("Occupancy plot for" + name).c_str(),
0038                                100, z_limit_min, z_limit_max, 200, r_limit_min, r_limit_max);
0039     occup_hist->SetDirectory(dir);
0040     m_hits_occup_hists.push_back(occup_hist);
0041   }
0042 }
0043
0044 void HitReconstructionAnalysis::process(const std::shared_ptr<const JEvent>& event) {
0045   for (std::size_t name_index = 0; name_index < m_data_names.size(); name_index++) {
0046     std::string data_name = m_data_names[name_index];
0047     auto& count_hist      = m_hits_count_hists[name_index];
0048     auto& occup_hist      = m_hits_occup_hists[name_index];
0049
0050     try {
0051       auto hits = event->Get<edm4eic::TrackerHit>(data_name);
0052       count_hist->Fill(hits.size());
0053       for (const auto* hit : hits) {
0054         float x = hit->getPosition().x;
0055         float y = hit->getPosition().y;
0056         float z = hit->getPosition().z;
0057         float r = sqrt(x * x + y * y);
0058         occup_hist->Fill(z, r);
0059         m_total_occup_th2->Fill(z, r);
0060       }
0061     } catch (std::exception& e) {
0062       // silently skip missing collections
0063     }
0064   }
0065 }