From 94abed0d493e9d354f98dcaa3c778529dda5ed3a Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Mon, 31 Mar 2025 14:58:47 -0500 Subject: [PATCH 01/24] first commit of NuGraphIcarus --- CMakeLists.txt | 1 + icaruscode/CMakeLists.txt | 1 + icaruscode/NuGraphIcarus/CMakeLists.txt | 63 ++ .../NuGraphIcarus/IcarusHDF5Maker_module.cc | 719 ++++++++++++++++++ .../NuGraphIcarus/IcarusNuGraphLoader_tool.cc | 160 ++++ .../IcarusNuSliceHitsProducer_module.cc | 173 +++++ .../IcarusTrueNuSliceHitsProducer_module.cc | 200 +++++ .../NuGraphIcarus/cafmaker_nugraph_test.fcl | 6 + .../NuGraphIcarus/hdf5maker_icarus_slice.fcl | 92 +++ .../testinference_slice_icarus.fcl | 115 +++ 10 files changed, 1530 insertions(+) create mode 100644 icaruscode/NuGraphIcarus/CMakeLists.txt create mode 100644 icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc create mode 100644 icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc create mode 100644 icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc create mode 100644 icaruscode/NuGraphIcarus/IcarusTrueNuSliceHitsProducer_module.cc create mode 100644 icaruscode/NuGraphIcarus/cafmaker_nugraph_test.fcl create mode 100644 icaruscode/NuGraphIcarus/hdf5maker_icarus_slice.fcl create mode 100644 icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl diff --git a/CMakeLists.txt b/CMakeLists.txt index a556fdd8a..f0c717daa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -66,6 +66,7 @@ find_package(lardata REQUIRED ) find_package(larsim REQUIRED ) find_package(larevt REQUIRED ) find_package(larreco REQUIRED ) +find_package(larrecodnn REQUIRED ) find_package(larana REQUIRED ) find_package(larpandora REQUIRED ) find_package(larpandoracontent REQUIRED ) diff --git a/icaruscode/CMakeLists.txt b/icaruscode/CMakeLists.txt index a87cbffd7..f603b9238 100644 --- a/icaruscode/CMakeLists.txt +++ b/icaruscode/CMakeLists.txt @@ -10,6 +10,7 @@ add_subdirectory(Filters) add_subdirectory(Generators) add_subdirectory(Geometry) add_subdirectory(LArG4) +add_subdirectory(NuGraphIcarus) add_subdirectory(PMT) add_subdirectory(RecoUtils) add_subdirectory(TPC) diff --git a/icaruscode/NuGraphIcarus/CMakeLists.txt b/icaruscode/NuGraphIcarus/CMakeLists.txt new file mode 100644 index 000000000..4c269b433 --- /dev/null +++ b/icaruscode/NuGraphIcarus/CMakeLists.txt @@ -0,0 +1,63 @@ +cet_enable_asserts() + +set( nugraph_tool_lib_list + lardataobj::RecoBase + lardataobj::AnalysisBase + larrecodnn::NuGraphBaseTools + TorchScatter::TorchScatter + art::Framework_Core + hep_concurrency::hep_concurrency + art_plugin_types::tool +) + +include_directories($ENV{HEP_HPC_INC}) + +cet_build_plugin(IcarusNuGraphLoader art::tool + LIBRARIES PRIVATE + ${nugraph_tool_lib_list} +) + +simple_plugin(IcarusNuSliceHitsProducer "module" + LIBRARIES PRIVATE + art::Framework_Core + nusimdata::SimulationBase + lardataobj::RecoBase + lardataalg::DetectorInfo + lardata::DetectorInfoServices_DetectorClocksServiceStandard_service + larcore::Geometry_Geometry_service + sbnobj::Common_Reco +) + +simple_plugin(IcarusTrueNuSliceHitsProducer "module" + LIBRARIES PRIVATE + art::Framework_Core + nusimdata::SimulationBase + lardataobj::RecoBase + lardataalg::DetectorInfo + lardata::DetectorInfoServices_DetectorClocksServiceStandard_service + larcore::Geometry_Geometry_service + sbnobj::Common_Reco + larsim::MCCheater_BackTrackerService_service + larsim::MCCheater_ParticleInventoryService_service +) + +simple_plugin(IcarusHDF5Maker "module" + LIBRARIES PRIVATE + art::Framework_Core + nusimdata::SimulationBase + larsim::MCCheater_BackTrackerService_service + larsim::MCCheater_ParticleInventoryService_service + hep_hpc_hdf5 + #hep::hpc_Utilities + lardata::RecoBaseProxy + lardataobj::RecoBase + lardataalg::DetectorInfo + lardata::DetectorInfoServices_DetectorClocksServiceStandard_service + larcore::Geometry_Geometry_service + ${HDF5_LIBRARIES} +) + +install_headers() +install_fhicl() +install_source() + diff --git a/icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc b/icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc new file mode 100644 index 000000000..a2eeb0515 --- /dev/null +++ b/icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc @@ -0,0 +1,719 @@ +//////////////////////////////////////////////////////////////////////// +//// Class: IcarusHDF5Maker +//// Plugin Type: analyzer (art v3_06_03) +//// File: IcarusHDF5Maker_module.cc +//// +//// Generated at Wed May 5 08:23:31 2021 by V Hewes using cetskelgen +/// from cetlib version v3_11_01. +////////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "larcore/Geometry/Geometry.h" +#include "larcorealg/Geometry/GeometryCore.h" +#include "larcorealg/Geometry/TPCGeo.h" +#include "larcore/Geometry/WireReadout.h" + +#include "nusimdata/SimulationBase/MCTruth.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/SpacePoint.h" +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/RecoBase/OpHit.h" +#include "lardataobj/RecoBase/OpFlash.h" +#include "lardata/RecoBaseProxy/ProxyBase.h" +#include "lardataobj/RecoBase/Wire.h" +#include "larcorealg/Geometry/Exceptions.h" +#include "lardataobj/RecoBase/Slice.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardata/RecoBaseProxy/ProxyBase.h" +#include "lardataobj/RecoBase/Vertex.h" +#include "lardataobj/RecoBase/Cluster.h" +#include "lardataobj/AnalysisBase/T0.h" +#include "lardataobj/RecoBase/PFParticleMetadata.h" + +#include "larsim/MCCheater/BackTrackerService.h" +#include "larsim/MCCheater/ParticleInventoryService.h" +#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "lardata/DetectorInfoServices/DetectorClocksService.h" + +#include "lardataobj/AnalysisBase/BackTrackerMatchingData.h" + +#include "hep_hpc/hdf5/make_ntuple.hpp" +#include "art/Framework/Services/Registry/ServiceHandle.h" + +class IcarusHDF5Maker : public art::EDAnalyzer { +public: + explicit IcarusHDF5Maker(fhicl::ParameterSet const& p); + ~IcarusHDF5Maker() noexcept {}; // bare pointers are cleaned up by endSubRun + + IcarusHDF5Maker(IcarusHDF5Maker const&) = delete; + IcarusHDF5Maker(IcarusHDF5Maker&&) = delete; + IcarusHDF5Maker& operator=(IcarusHDF5Maker const&) = delete; + IcarusHDF5Maker& operator=(IcarusHDF5Maker&&) = delete; + + void beginSubRun(art::SubRun const& sr) override; + void endSubRun(art::SubRun const& sr) override; + void analyze(art::Event const& e) override; + +private: + + std::string fTruthLabel; + std::string fHitLabel; + std::string fHitTruthLabel; + std::string fSPLabel; + std::string fOpHitLabel; + std::string fOpFlashLabel; + std::string fPFParticleLabel; + + bool fUseMap; + std::string fEventInfo; + std::string fOutputName; + + hep_hpc::hdf5::File fFile; ///< output HDF5 file + + hep_hpc::hdf5::Ntuple // event id (run, subrun, event) + >* fEventNtuple; ///< event ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // is cc + hep_hpc::hdf5::Column, // nu pdg + hep_hpc::hdf5::Column, // nu energy + hep_hpc::hdf5::Column, // lep energy + hep_hpc::hdf5::Column // nu dir (x, y, z) + >* fEventNtupleNu; ///< event ntuple with neutrino information + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // spacepoint id + hep_hpc::hdf5::Column, // 3d position (x, y, z) + hep_hpc::hdf5::Column, // 2d hit (u, v, y) + hep_hpc::hdf5::Column //ChiSquared of the hit + >* fSpacePointNtuple; ///< spacepoint ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // hit integral + hep_hpc::hdf5::Column, // hit rms + hep_hpc::hdf5::Column, // tpc id + hep_hpc::hdf5::Column, // global plane + hep_hpc::hdf5::Column, // global wire + hep_hpc::hdf5::Column, // global time + hep_hpc::hdf5::Column, // raw plane + hep_hpc::hdf5::Column, // raw wire + hep_hpc::hdf5::Column, // raw time + hep_hpc::hdf5::Column //cryostat + >* fHitNtuple; ///< hit ntuple + + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // g4 id + hep_hpc::hdf5::Column, // particle type + hep_hpc::hdf5::Column, // parent g4 id + hep_hpc::hdf5::Column, // is from nu + hep_hpc::hdf5::Column, // momentum + hep_hpc::hdf5::Column, // start position (x, y, z) + hep_hpc::hdf5::Column, // end position (x, y, z) + hep_hpc::hdf5::Column, // start process + hep_hpc::hdf5::Column // end process + >* fParticleNtuple; ///< particle ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // g4 id + hep_hpc::hdf5::Column, // deposited energy [ MeV ] + hep_hpc::hdf5::Column, // x position + hep_hpc::hdf5::Column, // y position + hep_hpc::hdf5::Column // z position + >* fEnergyDepNtuple; ///< energy deposition ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // hit_channel + hep_hpc::hdf5::Column, // wire pos + hep_hpc::hdf5::Column, // peaktime + hep_hpc::hdf5::Column, // width + hep_hpc::hdf5::Column, // area + hep_hpc::hdf5::Column, // amplitude + hep_hpc::hdf5::Column, // pe + hep_hpc::hdf5::Column // usedInFlash + >* fOpHitNtuple; ///< PMT hit ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // flash id + hep_hpc::hdf5::Column, // wire pos + hep_hpc::hdf5::Column, // time + hep_hpc::hdf5::Column, // time width + hep_hpc::hdf5::Column, // Y center + hep_hpc::hdf5::Column, // Y width + hep_hpc::hdf5::Column, // Z center + hep_hpc::hdf5::Column, // Z width + hep_hpc::hdf5::Column // totalpe + >* fOpFlashNtuple; ///< Flash ntuple + +hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // sumpe id + hep_hpc::hdf5::Column, // flash id + hep_hpc::hdf5::Column, // PMT channel + hep_hpc::hdf5::Column // pe + >* fOpFlashSumPENtuple; ///< Flash SumPE ntuple + + + + + + using ProxyPfpColl_t = decltype(proxy::getCollection >( + std::declval(),std::declval(), + proxy::withAssociated(std::declval()), + proxy::withAssociated(std::declval()), + proxy::withAssociated(std::declval()), + proxy::withAssociated(std::declval()) )); + using ProxyPfpElem_t = ProxyPfpColl_t::element_proxy_t; + + // proxy to connect cluster to hit + using ProxyClusColl_t = decltype(proxy::getCollection>( + std::declval(), std::declval(), + proxy::withAssociated(std::declval()))); + using ProxyClusElem_t = ProxyClusColl_t::element_proxy_t; + + std::vector nearwires; + int NearWire(const geo::WireReadoutGeom& geo, const geo::PlaneID &ip, const float x, const float y, const float z); + + int stitchedPlane(const geo::WireID wid) const { + int plane = wid.Plane; + if(wid.TPC==2 || wid.TPC==3) { + if(wid.Plane==1) plane=2; + else if(wid.Plane==2) plane=1; + } + return plane; + } + float stitchedTime(const geo::WireID wid, float timein) const { + float time = timein; + if(wid.TPC==2 || wid.TPC==3) { + //correction = 2*(tpcgeo.DriftDistance()/detProp.DriftVelocity()-clockData.TriggerOffsetTPC())/clockData.TPCClock().TickPeriod() = 6442.15 + time = 6442.15 - timein; + } + return time; + } + size_t stitchedWire(const geo::WireID wid) const { + size_t wire = wid.Wire; + int plane = stitchedPlane(wid); + if(wid.TPC==1 || wid.TPC==3) { + if(plane==1 || plane == 2) { + wire = wid.Wire + 2535; //2535 is the last part of the wires in cryos 0 an 2 before the cut in z=0 + } + } + return wire; + } + +}; + +IcarusHDF5Maker::IcarusHDF5Maker(fhicl::ParameterSet const& p) + : EDAnalyzer{p}, + fTruthLabel(p.get("TruthLabel")), + fHitLabel( p.get("HitLabel")), + fHitTruthLabel( p.get("HitTruthLabel","")), + fSPLabel( p.get("SPLabel")), + fOpHitLabel( p.get("OpHitLabel")), + fOpFlashLabel( p.get("OpFlashLabel")), + fPFParticleLabel( p.get("PFParticleLabel")), + fUseMap( p.get("UseMap", false)), + fEventInfo( p.get("EventInfo")), + fOutputName(p.get("OutputName")) +{ + if (fEventInfo != "none" && fEventInfo != "nu") + throw art::Exception(art::errors::Configuration) + << "EventInfo must be \"none\" or \"nu\", not " << fEventInfo; +} +std::vector tpc_ids_checked = {}; // add the TPC ids here so we only check them once + +void IcarusHDF5Maker::analyze(art::Event const& e) { + + const cheat::BackTrackerService* bt = 0; + if (!fUseMap) { + art::ServiceHandle bt_h; + bt = bt_h.get(); + } + geo::WireReadoutGeom const& geo = art::ServiceHandle()->Get(); + + std::set g4id; + art::ServiceHandle pi; + auto const clockData = art::ServiceHandle()->DataFor(e); + auto const detProp = art::ServiceHandle()->DataFor(e, clockData); + int run = e.id().run(); + int subrun = e.id().subRun(); + int event = e.id().event(); + + // auto tpcgeo = art::ServiceHandle()->TPC(); + // std::cout << "tpc drift distance=" << tpcgeo.DriftDistance() << std::endl; + // std::cout << "drift velocity=" << detProp.DriftVelocity() << std::endl; + // std::cout << "drift time=" << tpcgeo.DriftDistance()/detProp.DriftVelocity() << " 2*dt/0.4=" << 2*tpcgeo.DriftDistance()/detProp.DriftVelocity()/clockData.TPCClock().TickPeriod() << std::endl; + // std::cout << "tick period=" << clockData.TPCClock().TickPeriod() << std::endl; + // std::cout << "tpcTime=" << clockData.TPCTime() << std::endl; + // std::cout << "triggerOffsetTPC=" << clockData.TriggerOffsetTPC() << std::endl; + // std::cout << "triggerTime=" << clockData.TriggerTime() << std::endl; + // std::cout << "correction=" << 2*(tpcgeo.DriftDistance()/detProp.DriftVelocity()-clockData.TriggerOffsetTPC())/clockData.TPCClock().TickPeriod() << std::endl; + + std::cout << "NEW EVENT: " << run <<" "<< subrun << " "<< event< evtID { run, subrun, event }; + // Fill event table + if (fEventInfo == "none") { + fEventNtuple->insert( evtID.data() ); + mf::LogInfo("IcarusHDF5Maker") << "Filling event table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2]; + } + + if (fEventInfo == "nu") { + // Get MC truth + art::Handle> truthHandle; + e.getByLabel(fTruthLabel, truthHandle); + if (!truthHandle.isValid() || truthHandle->size() == 0) { + throw art::Exception(art::errors::LogicError) + << "Expected to find exactly one MC truth object!"; + } + simb::MCNeutrino nutruth = truthHandle->at(0).GetNeutrino(); + + auto up = nutruth.Nu().Momentum().Vect().Unit(); + std::array nuMomentum {(float)up.X(),(float)up.Y(),(float)up.Z()}; + + fEventNtupleNu->insert( evtID.data(), + nutruth.CCNC() == simb::kCC, + nutruth.Nu().PdgCode(), + nutruth.Nu().E(), + nutruth.Lepton().E(), + nuMomentum.data() + ); + + // for (int ip=0;ipat(0).NParticles();ip++) { + // std::cout << "mcp tkid=" << truthHandle->at(0).GetParticle(ip).TrackId() << " pdg=" << truthHandle->at(0).GetParticle(ip).PdgCode() + // << " mother=" << truthHandle->at(0).GetParticle(ip).Mother() + // << " vtx=" << truthHandle->at(0).GetParticle(ip).Vx() << " " << truthHandle->at(0).GetParticle(ip).Vy() << " " << truthHandle->at(0).GetParticle(ip).Vz() + // << std::endl; + // } + + mf::LogInfo("IcarusHDF5Maker") << "Filling event table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nis cc? " << (nutruth.CCNC() == simb::kCC) + << ", nu energy " << nutruth.Nu().E() + << ", lepton energy " << nutruth.Lepton().E() + << "\nnu momentum x " << nuMomentum[0] << ", y " + << nuMomentum[1] << ", z " << nuMomentum[2]; + } // if nu event info + + // Get spacepoints from the event record + art::Handle> spListHandle; + std::vector> splist; + if (e.getByLabel(fSPLabel, spListHandle)) art::fill_ptr_vector(splist, spListHandle); + + // Get hits from the event record + art::Handle> hitListHandle; + std::vector> hitlist; + if (e.getByLabel(fHitLabel, hitListHandle)) art::fill_ptr_vector(hitlist, hitListHandle); + + // Get assocations from spacepoints to hits + art::FindManyP fmp(spListHandle, e, fSPLabel); + + // Fill spacepoint table + for (size_t i = 0; i < spListHandle->size(); ++i) { + + art::Ptr spacePointPtr(spListHandle,i); + std::vector> associatedHits(fmp.at(spacePointPtr.key())); + if (associatedHits.size() < 3) { + std::cout << "I am certain this cannot happen... but here you go, space point with " << associatedHits.size() << " hits" << std::endl; + exit(1); + continue; + } + + std::array pos {(float)splist[i]->XYZ()[0],(float)splist[i]->XYZ()[1],(float)splist[i]->XYZ()[2]}; + + std::array hitID { -1, -1, -1 }; + for (size_t j = 0; j < associatedHits.size(); ++j) { + int plane = stitchedPlane(associatedHits[j]->WireID()); + //std::cout << "j=" << j << " tpc=" << associatedHits[j]->WireID().TPC << " plane=" << associatedHits[j]->WireID().Plane + // << " plane2=" << plane << " key=" << associatedHits[j].key() << std::endl; + hitID[plane] = associatedHits[j].key(); + } + if (hitID[0]==0 && hitID[1]==-1 && hitID[2]==-1) { + std::cout << "THIS SHOULD NOT HAPPEN -- sp i=" << i << " hit ids=" << hitID[0] << " " << hitID[1] << " " << hitID[2] << std::endl; + exit(1); + } + fSpacePointNtuple->insert(evtID.data(),splist[i]->ID(), pos.data(), hitID.data(),splist[i]->Chisq() ); + + mf::LogInfo("IcarusHDF5Maker") << "Filling spacepoint table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nspacepoint id " << splist[i]->ID() + << "\nposition x " << pos[0] << ", y " << pos[1] + << ", z " << pos[2] + << "\nhit ids " << hitID[0] << ", " << hitID[1] + << ", " << hitID[2] + << "Chi_squared " << splist[i]->Chisq(); + } // for spacepoint + + std::unique_ptr> hittruth; + if (fUseMap) { + hittruth = std::unique_ptr >(new art::FindManyP(hitListHandle, e, fHitTruthLabel)); + } + + // Loop over hits + for (art::Ptr hit : hitlist) { + + // Fill hit table + geo::WireID wireid = hit->WireID(); + size_t plane = stitchedPlane(wireid); + double time = stitchedTime(wireid,hit->PeakTime()); + size_t wire = stitchedWire(wireid); + fHitNtuple->insert(evtID.data(), + hit.key(), hit->Integral(), hit->RMS(), wireid.TPC, + plane, wire, time, + wireid.Plane, wireid.Wire, hit->PeakTime(),wireid.Cryostat + ); + + mf::LogInfo("IcarusHDF5Maker") << "Filling hit table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nhit id " << hit.key() << ", integral " + << hit->Integral() << ", RMS " << hit->RMS() + << ", TPC " << wireid.TPC + << "\nglobal plane " << plane << ", global wire " + << wire << ", global time " << time + << "\nlocal plane " << wireid.Plane + << ", local wire " << wireid.Wire + << ", local time " << hit->PeakTime() + << ", Cryostat " << wireid.Cryostat; + + // Fill energy deposit table + if (fUseMap) { + //not supported in icarus at the moment + exit(1); + /* + std::vector> particle_vec = hittruth->at(hit.key()); + std::vector match_vec = hittruth->data(hit.key()); + //loop over particles + for (size_t i_p = 0; i_p < particle_vec.size(); ++i_p) { + g4id.insert(particle_vec[i_p]->TrackId()); + fEnergyDepNtuple->insert(evtID.data(), hit.key(), + particle_vec[i_p]->TrackId(), + match_vec[i_p]->ideFraction, + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()); + mf::LogInfo("IcarusHDF5Maker") << "Filling energy deposit table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nhit id " << hit.key() << ", g4 id " + << particle_vec[i_p]->TrackId() << ", energy fraction " + << match_vec[i_p]->ideFraction << ", position (" + << std::nan << ", " << std::nan << ", " << std::nan << ")";; + } // for particle + */ + } else { + + std::vector h_ides = bt->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick()); + for (auto& tide : h_ides) { + int tid = tide.trackID; + if (pi->TrackIdToParticle_P(abs(tid))) tid = abs(tid); + g4id.insert(tid); + fEnergyDepNtuple->insert(evtID.data(),hit.key(), tid, tide.energyFrac, -1., -1., -1.); + } + + } // if using microboone map method or not + } // for hit + + std::set allIDs = g4id; // Copy original so we can safely modify it + + // Add invisible particles to hierarchy + for (int id : g4id) { + const simb::MCParticle* p = pi->TrackIdToParticle_P(abs(id)); + while (p!= NULL && p->Mother() != 0 ) { + allIDs.insert(abs(p->Mother())); + p = pi->TrackIdToParticle_P(abs(p->Mother())); + } + } + // Loop over true particles and fill table + for (int id : allIDs) { + const simb::MCParticle* p = pi->TrackIdToParticle_P(abs(id)); + if(p==NULL) std::cout<< "we have a problem p is null"<TrackIdToMCTruth_P(abs(id)); + int fromNu = (mct.isAvailable() ? mct->NeutrinoSet() : 0); + // std::cout << "all mcp tkid=" << p->TrackId() << " pdg=" << p->PdgCode() + // << " mother=" << p->Mother() + // << " vtx=" << p->Vx() << " " << p->Vy() << " " << p->Vz() + // << " mctruth=" << mct->NeutrinoSet() + // << std::endl; + std::array particleStart { (float)p->Vx(), (float)p->Vy(), (float)p->Vz() }; + std::array particleEnd { (float)p->EndX(), (float)p->EndY(), (float)p->EndZ() }; + fParticleNtuple->insert(evtID.data(), + abs(id), p->PdgCode(), p->Mother(), fromNu, (float)p->P(), + particleStart.data(), particleEnd.data(), + p->Process(), p->EndProcess() + ); + mf::LogInfo("IcarusHDF5Maker") << "Filling particle table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\ng4 id " << abs(id) << ", pdg code " + << p->PdgCode() << ", parent " << p->Mother() + << ", momentum " << p->P() + << "\nparticle start x " << particleStart[0] + << ", y " << particleStart[1] + << ", z " << particleStart[2] + << "\nparticle end x " << particleEnd[0] << ", y " + << particleEnd[1] << ", z " << particleEnd[2] + << "\nstart process " << p->Process() + << ", end process " << p->EndProcess(); + } + + // Get OpFlashs from the event record + art::Handle< std::vector< recob::OpFlash > > opFlashListHandle; + std::unique_ptr > assocFlashHit; + std::vector< art::Ptr< recob::OpFlash > > opflashlist; + if (fOpFlashLabel != "") { + if (e.getByLabel(fOpFlashLabel, opFlashListHandle)) art::fill_ptr_vector(opflashlist, opFlashListHandle); + assocFlashHit = std::unique_ptr >(new art::FindManyP(opFlashListHandle, e, fOpFlashLabel)); + } + + // get the flash matching the slice we selected + art::Handle< art::Assns > sliceOpFlashAssnsHandle; + e.getByLabel(fHitLabel, sliceOpFlashAssnsHandle); + int flkey = -1; + if (sliceOpFlashAssnsHandle->size()) { + art::Ptr fp = sliceOpFlashAssnsHandle->at(0).second; + flkey = fp.key(); + } + + std::vector > flashsumpepmtmap;//this has at most one element, due to the check on flkey + + // Loop over opflashs, pick only the one we care about + for (art::Ptr< recob::OpFlash > opflash : opflashlist) { + if (int(opflash.key())!=flkey) continue; + std::vector pes; + for (auto pe : opflash->PEs()) pes.push_back(pe); + std::vector nearwires; + double xyz[3] = {0.,opflash->YCenter(),opflash->ZCenter()}; + for (auto const& p : geo.Iterate()) nearwires.push_back(NearWire(geo,p,xyz[0],xyz[1],xyz[2])); + + // Fill opflash table + fOpFlashNtuple->insert(evtID.data(), + opflash.key(), + nearwires.data(), + opflash->Time(),opflash->TimeWidth(), + opflash->YCenter(),opflash->YWidth(),opflash->ZCenter(),opflash->ZWidth(), + opflash->TotalPE() + ); + size_t count = 0; + std::vector sumpepmtmap(art::ServiceHandle()->NOpDets(),-1); + for (size_t ipmt=0;ipmtinsert(evtID.data(),count,opflash.key(),ipmt,pes[ipmt]); + sumpepmtmap[ipmt] = count; + count++; + } + flashsumpepmtmap.push_back(sumpepmtmap); + mf::LogInfo("IcarusHDF5Maker") << "Filling opflash table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nflash id " << opflash.key() << ", Time " << opflash->Time() + << ", TotalPE " << opflash->TotalPE()//; + << "\nWireCenters size0 " << opflash->WireCenters().size()//; + << "\nYCenter " << opflash->YCenter()<< " ZCenter " << opflash->ZCenter() + << "\nYWidth " << opflash->YWidth()<< " ZWidth " << opflash->ZWidth() + << "\nInBeamFrame " << opflash->InBeamFrame()<< " OnBeamTime " << opflash->OnBeamTime() + << "\nAbsTime " << opflash->AbsTime() << " TimeWidth " << opflash->TimeWidth() << " FastToTotal " << opflash->FastToTotal(); + } + + // Get OpHits from the event record + art::Handle< std::vector< recob::OpHit > > opHitListHandle; + std::vector< art::Ptr< recob::OpHit > > ophitlist; + if (fOpHitLabel != "") { + if (e.getByLabel(fOpHitLabel, opHitListHandle)) + art::fill_ptr_vector(ophitlist, opHitListHandle); + } + + // Loop over ophits + for (art::Ptr< recob::OpHit > ophit : ophitlist) { + std::vector nearwires; + auto xyz = geo.OpDetGeoFromOpChannel(ophit->OpChannel()).GetCenter(); + for (auto const& p : geo.Iterate()) nearwires.push_back(NearWire(geo,p,xyz.X(),xyz.Y(),xyz.Z())); + + bool isInFlash = false; + if (flkey>=0) { + auto ophv = assocFlashHit->at(flkey); + isInFlash = (std::find(ophv.begin(),ophv.end(),ophit)!=ophv.end()); + } + // Fill ophit table + fOpHitNtuple->insert(evtID.data(),ophit.key(), + ophit->OpChannel(),nearwires.data(), + ophit->PeakTime(),ophit->Width(), + ophit->Area(),ophit->Amplitude(),ophit->PE(), + (isInFlash ? flashsumpepmtmap[0][ophit->OpChannel()] : -1) + ); + mf::LogInfo("IcarusHDF5Maker") << "\nFilling ophit table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nhit id " << ophit.key() << ", channel " + << ophit->OpChannel() << ", PeakTime " << ophit->PeakTime() + << ", Width " << ophit->Width() + << "\narea " << ophit->Area() << ", amplitude " + << ophit->Amplitude() << ", PE " << ophit->PE(); + } + //End optical analyzer + +} // function IcarusHDF5Maker::analyze + +void IcarusHDF5Maker::beginSubRun(art::SubRun const& sr) { + + struct timeval now; + gettimeofday(&now, NULL); + // Open HDF5 output + std::ostringstream fileName; + fileName << fOutputName << "_r" << std::setfill('0') << std::setw(5) << sr.run() + << "_s" << std::setfill('0') << std::setw(5) << sr.subRun() << "_ts" << std::setw(6) << now.tv_usec << ".h5"; + + fFile = hep_hpc::hdf5::File(fileName.str(), H5F_ACC_TRUNC); + + if (fEventInfo == "none") + fEventNtuple = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "event_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3) + )); + if (fEventInfo == "nu") + fEventNtupleNu = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "event_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("is_cc"), + hep_hpc::hdf5::make_scalar_column("nu_pdg"), + hep_hpc::hdf5::make_scalar_column("nu_energy"), + hep_hpc::hdf5::make_scalar_column("lep_energy"), + hep_hpc::hdf5::make_column("nu_dir", 3) + )); + + fSpacePointNtuple = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "spacepoint_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("spacepoint_id"), + hep_hpc::hdf5::make_column("position", 3), + hep_hpc::hdf5::make_column("hit_id", 3), + hep_hpc::hdf5::make_scalar_column("Chi_squared") + + )); + + fHitNtuple = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "hit_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("hit_id"), + hep_hpc::hdf5::make_scalar_column("integral"), + hep_hpc::hdf5::make_scalar_column("rms"), + hep_hpc::hdf5::make_scalar_column("tpc"), + hep_hpc::hdf5::make_scalar_column("global_plane"), + hep_hpc::hdf5::make_scalar_column("global_wire"), + hep_hpc::hdf5::make_scalar_column("global_time"), + hep_hpc::hdf5::make_scalar_column("local_plane"), + hep_hpc::hdf5::make_scalar_column("local_wire"), + hep_hpc::hdf5::make_scalar_column("local_time"), + hep_hpc::hdf5::make_scalar_column("Cryostat") + + )); + + fParticleNtuple = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "particle_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("g4_id"), + hep_hpc::hdf5::make_scalar_column("type"), + hep_hpc::hdf5::make_scalar_column("parent_id"), + hep_hpc::hdf5::make_scalar_column("from_nu"), + hep_hpc::hdf5::make_scalar_column("momentum"), + hep_hpc::hdf5::make_column("start_position", 3), + hep_hpc::hdf5::make_column("end_position", 3), + hep_hpc::hdf5::make_scalar_column("start_process"), + hep_hpc::hdf5::make_scalar_column("end_process") + )); + + fEnergyDepNtuple = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "edep_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("hit_id"), + hep_hpc::hdf5::make_scalar_column("g4_id"), + hep_hpc::hdf5::make_scalar_column("energy"), + hep_hpc::hdf5::make_scalar_column("x_position"), + hep_hpc::hdf5::make_scalar_column("y_position"), + hep_hpc::hdf5::make_scalar_column("z_position") + )); + + fOpHitNtuple = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "ophit_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("hit_id"), + hep_hpc::hdf5::make_scalar_column("hit_channel"), + hep_hpc::hdf5::make_column("wire_pos", 3),//3 views + hep_hpc::hdf5::make_scalar_column("peaktime"), + hep_hpc::hdf5::make_scalar_column("width"), + hep_hpc::hdf5::make_scalar_column("area"), + hep_hpc::hdf5:: make_scalar_column("amplitude"), + hep_hpc::hdf5::make_scalar_column("pe"), + hep_hpc::hdf5::make_scalar_column("sumpe_id") + )); + + fOpFlashNtuple = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "opflash_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("flash_id"), + hep_hpc::hdf5::make_column("wire_pos", 3),//3 views + hep_hpc::hdf5::make_scalar_column("time"), + hep_hpc::hdf5::make_scalar_column("time_width"), + hep_hpc::hdf5::make_scalar_column("y_center"), + hep_hpc::hdf5::make_scalar_column("y_width"), + hep_hpc::hdf5::make_scalar_column("z_center"), + hep_hpc::hdf5::make_scalar_column("z_width"), + hep_hpc::hdf5::make_scalar_column("totalpe") + )); + + fOpFlashSumPENtuple = new hep_hpc::hdf5::Ntuple( + hep_hpc::hdf5::make_ntuple({fFile, "opflashsumpe_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("sumpe_id"), + hep_hpc::hdf5::make_scalar_column("flash_id"), + hep_hpc::hdf5::make_scalar_column("pmt_channel"), + hep_hpc::hdf5::make_scalar_column("sumpe") + )); +} + +void IcarusHDF5Maker::endSubRun(art::SubRun const& sr) { + if (fEventInfo == "none") delete fEventNtuple; + if (fEventInfo == "nu") delete fEventNtupleNu; + delete fSpacePointNtuple; + delete fHitNtuple; + delete fParticleNtuple; + delete fEnergyDepNtuple; + delete fOpHitNtuple; + delete fOpFlashNtuple; + delete fOpFlashSumPENtuple; + fFile.close(); +} + +int IcarusHDF5Maker::NearWire(const geo::WireReadoutGeom& geo, const geo::PlaneID &ip, const float x, const float y, const float z) +//art::ServiceHandle geo; +{ + geo::PlaneGeo const& plane = geo.Plane(ip); + geo::WireID wireID; + try { + wireID = plane.NearestWireID(geo::Point_t(x,y,z)); + } + catch (geo::InvalidWireError const& e) { + if (!e.hasSuggestedWire()) throw; + wireID = plane.ClosestWireID(e.suggestedWireID()); + } + return wireID.Wire; +} + +DEFINE_ART_MODULE(IcarusHDF5Maker) diff --git a/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc b/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc new file mode 100644 index 000000000..8062ba80d --- /dev/null +++ b/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc @@ -0,0 +1,160 @@ +#include "larrecodnn/NuGraph/Tools/LoaderToolBase.h" + +#include "art/Utilities/ToolMacros.h" + +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "canvas/Utilities/InputTag.h" +#include "lardataobj/RecoBase/SpacePoint.h" +#include + +class IcarusNuGraphLoader : public LoaderToolBase { + +public: + /** + * @brief Constructor + * + * @param pset + */ + IcarusNuGraphLoader(const fhicl::ParameterSet& pset); + + /** + * @brief Virtual Destructor + */ + virtual ~IcarusNuGraphLoader() noexcept = default; + + /** + * @brief loadData function + * + * @param art::Event event record, list of input, idsmap + */ + void loadData(art::Event& e, + vector>& hitlist, + vector& inputs, + vector>& idsmap) override; + + int stitchedPlane(const geo::WireID wid) const { + int plane = wid.Plane; + if(wid.TPC==2 || wid.TPC==3) { + if(wid.Plane==1) plane=2; + else if(wid.Plane==2) plane=1; + } + return plane; + } + float stitchedTime(const geo::WireID wid, float timein) const { + float time = timein; + if(wid.TPC==2 || wid.TPC==3) { + //correction = 2*(tpcgeo.DriftDistance()/detProp.DriftVelocity()-clockData.TriggerOffsetTPC())/clockData.TPCClock().TickPeriod() = 6442.15 + time = 6442.15 - timein; + } + return time; + } + size_t stitchedWire(const geo::WireID wid) const { + size_t wire = wid.Wire; + int plane = stitchedPlane(wid); + if(wid.TPC==1 || wid.TPC==3) { + if(plane==1 || plane == 2) { + wire = wid.Wire + 2535; //2535 is the last part of the wires in cryos 0 an 2 before the cut in z=0 + } + } + return wire; + } + +private: + art::InputTag hitInput; + art::InputTag spsInput; +}; + +IcarusNuGraphLoader::IcarusNuGraphLoader(const fhicl::ParameterSet& p) + : hitInput{p.get("hitInput")}, spsInput{p.get("spsInput")} +{} + +void IcarusNuGraphLoader::loadData(art::Event& e, + vector>& hitlist, + vector& inputs, + vector>& idsmap) +{ + // + art::Handle> hitListHandle; + if (e.getByLabel(hitInput, hitListHandle)) { art::fill_ptr_vector(hitlist, hitListHandle); } + // + idsmap = std::vector>(planes.size(), std::vector()); + for (auto h : hitlist) { + // + size_t plane = stitchedPlane(h->WireID()); + idsmap[plane].push_back(h.key()); + } + + vector hit_table_hit_id_data; + vector hit_table_local_plane_data; + vector hit_table_local_time_data; + vector hit_table_local_wire_data; + vector hit_table_integral_data; + vector hit_table_rms_data; + vector spacepoint_table_spacepoint_id_data; + vector spacepoint_table_hit_id_u_data; + vector spacepoint_table_hit_id_v_data; + vector spacepoint_table_hit_id_y_data; + + // hit table + for (auto h : hitlist) { + + geo::WireID wireid = h->WireID(); + size_t plane = stitchedPlane(wireid); + double time = stitchedTime(wireid,h->PeakTime()); + size_t wire = stitchedWire(wireid); + + hit_table_hit_id_data.push_back(h.key()); + hit_table_local_plane_data.push_back(plane); + hit_table_local_time_data.push_back(time); + hit_table_local_wire_data.push_back(wire); + hit_table_integral_data.push_back(h->Integral()); + hit_table_rms_data.push_back(h->RMS()); + } + std::cout << "loader has nhits=" << hit_table_hit_id_data.size() << std::endl; + + // Get spacepoints from the event record + art::Handle> spListHandle; + vector> splist; + if (e.getByLabel(spsInput, spListHandle)) { art::fill_ptr_vector(splist, spListHandle); } + // Get assocations from spacepoints to hits + vector>> sp2Hit(splist.size()); + if (splist.size() > 0) { + art::FindManyP fmp(spListHandle, e, spsInput); + for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) { + if (splist[spIdx]->Chisq()>0.5) continue; + if (fmp.at(spIdx).size()!=3) continue; + sp2Hit[spIdx] = fmp.at(spIdx); + } + } + + // space point table + for (size_t i = 0; i < splist.size(); ++i) { + if (sp2Hit[i].size()==0) continue; + spacepoint_table_spacepoint_id_data.push_back(i); + spacepoint_table_hit_id_u_data.push_back(-1); + spacepoint_table_hit_id_v_data.push_back(-1); + spacepoint_table_hit_id_y_data.push_back(-1); + for (size_t j = 0; j < sp2Hit[i].size(); ++j) { + // + size_t plane = stitchedPlane(sp2Hit[i][j]->WireID()); + if (plane == 0) spacepoint_table_hit_id_u_data.back() = sp2Hit[i][j].key(); + if (plane == 1) spacepoint_table_hit_id_v_data.back() = sp2Hit[i][j].key(); + if (plane == 2) spacepoint_table_hit_id_y_data.back() = sp2Hit[i][j].key(); + } + } + std::cout << "loader has nsps=" << spacepoint_table_hit_id_u_data.size() << std::endl; + + inputs.emplace_back("hit_table_hit_id", hit_table_hit_id_data); + inputs.emplace_back("hit_table_local_plane", hit_table_local_plane_data); + inputs.emplace_back("hit_table_local_time", hit_table_local_time_data); + inputs.emplace_back("hit_table_local_wire", hit_table_local_wire_data); + inputs.emplace_back("hit_table_integral", hit_table_integral_data); + inputs.emplace_back("hit_table_rms", hit_table_rms_data); + + inputs.emplace_back("spacepoint_table_spacepoint_id", spacepoint_table_spacepoint_id_data); + inputs.emplace_back("spacepoint_table_hit_id_u", spacepoint_table_hit_id_u_data); + inputs.emplace_back("spacepoint_table_hit_id_v", spacepoint_table_hit_id_v_data); + inputs.emplace_back("spacepoint_table_hit_id_y", spacepoint_table_hit_id_y_data); +} +DEFINE_ART_CLASS_TOOL(IcarusNuGraphLoader) diff --git a/icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc b/icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc new file mode 100644 index 000000000..7624ac35b --- /dev/null +++ b/icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc @@ -0,0 +1,173 @@ +//////////////////////////////////////////////////////////////////////// +// Class: IcarusNuSliceHitsProducer +// Plugin Type: producer (art v3_06_03) +// File: IcarusNuSliceHitsProducer_module.cc +// +// Generated at Tue May 25 10:39:19 2021 by Giuseppe Cerati using cetskelgen +// from cetlib version v3_11_01. +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "canvas/Persistency/Common/FindManyP.h" + +#include + +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/RecoBase/OpFlash.h" +#include "canvas/Persistency/Common/Assns.h" +#include "art/Persistency/Common/PtrMaker.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/Slice.h" +#include "lardataobj/RecoBase/SpacePoint.h" + +#include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h" + +class IcarusNuSliceHitsProducer; +class IcarusNuSliceHitsProducer : public art::EDProducer { +public: + explicit IcarusNuSliceHitsProducer(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + IcarusNuSliceHitsProducer(IcarusNuSliceHitsProducer const&) = delete; + IcarusNuSliceHitsProducer(IcarusNuSliceHitsProducer&&) = delete; + IcarusNuSliceHitsProducer& operator=(IcarusNuSliceHitsProducer const&) = delete; + IcarusNuSliceHitsProducer& operator=(IcarusNuSliceHitsProducer&&) = delete; + + +private: + + // Declare member data here. + std::string fBaryMatchLabel; + std::string fSliceLabel; + std::string fSpsLabel; + bool doFlashAssn; + + std::vector nearwires; + + // Required functions. + void produce(art::Event& e) override; + +}; + + +IcarusNuSliceHitsProducer::IcarusNuSliceHitsProducer(fhicl::ParameterSet const& p) + : EDProducer{p}, + fBaryMatchLabel( p.get("BaryMatchLabel","tpcpmtbarycentermatchCryoE")), + fSliceLabel(p.get("SliceLabel","pandoraGausCryoE")), + fSpsLabel(p.get("SpsLabel","cluster3DCryoE")), + doFlashAssn(p.get("DoFlashAssn",false)) + // More initializers here. +{ + // Call appropriate produces<>() functions here. + produces>(); + produces>(); + produces>(); + produces>(); + if (doFlashAssn) produces>(); + + // Call appropriate consumes<>() for any products to be retrieved by this module. +} + +void IcarusNuSliceHitsProducer::produce(art::Event& e) +{ + + auto outputHits = std::make_unique >(); + auto assocSliceHitKeys = std::make_unique >(); + + art::PtrMaker hitPtrMaker(e); + art::PtrMaker spsPtrMaker(e); + + auto outputSpacepoints = std::make_unique >(); + auto outputSpsHitAssns = std::make_unique>(); + auto outputSlcFlhAssns = std::make_unique>(); + + art::ValidHandle< std::vector< sbn::TPCPMTBarycenterMatch > > baryMatchListHandle = e.getValidHandle >(fBaryMatchLabel); + art::FindManyP msl(baryMatchListHandle, e, fBaryMatchLabel); + art::FindManyP mfl(baryMatchListHandle, e, fBaryMatchLabel); + float minRad = 99999.; + int slkey = -1; + int flkey = -1; + int bmkey = -1; + for (size_t ibm=0; ibmsize(); ++ibm) { + art::Ptr bm(baryMatchListHandle,ibm); + if (msl.at(ibm).size()>0 && mfl.at(ibm).size()>0) { + std::cout << "ibm=" << ibm << " radius=" << bm->radius << " radius_Trigger=" << bm->radius_Trigger << " flashTime=" << bm->flashTime + << " slkey=" << msl.at(ibm).at(0).key() << " center=" << bm->chargeCenter + << " flkey=" << mfl.at(ibm).at(0).key() << " center=" << bm->flashCenter + << std::endl; + if (bm->radius_Trigger>0 && bm->radius_Triggerradius-bm->radius_Trigger)<0.00001) { + minRad = bm->radius_Trigger; + slkey = msl.at(ibm).at(0).key(); + flkey = mfl.at(ibm).at(0).key(); + bmkey = ibm; + } + } + } + std::cout << "keys: slice=" << slkey << " flash=" << flkey << " bm=" << bmkey << std::endl; + + if (doFlashAssn) { + if (slkey>=0 && bmkey>=0) { + outputSlcFlhAssns->addSingle(msl.at(bmkey).at(0),mfl.at(bmkey).at(0)); + } + } + + art::ValidHandle > inputSlice = e.getValidHandle >(fSliceLabel); + auto assocSliceHit = std::unique_ptr >(new art::FindManyP(inputSlice, e, fSliceLabel)); + + if (slkey>=0) { + std::cout << "slice hits=" << assocSliceHit->at(slkey).size() << std::endl; + for (auto h : assocSliceHit->at(slkey)) { + assocSliceHitKeys->push_back(h.key()); + outputHits->emplace_back(*h); + } + } + + // Get spacepoints from the event record + art::Handle> spListHandle; + std::vector> splist; + if (e.getByLabel(fSpsLabel, spListHandle)) art::fill_ptr_vector(splist, spListHandle); + size_t countsps = 0; + // Get assocations from spacepoints to hits + art::FindManyP fmp(spListHandle, e, fSpsLabel); + for (size_t spIdx = 0; spIdx < splist.size(); ++spIdx) { + auto assochits = fmp.at(spIdx); + if (assochits.size()<3) continue; + //if (splist[spIdx]->Chisq()>0.5) continue; + // + std::vector hitpos; + for (size_t j = 0; j < assochits.size(); ++j) { + auto pos = std::find(assocSliceHitKeys->begin(),assocSliceHitKeys->end(),assochits[j].key()); + if (pos==assocSliceHitKeys->end()) break; + hitpos.push_back(pos-assocSliceHitKeys->begin()); + } + if (hitpos.size()<3) continue; + outputSpacepoints->emplace_back(*splist[spIdx]); + // + const art::Ptr sps = spsPtrMaker(outputSpacepoints->size()-1); + for (size_t j = 0; j < hitpos.size(); ++j) { + const art::Ptr ahp = hitPtrMaker(hitpos[j]); + outputSpsHitAssns->addSingle(ahp,sps); + } + countsps++; + } // for spacepoint + std::cout << "sps count=" << countsps << std::endl; + + e.put(std::move(outputHits)); + e.put(std::move(assocSliceHitKeys)); + e.put(std::move(outputSpacepoints)); + e.put(std::move(outputSpsHitAssns)); + if (doFlashAssn) e.put(std::move(outputSlcFlhAssns)); + +} + +DEFINE_ART_MODULE(IcarusNuSliceHitsProducer) diff --git a/icaruscode/NuGraphIcarus/IcarusTrueNuSliceHitsProducer_module.cc b/icaruscode/NuGraphIcarus/IcarusTrueNuSliceHitsProducer_module.cc new file mode 100644 index 000000000..a81a0de64 --- /dev/null +++ b/icaruscode/NuGraphIcarus/IcarusTrueNuSliceHitsProducer_module.cc @@ -0,0 +1,200 @@ +//////////////////////////////////////////////////////////////////////// +// Class: IcarusTrueNuSliceHitsProducer +// Plugin Type: producer (art v3_06_03) +// File: IcarusTrueNuSliceHitsProducer_module.cc +// +// Generated at Tue May 25 10:39:19 2021 by Giuseppe Cerati using cetskelgen +// from cetlib version v3_11_01. +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "canvas/Persistency/Common/FindManyP.h" + +#include + +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/RecoBase/OpFlash.h" +#include "canvas/Persistency/Common/Assns.h" +#include "art/Persistency/Common/PtrMaker.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/Slice.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/RecoBase/SpacePoint.h" + +#include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h" +#include "larsim/MCCheater/ParticleInventoryService.h" +#include "larsim/MCCheater/BackTrackerService.h" +#include "lardata/DetectorInfoServices/DetectorClocksService.h" + +class IcarusTrueNuSliceHitsProducer; +class IcarusTrueNuSliceHitsProducer : public art::EDProducer { +public: + explicit IcarusTrueNuSliceHitsProducer(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + IcarusTrueNuSliceHitsProducer(IcarusTrueNuSliceHitsProducer const&) = delete; + IcarusTrueNuSliceHitsProducer(IcarusTrueNuSliceHitsProducer&&) = delete; + IcarusTrueNuSliceHitsProducer& operator=(IcarusTrueNuSliceHitsProducer const&) = delete; + IcarusTrueNuSliceHitsProducer& operator=(IcarusTrueNuSliceHitsProducer&&) = delete; + + +private: + + // Declare member data here. + std::string fBaryMatchLabel; + std::string fSliceLabel; + std::string fSpsLabel; + bool doFlashAssn; + + std::vector nearwires; + + // Required functions. + void produce(art::Event& e) override; + +}; + + +IcarusTrueNuSliceHitsProducer::IcarusTrueNuSliceHitsProducer(fhicl::ParameterSet const& p) + : EDProducer{p}, + fBaryMatchLabel( p.get("BaryMatchLabel","tpcpmtbarycentermatchCryoE")), + fSliceLabel(p.get("SliceLabel","pandoraGausCryoE")), + fSpsLabel(p.get("SpsLabel","cluster3DCryoE")), + doFlashAssn(p.get("DoFlashAssn",false)) + // More initializers here. +{ + // Call appropriate produces<>() functions here. + produces>(); + produces>(); + produces>(); + produces>(); + if (doFlashAssn) produces>(); + + // Call appropriate consumes<>() for any products to be retrieved by this module. +} + +void IcarusTrueNuSliceHitsProducer::produce(art::Event& e) +{ + + auto outputHits = std::make_unique >(); + auto assocSliceHitKeys = std::make_unique >(); + + art::PtrMaker hitPtrMaker(e); + art::PtrMaker spsPtrMaker(e); + + auto outputSpacepoints = std::make_unique >(); + auto outputSpsHitAssns = std::make_unique>(); + auto outputSlcFlhAssns = std::make_unique>(); + + art::ServiceHandle pi; + art::ServiceHandle bt_h; + auto const clockData = art::ServiceHandle()->DataFor(e); + + art::ValidHandle > inputSlice = e.getValidHandle >(fSliceLabel); + auto assocSliceHit = std::unique_ptr >(new art::FindManyP(inputSlice, e, fSliceLabel)); + int slkey = -1; + int maxCnt = 0; + for (size_t sk=0; sksize(); sk++) { + int slcnt = 0; + for (auto hit : assocSliceHit->at(sk)) { + std::vector h_ides = bt_h->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick()); + for (auto& tide : h_ides) { + int tid = tide.trackID; + auto mct = pi->TrackIdToMCTruth_P(abs(tid)); + int fromNu = (mct.isAvailable() ? mct->NeutrinoSet() : 0); + if (fromNu==1) { + slcnt++; + break; + } + } + if (slcnt>maxCnt) { + slkey = sk; + maxCnt = slcnt; + } + } + } + std::cout << "keys: slice=" << slkey << std::endl; + + if (slkey>=0) { + std::cout << "slice hits=" << assocSliceHit->at(slkey).size() << std::endl; + for (auto h : assocSliceHit->at(slkey)) { + assocSliceHitKeys->push_back(h.key()); + outputHits->emplace_back(*h); + } + } + + // Get spacepoints from the event record + art::Handle> spListHandle; + std::vector> splist; + if (e.getByLabel(fSpsLabel, spListHandle)) art::fill_ptr_vector(splist, spListHandle); + size_t countsps = 0; + // Get assocations from spacepoints to hits + art::FindManyP fmp(spListHandle, e, fSpsLabel); + for (size_t spIdx = 0; spIdx < splist.size(); ++spIdx) { + // if (splist[spIdx]->Chisq()>0.5) continue; + auto assochits = fmp.at(spIdx); + if (assochits.size()<3) continue; + // + std::vector hitpos; + for (size_t j = 0; j < assochits.size(); ++j) { + auto pos = std::find(assocSliceHitKeys->begin(),assocSliceHitKeys->end(),assochits[j].key()); + if (pos==assocSliceHitKeys->end()) break; + hitpos.push_back(pos-assocSliceHitKeys->begin()); + } + if (hitpos.size()<3) continue; + outputSpacepoints->emplace_back(*splist[spIdx]); + // + const art::Ptr sps = spsPtrMaker(outputSpacepoints->size()-1); + for (size_t j = 0; j < hitpos.size(); ++j) { + const art::Ptr ahp = hitPtrMaker(hitpos[j]); + outputSpsHitAssns->addSingle(ahp,sps); + } + countsps++; + } // for spacepoint + std::cout << "sps count=" << countsps << std::endl; + + art::ValidHandle< std::vector< sbn::TPCPMTBarycenterMatch > > baryMatchListHandle = e.getValidHandle >(fBaryMatchLabel); + art::FindManyP msl(baryMatchListHandle, e, fBaryMatchLabel); + art::FindManyP mfl(baryMatchListHandle, e, fBaryMatchLabel); + float minRad = 99999.; + int mbkey = -1; + for (size_t ibm=0; ibmsize(); ++ibm) { + art::Ptr bm(baryMatchListHandle,ibm); + if (int(msl.at(ibm).at(0).key())!=slkey) continue; + if (mfl.at(ibm).size()>0) { + std::cout << "ibm=" << ibm << " radius=" << bm->radius << " radius_Trigger=" << bm->radius_Trigger << " flashTime=" << bm->flashTime + << " slkey=" << msl.at(ibm).at(0).key() << " center=" << bm->chargeCenter + << " flkey=" << mfl.at(ibm).at(0).key() << " center=" << bm->flashCenter + << std::endl; + if (bm->radiusradius; + } + } + } + std::cout << "mbkey=" << mbkey << std::endl; + if (doFlashAssn){ + if (slkey>=0 && mbkey>=0) { + art::Ptr slp(inputSlice,slkey); + outputSlcFlhAssns->addSingle(slp,mfl.at(mbkey).at(0)); + } + } + + e.put(std::move(outputHits)); + e.put(std::move(assocSliceHitKeys)); + e.put(std::move(outputSpacepoints)); + e.put(std::move(outputSpsHitAssns)); + if (doFlashAssn) e.put(std::move(outputSlcFlhAssns)); +} + +DEFINE_ART_MODULE(IcarusTrueNuSliceHitsProducer) diff --git a/icaruscode/NuGraphIcarus/cafmaker_nugraph_test.fcl b/icaruscode/NuGraphIcarus/cafmaker_nugraph_test.fcl new file mode 100644 index 000000000..e217b591e --- /dev/null +++ b/icaruscode/NuGraphIcarus/cafmaker_nugraph_test.fcl @@ -0,0 +1,6 @@ +#include "cafmakerjob_icarus_detsim2d.fcl" +#### _systtools_and_fluxwgt + +physics.producers.cafmaker.NuGraphSliceHitLabel: "nuslhits" +physics.producers.cafmaker.NuGraphFilterLabel: "NuGraph:filter" +physics.producers.cafmaker.NuGraphSemanticLabel: "NuGraph:semantic" diff --git a/icaruscode/NuGraphIcarus/hdf5maker_icarus_slice.fcl b/icaruscode/NuGraphIcarus/hdf5maker_icarus_slice.fcl new file mode 100644 index 000000000..0a8a137cf --- /dev/null +++ b/icaruscode/NuGraphIcarus/hdf5maker_icarus_slice.fcl @@ -0,0 +1,92 @@ +#include "services_common_icarus.fcl" +#include "services_icarus_simulation.fcl" + +process_name: trueenergy + +services: +{ + TFileService: { fileName: "reco_hist.root" } + TimeTracker: {} + MemoryTracker: {} + RandomNumberGenerator: {} + FileCatalogMetadata: @local::art_file_catalog_mc + @table::icarus_geometry_services + @table::icarus_common_services + @table::icarus_backtracking_services + message: @local::icarus_message_services_prod_debug + ParticleInventoryService: { + ParticleInventory: { + EveIdCalculator: "EmEveIdCalculator" + G4ModuleLabel: "largeant" + OverrideRealData: true + } + } + + +} # services + +source: +{ + module_type: RootInput + maxEvents: -1 +} + +physics: +{ + + analyzers: + { + hdf5maker: { + module_type: "IcarusHDF5Maker" + TruthLabel: "generator" + EventInfo: "nu" + SPLabel: "nuslhits" + HitLabel: "nuslhits" + PFParticleLabel:"pandoraGausCryoE" + OpFlashLabel:"opflashCryoE" + OpHitLabel:"ophit" + UseMap: false + OutputName: "NeutrinoML" + } + } + + producers: + { + #nuslhitsE: { + # module_type: "IcarusNuSliceHitsProducer" + # BaryMatchLabel: "tpcpmtbarycentermatchCryoE" + # SliceLabel: "pandoraGausCryoE" + # SpsLabel: "cluster3DCryoE" + #} + #nuslhitsW: { + # module_type: "IcarusNuSliceHitsProducer" + # BaryMatchLabel: "tpcpmtbarycentermatchCryoW" + # SliceLabel: "pandoraGausCryoW" + # SpsLabel: "cluster3DCryoW" + #} + #nuslhits: { + # module_type: "HitMerger" + # HitProducerLabelVec: ["nuslhitsE","nuslhitsW"] + #} + nuslhits: { + #module_type: "IcarusNuSliceHitsProducer" + module_type: "IcarusTrueNuSliceHitsProducer" + BaryMatchLabel: "tpcpmtbarycentermatchCryoE" + SliceLabel: "pandoraGausCryoE" + SpsLabel: "cluster3DCryoE" + DoFlashAssn: true + } + # sps: cluster3DCryoE + } + + #reco: [nuslhitsE,nuslhitsW,nuslhits] + reco: [nuslhits] + ana: [ hdf5maker ] + + trigger_paths: [ reco ] + end_paths: [ ana ] + +} + +services.BackTrackerService.BackTracker.G4ModuleLabel: "largeant" +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "daq:simpleSC" diff --git a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl new file mode 100644 index 000000000..8a288efe4 --- /dev/null +++ b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl @@ -0,0 +1,115 @@ +#include "nugraph.fcl" +#include "services_common_icarus.fcl" +#include "services_icarus_simulation.fcl" + +process_name: testinference + +services: +{ + TFileService: { fileName: "reco_hist.root" } + TimeTracker: {} + MemoryTracker: {} + RandomNumberGenerator: {} + @table::icarus_geometry_services + @table::icarus_common_services + @table::icarus_backtracking_services + message: @local::icarus_message_services_prod_debug + ParticleInventoryService: { + ParticleInventory: { + EveIdCalculator: "EmEveIdCalculator" + G4ModuleLabel: "largeant" + OverrideRealData: true + } + } +} + +source: +{ + module_type: RootInput + maxEvents: -1 +} + +outputs: +{ + rootOutput: + { + module_type: RootOutput + dataTier: "reconstructed" + compressionLevel: 1 + saveMemoryObjectThreshold: 0 + fileName: "%ifb_%tc-%p.root" + fileProperties: {maxInputFiles: 1} + checkFileName: false + SelectEvents: [] + outputCommands: [ + "keep *_*_*_*" + #"drop *_*_*_*", + #"keep *_*_*_testinference" + ] + } +} + +physics: +{ + + producers: + { + nuslhitsCryoE: { + #module_type: "IcarusNuSliceHitsProducer" + module_type: "IcarusTrueNuSliceHitsProducer" + BaryMatchLabel: "tpcpmtbarycentermatchCryoE" + SliceLabel: "pandoraGausCryoE" + SpsLabel: "cluster3DCryoE" + } + nuslhitsCryoW: { + #module_type: "IcarusNuSliceHitsProducer" + module_type: "IcarusTrueNuSliceHitsProducer" + BaryMatchLabel: "tpcpmtbarycentermatchCryoW" + SliceLabel: "pandoraGausCryoW" + SpsLabel: "cluster3DCryoW" + } + nuslhits: { + module_type: "HitMerger" + HitProducerLabelVec: ["nuslhitsCryoE","nuslhitsCryoW"] + } + NuGraphCryoE: @local::ApptainerNuGraphNuSonicTriton + NuGraphCryoW: @local::ApptainerNuGraphNuSonicTriton + } + analyzers: + { + NuGraphAnaCryoE: { + module_type: "NuGraphAnalyzer" + NuGraphLabel: "NuGraphCryoE" + } + NuGraphAnaCryoW: { + module_type: "NuGraphAnalyzer" + NuGraphLabel: "NuGraphCryoW" + } + } + + reco: [nuslhitsCryoE, nuslhitsCryoW, NuGraphCryoE, NuGraphCryoW] + #reco: [ nuslhitsCryoE, NuGraphCryoE] + trigger_paths: [ reco ] + + ana: [NuGraphAnaCryoE,NuGraphAnaCryoW] + #ana: [NuGraphAnaCryoE] + streamROOT: [ rootOutput ] + end_paths: [ ana, streamROOT ] +} + +physics.producers.NuGraphCryoE.LoaderTool.tool_type: "IcarusNuGraphLoader" +physics.producers.NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE" +physics.producers.NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE" +physics.producers.NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE" +physics.producers.NuGraphCryoE.TritonConfig.modelName: "nugraph2_icarus_mpvmprbnb" +#physics.producers.NuGraphCryoE.debug: true + +physics.producers.NuGraphCryoW.LoaderTool.tool_type: "IcarusNuGraphLoader" +physics.producers.NuGraphCryoW.LoaderTool.hitInput: "nuslhitsCryoW" +physics.producers.NuGraphCryoW.LoaderTool.spsInput: "nuslhitsCryoW" +physics.producers.NuGraphCryoW.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoW" +physics.producers.NuGraphCryoW.TritonConfig.modelName: "nugraph2_icarus_mpvmprbnb" +#physics.producers.NuGraphCryoW.debug: true + +services.BackTrackerService.BackTracker.G4ModuleLabel: "largeant" +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "daq:simpleSC" From d728cdbeab5b2ae936e48782c9dc5e3e623eecc4 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Thu, 10 Apr 2025 13:38:43 -0500 Subject: [PATCH 02/24] use IcarusNuSliceHitsProducer --- icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl index 8a288efe4..a41d7cfc3 100644 --- a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl +++ b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl @@ -55,15 +55,15 @@ physics: producers: { nuslhitsCryoE: { - #module_type: "IcarusNuSliceHitsProducer" - module_type: "IcarusTrueNuSliceHitsProducer" + module_type: "IcarusNuSliceHitsProducer" + #module_type: "IcarusTrueNuSliceHitsProducer" BaryMatchLabel: "tpcpmtbarycentermatchCryoE" SliceLabel: "pandoraGausCryoE" SpsLabel: "cluster3DCryoE" } nuslhitsCryoW: { - #module_type: "IcarusNuSliceHitsProducer" - module_type: "IcarusTrueNuSliceHitsProducer" + module_type: "IcarusNuSliceHitsProducer" + #module_type: "IcarusTrueNuSliceHitsProducer" BaryMatchLabel: "tpcpmtbarycentermatchCryoW" SliceLabel: "pandoraGausCryoW" SpsLabel: "cluster3DCryoW" From 6f211ceeb85bd9b33e65adfa49d696fa8bb7a3ee Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Thu, 10 Apr 2025 17:15:05 -0500 Subject: [PATCH 03/24] add script to start local triton server --- icaruscode/NuGraphIcarus/CMakeLists.txt | 2 ++ icaruscode/NuGraphIcarus/scripts/CMakeLists.txt | 3 +++ .../scripts/setup_tritonservier-nugraph-v0.sh | 7 +++++++ icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl | 6 ++++++ 4 files changed, 18 insertions(+) create mode 100644 icaruscode/NuGraphIcarus/scripts/CMakeLists.txt create mode 100755 icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh diff --git a/icaruscode/NuGraphIcarus/CMakeLists.txt b/icaruscode/NuGraphIcarus/CMakeLists.txt index 4c269b433..fdfabcdc8 100644 --- a/icaruscode/NuGraphIcarus/CMakeLists.txt +++ b/icaruscode/NuGraphIcarus/CMakeLists.txt @@ -57,6 +57,8 @@ simple_plugin(IcarusHDF5Maker "module" ${HDF5_LIBRARIES} ) +add_subdirectory(scripts) + install_headers() install_fhicl() install_source() diff --git a/icaruscode/NuGraphIcarus/scripts/CMakeLists.txt b/icaruscode/NuGraphIcarus/scripts/CMakeLists.txt new file mode 100644 index 000000000..7b81c7840 --- /dev/null +++ b/icaruscode/NuGraphIcarus/scripts/CMakeLists.txt @@ -0,0 +1,3 @@ +install_scripts( + "setup_tritonservier-nugraph-v0.sh" + ) diff --git a/icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh b/icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh new file mode 100755 index 000000000..31df7e5fc --- /dev/null +++ b/icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh @@ -0,0 +1,7 @@ +#!/bin/bash + +unset PYTHONHOME +unset PYTHONPATH +rm -f tritonserver_nugraph-v0.log +export APPTAINER_BIND=/etc/hosts,/tmp,/cvmfs +/cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/bin/apptainer run --pid --ipc --home ~/:${HOME} --pwd ${PWD} /cvmfs/icarus.opensciencegrid.org/containers/tritonserver/nugraph-v0/ >& tritonserver_nugraph-v0.log & diff --git a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl index a41d7cfc3..27fd2f952 100644 --- a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl +++ b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl @@ -2,6 +2,12 @@ #include "services_common_icarus.fcl" #include "services_icarus_simulation.fcl" +# in order to setup the triton server locally, please do: +# > mrbslp +# > setup_nugraph-v0.sh +# then after the job you can kill the server process with: +# > killall -9 /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/x86_64/libexec/apptainer/libexec/starter + process_name: testinference services: From dc9e7065f56d564e0e3d8e098a48c770e5673f12 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Thu, 10 Apr 2025 22:47:53 -0500 Subject: [PATCH 04/24] fix script name --- icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl index 27fd2f952..95eaa3789 100644 --- a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl +++ b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl @@ -4,7 +4,7 @@ # in order to setup the triton server locally, please do: # > mrbslp -# > setup_nugraph-v0.sh +# > setup_tritonservier-nugraph-v0.sh # then after the job you can kill the server process with: # > killall -9 /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/x86_64/libexec/apptainer/libexec/starter From 81247b49d4a35621493df1f557fdf7710d85e6c7 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Fri, 11 Apr 2025 13:39:00 -0500 Subject: [PATCH 05/24] take a quick nap to avoid starting a lar job before the server is ready --- .../NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh b/icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh index 31df7e5fc..a971b9654 100755 --- a/icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh +++ b/icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh @@ -5,3 +5,4 @@ unset PYTHONPATH rm -f tritonserver_nugraph-v0.log export APPTAINER_BIND=/etc/hosts,/tmp,/cvmfs /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/bin/apptainer run --pid --ipc --home ~/:${HOME} --pwd ${PWD} /cvmfs/icarus.opensciencegrid.org/containers/tritonserver/nugraph-v0/ >& tritonserver_nugraph-v0.log & +sleep 90 From 23d0dfd14f66d29cc7af45191a4c37c045fd04aa Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Thu, 17 Apr 2025 17:09:48 -0500 Subject: [PATCH 06/24] fix script names, add version for grid processing --- .../NuGraphIcarus/scripts/CMakeLists.txt | 3 +- ...v0.sh => setup_tritonserver-nugraph-v0.sh} | 0 .../setup_tritonserver-nugraph-v0_grid.sh | 31 +++++++++++++++++++ 3 files changed, 33 insertions(+), 1 deletion(-) rename icaruscode/NuGraphIcarus/scripts/{setup_tritonservier-nugraph-v0.sh => setup_tritonserver-nugraph-v0.sh} (100%) create mode 100755 icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh diff --git a/icaruscode/NuGraphIcarus/scripts/CMakeLists.txt b/icaruscode/NuGraphIcarus/scripts/CMakeLists.txt index 7b81c7840..e0f7a8d2b 100644 --- a/icaruscode/NuGraphIcarus/scripts/CMakeLists.txt +++ b/icaruscode/NuGraphIcarus/scripts/CMakeLists.txt @@ -1,3 +1,4 @@ install_scripts( - "setup_tritonservier-nugraph-v0.sh" + "setup_tritonserver-nugraph-v0.sh" + "setup_tritonserver-nugraph-v0_grid.sh" ) diff --git a/icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh b/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh similarity index 100% rename from icaruscode/NuGraphIcarus/scripts/setup_tritonservier-nugraph-v0.sh rename to icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh diff --git a/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh b/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh new file mode 100755 index 000000000..af6f1392a --- /dev/null +++ b/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh @@ -0,0 +1,31 @@ +#!/bin/bash +# +unset PYTHONHOME +unset PYTHONPATH +rm -f tritonserver_nugraph-v0.log +export APPTAINER_BIND=/etc/hosts,/tmp,/cvmfs +# +export BASEPORT=8000 +export NPORTS=3 +# +export HTTPPORT=$BASEPORT +export GRPCPORT=$((BASEPORT+1)) +export METRPORT=$((BASEPORT+2)) +# +if [ -z "$FCL" ]; then + export FCL=$1 +fi +# +while 2>/dev/null >"/dev/tcp/0.0.0.0/$BASEPORT"; do + BASEPORT=$((BASEPORT+NPORTS)) + export HTTPPORT=$BASEPORT + export GRPCPORT=$((BASEPORT+1)) + export METRPORT=$((BASEPORT+2)) +done +# +echo "physics.producers.NuGraphCryoE.TritonConfig.serverURL: \"localhost:"$GRPCPORT"\"" >> $FCL +echo "physics.producers.NuGraphCryoW.TritonConfig.serverURL: \"localhost:"$GRPCPORT"\"" >> $FCL +# +/cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/bin/apptainer exec --pid --ipc --home ~/:${HOME} --pwd ${PWD} /cvmfs/icarus.opensciencegrid.org/containers/tritonserver/nugraph-v0/ tritonserver --model-repository /triton-server-config/models --http-port=${HTTPPORT} --grpc-port=${GRPCPORT} --metrics-port=${METRPORT} >& tritonserver_nugraph-v0.log & +# +sleep 90 From 17f1e7b78b6f0b3c8f753116e3516eda0a05b010 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Thu, 17 Apr 2025 17:31:11 -0500 Subject: [PATCH 07/24] no sleep better wait --- .../NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh | 2 +- .../NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh b/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh index a971b9654..31506929f 100755 --- a/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh +++ b/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh @@ -5,4 +5,4 @@ unset PYTHONPATH rm -f tritonserver_nugraph-v0.log export APPTAINER_BIND=/etc/hosts,/tmp,/cvmfs /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/bin/apptainer run --pid --ipc --home ~/:${HOME} --pwd ${PWD} /cvmfs/icarus.opensciencegrid.org/containers/tritonserver/nugraph-v0/ >& tritonserver_nugraph-v0.log & -sleep 90 +grep -q "Started" <(tail -f tritonserver_nugraph-v0.log) diff --git a/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh b/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh index af6f1392a..d6ed36a09 100755 --- a/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh +++ b/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh @@ -28,4 +28,4 @@ echo "physics.producers.NuGraphCryoW.TritonConfig.serverURL: \"localhost:"$GRPCP # /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/bin/apptainer exec --pid --ipc --home ~/:${HOME} --pwd ${PWD} /cvmfs/icarus.opensciencegrid.org/containers/tritonserver/nugraph-v0/ tritonserver --model-repository /triton-server-config/models --http-port=${HTTPPORT} --grpc-port=${GRPCPORT} --metrics-port=${METRPORT} >& tritonserver_nugraph-v0.log & # -sleep 90 +grep -q "Started" <(tail -f tritonserver_nugraph-v0.log) From d9bed1cdf91f1eb54234b4f02404bd164799896f Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Fri, 2 May 2025 12:59:37 -0500 Subject: [PATCH 08/24] add NuGraph to reco fcl, switch to libtorch inference by default --- fcl/caf/cafmaker_defs.fcl | 3 ++ fcl/reco/Definitions/stage1_icarus_defs.fcl | 15 +++++++ fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl | 1 + .../Stage1/Run2/stage1_run2_icarus_MC.fcl | 1 + .../NuGraphIcarus/IcarusNuGraphLoader_tool.cc | 3 +- icaruscode/NuGraphIcarus/nugraph_icarus.fcl | 45 +++++++++++++++++++ .../testinference_slice_icarus.fcl | 44 +++++------------- 7 files changed, 79 insertions(+), 33 deletions(-) create mode 100644 icaruscode/NuGraphIcarus/nugraph_icarus.fcl diff --git a/fcl/caf/cafmaker_defs.fcl b/fcl/caf/cafmaker_defs.fcl index 6c676dc26..6b58e7c6f 100644 --- a/fcl/caf/cafmaker_defs.fcl +++ b/fcl/caf/cafmaker_defs.fcl @@ -334,6 +334,9 @@ cafmaker.SystWeightLabels: ["genieweight", "fluxweight"] cafmaker.SaveGENIEEventRecord: true # save GENIE event record by default. Turn this off for data cafmaker fcl cafmaker.TPCPMTBarycenterMatchLabel: "tpcpmtbarycentermatch" cafmaker.TrackHitFillRREndCut: 30 # include entire PID region +cafmaker.NuGraphSliceHitLabel: "nuslhits" +cafmaker.NuGraphFilterLabel: "NuGraph:filter" +cafmaker.NuGraphSemanticLabel: "NuGraph:semantic" # Add CAFMaker to the list of producers caf_preprocess_producers.cafmaker: @local::cafmaker diff --git a/fcl/reco/Definitions/stage1_icarus_defs.fcl b/fcl/reco/Definitions/stage1_icarus_defs.fcl index dd82375cc..d3ae779fe 100644 --- a/fcl/reco/Definitions/stage1_icarus_defs.fcl +++ b/fcl/reco/Definitions/stage1_icarus_defs.fcl @@ -22,6 +22,7 @@ #include "supera_modules.fcl" #include "crtpmtmatching_parameters.fcl" #include "tpcpmtbarycentermatch_config.fcl" +#include "nugraph_icarus.fcl" BEGIN_PROLOG @@ -71,6 +72,13 @@ icarus_stage1_producers: tpcpmtbarycentermatchCryoE: @local::data_tpcpmtbarycentermatchproducer_east tpcpmtbarycentermatchCryoW: @local::data_tpcpmtbarycentermatchproducer_west + + ## NuGraph + nuslhitsCryoE: @local::nuslhitsCryoE + nuslhitsCryoW: @local::nuslhitsCryoW + NuGraphCryoE: @local::NuGraphCryoE + NuGraphCryoW: @local::NuGraphCryoW + } icarus_stage1_filters: @@ -217,6 +225,13 @@ icarus_tpcpmtbarycentermatch: [ tpcpmtbarycentermatchCryoW ] +icarus_nugraph: [ + nuslhitsCryoE, + nuslhitsCryoW, + NuGraphCryoE, + NuGraphCryoW + ] + icarus_crttrack: [crttrack] icarus_crtt0match: [CRTT0Matching] diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl index 40ea40c2a..3949690af 100644 --- a/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl +++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl @@ -6,6 +6,7 @@ physics.reco: [ @sequence::icarus_filter_cluster3D, @sequence::icarus_pandora_Gauss, @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, + @sequence::icarus_nugraph, @sequence::icarus_crttrack, @sequence::icarus_crtt0match, caloskimCalorimetryCryoE, caloskimCalorimetryCryoW] diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl index 8401564f1..86f7d09a7 100644 --- a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl +++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl @@ -52,6 +52,7 @@ physics.reco: [ @sequence::icarus_reco_Gauss_CryoW , @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, + @sequence::icarus_nugraph, caloskimCalorimetryCryoE, caloskimCalorimetryCryoW, mcassociationsGausCryoE, mcassociationsGausCryoW ] diff --git a/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc b/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc index 8062ba80d..dedb66b8d 100644 --- a/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc +++ b/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc @@ -129,9 +129,10 @@ void IcarusNuGraphLoader::loadData(art::Event& e, } // space point table + size_t spidx = 0; for (size_t i = 0; i < splist.size(); ++i) { if (sp2Hit[i].size()==0) continue; - spacepoint_table_spacepoint_id_data.push_back(i); + spacepoint_table_spacepoint_id_data.push_back(spidx++); spacepoint_table_hit_id_u_data.push_back(-1); spacepoint_table_hit_id_v_data.push_back(-1); spacepoint_table_hit_id_y_data.push_back(-1); diff --git a/icaruscode/NuGraphIcarus/nugraph_icarus.fcl b/icaruscode/NuGraphIcarus/nugraph_icarus.fcl new file mode 100644 index 000000000..24d2b39c1 --- /dev/null +++ b/icaruscode/NuGraphIcarus/nugraph_icarus.fcl @@ -0,0 +1,45 @@ +#include "nugraph.fcl" + +BEGIN_PROLOG + +nuslhitsCryoE: { + module_type: "IcarusNuSliceHitsProducer" + BaryMatchLabel: "tpcpmtbarycentermatchCryoE" + SliceLabel: "pandoraGausCryoE" + SpsLabel: "cluster3DCryoE" +} + +nuslhitsCryoW: { + module_type: "IcarusNuSliceHitsProducer" + BaryMatchLabel: "tpcpmtbarycentermatchCryoW" + SliceLabel: "pandoraGausCryoW" + SpsLabel: "cluster3DCryoW" +} + +#NuGraphCryoE: @local::ApptainerNuGraphNuSonicTriton +#NuGraphCryoE.LoaderTool.tool_type: "IcarusNuGraphLoader" +#NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE" +#NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE" +#NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE" +#NuGraphCryoE.TritonConfig.modelName: "nugraph2_icarus_mpvmprbnb" + +NuGraphCryoE: @local::NuGraphLibTorch +NuGraphCryoE.LoaderTool.tool_type: "IcarusNuGraphLoader" +NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE" +NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE" +NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE" +NuGraphCryoE.avgs_u: [168.04594 , 178.3245 , 266.6149 , 3.857218 ] +NuGraphCryoE.avgs_v: [1245.3547 , 176.54117 , 323.52786, 4.3267984] +NuGraphCryoE.avgs_y: [1225.5012 , 183.58075 , 310.83493, 4.3409133] +NuGraphCryoE.devs_u: [ 82.80644 , 67.60649 , 274.32666, 1.2912455] +NuGraphCryoE.devs_v: [ 293.06314, 66.8194 , 322.11386, 1.4249923] +NuGraphCryoE.devs_y: [ 307.1943 , 67.063324, 312.461 , 1.4532351] +NuGraphCryoE.modelFileName: "model_mpvmpr_bnb_numu_cos.pt" +#NuGraphCryoE.debug: true + +NuGraphCryoW: @local::NuGraphCryoE +NuGraphCryoW.LoaderTool.hitInput: "nuslhitsCryoW" +NuGraphCryoW.LoaderTool.spsInput: "nuslhitsCryoW" +NuGraphCryoW.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoW" + +END_PROLOG diff --git a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl index 95eaa3789..66902fcd0 100644 --- a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl +++ b/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl @@ -1,4 +1,4 @@ -#include "nugraph.fcl" +#include "nugraph_icarus.fcl" #include "services_common_icarus.fcl" #include "services_icarus_simulation.fcl" @@ -60,26 +60,14 @@ physics: producers: { - nuslhitsCryoE: { - module_type: "IcarusNuSliceHitsProducer" - #module_type: "IcarusTrueNuSliceHitsProducer" - BaryMatchLabel: "tpcpmtbarycentermatchCryoE" - SliceLabel: "pandoraGausCryoE" - SpsLabel: "cluster3DCryoE" - } - nuslhitsCryoW: { - module_type: "IcarusNuSliceHitsProducer" - #module_type: "IcarusTrueNuSliceHitsProducer" - BaryMatchLabel: "tpcpmtbarycentermatchCryoW" - SliceLabel: "pandoraGausCryoW" - SpsLabel: "cluster3DCryoW" - } - nuslhits: { - module_type: "HitMerger" - HitProducerLabelVec: ["nuslhitsCryoE","nuslhitsCryoW"] - } - NuGraphCryoE: @local::ApptainerNuGraphNuSonicTriton - NuGraphCryoW: @local::ApptainerNuGraphNuSonicTriton + nuslhitsCryoE: @local::nuslhitsCryoE + nuslhitsCryoW: @local::nuslhitsCryoW + #nuslhits: { + # module_type: "HitMerger" + # HitProducerLabelVec: ["nuslhitsCryoE","nuslhitsCryoW"] + #} + NuGraphCryoE: @local::NuGraphCryoE + NuGraphCryoW: @local::NuGraphCryoW } analyzers: { @@ -103,18 +91,10 @@ physics: end_paths: [ ana, streamROOT ] } -physics.producers.NuGraphCryoE.LoaderTool.tool_type: "IcarusNuGraphLoader" -physics.producers.NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE" -physics.producers.NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE" -physics.producers.NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE" -physics.producers.NuGraphCryoE.TritonConfig.modelName: "nugraph2_icarus_mpvmprbnb" -#physics.producers.NuGraphCryoE.debug: true +#physics.producers.nuslhitsCryoE.module_type: "IcarusTrueNuSliceHitsProducer" +#physics.producers.nuslhitsCryoW.module_type: "IcarusTrueNuSliceHitsProducer" -physics.producers.NuGraphCryoW.LoaderTool.tool_type: "IcarusNuGraphLoader" -physics.producers.NuGraphCryoW.LoaderTool.hitInput: "nuslhitsCryoW" -physics.producers.NuGraphCryoW.LoaderTool.spsInput: "nuslhitsCryoW" -physics.producers.NuGraphCryoW.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoW" -physics.producers.NuGraphCryoW.TritonConfig.modelName: "nugraph2_icarus_mpvmprbnb" +#physics.producers.NuGraphCryoE.debug: true #physics.producers.NuGraphCryoW.debug: true services.BackTrackerService.BackTracker.G4ModuleLabel: "largeant" From 8deda1bbdfb9d88d347f01f1cf21db4922c7ffd3 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Fri, 2 May 2025 14:25:29 -0500 Subject: [PATCH 09/24] make nugraph part of the pandora sequence --- fcl/reco/Definitions/stage1_icarus_defs.fcl | 11 ++++------- fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl | 1 - fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl | 1 - 3 files changed, 4 insertions(+), 9 deletions(-) diff --git a/fcl/reco/Definitions/stage1_icarus_defs.fcl b/fcl/reco/Definitions/stage1_icarus_defs.fcl index d3ae779fe..bda3ddd64 100644 --- a/fcl/reco/Definitions/stage1_icarus_defs.fcl +++ b/fcl/reco/Definitions/stage1_icarus_defs.fcl @@ -188,12 +188,16 @@ icarus_reco_cluster3DCryoW: [ cluster3DCryoW ] icarus_reco_cluster3DCryoE: [ cluster3DCryoE ] icarus_reco_pandoraGausCryoW: [ pandoraGausCryoW, + nuslhitsCryoW, + NuGraphCryoW, pandoraTrackGausCryoW, pandoraKalmanTrackGausCryoW, SBNShowerGausCryoW ] icarus_reco_pandoraGausCryoE: [ pandoraGausCryoE, + nuslhitsCryoE, + NuGraphCryoE, pandoraTrackGausCryoE, pandoraKalmanTrackGausCryoE, SBNShowerGausCryoE @@ -225,13 +229,6 @@ icarus_tpcpmtbarycentermatch: [ tpcpmtbarycentermatchCryoW ] -icarus_nugraph: [ - nuslhitsCryoE, - nuslhitsCryoW, - NuGraphCryoE, - NuGraphCryoW - ] - icarus_crttrack: [crttrack] icarus_crtt0match: [CRTT0Matching] diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl index 3949690af..40ea40c2a 100644 --- a/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl +++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl @@ -6,7 +6,6 @@ physics.reco: [ @sequence::icarus_filter_cluster3D, @sequence::icarus_pandora_Gauss, @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, - @sequence::icarus_nugraph, @sequence::icarus_crttrack, @sequence::icarus_crtt0match, caloskimCalorimetryCryoE, caloskimCalorimetryCryoW] diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl index 86f7d09a7..8401564f1 100644 --- a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl +++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl @@ -52,7 +52,6 @@ physics.reco: [ @sequence::icarus_reco_Gauss_CryoW , @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, - @sequence::icarus_nugraph, caloskimCalorimetryCryoE, caloskimCalorimetryCryoW, mcassociationsGausCryoE, mcassociationsGausCryoW ] From e87380b5bff2d7a509f5ef2bb421c8cbe1cfb83d Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Mon, 5 May 2025 17:40:16 -0500 Subject: [PATCH 10/24] add icarus_postPandora_Gauss --- fcl/reco/Definitions/stage1_icarus_defs.fcl | 17 ++++++++++++----- fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl | 1 + fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl | 4 ++-- 3 files changed, 15 insertions(+), 7 deletions(-) diff --git a/fcl/reco/Definitions/stage1_icarus_defs.fcl b/fcl/reco/Definitions/stage1_icarus_defs.fcl index bda3ddd64..446c6072f 100644 --- a/fcl/reco/Definitions/stage1_icarus_defs.fcl +++ b/fcl/reco/Definitions/stage1_icarus_defs.fcl @@ -187,16 +187,18 @@ icarus_reco_cluster3DCryoW: [ cluster3DCryoW ] icarus_reco_cluster3DCryoE: [ cluster3DCryoE ] -icarus_reco_pandoraGausCryoW: [ pandoraGausCryoW, - nuslhitsCryoW, +icarus_reco_pandoraGausCryoW: [ pandoraGausCryoW ] + +icarus_reco_pandoraGausCryoE: [ pandoraGausCryoE ] + +icarus_reco_postPandoraGausCryoW: [ nuslhitsCryoW, NuGraphCryoW, pandoraTrackGausCryoW, pandoraKalmanTrackGausCryoW, SBNShowerGausCryoW ] -icarus_reco_pandoraGausCryoE: [ pandoraGausCryoE, - nuslhitsCryoE, +icarus_reco_postPandoraGausCryoE: [ nuslhitsCryoE, NuGraphCryoE, pandoraTrackGausCryoE, pandoraKalmanTrackGausCryoE, @@ -218,6 +220,11 @@ icarus_pandora_Gauss: [ @sequence::icarus_reco_pandoraGausCryoW ] +icarus_postPandora_Gauss: [ + @sequence::icarus_reco_postPandoraGausCryoE, + @sequence::icarus_reco_postPandoraGausCryoW + ] + #Add flash matching icarus_reco_fm: [ fmatchCryoE, fmatchCryoW, @@ -236,7 +243,7 @@ icarus_crtt0match: [CRTT0Matching] icarus_crtt0match_eff: [CRTT0MatchingE, CRTT0MatchingW] ### Below we include overrides for the modules above - +q ## Overrides for filtering of cluster3D hits icarus_stage1_filters.TPCHitFilterCryoW.HitDataLabelVec: ["cluster3DCryoW"] icarus_stage1_filters.TPCHitFilterCryoW.MaximumHits: 60000 diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl index 40ea40c2a..fe416f761 100644 --- a/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl +++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus.fcl @@ -6,6 +6,7 @@ physics.reco: [ @sequence::icarus_filter_cluster3D, @sequence::icarus_pandora_Gauss, @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, + @sequence::icarus_postPandora_Gauss, @sequence::icarus_crttrack, @sequence::icarus_crtt0match, caloskimCalorimetryCryoE, caloskimCalorimetryCryoW] diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl index 8401564f1..4d111f063 100644 --- a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl +++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl @@ -48,10 +48,10 @@ physics.producers: { } physics.reco: [ - @sequence::icarus_reco_Gauss_CryoE , - @sequence::icarus_reco_Gauss_CryoW , + @sequence::icarus_pandora_Gauss , @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, + @sequence::icarus_postPandora_Gauss , caloskimCalorimetryCryoE, caloskimCalorimetryCryoW, mcassociationsGausCryoE, mcassociationsGausCryoW ] From 0aa9cdcc03c7b78f70c68fd89357e0ddc7c46bb3 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Mon, 5 May 2025 18:33:12 -0500 Subject: [PATCH 11/24] remove spurious character --- fcl/reco/Definitions/stage1_icarus_defs.fcl | 1 - 1 file changed, 1 deletion(-) diff --git a/fcl/reco/Definitions/stage1_icarus_defs.fcl b/fcl/reco/Definitions/stage1_icarus_defs.fcl index 446c6072f..3a6d9d394 100644 --- a/fcl/reco/Definitions/stage1_icarus_defs.fcl +++ b/fcl/reco/Definitions/stage1_icarus_defs.fcl @@ -243,7 +243,6 @@ icarus_crtt0match: [CRTT0Matching] icarus_crtt0match_eff: [CRTT0MatchingE, CRTT0MatchingW] ### Below we include overrides for the modules above -q ## Overrides for filtering of cluster3D hits icarus_stage1_filters.TPCHitFilterCryoW.HitDataLabelVec: ["cluster3DCryoW"] icarus_stage1_filters.TPCHitFilterCryoW.MaximumHits: 60000 From cb5276c1c0e3d35648e0340f338e58759cba9f2c Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Thu, 8 May 2025 12:09:17 -0500 Subject: [PATCH 12/24] restore cluster3D --- fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl index 4d111f063..fb445b64e 100644 --- a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl +++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl @@ -48,10 +48,11 @@ physics.producers: { } physics.reco: [ - @sequence::icarus_pandora_Gauss , + @sequence::icarus_filter_cluster3D, + @sequence::icarus_pandora_Gauss, @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, - @sequence::icarus_postPandora_Gauss , + @sequence::icarus_postPandora_Gauss, caloskimCalorimetryCryoE, caloskimCalorimetryCryoW, mcassociationsGausCryoE, mcassociationsGausCryoW ] From 70b8af6bbf932a60db67405e077434e765936c89 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Thu, 29 May 2025 16:33:43 -0500 Subject: [PATCH 13/24] add icarus_postPandora_Gauss to 1d fcl as well --- fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl | 1 + 1 file changed, 1 insertion(+) diff --git a/fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl b/fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl index 85b155397..cbf9f4e40 100644 --- a/fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl +++ b/fcl/reco/Stage1/mc/stage1_run2_1d_icarus_MC.fcl @@ -6,6 +6,7 @@ physics.reco: [ @sequence::icarus_reco_Gauss1D_CryoW , @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, + @sequence::icarus_postPandora_Gauss, @sequence::icarus_crttrack, @sequence::icarus_crtt0tagging, caloskimCalorimetryCryoE, caloskimCalorimetryCryoW, From 6f1b8b7b6beb368866bb9f273cfa3f8a439ea92e Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Thu, 5 Jun 2025 19:17:41 -0500 Subject: [PATCH 14/24] address G.P.'s comments --- icaruscode/NuGraphIcarus/CMakeLists.txt | 15 +- .../NuGraphIcarus/IcarusHDF5Maker_module.cc | 834 ++++++++---------- .../NuGraphIcarus/IcarusNuGraphLoader_tool.cc | 104 +-- .../IcarusNuSliceHitsProducer_module.cc | 120 +-- .../IcarusTrueNuSliceHitsProducer_module.cc | 116 +-- icaruscode/NuGraphIcarus/StitchingUtils.h | 50 ++ .../NuGraphIcarus/hdf5maker_icarus_slice.fcl | 2 - .../testinference_slice_icarus.fcl | 15 +- 8 files changed, 613 insertions(+), 643 deletions(-) create mode 100644 icaruscode/NuGraphIcarus/StitchingUtils.h diff --git a/icaruscode/NuGraphIcarus/CMakeLists.txt b/icaruscode/NuGraphIcarus/CMakeLists.txt index fdfabcdc8..0392247c3 100644 --- a/icaruscode/NuGraphIcarus/CMakeLists.txt +++ b/icaruscode/NuGraphIcarus/CMakeLists.txt @@ -2,7 +2,6 @@ cet_enable_asserts() set( nugraph_tool_lib_list lardataobj::RecoBase - lardataobj::AnalysisBase larrecodnn::NuGraphBaseTools TorchScatter::TorchScatter art::Framework_Core @@ -20,11 +19,7 @@ cet_build_plugin(IcarusNuGraphLoader art::tool simple_plugin(IcarusNuSliceHitsProducer "module" LIBRARIES PRIVATE art::Framework_Core - nusimdata::SimulationBase lardataobj::RecoBase - lardataalg::DetectorInfo - lardata::DetectorInfoServices_DetectorClocksServiceStandard_service - larcore::Geometry_Geometry_service sbnobj::Common_Reco ) @@ -34,11 +29,10 @@ simple_plugin(IcarusTrueNuSliceHitsProducer "module" nusimdata::SimulationBase lardataobj::RecoBase lardataalg::DetectorInfo - lardata::DetectorInfoServices_DetectorClocksServiceStandard_service - larcore::Geometry_Geometry_service sbnobj::Common_Reco larsim::MCCheater_BackTrackerService_service larsim::MCCheater_ParticleInventoryService_service + lardata::DetectorClocksService ) simple_plugin(IcarusHDF5Maker "module" @@ -48,12 +42,11 @@ simple_plugin(IcarusHDF5Maker "module" larsim::MCCheater_BackTrackerService_service larsim::MCCheater_ParticleInventoryService_service hep_hpc_hdf5 - #hep::hpc_Utilities - lardata::RecoBaseProxy lardataobj::RecoBase lardataalg::DetectorInfo - lardata::DetectorInfoServices_DetectorClocksServiceStandard_service - larcore::Geometry_Geometry_service + lardata::DetectorClocksService + larcorealg::Geometry + larcore::WireReadout ${HDF5_LIBRARIES} ) diff --git a/icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc b/icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc index a2eeb0515..a6ff7e9ee 100644 --- a/icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc +++ b/icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc @@ -1,17 +1,13 @@ -//////////////////////////////////////////////////////////////////////// -//// Class: IcarusHDF5Maker -//// Plugin Type: analyzer (art v3_06_03) -//// File: IcarusHDF5Maker_module.cc -//// -//// Generated at Wed May 5 08:23:31 2021 by V Hewes using cetskelgen -/// from cetlib version v3_11_01. -////////////////////////////////////////////////////////////////////////// +/** + * @file icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc + * @brief Implementation of `IcarusHDF5Maker` _art_ module. + * @author Giuseppe Cerati (cerati@fnal.gov), V Hewes + */ #include "art/Framework/Core/EDAnalyzer.h" #include "art/Framework/Core/ModuleMacros.h" #include "art/Framework/Principal/Event.h" #include "art/Framework/Principal/Handle.h" -#include "art/Framework/Principal/Run.h" #include "art/Framework/Principal/SubRun.h" #include "canvas/Utilities/InputTag.h" #include "fhiclcpp/ParameterSet.h" @@ -20,25 +16,17 @@ #include "larcore/Geometry/Geometry.h" #include "larcorealg/Geometry/GeometryCore.h" -#include "larcorealg/Geometry/TPCGeo.h" +#include "larcorealg/Geometry/PlaneGeo.h" #include "larcore/Geometry/WireReadout.h" +#include "larcorealg/Geometry/WireReadoutGeom.h" #include "nusimdata/SimulationBase/MCTruth.h" #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/SpacePoint.h" -#include "lardataobj/RecoBase/Wire.h" #include "lardataobj/RecoBase/OpHit.h" #include "lardataobj/RecoBase/OpFlash.h" -#include "lardata/RecoBaseProxy/ProxyBase.h" -#include "lardataobj/RecoBase/Wire.h" #include "larcorealg/Geometry/Exceptions.h" #include "lardataobj/RecoBase/Slice.h" -#include "lardataobj/RecoBase/PFParticle.h" -#include "lardata/RecoBaseProxy/ProxyBase.h" -#include "lardataobj/RecoBase/Vertex.h" -#include "lardataobj/RecoBase/Cluster.h" -#include "lardataobj/AnalysisBase/T0.h" -#include "lardataobj/RecoBase/PFParticleMetadata.h" #include "larsim/MCCheater/BackTrackerService.h" #include "larsim/MCCheater/ParticleInventoryService.h" @@ -49,6 +37,12 @@ #include "hep_hpc/hdf5/make_ntuple.hpp" #include "art/Framework/Services/Registry/ServiceHandle.h" +#include "larcore/CoreUtils/ServiceUtil.h" + +#include +#include + +#include "StitchingUtils.h" class IcarusHDF5Maker : public art::EDAnalyzer { public: @@ -66,164 +60,237 @@ class IcarusHDF5Maker : public art::EDAnalyzer { private: - std::string fTruthLabel; - std::string fHitLabel; - std::string fHitTruthLabel; - std::string fSPLabel; - std::string fOpHitLabel; - std::string fOpFlashLabel; - std::string fPFParticleLabel; + art::InputTag fTruthLabel; + art::InputTag fHitLabel; + art::InputTag fHitTruthLabel; + art::InputTag fSPLabel; + art::InputTag fOpHitLabel; + art::InputTag fOpFlashLabel; bool fUseMap; std::string fEventInfo; std::string fOutputName; - hep_hpc::hdf5::File fFile; ///< output HDF5 file - - hep_hpc::hdf5::Ntuple // event id (run, subrun, event) - >* fEventNtuple; ///< event ntuple - - hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // is cc - hep_hpc::hdf5::Column, // nu pdg - hep_hpc::hdf5::Column, // nu energy - hep_hpc::hdf5::Column, // lep energy - hep_hpc::hdf5::Column // nu dir (x, y, z) - >* fEventNtupleNu; ///< event ntuple with neutrino information - - hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // spacepoint id - hep_hpc::hdf5::Column, // 3d position (x, y, z) - hep_hpc::hdf5::Column, // 2d hit (u, v, y) - hep_hpc::hdf5::Column //ChiSquared of the hit - >* fSpacePointNtuple; ///< spacepoint ntuple - - hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // hit id - hep_hpc::hdf5::Column, // hit integral - hep_hpc::hdf5::Column, // hit rms - hep_hpc::hdf5::Column, // tpc id - hep_hpc::hdf5::Column, // global plane - hep_hpc::hdf5::Column, // global wire - hep_hpc::hdf5::Column, // global time - hep_hpc::hdf5::Column, // raw plane - hep_hpc::hdf5::Column, // raw wire - hep_hpc::hdf5::Column, // raw time - hep_hpc::hdf5::Column //cryostat - >* fHitNtuple; ///< hit ntuple - - - hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // g4 id - hep_hpc::hdf5::Column, // particle type - hep_hpc::hdf5::Column, // parent g4 id - hep_hpc::hdf5::Column, // is from nu - hep_hpc::hdf5::Column, // momentum - hep_hpc::hdf5::Column, // start position (x, y, z) - hep_hpc::hdf5::Column, // end position (x, y, z) - hep_hpc::hdf5::Column, // start process - hep_hpc::hdf5::Column // end process - >* fParticleNtuple; ///< particle ntuple - - hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // hit id - hep_hpc::hdf5::Column, // g4 id - hep_hpc::hdf5::Column, // deposited energy [ MeV ] - hep_hpc::hdf5::Column, // x position - hep_hpc::hdf5::Column, // y position - hep_hpc::hdf5::Column // z position - >* fEnergyDepNtuple; ///< energy deposition ntuple - - hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // hit id - hep_hpc::hdf5::Column, // hit_channel - hep_hpc::hdf5::Column, // wire pos - hep_hpc::hdf5::Column, // peaktime - hep_hpc::hdf5::Column, // width - hep_hpc::hdf5::Column, // area - hep_hpc::hdf5::Column, // amplitude - hep_hpc::hdf5::Column, // pe - hep_hpc::hdf5::Column // usedInFlash - >* fOpHitNtuple; ///< PMT hit ntuple - - hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // flash id - hep_hpc::hdf5::Column, // wire pos - hep_hpc::hdf5::Column, // time - hep_hpc::hdf5::Column, // time width - hep_hpc::hdf5::Column, // Y center - hep_hpc::hdf5::Column, // Y width - hep_hpc::hdf5::Column, // Z center - hep_hpc::hdf5::Column, // Z width - hep_hpc::hdf5::Column // totalpe - >* fOpFlashNtuple; ///< Flash ntuple - -hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // sumpe id - hep_hpc::hdf5::Column, // flash id - hep_hpc::hdf5::Column, // PMT channel - hep_hpc::hdf5::Column // pe - >* fOpFlashSumPENtuple; ///< Flash SumPE ntuple - - - - - - using ProxyPfpColl_t = decltype(proxy::getCollection >( - std::declval(),std::declval(), - proxy::withAssociated(std::declval()), - proxy::withAssociated(std::declval()), - proxy::withAssociated(std::declval()), - proxy::withAssociated(std::declval()) )); - using ProxyPfpElem_t = ProxyPfpColl_t::element_proxy_t; - - // proxy to connect cluster to hit - using ProxyClusColl_t = decltype(proxy::getCollection>( - std::declval(), std::declval(), - proxy::withAssociated(std::declval()))); - using ProxyClusElem_t = ProxyClusColl_t::element_proxy_t; - - std::vector nearwires; - int NearWire(const geo::WireReadoutGeom& geo, const geo::PlaneID &ip, const float x, const float y, const float z); - - int stitchedPlane(const geo::WireID wid) const { - int plane = wid.Plane; - if(wid.TPC==2 || wid.TPC==3) { - if(wid.Plane==1) plane=2; - else if(wid.Plane==2) plane=1; - } - return plane; - } - float stitchedTime(const geo::WireID wid, float timein) const { - float time = timein; - if(wid.TPC==2 || wid.TPC==3) { - //correction = 2*(tpcgeo.DriftDistance()/detProp.DriftVelocity()-clockData.TriggerOffsetTPC())/clockData.TPCClock().TickPeriod() = 6442.15 - time = 6442.15 - timein; - } - return time; - } - size_t stitchedWire(const geo::WireID wid) const { - size_t wire = wid.Wire; - int plane = stitchedPlane(wid); - if(wid.TPC==1 || wid.TPC==3) { - if(plane==1 || plane == 2) { - wire = wid.Wire + 2535; //2535 is the last part of the wires in cryos 0 an 2 before the cut in z=0 - } + std::vector tpc_ids_checked = {}; // add the TPC ids here so we only check them once + + struct HDFDataFile { + hep_hpc::hdf5::File file; ///< output HDF5 file + + std::optional // event id (run, subrun, event) + >> eventNtuple; ///< event ntuple + + std::optional, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // is cc + hep_hpc::hdf5::Column, // nu pdg + hep_hpc::hdf5::Column, // nu energy + hep_hpc::hdf5::Column, // lep energy + hep_hpc::hdf5::Column // nu dir (x, y, z) + >> eventNtupleNu; ///< event ntuple with neutrino information + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // spacepoint id + hep_hpc::hdf5::Column, // 3d position (x, y, z) + hep_hpc::hdf5::Column, // 2d hit (u, v, y) + hep_hpc::hdf5::Column //ChiSquared of the hit + > spacePointNtuple; ///< spacepoint ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // hit integral + hep_hpc::hdf5::Column, // hit rms + hep_hpc::hdf5::Column, // tpc id + hep_hpc::hdf5::Column, // global plane + hep_hpc::hdf5::Column, // global wire + hep_hpc::hdf5::Column, // global time + hep_hpc::hdf5::Column, // raw plane + hep_hpc::hdf5::Column, // raw wire + hep_hpc::hdf5::Column, // raw time + hep_hpc::hdf5::Column //cryostat + > hitNtuple; ///< hit ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // g4 id + hep_hpc::hdf5::Column, // particle type + hep_hpc::hdf5::Column, // parent g4 id + hep_hpc::hdf5::Column, // is from nu + hep_hpc::hdf5::Column, // momentum + hep_hpc::hdf5::Column, // start position (x, y, z) + hep_hpc::hdf5::Column, // end position (x, y, z) + hep_hpc::hdf5::Column, // start process + hep_hpc::hdf5::Column // end process + > particleNtuple; ///< particle ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // g4 id + hep_hpc::hdf5::Column, // deposited energy [ MeV ] + hep_hpc::hdf5::Column, // x position + hep_hpc::hdf5::Column, // y position + hep_hpc::hdf5::Column // z position + > energyDepNtuple; ///< energy deposition ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // hit_channel + hep_hpc::hdf5::Column, // wire pos + hep_hpc::hdf5::Column, // peaktime + hep_hpc::hdf5::Column, // width + hep_hpc::hdf5::Column, // area + hep_hpc::hdf5::Column, // amplitude + hep_hpc::hdf5::Column, // pe + hep_hpc::hdf5::Column // usedInFlash + > opHitNtuple; ///< PMT hit ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // flash id + hep_hpc::hdf5::Column, // wire pos + hep_hpc::hdf5::Column, // time + hep_hpc::hdf5::Column, // time width + hep_hpc::hdf5::Column, // Y center + hep_hpc::hdf5::Column, // Y width + hep_hpc::hdf5::Column, // Z center + hep_hpc::hdf5::Column, // Z width + hep_hpc::hdf5::Column // totalpe + > opFlashNtuple; ///< Flash ntuple + + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // sumpe id + hep_hpc::hdf5::Column, // flash id + hep_hpc::hdf5::Column, // PMT channel + hep_hpc::hdf5::Column // pe + > opFlashSumPENtuple; ///< Flash SumPE ntuple + + static std::string makeOutputFileName(std::string const& outputName, art::SubRunID const& sr) + { + struct timeval now; + gettimeofday(&now, NULL); + std::ostringstream fileName; + fileName << outputName + << "_r" << std::setfill('0') << std::setw(5) << sr.run() + << "_s" << std::setfill('0') << std::setw(5) << sr.subRun() + << "_ts" << std::setw(6) << now.tv_usec << ".h5"; + return fileName.str(); } - return wire; - } - + + HDFDataFile(std::string const& outputName, std::string const& eventInfo, art::SubRunID const& sr) + : file{ makeOutputFileName(outputName, sr), H5F_ACC_TRUNC } + , eventNtuple{ (eventInfo == "none") + ? std::optional > + >(hep_hpc::hdf5::make_ntuple({file, "event_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3)) + ) : std::nullopt + } + , eventNtupleNu{ (eventInfo == "nu") + ? std::optional, + hep_hpc::hdf5::Column, + hep_hpc::hdf5::Column, + hep_hpc::hdf5::Column, + hep_hpc::hdf5::Column, + hep_hpc::hdf5::Column > + >(hep_hpc::hdf5::make_ntuple({file, "event_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("is_cc"), + hep_hpc::hdf5::make_scalar_column("nu_pdg"), + hep_hpc::hdf5::make_scalar_column("nu_energy"), + hep_hpc::hdf5::make_scalar_column("lep_energy"), + hep_hpc::hdf5::make_column("nu_dir", 3)) + ) : std::nullopt + } + , spacePointNtuple{ + hep_hpc::hdf5::make_ntuple({file, "spacepoint_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("spacepoint_id"), + hep_hpc::hdf5::make_column("position", 3), + hep_hpc::hdf5::make_column("hit_id", 3), + hep_hpc::hdf5::make_scalar_column("Chi_squared")) + } + , hitNtuple{ + hep_hpc::hdf5::make_ntuple({file, "hit_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("hit_id"), + hep_hpc::hdf5::make_scalar_column("integral"), + hep_hpc::hdf5::make_scalar_column("rms"), + hep_hpc::hdf5::make_scalar_column("tpc"), + hep_hpc::hdf5::make_scalar_column("global_plane"), + hep_hpc::hdf5::make_scalar_column("global_wire"), + hep_hpc::hdf5::make_scalar_column("global_time"), + hep_hpc::hdf5::make_scalar_column("local_plane"), + hep_hpc::hdf5::make_scalar_column("local_wire"), + hep_hpc::hdf5::make_scalar_column("local_time"), + hep_hpc::hdf5::make_scalar_column("Cryostat")) + } + , particleNtuple{ + hep_hpc::hdf5::make_ntuple({file, "particle_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("g4_id"), + hep_hpc::hdf5::make_scalar_column("type"), + hep_hpc::hdf5::make_scalar_column("parent_id"), + hep_hpc::hdf5::make_scalar_column("from_nu"), + hep_hpc::hdf5::make_scalar_column("momentum"), + hep_hpc::hdf5::make_column("start_position", 3), + hep_hpc::hdf5::make_column("end_position", 3), + hep_hpc::hdf5::make_scalar_column("start_process"), + hep_hpc::hdf5::make_scalar_column("end_process")) + } + , energyDepNtuple{ + hep_hpc::hdf5::make_ntuple({file, "edep_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("hit_id"), + hep_hpc::hdf5::make_scalar_column("g4_id"), + hep_hpc::hdf5::make_scalar_column("energy"), + hep_hpc::hdf5::make_scalar_column("x_position"), + hep_hpc::hdf5::make_scalar_column("y_position"), + hep_hpc::hdf5::make_scalar_column("z_position")) + } + , opHitNtuple{ + hep_hpc::hdf5::make_ntuple({file, "ophit_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("hit_id"), + hep_hpc::hdf5::make_scalar_column("hit_channel"), + hep_hpc::hdf5::make_column("wire_pos", 3),//3 views + hep_hpc::hdf5::make_scalar_column("peaktime"), + hep_hpc::hdf5::make_scalar_column("width"), + hep_hpc::hdf5::make_scalar_column("area"), + hep_hpc::hdf5:: make_scalar_column("amplitude"), + hep_hpc::hdf5::make_scalar_column("pe"), + hep_hpc::hdf5::make_scalar_column("sumpe_id")) + } + , opFlashNtuple{ + hep_hpc::hdf5::make_ntuple({file, "opflash_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("flash_id"), + hep_hpc::hdf5::make_column("wire_pos", 3),//3 views + hep_hpc::hdf5::make_scalar_column("time"), + hep_hpc::hdf5::make_scalar_column("time_width"), + hep_hpc::hdf5::make_scalar_column("y_center"), + hep_hpc::hdf5::make_scalar_column("y_width"), + hep_hpc::hdf5::make_scalar_column("z_center"), + hep_hpc::hdf5::make_scalar_column("z_width"), + hep_hpc::hdf5::make_scalar_column("totalpe")) + } + , opFlashSumPENtuple{ + hep_hpc::hdf5::make_ntuple({file, "opflashsumpe_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("sumpe_id"), + hep_hpc::hdf5::make_scalar_column("flash_id"), + hep_hpc::hdf5::make_scalar_column("pmt_channel"), + hep_hpc::hdf5::make_scalar_column("sumpe")) + } + {} + }; + std::unique_ptr fHDFData; + + int NearWire(const geo::PlaneGeo &plane, float x, float y, float z) const; }; IcarusHDF5Maker::IcarusHDF5Maker(fhicl::ParameterSet const& p) : EDAnalyzer{p}, - fTruthLabel(p.get("TruthLabel")), - fHitLabel( p.get("HitLabel")), - fHitTruthLabel( p.get("HitTruthLabel","")), - fSPLabel( p.get("SPLabel")), - fOpHitLabel( p.get("OpHitLabel")), - fOpFlashLabel( p.get("OpFlashLabel")), - fPFParticleLabel( p.get("PFParticleLabel")), + fTruthLabel(p.get("TruthLabel")), + fHitLabel( p.get("HitLabel")), + fHitTruthLabel( p.get("HitTruthLabel","")), + fSPLabel( p.get("SPLabel")), + fOpHitLabel( p.get("OpHitLabel")), + fOpFlashLabel( p.get("OpFlashLabel")), fUseMap( p.get("UseMap", false)), fEventInfo( p.get("EventInfo")), fOutputName(p.get("OutputName")) @@ -232,19 +299,14 @@ IcarusHDF5Maker::IcarusHDF5Maker(fhicl::ParameterSet const& p) throw art::Exception(art::errors::Configuration) << "EventInfo must be \"none\" or \"nu\", not " << fEventInfo; } -std::vector tpc_ids_checked = {}; // add the TPC ids here so we only check them once void IcarusHDF5Maker::analyze(art::Event const& e) { - const cheat::BackTrackerService* bt = 0; - if (!fUseMap) { - art::ServiceHandle bt_h; - bt = bt_h.get(); - } - geo::WireReadoutGeom const& geo = art::ServiceHandle()->Get(); + const cheat::BackTrackerService* bt = fUseMap? nullptr: art::ServiceHandle().get(); + geo::WireReadoutGeom const& geom = art::ServiceHandle()->Get(); std::set g4id; - art::ServiceHandle pi; + auto const* pi = lar::providerFrom(); auto const clockData = art::ServiceHandle()->DataFor(e); auto const detProp = art::ServiceHandle()->DataFor(e, clockData); int run = e.id().run(); @@ -261,12 +323,12 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // std::cout << "triggerTime=" << clockData.TriggerTime() << std::endl; // std::cout << "correction=" << 2*(tpcgeo.DriftDistance()/detProp.DriftVelocity()-clockData.TriggerOffsetTPC())/clockData.TPCClock().TickPeriod() << std::endl; - std::cout << "NEW EVENT: " << run <<" "<< subrun << " "<< event< evtID { run, subrun, event }; // Fill event table if (fEventInfo == "none") { - fEventNtuple->insert( evtID.data() ); + fHDFData->eventNtuple->insert( evtID.data() ); mf::LogInfo("IcarusHDF5Maker") << "Filling event table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2]; @@ -274,24 +336,23 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { if (fEventInfo == "nu") { // Get MC truth - art::Handle> truthHandle; - e.getByLabel(fTruthLabel, truthHandle); - if (!truthHandle.isValid() || truthHandle->size() == 0) { + auto truthHandle = e.getValidHandle>(fTruthLabel); + if (truthHandle->size() != 1) { throw art::Exception(art::errors::LogicError) << "Expected to find exactly one MC truth object!"; } - simb::MCNeutrino nutruth = truthHandle->at(0).GetNeutrino(); + simb::MCNeutrino const& nutruth = truthHandle->at(0).GetNeutrino(); auto up = nutruth.Nu().Momentum().Vect().Unit(); std::array nuMomentum {(float)up.X(),(float)up.Y(),(float)up.Z()}; - fEventNtupleNu->insert( evtID.data(), - nutruth.CCNC() == simb::kCC, - nutruth.Nu().PdgCode(), - nutruth.Nu().E(), - nutruth.Lepton().E(), - nuMomentum.data() - ); + fHDFData->eventNtupleNu->insert( evtID.data(), + nutruth.CCNC() == simb::kCC, + nutruth.Nu().PdgCode(), + nutruth.Nu().E(), + nutruth.Lepton().E(), + nuMomentum.data() + ); // for (int ip=0;ipat(0).NParticles();ip++) { // std::cout << "mcp tkid=" << truthHandle->at(0).GetParticle(ip).TrackId() << " pdg=" << truthHandle->at(0).GetParticle(ip).PdgCode() @@ -300,25 +361,24 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // << std::endl; // } - mf::LogInfo("IcarusHDF5Maker") << "Filling event table" - << "\nrun " << evtID[0] << ", subrun " << evtID[1] - << ", event " << evtID[2] - << "\nis cc? " << (nutruth.CCNC() == simb::kCC) - << ", nu energy " << nutruth.Nu().E() - << ", lepton energy " << nutruth.Lepton().E() - << "\nnu momentum x " << nuMomentum[0] << ", y " - << nuMomentum[1] << ", z " << nuMomentum[2]; + mf::LogDebug("IcarusHDF5Maker") << "Filling event table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nis cc? " << (nutruth.CCNC() == simb::kCC) + << ", nu energy " << nutruth.Nu().E() + << ", lepton energy " << nutruth.Lepton().E() + << "\nnu momentum x " << nuMomentum[0] << ", y " + << nuMomentum[1] << ", z " << nuMomentum[2]; } // if nu event info // Get spacepoints from the event record - art::Handle> spListHandle; - std::vector> splist; - if (e.getByLabel(fSPLabel, spListHandle)) art::fill_ptr_vector(splist, spListHandle); + auto spListHandle = e.getValidHandle>(fSPLabel); + std::vector const& splist = *spListHandle; // Get hits from the event record - art::Handle> hitListHandle; + auto hitListHandle = e.getValidHandle>(fHitLabel); std::vector> hitlist; - if (e.getByLabel(fHitLabel, hitListHandle)) art::fill_ptr_vector(hitlist, hitListHandle); + art::fill_ptr_vector(hitlist, hitListHandle); // Get assocations from spacepoints to hits art::FindManyP fmp(spListHandle, e, fSPLabel); @@ -326,38 +386,34 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // Fill spacepoint table for (size_t i = 0; i < spListHandle->size(); ++i) { - art::Ptr spacePointPtr(spListHandle,i); - std::vector> associatedHits(fmp.at(spacePointPtr.key())); + std::vector> const& associatedHits(fmp.at(i)); if (associatedHits.size() < 3) { - std::cout << "I am certain this cannot happen... but here you go, space point with " << associatedHits.size() << " hits" << std::endl; - exit(1); - continue; + throw art::Exception(art::errors::LogicError) << "I am certain this cannot happen... but here you go, space point with " << associatedHits.size() << " hits"; } - std::array pos {(float)splist[i]->XYZ()[0],(float)splist[i]->XYZ()[1],(float)splist[i]->XYZ()[2]}; + std::array pos {(float)splist[i].XYZ()[0],(float)splist[i].XYZ()[1],(float)splist[i].XYZ()[2]}; std::array hitID { -1, -1, -1 }; - for (size_t j = 0; j < associatedHits.size(); ++j) { - int plane = stitchedPlane(associatedHits[j]->WireID()); - //std::cout << "j=" << j << " tpc=" << associatedHits[j]->WireID().TPC << " plane=" << associatedHits[j]->WireID().Plane - // << " plane2=" << plane << " key=" << associatedHits[j].key() << std::endl; - hitID[plane] = associatedHits[j].key(); + for (art::Ptr const& hit: associatedHits) { + int plane = util::stitchedPlane(hit->WireID()); + //std::cout << " tpc=" << hit->WireID().TPC << " plane=" << hit->WireID().Plane + // << " plane2=" << plane << " key=" << hit.key() << std::endl; + hitID[plane] = hit.key(); } if (hitID[0]==0 && hitID[1]==-1 && hitID[2]==-1) { - std::cout << "THIS SHOULD NOT HAPPEN -- sp i=" << i << " hit ids=" << hitID[0] << " " << hitID[1] << " " << hitID[2] << std::endl; - exit(1); + throw art::Exception(art::errors::LogicError) << "THIS SHOULD NOT HAPPEN -- sp i=" << i << " hit ids=" << hitID[0] << " " << hitID[1] << " " << hitID[2]; } - fSpacePointNtuple->insert(evtID.data(),splist[i]->ID(), pos.data(), hitID.data(),splist[i]->Chisq() ); - - mf::LogInfo("IcarusHDF5Maker") << "Filling spacepoint table" - << "\nrun " << evtID[0] << ", subrun " << evtID[1] - << ", event " << evtID[2] - << "\nspacepoint id " << splist[i]->ID() - << "\nposition x " << pos[0] << ", y " << pos[1] - << ", z " << pos[2] - << "\nhit ids " << hitID[0] << ", " << hitID[1] - << ", " << hitID[2] - << "Chi_squared " << splist[i]->Chisq(); + fHDFData->spacePointNtuple.insert(evtID.data(),splist[i].ID(), pos.data(), hitID.data(),splist[i].Chisq() ); + + mf::LogDebug("IcarusHDF5Maker") << "Filling spacepoint table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nspacepoint id " << splist[i].ID() + << "\nposition x " << pos[0] << ", y " << pos[1] + << ", z " << pos[2] + << "\nhit ids " << hitID[0] << ", " << hitID[1] + << ", " << hitID[2] + << "Chi_squared " << splist[i].Chisq(); } // for spacepoint std::unique_ptr> hittruth; @@ -370,14 +426,14 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // Fill hit table geo::WireID wireid = hit->WireID(); - size_t plane = stitchedPlane(wireid); - double time = stitchedTime(wireid,hit->PeakTime()); - size_t wire = stitchedWire(wireid); - fHitNtuple->insert(evtID.data(), - hit.key(), hit->Integral(), hit->RMS(), wireid.TPC, - plane, wire, time, - wireid.Plane, wireid.Wire, hit->PeakTime(),wireid.Cryostat - ); + size_t plane = util::stitchedPlane(wireid); + double time = util::stitchedTime(wireid,hit->PeakTime()); + size_t wire = util::stitchedWire(wireid); + fHDFData->hitNtuple.insert(evtID.data(), + hit.key(), hit->Integral(), hit->RMS(), wireid.TPC, + plane, wire, time, + wireid.Plane, wireid.Wire, hit->PeakTime(),wireid.Cryostat + ); mf::LogInfo("IcarusHDF5Maker") << "Filling hit table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] @@ -395,7 +451,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // Fill energy deposit table if (fUseMap) { //not supported in icarus at the moment - exit(1); + throw art::Exception{ art::errors::UnimplementedFeature } << "The use of map is currently not supported in ICARUS"; /* std::vector> particle_vec = hittruth->at(hit.key()); std::vector match_vec = hittruth->data(hit.key()); @@ -419,32 +475,35 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { */ } else { - std::vector h_ides = bt->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick()); - for (auto& tide : h_ides) { + std::vector const& h_ides = bt->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick()); + for (auto const& tide : h_ides) { int tid = tide.trackID; + // if negative id, make sure there is a corresponding particle to look for before taking the abs. This way negative means no asociated particle. if (pi->TrackIdToParticle_P(abs(tid))) tid = abs(tid); g4id.insert(tid); - fEnergyDepNtuple->insert(evtID.data(),hit.key(), tid, tide.energyFrac, -1., -1., -1.); + fHDFData->energyDepNtuple.insert(evtID.data(),hit.key(), tid, tide.energyFrac, -1., -1., -1.); } } // if using microboone map method or not } // for hit - std::set allIDs = g4id; // Copy original so we can safely modify it + std::unordered_map allIDs; // Add invisible particles to hierarchy for (int id : g4id) { const simb::MCParticle* p = pi->TrackIdToParticle_P(abs(id)); - while (p!= NULL && p->Mother() != 0 ) { - allIDs.insert(abs(p->Mother())); + allIDs.emplace(abs(id), p); + while (p && p->Mother() != 0 ) { p = pi->TrackIdToParticle_P(abs(p->Mother())); + allIDs.emplace(abs(p->Mother()), p); } } // Loop over true particles and fill table - for (int id : allIDs) { - const simb::MCParticle* p = pi->TrackIdToParticle_P(abs(id)); - if(p==NULL) std::cout<< "we have a problem p is null"<TrackIdToMCTruth_P(abs(id)); int fromNu = (mct.isAvailable() ? mct->NeutrinoSet() : 0); // std::cout << "all mcp tkid=" << p->TrackId() << " pdg=" << p->PdgCode() @@ -454,257 +513,130 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // << std::endl; std::array particleStart { (float)p->Vx(), (float)p->Vy(), (float)p->Vz() }; std::array particleEnd { (float)p->EndX(), (float)p->EndY(), (float)p->EndZ() }; - fParticleNtuple->insert(evtID.data(), - abs(id), p->PdgCode(), p->Mother(), fromNu, (float)p->P(), - particleStart.data(), particleEnd.data(), - p->Process(), p->EndProcess() - ); - mf::LogInfo("IcarusHDF5Maker") << "Filling particle table" - << "\nrun " << evtID[0] << ", subrun " << evtID[1] - << ", event " << evtID[2] - << "\ng4 id " << abs(id) << ", pdg code " - << p->PdgCode() << ", parent " << p->Mother() - << ", momentum " << p->P() - << "\nparticle start x " << particleStart[0] - << ", y " << particleStart[1] - << ", z " << particleStart[2] - << "\nparticle end x " << particleEnd[0] << ", y " - << particleEnd[1] << ", z " << particleEnd[2] - << "\nstart process " << p->Process() - << ", end process " << p->EndProcess(); + fHDFData->particleNtuple.insert(evtID.data(), + abs(id), p->PdgCode(), p->Mother(), fromNu, (float)p->P(), + particleStart.data(), particleEnd.data(), + p->Process(), p->EndProcess() + ); + mf::LogDebug("IcarusHDF5Maker") << "Filling particle table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\ng4 id " << abs(id) << ", pdg code " + << p->PdgCode() << ", parent " << p->Mother() + << ", momentum " << p->P() + << "\nparticle start x " << particleStart[0] + << ", y " << particleStart[1] + << ", z " << particleStart[2] + << "\nparticle end x " << particleEnd[0] << ", y " + << particleEnd[1] << ", z " << particleEnd[2] + << "\nstart process " << p->Process() + << ", end process " << p->EndProcess(); } // Get OpFlashs from the event record art::Handle< std::vector< recob::OpFlash > > opFlashListHandle; - std::unique_ptr > assocFlashHit; - std::vector< art::Ptr< recob::OpFlash > > opflashlist; - if (fOpFlashLabel != "") { - if (e.getByLabel(fOpFlashLabel, opFlashListHandle)) art::fill_ptr_vector(opflashlist, opFlashListHandle); - assocFlashHit = std::unique_ptr >(new art::FindManyP(opFlashListHandle, e, fOpFlashLabel)); + std::optional > assocFlashHit; + std::vector const* opflashlist = nullptr; + if (!fOpFlashLabel.empty()) { + e.getByLabel(fOpFlashLabel, opFlashListHandle); + opflashlist = opFlashListHandle.product(); + assocFlashHit.emplace(opFlashListHandle, e, fOpFlashLabel); } // get the flash matching the slice we selected - art::Handle< art::Assns > sliceOpFlashAssnsHandle; - e.getByLabel(fHitLabel, sliceOpFlashAssnsHandle); - int flkey = -1; - if (sliceOpFlashAssnsHandle->size()) { - art::Ptr fp = sliceOpFlashAssnsHandle->at(0).second; - flkey = fp.key(); - } + auto const& sliceOpFlashAssns = e.getProduct>(fHitLabel); + int flkey = sliceOpFlashAssns.size()==0 ? -1: sliceOpFlashAssns.at(0).second.key(); std::vector > flashsumpepmtmap;//this has at most one element, due to the check on flkey - // Loop over opflashs, pick only the one we care about - for (art::Ptr< recob::OpFlash > opflash : opflashlist) { - if (int(opflash.key())!=flkey) continue; + // Pick only the flash we care about + if (flkey>=0) { + recob::OpFlash const& opflash = opflashlist->at(flkey); std::vector pes; - for (auto pe : opflash->PEs()) pes.push_back(pe); + for (auto pe : opflash.PEs()) pes.push_back(pe); std::vector nearwires; - double xyz[3] = {0.,opflash->YCenter(),opflash->ZCenter()}; - for (auto const& p : geo.Iterate()) nearwires.push_back(NearWire(geo,p,xyz[0],xyz[1],xyz[2])); + double xyz[3] = {0.,opflash.YCenter(),opflash.ZCenter()}; + for (auto const& p : geom.Iterate()) nearwires.push_back(NearWire(p,xyz[0],xyz[1],xyz[2])); // Fill opflash table - fOpFlashNtuple->insert(evtID.data(), - opflash.key(), - nearwires.data(), - opflash->Time(),opflash->TimeWidth(), - opflash->YCenter(),opflash->YWidth(),opflash->ZCenter(),opflash->ZWidth(), - opflash->TotalPE() - ); + fHDFData->opFlashNtuple.insert(evtID.data(), + flkey, + nearwires.data(), + opflash.Time(),opflash.TimeWidth(), + opflash.YCenter(),opflash.YWidth(),opflash.ZCenter(),opflash.ZWidth(), + opflash.TotalPE() + ); size_t count = 0; - std::vector sumpepmtmap(art::ServiceHandle()->NOpDets(),-1); + std::vector sumpepmtmap(art::ServiceHandle()->NOpDets(),0); for (size_t ipmt=0;ipmtinsert(evtID.data(),count,opflash.key(),ipmt,pes[ipmt]); + fHDFData->opFlashSumPENtuple.insert(evtID.data(),count,flkey,ipmt,pes[ipmt]); sumpepmtmap[ipmt] = count; count++; } flashsumpepmtmap.push_back(sumpepmtmap); - mf::LogInfo("IcarusHDF5Maker") << "Filling opflash table" - << "\nrun " << evtID[0] << ", subrun " << evtID[1] - << ", event " << evtID[2] - << "\nflash id " << opflash.key() << ", Time " << opflash->Time() - << ", TotalPE " << opflash->TotalPE()//; - << "\nWireCenters size0 " << opflash->WireCenters().size()//; - << "\nYCenter " << opflash->YCenter()<< " ZCenter " << opflash->ZCenter() - << "\nYWidth " << opflash->YWidth()<< " ZWidth " << opflash->ZWidth() - << "\nInBeamFrame " << opflash->InBeamFrame()<< " OnBeamTime " << opflash->OnBeamTime() - << "\nAbsTime " << opflash->AbsTime() << " TimeWidth " << opflash->TimeWidth() << " FastToTotal " << opflash->FastToTotal(); + mf::LogDebug("IcarusHDF5Maker") << "Filling opflash table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nflash id " << flkey << ", Time " << opflash.Time() + << ", TotalPE " << opflash.TotalPE()//; + << "\nWireCenters size0 " << opflash.WireCenters().size()//; + << "\nYCenter " << opflash.YCenter()<< " ZCenter " << opflash.ZCenter() + << "\nYWidth " << opflash.YWidth()<< " ZWidth " << opflash.ZWidth() + << "\nInBeamFrame " << opflash.InBeamFrame()<< " OnBeamTime " << opflash.OnBeamTime() + << "\nAbsTime " << opflash.AbsTime() << " TimeWidth " << opflash.TimeWidth() << " FastToTotal " << opflash.FastToTotal(); } // Get OpHits from the event record - art::Handle< std::vector< recob::OpHit > > opHitListHandle; std::vector< art::Ptr< recob::OpHit > > ophitlist; - if (fOpHitLabel != "") { - if (e.getByLabel(fOpHitLabel, opHitListHandle)) - art::fill_ptr_vector(ophitlist, opHitListHandle); + if (!fOpHitLabel.empty()) { + auto opHitListHandle = e.getValidHandle< std::vector< recob::OpHit > >(fOpHitLabel); + art::fill_ptr_vector(ophitlist, opHitListHandle); } + std::set> ophv; + if (flkey >= 0) { + auto const& flashHits = assocFlashHit->at(flkey); + ophv.insert(begin(flashHits), end(flashHits)); + } // Loop over ophits for (art::Ptr< recob::OpHit > ophit : ophitlist) { std::vector nearwires; - auto xyz = geo.OpDetGeoFromOpChannel(ophit->OpChannel()).GetCenter(); - for (auto const& p : geo.Iterate()) nearwires.push_back(NearWire(geo,p,xyz.X(),xyz.Y(),xyz.Z())); + auto xyz = geom.OpDetGeoFromOpChannel(ophit->OpChannel()).GetCenter(); + for (auto const& p : geom.Iterate()) nearwires.push_back(NearWire(p,xyz.X(),xyz.Y(),xyz.Z())); + + bool isInFlash = (ophv.count(ophit) > 0); - bool isInFlash = false; - if (flkey>=0) { - auto ophv = assocFlashHit->at(flkey); - isInFlash = (std::find(ophv.begin(),ophv.end(),ophit)!=ophv.end()); - } // Fill ophit table - fOpHitNtuple->insert(evtID.data(),ophit.key(), - ophit->OpChannel(),nearwires.data(), - ophit->PeakTime(),ophit->Width(), - ophit->Area(),ophit->Amplitude(),ophit->PE(), - (isInFlash ? flashsumpepmtmap[0][ophit->OpChannel()] : -1) - ); - mf::LogInfo("IcarusHDF5Maker") << "\nFilling ophit table" - << "\nrun " << evtID[0] << ", subrun " << evtID[1] - << ", event " << evtID[2] - << "\nhit id " << ophit.key() << ", channel " - << ophit->OpChannel() << ", PeakTime " << ophit->PeakTime() - << ", Width " << ophit->Width() - << "\narea " << ophit->Area() << ", amplitude " - << ophit->Amplitude() << ", PE " << ophit->PE(); + fHDFData->opHitNtuple.insert(evtID.data(),ophit.key(), + ophit->OpChannel(),nearwires.data(), + ophit->PeakTime(),ophit->Width(), + ophit->Area(),ophit->Amplitude(),ophit->PE(), + (isInFlash ? flashsumpepmtmap[0][ophit->OpChannel()] : -1) + ); + mf::LogDebug("IcarusHDF5Maker") << "\nFilling ophit table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nhit id " << ophit.key() << ", channel " + << ophit->OpChannel() << ", PeakTime " << ophit->PeakTime() + << ", Width " << ophit->Width() + << "\narea " << ophit->Area() << ", amplitude " + << ophit->Amplitude() << ", PE " << ophit->PE(); } //End optical analyzer } // function IcarusHDF5Maker::analyze void IcarusHDF5Maker::beginSubRun(art::SubRun const& sr) { - - struct timeval now; - gettimeofday(&now, NULL); - // Open HDF5 output - std::ostringstream fileName; - fileName << fOutputName << "_r" << std::setfill('0') << std::setw(5) << sr.run() - << "_s" << std::setfill('0') << std::setw(5) << sr.subRun() << "_ts" << std::setw(6) << now.tv_usec << ".h5"; - - fFile = hep_hpc::hdf5::File(fileName.str(), H5F_ACC_TRUNC); - - if (fEventInfo == "none") - fEventNtuple = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "event_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3) - )); - if (fEventInfo == "nu") - fEventNtupleNu = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "event_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("is_cc"), - hep_hpc::hdf5::make_scalar_column("nu_pdg"), - hep_hpc::hdf5::make_scalar_column("nu_energy"), - hep_hpc::hdf5::make_scalar_column("lep_energy"), - hep_hpc::hdf5::make_column("nu_dir", 3) - )); - - fSpacePointNtuple = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "spacepoint_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("spacepoint_id"), - hep_hpc::hdf5::make_column("position", 3), - hep_hpc::hdf5::make_column("hit_id", 3), - hep_hpc::hdf5::make_scalar_column("Chi_squared") - - )); - - fHitNtuple = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "hit_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("hit_id"), - hep_hpc::hdf5::make_scalar_column("integral"), - hep_hpc::hdf5::make_scalar_column("rms"), - hep_hpc::hdf5::make_scalar_column("tpc"), - hep_hpc::hdf5::make_scalar_column("global_plane"), - hep_hpc::hdf5::make_scalar_column("global_wire"), - hep_hpc::hdf5::make_scalar_column("global_time"), - hep_hpc::hdf5::make_scalar_column("local_plane"), - hep_hpc::hdf5::make_scalar_column("local_wire"), - hep_hpc::hdf5::make_scalar_column("local_time"), - hep_hpc::hdf5::make_scalar_column("Cryostat") - - )); - - fParticleNtuple = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "particle_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("g4_id"), - hep_hpc::hdf5::make_scalar_column("type"), - hep_hpc::hdf5::make_scalar_column("parent_id"), - hep_hpc::hdf5::make_scalar_column("from_nu"), - hep_hpc::hdf5::make_scalar_column("momentum"), - hep_hpc::hdf5::make_column("start_position", 3), - hep_hpc::hdf5::make_column("end_position", 3), - hep_hpc::hdf5::make_scalar_column("start_process"), - hep_hpc::hdf5::make_scalar_column("end_process") - )); - - fEnergyDepNtuple = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "edep_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("hit_id"), - hep_hpc::hdf5::make_scalar_column("g4_id"), - hep_hpc::hdf5::make_scalar_column("energy"), - hep_hpc::hdf5::make_scalar_column("x_position"), - hep_hpc::hdf5::make_scalar_column("y_position"), - hep_hpc::hdf5::make_scalar_column("z_position") - )); - - fOpHitNtuple = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "ophit_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("hit_id"), - hep_hpc::hdf5::make_scalar_column("hit_channel"), - hep_hpc::hdf5::make_column("wire_pos", 3),//3 views - hep_hpc::hdf5::make_scalar_column("peaktime"), - hep_hpc::hdf5::make_scalar_column("width"), - hep_hpc::hdf5::make_scalar_column("area"), - hep_hpc::hdf5:: make_scalar_column("amplitude"), - hep_hpc::hdf5::make_scalar_column("pe"), - hep_hpc::hdf5::make_scalar_column("sumpe_id") - )); - - fOpFlashNtuple = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "opflash_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("flash_id"), - hep_hpc::hdf5::make_column("wire_pos", 3),//3 views - hep_hpc::hdf5::make_scalar_column("time"), - hep_hpc::hdf5::make_scalar_column("time_width"), - hep_hpc::hdf5::make_scalar_column("y_center"), - hep_hpc::hdf5::make_scalar_column("y_width"), - hep_hpc::hdf5::make_scalar_column("z_center"), - hep_hpc::hdf5::make_scalar_column("z_width"), - hep_hpc::hdf5::make_scalar_column("totalpe") - )); - - fOpFlashSumPENtuple = new hep_hpc::hdf5::Ntuple( - hep_hpc::hdf5::make_ntuple({fFile, "opflashsumpe_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("sumpe_id"), - hep_hpc::hdf5::make_scalar_column("flash_id"), - hep_hpc::hdf5::make_scalar_column("pmt_channel"), - hep_hpc::hdf5::make_scalar_column("sumpe") - )); + fHDFData = std::make_unique(fOutputName, fEventInfo, sr.id()); } void IcarusHDF5Maker::endSubRun(art::SubRun const& sr) { - if (fEventInfo == "none") delete fEventNtuple; - if (fEventInfo == "nu") delete fEventNtupleNu; - delete fSpacePointNtuple; - delete fHitNtuple; - delete fParticleNtuple; - delete fEnergyDepNtuple; - delete fOpHitNtuple; - delete fOpFlashNtuple; - delete fOpFlashSumPENtuple; - fFile.close(); + fHDFData.reset(); } -int IcarusHDF5Maker::NearWire(const geo::WireReadoutGeom& geo, const geo::PlaneID &ip, const float x, const float y, const float z) -//art::ServiceHandle geo; +int IcarusHDF5Maker::NearWire(const geo::PlaneGeo &plane, float x, float y, float z) const { - geo::PlaneGeo const& plane = geo.Plane(ip); geo::WireID wireID; try { wireID = plane.NearestWireID(geo::Point_t(x,y,z)); diff --git a/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc b/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc index dedb66b8d..e3fea33f6 100644 --- a/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc +++ b/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc @@ -1,13 +1,24 @@ +/** + * @file icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc + * @author Giuseppe Cerati (cerati@fnal.gov) + */ + #include "larrecodnn/NuGraph/Tools/LoaderToolBase.h" #include "art/Utilities/ToolMacros.h" +#include "messagefacility/MessageLogger/MessageLogger.h" #include "canvas/Persistency/Common/FindManyP.h" #include "canvas/Persistency/Common/Ptr.h" #include "canvas/Utilities/InputTag.h" #include "lardataobj/RecoBase/SpacePoint.h" +#include "lardataobj/RecoBase/Hit.h" +#include "larcoreobj/SimpleTypesAndConstants/geo_types.h" // geo::WireID +#include // std::move() #include +#include "StitchingUtils.h" + class IcarusNuGraphLoader : public LoaderToolBase { public: @@ -18,11 +29,6 @@ class IcarusNuGraphLoader : public LoaderToolBase { */ IcarusNuGraphLoader(const fhicl::ParameterSet& pset); - /** - * @brief Virtual Destructor - */ - virtual ~IcarusNuGraphLoader() noexcept = default; - /** * @brief loadData function * @@ -33,40 +39,14 @@ class IcarusNuGraphLoader : public LoaderToolBase { vector& inputs, vector>& idsmap) override; - int stitchedPlane(const geo::WireID wid) const { - int plane = wid.Plane; - if(wid.TPC==2 || wid.TPC==3) { - if(wid.Plane==1) plane=2; - else if(wid.Plane==2) plane=1; - } - return plane; - } - float stitchedTime(const geo::WireID wid, float timein) const { - float time = timein; - if(wid.TPC==2 || wid.TPC==3) { - //correction = 2*(tpcgeo.DriftDistance()/detProp.DriftVelocity()-clockData.TriggerOffsetTPC())/clockData.TPCClock().TickPeriod() = 6442.15 - time = 6442.15 - timein; - } - return time; - } - size_t stitchedWire(const geo::WireID wid) const { - size_t wire = wid.Wire; - int plane = stitchedPlane(wid); - if(wid.TPC==1 || wid.TPC==3) { - if(plane==1 || plane == 2) { - wire = wid.Wire + 2535; //2535 is the last part of the wires in cryos 0 an 2 before the cut in z=0 - } - } - return wire; - } - private: - art::InputTag hitInput; - art::InputTag spsInput; + + art::InputTag const fHitInput; + art::InputTag const fSpsInput; }; IcarusNuGraphLoader::IcarusNuGraphLoader(const fhicl::ParameterSet& p) - : hitInput{p.get("hitInput")}, spsInput{p.get("spsInput")} + : fHitInput{p.get("hitInput")}, fSpsInput{p.get("spsInput")} {} void IcarusNuGraphLoader::loadData(art::Event& e, @@ -75,13 +55,13 @@ void IcarusNuGraphLoader::loadData(art::Event& e, vector>& idsmap) { // - art::Handle> hitListHandle; - if (e.getByLabel(hitInput, hitListHandle)) { art::fill_ptr_vector(hitlist, hitListHandle); } + auto hitListHandle = e.getValidHandle>(fHitInput); + art::fill_ptr_vector(hitlist, hitListHandle); // idsmap = std::vector>(planes.size(), std::vector()); for (auto h : hitlist) { // - size_t plane = stitchedPlane(h->WireID()); + size_t plane = util::stitchedPlane(h->WireID()); idsmap[plane].push_back(h.key()); } @@ -100,9 +80,9 @@ void IcarusNuGraphLoader::loadData(art::Event& e, for (auto h : hitlist) { geo::WireID wireid = h->WireID(); - size_t plane = stitchedPlane(wireid); - double time = stitchedTime(wireid,h->PeakTime()); - size_t wire = stitchedWire(wireid); + size_t plane = util::stitchedPlane(wireid); + double time = util::stitchedTime(wireid,h->PeakTime()); + size_t wire = util::stitchedWire(wireid); hit_table_hit_id_data.push_back(h.key()); hit_table_local_plane_data.push_back(plane); @@ -111,18 +91,20 @@ void IcarusNuGraphLoader::loadData(art::Event& e, hit_table_integral_data.push_back(h->Integral()); hit_table_rms_data.push_back(h->RMS()); } - std::cout << "loader has nhits=" << hit_table_hit_id_data.size() << std::endl; + mf::LogDebug{ "IcarusNuGraphLoader" } << "loader has nhits=" << hit_table_hit_id_data.size(); // Get spacepoints from the event record art::Handle> spListHandle; vector> splist; - if (e.getByLabel(spsInput, spListHandle)) { art::fill_ptr_vector(splist, spListHandle); } + if (e.getByLabel(fSpsInput, spListHandle)) { art::fill_ptr_vector(splist, spListHandle); } // Get assocations from spacepoints to hits vector>> sp2Hit(splist.size()); - if (splist.size() > 0) { - art::FindManyP fmp(spListHandle, e, spsInput); + if (!splist.empty()) { + art::FindManyP fmp(spListHandle, e, fSpsInput); for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) { - if (splist[spIdx]->Chisq()>0.5) continue; + static constexpr double MinChiSq = 0.5; ///< Threshold to consider a space point good. + if (splist[spIdx]->Chisq()>MinChiSq) continue; + // only space points with hits on all planes are enough for NuGraph if (fmp.at(spIdx).size()!=3) continue; sp2Hit[spIdx] = fmp.at(spIdx); } @@ -131,31 +113,31 @@ void IcarusNuGraphLoader::loadData(art::Event& e, // space point table size_t spidx = 0; for (size_t i = 0; i < splist.size(); ++i) { - if (sp2Hit[i].size()==0) continue; + if (sp2Hit[i].empty()) continue; spacepoint_table_spacepoint_id_data.push_back(spidx++); spacepoint_table_hit_id_u_data.push_back(-1); spacepoint_table_hit_id_v_data.push_back(-1); spacepoint_table_hit_id_y_data.push_back(-1); for (size_t j = 0; j < sp2Hit[i].size(); ++j) { // - size_t plane = stitchedPlane(sp2Hit[i][j]->WireID()); + size_t plane = util::stitchedPlane(sp2Hit[i][j]->WireID()); if (plane == 0) spacepoint_table_hit_id_u_data.back() = sp2Hit[i][j].key(); if (plane == 1) spacepoint_table_hit_id_v_data.back() = sp2Hit[i][j].key(); if (plane == 2) spacepoint_table_hit_id_y_data.back() = sp2Hit[i][j].key(); } } - std::cout << "loader has nsps=" << spacepoint_table_hit_id_u_data.size() << std::endl; - - inputs.emplace_back("hit_table_hit_id", hit_table_hit_id_data); - inputs.emplace_back("hit_table_local_plane", hit_table_local_plane_data); - inputs.emplace_back("hit_table_local_time", hit_table_local_time_data); - inputs.emplace_back("hit_table_local_wire", hit_table_local_wire_data); - inputs.emplace_back("hit_table_integral", hit_table_integral_data); - inputs.emplace_back("hit_table_rms", hit_table_rms_data); - - inputs.emplace_back("spacepoint_table_spacepoint_id", spacepoint_table_spacepoint_id_data); - inputs.emplace_back("spacepoint_table_hit_id_u", spacepoint_table_hit_id_u_data); - inputs.emplace_back("spacepoint_table_hit_id_v", spacepoint_table_hit_id_v_data); - inputs.emplace_back("spacepoint_table_hit_id_y", spacepoint_table_hit_id_y_data); + mf::LogDebug{ "IcarusNuGraphLoader" } << "loader has nsps=" << spacepoint_table_hit_id_u_data.size(); + + inputs.emplace_back("hit_table_hit_id", std::move(hit_table_hit_id_data)); + inputs.emplace_back("hit_table_local_plane", std::move(hit_table_local_plane_data)); + inputs.emplace_back("hit_table_local_time", std::move(hit_table_local_time_data)); + inputs.emplace_back("hit_table_local_wire", std::move(hit_table_local_wire_data)); + inputs.emplace_back("hit_table_integral", std::move(hit_table_integral_data)); + inputs.emplace_back("hit_table_rms", std::move(hit_table_rms_data)); + + inputs.emplace_back("spacepoint_table_spacepoint_id", std::move(spacepoint_table_spacepoint_id_data)); + inputs.emplace_back("spacepoint_table_hit_id_u", std::move(spacepoint_table_hit_id_u_data)); + inputs.emplace_back("spacepoint_table_hit_id_v", std::move(spacepoint_table_hit_id_v_data)); + inputs.emplace_back("spacepoint_table_hit_id_y", std::move(spacepoint_table_hit_id_y_data)); } DEFINE_ART_CLASS_TOOL(IcarusNuGraphLoader) diff --git a/icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc b/icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc index 7624ac35b..7c521d5e6 100644 --- a/icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc +++ b/icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc @@ -1,26 +1,26 @@ -//////////////////////////////////////////////////////////////////////// -// Class: IcarusNuSliceHitsProducer -// Plugin Type: producer (art v3_06_03) -// File: IcarusNuSliceHitsProducer_module.cc -// -// Generated at Tue May 25 10:39:19 2021 by Giuseppe Cerati using cetskelgen -// from cetlib version v3_11_01. -//////////////////////////////////////////////////////////////////////// +/** + * @file icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc + * @brief Implementation of `IcarusNuSliceHitsProducer` _art_ module. + * @author Giuseppe Cerati (cerati@fnal.gov) + * @date May 25, 2021 + */ #include "art/Framework/Core/EDProducer.h" #include "art/Framework/Core/ModuleMacros.h" #include "art/Framework/Principal/Event.h" #include "art/Framework/Principal/Handle.h" -#include "art/Framework/Principal/Run.h" -#include "art/Framework/Principal/SubRun.h" #include "canvas/Utilities/InputTag.h" #include "fhiclcpp/ParameterSet.h" #include "messagefacility/MessageLogger/MessageLogger.h" +#include "canvas/Persistency/Common/FindOneP.h" #include "canvas/Persistency/Common/FindManyP.h" +#include // std::find() +#include // std::abs() #include +#include // std::move() +#include -#include "lardataobj/RecoBase/Wire.h" #include "lardataobj/RecoBase/OpFlash.h" #include "canvas/Persistency/Common/Assns.h" #include "art/Persistency/Common/PtrMaker.h" @@ -30,8 +30,17 @@ #include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h" -class IcarusNuSliceHitsProducer; class IcarusNuSliceHitsProducer : public art::EDProducer { + +/** + * @brief Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction. + * + * Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction. + * It also produces a vector of ints, that maps the new hit collection to the original one. + * Optionally, it can store the association to the trigger flash used to pick the slice. + * + */ + public: explicit IcarusNuSliceHitsProducer(fhicl::ParameterSet const& p); // The compiler-generated destructor is fine for non-base @@ -47,13 +56,11 @@ class IcarusNuSliceHitsProducer : public art::EDProducer { private: // Declare member data here. - std::string fBaryMatchLabel; - std::string fSliceLabel; - std::string fSpsLabel; + art::InputTag const fBaryMatchLabel; + art::InputTag const fSliceLabel; + art::InputTag const fSpsLabel; bool doFlashAssn; - std::vector nearwires; - // Required functions. void produce(art::Event& e) override; @@ -62,9 +69,9 @@ class IcarusNuSliceHitsProducer : public art::EDProducer { IcarusNuSliceHitsProducer::IcarusNuSliceHitsProducer(fhicl::ParameterSet const& p) : EDProducer{p}, - fBaryMatchLabel( p.get("BaryMatchLabel","tpcpmtbarycentermatchCryoE")), - fSliceLabel(p.get("SliceLabel","pandoraGausCryoE")), - fSpsLabel(p.get("SpsLabel","cluster3DCryoE")), + fBaryMatchLabel( p.get("BaryMatchLabel","tpcpmtbarycentermatchCryoE")), + fSliceLabel(p.get("SliceLabel","pandoraGausCryoE")), + fSpsLabel(p.get("SpsLabel","cluster3DCryoE")), doFlashAssn(p.get("DoFlashAssn",false)) // More initializers here. { @@ -92,66 +99,69 @@ void IcarusNuSliceHitsProducer::produce(art::Event& e) auto outputSlcFlhAssns = std::make_unique>(); art::ValidHandle< std::vector< sbn::TPCPMTBarycenterMatch > > baryMatchListHandle = e.getValidHandle >(fBaryMatchLabel); - art::FindManyP msl(baryMatchListHandle, e, fBaryMatchLabel); - art::FindManyP mfl(baryMatchListHandle, e, fBaryMatchLabel); + art::FindOneP msl(baryMatchListHandle, e, fBaryMatchLabel); + art::FindOneP mfl(baryMatchListHandle, e, fBaryMatchLabel); float minRad = 99999.; - int slkey = -1; - int flkey = -1; - int bmkey = -1; + art::Ptr triggeringSlice; + art::Ptr triggeringFlash; for (size_t ibm=0; ibmsize(); ++ibm) { - art::Ptr bm(baryMatchListHandle,ibm); - if (msl.at(ibm).size()>0 && mfl.at(ibm).size()>0) { - std::cout << "ibm=" << ibm << " radius=" << bm->radius << " radius_Trigger=" << bm->radius_Trigger << " flashTime=" << bm->flashTime - << " slkey=" << msl.at(ibm).at(0).key() << " center=" << bm->chargeCenter - << " flkey=" << mfl.at(ibm).at(0).key() << " center=" << bm->flashCenter - << std::endl; - if (bm->radius_Trigger>0 && bm->radius_Triggerradius-bm->radius_Trigger)<0.00001) { - minRad = bm->radius_Trigger; - slkey = msl.at(ibm).at(0).key(); - flkey = mfl.at(ibm).at(0).key(); - bmkey = ibm; - } + sbn::TPCPMTBarycenterMatch const& bm = (*baryMatchListHandle)[ibm]; + // require the radius to be set (i.e. that there is a best match) and to match the one of the trigger within an arbitrarily small enough margin + if (bm.radius_Trigger<=0 || std::abs(bm.radius-bm.radius_Trigger)>=0.00001) continue; + art::Ptr const& slice = msl.at(ibm); + art::Ptr const& flash = mfl.at(ibm); + if (!slice || !flash) throw std::logic_error("Unexpected missing flash or slice in barycenter matching associations"); + mf::LogTrace{ "IcarusNuSliceHitsProducer" } + << "ibm=" << ibm << " radius=" << bm.radius << " radius_Trigger=" << bm.radius_Trigger << " flashTime=" << bm.flashTime + << " slkey=" << msl.at(ibm).key() << " center=" << bm.chargeCenter + << " flkey=" << mfl.at(ibm).key() << " center=" << bm.flashCenter; + if (bm.radius_Trigger +#include // std::find() +#include // std::move() +#include -#include "lardataobj/RecoBase/Wire.h" #include "lardataobj/RecoBase/OpFlash.h" #include "canvas/Persistency/Common/Assns.h" #include "art/Persistency/Common/PtrMaker.h" @@ -33,9 +30,19 @@ #include "larsim/MCCheater/ParticleInventoryService.h" #include "larsim/MCCheater/BackTrackerService.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" +#include "larcore/CoreUtils/ServiceUtil.h" -class IcarusTrueNuSliceHitsProducer; class IcarusTrueNuSliceHitsProducer : public art::EDProducer { + +/** + * @brief Produce separate collections of hits and space point from the slice with most true hits from the neutrino interaction. + * + * Produce separate collections of hits and space point from the slice with most true hits from the neutrino interaction. + * It also produces a vector of ints, that maps the new hit collection to the original one. + * Optionally, it can store the association to the trigger flash matched to the slice. + * + */ + public: explicit IcarusTrueNuSliceHitsProducer(fhicl::ParameterSet const& p); // The compiler-generated destructor is fine for non-base @@ -51,12 +58,10 @@ class IcarusTrueNuSliceHitsProducer : public art::EDProducer { private: // Declare member data here. - std::string fBaryMatchLabel; - std::string fSliceLabel; - std::string fSpsLabel; - bool doFlashAssn; - - std::vector nearwires; + art::InputTag const fBaryMatchLabel; + art::InputTag const fSliceLabel; + art::InputTag const fSpsLabel; + bool const doFlashAssn; // Required functions. void produce(art::Event& e) override; @@ -66,9 +71,9 @@ class IcarusTrueNuSliceHitsProducer : public art::EDProducer { IcarusTrueNuSliceHitsProducer::IcarusTrueNuSliceHitsProducer(fhicl::ParameterSet const& p) : EDProducer{p}, - fBaryMatchLabel( p.get("BaryMatchLabel","tpcpmtbarycentermatchCryoE")), - fSliceLabel(p.get("SliceLabel","pandoraGausCryoE")), - fSpsLabel(p.get("SpsLabel","cluster3DCryoE")), + fBaryMatchLabel( p.get("BaryMatchLabel") ), + fSliceLabel(p.get("SliceLabel")), + fSpsLabel(p.get("SpsLabel")), doFlashAssn(p.get("DoFlashAssn",false)) // More initializers here. { @@ -95,63 +100,65 @@ void IcarusTrueNuSliceHitsProducer::produce(art::Event& e) auto outputSpsHitAssns = std::make_unique>(); auto outputSlcFlhAssns = std::make_unique>(); - art::ServiceHandle pi; - art::ServiceHandle bt_h; + auto const* pi = lar::providerFrom(); + auto const* bt = lar::providerFrom(); auto const clockData = art::ServiceHandle()->DataFor(e); art::ValidHandle > inputSlice = e.getValidHandle >(fSliceLabel); - auto assocSliceHit = std::unique_ptr >(new art::FindManyP(inputSlice, e, fSliceLabel)); + art::FindManyP assocSliceHit(inputSlice, e, fSliceLabel); int slkey = -1; int maxCnt = 0; + // find the slice with the largest number of hits contributed by neutrino interactions for (size_t sk=0; sksize(); sk++) { - int slcnt = 0; - for (auto hit : assocSliceHit->at(sk)) { - std::vector h_ides = bt_h->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick()); - for (auto& tide : h_ides) { + int slhcnt = 0; + for (auto hit : assocSliceHit.at(sk)) { + std::vector const& h_ides = bt->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick()); + for (auto const& tide : h_ides) { int tid = tide.trackID; auto mct = pi->TrackIdToMCTruth_P(abs(tid)); - int fromNu = (mct.isAvailable() ? mct->NeutrinoSet() : 0); - if (fromNu==1) { - slcnt++; + bool fromNu = mct.isAvailable() && mct->NeutrinoSet(); + if (fromNu) { + slhcnt++; break; } } - if (slcnt>maxCnt) { - slkey = sk; - maxCnt = slcnt; - } + } + if (slhcnt>maxCnt) { + slkey = sk; + maxCnt = slhcnt; } } - std::cout << "keys: slice=" << slkey << std::endl; + mf::LogTrace{ "IcarusTrueNuSliceHitsProducer" } << "keys: slice=" << slkey; if (slkey>=0) { - std::cout << "slice hits=" << assocSliceHit->at(slkey).size() << std::endl; - for (auto h : assocSliceHit->at(slkey)) { + auto const& hits = assocSliceHit.at(slkey); + mf::LogTrace{ "IcarusTrueNuSliceHitsProducer" } << "slice hits=" << hits.size(); + for (auto const& h : hits) { assocSliceHitKeys->push_back(h.key()); - outputHits->emplace_back(*h); } + std::sort(begin(*assocSliceHitKeys), end(*assocSliceHitKeys)); + outputHits->reserve(assocSliceHitKeys->size()); + for (auto const& hit: *assocSliceHitKeys) outputHits->push_back(*hits[hit]); } // Get spacepoints from the event record - art::Handle> spListHandle; - std::vector> splist; - if (e.getByLabel(fSpsLabel, spListHandle)) art::fill_ptr_vector(splist, spListHandle); + auto splist = e.getValidHandle>(fSpsLabel); size_t countsps = 0; // Get assocations from spacepoints to hits - art::FindManyP fmp(spListHandle, e, fSpsLabel); - for (size_t spIdx = 0; spIdx < splist.size(); ++spIdx) { - // if (splist[spIdx]->Chisq()>0.5) continue; - auto assochits = fmp.at(spIdx); + art::FindManyP spacePointToHits(splist, e, fSpsLabel); + for (size_t spIdx = 0; spIdx < splist->size(); ++spIdx) { + auto const& assochits = spacePointToHits.at(spIdx); + // Consider only space points with hits associated on all planes. This is enough for NuGraph. if (assochits.size()<3) continue; // std::vector hitpos; for (size_t j = 0; j < assochits.size(); ++j) { - auto pos = std::find(assocSliceHitKeys->begin(),assocSliceHitKeys->end(),assochits[j].key()); - if (pos==assocSliceHitKeys->end()) break; + auto pos = std::lower_bound(assocSliceHitKeys->begin(),assocSliceHitKeys->end(),assochits[j].key()); + if ( (pos == assocSliceHitKeys->end()) || (*pos != assochits[j].key()) ) break; hitpos.push_back(pos-assocSliceHitKeys->begin()); } if (hitpos.size()<3) continue; - outputSpacepoints->emplace_back(*splist[spIdx]); + outputSpacepoints->emplace_back((*splist)[spIdx]); // const art::Ptr sps = spsPtrMaker(outputSpacepoints->size()-1); for (size_t j = 0; j < hitpos.size(); ++j) { @@ -160,7 +167,7 @@ void IcarusTrueNuSliceHitsProducer::produce(art::Event& e) } countsps++; } // for spacepoint - std::cout << "sps count=" << countsps << std::endl; + mf::LogTrace{ "IcarusTrueNuSliceHitsProducer" } << "sps count=" << countsps; art::ValidHandle< std::vector< sbn::TPCPMTBarycenterMatch > > baryMatchListHandle = e.getValidHandle >(fBaryMatchLabel); art::FindManyP msl(baryMatchListHandle, e, fBaryMatchLabel); @@ -171,10 +178,9 @@ void IcarusTrueNuSliceHitsProducer::produce(art::Event& e) art::Ptr bm(baryMatchListHandle,ibm); if (int(msl.at(ibm).at(0).key())!=slkey) continue; if (mfl.at(ibm).size()>0) { - std::cout << "ibm=" << ibm << " radius=" << bm->radius << " radius_Trigger=" << bm->radius_Trigger << " flashTime=" << bm->flashTime - << " slkey=" << msl.at(ibm).at(0).key() << " center=" << bm->chargeCenter - << " flkey=" << mfl.at(ibm).at(0).key() << " center=" << bm->flashCenter - << std::endl; + mf::LogTrace{ "IcarusTrueNuSliceHitsProducer" } << "ibm=" << ibm << " radius=" << bm->radius << " radius_Trigger=" << bm->radius_Trigger << " flashTime=" << bm->flashTime + << " slkey=" << msl.at(ibm).at(0).key() << " center=" << bm->chargeCenter + << " flkey=" << mfl.at(ibm).at(0).key() << " center=" << bm->flashCenter; if (bm->radius Date: Mon, 9 Jun 2025 13:57:52 -0500 Subject: [PATCH 15/24] moved NuGraph folder, improved fcl documentation --- icaruscode/CMakeLists.txt | 1 - icaruscode/TPC/CMakeLists.txt | 1 + .../NuGraph}/CMakeLists.txt | 0 .../NuGraph}/IcarusHDF5Maker_module.cc | 0 .../NuGraph}/IcarusNuGraphLoader_tool.cc | 0 .../IcarusNuSliceHitsProducer_module.cc | 0 .../IcarusTrueNuSliceHitsProducer_module.cc | 0 .../NuGraph}/StitchingUtils.h | 0 .../NuGraph}/cafmaker_nugraph_test.fcl | 0 .../NuGraph}/hdf5maker_icarus_slice.fcl | 7 ++++- .../NuGraph}/nugraph_icarus.fcl | 7 ----- .../NuGraph}/scripts/CMakeLists.txt | 0 .../scripts/setup_tritonserver-nugraph-v0.sh | 0 .../setup_tritonserver-nugraph-v0_grid.sh | 0 .../NuGraph}/testinference_slice_icarus.fcl | 26 ++++++++++++++----- 15 files changed, 27 insertions(+), 15 deletions(-) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/CMakeLists.txt (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/IcarusHDF5Maker_module.cc (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/IcarusNuGraphLoader_tool.cc (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/IcarusNuSliceHitsProducer_module.cc (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/IcarusTrueNuSliceHitsProducer_module.cc (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/StitchingUtils.h (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/cafmaker_nugraph_test.fcl (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/hdf5maker_icarus_slice.fcl (93%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/nugraph_icarus.fcl (79%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/scripts/CMakeLists.txt (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/scripts/setup_tritonserver-nugraph-v0.sh (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/scripts/setup_tritonserver-nugraph-v0_grid.sh (100%) rename icaruscode/{NuGraphIcarus => TPC/NuGraph}/testinference_slice_icarus.fcl (79%) diff --git a/icaruscode/CMakeLists.txt b/icaruscode/CMakeLists.txt index a2d709313..606f87f42 100644 --- a/icaruscode/CMakeLists.txt +++ b/icaruscode/CMakeLists.txt @@ -10,7 +10,6 @@ add_subdirectory(Filters) add_subdirectory(Generators) add_subdirectory(Geometry) add_subdirectory(LArG4) -add_subdirectory(NuGraphIcarus) add_subdirectory(PMT) add_subdirectory(RecoUtils) add_subdirectory(TPC) diff --git a/icaruscode/TPC/CMakeLists.txt b/icaruscode/TPC/CMakeLists.txt index 3bca0be5a..2652769c2 100644 --- a/icaruscode/TPC/CMakeLists.txt +++ b/icaruscode/TPC/CMakeLists.txt @@ -9,3 +9,4 @@ add_subdirectory(Simulation) add_subdirectory(Tracking) add_subdirectory(ICARUSWireCell) add_subdirectory(SPAna) +add_subdirectory(NuGraph) diff --git a/icaruscode/NuGraphIcarus/CMakeLists.txt b/icaruscode/TPC/NuGraph/CMakeLists.txt similarity index 100% rename from icaruscode/NuGraphIcarus/CMakeLists.txt rename to icaruscode/TPC/NuGraph/CMakeLists.txt diff --git a/icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc b/icaruscode/TPC/NuGraph/IcarusHDF5Maker_module.cc similarity index 100% rename from icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc rename to icaruscode/TPC/NuGraph/IcarusHDF5Maker_module.cc diff --git a/icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc b/icaruscode/TPC/NuGraph/IcarusNuGraphLoader_tool.cc similarity index 100% rename from icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc rename to icaruscode/TPC/NuGraph/IcarusNuGraphLoader_tool.cc diff --git a/icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/IcarusNuSliceHitsProducer_module.cc similarity index 100% rename from icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc rename to icaruscode/TPC/NuGraph/IcarusNuSliceHitsProducer_module.cc diff --git a/icaruscode/NuGraphIcarus/IcarusTrueNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/IcarusTrueNuSliceHitsProducer_module.cc similarity index 100% rename from icaruscode/NuGraphIcarus/IcarusTrueNuSliceHitsProducer_module.cc rename to icaruscode/TPC/NuGraph/IcarusTrueNuSliceHitsProducer_module.cc diff --git a/icaruscode/NuGraphIcarus/StitchingUtils.h b/icaruscode/TPC/NuGraph/StitchingUtils.h similarity index 100% rename from icaruscode/NuGraphIcarus/StitchingUtils.h rename to icaruscode/TPC/NuGraph/StitchingUtils.h diff --git a/icaruscode/NuGraphIcarus/cafmaker_nugraph_test.fcl b/icaruscode/TPC/NuGraph/cafmaker_nugraph_test.fcl similarity index 100% rename from icaruscode/NuGraphIcarus/cafmaker_nugraph_test.fcl rename to icaruscode/TPC/NuGraph/cafmaker_nugraph_test.fcl diff --git a/icaruscode/NuGraphIcarus/hdf5maker_icarus_slice.fcl b/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl similarity index 93% rename from icaruscode/NuGraphIcarus/hdf5maker_icarus_slice.fcl rename to icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl index 6cd27b3af..71dae04a2 100644 --- a/icaruscode/NuGraphIcarus/hdf5maker_icarus_slice.fcl +++ b/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl @@ -1,7 +1,12 @@ +### +### fcl file for producing the hdf5 files used as input for NuGraph training. +### Author: G. Cerati (FNAL) +### + #include "services_common_icarus.fcl" #include "services_icarus_simulation.fcl" -process_name: trueenergy +process_name: makehdf5 services: { diff --git a/icaruscode/NuGraphIcarus/nugraph_icarus.fcl b/icaruscode/TPC/NuGraph/nugraph_icarus.fcl similarity index 79% rename from icaruscode/NuGraphIcarus/nugraph_icarus.fcl rename to icaruscode/TPC/NuGraph/nugraph_icarus.fcl index 24d2b39c1..071aa2801 100644 --- a/icaruscode/NuGraphIcarus/nugraph_icarus.fcl +++ b/icaruscode/TPC/NuGraph/nugraph_icarus.fcl @@ -16,13 +16,6 @@ nuslhitsCryoW: { SpsLabel: "cluster3DCryoW" } -#NuGraphCryoE: @local::ApptainerNuGraphNuSonicTriton -#NuGraphCryoE.LoaderTool.tool_type: "IcarusNuGraphLoader" -#NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE" -#NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE" -#NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE" -#NuGraphCryoE.TritonConfig.modelName: "nugraph2_icarus_mpvmprbnb" - NuGraphCryoE: @local::NuGraphLibTorch NuGraphCryoE.LoaderTool.tool_type: "IcarusNuGraphLoader" NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE" diff --git a/icaruscode/NuGraphIcarus/scripts/CMakeLists.txt b/icaruscode/TPC/NuGraph/scripts/CMakeLists.txt similarity index 100% rename from icaruscode/NuGraphIcarus/scripts/CMakeLists.txt rename to icaruscode/TPC/NuGraph/scripts/CMakeLists.txt diff --git a/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh b/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0.sh similarity index 100% rename from icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0.sh rename to icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0.sh diff --git a/icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh b/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh similarity index 100% rename from icaruscode/NuGraphIcarus/scripts/setup_tritonserver-nugraph-v0_grid.sh rename to icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh diff --git a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl b/icaruscode/TPC/NuGraph/testinference_slice_icarus.fcl similarity index 79% rename from icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl rename to icaruscode/TPC/NuGraph/testinference_slice_icarus.fcl index b88e8e8a6..9b64ee70f 100644 --- a/icaruscode/NuGraphIcarus/testinference_slice_icarus.fcl +++ b/icaruscode/TPC/NuGraph/testinference_slice_icarus.fcl @@ -1,13 +1,12 @@ +### +### Test fcl file for NuGraph inference. For each cryostat it runs two producers (IcarusNuSliceHitsProducer and NuGraphInference) plus one analyzer (NuGraphAnalyzer). +### Author: G. Cerati (FNAL) +### + #include "nugraph_icarus.fcl" #include "services_common_icarus.fcl" #include "services_icarus_simulation.fcl" -# in order to setup the triton server locally, please do: -# > mrbslp -# > setup_tritonservier-nugraph-v0.sh -# then after the job you can kill the server process with: -# > killall -9 /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/x86_64/libexec/apptainer/libexec/starter - process_name: testinference services: @@ -98,3 +97,18 @@ physics: services.BackTrackerService.BackTracker.G4ModuleLabel: "largeant" services.BackTrackerService.BackTracker.SimChannelModuleLabel: "daq:simpleSC" + + +### Parameters and instructions to test the NuSonic inference model (still under development) +# in order to setup the triton server locally, please do: +# > mrbslp +# > setup_tritonservier-nugraph-v0.sh +# then after the job you can kill the server process with: +# > killall -9 /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/x86_64/libexec/apptainer/libexec/starter +#NuGraphCryoE: @local::ApptainerNuGraphNuSonicTriton +#NuGraphCryoE.LoaderTool.tool_type: "IcarusNuGraphLoader" +#NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE" +#NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE" +#NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE" +#NuGraphCryoE.TritonConfig.modelName: "nugraph2_icarus_mpvmprbnb" + From ac05da7b2c1bad14940d25f323b4a7b554c4e7f3 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Mon, 9 Jun 2025 13:58:19 -0500 Subject: [PATCH 16/24] remove cafmaker_nugraph_test.fcl --- icaruscode/TPC/NuGraph/cafmaker_nugraph_test.fcl | 6 ------ 1 file changed, 6 deletions(-) delete mode 100644 icaruscode/TPC/NuGraph/cafmaker_nugraph_test.fcl diff --git a/icaruscode/TPC/NuGraph/cafmaker_nugraph_test.fcl b/icaruscode/TPC/NuGraph/cafmaker_nugraph_test.fcl deleted file mode 100644 index e217b591e..000000000 --- a/icaruscode/TPC/NuGraph/cafmaker_nugraph_test.fcl +++ /dev/null @@ -1,6 +0,0 @@ -#include "cafmakerjob_icarus_detsim2d.fcl" -#### _systtools_and_fluxwgt - -physics.producers.cafmaker.NuGraphSliceHitLabel: "nuslhits" -physics.producers.cafmaker.NuGraphFilterLabel: "NuGraph:filter" -physics.producers.cafmaker.NuGraphSemanticLabel: "NuGraph:semantic" From 5297b0c9059ab9fad3a0c536ad3bbb41cf34a19c Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Mon, 9 Jun 2025 14:35:57 -0500 Subject: [PATCH 17/24] rename Icarus to ICARUS --- icaruscode/TPC/NuGraph/CMakeLists.txt | 8 +-- ...er_module.cc => ICARUSHDF5Maker_module.cc} | 50 +++++++++---------- ...er_tool.cc => ICARUSNuGraphLoader_tool.cc} | 16 +++--- ...cc => ICARUSNuSliceHitsProducer_module.cc} | 30 +++++------ ...> ICARUSTrueNuSliceHitsProducer_module.cc} | 32 ++++++------ .../TPC/NuGraph/hdf5maker_icarus_slice.fcl | 10 ++-- icaruscode/TPC/NuGraph/nugraph_icarus.fcl | 6 +-- .../NuGraph/testinference_slice_icarus.fcl | 8 +-- 8 files changed, 80 insertions(+), 80 deletions(-) rename icaruscode/TPC/NuGraph/{IcarusHDF5Maker_module.cc => ICARUSHDF5Maker_module.cc} (95%) rename icaruscode/TPC/NuGraph/{IcarusNuGraphLoader_tool.cc => ICARUSNuGraphLoader_tool.cc} (91%) rename icaruscode/TPC/NuGraph/{IcarusNuSliceHitsProducer_module.cc => ICARUSNuSliceHitsProducer_module.cc} (88%) rename icaruscode/TPC/NuGraph/{IcarusTrueNuSliceHitsProducer_module.cc => ICARUSTrueNuSliceHitsProducer_module.cc} (87%) diff --git a/icaruscode/TPC/NuGraph/CMakeLists.txt b/icaruscode/TPC/NuGraph/CMakeLists.txt index 0392247c3..7fdd40479 100644 --- a/icaruscode/TPC/NuGraph/CMakeLists.txt +++ b/icaruscode/TPC/NuGraph/CMakeLists.txt @@ -11,19 +11,19 @@ set( nugraph_tool_lib_list include_directories($ENV{HEP_HPC_INC}) -cet_build_plugin(IcarusNuGraphLoader art::tool +cet_build_plugin(ICARUSNuGraphLoader art::tool LIBRARIES PRIVATE ${nugraph_tool_lib_list} ) -simple_plugin(IcarusNuSliceHitsProducer "module" +simple_plugin(ICARUSNuSliceHitsProducer "module" LIBRARIES PRIVATE art::Framework_Core lardataobj::RecoBase sbnobj::Common_Reco ) -simple_plugin(IcarusTrueNuSliceHitsProducer "module" +simple_plugin(ICARUSTrueNuSliceHitsProducer "module" LIBRARIES PRIVATE art::Framework_Core nusimdata::SimulationBase @@ -35,7 +35,7 @@ simple_plugin(IcarusTrueNuSliceHitsProducer "module" lardata::DetectorClocksService ) -simple_plugin(IcarusHDF5Maker "module" +simple_plugin(ICARUSHDF5Maker "module" LIBRARIES PRIVATE art::Framework_Core nusimdata::SimulationBase diff --git a/icaruscode/TPC/NuGraph/IcarusHDF5Maker_module.cc b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc similarity index 95% rename from icaruscode/TPC/NuGraph/IcarusHDF5Maker_module.cc rename to icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc index a6ff7e9ee..5025ba1fb 100644 --- a/icaruscode/TPC/NuGraph/IcarusHDF5Maker_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc @@ -1,6 +1,6 @@ /** - * @file icaruscode/NuGraphIcarus/IcarusHDF5Maker_module.cc - * @brief Implementation of `IcarusHDF5Maker` _art_ module. + * @file icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc + * @brief Implementation of `ICARUSHDF5Maker` _art_ module. * @author Giuseppe Cerati (cerati@fnal.gov), V Hewes */ @@ -44,15 +44,15 @@ #include "StitchingUtils.h" -class IcarusHDF5Maker : public art::EDAnalyzer { +class ICARUSHDF5Maker : public art::EDAnalyzer { public: - explicit IcarusHDF5Maker(fhicl::ParameterSet const& p); - ~IcarusHDF5Maker() noexcept {}; // bare pointers are cleaned up by endSubRun + explicit ICARUSHDF5Maker(fhicl::ParameterSet const& p); + ~ICARUSHDF5Maker() noexcept {}; // bare pointers are cleaned up by endSubRun - IcarusHDF5Maker(IcarusHDF5Maker const&) = delete; - IcarusHDF5Maker(IcarusHDF5Maker&&) = delete; - IcarusHDF5Maker& operator=(IcarusHDF5Maker const&) = delete; - IcarusHDF5Maker& operator=(IcarusHDF5Maker&&) = delete; + ICARUSHDF5Maker(ICARUSHDF5Maker const&) = delete; + ICARUSHDF5Maker(ICARUSHDF5Maker&&) = delete; + ICARUSHDF5Maker& operator=(ICARUSHDF5Maker const&) = delete; + ICARUSHDF5Maker& operator=(ICARUSHDF5Maker&&) = delete; void beginSubRun(art::SubRun const& sr) override; void endSubRun(art::SubRun const& sr) override; @@ -283,7 +283,7 @@ class IcarusHDF5Maker : public art::EDAnalyzer { int NearWire(const geo::PlaneGeo &plane, float x, float y, float z) const; }; -IcarusHDF5Maker::IcarusHDF5Maker(fhicl::ParameterSet const& p) +ICARUSHDF5Maker::ICARUSHDF5Maker(fhicl::ParameterSet const& p) : EDAnalyzer{p}, fTruthLabel(p.get("TruthLabel")), fHitLabel( p.get("HitLabel")), @@ -300,7 +300,7 @@ IcarusHDF5Maker::IcarusHDF5Maker(fhicl::ParameterSet const& p) << "EventInfo must be \"none\" or \"nu\", not " << fEventInfo; } -void IcarusHDF5Maker::analyze(art::Event const& e) { +void ICARUSHDF5Maker::analyze(art::Event const& e) { const cheat::BackTrackerService* bt = fUseMap? nullptr: art::ServiceHandle().get(); geo::WireReadoutGeom const& geom = art::ServiceHandle()->Get(); @@ -329,7 +329,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // Fill event table if (fEventInfo == "none") { fHDFData->eventNtuple->insert( evtID.data() ); - mf::LogInfo("IcarusHDF5Maker") << "Filling event table" + mf::LogInfo("ICARUSHDF5Maker") << "Filling event table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2]; } @@ -361,7 +361,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // << std::endl; // } - mf::LogDebug("IcarusHDF5Maker") << "Filling event table" + mf::LogDebug("ICARUSHDF5Maker") << "Filling event table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2] << "\nis cc? " << (nutruth.CCNC() == simb::kCC) @@ -405,7 +405,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { } fHDFData->spacePointNtuple.insert(evtID.data(),splist[i].ID(), pos.data(), hitID.data(),splist[i].Chisq() ); - mf::LogDebug("IcarusHDF5Maker") << "Filling spacepoint table" + mf::LogDebug("ICARUSHDF5Maker") << "Filling spacepoint table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2] << "\nspacepoint id " << splist[i].ID() @@ -435,7 +435,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { wireid.Plane, wireid.Wire, hit->PeakTime(),wireid.Cryostat ); - mf::LogInfo("IcarusHDF5Maker") << "Filling hit table" + mf::LogInfo("ICARUSHDF5Maker") << "Filling hit table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2] << "\nhit id " << hit.key() << ", integral " @@ -464,7 +464,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()); - mf::LogInfo("IcarusHDF5Maker") << "Filling energy deposit table" + mf::LogInfo("ICARUSHDF5Maker") << "Filling energy deposit table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2] << "\nhit id " << hit.key() << ", g4 id " @@ -501,7 +501,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { // Loop over true particles and fill table for (auto [ id, p ] : allIDs) { if (!p) { - mf::LogError("IcarusHDF5Maker") << "Track ID=" << id << " not tracked back to any particle."; + mf::LogError("ICARUSHDF5Maker") << "Track ID=" << id << " not tracked back to any particle."; continue; } auto mct = pi->TrackIdToMCTruth_P(abs(id)); @@ -518,7 +518,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { particleStart.data(), particleEnd.data(), p->Process(), p->EndProcess() ); - mf::LogDebug("IcarusHDF5Maker") << "Filling particle table" + mf::LogDebug("ICARUSHDF5Maker") << "Filling particle table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2] << "\ng4 id " << abs(id) << ", pdg code " @@ -575,7 +575,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { count++; } flashsumpepmtmap.push_back(sumpepmtmap); - mf::LogDebug("IcarusHDF5Maker") << "Filling opflash table" + mf::LogDebug("ICARUSHDF5Maker") << "Filling opflash table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2] << "\nflash id " << flkey << ", Time " << opflash.Time() @@ -614,7 +614,7 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { ophit->Area(),ophit->Amplitude(),ophit->PE(), (isInFlash ? flashsumpepmtmap[0][ophit->OpChannel()] : -1) ); - mf::LogDebug("IcarusHDF5Maker") << "\nFilling ophit table" + mf::LogDebug("ICARUSHDF5Maker") << "\nFilling ophit table" << "\nrun " << evtID[0] << ", subrun " << evtID[1] << ", event " << evtID[2] << "\nhit id " << ophit.key() << ", channel " @@ -625,17 +625,17 @@ void IcarusHDF5Maker::analyze(art::Event const& e) { } //End optical analyzer -} // function IcarusHDF5Maker::analyze +} // function ICARUSHDF5Maker::analyze -void IcarusHDF5Maker::beginSubRun(art::SubRun const& sr) { +void ICARUSHDF5Maker::beginSubRun(art::SubRun const& sr) { fHDFData = std::make_unique(fOutputName, fEventInfo, sr.id()); } -void IcarusHDF5Maker::endSubRun(art::SubRun const& sr) { +void ICARUSHDF5Maker::endSubRun(art::SubRun const& sr) { fHDFData.reset(); } -int IcarusHDF5Maker::NearWire(const geo::PlaneGeo &plane, float x, float y, float z) const +int ICARUSHDF5Maker::NearWire(const geo::PlaneGeo &plane, float x, float y, float z) const { geo::WireID wireID; try { @@ -648,4 +648,4 @@ int IcarusHDF5Maker::NearWire(const geo::PlaneGeo &plane, float x, float y, floa return wireID.Wire; } -DEFINE_ART_MODULE(IcarusHDF5Maker) +DEFINE_ART_MODULE(ICARUSHDF5Maker) diff --git a/icaruscode/TPC/NuGraph/IcarusNuGraphLoader_tool.cc b/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc similarity index 91% rename from icaruscode/TPC/NuGraph/IcarusNuGraphLoader_tool.cc rename to icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc index e3fea33f6..a09b3b253 100644 --- a/icaruscode/TPC/NuGraph/IcarusNuGraphLoader_tool.cc +++ b/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc @@ -1,5 +1,5 @@ /** - * @file icaruscode/NuGraphIcarus/IcarusNuGraphLoader_tool.cc + * @file icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc * @author Giuseppe Cerati (cerati@fnal.gov) */ @@ -19,7 +19,7 @@ #include "StitchingUtils.h" -class IcarusNuGraphLoader : public LoaderToolBase { +class ICARUSNuGraphLoader : public LoaderToolBase { public: /** @@ -27,7 +27,7 @@ class IcarusNuGraphLoader : public LoaderToolBase { * * @param pset */ - IcarusNuGraphLoader(const fhicl::ParameterSet& pset); + ICARUSNuGraphLoader(const fhicl::ParameterSet& pset); /** * @brief loadData function @@ -45,11 +45,11 @@ class IcarusNuGraphLoader : public LoaderToolBase { art::InputTag const fSpsInput; }; -IcarusNuGraphLoader::IcarusNuGraphLoader(const fhicl::ParameterSet& p) +ICARUSNuGraphLoader::ICARUSNuGraphLoader(const fhicl::ParameterSet& p) : fHitInput{p.get("hitInput")}, fSpsInput{p.get("spsInput")} {} -void IcarusNuGraphLoader::loadData(art::Event& e, +void ICARUSNuGraphLoader::loadData(art::Event& e, vector>& hitlist, vector& inputs, vector>& idsmap) @@ -91,7 +91,7 @@ void IcarusNuGraphLoader::loadData(art::Event& e, hit_table_integral_data.push_back(h->Integral()); hit_table_rms_data.push_back(h->RMS()); } - mf::LogDebug{ "IcarusNuGraphLoader" } << "loader has nhits=" << hit_table_hit_id_data.size(); + mf::LogDebug{ "ICARUSNuGraphLoader" } << "loader has nhits=" << hit_table_hit_id_data.size(); // Get spacepoints from the event record art::Handle> spListHandle; @@ -126,7 +126,7 @@ void IcarusNuGraphLoader::loadData(art::Event& e, if (plane == 2) spacepoint_table_hit_id_y_data.back() = sp2Hit[i][j].key(); } } - mf::LogDebug{ "IcarusNuGraphLoader" } << "loader has nsps=" << spacepoint_table_hit_id_u_data.size(); + mf::LogDebug{ "ICARUSNuGraphLoader" } << "loader has nsps=" << spacepoint_table_hit_id_u_data.size(); inputs.emplace_back("hit_table_hit_id", std::move(hit_table_hit_id_data)); inputs.emplace_back("hit_table_local_plane", std::move(hit_table_local_plane_data)); @@ -140,4 +140,4 @@ void IcarusNuGraphLoader::loadData(art::Event& e, inputs.emplace_back("spacepoint_table_hit_id_v", std::move(spacepoint_table_hit_id_v_data)); inputs.emplace_back("spacepoint_table_hit_id_y", std::move(spacepoint_table_hit_id_y_data)); } -DEFINE_ART_CLASS_TOOL(IcarusNuGraphLoader) +DEFINE_ART_CLASS_TOOL(ICARUSNuGraphLoader) diff --git a/icaruscode/TPC/NuGraph/IcarusNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc similarity index 88% rename from icaruscode/TPC/NuGraph/IcarusNuSliceHitsProducer_module.cc rename to icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc index 7c521d5e6..73c5a6a96 100644 --- a/icaruscode/TPC/NuGraph/IcarusNuSliceHitsProducer_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc @@ -1,6 +1,6 @@ /** - * @file icaruscode/NuGraphIcarus/IcarusNuSliceHitsProducer_module.cc - * @brief Implementation of `IcarusNuSliceHitsProducer` _art_ module. + * @file icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc + * @brief Implementation of `ICARUSNuSliceHitsProducer` _art_ module. * @author Giuseppe Cerati (cerati@fnal.gov) * @date May 25, 2021 */ @@ -30,7 +30,7 @@ #include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h" -class IcarusNuSliceHitsProducer : public art::EDProducer { +class ICARUSNuSliceHitsProducer : public art::EDProducer { /** * @brief Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction. @@ -42,15 +42,15 @@ class IcarusNuSliceHitsProducer : public art::EDProducer { */ public: - explicit IcarusNuSliceHitsProducer(fhicl::ParameterSet const& p); + explicit ICARUSNuSliceHitsProducer(fhicl::ParameterSet const& p); // The compiler-generated destructor is fine for non-base // classes without bare pointers or other resource use. // Plugins should not be copied or assigned. - IcarusNuSliceHitsProducer(IcarusNuSliceHitsProducer const&) = delete; - IcarusNuSliceHitsProducer(IcarusNuSliceHitsProducer&&) = delete; - IcarusNuSliceHitsProducer& operator=(IcarusNuSliceHitsProducer const&) = delete; - IcarusNuSliceHitsProducer& operator=(IcarusNuSliceHitsProducer&&) = delete; + ICARUSNuSliceHitsProducer(ICARUSNuSliceHitsProducer const&) = delete; + ICARUSNuSliceHitsProducer(ICARUSNuSliceHitsProducer&&) = delete; + ICARUSNuSliceHitsProducer& operator=(ICARUSNuSliceHitsProducer const&) = delete; + ICARUSNuSliceHitsProducer& operator=(ICARUSNuSliceHitsProducer&&) = delete; private: @@ -67,7 +67,7 @@ class IcarusNuSliceHitsProducer : public art::EDProducer { }; -IcarusNuSliceHitsProducer::IcarusNuSliceHitsProducer(fhicl::ParameterSet const& p) +ICARUSNuSliceHitsProducer::ICARUSNuSliceHitsProducer(fhicl::ParameterSet const& p) : EDProducer{p}, fBaryMatchLabel( p.get("BaryMatchLabel","tpcpmtbarycentermatchCryoE")), fSliceLabel(p.get("SliceLabel","pandoraGausCryoE")), @@ -85,7 +85,7 @@ IcarusNuSliceHitsProducer::IcarusNuSliceHitsProducer(fhicl::ParameterSet const& // Call appropriate consumes<>() for any products to be retrieved by this module. } -void IcarusNuSliceHitsProducer::produce(art::Event& e) +void ICARUSNuSliceHitsProducer::produce(art::Event& e) { auto outputHits = std::make_unique >(); @@ -111,7 +111,7 @@ void IcarusNuSliceHitsProducer::produce(art::Event& e) art::Ptr const& slice = msl.at(ibm); art::Ptr const& flash = mfl.at(ibm); if (!slice || !flash) throw std::logic_error("Unexpected missing flash or slice in barycenter matching associations"); - mf::LogTrace{ "IcarusNuSliceHitsProducer" } + mf::LogTrace{ "ICARUSNuSliceHitsProducer" } << "ibm=" << ibm << " radius=" << bm.radius << " radius_Trigger=" << bm.radius_Trigger << " flashTime=" << bm.flashTime << " slkey=" << msl.at(ibm).key() << " center=" << bm.chargeCenter << " flkey=" << mfl.at(ibm).key() << " center=" << bm.flashCenter; @@ -121,7 +121,7 @@ void IcarusNuSliceHitsProducer::produce(art::Event& e) triggeringFlash = flash; } } - mf::LogTrace{ "IcarusNuSliceHitsProducer" } << "keys: slice=" << triggeringSlice.key() << " flash=" << triggeringFlash.key(); + mf::LogTrace{ "ICARUSNuSliceHitsProducer" } << "keys: slice=" << triggeringSlice.key() << " flash=" << triggeringFlash.key(); if (doFlashAssn) { if (triggeringSlice && triggeringFlash) { @@ -133,7 +133,7 @@ void IcarusNuSliceHitsProducer::produce(art::Event& e) art::FindManyP assocSliceHit(inputSlice, e, fSliceLabel); if (triggeringSlice) { - mf::LogTrace{ "IcarusNuSliceHitsProducer" } << "slice hits=" << assocSliceHit.at(triggeringSlice.key()).size(); + mf::LogTrace{ "ICARUSNuSliceHitsProducer" } << "slice hits=" << assocSliceHit.at(triggeringSlice.key()).size(); auto const& hits = assocSliceHit.at(triggeringSlice.key()); for (auto h : hits) { assocSliceHitKeys->push_back(h.key()); @@ -170,7 +170,7 @@ void IcarusNuSliceHitsProducer::produce(art::Event& e) } countsps++; } // for spacepoint - mf::LogTrace{ "IcarusNuSliceHitsProducer" } << "sps count=" << countsps; + mf::LogTrace{ "ICARUSNuSliceHitsProducer" } << "sps count=" << countsps; e.put(std::move(outputHits)); e.put(std::move(assocSliceHitKeys)); @@ -180,4 +180,4 @@ void IcarusNuSliceHitsProducer::produce(art::Event& e) } -DEFINE_ART_MODULE(IcarusNuSliceHitsProducer) +DEFINE_ART_MODULE(ICARUSNuSliceHitsProducer) diff --git a/icaruscode/TPC/NuGraph/IcarusTrueNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc similarity index 87% rename from icaruscode/TPC/NuGraph/IcarusTrueNuSliceHitsProducer_module.cc rename to icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc index ac6fc8d98..6a22493f0 100644 --- a/icaruscode/TPC/NuGraph/IcarusTrueNuSliceHitsProducer_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc @@ -1,6 +1,6 @@ /** - * @file icaruscode/NuGraphIcarus/IcarusTrueNuSliceHitsProducer_module.cc - * @brief Implementation of `IcarusTrueNuSliceHitsProducer` module + * @file icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc + * @brief Implementation of `ICARUSTrueNuSliceHitsProducer` module * @author Giuseppe Cerati (cerati@fnal.gov) */ @@ -32,7 +32,7 @@ #include "lardata/DetectorInfoServices/DetectorClocksService.h" #include "larcore/CoreUtils/ServiceUtil.h" -class IcarusTrueNuSliceHitsProducer : public art::EDProducer { +class ICARUSTrueNuSliceHitsProducer : public art::EDProducer { /** * @brief Produce separate collections of hits and space point from the slice with most true hits from the neutrino interaction. @@ -44,15 +44,15 @@ class IcarusTrueNuSliceHitsProducer : public art::EDProducer { */ public: - explicit IcarusTrueNuSliceHitsProducer(fhicl::ParameterSet const& p); + explicit ICARUSTrueNuSliceHitsProducer(fhicl::ParameterSet const& p); // The compiler-generated destructor is fine for non-base // classes without bare pointers or other resource use. // Plugins should not be copied or assigned. - IcarusTrueNuSliceHitsProducer(IcarusTrueNuSliceHitsProducer const&) = delete; - IcarusTrueNuSliceHitsProducer(IcarusTrueNuSliceHitsProducer&&) = delete; - IcarusTrueNuSliceHitsProducer& operator=(IcarusTrueNuSliceHitsProducer const&) = delete; - IcarusTrueNuSliceHitsProducer& operator=(IcarusTrueNuSliceHitsProducer&&) = delete; + ICARUSTrueNuSliceHitsProducer(ICARUSTrueNuSliceHitsProducer const&) = delete; + ICARUSTrueNuSliceHitsProducer(ICARUSTrueNuSliceHitsProducer&&) = delete; + ICARUSTrueNuSliceHitsProducer& operator=(ICARUSTrueNuSliceHitsProducer const&) = delete; + ICARUSTrueNuSliceHitsProducer& operator=(ICARUSTrueNuSliceHitsProducer&&) = delete; private: @@ -69,7 +69,7 @@ class IcarusTrueNuSliceHitsProducer : public art::EDProducer { }; -IcarusTrueNuSliceHitsProducer::IcarusTrueNuSliceHitsProducer(fhicl::ParameterSet const& p) +ICARUSTrueNuSliceHitsProducer::ICARUSTrueNuSliceHitsProducer(fhicl::ParameterSet const& p) : EDProducer{p}, fBaryMatchLabel( p.get("BaryMatchLabel") ), fSliceLabel(p.get("SliceLabel")), @@ -87,7 +87,7 @@ IcarusTrueNuSliceHitsProducer::IcarusTrueNuSliceHitsProducer(fhicl::ParameterSet // Call appropriate consumes<>() for any products to be retrieved by this module. } -void IcarusTrueNuSliceHitsProducer::produce(art::Event& e) +void ICARUSTrueNuSliceHitsProducer::produce(art::Event& e) { auto outputHits = std::make_unique >(); @@ -128,11 +128,11 @@ void IcarusTrueNuSliceHitsProducer::produce(art::Event& e) maxCnt = slhcnt; } } - mf::LogTrace{ "IcarusTrueNuSliceHitsProducer" } << "keys: slice=" << slkey; + mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "keys: slice=" << slkey; if (slkey>=0) { auto const& hits = assocSliceHit.at(slkey); - mf::LogTrace{ "IcarusTrueNuSliceHitsProducer" } << "slice hits=" << hits.size(); + mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "slice hits=" << hits.size(); for (auto const& h : hits) { assocSliceHitKeys->push_back(h.key()); } @@ -167,7 +167,7 @@ void IcarusTrueNuSliceHitsProducer::produce(art::Event& e) } countsps++; } // for spacepoint - mf::LogTrace{ "IcarusTrueNuSliceHitsProducer" } << "sps count=" << countsps; + mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "sps count=" << countsps; art::ValidHandle< std::vector< sbn::TPCPMTBarycenterMatch > > baryMatchListHandle = e.getValidHandle >(fBaryMatchLabel); art::FindManyP msl(baryMatchListHandle, e, fBaryMatchLabel); @@ -178,7 +178,7 @@ void IcarusTrueNuSliceHitsProducer::produce(art::Event& e) art::Ptr bm(baryMatchListHandle,ibm); if (int(msl.at(ibm).at(0).key())!=slkey) continue; if (mfl.at(ibm).size()>0) { - mf::LogTrace{ "IcarusTrueNuSliceHitsProducer" } << "ibm=" << ibm << " radius=" << bm->radius << " radius_Trigger=" << bm->radius_Trigger << " flashTime=" << bm->flashTime + mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "ibm=" << ibm << " radius=" << bm->radius << " radius_Trigger=" << bm->radius_Trigger << " flashTime=" << bm->flashTime << " slkey=" << msl.at(ibm).at(0).key() << " center=" << bm->chargeCenter << " flkey=" << mfl.at(ibm).at(0).key() << " center=" << bm->flashCenter; if (bm->radius killall -9 /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/x86_64/libexec/apptainer/libexec/starter #NuGraphCryoE: @local::ApptainerNuGraphNuSonicTriton -#NuGraphCryoE.LoaderTool.tool_type: "IcarusNuGraphLoader" +#NuGraphCryoE.LoaderTool.tool_type: "ICARUSNuGraphLoader" #NuGraphCryoE.LoaderTool.hitInput: "nuslhitsCryoE" #NuGraphCryoE.LoaderTool.spsInput: "nuslhitsCryoE" #NuGraphCryoE.DecoderTools.SemanticDecoderTool.hitInput: "nuslhitsCryoE" From a6a3c8fc77fc6254474f03414344ac62e3cb0bbb Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Mon, 9 Jun 2025 17:22:46 -0500 Subject: [PATCH 18/24] remove fEventInfo --- .../TPC/NuGraph/ICARUSHDF5Maker_module.cc | 244 ++++++++---------- .../TPC/NuGraph/hdf5maker_icarus_slice.fcl | 1 - 2 files changed, 107 insertions(+), 138 deletions(-) diff --git a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc index 5025ba1fb..a77e966c8 100644 --- a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc @@ -68,7 +68,6 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { art::InputTag fOpFlashLabel; bool fUseMap; - std::string fEventInfo; std::string fOutputName; std::vector tpc_ids_checked = {}; // add the TPC ids here so we only check them once @@ -76,16 +75,13 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { struct HDFDataFile { hep_hpc::hdf5::File file; ///< output HDF5 file - std::optional // event id (run, subrun, event) - >> eventNtuple; ///< event ntuple - - std::optional, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // is cc - hep_hpc::hdf5::Column, // nu pdg - hep_hpc::hdf5::Column, // nu energy - hep_hpc::hdf5::Column, // lep energy - hep_hpc::hdf5::Column // nu dir (x, y, z) - >> eventNtupleNu; ///< event ntuple with neutrino information + hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) + hep_hpc::hdf5::Column, // is cc + hep_hpc::hdf5::Column, // nu pdg + hep_hpc::hdf5::Column, // nu energy + hep_hpc::hdf5::Column, // lep energy + hep_hpc::hdf5::Column // nu dir (x, y, z) + > eventNtupleNu; ///< event ntuple with neutrino information hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) hep_hpc::hdf5::Column, // spacepoint id @@ -95,69 +91,69 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { > spacePointNtuple; ///< spacepoint ntuple hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // hit id - hep_hpc::hdf5::Column, // hit integral - hep_hpc::hdf5::Column, // hit rms - hep_hpc::hdf5::Column, // tpc id - hep_hpc::hdf5::Column, // global plane - hep_hpc::hdf5::Column, // global wire - hep_hpc::hdf5::Column, // global time - hep_hpc::hdf5::Column, // raw plane - hep_hpc::hdf5::Column, // raw wire - hep_hpc::hdf5::Column, // raw time - hep_hpc::hdf5::Column //cryostat + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // hit integral + hep_hpc::hdf5::Column, // hit rms + hep_hpc::hdf5::Column, // tpc id + hep_hpc::hdf5::Column, // global plane + hep_hpc::hdf5::Column, // global wire + hep_hpc::hdf5::Column, // global time + hep_hpc::hdf5::Column, // raw plane + hep_hpc::hdf5::Column, // raw wire + hep_hpc::hdf5::Column, // raw time + hep_hpc::hdf5::Column //cryostat > hitNtuple; ///< hit ntuple hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // g4 id - hep_hpc::hdf5::Column, // particle type - hep_hpc::hdf5::Column, // parent g4 id - hep_hpc::hdf5::Column, // is from nu - hep_hpc::hdf5::Column, // momentum - hep_hpc::hdf5::Column, // start position (x, y, z) - hep_hpc::hdf5::Column, // end position (x, y, z) - hep_hpc::hdf5::Column, // start process - hep_hpc::hdf5::Column // end process + hep_hpc::hdf5::Column, // g4 id + hep_hpc::hdf5::Column, // particle type + hep_hpc::hdf5::Column, // parent g4 id + hep_hpc::hdf5::Column, // is from nu + hep_hpc::hdf5::Column, // momentum + hep_hpc::hdf5::Column, // start position (x, y, z) + hep_hpc::hdf5::Column, // end position (x, y, z) + hep_hpc::hdf5::Column, // start process + hep_hpc::hdf5::Column // end process > particleNtuple; ///< particle ntuple hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // hit id - hep_hpc::hdf5::Column, // g4 id - hep_hpc::hdf5::Column, // deposited energy [ MeV ] - hep_hpc::hdf5::Column, // x position - hep_hpc::hdf5::Column, // y position - hep_hpc::hdf5::Column // z position + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // g4 id + hep_hpc::hdf5::Column, // deposited energy [ MeV ] + hep_hpc::hdf5::Column, // x position + hep_hpc::hdf5::Column, // y position + hep_hpc::hdf5::Column // z position > energyDepNtuple; ///< energy deposition ntuple hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // hit id - hep_hpc::hdf5::Column, // hit_channel - hep_hpc::hdf5::Column, // wire pos - hep_hpc::hdf5::Column, // peaktime - hep_hpc::hdf5::Column, // width - hep_hpc::hdf5::Column, // area - hep_hpc::hdf5::Column, // amplitude - hep_hpc::hdf5::Column, // pe - hep_hpc::hdf5::Column // usedInFlash + hep_hpc::hdf5::Column, // hit id + hep_hpc::hdf5::Column, // hit_channel + hep_hpc::hdf5::Column, // wire pos + hep_hpc::hdf5::Column, // peaktime + hep_hpc::hdf5::Column, // width + hep_hpc::hdf5::Column, // area + hep_hpc::hdf5::Column, // amplitude + hep_hpc::hdf5::Column, // pe + hep_hpc::hdf5::Column // usedInFlash > opHitNtuple; ///< PMT hit ntuple hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // flash id - hep_hpc::hdf5::Column, // wire pos - hep_hpc::hdf5::Column, // time - hep_hpc::hdf5::Column, // time width - hep_hpc::hdf5::Column, // Y center - hep_hpc::hdf5::Column, // Y width - hep_hpc::hdf5::Column, // Z center - hep_hpc::hdf5::Column, // Z width - hep_hpc::hdf5::Column // totalpe + hep_hpc::hdf5::Column, // flash id + hep_hpc::hdf5::Column, // wire pos + hep_hpc::hdf5::Column, // time + hep_hpc::hdf5::Column, // time width + hep_hpc::hdf5::Column, // Y center + hep_hpc::hdf5::Column, // Y width + hep_hpc::hdf5::Column, // Z center + hep_hpc::hdf5::Column, // Z width + hep_hpc::hdf5::Column // totalpe > opFlashNtuple; ///< Flash ntuple hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) - hep_hpc::hdf5::Column, // sumpe id - hep_hpc::hdf5::Column, // flash id - hep_hpc::hdf5::Column, // PMT channel - hep_hpc::hdf5::Column // pe + hep_hpc::hdf5::Column, // sumpe id + hep_hpc::hdf5::Column, // flash id + hep_hpc::hdf5::Column, // PMT channel + hep_hpc::hdf5::Column // pe > opFlashSumPENtuple; ///< Flash SumPE ntuple static std::string makeOutputFileName(std::string const& outputName, art::SubRunID const& sr) @@ -169,32 +165,20 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { << "_r" << std::setfill('0') << std::setw(5) << sr.run() << "_s" << std::setfill('0') << std::setw(5) << sr.subRun() << "_ts" << std::setw(6) << now.tv_usec << ".h5"; + std::cout << fileName.str() << std::endl; return fileName.str(); } - HDFDataFile(std::string const& outputName, std::string const& eventInfo, art::SubRunID const& sr) + HDFDataFile(std::string const& outputName, art::SubRunID const& sr) : file{ makeOutputFileName(outputName, sr), H5F_ACC_TRUNC } - , eventNtuple{ (eventInfo == "none") - ? std::optional > - >(hep_hpc::hdf5::make_ntuple({file, "event_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3)) - ) : std::nullopt - } - , eventNtupleNu{ (eventInfo == "nu") - ? std::optional, - hep_hpc::hdf5::Column, - hep_hpc::hdf5::Column, - hep_hpc::hdf5::Column, - hep_hpc::hdf5::Column, - hep_hpc::hdf5::Column > - >(hep_hpc::hdf5::make_ntuple({file, "event_table", 1000}, - hep_hpc::hdf5::make_column("event_id", 3), - hep_hpc::hdf5::make_scalar_column("is_cc"), - hep_hpc::hdf5::make_scalar_column("nu_pdg"), - hep_hpc::hdf5::make_scalar_column("nu_energy"), - hep_hpc::hdf5::make_scalar_column("lep_energy"), - hep_hpc::hdf5::make_column("nu_dir", 3)) - ) : std::nullopt + , eventNtupleNu{ + hep_hpc::hdf5::make_ntuple({file, "event_table", 1000}, + hep_hpc::hdf5::make_column("event_id", 3), + hep_hpc::hdf5::make_scalar_column("is_cc"), + hep_hpc::hdf5::make_scalar_column("nu_pdg"), + hep_hpc::hdf5::make_scalar_column("nu_energy"), + hep_hpc::hdf5::make_scalar_column("lep_energy"), + hep_hpc::hdf5::make_column("nu_dir", 3)) } , spacePointNtuple{ hep_hpc::hdf5::make_ntuple({file, "spacepoint_table", 1000}, @@ -217,7 +201,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::make_scalar_column("local_plane"), hep_hpc::hdf5::make_scalar_column("local_wire"), hep_hpc::hdf5::make_scalar_column("local_time"), - hep_hpc::hdf5::make_scalar_column("Cryostat")) + hep_hpc::hdf5::make_scalar_column("Cryostat")) } , particleNtuple{ hep_hpc::hdf5::make_ntuple({file, "particle_table", 1000}, @@ -230,7 +214,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::make_column("start_position", 3), hep_hpc::hdf5::make_column("end_position", 3), hep_hpc::hdf5::make_scalar_column("start_process"), - hep_hpc::hdf5::make_scalar_column("end_process")) + hep_hpc::hdf5::make_scalar_column("end_process")) } , energyDepNtuple{ hep_hpc::hdf5::make_ntuple({file, "edep_table", 1000}, @@ -240,7 +224,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::make_scalar_column("energy"), hep_hpc::hdf5::make_scalar_column("x_position"), hep_hpc::hdf5::make_scalar_column("y_position"), - hep_hpc::hdf5::make_scalar_column("z_position")) + hep_hpc::hdf5::make_scalar_column("z_position")) } , opHitNtuple{ hep_hpc::hdf5::make_ntuple({file, "ophit_table", 1000}, @@ -253,7 +237,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::make_scalar_column("area"), hep_hpc::hdf5:: make_scalar_column("amplitude"), hep_hpc::hdf5::make_scalar_column("pe"), - hep_hpc::hdf5::make_scalar_column("sumpe_id")) + hep_hpc::hdf5::make_scalar_column("sumpe_id")) } , opFlashNtuple{ hep_hpc::hdf5::make_ntuple({file, "opflash_table", 1000}, @@ -266,7 +250,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::make_scalar_column("y_width"), hep_hpc::hdf5::make_scalar_column("z_center"), hep_hpc::hdf5::make_scalar_column("z_width"), - hep_hpc::hdf5::make_scalar_column("totalpe")) + hep_hpc::hdf5::make_scalar_column("totalpe")) } , opFlashSumPENtuple{ hep_hpc::hdf5::make_ntuple({file, "opflashsumpe_table", 1000}, @@ -274,9 +258,9 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::make_scalar_column("sumpe_id"), hep_hpc::hdf5::make_scalar_column("flash_id"), hep_hpc::hdf5::make_scalar_column("pmt_channel"), - hep_hpc::hdf5::make_scalar_column("sumpe")) + hep_hpc::hdf5::make_scalar_column("sumpe")) } - {} + { } }; std::unique_ptr fHDFData; @@ -292,13 +276,8 @@ ICARUSHDF5Maker::ICARUSHDF5Maker(fhicl::ParameterSet const& p) fOpHitLabel( p.get("OpHitLabel")), fOpFlashLabel( p.get("OpFlashLabel")), fUseMap( p.get("UseMap", false)), - fEventInfo( p.get("EventInfo")), fOutputName(p.get("OutputName")) -{ - if (fEventInfo != "none" && fEventInfo != "nu") - throw art::Exception(art::errors::Configuration) - << "EventInfo must be \"none\" or \"nu\", not " << fEventInfo; -} +{ } void ICARUSHDF5Maker::analyze(art::Event const& e) { @@ -327,49 +306,40 @@ void ICARUSHDF5Maker::analyze(art::Event const& e) { std::array evtID { run, subrun, event }; // Fill event table - if (fEventInfo == "none") { - fHDFData->eventNtuple->insert( evtID.data() ); - mf::LogInfo("ICARUSHDF5Maker") << "Filling event table" - << "\nrun " << evtID[0] << ", subrun " << evtID[1] - << ", event " << evtID[2]; + // Get MC truth + auto truthHandle = e.getValidHandle>(fTruthLabel); + if (truthHandle->size() != 1) { + //avoid pile-up, which is not handled downstream + return; } - - if (fEventInfo == "nu") { - // Get MC truth - auto truthHandle = e.getValidHandle>(fTruthLabel); - if (truthHandle->size() != 1) { - throw art::Exception(art::errors::LogicError) - << "Expected to find exactly one MC truth object!"; - } - simb::MCNeutrino const& nutruth = truthHandle->at(0).GetNeutrino(); - - auto up = nutruth.Nu().Momentum().Vect().Unit(); - std::array nuMomentum {(float)up.X(),(float)up.Y(),(float)up.Z()}; - - fHDFData->eventNtupleNu->insert( evtID.data(), - nutruth.CCNC() == simb::kCC, - nutruth.Nu().PdgCode(), - nutruth.Nu().E(), - nutruth.Lepton().E(), - nuMomentum.data() - ); - - // for (int ip=0;ipat(0).NParticles();ip++) { - // std::cout << "mcp tkid=" << truthHandle->at(0).GetParticle(ip).TrackId() << " pdg=" << truthHandle->at(0).GetParticle(ip).PdgCode() - // << " mother=" << truthHandle->at(0).GetParticle(ip).Mother() - // << " vtx=" << truthHandle->at(0).GetParticle(ip).Vx() << " " << truthHandle->at(0).GetParticle(ip).Vy() << " " << truthHandle->at(0).GetParticle(ip).Vz() - // << std::endl; - // } - - mf::LogDebug("ICARUSHDF5Maker") << "Filling event table" - << "\nrun " << evtID[0] << ", subrun " << evtID[1] - << ", event " << evtID[2] - << "\nis cc? " << (nutruth.CCNC() == simb::kCC) - << ", nu energy " << nutruth.Nu().E() - << ", lepton energy " << nutruth.Lepton().E() - << "\nnu momentum x " << nuMomentum[0] << ", y " - << nuMomentum[1] << ", z " << nuMomentum[2]; - } // if nu event info + simb::MCNeutrino const& nutruth = truthHandle->at(0).GetNeutrino(); + + auto up = nutruth.Nu().Momentum().Vect().Unit(); + std::array nuMomentum {(float)up.X(),(float)up.Y(),(float)up.Z()}; + + fHDFData->eventNtupleNu.insert( evtID.data(), + nutruth.CCNC() == simb::kCC, + nutruth.Nu().PdgCode(), + nutruth.Nu().E(), + nutruth.Lepton().E(), + nuMomentum.data() + ); + + // for (int ip=0;ipat(0).NParticles();ip++) { + // std::cout << "mcp tkid=" << truthHandle->at(0).GetParticle(ip).TrackId() << " pdg=" << truthHandle->at(0).GetParticle(ip).PdgCode() + // << " mother=" << truthHandle->at(0).GetParticle(ip).Mother() + // << " vtx=" << truthHandle->at(0).GetParticle(ip).Vx() << " " << truthHandle->at(0).GetParticle(ip).Vy() << " " << truthHandle->at(0).GetParticle(ip).Vz() + // << std::endl; + // } + + mf::LogDebug("ICARUSHDF5Maker") << "Filling event table" + << "\nrun " << evtID[0] << ", subrun " << evtID[1] + << ", event " << evtID[2] + << "\nis cc? " << (nutruth.CCNC() == simb::kCC) + << ", nu energy " << nutruth.Nu().E() + << ", lepton energy " << nutruth.Lepton().E() + << "\nnu momentum x " << nuMomentum[0] << ", y " + << nuMomentum[1] << ", z " << nuMomentum[2]; // Get spacepoints from the event record auto spListHandle = e.getValidHandle>(fSPLabel); @@ -628,7 +598,7 @@ void ICARUSHDF5Maker::analyze(art::Event const& e) { } // function ICARUSHDF5Maker::analyze void ICARUSHDF5Maker::beginSubRun(art::SubRun const& sr) { - fHDFData = std::make_unique(fOutputName, fEventInfo, sr.id()); + fHDFData = std::make_unique(fOutputName, sr.id()); } void ICARUSHDF5Maker::endSubRun(art::SubRun const& sr) { diff --git a/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl b/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl index 6f4ac1047..decac0c63 100644 --- a/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl +++ b/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl @@ -43,7 +43,6 @@ physics: hdf5maker: { module_type: "ICARUSHDF5Maker" TruthLabel: "generator" - EventInfo: "nu" SPLabel: "nuslhits" HitLabel: "nuslhits" OpFlashLabel:"opflashCryoE" From 2f14f07403d892f0a7bcc0d98bdaf8b9735acf60 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Mon, 9 Jun 2025 22:24:44 -0500 Subject: [PATCH 19/24] fix assocSliceHitKeys usage --- .../TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc | 9 ++++++--- .../NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc | 12 ++++++++---- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc index 73c5a6a96..0eaf86b66 100644 --- a/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc @@ -18,6 +18,7 @@ #include // std::find() #include // std::abs() #include +#include #include // std::move() #include @@ -135,12 +136,14 @@ void ICARUSNuSliceHitsProducer::produce(art::Event& e) if (triggeringSlice) { mf::LogTrace{ "ICARUSNuSliceHitsProducer" } << "slice hits=" << assocSliceHit.at(triggeringSlice.key()).size(); auto const& hits = assocSliceHit.at(triggeringSlice.key()); - for (auto h : hits) { - assocSliceHitKeys->push_back(h.key()); + std::unordered_map keyIdxMap; + for (size_t ih=0; ihpush_back(hits[ih].key()); + keyIdxMap.insert({hits[ih].key(), ih}); } std::sort(begin(*assocSliceHitKeys), end(*assocSliceHitKeys)); outputHits->reserve(assocSliceHitKeys->size()); - for (auto const& hit: *assocSliceHitKeys) outputHits->push_back(*hits[hit]); + for (auto const& key: *assocSliceHitKeys) outputHits->push_back(*hits[keyIdxMap.at(key)]); } // Get spacepoints from the event record diff --git a/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc index 6a22493f0..1b85b3db1 100644 --- a/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSTrueNuSliceHitsProducer_module.cc @@ -13,8 +13,10 @@ #include "messagefacility/MessageLogger/MessageLogger.h" #include "canvas/Persistency/Common/FindManyP.h" -#include #include // std::find() +#include // std::abs() +#include +#include #include // std::move() #include @@ -133,12 +135,14 @@ void ICARUSTrueNuSliceHitsProducer::produce(art::Event& e) if (slkey>=0) { auto const& hits = assocSliceHit.at(slkey); mf::LogTrace{ "ICARUSTrueNuSliceHitsProducer" } << "slice hits=" << hits.size(); - for (auto const& h : hits) { - assocSliceHitKeys->push_back(h.key()); + std::unordered_map keyIdxMap; + for (size_t ih=0; ihpush_back(hits[ih].key()); + keyIdxMap.insert({hits[ih].key(), ih}); } std::sort(begin(*assocSliceHitKeys), end(*assocSliceHitKeys)); outputHits->reserve(assocSliceHitKeys->size()); - for (auto const& hit: *assocSliceHitKeys) outputHits->push_back(*hits[hit]); + for (auto const& key: *assocSliceHitKeys) outputHits->push_back(*hits[keyIdxMap.at(key)]); } // Get spacepoints from the event record From c3b784bf7520bf577ac691f897feaee17fe39f61 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Tue, 10 Jun 2025 17:37:36 -0500 Subject: [PATCH 20/24] restore correct behavior for allIDs, add PMT positions, fix SimChannelModuleLabel --- icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc | 15 ++++++++++----- icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl | 2 +- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc index a77e966c8..51cd52a20 100644 --- a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc @@ -139,7 +139,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::Ntuple, // event id (run, subrun, event) hep_hpc::hdf5::Column, // flash id - hep_hpc::hdf5::Column, // wire pos + hep_hpc::hdf5::Column, // wire pos hep_hpc::hdf5::Column, // time hep_hpc::hdf5::Column, // time width hep_hpc::hdf5::Column, // Y center @@ -153,6 +153,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::Column, // sumpe id hep_hpc::hdf5::Column, // flash id hep_hpc::hdf5::Column, // PMT channel + hep_hpc::hdf5::Column, // YZ pos hep_hpc::hdf5::Column // pe > opFlashSumPENtuple; ///< Flash SumPE ntuple @@ -258,6 +259,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::make_scalar_column("sumpe_id"), hep_hpc::hdf5::make_scalar_column("flash_id"), hep_hpc::hdf5::make_scalar_column("pmt_channel"), + hep_hpc::hdf5::make_column("yz_pos", 2), hep_hpc::hdf5::make_scalar_column("sumpe")) } { } @@ -349,7 +351,7 @@ void ICARUSHDF5Maker::analyze(art::Event const& e) { auto hitListHandle = e.getValidHandle>(fHitLabel); std::vector> hitlist; art::fill_ptr_vector(hitlist, hitListHandle); - + // Get assocations from spacepoints to hits art::FindManyP fmp(spListHandle, e, fSPLabel); @@ -464,8 +466,9 @@ void ICARUSHDF5Maker::analyze(art::Event const& e) { const simb::MCParticle* p = pi->TrackIdToParticle_P(abs(id)); allIDs.emplace(abs(id), p); while (p && p->Mother() != 0 ) { - p = pi->TrackIdToParticle_P(abs(p->Mother())); - allIDs.emplace(abs(p->Mother()), p); + auto mid = abs(p->Mother()); + p = pi->TrackIdToParticle_P(mid); + allIDs.emplace(mid, p); } } // Loop over true particles and fill table @@ -540,7 +543,9 @@ void ICARUSHDF5Maker::analyze(art::Event const& e) { std::vector sumpepmtmap(art::ServiceHandle()->NOpDets(),0); for (size_t ipmt=0;ipmtopFlashSumPENtuple.insert(evtID.data(),count,flkey,ipmt,pes[ipmt]); + auto xyz = geom.OpDetGeoFromOpChannel(ipmt).GetCenter(); + std::vector yzpos = {float(xyz.Y()),float(xyz.Z())}; + fHDFData->opFlashSumPENtuple.insert(evtID.data(),count,flkey,ipmt,yzpos.data(),pes[ipmt]); sumpepmtmap[ipmt] = count; count++; } diff --git a/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl b/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl index decac0c63..399d72c6a 100644 --- a/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl +++ b/icaruscode/TPC/NuGraph/hdf5maker_icarus_slice.fcl @@ -91,4 +91,4 @@ physics: } services.BackTrackerService.BackTracker.G4ModuleLabel: "largeant" -services.BackTrackerService.BackTracker.SimChannelModuleLabel: "daq:simpleSC" +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" From 438a62defe16e84e2f4d332560e117d2156d06e8 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Wed, 11 Jun 2025 17:34:32 -0500 Subject: [PATCH 21/24] addressing comments --- icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc | 10 ++++------ icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc | 11 ++++++++++- .../TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc | 3 ++- .../scripts/setup_tritonserver-nugraph-v0_grid.sh | 4 ++-- 4 files changed, 18 insertions(+), 10 deletions(-) diff --git a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc index 51cd52a20..18f9d480c 100644 --- a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc @@ -47,7 +47,6 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { public: explicit ICARUSHDF5Maker(fhicl::ParameterSet const& p); - ~ICARUSHDF5Maker() noexcept {}; // bare pointers are cleaned up by endSubRun ICARUSHDF5Maker(ICARUSHDF5Maker const&) = delete; ICARUSHDF5Maker(ICARUSHDF5Maker&&) = delete; @@ -70,8 +69,6 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { bool fUseMap; std::string fOutputName; - std::vector tpc_ids_checked = {}; // add the TPC ids here so we only check them once - struct HDFDataFile { hep_hpc::hdf5::File file; ///< output HDF5 file @@ -179,7 +176,7 @@ class ICARUSHDF5Maker : public art::EDAnalyzer { hep_hpc::hdf5::make_scalar_column("nu_pdg"), hep_hpc::hdf5::make_scalar_column("nu_energy"), hep_hpc::hdf5::make_scalar_column("lep_energy"), - hep_hpc::hdf5::make_column("nu_dir", 3)) + hep_hpc::hdf5::make_column("nu_dir", 3)) } , spacePointNtuple{ hep_hpc::hdf5::make_ntuple({file, "spacepoint_table", 1000}, @@ -358,7 +355,7 @@ void ICARUSHDF5Maker::analyze(art::Event const& e) { // Fill spacepoint table for (size_t i = 0; i < spListHandle->size(); ++i) { - std::vector> const& associatedHits(fmp.at(i)); + std::vector> const& associatedHits = fmp.at(i); if (associatedHits.size() < 3) { throw art::Exception(art::errors::LogicError) << "I am certain this cannot happen... but here you go, space point with " << associatedHits.size() << " hits"; } @@ -450,7 +447,8 @@ void ICARUSHDF5Maker::analyze(art::Event const& e) { std::vector const& h_ides = bt->ChannelToTrackIDEs(clockData, hit->Channel(), hit->StartTick(), hit->EndTick()); for (auto const& tide : h_ides) { int tid = tide.trackID; - // if negative id, make sure there is a corresponding particle to look for before taking the abs. This way negative means no asociated particle. + // if negative id, make sure there is a corresponding particle to look for before taking the abs. + // This way negative means no associated particle (a convention that can be used in the code that processes the ntuple). if (pi->TrackIdToParticle_P(abs(tid))) tid = abs(tid); g4id.insert(tid); fHDFData->energyDepNtuple.insert(evtID.data(),hit.key(), tid, tide.energyFrac, -1., -1., -1.); diff --git a/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc b/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc index a09b3b253..2cf0ad499 100644 --- a/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc +++ b/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc @@ -21,6 +21,15 @@ class ICARUSNuGraphLoader : public LoaderToolBase { +/** + * @brief Tool to collect the inputs needed by NuGraph from ICARUS data products. + * + * Read hit and space points from the event record, and package them for usage in NuGraph. This tool is called from larrecodnn/NuGraph/NuGraphInference_module.cc. + * Hits are stitched so that the 4 ICARUS TPCs in a cryostat are viewed in a single time vs wire plane. + * Only space points with chi2 fmp(spListHandle, e, fSpsInput); for (size_t spIdx = 0; spIdx < sp2Hit.size(); ++spIdx) { - static constexpr double MinChiSq = 0.5; ///< Threshold to consider a space point good. if (splist[spIdx]->Chisq()>MinChiSq) continue; // only space points with hits on all planes are enough for NuGraph if (fmp.at(spIdx).size()!=3) continue; diff --git a/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc index 0eaf86b66..9f669d589 100644 --- a/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc @@ -37,6 +37,7 @@ class ICARUSNuSliceHitsProducer : public art::EDProducer { * @brief Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction. * * Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction. + * The slice is the one idenfied as the one producing the trigger flash. Only space points made out of 3 hits are saved. * It also produces a vector of ints, that maps the new hit collection to the original one. * Optionally, it can store the association to the trigger flash used to pick the slice. * @@ -60,7 +61,7 @@ class ICARUSNuSliceHitsProducer : public art::EDProducer { art::InputTag const fBaryMatchLabel; art::InputTag const fSliceLabel; art::InputTag const fSpsLabel; - bool doFlashAssn; + bool const doFlashAssn; // Required functions. void produce(art::Event& e) override; diff --git a/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh b/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh index d6ed36a09..174a461d6 100755 --- a/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh +++ b/icaruscode/TPC/NuGraph/scripts/setup_tritonserver-nugraph-v0_grid.sh @@ -23,8 +23,8 @@ while 2>/dev/null >"/dev/tcp/0.0.0.0/$BASEPORT"; do export METRPORT=$((BASEPORT+2)) done # -echo "physics.producers.NuGraphCryoE.TritonConfig.serverURL: \"localhost:"$GRPCPORT"\"" >> $FCL -echo "physics.producers.NuGraphCryoW.TritonConfig.serverURL: \"localhost:"$GRPCPORT"\"" >> $FCL +echo "physics.producers.NuGraphCryoE.TritonConfig.serverURL: 'localhost:${GRPCPORT}'" >> "$FCL" +echo "physics.producers.NuGraphCryoW.TritonConfig.serverURL: 'localhost:${GRPCPORT}'" >> "$FCL" # /cvmfs/oasis.opensciencegrid.org/mis/apptainer/1.3.2/bin/apptainer exec --pid --ipc --home ~/:${HOME} --pwd ${PWD} /cvmfs/icarus.opensciencegrid.org/containers/tritonserver/nugraph-v0/ tritonserver --model-repository /triton-server-config/models --http-port=${HTTPPORT} --grpc-port=${GRPCPORT} --metrics-port=${METRPORT} >& tritonserver_nugraph-v0.log & # From 55e008729b165b5c393b6c17b56667605deefe0a Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Wed, 11 Jun 2025 17:41:06 -0500 Subject: [PATCH 22/24] addressing comments --- icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc | 3 --- 1 file changed, 3 deletions(-) diff --git a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc index 18f9d480c..c6f80bb6d 100644 --- a/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSHDF5Maker_module.cc @@ -369,9 +369,6 @@ void ICARUSHDF5Maker::analyze(art::Event const& e) { // << " plane2=" << plane << " key=" << hit.key() << std::endl; hitID[plane] = hit.key(); } - if (hitID[0]==0 && hitID[1]==-1 && hitID[2]==-1) { - throw art::Exception(art::errors::LogicError) << "THIS SHOULD NOT HAPPEN -- sp i=" << i << " hit ids=" << hitID[0] << " " << hitID[1] << " " << hitID[2]; - } fHDFData->spacePointNtuple.insert(evtID.data(),splist[i].ID(), pos.data(), hitID.data(),splist[i].Chisq() ); mf::LogDebug("ICARUSHDF5Maker") << "Filling spacepoint table" From d789fe9a0f46287437f76698ecaa7c62dafed949 Mon Sep 17 00:00:00 2001 From: Giuseppe Cerati Date: Wed, 11 Jun 2025 18:04:52 -0500 Subject: [PATCH 23/24] improve documentation --- icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc b/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc index 2cf0ad499..0774c17b5 100644 --- a/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc +++ b/icaruscode/TPC/NuGraph/ICARUSNuGraphLoader_tool.cc @@ -26,7 +26,7 @@ class ICARUSNuGraphLoader : public LoaderToolBase { * * Read hit and space points from the event record, and package them for usage in NuGraph. This tool is called from larrecodnn/NuGraph/NuGraphInference_module.cc. * Hits are stitched so that the 4 ICARUS TPCs in a cryostat are viewed in a single time vs wire plane. - * Only space points with chi2 Date: Wed, 11 Jun 2025 20:44:09 -0500 Subject: [PATCH 24/24] fix typo in description --- icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc index 9f669d589..3b48cb6d6 100644 --- a/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc +++ b/icaruscode/TPC/NuGraph/ICARUSNuSliceHitsProducer_module.cc @@ -37,7 +37,7 @@ class ICARUSNuSliceHitsProducer : public art::EDProducer { * @brief Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction. * * Produce separate collections of hits and space point from the slice identified as most likely from the neutrino interaction. - * The slice is the one idenfied as the one producing the trigger flash. Only space points made out of 3 hits are saved. + * The slice is the one identified as the one producing the trigger flash. Only space points made out of 3 hits are saved. * It also produces a vector of ints, that maps the new hit collection to the original one. * Optionally, it can store the association to the trigger flash used to pick the slice. *