Merge branch 'release/v10_10_05'

Author: kjplows

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.10.03.02 LANGUAGES CXX)
+project(sbncode VERSION 10.10.05 LANGUAGES CXX)
 
 message(STATUS "\n\n  ==========================   ${PROJECT_NAME}   ==========================")
 

sbncode/CAFMaker/CAFMakerParams.h

@@ -321,6 +321,12 @@ namespace caf
       "crthit" // icarus
     };
 
+    Atom<string> CRTSimChanLabel {
+      Name("CRTSimChanLabel"),
+      Comment("Label of AuxDetSimChannels."),
+      "genericcrt" // icarus
+    };
+
     Atom<string> CRTTrackLabel {
       Name("CRTTrackLabel"),
       Comment("Label of sbn CRT tracks."),
@@ -417,6 +423,30 @@ namespace caf
       "emuTriggerUnshifted"
     };
 
+    Atom<art::InputTag> MonPulsesTriggerLabel {
+      Name("MonPulsesTriggerLabel"),
+      Comment("Label of trigger emulation product MonPulses (number of PMT pairs above threshold for all channels) for all flashes."),
+      art::InputTag("opdaq", "MonPulses", "DetSim")
+    };
+
+    Atom<art::InputTag> MonPulseSizesTriggerLabel {
+      Name("MonPulseSizesTriggerLabel"),
+      Comment("Label of trigger emulation product MonPulses Sizes, which gives the length of each trigger response in MonPulses."),
+      art::InputTag("opdaq", "MonPulseSizes", "DetSim")
+    };
+
+    Atom<art::InputTag> PairsTriggerLabel {
+      Name("PairsTriggerLabel"),
+      Comment("Label of number of PMT pairs over threshold."),
+      art::InputTag("opdaq", "pairsOverThreshold", "DetSim")
+    };
+
+    Atom<art::InputTag> EmulatedTriggerLabel {
+      Name("EmulatedTriggerLabel"),
+      Comment("Label of bool of passing the trigger."),
+      art::InputTag("opdaq", "triggerEmulation", "DetSim")
+    };
+
     Atom<string> FlashTrigLabel {
       Name("FlashTrigLabel"),
       Comment("Label of bool of passing flash trigger."),

sbncode/CAFMaker/CAFMaker_module.cc

@@ -141,6 +141,7 @@
 #include "sbnanaobj/StandardRecord/Flat/FlatRecord.h"
 #include "lardataobj/RawData/ExternalTrigger.h"
 #include "lardataobj/RawData/TriggerData.h"
+#include "lardataobj/Simulation/AuxDetSimChannel.h"
 
 // // CAFMaker
 #include "sbncode/CAFMaker/AssociationUtil.h"
@@ -288,6 +289,11 @@ class CAFMaker : public art::EDProducer {
   std::vector<std::vector<geo::BoxBoundedGeo>> fTPCVolumes;
   std::vector<geo::BoxBoundedGeo> fActiveVolumes;
 
+  // CRT geometry info
+  //
+  // ICARUS
+  std::map<std::pair<int, int>, int> fFEBChannel2AuxDetID;
+
   // random number generator for fake reco
   CLHEP::HepRandomEngine& fFakeRecoRandomEngine;
 
@@ -316,6 +322,7 @@ class CAFMaker : public art::EDProducer {
   double GetBlindPOTScale() const;
 
   void InitVolumes(); ///< Initialize volumes from Gemotry service
+  void InitCRTMapping(); ///< Initialize CRT mapping
 
   void FixPMTReferenceTimes(StandardRecord &rec, double PMT_reference_time);
   void FixCRTReferenceTimes(StandardRecord &rec, double CRTT0_reference_time, double CRTT1_reference_time);
@@ -424,6 +431,9 @@ class CAFMaker : public art::EDProducer {
   // setup volume definitions
   InitVolumes();
 
+  // setup CRT mapping
+  InitCRTMapping();
+
   fSaveGENIEEventRecord = fParams.SaveGENIEEventRecord();
 
 }
@@ -597,6 +607,36 @@ void CAFMaker::FixCRTReferenceTimes(StandardRecord &rec, double CRTT0_reference_
 
 }
 
+void CAFMaker::InitCRTMapping() {
+
+  // ICARUS
+  cet::search_path searchPath("FW_SEARCH_PATH");
+  std::string fullFileName;
+  searchPath.find_file("feb_map.txt",fullFileName);
+  std::ifstream fin;
+  fin.open(fullFileName, std::ios::in);
+
+  if (fin.good()) {
+    std::string line;
+
+    while(std::getline(fin,line)) {
+      std::vector<std::string> row;
+      std::string word;
+
+      std::stringstream s(line);
+      while (std::getline(s, word, ',')) {
+        row.push_back(word);
+      }
+      int auxDetID = std::stoi(row[0]); // auxDetID
+      int mac5 = std::stoi(row[1]); // FEB ID
+      int chan = std::stoi(row[2]); // FEB channel
+      fFEBChannel2AuxDetID[std::make_pair(mac5, chan)] = auxDetID;
+    }
+
+    fin.close();
+  }
+}
+
 void CAFMaker::InitVolumes() {
   const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
 
@@ -1600,7 +1640,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   //#######################################################
   // Fill detector & reco
   //#######################################################
-
+  
   //Beam gate and Trigger info
   art::Handle<sbn::ExtraTriggerInfo> extratrig_handle;
   GetByLabelStrict(evt, fParams.TriggerLabel().encode(), extratrig_handle);
@@ -1612,16 +1652,40 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   if (!isRealData)
     GetByLabelStrict(evt, fParams.UnshiftedTriggerLabel().encode(), unshifted_trig_handle);
 
+  //Trigger emulation handles
+  art::Handle<std::vector<int>> monpulses_handle;
+  GetByLabelStrict(evt, fParams.MonPulsesTriggerLabel().encode(), monpulses_handle);
+
+  art::Handle<std::vector<int>> monpulse_sizes_handle;
+  GetByLabelStrict(evt, fParams.MonPulseSizesTriggerLabel().encode(), monpulse_sizes_handle);
+
+  art::Handle<int> pairs_handle;
+  GetByLabelStrict(evt, fParams.PairsTriggerLabel().encode(), pairs_handle);
+
+  art::Handle<bool> trigemu_handle;
+  GetByLabelStrict(evt, fParams.EmulatedTriggerLabel().encode(), trigemu_handle);
+
+  // Check trigger handles
   const bool isValidTrigger = extratrig_handle.isValid() && trig_handle.isValid() && trig_handle->size() == 1;
   const bool isValidUnshiftedTrigger = unshifted_trig_handle.isValid() && unshifted_trig_handle->size() == 1;
+  const bool isValidEmulationTrigger = monpulses_handle.isValid() && monpulse_sizes_handle.isValid() && pairs_handle.isValid() && trigemu_handle.isValid();
 
   const double triggerShift = (isValidUnshiftedTrigger && isValidTrigger)?
     unshifted_trig_handle->at(0).TriggerTime() - trig_handle->at(0).TriggerTime() : 0.;
 
+  // Fill local ExtraTriggerInfo struct
+  sbn::ExtraTriggerInfo extratrig;
+  if (extratrig_handle.isValid()) extratrig = *extratrig_handle;
+
   caf::SRTrigger srtrigger;
   if (isValidTrigger) {
       FillTrigger(*extratrig_handle, trig_handle->at(0), srtrigger, triggerShift);
   }
+  // Fill trigger emulation information
+  if (isValidEmulationTrigger) { 
+      FillTriggerEmulation(monpulses_handle, monpulse_sizes_handle, pairs_handle, trigemu_handle, srtrigger);
+  }
+
   // 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.
@@ -1652,16 +1716,36 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   caf::SRSBNDFrameShiftInfo srsbndframeshiftinfo;
   caf::SRSBNDTimingInfo srsbndtiminginfo;
 
+  // Mapping of (feb, channel) to truth information (AuxDetSimChannel) -- filled for ICARUS
+  std::map<std::pair<int, int>, sim::AuxDetSimChannel> crtsimchanmap;
+
   if(fDet == kICARUS)
     {
       art::Handle<std::vector<sbn::crt::CRTHit>> crthits_handle;
       GetByLabelStrict(evt, fParams.CRTHitLabel(), crthits_handle);
+
+      art::Handle<std::vector<sim::AuxDetSimChannel>> auxdetsimchan_handle;
+      GetByLabelStrict(evt, fParams.CRTSimChanLabel(), auxdetsimchan_handle);
+
       // fill into event
       if (crthits_handle.isValid()) {
         const std::vector<sbn::crt::CRTHit> &crthits = *crthits_handle;
+
+        std::vector<sim::AuxDetSimChannel> empty;
+        const std::vector<sim::AuxDetSimChannel> &crtsimchanvec = (auxdetsimchan_handle.isValid()) ? *auxdetsimchan_handle : empty;
+        // Turn the AuxDetSimChannel's into a map that can be looked up from a CRT Hit
+        for (const sim::AuxDetSimChannel &crtsimchan: crtsimchanvec) {
+          for (auto const &map_pair: fFEBChannel2AuxDetID) {
+            if (map_pair.second == (int)crtsimchan.AuxDetID()) {
+              crtsimchanmap[map_pair.first] = crtsimchan;
+              break;
+            }
+          }
+        }
+
         for (unsigned i = 0; i < crthits.size(); i++) {
           srcrthits.emplace_back();
-          FillCRTHit(crthits[i], fParams.CRTUseTS0(), CRT_T0_reference_time, CRT_T1_reference_time, srcrthits.back());
+          FillCRTHit(crthits[i], fParams.CRTUseTS0(), CRT_T0_reference_time, CRT_T1_reference_time, crtsimchanmap, srcrthits.back());
         }
       }
 
@@ -2359,7 +2443,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
             crthittagginginfo = fmCRTHitMatchInfo.at(iPart);
           }          
 
-          FillTrackCRTHit(fmCRTHitMatch.at(iPart), crthitmatch, crthittagginginfo, fParams.CRTUseTS0(), CRT_T0_reference_time, CRT_T1_reference_time, trk);
+          FillTrackCRTHit(fmCRTHitMatch.at(iPart), crthitmatch, crthittagginginfo, fParams.CRTUseTS0(), CRT_T0_reference_time, CRT_T1_reference_time, crtsimchanmap, trk);
         }
         // NOTE: SEE TODO AT fmCRTTrackMatch
         if (fmCRTTrackMatch.isValid() && fDet == kICARUS) {

sbncode/CAFMaker/CMakeLists.txt

@@ -36,6 +36,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker
                   sbnobj::Common_Reco
                   sbnobj::Common_Analysis
                   sbnobj::Common_PMT_Data
+                  sbnobj::Common_Trigger
                   sbnobj::SBND_CRT
                   sbnobj::SBND_Timing
                   lardataalg::DetectorInfo

sbncode/CAFMaker/FillReco.cxx

@@ -68,6 +68,7 @@ namespace caf
                   bool use_ts0,
                   int64_t CRT_T0_reference_time, // ns, signed
                   double CRT_T1_reference_time, // us
+                  const std::map<std::pair<int, int>, sim::AuxDetSimChannel> &crtsimchanmap,
                   caf::SRCRTHit &srhit,
                   bool allowEmpty) {
 
@@ -86,6 +87,52 @@ namespace caf
     srhit.pe = hit.peshit;
     srhit.plane = hit.plane;
 
+    // lookup truth matching information
+    std::map<int, float> true_energy;
+    std::set<std::pair<int, int>> checked; // keep track of what AuxDetSimChannels we have looked at
+    for (auto const &mac_to_pes: hit.pesmap) {
+      int mac = (int)mac_to_pes.first;
+      for (const std::pair<int,float> &chan_pe: mac_to_pes.second) {
+        int chan = chan_pe.first;
+
+        // check for the 3-pos Minos ("m") modules, or the 1-pos other modules
+        bool is_minos_module = mac < 94;
+        int pos = is_minos_module ? (chan/10 + 1) : 1;
+
+        std::pair<int, int> key {mac, pos};
+        if (checked.count(key)) continue; // already looked here
+
+        checked.insert(key);
+   
+        if (crtsimchanmap.count(key)) {
+          const sim::AuxDetSimChannel &crtsimchan = crtsimchanmap.at(key);
+          for (const sim::AuxDetIDE &ide: crtsimchan.AuxDetIDEs()) {
+            double hit_time_us = srhit.time;
+            double tru_time_us = ide.entryT/1e3; // ns -> us 
+
+            // 1us cut on truth matching
+            if (abs(hit_time_us - tru_time_us) < 1) {
+            // std::cout << "Hit at time: " << (srhit.time/1e3) << " " << (srhit.t0/1e3) << " " << (srhit.t1/1e3) << " pes: " << chan_pe.second << " of " << srhit.pe << " on FEB: " << mac << " channel: " << chan << " matched to AuxDetID: " << crtsimchan.AuxDetID() << " G4ID: " << ide.trackID << " E: " << ide.energyDeposited << " T: " << (ide.entryT/1e6) << std::endl;
+              true_energy[ide.trackID] += ide.energyDeposited;
+            }
+
+          }
+        }
+      }
+    }
+
+    // sort results by energy
+    std::vector<std::pair<float, int>> true_energy_list;
+    for (auto const &pair: true_energy) true_energy_list.push_back(std::make_pair(pair.second, pair.first));
+    std::sort(true_energy_list.begin(), true_energy_list.end(), std::greater<>());
+
+    // Save to the SRCRTHit truth
+    for (auto const &pair: true_energy_list) {
+      srhit.truth.match_e.push_back(pair.first);
+      srhit.truth.match_id.push_back(pair.second);
+    }
+    if (true_energy_list.size() > 0) srhit.truth.bestmatch_id = true_energy_list[0].second;
+
   }
 
   void FillCRTTrack(const sbn::crt::CRTTrack &track,
@@ -658,6 +705,7 @@ namespace caf
                        bool use_ts0,
                        int64_t CRT_T0_reference_time, // ns, signed
                        double CRT_T1_reference_time, // us
+                       const std::map<std::pair<int, int>, sim::AuxDetSimChannel> &crtsimchanmap,
                        caf::SRTrack &srtrack,
                        bool allowEmpty)
   {
@@ -684,7 +732,7 @@ namespace caf
       srtrack.crthit.hit.plane = t0match[0]->fID;
     }
     if (hitmatch.size()) {
-      FillCRTHit(*hitmatch[0], use_ts0, CRT_T0_reference_time, CRT_T1_reference_time, srtrack.crthit.hit, allowEmpty);
+      FillCRTHit(*hitmatch[0], use_ts0, CRT_T0_reference_time, CRT_T1_reference_time, crtsimchanmap, srtrack.crthit.hit, allowEmpty);
     }
   }
 

sbncode/CAFMaker/FillReco.h

@@ -25,6 +25,7 @@
 #include "lardataobj/AnalysisBase/T0.h"
 #include "lardataobj/RecoBase/PFParticleMetadata.h"
 #include "lardataobj/RecoBase/MCSFitResult.h"
+#include "lardataobj/Simulation/AuxDetSimChannel.h"
 #include "larrecodnn/CVN/func/Result.h"
 #include "sbnobj/Common/Reco/Stub.h"
 #include "sbnobj/Common/Reco/RangeP.h"
@@ -193,6 +194,7 @@ namespace caf
                        bool use_ts0,
                        int64_t CRT_T0_reference_time, // ns, signed
                        double CRT_T1_reference_time, // us
+                       const std::map<std::pair<int, int>, sim::AuxDetSimChannel> &crtsimchanmap,
                        caf::SRTrack &srtrack,
                        bool allowEmpty = false);
 
@@ -260,6 +262,7 @@ namespace caf
                   bool use_ts0,
                   int64_t CRT_T0_reference_time, // ns, signed
                   double CRT_T1_reference_time, // us
+                  const std::map<std::pair<int, int>, sim::AuxDetSimChannel> &crtsimchanmap,
                   caf::SRCRTHit &srhit,
                   bool allowEmpty = false);
   void FillCRTTrack(const sbn::crt::CRTTrack &track,

sbncode/CAFMaker/FillTrigger.cxx

@@ -26,6 +26,7 @@ namespace caf
     triggerInfo.trigger_count = addltrig_info.triggerCount;
     triggerInfo.gate_count = addltrig_info.gateCount;
     triggerInfo.gate_delta = addltrig_info.gateCountFromPreviousTrigger;
+
   }
 
   void FillTriggerMC(double absolute_time, caf::SRTrigger& triggerInfo) {
@@ -35,7 +36,6 @@ namespace caf
     // Set this to 0 since the "MC" trigger is (for now) always at the spill time
     triggerInfo.trigger_within_gate = 0.;
 
-    // TODO: fill others?
   }
 
   void FillTriggerSBND(caf::SRSBNDTimingInfo& timingInfo, caf::SRTrigger& triggerInfo){
@@ -49,4 +49,17 @@ namespace caf
     triggerInfo.trigger_within_gate = diff_ts;
   }
 
