Merge tag 'v10_20_05' into feature/frame_shift_refactor_v10_14_02

v10_20_05

Author: VCLanNguyen

CMakeLists.txt

@@ -16,7 +16,7 @@
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
 find_package(cetmodules 3.20.00 REQUIRED)
-project(sbncode VERSION 10.14.02 LANGUAGES CXX)
+project(sbncode VERSION 10.20.05 LANGUAGES CXX)
 
 message(STATUS "\n\n  ==========================   ${PROJECT_NAME}   ==========================")
 

sbncode/BeamSpillInfoRetriever/POTTools.cpp

@@ -20,11 +20,11 @@ namespace sbn::pot{
             if (ctb_frag.Trigger(word_i)->IsHLT() && ctb_frag.Trigger(word_i)->IsTrigger(HLT))
             {
               foundHLT = true;
-              uint64_t RawprevPTBTimeStamp = ctb_frag.PTBWord(word_i)->prevTS * 20;
-              uint64_t RawcurrPTBTimeStamp = ctb_frag.Trigger(word_i)->timestamp * 20;
-              double currTS_candidate = std::bitset<64>(RawcurrPTBTimeStamp/20).to_ullong()/50e6;
+              uint64_t RawprevPTBTimeStamp = ctb_frag.PTBWord(word_i)->prevTS;
+              uint64_t RawcurrPTBTimeStamp = ctb_frag.Trigger(word_i)->timestamp;
+              std::uint64_t currTS_candidate = std::bitset<64>(RawcurrPTBTimeStamp).to_ullong() * 20;
               if(currTS_candidate < PTBInfo.currPTBTimeStamp){
-                PTBInfo.prevPTBTimeStamp = std::bitset<64>(RawprevPTBTimeStamp / 20).to_ullong()/50e6;
+                PTBInfo.prevPTBTimeStamp = std::bitset<64>(RawprevPTBTimeStamp).to_ullong() * 20;
                 PTBInfo.currPTBTimeStamp = currTS_candidate;
                 PTBInfo.GateCounter = ctb_frag.Trigger(word_i)->gate_counter;
               }
@@ -43,6 +43,35 @@ namespace sbn::pot{
     }
   }
 
+  std::vector<PTBInfo_t> extractAllPTBInfo(std::vector<artdaq::Fragment> const& cont_frags) {
+    std::vector<PTBInfo_t> PTBInfoVec;
+    for (auto const& cont : cont_frags)
+    {
+      artdaq::ContainerFragment cont_frag(cont);
+      for (size_t fragi = 0; fragi < cont_frag.block_count(); ++fragi)
+      {
+        artdaq::Fragment frag = *cont_frag[fragi];
+        sbndaq::CTBFragment ctb_frag(frag);   
+        for(size_t word_i = 0; word_i < ctb_frag.NWords(); ++word_i)
+        {
+          if(ctb_frag.Trigger(word_i)){
+            PTBInfo_t PTBInfo;
+            uint64_t RawprevPTBTimeStamp = ctb_frag.PTBWord(word_i)->prevTS;
+            uint64_t RawcurrPTBTimeStamp = ctb_frag.Trigger(word_i)->timestamp;
+            PTBInfo.prevPTBTimeStamp = std::bitset<64>(RawprevPTBTimeStamp).to_ullong() * 20;
+            PTBInfo.currPTBTimeStamp = std::bitset<64>(RawcurrPTBTimeStamp).to_ullong() * 20;
+            PTBInfo.GateCounter = ctb_frag.Trigger(word_i)->gate_counter;
+            PTBInfo.isHLT = ctb_frag.Trigger(word_i)->IsHLT();
+            PTBInfo.triggerWord = ctb_frag.Trigger(word_i)->trigger_word;
+            PTBInfoVec.push_back(PTBInfo);
+          }
+        } 
+      } 
+    } 
+  
+    return PTBInfoVec;
+  }
+
   double extractTDCTimeStamp(art::Handle<std::vector<artdaq::Fragment> > cont_frags) {
 
     double TDCTimeStamp = 0;

sbncode/BeamSpillInfoRetriever/POTTools.h

@@ -13,18 +13,21 @@
 
 #include "sbncode/BeamSpillInfoRetriever/MWRData.h"
 #include "larcorealg/CoreUtils/counter.h"
+#include <vector>
 
 namespace sbn::pot{
 
   typedef struct PTBInfo_t {
-    double currPTBTimeStamp  = 1e20;
-    double prevPTBTimeStamp  = 0;
+    std::uint64_t currPTBTimeStamp  = UINT64_MAX;  ///< Timestamp in UTC nanoseconds since Unix epoch (converted from 20ns clock ticks)
+    std::uint64_t prevPTBTimeStamp  = 0; 
     unsigned int GateCounter = 0;
+    bool isHLT = false;            
+    uint64_t triggerWord = 0;  ///< Timestamp in s since beam extraction signal
   } PTBInfo_t;
 
   typedef struct TriggerInfo_t {
     int gate_type = 0; ///< Source of the spill: `1`: BNB, `2`: NuMI
-    double t_current_event  = 0;
+    double t_current_event  = 0;  ///< Timestamp in UTC seconds since Unix epoch (converted from 20ns clock ticks)
     double t_previous_event = 0;
     unsigned int number_of_gates_since_previous_event = 0;
     std::int64_t WR_to_Spill_conversion = 0;
@@ -36,13 +39,21 @@ namespace sbn::pot{
   } MWRdata_t;
 
   /**
-   * @brief Extracts information from PTB for use in SBND POT accounting.
+   * @brief Extracts information from PTB for a single HLT for use in SBND POT accounting.
    * 
    * @param cont_frags The PTB fragments to examine.
    * @param HLT The high level trigger we are searching for.
    */
   PTBInfo_t extractPTBInfo(art::Handle<std::vector<artdaq::Fragment> > cont_frags, int HLT);
 
+  /**
+   * @brief Extracts ALL PTB information for using in SBND CAF files.
+   * 
+   * @param cont_frags The PTB fragments to examine.
+   * @return Vector of PTBInfo_t containing all triggers found.
+   */
+  std::vector<PTBInfo_t> extractAllPTBInfo(std::vector<artdaq::Fragment> const& cont_frags);
+
   /**
    * @brief Extracts information from TDC for use in SBND POT accounting.
    * 

sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc

@@ -123,13 +123,13 @@ sbn::pot::TriggerInfo_t sbn::SBNDBNBRetriever::extractTriggerInfo(art::Event con
   else{
     // If missing TDC, use PTB instead
     mf::LogDebug("SBNDBNBRetriever") << " Missing TDC Container Fragments!!! " << std::endl;
-    triggerInfo.t_current_event = PTBInfo.currPTBTimeStamp - fBESOffset;
+    triggerInfo.t_current_event = PTBInfo.currPTBTimeStamp / 1e9 - fBESOffset;
   }
 
-  triggerInfo.t_previous_event = PTBInfo.prevPTBTimeStamp - fBESOffset;
+  triggerInfo.t_previous_event = PTBInfo.prevPTBTimeStamp / 1e9 - fBESOffset;
   triggerInfo.number_of_gates_since_previous_event = PTBInfo.GateCounter;
 
-  double PTBandCurrOffset = PTBInfo.currPTBTimeStamp - triggerInfo.t_current_event - fBESOffset;
+  double PTBandCurrOffset = PTBInfo.currPTBTimeStamp / 1e9 - triggerInfo.t_current_event - fBESOffset;
 
   // Catch for an issue seen a few times where PTB off by a second.
   // Only need to correct prevTS because either currTS is from TDC

sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc

@@ -126,13 +126,13 @@ sbn::pot::TriggerInfo_t sbn::SBNDBNBZEROBIASRetriever::extractTriggerInfo(art::E
   else{
     // If missing TDC, use PTB instead
     mf::LogDebug("SBNDBNBZEROBIASRetriever") << " Missing TDC Container Fragments!!!" << std::endl;
-    triggerInfo.t_current_event = PTBInfo.currPTBTimeStamp - fBESOffset;
+    triggerInfo.t_current_event = PTBInfo.currPTBTimeStamp / 1e9 - fBESOffset;
   }
 
-  triggerInfo.t_previous_event = PTBInfo.prevPTBTimeStamp - fBESOffset;
+  triggerInfo.t_previous_event = PTBInfo.prevPTBTimeStamp / 1e9 - fBESOffset;
   triggerInfo.number_of_gates_since_previous_event = PTBInfo.GateCounter;
 
-  double PTBandCurrOffset = PTBInfo.currPTBTimeStamp - triggerInfo.t_current_event - fBESOffset;
+  double PTBandCurrOffset = PTBInfo.currPTBTimeStamp / 1e9 - triggerInfo.t_current_event - fBESOffset;
 
   // Catch for an issue seen a few times where PTB off by a second.
   // Only need to correct prevTS because either currTS is from TDC

sbncode/CAFMaker/CAFMakerParams.h

@@ -483,6 +483,12 @@ namespace caf
       "simdrift"
     };
 
+    Atom<art::InputTag> SimEnergyDepositLabel {
+      Name("SimEnergyDepositLabel"),
+      Comment("Label of input sim::SimEnergyDeposit objects."),
+      art::InputTag("ionandscint", "priorSCE","G4")
+    };
+
     Atom<bool> FillTrueParticles {
       Name("FillTrueParticles"),
       Comment("Whether to fill the rec.true_particles branch. The information on true particles"
@@ -634,7 +640,18 @@ namespace caf
       Comment("Label of CVN scores."),
       "cvn" 
     };
+
+    Atom<string> LightCaloLabel {
+      Name("LightCaloLabel"),
+      Comment("Label of light calorimetry producer"),
+      "lightcalo"
+    };
 
+    Atom<std::string> fBlipTag { 
+      Name("BlipTag"),
+      Comment("Provide a string to label the blip input"), 
+      "blipreco"
+    };
   };
 }
 

sbncode/CAFMaker/CAFMaker_module.cc

@@ -20,6 +20,7 @@
 #include "sbncode/CAFMaker/FillExposure.h"
 #include "sbncode/CAFMaker/FillTrigger.h"
 #include "sbncode/CAFMaker/Utils.h"
+#include "sbncode/CAFMaker/FillBlip.h"
 
 // C/C++ includes
 #include <fenv.h>
@@ -97,6 +98,7 @@
 #include "lardataobj/RecoBase/MCSFitResult.h"
 #include "lardataobj/RecoBase/Cluster.h"
 #include "lardataobj/AnalysisBase/MVAOutput.h"
+#include "lardataobj/Simulation/SimEnergyDeposit.h"
 
 #include "nusimdata/SimulationBase/MCFlux.h"
 #include "nusimdata/SimulationBase/MCTruth.h"
@@ -115,11 +117,14 @@
 #include "sbnobj/Common/POTAccounting/BNBSpillInfo.h"
 #include "sbnobj/Common/POTAccounting/EXTCountInfo.h"
 #include "sbnobj/Common/POTAccounting/NuMISpillInfo.h"
+#include "sbncode/BeamSpillInfoRetriever/POTTools.h"
+#include "artdaq-core/Data/ContainerFragment.hh"
 #include "sbnobj/Common/Trigger/ExtraTriggerInfo.h"
 #include "sbnobj/Common/Reco/CRUMBSResult.h"
 #include "sbnobj/Common/Reco/OpT0FinderResult.h"
 #include "sbnobj/SBND/CRT/CRTVeto.hh"
 #include "sbnobj/Common/Reco/CorrectedOpFlashTiming.h"
+#include "sbnobj/Common/Reco/LightCalo.h"
 #include "sbnobj/SBND/Timing/TimingInfo.hh"
 #include "sbnobj/SBND/Timing/FrameShiftInfo.hh"
 
@@ -352,6 +357,10 @@ class CAFMaker : public art::EDProducer {
   art::FindOneP<T> FindOnePStrict(const U& from, const art::Event& evt,
 				  const art::InputTag& label) const;
 
+  template <class T, class U>
+  art::FindOneP<T> FindOnePStrictSingle(const U& from, const art::Event& evt,
+                                        const art::InputTag& label) const;
+
   template <class T, class D, class U>
   art::FindOneP<T, D> FindOnePDStrict(const U& from,
                                       const art::Event& evt,
@@ -1299,6 +1308,14 @@ art::FindOneP<T> CAFMaker::FindOnePStrict(const U& from,
   return ret;
 }
 
+//......................................................................
+template <class T, class U>
+art::FindOneP<T> CAFMaker::FindOnePStrictSingle(const U& from,
+                                                const art::Event& evt,
+                                                const art::InputTag& tag) const {
+  return FindOnePStrict<T>(std::vector{ from }, evt, tag);
+}
+
 //......................................................................
 template <class T, class D, class U>
 art::FindOneP<T, D> CAFMaker::FindOnePDStrict(const U& from,
@@ -1378,6 +1395,7 @@ bool CAFMaker::GetPsetParameter(const fhicl::ParameterSet& pset,
 
 //......................................................................
 void CAFMaker::produce(art::Event& evt) noexcept {
+  mf::LogInfo("CAFMaker") << "CAFMaker::produce called for event: " << evt.id();
 
   bool const firstInFile = (fIndexInFile++ == 0);
 
@@ -1637,6 +1655,36 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     } // end for fm
   } // end for i (mctruths)
 
+  // get sim energy deposits if they're there
+  ::art::Handle<std::vector<sim::SimEnergyDeposit>> sed_handle;
+  GetByLabelStrict(evt, fParams.SimEnergyDepositLabel().encode(), sed_handle);
+
+  if (!isRealData && sed_handle.isValid()){
+    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
+
+    srtruthbranch.dep.reserve(mctruths.size());
+    for (size_t n=0; n<mctruths.size();n++){
+      SRTrueDeposit init;
+      init.electrons = 0;
+      init.photons = 0;
+      init.energy = 0; 
+      srtruthbranch.dep.push_back(init);
+    }
+
+    for (sim::SimEnergyDeposit const& sed: *sed_handle){
+      const auto trackID = sed.TrackID();
+
+      art::Ptr<simb::MCTruth> mctruth = pi_serv->TrackIdToMCTruth_P(trackID);
+      auto it = std::find(mctruths.begin(), mctruths.end(), mctruth);
+      if (it == mctruths.end()) continue; 
+
+      auto idx = std::distance(mctruths.begin(), it);
+      srtruthbranch.dep.at(idx).energy    += sed.Energy()*1e-3; // GeV
+      srtruthbranch.dep.at(idx).photons   += sed.NumPhotons();
+      srtruthbranch.dep.at(idx).electrons += sed.NumElectrons();
+    }
+  }
+
   // get the number of events generated in the gen stage
   unsigned n_gen_evt = 0;
   for (const art::ProcessConfiguration &process: evt.processHistory()) {
@@ -1711,6 +1759,24 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       FillTriggerEmulation(monpulses_handle, monpulse_sizes_handle, pairs_handle, trigemu_handle, srtrigger);
   }
 
+  
+  // Fill PTB (Penn Trigger Board) information for SBND real data
+  if (isRealData && fDet == kSBND) {
+    art::InputTag PTB_itag("daq", "ContainerPTB");
+    art::Handle<artdaq::Fragments> ptb_frags_handle;
+    evt.getByLabel(PTB_itag, ptb_frags_handle);
+    if (ptb_frags_handle.isValid()) {
+      mf::LogDebug("CAFMaker") << "Found ContainerPTB, extracting PTB triggers...";
+      std::vector<sbn::pot::PTBInfo_t> ptb_triggers = sbn::pot::extractAllPTBInfo(*ptb_frags_handle);
+      mf::LogDebug("CAFMaker") << "Extracted " << ptb_triggers.size() << " PTB triggers";
+      FillPTBTriggersSBND(ptb_triggers, srtrigger);
+      mf::LogDebug("CAFMaker") << "PTB HLT triggers: " << srtrigger.ptb_hlt_timestamp.size() 
+                                << ", LLT triggers: " << srtrigger.ptb_llt_timestamp.size();
+    } else {
+      mf::LogDebug("CAFMaker") << "ContainerPTB not found for event " << evt.id();
+    }
+  }
+
   // If not real data, fill in enough of the SRTrigger to make (e.g.) the CRT
   // time referencing work. TODO: add more stuff to a "MC"-Trigger?
   // No longer needed with incorporation of trigger emulation in the MC.
@@ -1791,12 +1857,15 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     {
       art::Handle<std::vector<sbnd::crt::CRTSpacePoint>> crtspacepoints_handle;
       GetByLabelStrict(evt, fParams.CRTSpacePointLabel(), crtspacepoints_handle);
+      art::FindOneP<sbnd::crt::CRTCluster> foCRTCluster =
+        FindOnePStrict<sbnd::crt::CRTCluster>(crtspacepoints_handle, evt, fParams.CRTSpacePointLabel());
 
       if (crtspacepoints_handle.isValid()) {
         const std::vector<sbnd::crt::CRTSpacePoint> &crtspacepoints = *crtspacepoints_handle;
         for (unsigned i = 0; i < crtspacepoints.size(); i++) {
           srcrtspacepoints.emplace_back();
-          FillCRTSpacePoint(crtspacepoints[i], srcrtspacepoints.back());
+          const art::Ptr<sbnd::crt::CRTCluster> crtcluster = foCRTCluster.at(i);
+          FillCRTSpacePoint(crtspacepoints[i], *crtcluster, srcrtspacepoints.back());
         }
       }
 
@@ -1928,6 +1997,13 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     }
   }
 
+  //Fill blips. art::handle for blips and then call fill blips for each one. Make a vector to hold all of them. I handle for loop in Fill blip
+  art::Handle<std::vector<blip::Blip>> blipHandle;
+  std::vector<caf::SRBlip> srblips;
+  if(evt.getByLabel( fParams.fBlipTag(), blipHandle)) //fill SR blips
+  {
+    FillBlip( *blipHandle, srblips);
+  }
   // collect the TPC slices
   std::vector<art::Ptr<recob::Slice>> slices;
   std::vector<std::string> slice_tag_suffixes;
@@ -2049,6 +2125,10 @@ void CAFMaker::produce(art::Event& evt) noexcept {
      if (fmCorrectedOpFlash.isValid())
       slcCorrectedOpFlash = fmCorrectedOpFlash.at(0);
 
+    art::FindOneP<sbn::LightCalo> foLightCalo =
+      FindOnePStrict<sbn::LightCalo>(sliceList,evt,
+        fParams.LightCaloLabel() + slice_tag_suff);
+    const sbn::LightCalo *slcLightCalo = foLightCalo.isValid()? foLightCalo.at(0).get() : nullptr;
 
     art::FindOneP<lcvn::Result> foCVNResult =
       FindOnePStrict<lcvn::Result>(sliceList, evt,
@@ -2239,6 +2319,9 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       FindOnePDStrict<sbnd::crt::CRTSpacePoint, anab::T0>(slcTracks, evt,
                fParams.CRTSpacePointMatchLabel() + slice_tag_suff);
 
+    art::Handle<std::vector<sbnd::crt::CRTSpacePoint>> crtspacepoints_handle;
+    GetByLabelStrict(evt, fParams.CRTSpacePointLabel(), crtspacepoints_handle);
+
     art::FindOneP<sbnd::crt::CRTTrack, anab::T0> foSBNDCRTTrackMatch =
       FindOnePDStrict<sbnd::crt::CRTTrack, anab::T0>(slcTracks, evt,
                fParams.SBNDCRTTrackMatchLabel() + slice_tag_suff);
@@ -2308,6 +2391,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     FillSliceCRUMBS(slcCRUMBS, recslc);
     FillSliceOpT0Finder(slcOpT0, recslc);
     FillSliceBarycenter(slcHits, slcSpacePoints, recslc);
+    FillSliceLightCalo(slcLightCalo, recslc);
     FillTPCPMTBarycenterMatch(barycenterMatch, recslc);
     FillCorrectedOpFlashTiming(slcCorrectedOpFlash, recslc);
     FillCVNScores(cvnResult, recslc);
@@ -2515,8 +2599,12 @@ void CAFMaker::produce(art::Event& evt) noexcept {
           {
             const art::Ptr<sbnd::crt::CRTSpacePoint> crtspacepoint = foCRTSpacePointMatch.at(iPart);
 
+            art::FindOneP<sbnd::crt::CRTCluster> foCRTCluster =
+              FindOnePStrictSingle<sbnd::crt::CRTCluster>(crtspacepoint, evt, fParams.CRTSpacePointLabel());
+            const art::Ptr<sbnd::crt::CRTCluster>& crtcluster = foCRTCluster.at(0);
+
             if(crtspacepoint.isNonnull())
-              FillTrackCRTSpacePoint(foCRTSpacePointMatch.data(iPart).ref(), crtspacepoint, trk);
+              FillTrackCRTSpacePoint(foCRTSpacePointMatch.data(iPart).ref(), *crtspacepoint, *crtcluster, trk);
           }
         if(foSBNDCRTTrackMatch.isValid() && fDet == kSBND)
           {
@@ -2605,6 +2693,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   rec.sbnd_crt_veto    = srsbndcrtveto;
   rec.opflashes        = srflashes;
   rec.nopflashes       = srflashes.size();
+  rec.blips            = srblips;
   rec.sbnd_frames      = srsbndframeshiftinfo;
   rec.sbnd_timings     = srsbndtiminginfo;
   rec.soft_trig        = srsbndsofttrig;