Merge pull request #641 from SBNSoftware/feature/hlay_crt_channel_status

CRT Channel Status

Author: bear-is-asleep

sbndcode/CRT/CRTReco/CRTStripHitProducer_module.cc

@@ -222,7 +222,7 @@ std::vector<sbnd::crt::CRTStripHit> sbnd::crt::CRTStripHitProducer::CreateStripH
 
       if(unix_diff < -1 || unix_diff > 1)
         {
-          throw std::runtime_error("Unix timestamps differ by more than 1" + unix_diff);
+          throw std::runtime_error(Form("Unix timestamps differ by more than 1 (%li)", unix_diff));
         }
 
       if(unix_diff == 1)
@@ -254,6 +254,9 @@ std::vector<sbnd::crt::CRTStripHit> sbnd::crt::CRTStripHitProducer::CreateStripH
       const CRTSiPMGeo sipm1  = fCRTGeoAlg.GetSiPM(channel);
       const CRTSiPMGeo sipm2  = fCRTGeoAlg.GetSiPM(channel+1);
 
+      if(sipm1.status == CRTChannelStatus::kDeadChannel || sipm2.status == CRTChannelStatus::kDeadChannel)
+        continue;
+
       // Subtract channel pedestals
       const uint16_t adc1 = sipm1.pedestal < sipm_adcs[adc_i]   ? sipm_adcs[adc_i] - sipm1.pedestal   : 0;
       const uint16_t adc2 = sipm2.pedestal < sipm_adcs[adc_i+1] ? sipm_adcs[adc_i+1] - sipm2.pedestal : 0;

sbndcode/Calibration/CRT/CalibService/CMakeLists.txt

@@ -12,4 +12,5 @@ art_make(SERVICE_LIBRARIES
                            art::Persistency_Provenance
                            messagefacility::MF_MessageLogger
                            ROOT::Core
+                           sbnobj::SBND_CRT
 )

sbndcode/Calibration/CRT/CalibService/CRTCalibService.h

@@ -17,6 +17,7 @@
 #include "fhiclcpp/ParameterSet.h"
 #include "art/Framework/Core/ModuleMacros.h"
 #include "art/Framework/Services/Registry/ServiceMacros.h"
