Gen2 Production PR: Manual Fix for GENIE timing bug in CAFMaker + Extend existing-data shift in CAFMaker

sbncode/CAFMaker/CAFMaker_module.cc

@@ -42,6 +42,8 @@
 
 #include "ifdh_art/IFDHService/IFDH_service.h"
 
+#include "dk2nu/tree/dk2nu.h"
+
 // ROOT includes
 #include "TFile.h"
 #include "TH1D.h"
@@ -136,6 +138,7 @@
 
 #include "canvas/Persistency/Provenance/ProcessConfiguration.h"
 #include "larcoreobj/SummaryData/POTSummary.h"
+#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"
 
 // StandardRecord
 #include "sbnanaobj/StandardRecord/StandardRecord.h"
@@ -338,6 +341,10 @@ class CAFMaker : public art::EDProducer {
   void SBNDShiftCRTReference(StandardRecord &rec, double SBNDFrame) const;
   void SBNDShiftPMTReference(StandardRecord &rec, double SBNDFrame) const;
 
+  /// Corrects the dk2nu v01_11_00 vertex-time bug (fixed in dk2nu v01_12_00, larsoft v10_23_00)
+  /// by adding each GENIE neutrino's `prod_time`/`dk2gen` offset to its truth time and matched slice.
+  void CorrectMCTiming(StandardRecord &rec) const;
+
   /// Equivalent of FindManyP except a return that is !isValid() prints a
   /// messsage and aborts if StrictMode is true.
   template <class T, class U>
@@ -536,14 +543,103 @@ void CAFMaker::SBNDShiftCRTReference(StandardRecord &rec, double SBNDFrame) cons
     trk.time += SBNDFrame; //ns
   }
 
-  //TODO: CRT Space Point and Track Match
   for (SRSlice &slc: rec.slc){
     for (SRPFP &pfp: slc.reco.pfp){
+      //PFP track-matched CRT
       if(!std::isnan(pfp.trk.crtspacepoint.score)) pfp.trk.crtspacepoint.spacepoint.time += SBNDFrame;
-
       if(!std::isnan(pfp.trk.crtsbndtrack.score)) pfp.trk.crtsbndtrack.track.time += SBNDFrame;
     }
+
+    // CRUMBS CRT input features [us]; unfilled default is NaN, not -9999.
+    if (!std::isnan(slc.crumbs_result.crt.tracktime)) slc.crumbs_result.crt.tracktime += SBNDFrame / 1000.;
+    if (!std::isnan(slc.crumbs_result.crt.sptime))    slc.crumbs_result.crt.sptime    += SBNDFrame / 1000.;
   }
