reconstruction/tof_efficiency/TofEfficiency_processor.cc

0001 #include "TofEfficiency_processor.h"
0002
0003 #include <JANA/JApplication.h>
0004 #include <JANA/Services/JGlobalRootLock.h>
0005 #include <edm4eic/TrackPoint.h>
0006 #include <edm4eic/TrackSegmentCollection.h>
0007 #include <edm4eic/TrackerHitCollection.h>
0008 #include <edm4hep/MCParticleCollection.h>
0009 #include <edm4hep/Vector3f.h>
0010 #include <fmt/core.h>
0011 #include <podio/RelationRange.h>
0012 #include <spdlog/logger.h>
0013 #include <cmath>
0014 #include <vector>
0015
0016 #include "services/rootfile/RootFile_service.h"
0017
0018 //-------------------------------------------
0019 // InitWithGlobalRootLock
0020 //-------------------------------------------
0021 void TofEfficiency_processor::InitWithGlobalRootLock(){
0022     std::string plugin_name=("tof_efficiency");
0023
0024     InitLogger(GetApplication(), plugin_name);
0025
0026     // Get JANA application
0027     auto *app = GetApplication();
0028
0029     // Ask service locator a file to write histograms to
0030     auto root_file_service = app->GetService<RootFile_service>();
0031
0032     // Get TDirectory for histograms root file
0033     auto globalRootLock = app->GetService<JGlobalRootLock>();
0034     globalRootLock->acquire_write_lock();
0035     auto *file = root_file_service->GetHistFile();
0036     globalRootLock->release_lock();
0037
0038     // Create a directory for this plugin. And subdirectories for series of histograms
0039     m_dir_main = file->mkdir(plugin_name.c_str());
0040
0041     auto phi_limit_min = 0;
0042     auto phi_limit_max = 6.29;
0043     auto z_limit_min = -1250;
0044     auto z_limit_max = 1250;
0045     auto r_limit_min = 50;
0046     auto r_limit_max = 675;
0047
0048     m_th2_btof_phiz = new TH2F("btof_phiz", "Hit position for Barrel TOF", 100, phi_limit_min, phi_limit_max, 100, z_limit_min, z_limit_max);
0049     m_th2_ftof_rphi = new TH2F("ftof_rphi", "Hit position for Forward TOF", 100, r_limit_min, r_limit_max, 100, phi_limit_min, phi_limit_max);
0050     m_th2_btof_phiz->SetDirectory(m_dir_main);
0051     m_th2_ftof_rphi->SetDirectory(m_dir_main);
0052
0053     m_tntuple_track = new TNtuple("track","track with tof","det:proj_x:proj_y:proj_z:proj_pathlength:tofhit_x:tofhit_y:tofhit_z:tofhit_t:tofhit_dca");
0054     m_tntuple_track->SetDirectory(m_dir_main);
0055 }
0056
0057 //-------------------------------------------
0058 // ProcessSequential
0059 //-------------------------------------------
0060 void TofEfficiency_processor::ProcessSequential(const std::shared_ptr<const JEvent>& event) {
0061     const auto &mcParticles   = *(event->GetCollection<edm4hep::MCParticle>("MCParticles"));
0062     const auto &trackSegments = *(event->GetCollection<edm4eic::TrackSegment>("CentralTrackSegments"));
0063     const auto &barrelHits    = *(event->GetCollection<edm4eic::TrackerHit>("TOFBarrelRecHit"));
0064     const auto &endcapHits    = *(event->GetCollection<edm4eic::TrackerHit>("TOFEndcapRecHits"));
0065
0066     // List TOF Barrel hits from barrel
0067     logger()->trace("TOF barrel hits:");
0068     m_log->trace("   {:>10} {:>10} {:>10} {:>10}", "[x]", "[y]", "[z]", "[time]");
0069     for (const auto hit: barrelHits) {
0070         const auto& pos = hit.getPosition();
0071         float r=sqrt(pos.x*pos.x+pos.y*pos.y);
0072         float phi=acos(pos.x/r); if(pos.y<0) phi+=3.1415927;
0073         m_th2_btof_phiz->Fill(phi, pos.z);
0074         m_log->trace("   {:>10.2f} {:>10.2f} {:>10.2f} {:>10.4f}", pos.x, pos.y, pos.z, hit.getTime());
0075     }
0076
0077     // List TOF endcap hits
0078     logger()->trace("TOF endcap hits:");
0079     m_log->trace("   {:>10} {:>10} {:>10} {:>10}", "[x]", "[y]", "[z]", "[time]");
0080     for (const auto hit: endcapHits) {
0081         const auto& pos = hit.getPosition();
0082         float r=sqrt(pos.x*pos.x+pos.y*pos.y);
0083         float phi=acos(pos.x/r); if(pos.y<0) phi+=3.1415927;
0084         m_th2_ftof_rphi->Fill(r, phi);
0085         m_log->trace("   {:>10.2f} {:>10.2f} {:>10.2f} {:>10.4f}", pos.x, pos.y, pos.z, hit.getTime());
0086     }
0087
0088     // Now go through reconstructed tracks points
0089     logger()->trace("Going over tracks:");
0090     m_log->trace("   {:>10} {:>10} {:>10} {:>10}", "[x]", "[y]", "[z]", "[length]");
0091     for( const auto track_segment : trackSegments ){
0092         logger()->trace(" Track trajectory");
0093
0094         for(const auto point: track_segment.getPoints()) {
0095             auto &pos = point.position;
0096             m_log->trace("   {:>10.2f} {:>10.2f} {:>10.2f} {:>10.2f}", pos.x, pos.y, pos.z, point.pathlength);
0097
0098             int det=IsTOFHit(pos.x, pos.y, pos.z);
0099             float distance_closest=1e6;
0100             float hit_x=-1000, hit_y=-1000, hit_z=-1000, hit_t=-1000;
0101             float hit_px=-1000, hit_py=-1000, hit_pz=-1000, hit_e=-1000;
0102             if(det==1) {
0103                 for (const auto hit: barrelHits) {
0104                     const auto& hitpos = hit.getPosition();
0105                     float distance=sqrt((hitpos.x-pos.x)*(hitpos.x-pos.x)+(hitpos.y-pos.y)*(hitpos.y-pos.y)+(hitpos.z-pos.z)*(hitpos.z-pos.z));
0106                     if(distance<distance_closest) {
0107                         distance_closest=distance;
0108                         hit_x=hitpos.x;
0109                         hit_y=hitpos.y;
0110                         hit_z=hitpos.z;
0111                         hit_t=hit.getTime();
0112                     }
0113                 }
0114             }
0115             if(det==2) {
0116                 for (const auto hit: endcapHits) {
0117                     const auto& hitpos = hit.getPosition();
0118                     float distance=sqrt((hitpos.x-pos.x)*(hitpos.x-pos.x)+(hitpos.y-pos.y)*(hitpos.y-pos.y)+(hitpos.z-pos.z)*(hitpos.z-pos.z));
0119                     if(distance<distance_closest) {
0120                         distance_closest=distance;
0121                         hit_x=hitpos.x;
0122                         hit_y=hitpos.y;
0123                         hit_z=hitpos.z;
0124                         hit_t=hit.getTime();
0125                     }
0126                 }
0127             }
0128             if(det!=0) m_tntuple_track->Fill(det, pos.x, pos.y, pos.z, point.pathlength, hit_x, hit_y, hit_z, hit_t, distance_closest);
0129         }
0130     }
0131 }
0132
0133 //-------------------------------------------
0134 // FinishWithGlobalRootLock
0135 //-------------------------------------------
0136 void TofEfficiency_processor::FinishWithGlobalRootLock() {
0137
0138     // Do any final calculations here.
0139
0140 }
0141
0142 int TofEfficiency_processor::IsTOFHit(float x, float y, float z) {
0143     const float btof_rmin=630;
0144     const float btof_rmax=642;
0145     const float btof_zmin=-1210;
0146     const float btof_zmax=+1210;
0147
0148     const float ftof_rmin=65;
0149     const float ftof_rmax=680;
0150     const float ftof_zmin=1900;
0151     const float ftof_zmax=1940;
0152
0153     float r=sqrt(x*x+y*y);
0154
0155     if(r>btof_rmin&&r<btof_rmax&&z>btof_zmin&&z<btof_zmax) return 1;
0156     else if(r>ftof_rmin&&r<ftof_rmax&&z>ftof_zmin&&z<ftof_zmax) return 2;
0157     else return 0;
0158 }