+#include "sbnobj/SBND/CRT/CRTEnums.hh"
 
 namespace SBND {
   class CRTCalibService;
@@ -34,12 +35,17 @@ class SBND::CRTCalibService {
 
   double GetPedestalFromFEBMAC5AndChannel(unsigned int feb_mac5, unsigned int ch) const;
 
+  enum sbnd::crt::CRTChannelStatus GetChannelStatusFromFEBMAC5AndChannel(unsigned int feb_mac5,
+                                                                         unsigned int ch) const;
+
 private:
 
   std::unordered_map<unsigned int, double> fTimingOffsetFromFEBMAC5;
 
   std::unordered_map<unsigned int, std::unordered_map<unsigned int, double>> fPedestalFromFEBMAC5AndChannel;
 
+  std::unordered_map<unsigned int, std::unordered_map<unsigned int, sbnd::crt::CRTChannelStatus>> fChannelStatusFromFEBMAC5AndChannel;
+
 };
 
 DECLARE_ART_SERVICE(SBND::CRTCalibService, LEGACY)

sbndcode/Calibration/CRT/CalibService/CRTCalibService_service.cc

@@ -19,11 +19,13 @@ SBND::CRTCalibService::CRTCalibService(fhicl::ParameterSet const& pset)
 {
   const std::string timingOffsetFile = pset.get<std::string>("TimingOffsetFileName");
   const std::string pedestalFile     = pset.get<std::string>("PedestalFileName");
+  const std::string badChannelsFile  = pset.get<std::string>("BadChannelsFileName");
 
-  std::string timingOffsetPath, pedestalPath;
+  std::string timingOffsetPath, pedestalPath, badChannelsPath;
   cet::search_path sp("FW_SEARCH_PATH");
   sp.find_file(timingOffsetFile, timingOffsetPath);
   sp.find_file(pedestalFile, pedestalPath);
+  sp.find_file(badChannelsFile, badChannelsPath);
 
   if(timingOffsetPath.empty())
     {
@@ -37,9 +39,16 @@ SBND::CRTCalibService::CRTCalibService(fhicl::ParameterSet const& pset)
       throw cet::exception("File not found");
     }
 
+  if(badChannelsPath.empty())
+    {
+      std::cout << "SBND::CRTCalibService Input file " << badChannelsFile << " not found" << std::endl;
+      throw cet::exception("File not found");
+    }
+
   std::cout << "SBND CRT Channel Map: Building map from files...\n"
             << "\t Timing Offsets: " << timingOffsetFile << '\n'
-            << "\t Pedestals: " << pedestalFile << std::endl;
+            << "\t Pedestals: " << pedestalFile << '\n'
+            << "\t Bad Channels: " << badChannelsFile << std::endl;
 
   std::ifstream timingOffsetStream(timingOffsetPath, std::ios::in);
   std::string line;
@@ -72,8 +81,49 @@ SBND::CRTCalibService::CRTCalibService(fhicl::ParameterSet const& pset)
 
       fPedestalFromFEBMAC5AndChannel[mac5][ch] = pedestal;
     }
-  
+
   pedestalStream.close();
+
+  std::ifstream badChannelsStream(badChannelsPath, std::ios::in);
+
+  while(std::getline(badChannelsStream, line))
+    {
+      std::stringstream linestream(line);
+
+      unsigned int mac5, ch;
+      std::string status_string;
+
+      linestream
+        >> mac5
+        >> ch
+        >> status_string;
+
+      sbnd::crt::CRTChannelStatus status      = sbnd::crt::CRTChannelStatus::kGoodChannel;
+      sbnd::crt::CRTChannelStatus status_pair = sbnd::crt::CRTChannelStatus::kGoodChannel;
+
+      if(status_string.compare("kDeadChannel") == 0)
+        {
+          status      = sbnd::crt::CRTChannelStatus::kDeadChannel;
+          status_pair = sbnd::crt::CRTChannelStatus::kDeadNeighbourChannel;
+        }
+      else if(status_string.compare("kQuietChannel") == 0)
+        {
+          status      = sbnd::crt::CRTChannelStatus::kQuietChannel;
+          status_pair = sbnd::crt::CRTChannelStatus::kQuietNeighbourChannel;
+        }
+      else
+        {
+          std::cout << "SBND::CRTCalibService unknown channel status " << status_string << std::endl;
+          throw cet::exception("Unknown status");
+        }
+
+      fChannelStatusFromFEBMAC5AndChannel[mac5][ch] = status;
+
+      unsigned int ch_pair = ch % 2 ? ch - 1 : ch + 1;
+      fChannelStatusFromFEBMAC5AndChannel[mac5][ch_pair] = status_pair;
+    }
+
+  badChannelsStream.close();
 }
 
 SBND::CRTCalibService::CRTCalibService(fhicl::ParameterSet const& pset, art::ActivityRegistry&)
@@ -125,4 +175,20 @@ double SBND::CRTCalibService::GetPedestalFromFEBMAC5AndChannel(unsigned int feb_
   return subIter->second;
 }
 
+enum sbnd::crt::CRTChannelStatus SBND::CRTCalibService::GetChannelStatusFromFEBMAC5AndChannel(unsigned int feb_mac5,
+                                                                                              unsigned int ch) const
+{
+  auto iter = fChannelStatusFromFEBMAC5AndChannel.find(feb_mac5);
+
+  if(iter == fChannelStatusFromFEBMAC5AndChannel.end())
+    return sbnd::crt::CRTChannelStatus::kGoodChannel;
+
+  auto subIter = iter->second.find(ch);
+
+  if(subIter == iter->second.end())
+    return sbnd::crt::CRTChannelStatus::kGoodChannel;
+
+  return subIter->second;
+}
+
 DEFINE_ART_SERVICE(SBND::CRTCalibService)

sbndcode/Calibration/CRT/CalibService/SBNDCRTBadChannels_v1.txt

