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 <stddef.h>
0010 #include <cmath>
0011 #include <exception>
0012 #include <map>
0013 #include <memory>
0014
0015 #include "TrackingOccupancyAnalysis.h"
0016
0017
0018 void TrackingOccupancyAnalysis::init(JApplication *app, TDirectory *plugin_tdir) {
0019     auto *dir = plugin_tdir->mkdir("SimOccupancies");     // 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
0027     m_total_occup_th2 = new TH2F("total_occup", "Occupancy plot for all readouts", 200, z_limit_min, +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(), ("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 = std::make_shared<TH2F>(("occup_" + name).c_str(), ("Occupancy plot for" + name).c_str(), 100, z_limit_min, z_limit_max, 200, r_limit_min, r_limit_max);
0036         occup_hist->SetDirectory(dir);
0037         m_hits_occup_hists.push_back(occup_hist);
0038     }
0039 }
0040
0041
0042 void TrackingOccupancyAnalysis::process(const std::shared_ptr<const JEvent> &event) {
0043
0044     for(size_t name_index = 0; name_index < m_data_names.size(); name_index++ ) {
0045         std::string data_name = m_data_names[name_index];
0046         auto &count_hist = m_hits_count_hists[name_index];
0047         auto &occup_hist = m_hits_occup_hists[name_index];
0048
0049         try {
0050             auto hits = event->Get<edm4hep::SimTrackerHit>(data_name);
0051             count_hist->Fill(hits.size());
0052             for(const auto *hit: hits) {
0053                 float x = hit->getPosition().x;
0054                 float y = hit->getPosition().y;
0055                 float z = hit->getPosition().z;
0056                 float r = sqrt(x*x + y*y);
0057                 occup_hist->Fill(z, r);
0058                 m_total_occup_th2->Fill(z, r);
0059             }
0060         } catch(std::exception& e) {
0061             // silently skip missing collections
0062         }
0063     }
0064 }