Merge pull request #920 from SBNSoftware/feature/hlay_crt_ana

CRT Analysis Tools

Author: nathanielerowe

sbndcode/CRT/CRTAna/ADRIFT_module.cc

@@ -1,12 +1,17 @@
 art_make(
          MODULE_LIBRARIES
          ROOT::Tree
+         ROOT::Graf3d
          art_root_io::TFileService_service
          larsim::Utils
          lardata::DetectorClocksService
+         sbndaq_artdaq_core::sbndaq-artdaq-core_Obj_SBND
          sbnobj::SBND_CRT
          sbnobj::SBND_Timing
+         sbnobj::Common_Reco
          sbndcode_CRT_CRTBackTracker
+         sbndcode_CRT_CRTReco
+         sbndcode_ChannelMaps_CRT_CRTChannelMapService_service
 )
 
 install_fhicl()

sbndcode/CRT/CRTAna/CMakeLists.txt

@@ -37,13 +37,13 @@
 #include "sbnobj/SBND/Timing/DAQTimestamp.hh"
 
 #include "sbndcode/Geometry/GeometryWrappers/CRTGeoService.h"
-#include "sbndcode/Geometry/GeometryWrappers/TPCGeoAlg.h"
 #include "sbndcode/ChannelMaps/CRT/CRTChannelMapService.h"
 #include "sbndcode/CRT/CRTBackTracker/CRTBackTrackerAlg.h"
 #include "sbndcode/CRT/CRTUtils/CRTCommonUtils.h"
 #include "sbndcode/CRT/CRTUtils/TPCGeoUtil.h"
 #include "sbndcode/Decoders/PTB/sbndptb.h"
 #include "sbndcode/Timing/SBNDRawTimingObj.h"
