Merge branch 'release/v10_04_01_01'

Author: Brinden Carlson

CMakeLists.txt

@@ -15,12 +15,12 @@
 
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
-set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.04.01)
+set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.04.01.01)
 find_package(cetmodules REQUIRED)
 project(sbndcode LANGUAGES CXX)
 
 # for CI:
-# project(sbndcode VERSION 10.04.01)
+# project(sbndcode VERSION 10.04.01.01)
 
 message(STATUS "\n")
 message(STATUS "=================================   ${PROJECT_NAME}   =================================")

sbndcode/CMakeLists.txt

@@ -44,5 +44,4 @@ add_subdirectory(Decoders)
 add_subdirectory(ChannelMaps)
 
 # for CVN
-# Disable for now
-#add_subdirectory(SBNDCVN)
+add_subdirectory(SBNDCVN)

sbndcode/CRT/CRTAna/CRTAnalysis_module.cc

@@ -36,6 +36,7 @@
 
 #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"
@@ -97,6 +98,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
     fCRTSpacePointMatchingModuleLabel, fCRTTrackMatchingModuleLabel, fPFPModuleLabel, fPTBModuleLabel,
     fTDCModuleLabel, fTimingReferenceModuleLabel;
   bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC;
+  //! Adding some of the reco parameters to save corrections
+  double fPEAttenuation, fTimeWalkNorm, fTimeWalkScale, fPropDelay;
 
   TTree* fTree;
 
@@ -181,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;
@@ -340,6 +348,11 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   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()));
@@ -427,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);
@@ -1093,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);
@@ -1123,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];
@@ -1134,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)
@@ -1167,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
     }
   }
 

sbndcode/Commissioning/CMakeLists.txt

