diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index b38ed1f6a..fe35bad7a 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/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 @@ -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 offset_vector_ns(rec.mc.nu.size(), std::numeric_limits::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> mcflux_handle; - GetByLabelStrict(evt, std::string("generator"), mcflux_handle); + GetByLabelStrict(evt, fParams.GenLabel(), mcflux_handle); std::vector> mcfluxes; if (mcflux_handle.isValid()) { art::fill_ptr_vector(mcfluxes, mcflux_handle); } + art::Handle> dk2nu_handle; + GetByLabelStrict(evt, fParams.GenLabel(), dk2nu_handle); + + std::vector const* dk2nus = + dk2nu_handle.isValid() ? dk2nu_handle.product() : nullptr; + // get the MCReco for the fake-reco art::Handle> 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; iat(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,14 +2821,23 @@ 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); @@ -2709,11 +2845,11 @@ void CAFMaker::produce(art::Event& evt) noexcept { 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); } } diff --git a/sbncode/CAFMaker/CMakeLists.txt b/sbncode/CAFMaker/CMakeLists.txt index 7a5162082..a53793dd5 100644 --- a/sbncode/CAFMaker/CMakeLists.txt +++ b/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 diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 4d1d794d0..3cf135a04 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/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 mctruth, const simb::MCFlux &mcflux, + const bsim::Dk2Nu &dk2nu, const simb::GTruth& gtruth, const std::vector &srparticles, const std::map>> &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::signaling_NaN(); srneutrino.parent_pdg = mcflux.fptype; srneutrino.parent_dcy_mode = mcflux.fndecay; srneutrino.prod_vtx.x = mcflux.fvx; diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index 419ed6ee9..e7579d799 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/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 mctruth, const simb::MCFlux &mcflux, + const bsim::Dk2Nu &dk2nu, const simb::GTruth& gtruth, const std::vector &srparticles, const std::map>> &id_to_truehit_map,