+#include "sbndcode/CRT/CRTReco/CRTClusterCharacterisationAlg.h"
 
 namespace sbnd::crt {
   class CRTAnalysis;
@@ -103,16 +103,14 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   art::ServiceHandle<CRTGeoService>              fCRTGeoService;
   art::ServiceHandle<SBND::CRTChannelMapService> fCRTChannelMapService;
 
-  TPCGeoAlg         fTPCGeoAlg;
-  CRTBackTrackerAlg fCRTBackTrackerAlg;
+  CRTBackTrackerAlg             fCRTBackTrackerAlg;
+  CRTClusterCharacterisationAlg fCRTClusterCharacAlg;
 
   std::string fMCParticleModuleLabel, fSimDepositModuleLabel, fFEBDataModuleLabel, fCRTStripHitModuleLabel,
     fCRTClusterModuleLabel, fCRTSpacePointModuleLabel, fCRTTrackModuleLabel, fCRTBlobModuleLabel, fTPCTrackModuleLabel,
     fCRTSpacePointMatchingModuleLabel, fCRTTrackMatchingModuleLabel, fPFPModuleLabel, fPTBModuleLabel,
     fTDCModuleLabel, fTimingReferenceModuleLabel;
   bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC, fHasBlobs, fTruthMatch;
-  //! Adding some of the reco parameters to save corrections
-  double fPEAttenuation, fTimeWalkNorm, fTimeWalkScale, fPropDelay;
 
   TTree* fTree;
 
@@ -124,7 +122,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   int _crt_timing_reference_type;
   int _crt_timing_reference_channel;
 
-  //mc truth
+  // True Particles (from G4)
   std::vector<int16_t>              _mc_trackid;
   std::vector<int16_t>              _mc_pdg;
   std::vector<int16_t>              _mc_status;
@@ -147,7 +145,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>               _mc_endpz;
   std::vector<double>               _mc_ende;
 
-  //G4 detector id
+  // True Energy Depositions
   std::vector<int16_t> _ide_trackid;
   std::vector<float>   _ide_e;
   std::vector<float>   _ide_entryx;
@@ -159,7 +157,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<float>   _ide_exitz;
   std::vector<float>   _ide_exitt;
 
-  //front end mother board
+  // Raw Readouts (FEBDatas)
   std::vector<uint16_t>              _feb_mac5;
   std::vector<int16_t>               _feb_tagger;
   std::vector<uint16_t>              _feb_flags;
@@ -169,7 +167,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<std::vector<uint16_t>> _feb_adc;
   std::vector<uint32_t>              _feb_coinc;
 
-  //strip hit to select the strip which has ADC above threshold
+  // Strip Hits
   std::vector<uint32_t> _sh_channel;
   std::vector<int16_t>  _sh_tagger;
   std::vector<double>   _sh_ts0;
@@ -188,7 +186,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>   _sh_truth_energy;
   std::vector<double>   _sh_truth_time;
 
-  //cluster from x-y coincidence for CRTSpacePoint , this is what we normally call a CRT hit
+  // Clusters and their corresponding SpacePoints
   std::vector<double>                _cl_ts0;
   std::vector<double>                _cl_ts1;
   std::vector<uint32_t>              _cl_unixs;
@@ -197,11 +195,11 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<uint8_t>               _cl_composition;
   std::vector<std::vector<uint32_t>> _cl_channel_set;
   std::vector<std::vector<uint16_t>> _cl_adc_set;
-  std::vector<std::vector<double>>   _cl_sh_ts0_set; //! To store t0 from x-y coincidences
-  std::vector<std::vector<double>>   _cl_sh_ts1_set; //! To store t1 from x-y coincidences
-  std::vector<std::vector<uint16_t>> _cl_sh_feb_mac5_set; //! MAC5 addresses of StripHit FEBs
-  std::vector<std::vector<double>>   _cl_sh_time_walk_set; //! Time walk correction
-  std::vector<std::vector<double>>   _cl_sh_prop_delay_set; //! Light propagation correction
+  std::vector<std::vector<double>>   _cl_sh_ts0_set;
+  std::vector<std::vector<double>>   _cl_sh_ts1_set;
+  std::vector<std::vector<uint16_t>> _cl_sh_feb_mac5_set;
+  std::vector<std::vector<double>>   _cl_sh_time_walk_set;
+  std::vector<std::vector<double>>   _cl_sh_prop_delay_set;
   std::vector<int>                   _cl_truth_trackid;
   std::vector<double>                _cl_truth_completeness;
   std::vector<double>                _cl_truth_purity;
@@ -234,7 +232,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>                _cl_sp_dts1;
   std::vector<bool>                  _cl_sp_complete;
 
-  //backtrack truth information from reco level
+  // True Deposits per particle per tagger
   std::vector<int>     _td_tag_trackid;
   std::vector<int>     _td_tag_pdg;
   std::vector<int16_t> _td_tag_tagger;
@@ -245,6 +243,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>  _td_tag_z;
   std::vector<bool>    _td_tag_reco_status;
 
+  // True Deposits per particle
   std::vector<int>     _td_trackid;
   std::vector<int>     _td_pdg;
   std::vector<double>  _td_energy;
@@ -253,7 +252,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<bool>    _td_reco_status;
   std::vector<bool>    _td_reco_triple;
 
-  //track level information
+  // Tracks
   std::vector<double>                _tr_start_x;
   std::vector<double>                _tr_start_y;
   std::vector<double>                _tr_start_z;
@@ -299,7 +298,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>                _tr_truth_theta;
   std::vector<double>                _tr_truth_phi;
 
-  // crt blob information
+  // Blobs
   std::vector<double>           _bl_ts0;
   std::vector<double>           _bl_ets0;
   std::vector<double>           _bl_ts1;
@@ -308,7 +307,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<int>              _bl_nsps;
   std::vector<std::vector<int>> _bl_nsps_per_tagger;
 
-  // tpc track information (including crt matching)
+  // TPC Tracks and their CRT matches
   std::vector<double>                _tpc_start_x;
   std::vector<double>                _tpc_start_y;
   std::vector<double>                _tpc_start_z;
@@ -323,6 +322,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>                _tpc_end_dir_z;
   std::vector<double>                _tpc_length;
   std::vector<double>                _tpc_track_score;
+  std::vector<int>                   _tpc_whichtpc;
   std::vector<int>                   _tpc_truth_trackid;
   std::vector<int>                   _tpc_truth_pdg;
   std::vector<double>                _tpc_truth_energy;
@@ -343,6 +343,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<bool>                  _tpc_tr_matchable;
   std::vector<bool>                  _tpc_tr_matched;
   std::vector<bool>                  _tpc_tr_good_match;
+  std::vector<double>                _tpc_tr_xshift;
   std::vector<double>                _tpc_tr_ts0;
   std::vector<double>                _tpc_tr_ts1;
   std::vector<std::vector<int>>      _tpc_tr_taggers;
@@ -354,14 +355,15 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>                _tpc_tr_end_z;
   std::vector<double>                _tpc_tr_score;
 
-  // ptb information (trigger board)
+  // Penn Trigger Board HLTs
   std::vector<uint64_t> _ptb_hlt_trigger;
   std::vector<uint64_t> _ptb_hlt_timestamp;
 
+  // Penn Trigger Board LLTs
   std::vector<uint64_t> _ptb_llt_trigger;
   std::vector<uint64_t> _ptb_llt_timestamp;
 
-  // spec tdc information (timing board)
+  // SPEC TDC Timestamps
   std::vector<uint32_t>    _tdc_channel;
   std::vector<uint64_t>    _tdc_timestamp;
   std::vector<uint64_t>    _tdc_offset;
@@ -371,6 +373,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
 sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   : EDAnalyzer{p}
   , fCRTBackTrackerAlg(p.get<fhicl::ParameterSet>("CRTBackTrackerAlg", fhicl::ParameterSet()))
+  , fCRTClusterCharacAlg(p.get<fhicl::ParameterSet>("CRTClusterCharacterisationAlg", fhicl::ParameterSet()))
 {
   fMCParticleModuleLabel            = p.get<std::string>("MCParticleModuleLabel", "largeant");
   fSimDepositModuleLabel            = p.get<std::string>("SimDepositModuleLabel", "genericcrt");
@@ -394,10 +397,6 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   fHasTDC                           = p.get<bool>("HasTDC", false);
   fHasBlobs                         = p.get<bool>("HasBlobs", false);
   fTruthMatch                       = p.get<bool>("TruthMatch", true);
-  fPEAttenuation                    = p.get<double>("PEAttenuation", 1.0);
-  fTimeWalkNorm                     = p.get<double>("TimeWalkNorm",  0.0);
-  fTimeWalkScale                    = p.get<double>("TimeWalkScale", 0.0);
-  fPropDelay                        = p.get<double>("PropDelay",     0.0);
 
   if(!fDataMode && fTruthMatch)
     fCRTBackTrackerAlg = CRTBackTrackerAlg(p.get<fhicl::ParameterSet>("CRTBackTrackerAlg", fhicl::ParameterSet()));
@@ -621,6 +620,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
       fTree->Branch("tpc_end_dir_z", "std::vector<double>", &_tpc_end_dir_z);
       fTree->Branch("tpc_length", "std::vector<double>", &_tpc_length);
       fTree->Branch("tpc_track_score", "std::vector<double>", &_tpc_track_score);
+      fTree->Branch("tpc_whichtpc", "std::vector<int>", &_tpc_whichtpc);
       fTree->Branch("tpc_sp_matched", "std::vector<bool>", &_tpc_sp_matched);
       fTree->Branch("tpc_sp_channel_set", "std::vector<std::vector<uint32_t>>", &_tpc_sp_channel_set);
       fTree->Branch("tpc_sp_xshift", "std::vector<double>", &_tpc_sp_xshift);
@@ -633,6 +633,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
       fTree->Branch("tpc_sp_z", "std::vector<double>", &_tpc_sp_z);
       fTree->Branch("tpc_sp_score", "std::vector<double>", &_tpc_sp_score);
       fTree->Branch("tpc_tr_matched", "std::vector<bool>", &_tpc_tr_matched);
+      fTree->Branch("tpc_tr_xshift", "std::vector<double>", &_tpc_tr_xshift);
       fTree->Branch("tpc_tr_ts0", "std::vector<double>", &_tpc_tr_ts0);
       fTree->Branch("tpc_tr_ts1", "std::vector<double>", &_tpc_tr_ts1);
       fTree->Branch("tpc_tr_taggers", "std::vector<std::vector<int>>", &_tpc_tr_taggers);
@@ -1326,25 +1327,10 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
           _cl_sh_ts1_set[i][ii]      = striphit->Ts1();
           _cl_sh_feb_mac5_set[i][ii] = fCRTChannelMapService->GetMAC5FromOfflineChannelID(striphit->Channel());
 
-          /*
-           * The below segment reimplements the CorrectTime() method
-           * from CRTReco/CRTClusterCharacterisationAlg.cc .
-           * Because the Ts0(), Ts1() getters invoked in _cl_sp_ts*, _cl_sh_ts*_set are raw T0/1
-           * counters, the time walk and propagation delay are saved as explicit branches here.
-           */
           if(spacepoints.size() == 1)
             {
-              double pe0 = fCRTGeoService->GetSiPM( striphit->Channel() ).gain * striphit->ADC1();
-              double pe1 = fCRTGeoService->GetSiPM( striphit->Channel() + 1 ).gain * striphit->ADC2();
-              double pe  = pe0 + pe1;
-
-              double dist = fCRTGeoService->DistanceDownStrip( spacepoints[0]->Pos(), striphit->Channel() );
-
-              double corr = std::pow( dist - fPEAttenuation, 2.0 ) / std::pow( fPEAttenuation, 2.0 );
-              double tw_pe = pe * corr;
-
-              _cl_sh_time_walk_set[i][ii]  = fTimeWalkNorm * std::exp( -fTimeWalkScale * tw_pe );
-              _cl_sh_prop_delay_set[i][ii] = fPropDelay * dist;
+              _cl_sh_time_walk_set[i][ii]  = fCRTClusterCharacAlg.TimeWalk(striphit, spacepoints[0]->Pos());
+              _cl_sh_prop_delay_set[i][ii] = fCRTClusterCharacAlg.PropagationDelay(striphit, spacepoints[0]->Pos());
 
               ts0_set.push_back({_cl_sh_feb_mac5_set[i][ii], _cl_sh_ts0_set[i][ii] - _cl_sh_time_walk_set[i][ii] - _cl_sh_prop_delay_set[i][ii]});
               ts1_set.push_back({_cl_sh_feb_mac5_set[i][ii], _cl_sh_ts1_set[i][ii] - _cl_sh_time_walk_set[i][ii] - _cl_sh_prop_delay_set[i][ii]});
@@ -1682,6 +1668,7 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::
   _tpc_end_dir_z.resize(nTracks);
   _tpc_length.resize(nTracks);
   _tpc_track_score.resize(nTracks);
+  _tpc_whichtpc.resize(nTracks);
   _tpc_truth_trackid.resize(nTracks);
   _tpc_truth_pdg.resize(nTracks);
   _tpc_truth_energy.resize(nTracks);
@@ -1702,6 +1689,7 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::
   _tpc_tr_matchable.resize(nTracks);
   _tpc_tr_matched.resize(nTracks);
   _tpc_tr_good_match.resize(nTracks);
+  _tpc_tr_xshift.resize(nTracks);
   _tpc_tr_ts0.resize(nTracks);
   _tpc_tr_ts1.resize(nTracks);
   _tpc_tr_taggers.resize(nTracks);
@@ -1771,6 +1759,7 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::
       const art::Ptr<CRTTrack> crttrack = tracksToTrackMatches.at(track.key());
 
       const std::vector<art::Ptr<recob::Hit>> trackHits = tracksToHits.at(track.key());
+      _tpc_whichtpc[nActualTracks]                      = TPCGeoUtil::DetectedInTPC(trackHits);
 
       if(spacepoint.isNonnull())
         {
@@ -1813,7 +1802,11 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::
         {
           const anab::T0 trackMatch = tracksToTrackMatches.data(track.key()).ref();
 
+          const int driftDirection     = TPCGeoUtil::DriftDirectionFromHits(geometryService, trackHits);
+          const double crtShiftingTime = fDataMode ? crttrack->Ts0() * 1e-3 : crttrack->Ts1() * 1e-3;
+
           _tpc_tr_matched[nActualTracks] = true;
+          _tpc_tr_xshift[nActualTracks]  = driftDirection * crtShiftingTime * detProp.DriftVelocity();
           _tpc_tr_ts0[nActualTracks]     = crttrack->Ts0();
           _tpc_tr_ts1[nActualTracks]     = crttrack->Ts1();
           _tpc_tr_score[nActualTracks]   = trackMatch.TriggerConfidence();
@@ -1837,6 +1830,7 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::
       else
         {
           _tpc_tr_matched[nActualTracks] = false;
+          _tpc_tr_xshift[nActualTracks]  = -std::numeric_limits<double>::max();
           _tpc_tr_ts0[nActualTracks]     = -std::numeric_limits<double>::max();
           _tpc_tr_ts1[nActualTracks]     = -std::numeric_limits<double>::max();
           _tpc_tr_score[nActualTracks]   = -std::numeric_limits<double>::max();
@@ -1915,6 +1909,7 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::
   _tpc_end_dir_z.resize(nActualTracks);
   _tpc_length.resize(nActualTracks);
   _tpc_track_score.resize(nActualTracks);
+  _tpc_whichtpc.resize(nActualTracks);
   _tpc_truth_trackid.resize(nActualTracks);
   _tpc_truth_pdg.resize(nActualTracks);
   _tpc_truth_energy.resize(nActualTracks);