+
+  // Note: SRSBNDCRTVeto.sp_time data is already RWM-shifted in CRTVeto module
+}
+
+void CAFMaker::CorrectMCTiming(StandardRecord &rec) const {
+
+  mf::LogDebug("CAFMaker") << "CorrectMCTiming: running";
+
+  // Recovers the dk2nu v01_11_00 SetT(0) bug (fixed in v01_12_00) by adding back
+  // the missing meson decay time + tof from decay to flux window.
+
+  // Shifts a true particle's gen/start/end times, skipping any left at its unfilled default.
+  auto shiftMCParticleTimes = [](SRTrueParticle &part, double offset_us) {
+    if (part.genT   > -9998.) part.genT   += offset_us;
+    if (part.startT > -9998.) part.startT += offset_us;
+    if (part.endT   > -9998.) part.endT   += offset_us;
+  };
+  // Shifts a true interaction's own time plus all of its primary particles' times.
+  auto shiftTruthTimes = [&shiftMCParticleTimes](SRTrueInteraction &truth, double offset_us) {
+    truth.time += offset_us;
+    for (SRTrueParticle &part: truth.prim) shiftMCParticleTimes(part, offset_us);
+  };
+
+  // Per-neutrino offset, indexed the same as rec.mc.nu / SRSlice.tmatch.index.
+  std::vector<double> offset_vector_ns(rec.mc.nu.size(), std::numeric_limits<double>::signaling_NaN());
+
+  for (unsigned i = 0; i < rec.mc.nu.size(); i++) {
+    SRTrueInteraction &nu = rec.mc.nu[i];
+
+    if (nu.generator != caf::kGENIE) continue; //GENIE-dk2nu-specific correction;
+    if (std::isnan(nu.prod_time)) continue; // no dk2nu info for this sample/interaction
+    if (nu.dk2gen <= -998.) continue; // simb::MCFlux::fdk2gen unfilled default is -999.
+
+    offset_vector_ns[i] = nu.prod_time * 1000. + nu.dk2gen / util::kc; // prod_time is us, dk2gen/kc is ns
+    double const offset_us = offset_vector_ns[i] / 1000.;
+
+    shiftTruthTimes(nu, offset_us);
+  }
+
+  // Apply each slice's matched neutrino's own offset (SRSlice.tmatch.index)
+  // Unmatched/cosmic slices (tmatch.index < 0) are left untouched.
+  for (SRSlice &slc: rec.slc) {
+    if (slc.tmatch.index < 0 || std::isnan(offset_vector_ns[slc.tmatch.index])) continue;
+
+    double const offset_ns = offset_vector_ns[slc.tmatch.index];
+    double const offset_us = offset_ns / 1000.;
+
+    // slc.truth is a copy of rec.mc.nu[tmatch.index] taken earlier (FillTrue.cxx:372).
+    shiftTruthTimes(slc.truth, offset_us);
+
+    // Opt0 FM
+    if (!std::isnan(slc.opt0.time))     slc.opt0.time     += offset_us; // NaN default
+    if (!std::isnan(slc.opt0_sec.time)) slc.opt0_sec.time += offset_us; // NaN default
+
+    // BFM
+    if (slc.barycenterFM.flashTime     != -9999.) slc.barycenterFM.flashTime     += offset_us; // -9999. default
+    if (slc.barycenterFM.flashFirstHit != -9999.) slc.barycenterFM.flashFirstHit += offset_us; // -9999. default
+
+    // BFM + Light Propagation Module
+    if (slc.correctedOpFlash.OpFlashT0 != -9999.) { // -9999. default
+      slc.correctedOpFlash.OpFlashT0          += offset_us;
+      slc.correctedOpFlash.OpFlashT0Corrected += offset_us;
+    }
+
+    // Simple FM: PMT & XARAPUCA
+    if (slc.fmatch.time      != -9999.) slc.fmatch.time      += offset_us; // -9999. default
+    if (slc.fmatchop.time    != -9999.) slc.fmatchop.time    += offset_us; // -9999. default
+    if (slc.fmatchara.time   != -9999.) slc.fmatchara.time   += offset_us; // -9999. default
+    if (slc.fmatchopara.time != -9999.) slc.fmatchopara.time += offset_us; // -9999. default
+
+    //PFP track-matched CRT
+    for (SRPFP &pfp: slc.reco.pfp) {
+      if (!std::isnan(pfp.trk.crtspacepoint.score)) pfp.trk.crtspacepoint.spacepoint.time += offset_ns;
+      if (!std::isnan(pfp.trk.crtsbndtrack.score))  pfp.trk.crtsbndtrack.track.time      += offset_ns;
+    }
+
+    // CRUMBS BDT
+    if (!std::isnan(slc.crumbs_result.pds.fmtime))    slc.crumbs_result.pds.fmtime    += offset_us; // NaN default
+    if (!std::isnan(slc.crumbs_result.crt.tracktime)) slc.crumbs_result.crt.tracktime += offset_us; // NaN default
+    if (!std::isnan(slc.crumbs_result.crt.sptime))    slc.crumbs_result.crt.sptime    += offset_us; // NaN default
+  }
+
+  // TODO: Non-correction list
+  // 1. No slice-truth-matching path exists: rec.opflashes, rec.crt_spacepoints, rec.sbnd_crt_tracks
+  // 2. TPC timing variables: SRCaloPoint.t, SRHit.peakTime, SRPFP.t0, SRBlip.time/.timeTick, SRBlipHitClust.time/.startTime/.endTime/.timespan/.timeTick
+  // 3. Has their own timing shift in local modules: SRSBNDCRTVeto.sp_time, SRSoftwareTrigger.flash_peaktime
 }
 
 void CAFMaker::SBNDShiftPMTReference(StandardRecord &rec, double SBNDFrame) const {
@@ -556,10 +652,25 @@ void CAFMaker::SBNDShiftPMTReference(StandardRecord &rec, double SBNDFrame) cons
     opf.firsttime += SBNDFrame_us;
   }
 
