Merge branch 'develop' into feature/gputnam-TPCChannelInfo

Author: bear-is-asleep

CMakeLists.txt

@@ -15,12 +15,12 @@
 
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
-set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 09.93.01.01)
+set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.04.01)
 find_package(cetmodules REQUIRED)
 project(sbndcode LANGUAGES CXX)
 
 # for CI:
-# project(sbndcode VERSION 09.93.01.01)
+# project(sbndcode VERSION 10.04.01)
 
 message(STATUS "\n")
 message(STATUS "=================================   ${PROJECT_NAME}   =================================")
@@ -80,6 +80,7 @@ find_package( larcore REQUIRED )
 find_package( larcorealg REQUIRED )
 find_package( larcoreobj REQUIRED )
 find_package( larreco REQUIRED )
+find_package( larrecodnn REQUIRED )
 find_package( larana REQUIRED )
 find_package( PostgreSQL REQUIRED )
 find_package( hep_concurrency REQUIRED )

sbndcode/AnalysisTree/AnalysisTree_module.cc

@@ -3677,17 +3677,15 @@ bool sbnd::AnalysisTree::TrackIdToMCTruth( Int_t const trkID, art::Ptr<simb::MCT
 double sbnd::AnalysisTree::bdist(const recob::Track::Point_t& pos)
 {
   // Get geometry.
-  art::ServiceHandle<geo::Geometry> geom;
-
+  geo::TPCGeo const& tpc = art::ServiceHandle<geo::Geometry>()->TPC({0, 0});
   double d1 = pos.X();                             // Distance to right side (wires).
-  double d2 = 2.*geom->DetHalfWidth() - pos.X();   // Distance to left side (cathode).
-  double d3 = pos.Y() + geom->DetHalfHeight();     // Distance to bottom.
-  double d4 = geom->DetHalfHeight() - pos.Y();     // Distance to top.
+  double d2 = 2.*tpc.HalfWidth() - pos.X();   // Distance to left side (cathode).
+  double d3 = pos.Y() + tpc.HalfHeight();     // Distance to bottom.
+  double d4 = tpc.HalfHeight() - pos.Y();     // Distance to top.
   double d5 = pos.Z();                             // Distance to front.
-  double d6 = geom->DetLength() - pos.Z();           // Distance to back.
+  double d6 = tpc.Length() - pos.Z();           // Distance to back.
 
-  double result = std::min(std::min(std::min(std::min(std::min(d1, d2), d3), d4), d5), d6);
-  return result;
+  return std::min({d1, d2, d3, d4, d5, d6});
 }
 
 
@@ -3701,16 +3699,14 @@ double sbnd::AnalysisTree::length(const recob::Track& track)
 // Length of MC particle, trajectory by trajectory.
 double sbnd::AnalysisTree::length(const simb::MCParticle& part, TVector3& start, TVector3& end)
 {
-  // Get geometry.
-  art::ServiceHandle<geo::Geometry> geom;
-
   // Get active volume boundary.
-  double xmin = -2.0 * geom->DetHalfWidth() - 1e-8;
-  double xmax = 2.0 * geom->DetHalfWidth() + 1e-8;
-  double ymin = -geom->DetHalfHeight() -1e-8;
-  double ymax = geom->DetHalfHeight() + 1e-8;
+  geo::TPCGeo const& tpc = art::ServiceHandle<geo::Geometry>()->TPC({0, 0});
+  double xmin = -2.0 * tpc.HalfWidth() - 1e-8;
+  double xmax = 2.0 * tpc.HalfWidth() + 1e-8;
+  double ymin = -tpc.HalfHeight() -1e-8;
+  double ymax = tpc.HalfHeight() + 1e-8;
   double zmin = 0. -1e-8;
-  double zmax = geom->DetLength() + 1e-8;
+  double zmax = tpc.Length() + 1e-8;
 
   // Get number traj points
   int n = part.NumberTrajectoryPoints();

sbndcode/AnalysisTree/CMakeLists.txt

@@ -13,7 +13,6 @@ set(      MODULE_LIBRARIES larcorealg::Geometry
                            lardata::RecoObjects
                            larpandora::LArPandoraInterface
                            sbndcode_CRTUtils
-                           sbndcode_CRT
                            sbnobj::Common_CRT
                            nusimdata::SimulationBase
                            art::Framework_Core

sbndcode/CMakeLists.txt

@@ -21,6 +21,8 @@ add_subdirectory(OpT0Finder)
 add_subdirectory(ShowerAna)
 add_subdirectory(Commissioning)
 add_subdirectory(Trigger)
+add_subdirectory(Timing)
+add_subdirectory(PTBAna)
 
 # subdirectory for gallery stuff
 add_subdirectory(gallery)
@@ -40,3 +42,7 @@ add_subdirectory(WireCell)
 
 add_subdirectory(Decoders)
 add_subdirectory(ChannelMaps)
+
+# for CVN
+# Disable for now
+add_subdirectory(SBNDCVN)

sbndcode/CRT/CMakeLists.txt

@@ -6,48 +6,37 @@ add_subdirectory(CRTSimulation)
 add_subdirectory(CRTTPCMatching)
 add_subdirectory(CRTUtils)
 
-art_make(
-         EXCLUDE
-                CRTChannelMapAlg.cxx
-                CRTGeometryHelper_service.cc
-         LIBRARY_NAME
-                sbndcode_CRTData
-         LIB_LIBRARIES
-                larcorealg::Geometry
-                larcore::Geometry_Geometry_service
-                lardataobj::RecoBase
-                art::Framework_Core
-                art::Framework_IO_Sources
-                art::Framework_Principal
-                canvas::canvas
-                art::Framework_Services_Registry
-                messagefacility::MF_MessageLogger
-
-                fhiclcpp::fhiclcpp
-                ROOT::Core
-                ROOT::Geom
-                cetlib::cetlib
-)
-
 art_make_library(
-                 SOURCE       CRTChannelMapAlg.cxx
-                 LIBRARIES
-                        larcorealg::Geometry
-                        sbndcode_CRTData
-                        cetlib_except::cetlib_except
-                        messagefacility::MF_MessageLogger
-                        fhiclcpp::fhiclcpp
-                        ROOT::Core
-                        ROOT::Geom
-                        cetlib::cetlib
+    SOURCE
+        CRTGeoObjectSorter.cxx
+    LIBRARY_NAME
+        sbndcode_CRTData
+    LIBRARIES
+        larcorealg::Geometry
+        larcore::Geometry_Geometry_service
+        lardataobj::RecoBase
+        art::Framework_Core
+        art::Framework_IO_Sources
+        art::Framework_Principal
+        canvas::canvas
+        art::Framework_Services_Registry
+        messagefacility::MF_MessageLogger
+        fhiclcpp::fhiclcpp
+        ROOT::Core
+        ROOT::Geom
+        cetlib::cetlib
 )
-simple_plugin(
-              CRTGeometryHelper service
-              sbndcode_CRT
-              art::Framework_Core
-              art::Framework_IO_Sources
-              art::Framework_Principal
-              canvas::canvas
-              art::Framework_Services_Registry
-              messagefacility::MF_MessageLogger
+
+cet_build_plugin( CRTAuxDetInitializer art::tool
+    LIBRARIES REG
+        larcorealg::Geometry
+        sbndcode_CRTData
+        cetlib_except::cetlib_except
+        messagefacility::MF_MessageLogger
+        fhiclcpp::fhiclcpp
+        cetlib::cetlib
 )
+
+install_headers()
+install_fhicl()
+install_source()

sbndcode/CRT/CRTAna/CRTAnalysis_module.cc

@@ -36,9 +36,11 @@
 
 #include "sbndcode/Geometry/GeometryWrappers/CRTGeoAlg.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/Decoders/PTB/sbndptb.h"
+#include "sbndcode/Timing/SBNDRawTimingObj.h"
 
 namespace sbnd::crt {
   class CRTAnalysis;
@@ -94,8 +96,10 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::string fMCParticleModuleLabel, fSimDepositModuleLabel, fFEBDataModuleLabel, fCRTStripHitModuleLabel,
     fCRTClusterModuleLabel, fCRTSpacePointModuleLabel, fCRTTrackModuleLabel, fTPCTrackModuleLabel,
     fCRTSpacePointMatchingModuleLabel, fCRTTrackMatchingModuleLabel, fPFPModuleLabel, fPTBModuleLabel,
-    fTDCModuleLabel;
+    fTDCModuleLabel, fTimingReferenceModuleLabel;
   bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC;
+  //! Adding some of the reco parameters to save corrections
+  double fPEAttenuation, fTimeWalkNorm, fTimeWalkScale, fPropDelay;
 
   TTree* fTree;
 
@@ -104,6 +108,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   int _run;
   int _subrun;
   int _event;
+  int _crt_timing_reference_type;
+  int _crt_timing_reference_channel;
 
   //mc truth
   std::vector<int16_t>              _mc_trackid;
@@ -178,6 +184,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<uint32_t>> _cl_sh_ts0_set; //! To store t0 from x-y coincidences
+  std::vector<std::vector<uint32_t>> _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<int>                   _cl_truth_trackid;
   std::vector<double>                _cl_truth_completeness;
   std::vector<double>                _cl_truth_purity;
@@ -331,11 +342,17 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   fPFPModuleLabel                   = p.get<std::string>("PFPModuleLabel", "pandora");
   fPTBModuleLabel                   = p.get<std::string>("PTBModuleLabel", "ptbdecoder");
   fTDCModuleLabel                   = p.get<std::string>("TDCModuleLabel", "tdcdecoder");
+  fTimingReferenceModuleLabel       = p.get<std::string>("TimingReferenceModuleLabel", "crtstrips");
   fDebug                            = p.get<bool>("Debug", false);
   fDataMode                         = p.get<bool>("DataMode", false);
   fNoTPC                            = p.get<bool>("NoTPC", false);
   fHasPTB                           = p.get<bool>("HasPTB", false);
   fHasTDC                           = p.get<bool>("HasTDC", false);
+  //! Adding some of the reco parameters to save corrections
+  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)
     fCRTBackTrackerAlg = CRTBackTrackerAlg(p.get<fhicl::ParameterSet>("CRTBackTrackerAlg", fhicl::ParameterSet()));
@@ -346,6 +363,8 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   fTree->Branch("run", &_run);
   fTree->Branch("subrun", &_subrun);
   fTree->Branch("event", &_event);
+  fTree->Branch("crt_timing_reference_type", &_crt_timing_reference_type);
+  fTree->Branch("crt_timing_reference_channel", &_crt_timing_reference_channel);
 
   if(!fDataMode)
     {
@@ -421,6 +440,11 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   fTree->Branch("cl_composition", "std::vector<uint8_t>", &_cl_composition);
   fTree->Branch("cl_channel_set", "std::vector<std::vector<uint32_t>>", &_cl_channel_set);
   fTree->Branch("cl_adc_set", "std::vector<std::vector<uint16_t>>", &_cl_adc_set);
+  fTree->Branch("cl_sh_ts0_set", "std::vector<std::vector<uint32_t>>", &_cl_sh_ts0_set);
+  fTree->Branch("cl_sh_ts1_set", "std::vector<std::vector<uint32_t>>", &_cl_sh_ts1_set);
+  fTree->Branch("cl_sh_feb_mac5_set", "std::vector<std::vector<uint16_t>>", &_cl_sh_feb_mac5_set);
+  fTree->Branch("cl_sh_time_walk_set", "std::vector<std::vector<double>>", &_cl_sh_time_walk_set);
+  fTree->Branch("cl_sh_prop_delay_set", "std::vector<std::vector<double>>", &_cl_sh_prop_delay_set);
   if(!fDataMode)
     {
       fTree->Branch("cl_truth_trackid", "std::vector<int>", &_cl_truth_trackid);
@@ -621,6 +645,17 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e)
 
   if(fDebug) std::cout << "This is event " << _run << "-" << _subrun << "-" << _event << std::endl;
 
+  _crt_timing_reference_type    = -1;
+  _crt_timing_reference_channel = -1;
+
+  art::Handle<raw::TimingReferenceInfo> TimingReferenceHandle;
+  e.getByLabel(fTimingReferenceModuleLabel, TimingReferenceHandle);
+  if(TimingReferenceHandle.isValid())
+    {
+      _crt_timing_reference_type    = TimingReferenceHandle->timingType;
+      _crt_timing_reference_channel = TimingReferenceHandle->timingChannel;
+    }
+
   if(fHasPTB)
     {
       // Get PTBs
@@ -1076,6 +1111,11 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
   _cl_composition.resize(nClusters);
   _cl_channel_set.resize(nClusters);
   _cl_adc_set.resize(nClusters);
+  _cl_sh_ts0_set.resize(nClusters);
+  _cl_sh_ts1_set.resize(nClusters);
+  _cl_sh_feb_mac5_set.resize(nClusters);
+  _cl_sh_time_walk_set.resize(nClusters);
+  _cl_sh_prop_delay_set.resize(nClusters);
   _cl_truth_trackid.resize(nClusters);
   _cl_truth_completeness.resize(nClusters);
   _cl_truth_purity.resize(nClusters);
@@ -1106,6 +1146,8 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
   _cl_sp_ets1.resize(nClusters);
   _cl_sp_complete.resize(nClusters);
 
+  art::ServiceHandle<SBND::CRTChannelMapService> ChannelMapService;
+
   for(unsigned i = 0; i < nClusters; ++i)
     {
       const auto cluster = CRTClusterVec[i];
@@ -1117,16 +1159,52 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
       _cl_tagger[i]      = cluster->Tagger();
       _cl_composition[i] = cluster->Composition();
 
+      const auto spacepoints = clustersToSpacePoints.at(cluster.key());
+
       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_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]);
 
       for(unsigned ii = 0; ii < _cl_nhits[i]; ++ii)
         {
           const auto striphit = striphits[ii];
           _cl_channel_set[i][ii] = striphit->Channel();
           _cl_adc_set[i][2*ii]   = striphit->ADC1();
           _cl_adc_set[i][2*ii+1] = striphit->ADC2();
+	  _cl_sh_ts0_set[i][ii]  = striphit->Ts0();
+	  _cl_sh_ts1_set[i][ii]  = striphit->Ts1();
+	  SBND::CRTChannelMapService::ModuleInfo_t module_info = ChannelMapService->GetModuleInfoFromOfflineID( striphit->Channel() / 32 );
+	  _cl_sh_feb_mac5_set[i][ii] = ( module_info.valid ) ? module_info.feb_mac5 : 0;
+
+	  /*
+	   * 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) { // 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() );
+
+	    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;
+
+	  } else { // fill with nonsense
+	    _cl_sh_time_walk_set[i][ii]  = -999999.;
+	    _cl_sh_prop_delay_set[i][ii] = -999999.;
+	  }
+
         }
 
       if(!fDataMode)
@@ -1150,7 +1228,6 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
           _cl_truth_core_z[i]           = truthMatch.deposit.coreZ;
         }
 
-      const auto spacepoints = clustersToSpacePoints.at(cluster.key());
       if(spacepoints.size() == 1)
         { 
           const auto spacepoint = spacepoints[0];

sbndcode/CRT/CRTAna/run_crtana.fcl

@@ -1,5 +1,6 @@
 #include "crtgeoalg_sbnd.fcl"
 #include "crtbacktrackeralg_sbnd.fcl"
+#include "crtrecoproducers_sbnd.fcl"
 #include "simulationservices_sbnd.fcl"
 
 process_name: CRTAnalyzer
@@ -27,6 +28,10 @@ physics:
       module_type:       "CRTAnalysis"
       CRTGeoAlg:         @local::crtgeoalg_sbnd
       CRTBackTrackerAlg: @local::crtbacktrackeralg_sbnd
+      PEAttenuation:     @local::sbnd_crtsim.DetSimParams.NpeScaleShift	
+      PropDelay:         @local::sbnd_crtsim.DetSimParams.PropDelay	
+      TimeWalkNorm:      @local::sbnd_crtsim.DetSimParams.TDelayNorm	
+      TimeWalkScale:     @local::sbnd_crtsim.DetSimParams.TDelayScale
     }
   }