diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index b38ed1f6a..03afbd8f2 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" @@ -338,6 +340,8 @@ class CAFMaker : public art::EDProducer { void SBNDShiftCRTReference(StandardRecord &rec, double SBNDFrame) const; void SBNDShiftPMTReference(StandardRecord &rec, double SBNDFrame) const; + 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 +540,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 { + + std::cout << "CAFMaker::CorrectMCTiming: running" << std::endl; + + // 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. + constexpr double kSpeedOfLight = 29.9792458; // cm/ns + + // Per-neutrino offset, indexed the same as rec.mc.nu / SRSlice.tmatch.index. + std::vector offset_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_ns[i] = nu.prod_time + nu.dk2gen * 100. / kSpeedOfLight; + double const delta_us = offset_ns[i] / 1000.; + + nu.time += delta_us; + + for (SRTrueParticle &part: nu.prim) { + if (part.genT > -9998.) part.genT += delta_us; + if (part.startT > -9998.) part.startT += delta_us; + if (part.endT > -9998.) part.endT += delta_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_ns[slc.tmatch.index])) continue; + + double const off_ns = offset_ns[slc.tmatch.index]; + double const delta_us = off_ns / 1000.; + + // slc.truth is a copy of rec.mc.nu[tmatch.index] taken earlier (FillTrue.cxx:372). + slc.truth.time += delta_us; + for (SRTrueParticle &part: slc.truth.prim) { + if (part.genT > -9998.) part.genT += delta_us; + if (part.startT > -9998.) part.startT += delta_us; + if (part.endT > -9998.) part.endT += delta_us; + } + + // Opt0 FM + if (!std::isnan(slc.opt0.time)) slc.opt0.time += delta_us; // NaN default + if (!std::isnan(slc.opt0_sec.time)) slc.opt0_sec.time += delta_us; // NaN default + + // BFM + if (slc.barycenterFM.flashTime != -9999.) slc.barycenterFM.flashTime += delta_us; // -9999. default + if (slc.barycenterFM.flashFirstHit != -9999.) slc.barycenterFM.flashFirstHit += delta_us; // -9999. default + + // BFM + Light Propagation Module + if (slc.correctedOpFlash.OpFlashT0 != -9999.) { // -9999. default + slc.correctedOpFlash.OpFlashT0 += delta_us; + slc.correctedOpFlash.OpFlashT0Corrected += delta_us; + } + + // Simple FM: PMT & XARAPUCA + if (slc.fmatch.time != -9999.) slc.fmatch.time += delta_us; // -9999. default + if (slc.fmatchop.time != -9999.) slc.fmatchop.time += delta_us; // -9999. default + if (slc.fmatchara.time != -9999.) slc.fmatchara.time += delta_us; // -9999. default + if (slc.fmatchopara.time != -9999.) slc.fmatchopara.time += delta_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 += off_ns; + if (!std::isnan(pfp.trk.crtsbndtrack.score)) pfp.trk.crtsbndtrack.track.time += off_ns; + } + + // CRUMBS BDT + if (!std::isnan(slc.crumbs_result.pds.fmtime)) slc.crumbs_result.pds.fmtime += delta_us; // NaN default + if (!std::isnan(slc.crumbs_result.crt.tracktime)) slc.crumbs_result.crt.tracktime += delta_us; // NaN default + if (!std::isnan(slc.crumbs_result.crt.sptime)) slc.crumbs_result.crt.sptime += delta_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 +649,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) { @@ -1454,6 +1562,14 @@ void CAFMaker::produce(art::Event& evt) noexcept { art::fill_ptr_vector(mcfluxes, mcflux_handle); } + art::Handle> dk2nu_handle; + GetByLabelStrict(evt, std::string("generator"), dk2nu_handle); + + std::vector> dk2nus; + if (dk2nu_handle.isValid()) { + art::fill_ptr_vector(dk2nus, dk2nu_handle); + } + // get the MCReco for the fake-reco art::Handle> mctrack_handle; GetByLabelStrict(evt, std::string("mcreco"), mctrack_handle); @@ -1537,9 +1653,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 baddk2nu; + for (size_t i=0; i v10_20_09 if (isRealData && (fDet == kSBND) && fSubRunPOT > 0) { + std::cout << "CAFMaker: run " << evt.run() << " subRun " << evt.subRun() + << " fSubRunPOT=" << fSubRunPOT + << " frameApplyAtCaf=" << rec.sbnd_frames.frameApplyAtCaf << std::endl; + // Fill trigger info FillTriggerSBND(srsbndtiminginfo, srtrigger); @@ -2709,11 +2846,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..8be94a3f6 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -1,5 +1,7 @@ #include "FillTrue.h" +#include "dk2nu/tree/dk2nu.h" + #include "larcorealg/GeoAlgo/GeoAlgo.h" #include "larcorealg/Geometry/GeometryCore.h" #include "larcorealg/Geometry/WireReadoutGeom.h" @@ -385,6 +387,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 +460,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, @@ -549,6 +553,9 @@ namespace caf { // Set the MCFlux stuff srneutrino.initpdg = mcflux.fntype; srneutrino.baseline = mcflux.fdk2gen + mcflux.fgen2vtx; + srneutrino.dk2gen = mcflux.fdk2gen; + srneutrino.prod_time = (!dk2nu.ancestor.empty()) + ? dk2nu.ancestor.back().startt : 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,