reconstruction/tracking_occupancy/TrackingOccupancyAnalysis.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 <edm4hep/SimTrackerHitCollection.h>
0008 #include <edm4hep/Vector3d.h>
0009 #include <cmath>
0010 #include <cstddef>
0011 #include <exception>
0012 #include <map>
0013 #include <memory>
0014
0015 #include "TrackingOccupancyAnalysis.h"
0016
0017 void TrackingOccupancyAnalysis::init(JApplication* /* app */, TDirectory* plugin_tdir) {
0018   auto* dir = plugin_tdir->mkdir("SimOccupancies"); // TODO create directory for this analysis
0019
0020   auto z_limit_min = -2000;
0021   auto z_limit_max = 2000;
0022   auto r_limit_min = 0;
0023   auto r_limit_max = 1200;
0024
0025   m_total_occup_th2 = new TH2F("total_occup", "Occupancy plot for all readouts", 200, z_limit_min,
0026                                +z_limit_max, 100, r_limit_min, r_limit_max);
0027   m_total_occup_th2->SetDirectory(dir);
0028
0029   for (auto& name : m_data_names) {
0030     auto count_hist = std::make_shared<TH1F>(("count_" + name).c_str(),
0031                                              ("Count hits for " + name).c_str(), 100, 0, 30);
0032     count_hist->SetDirectory(dir);
0033     m_hits_count_hists.push_back(count_hist);
0034
0035     auto occup_hist =
0036         std::make_shared<TH2F>(("occup_" + name).c_str(), ("Occupancy plot for" + name).c_str(),
0037                                100, z_limit_min, z_limit_max, 200, r_limit_min, r_limit_max);
0038     occup_hist->SetDirectory(dir);
0039     m_hits_occup_hists.push_back(occup_hist);
0040   }
0041 }
0042
0043 void TrackingOccupancyAnalysis::process(const std::shared_ptr<const JEvent>& event) {
0044
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<edm4hep::SimTrackerHit>(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 }