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 <stddef.h>
0010 #include <cmath>
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, +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(), ("Count hits for " + name).c_str(), 100, 0, 30);
0031         count_hist->SetDirectory(dir);
0032         m_hits_count_hists.push_back(count_hist);
0033
0034         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);
0035         occup_hist->SetDirectory(dir);
0036         m_hits_occup_hists.push_back(occup_hist);
0037     }
0038 }
0039
0040 void HitReconstructionAnalysis::process(const std::shared_ptr<const JEvent> &event) {
0041     for(size_t name_index = 0; name_index < m_data_names.size(); name_index++ ) {
0042         std::string data_name = m_data_names[name_index];
0043         auto &count_hist = m_hits_count_hists[name_index];
0044         auto &occup_hist = m_hits_occup_hists[name_index];
0045
0046         try {
0047             auto hits = event->Get<edm4eic::TrackerHit>(data_name);
0048             count_hist->Fill(hits.size());
0049             for(const auto *hit: hits) {
0050                 float x = hit->getPosition().x;
0051                 float y = hit->getPosition().y;
0052                 float z = hit->getPosition().z;
0053                 float r = sqrt(x*x + y*y);
0054                 occup_hist->Fill(z, r);
0055                 m_total_occup_th2->Fill(z, r);
0056             }
0057         } catch(std::exception& e) {
0058             // silently skip missing collections
0059         }
0060     }
0061 }