@@ -0,0 +1,14 @@
+182	30	kDeadChannel
+136	0	kDeadChannel
+160	29	kDeadChannel
+160	31	kDeadChannel
+55	8	kDeadChannel
+130	3	kDeadChannel
+82	27	kDeadChannel
+88	27	kDeadChannel
+131	18	kDeadChannel
+117	12	kDeadChannel
+109	11	kDeadChannel
+111	31	kDeadChannel
+77	2	kQuietChannel
+216	11	kQuietChannel

sbndcode/Calibration/CRT/CalibService/crt_calib_service.fcl

@@ -4,6 +4,7 @@ crt_calib_service:
 {
   TimingOffsetFileName: "SBNDCRTTimingOffsets_v1.txt"
   PedestalFileName:     "SBNDCRTPedestals_v2.txt"
+  BadChannelsFileName:  "SBNDCRTBadChannels_v1.txt"
 }
 
 END_PROLOG

sbndcode/Geometry/GeometryWrappers/CRTGeoAlg.cc

@@ -130,22 +130,26 @@ namespace sbnd::crt {
             const uint32_t actualChannel0 = invert ? 31 - (2 * ads_i) : 2 * ads_i;
             const uint32_t actualChannel1 = invert ? actualChannel0 - 1 : actualChannel0 + 1;
 
-            uint32_t pedestal0 = 0;
-            uint32_t pedestal1 = 0;
+            uint32_t pedestal0       = 0;
+            uint32_t pedestal1       = 0;
+            CRTChannelStatus status0 = CRTChannelStatus::kGoodChannel;
+            CRTChannelStatus status1 = CRTChannelStatus::kGoodChannel;
 
             if(!fMC)
               {
                 art::ServiceHandle<SBND::CRTCalibService> CalibService;
                 pedestal0 = CalibService->GetPedestalFromFEBMAC5AndChannel(mac5, actualChannel0);
                 pedestal1 = CalibService->GetPedestalFromFEBMAC5AndChannel(mac5, actualChannel1);
+                status0   = CalibService->GetChannelStatusFromFEBMAC5AndChannel(mac5, actualChannel0);
+                status1   = CalibService->GetChannelStatusFromFEBMAC5AndChannel(mac5, actualChannel1);
               }
 
             const double gain0 = fDefaultGain;
             const double gain1 = fDefaultGain;
 
             // Fill SiPM information
-            CRTSiPMGeo sipm0 = CRTSiPMGeo(stripName, channel0, sipm0XYZWorld, pedestal0, gain0);
-            CRTSiPMGeo sipm1 = CRTSiPMGeo(stripName, channel1, sipm1XYZWorld, pedestal1, gain1);
+            CRTSiPMGeo sipm0 = CRTSiPMGeo(stripName, channel0, sipm0XYZWorld, pedestal0, gain0, status0);
+            CRTSiPMGeo sipm1 = CRTSiPMGeo(stripName, channel1, sipm1XYZWorld, pedestal1, gain1, status1);
             fSiPMs.insert(std::pair<uint16_t, CRTSiPMGeo>(channel0, sipm0));
             fSiPMs.insert(std::pair<uint16_t, CRTSiPMGeo>(channel1, sipm1));
           }

sbndcode/Geometry/GeometryWrappers/CRTGeoAlg.h

@@ -45,7 +45,7 @@ namespace sbnd::crt {
 
   struct CRTSiPMGeo{
     CRTSiPMGeo(const std::string &_stripName, const uint32_t _channel, const geo::Point_t location,
-               const uint32_t _pedestal, const double _gain)
+               const uint32_t _pedestal, const double _gain, const CRTChannelStatus _status)
     {
       stripName = _stripName;
       channel   = _channel;
@@ -54,16 +54,18 @@ namespace sbnd::crt {
       z         = location.Z();
       pedestal  = _pedestal;
       gain      = _gain;
+      status    = _status;
       null      = false;
     }
-    std::string stripName;
-    uint16_t    channel;
-    double      x;
-    double      y;
-    double      z;
-    bool        null;
-    uint32_t    pedestal;
-    double      gain;
+    std::string      stripName;
+    uint16_t         channel;
+    double           x;
+    double           y;
+    double           z;
+    bool             null;
+    uint32_t         pedestal;
+    CRTChannelStatus status;
+    double           gain;
   };
 
   // CRT strip geometry struct contains dimensions and mother module