@@ -20,6 +20,7 @@ set(
         sbndcode_CRTUtils
         sbnobj::Common_CRT
         sbnobj::SBND_CRT
+	sbnobj::SBND_Timing
         sbndaq_artdaq_core::sbndaq-artdaq-core_Obj_SBND
         nusimdata::SimulationBase
         art::Framework_Core

sbndcode/Commissioning/HitDumper_module.cc

@@ -58,6 +58,8 @@
 #include "sbndaq-artdaq-core/Obj/SBND/pmtSoftwareTrigger.hh"
 #include "sbndaq-artdaq-core/Obj/SBND/CRTmetric.hh"
 
+#include "sbnobj/SBND/Timing/DAQTimestamp.hh"
+
 // Truth includes
 //#include "larsim/MCCheater/BackTrackerService.h"
 //#include "larsim/MCCheater/ParticleInventoryService.h"
@@ -127,6 +129,9 @@ class Hitdumper : public art::EDAnalyzer {
 
   // Called at the beginning of every subrun
   virtual void beginSubRun(art::SubRun const& sr) override;
+
+  void AnalyseTDCs(std::vector<art::Ptr<sbnd::timing::DAQTimestamp>> &TDCVec);
+
 private:
 
   /// Resets the variables that are saved to the TTree
@@ -162,6 +167,11 @@ class Hitdumper : public art::EDAnalyzer {
   /// Resize the data structure for MCShowers
   void ResizeMCShower(int nShowers);
 
+  std::string fTDCModuleLabel;
+
+  bool fHasTDC;
+
+
   opdet::sbndPDMapAlg _pd_map;
 
   TTree* fTree;
@@ -438,6 +448,12 @@ class Hitdumper : public art::EDAnalyzer {
 
   const int MAX_INT = std::numeric_limits<int>::max();
   const long int TIME_CORRECTION = (long int) std::numeric_limits<int>::max() * 2;
+
+  std::vector<uint32_t>    _tdc_channel;
+  std::vector<uint64_t>    _tdc_timestamp;
+  std::vector<uint64_t>    _tdc_offset;
+  std::vector<std::string> _tdc_name;
+
 };
 
 
@@ -460,6 +476,9 @@ Hitdumper::Hitdumper(fhicl::ParameterSet const& pset)
 void Hitdumper::reconfigure(fhicl::ParameterSet const& p)
 {
 
+  fTDCModuleLabel                   = p.get<std::string>("TDCModuleLabel", "tdcdecoder");
+  fHasTDC                           = p.get<bool>("HasTDC", false);
+
   _max_hits = p.get<int>("MaxHits", 50000);
   _max_ophits = p.get<int>("MaxOpHits", 50000);
   _max_samples = p.get<int>("MaxSamples", 5001);
@@ -523,6 +542,24 @@ void Hitdumper::analyze(const art::Event& evt)
   _t0 = 0.;
   // t0 = detprop->TriggerOffset();  // units of TPC ticks
 
+
+  if(fHasTDC)
+    {
+      // Get TDCs                                                                               
+      art::Handle<std::vector<sbnd::timing::DAQTimestamp>> TDCHandle;
+      evt.getByLabel(fTDCModuleLabel, TDCHandle);
+      if(!TDCHandle.isValid()){
+	std::cout << "TDC product " << fTDCModuleLabel << " not found..." << std::endl;
+        throw std::exception();
+      }
+      std::vector<art::Ptr<sbnd::timing::DAQTimestamp>> TDCVec;
+      art::fill_ptr_vector(TDCVec, TDCHandle);
+
+      // Fill TDC variables                                                                     
+      AnalyseTDCs(TDCVec);
+    }
+
+
   //
   // Hits
   //
@@ -1190,6 +1227,15 @@ void Hitdumper::analyze(const art::Event& evt)
   fTree->Branch("evttime",&_evttime,"evttime/D");
   fTree->Branch("t0",&_t0,"t0/I");
 
+  if(fHasTDC)
+    {
+      fTree->Branch("tdc_channel", "std::vector<uint32_t>", &_tdc_channel);
+      fTree->Branch("tdc_timestamp", "std::vector<uint64_t>", &_tdc_timestamp);
+      fTree->Branch("tdc_offset", "std::vector<uint64_t>", &_tdc_offset);
+      fTree->Branch("tdc_name", "std::vector<std::string>", &_tdc_name);
+    }
+
+
   fTree->Branch("nhits", &_nhits, "nhits/I");
   fTree->Branch("hit_cryostat", &_hit_cryostat);
   fTree->Branch("hit_tpc", &_hit_tpc);
@@ -1410,6 +1456,29 @@ void Hitdumper::analyze(const art::Event& evt)
 }
 
 
+void Hitdumper::AnalyseTDCs(std::vector<art::Ptr<sbnd::timing::DAQTimestamp>> &TDCVec)
+{
+  const unsigned nTDCs = TDCVec.size();
+
+  _tdc_channel.resize(nTDCs);
+  _tdc_timestamp.resize(nTDCs);
+  _tdc_offset.resize(nTDCs);
+  _tdc_name.resize(nTDCs);
+
+  unsigned tdc_i = 0;
+
+  for(auto const& tdc : TDCVec)
+    {
+      _tdc_channel[tdc_i]   = tdc->Channel();
+      _tdc_timestamp[tdc_i] = tdc->Timestamp();
+      _tdc_offset[tdc_i]    = tdc->Offset();
+      _tdc_name[tdc_i]      = tdc->Name();
+
+      ++tdc_i;
+    }
+}
+
+
 void Hitdumper::ResetWireHitsVars(int n) {
   _hit_cryostat.assign(n, DEFAULT_VALUE);
   _hit_tpc.assign(n, DEFAULT_VALUE);

sbndcode/Commissioning/fcls/hitdumpermodule.fcl

@@ -27,6 +27,7 @@ hitdumper:
     MCParticleModuleLabel:    "largeant"
     MCTrackModuleLabel:       "mcreco"
     MCShowerModuleLabel:      "mcreco"
+    TDCModuleLabel:           "tdcdecoder"
 
     KeepCRTStripHits:         true
     KeepCRTSpacePoints:       true
@@ -41,6 +42,7 @@ hitdumper:
     readMCParticle:           false
     savePOTinfo:              true
     checkTransparency:        false
+    HasTDC:                   false
 
     window:                   100
 
@@ -55,5 +57,6 @@ hitdumper:
 hitdumper_data: @local::hitdumper
 hitdumper_data.OpHitsModuleLabel: ["ophitpmt"] #only have PMT data as of 05/30/2024
 hitdumper_data.readTruth: false
+hitdumper_data.HasTDC: true
 
 END_PROLOG