-  //OpT0 match to slice
+  // OpT0 match to slice; guards match each field's own unfilled convention (see CorrectMCTiming)
   for (SRSlice &s: rec.slc) {
-    s.opt0.time += SBNDFrame_us;
+    if (!std::isnan(s.opt0.time))     s.opt0.time     += SBNDFrame_us; // NaN default
+    if (!std::isnan(s.opt0_sec.time)) s.opt0_sec.time += SBNDFrame_us; // NaN default
+
+    if (s.barycenterFM.flashTime     != -9999.) s.barycenterFM.flashTime     += SBNDFrame_us; // -9999. default
+    if (s.barycenterFM.flashFirstHit != -9999.) s.barycenterFM.flashFirstHit += SBNDFrame_us; // -9999. default
+
+    if (s.fmatch.time      != -9999.) s.fmatch.time      += SBNDFrame_us; // -9999. default
+    if (s.fmatchop.time    != -9999.) s.fmatchop.time    += SBNDFrame_us; // -9999. default
+    if (s.fmatchara.time   != -9999.) s.fmatchara.time   += SBNDFrame_us; // -9999. default
+    if (s.fmatchopara.time != -9999.) s.fmatchopara.time += SBNDFrame_us; // -9999. default
+
+    if (!std::isnan(s.crumbs_result.pds.fmtime)) s.crumbs_result.pds.fmtime += SBNDFrame_us; // NaN default
   }
+
+  // Note:
+  // 1. SRSoftwareTrigger.flash_peaktime data: already RWM-shifted in SoftwareTrigger module
+  // 2. correctedOpFlash.OpFlashT0/.OpFlashT0Corrected: already RWM-shifted in LightPropagation module
 }
 
 void CAFMaker::FixPMTReferenceTimes(StandardRecord &rec, double PMT_reference_time) {
@@ -1447,13 +1558,19 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   }
 
   art::Handle<std::vector<simb::MCFlux>> mcflux_handle;
-  GetByLabelStrict(evt, std::string("generator"), mcflux_handle);
+  GetByLabelStrict(evt, fParams.GenLabel(), mcflux_handle);
 
   std::vector<art::Ptr<simb::MCFlux>> mcfluxes;
   if (mcflux_handle.isValid()) {
     art::fill_ptr_vector(mcfluxes, mcflux_handle);
   }
 
+  art::Handle<std::vector<bsim::Dk2Nu>> dk2nu_handle;
+  GetByLabelStrict(evt, fParams.GenLabel(), dk2nu_handle);
+
+  std::vector<bsim::Dk2Nu> const* dk2nus =
+    dk2nu_handle.isValid() ? dk2nu_handle.product() : nullptr;
+
   // get the MCReco for the fake-reco
   art::Handle<std::vector<sim::MCTrack>> mctrack_handle;
   GetByLabelStrict(evt, std::string("mcreco"), mctrack_handle);
@@ -1537,9 +1654,13 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   // holder for invalid MCFlux
   simb::MCFlux badflux; // default constructor gives nonsense values
 
