Add relevant value to CRTAnalysis module

Author: henrylay97

sbndcode/CRT/CRTAna/CRTAnalysis_module.cc

@@ -87,6 +87,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
                           const art::Handle<std::vector<CRTSpacePoint>> &CRTSpacePointModuleLabel, const art::Handle<std::vector<recob::PFParticle>> &PFPHandle,
                           const std::map<CRTBackTrackerAlg::Category, bool> &spacePointRecoStatusMap = {}, const std::map<int, std::pair<bool, bool>> &trackRecoStatusMap = {});
 
+  double SPTimeDelta(std::vector<std::pair<int, double>> &time_set);
+
 private:
 
   CRTGeoAlg fCRTGeoAlg;
@@ -217,6 +219,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>                _cl_sp_ets0;
   std::vector<double>                _cl_sp_ts1;
   std::vector<double>                _cl_sp_ets1;
+  std::vector<double>                _cl_sp_dts0;
+  std::vector<double>                _cl_sp_dts1;
   std::vector<bool>                  _cl_sp_complete;
 
   //backtrack truth information from reco level
@@ -476,6 +480,8 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   fTree->Branch("cl_sp_ets0", "std::vector<double>", &_cl_sp_ets0);
   fTree->Branch("cl_sp_ts1", "std::vector<double>", &_cl_sp_ts1);
   fTree->Branch("cl_sp_ets1", "std::vector<double>", &_cl_sp_ets1);
+  fTree->Branch("cl_sp_dts0", "std::vector<double>", &_cl_sp_dts0);
+  fTree->Branch("cl_sp_dts1", "std::vector<double>", &_cl_sp_dts1);
   fTree->Branch("cl_sp_complete", "std::vector<bool>", &_cl_sp_complete);
 
   if(!fDataMode)
@@ -1144,6 +1150,8 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
   _cl_sp_ets0.resize(nClusters);
   _cl_sp_ts1.resize(nClusters);
   _cl_sp_ets1.resize(nClusters);
+  _cl_sp_dts0.resize(nClusters);
+  _cl_sp_dts1.resize(nClusters);
   _cl_sp_complete.resize(nClusters);
 
   for(unsigned i = 0; i < nClusters; ++i)
@@ -1162,12 +1170,14 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
       const auto striphits = clustersToStripHits.at(cluster.key());
       _cl_channel_set[i].resize(_cl_nhits[i]);
       _cl_adc_set[i].resize(2 * _cl_nhits[i]);
-      _cl_sh_ts0_set[i].resize(_cl_nhits[i] );
-      _cl_sh_ts1_set[i].resize(_cl_nhits[i] );
+      _cl_sh_ts0_set[i].resize(_cl_nhits[i]);
+      _cl_sh_ts1_set[i].resize(_cl_nhits[i]);
       _cl_sh_feb_mac5_set[i].resize(_cl_nhits[i]);
       _cl_sh_time_walk_set[i].resize(_cl_nhits[i]);
       _cl_sh_prop_delay_set[i].resize(_cl_nhits[i]);
 
+      std::vector<std::pair<int, double>> ts0_set, ts1_set;
+
       for(unsigned ii = 0; ii < _cl_nhits[i]; ++ii)
         {
           const auto striphit = striphits[ii];
@@ -1192,24 +1202,28 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
            * 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) { // need unique position of spacepoint
-            double pe0 = fCRTGeoAlg.GetSiPM( striphit->Channel() ).gain * striphit->ADC1();
-            double pe1 = fCRTGeoAlg.GetSiPM( striphit->Channel() + 1 ).gain * striphit->ADC2();
-            double pe  = pe0 + pe1;
-
-            double dist = fCRTGeoAlg.DistanceDownStrip( spacepoints[0]->Pos(), striphit->Channel() );
+          if(spacepoints.size() == 1)
+            {
+              double pe0 = fCRTGeoAlg.GetSiPM( striphit->Channel() ).gain * striphit->ADC1();
+              double pe1 = fCRTGeoAlg.GetSiPM( striphit->Channel() + 1 ).gain * striphit->ADC2();
+              double pe  = pe0 + pe1;
 
-            double corr = std::pow( dist - fPEAttenuation, 2.0 ) / std::pow( fPEAttenuation, 2.0 );
-            double tw_pe = pe * corr;
+              double dist = fCRTGeoAlg.DistanceDownStrip( spacepoints[0]->Pos(), striphit->Channel() );
 
-            _cl_sh_time_walk_set[i][ii]  = fTimeWalkNorm * std::exp( -fTimeWalkScale * tw_pe );
-            _cl_sh_prop_delay_set[i][ii] = fPropDelay * dist;
+              double corr = std::pow( dist - fPEAttenuation, 2.0 ) / std::pow( fPEAttenuation, 2.0 );
+              double tw_pe = pe * corr;
 
-          } else { // fill with nonsense
-            _cl_sh_time_walk_set[i][ii]  = -999999.;
-            _cl_sh_prop_delay_set[i][ii] = -999999.;
-          }
+              _cl_sh_time_walk_set[i][ii]  = fTimeWalkNorm * std::exp( -fTimeWalkScale * tw_pe );
+              _cl_sh_prop_delay_set[i][ii] = fPropDelay * dist;
 
+              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]});
+            }
+          else
+            {
+              _cl_sh_time_walk_set[i][ii]  = -999999.;
+              _cl_sh_prop_delay_set[i][ii] = -999999.;
+            }
         }
 
       if(!fDataMode)
@@ -1249,6 +1263,8 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
           _cl_sp_ets0[i]     = spacepoint->Ts0Err();
           _cl_sp_ts1[i]      = spacepoint->Ts1();
           _cl_sp_ets1[i]     = spacepoint->Ts1Err();
+          _cl_sp_dts0[i]     = SPTimeDelta(ts0_set);
+          _cl_sp_dts1[i]     = SPTimeDelta(ts1_set);
           _cl_sp_complete[i] = spacepoint->Complete();
         }
       else
@@ -1265,6 +1281,8 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
           _cl_sp_ets0[i]     = -999999.;
           _cl_sp_ts1[i]      = -999999.;
           _cl_sp_ets1[i]     = -999999.;
+          _cl_sp_dts0[i]     = -999999.;
+          _cl_sp_dts1[i]     = -999999.;
           _cl_sp_complete[i] = false;
         }
     }
@@ -1665,4 +1683,22 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::
   _tpc_tr_score.resize(nActualTracks);
 }
 
+double sbnd::crt::CRTAnalysis::SPTimeDelta(std::vector<std::pair<int, double>> &time_set)
+{
+  std::sort(time_set.begin(), time_set.end(), [](auto &a, auto &b)
+  { return a.first < b.first; });
+
+  double sum = 0;
+
+  for(unsigned int i = 0; i < time_set.size(); ++i)
+    {
+      for(unsigned int ii = i + 1; ii < time_set.size(); ++ii)
+        sum += (time_set[i].second - time_set[ii].second);
+    }
+
+  sum /= (time_set.size() - 1);
+
+  return sum;
+}
+
 DEFINE_ART_MODULE(sbnd::crt::CRTAnalysis)