diff --git a/sbncode/CAFMaker/CAFMakerParams.h b/sbncode/CAFMaker/CAFMakerParams.h index 28d2d6d66..544ec72f4 100644 --- a/sbncode/CAFMaker/CAFMakerParams.h +++ b/sbncode/CAFMaker/CAFMakerParams.h @@ -320,6 +320,12 @@ namespace caf "crthit" // icarus }; + Atom CRTSimChanLabel { + Name("CRTSimChanLabel"), + Comment("Label of AuxDetSimChannels."), + "genericcrt" // icarus + }; + Atom CRTTrackLabel { Name("CRTTrackLabel"), Comment("Label of sbn CRT tracks."), diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index a2edf0349..985b8e127 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -140,6 +140,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" @@ -287,6 +288,11 @@ class CAFMaker : public art::EDProducer { std::vector> fTPCVolumes; std::vector fActiveVolumes; + // CRT geometry info + // + // ICARUS + std::map, int> fFEBChannel2AuxDetID; + // random number generator for fake reco CLHEP::HepRandomEngine& fFakeRecoRandomEngine; @@ -315,6 +321,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); @@ -423,6 +430,9 @@ class CAFMaker : public art::EDProducer { // setup volume definitions InitVolumes(); + // setup CRT mapping + InitCRTMapping(); + fSaveGENIEEventRecord = fParams.SaveGENIEEventRecord(); } @@ -596,6 +606,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 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(); @@ -1649,16 +1689,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, sim::AuxDetSimChannel> crtsimchanmap; + if(fDet == kICARUS) { art::Handle> crthits_handle; GetByLabelStrict(evt, fParams.CRTHitLabel(), crthits_handle); + + art::Handle> auxdetsimchan_handle; + GetByLabelStrict(evt, fParams.CRTSimChanLabel(), auxdetsimchan_handle); + // fill into event if (crthits_handle.isValid()) { const std::vector &crthits = *crthits_handle; + + std::vector empty; + const std::vector &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()); } } @@ -2318,7 +2378,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) { diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index a6eb34f1a..fe79f99e1 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/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, 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 true_energy; + std::set> 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 &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 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> 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, @@ -635,6 +682,7 @@ namespace caf bool use_ts0, int64_t CRT_T0_reference_time, // ns, signed double CRT_T1_reference_time, // us + const std::map, sim::AuxDetSimChannel> &crtsimchanmap, caf::SRTrack &srtrack, bool allowEmpty) { @@ -661,7 +709,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); } } diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h index 7eb7a7948..39bfd2318 100644 --- a/sbncode/CAFMaker/FillReco.h +++ b/sbncode/CAFMaker/FillReco.h @@ -28,6 +28,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" @@ -177,6 +178,7 @@ namespace caf bool use_ts0, int64_t CRT_T0_reference_time, // ns, signed double CRT_T1_reference_time, // us + const std::map, sim::AuxDetSimChannel> &crtsimchanmap, caf::SRTrack &srtrack, bool allowEmpty = false); @@ -250,6 +252,7 @@ namespace caf bool use_ts0, int64_t CRT_T0_reference_time, // ns, signed double CRT_T1_reference_time, // us + const std::map, sim::AuxDetSimChannel> &crtsimchanmap, caf::SRCRTHit &srhit, bool allowEmpty = false); void FillCRTTrack(const sbn::crt::CRTTrack &track,