+  void FillTriggerEmulation(art::Handle<std::vector<int>> const& monpulsesFlat,
+                             art::Handle<std::vector<int>> const& monpulseSizes,
+                             art::Handle<int> const& numPairs,
+                             art::Handle<bool> const& passedTrig,
+                             caf::SRTrigger& triggerInfo) {
+
+    triggerInfo.monpulses_flat = *monpulsesFlat;
+    triggerInfo.monpulse_sizes = *monpulseSizes;
+    triggerInfo.num_pairs_over_threshold = *numPairs;
+    triggerInfo.passed_trigger = *passedTrig;
+
+  }
+
 }

sbncode/CAFMaker/FillTrigger.h

@@ -6,6 +6,7 @@
 #include "sbnanaobj/StandardRecord/SRTrigger.h"
 #include "sbnanaobj/StandardRecord/SRSBNDTimingInfo.h"
 #include "lardataobj/RawData/TriggerData.h"
+#include "art/Framework/Principal/Handle.h"
 
 #include <vector>
 
@@ -20,6 +21,12 @@ namespace caf
   void FillTriggerMC(double absolute_time, caf::SRTrigger& triggerInfo);
 
   void FillTriggerSBND(caf::SRSBNDTimingInfo& timingInfo, caf::SRTrigger& triggerInfo);
+
+  void FillTriggerEmulation(art::Handle<std::vector<int>> const& monpulsesFlat,
+                             art::Handle<std::vector<int>> const& monpulseSizes,
+                             art::Handle<int> const& numPairs,
+                             art::Handle<bool> const& passedTrig,
+                             caf::SRTrigger& triggerInfo);
 }
 
 #endif