+  // holder for invalid Dk2Nu (default-constructed ancestor vector is empty -- treated as "absent")
+  bsim::Dk2Nu const baddk2nu;
+
   for (size_t i=0; i<mctruths.size(); i++) {
     auto const& mctruth = mctruths.at(i);
     const simb::MCFlux &mcflux = (mcfluxes.size()) ? *mcfluxes.at(i) : badflux;
+    const bsim::Dk2Nu &dk2nu = dk2nus ? dk2nus->at(i) : baddk2nu;
 
     simb::GTruth gtruth;
     bool ok = GetAssociatedProduct(fmp_gtruth, i, gtruth);
@@ -1554,7 +1675,13 @@ void CAFMaker::produce(art::Event& evt) noexcept {
 
     if ( !isRealData ){
 
-      FillTrueNeutrino(mctruth, mcflux, gtruth, true_particles, id_to_truehit_map, srtruthbranch.nu.back(), i, fActiveVolumes);
+      FillTrueNeutrino(mctruth, mcflux, dk2nu, gtruth, true_particles, id_to_truehit_map, srtruthbranch.nu.back(), i, fActiveVolumes);
+
+      SRTrueInteraction &nu = srtruthbranch.nu.back();
+      caf::generator_ &generator = nu.generator;
+      if (ok)                             generator = caf::kGENIE;
+      else if (i < mevprtl_truths.size()) generator = caf::kMeVPrtl;
+      else                                generator = caf::kUnknownGenerator;
 
       srtruthbranch.nu.back().genie_evtrec_idx = fGenieEventCounter;
 
@@ -2694,26 +2821,35 @@ void CAFMaker::produce(art::Event& evt) noexcept {
 
   FixPMTReferenceTimes(rec, PMT_reference_time);
 
+  // SBND MC: recover the dk2nu v01_11_00 SetT(0) vertex-time bug (gen2-tag production)
+  if (!isRealData && (fDet == kSBND)) {
+    CorrectMCTiming(rec);
+  }
+
   // TODO: TPC?
 
-  // SBND: Fix the Reference time in data depending on the stream
-  // For more information, see: 
-  // https://sbn-docdb.fnal.gov/cgi-bin/sso/RetrieveFile?docid=43090
+  // SBND: Fix the reference time in data exclusively for Gen2 Sample (v10_14_02)
+  // For information, see: https://sbn-docdb.fnal.gov/cgi-bin/sso/RetrieveFile?docid=43090
+  // Otherwise, this code block is legacy for sbncode > v10_20_09
 
   if (isRealData && (fDet == kSBND) && fSubRunPOT > 0)
   {
+    mf::LogDebug("CAFMaker") << "run " << evt.run() << " subRun " << evt.subRun()
+                             << " fSubRunPOT=" << fSubRunPOT
+                             << " frameApplyAtCaf=" << rec.sbnd_frames.frameApplyAtCaf;
+
     // Fill trigger info
     FillTriggerSBND(srsbndtiminginfo, srtrigger);
 
     // Shift timing reference frame
     if (!std::isnan(rec.sbnd_frames.frameApplyAtCaf) && (rec.sbnd_frames.frameApplyAtCaf != 0.0)){
       mf::LogInfo("CAFMaker") << "Setting Reference Timing for timing object in SBND \n"
                               << "    Shift Apply At Caf Level = " << rec.sbnd_frames.frameApplyAtCaf << " ns\n";
-      
-      //shift reference frame for CRT objects: crt trk, crt sp, crt sp match, crt trk match
+
+      //shift reference frame for CRT objects: crt trk, crt sp, crt sp match, crt trk match, CRUMBS CRT features
       SBNDShiftCRTReference(rec, rec.sbnd_frames.frameApplyAtCaf);
 
-      //shift reference frame for PMT objects: opflash, opt0
+      //shift reference frame for PMT objects: opflash, opt0/opt0_sec, barycenter flash match, simple flash matches, CRUMBS PDS feature
       SBNDShiftPMTReference(rec, rec.sbnd_frames.frameApplyAtCaf);
     }
   }

sbncode/CAFMaker/CMakeLists.txt

@@ -23,6 +23,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker
                   ROOT::EG
                   CLHEP::Random
                   caf_RecoUtils
+                  dk2nu::Tree
                   lardataobj::AnalysisBase
                   lardataobj::RecoBase
                   lardataobj::Simulation

sbncode/CAFMaker/FillTrue.cxx

@@ -1,10 +1,13 @@
 #include "FillTrue.h"
 
+#include "dk2nu/tree/dk2nu.h"
+
 #include "larcorealg/GeoAlgo/GeoAlgo.h"
 #include "larcorealg/Geometry/GeometryCore.h"
 #include "larcorealg/Geometry/WireReadoutGeom.h"
 #include "larevt/SpaceCharge/SpaceCharge.h"
 #include "larcore/CoreUtils/ServiceUtil.h"
+#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"
 #include "RecoUtils/RecoUtils.h"
 
 #include "CLHEP/Random/RandGauss.h"
@@ -385,6 +388,7 @@ namespace caf {
    srtruth.position.y = truth.decay_pos.Y();
    srtruth.position.z = truth.decay_pos.Z();
    srtruth.time = truth.decay_pos.T();
+   srtruth.prod_time = truth.meson_dpos_beamcoord.T(); //Meson time from flux file.
 
    srtruth.momentum.x = truth.mevprtl_mom.X();
    srtruth.momentum.y = truth.mevprtl_mom.Y();
@@ -457,6 +461,7 @@ namespace caf {
   //------------------------------------------------
   void FillTrueNeutrino(const art::Ptr<simb::MCTruth> mctruth,
       const simb::MCFlux &mcflux,
+      const bsim::Dk2Nu &dk2nu,
       const simb::GTruth& gtruth,
       const std::vector<caf::SRTrueParticle> &srparticles,
       const std::map<int, std::vector<art::Ptr<recob::Hit>>> &id_to_truehit_map,
@@ -548,7 +553,10 @@ namespace caf {
 
     // Set the MCFlux stuff
     srneutrino.initpdg = mcflux.fntype;
-    srneutrino.baseline = mcflux.fdk2gen + mcflux.fgen2vtx;
+    srneutrino.baseline = (mcflux.fdk2gen + mcflux.fgen2vtx) * util::kMeterToCentimeter;
+    srneutrino.dk2gen = mcflux.fdk2gen * util::kMeterToCentimeter;
+    srneutrino.prod_time = (!dk2nu.ancestor.empty())
+        ? dk2nu.ancestor.back().startt / 1000. /* ns -> us */ : std::numeric_limits<float>::signaling_NaN();
     srneutrino.parent_pdg = mcflux.fptype;
     srneutrino.parent_dcy_mode = mcflux.fndecay;
     srneutrino.prod_vtx.x = mcflux.fvx;

sbncode/CAFMaker/FillTrue.h

@@ -33,6 +33,8 @@
 #include "sbnanaobj/StandardRecord/StandardRecord.h"
 #include "sbnanaobj/StandardRecord/SRMeVPrtl.h"
 
+namespace bsim { class Dk2Nu; }
+
 namespace caf
 {
   struct HitsEnergy {
@@ -86,6 +88,7 @@ namespace caf
 
   void FillTrueNeutrino(const art::Ptr<simb::MCTruth> mctruth,
                         const simb::MCFlux &mcflux,
+                        const bsim::Dk2Nu &dk2nu,
                         const simb::GTruth& gtruth,
                         const std::vector<caf::SRTrueParticle> &srparticles,
                         const std::map<int, std::vector<art::Ptr<recob::Hit>>> &id_to_truehit_map,