diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 3ce38060a..4fb02d234 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -65,7 +65,7 @@ #include "art/Framework/Services/System/TriggerNamesService.h" #include "nurandom/RandomUtils/NuRandomService.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" -#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" #include "larcore/CoreUtils/ServiceUtil.h" #include "larevt/SpaceCharge/SpaceCharge.h" @@ -93,6 +93,7 @@ #include "lardataobj/RecoBase/PFParticle.h" #include "lardataobj/RecoBase/Slice.h" #include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/TrackHitMeta.h" #include "lardataobj/RecoBase/Vertex.h" #include "lardataobj/RecoBase/Shower.h" #include "lardataobj/RecoBase/MCSFitResult.h" @@ -1510,15 +1511,15 @@ void CAFMaker::produce(art::Event& evt) noexcept { } // Prep truth-to-reco-matching info - std::map>> id_to_ide_map; + std::map> id_to_ide_map; std::map>> id_to_truehit_map; std::map id_to_hit_energy_map; if ( !isRealData ) { art::ServiceHandle bt_serv; - id_to_ide_map = PrepSimChannels(simchannels, wireReadout); - id_to_truehit_map = PrepTrueHits(hits, clock_data, *bt_serv); + id_to_ide_map = sbn::PrepSimChannels(simchannels, wireReadout); + id_to_truehit_map = sbn::PrepTrueHits(hits, clock_data, *bt_serv); id_to_hit_energy_map = SetupIDHitEnergyMap(hits, clock_data, *bt_serv); } @@ -2291,8 +2292,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { FindManyPStrict(fmPFPart, evt, fParams.PFParticleLabel() + slice_tag_suff); - art::FindManyP fmTrackHit = - FindManyPStrict(slcTracks, evt, + art::FindManyP fmTrackHit = + FindManyPDStrict(slcTracks, evt, fParams.RecoTrackLabel() + slice_tag_suff); art::FindManyP fmShowerHit = @@ -2571,10 +2572,10 @@ void CAFMaker::produce(art::Event& evt) noexcept { FillTrackDazzle(fmTrackDazzle.at(iPart).front(), trk); } if (fmCalo.isValid()) { - FillTrackCalo(fmCalo.at(iPart), fmTrackHit.at(iPart), + FillTrackCalo(fmCalo.at(iPart), *thisTrack[0], fmTrackHit.at(iPart), fmTrackHit.data(iPart), (fParams.FillHitsNeutrinoSlices() && NeutrinoSlice) || fParams.FillHitsAllSlices(), fParams.TrackHitFillRRStartCut(), fParams.TrackHitFillRREndCut(), - dprop, trk); + dprop, *sce, trk); } if (fmCRTHitMatch.isValid() && fDet == kICARUS) { diff --git a/sbncode/CAFMaker/CMakeLists.txt b/sbncode/CAFMaker/CMakeLists.txt index 1dbde57ed..f0e59061a 100644 --- a/sbncode/CAFMaker/CMakeLists.txt +++ b/sbncode/CAFMaker/CMakeLists.txt @@ -30,6 +30,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker larcorealg::Geometry larcore::Geometry_Geometry_service larcorealg::GeoAlgo + lardataalg::headers larsim::MCCheater_BackTrackerService_service nusimdata::SimulationBase larsim::MCCheater_ParticleInventoryService_service @@ -72,6 +73,7 @@ cet_build_plugin ( CAFMaker art::module hep_concurrency::hep_concurrency lardataobj::RecoBase nurandom::RandomUtils_NuRandomService_service + lardata::DetectorPropertiesService lardata::DetectorClocksService BASENAME_ONLY ) diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index fe8c7cd78..2150395f0 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -1014,11 +1014,31 @@ namespace caf } } } + + // Helper function: get the e field + double GetEfield(const detinfo::DetectorPropertiesData& dprop, spacecharge::SpaceCharge const& sce, const geo::Point_t& loc) { + + double EField = dprop.Efield(); + if (sce.EnableSimEfieldSCE()) { + // Gets fractional E field distortions w.r.t. the nominal field on x-axis + geo::Vector_t EFieldOffsets = sce.GetEfieldOffsets(loc); + // Add 1 in X direction as this is the direction of the drift field, not caring if it is +x or -x direction, since we only want |E| + EFieldOffsets += geo::Vector_t{1, 0, 0}; + // Convert to Absolute E Field from relative + EFieldOffsets *= EField; + // We only care about the magnitude for recombination + EField = EFieldOffsets.r(); + } + return EField; + } void FillTrackPlaneCalo(const anab::Calorimetry &calo, + const recob::Track& track, const std::vector> &hits, + const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, + const spacecharge::SpaceCharge &sce, caf::SRTrackCalo &srcalo) { // Collect info from Calorimetry @@ -1057,7 +1077,8 @@ namespace caf // lookup the wire -- the Calorimery object makes this // __way__ harder than it should be - for (const art::Ptr &h: hits) { + for (unsigned i_hit = 0; i_hit < hits.size(); i_hit++) { + const art::Ptr &h = hits[i_hit]; if (h.key() == tps[i]) { p.wire = h->WireID().Wire; p.tpc = h->WireID().TPC; @@ -1069,6 +1090,19 @@ namespace caf p.mult = h->Multiplicity(); p.start = h->StartTick(); p.end = h->EndTick(); + + + // Get the trajectory point index from this hit. Again -- this is too hard. + // + // Use this to get the (SCE corrected) efield and the angle to the drift direction + unsigned traj_point_index = thms.at(i_hit)->Index(); + if (track.HasValidPoint(traj_point_index)) { + float costh_drift = track.DirectionAtPoint(traj_point_index).X(); + float phi = acos(abs(costh_drift)); + float efield = GetEfield(dprop, sce, track.LocationAtPoint(traj_point_index)); + p.efield = efield; + p.phi = phi; + } } } @@ -1122,9 +1156,12 @@ namespace caf } void FillTrackCalo(const std::vector> &calos, + const recob::Track& track, const std::vector> &hits, + const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, + const spacecharge::SpaceCharge &sce, caf::SRTrack& srtrack, bool allowEmpty) { @@ -1137,7 +1174,7 @@ namespace caf if (calo.PlaneID()) { unsigned plane_id = calo.PlaneID().Plane; assert(plane_id < 3); - FillTrackPlaneCalo(calo, hits, fill_calo_points, fillhit_rrstart, fillhit_rrend, dprop, srtrack.calo[plane_id]); + FillTrackPlaneCalo(calo, track, hits, thms, fill_calo_points, fillhit_rrstart, fillhit_rrend, dprop, sce, srtrack.calo[plane_id]); } } diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h index a837e0b56..786b5dd0f 100644 --- a/sbncode/CAFMaker/FillReco.h +++ b/sbncode/CAFMaker/FillReco.h @@ -5,15 +5,18 @@ #include "art/Framework/Services/Registry/ServiceHandle.h" #include "larsim/MCCheater/ParticleInventoryService.h" -#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" // LArSoft includes #include "larcorealg/Geometry/fwd.h" +#include "larevt/SpaceCharge/SpaceCharge.h" + #include "lardataobj/RecoBase/PFParticle.h" #include "lardataobj/RecoBase/Shower.h" #include "lardataobj/RecoBase/Slice.h" #include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/TrackHitMeta.h" #include "lardataobj/RecoBase/Vertex.h" #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/SpacePoint.h" @@ -239,9 +242,12 @@ namespace caf bool allowEmpty = false); void FillTrackPlaneCalo(const anab::Calorimetry &calo, + const recob::Track& track, const std::vector> &hits, + const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, + spacecharge::SpaceCharge const& sce, caf::SRTrackCalo &srcalo); void FillTrackScatterClosestApproach(const art::Ptr closestapproach, @@ -257,9 +263,12 @@ namespace caf bool allowEmpty = false); void FillTrackCalo(const std::vector> &calos, + const recob::Track& track, const std::vector> &hits, + const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, + spacecharge::SpaceCharge const& sce, caf::SRTrack& srtrack, bool allowEmpty = false); diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 4d1d794d0..4c411b6d0 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -5,7 +5,6 @@ #include "larcorealg/Geometry/WireReadoutGeom.h" #include "larevt/SpaceCharge/SpaceCharge.h" #include "larcore/CoreUtils/ServiceUtil.h" -#include "RecoUtils/RecoUtils.h" #include "CLHEP/Random/RandGauss.h" @@ -136,7 +135,7 @@ namespace caf { }//FillTrackTruth // Assumes truth matching and calo-points are filled - void FillTrackCaloTruth(const std::map>> &id_to_ide_map, + void FillTrackCaloTruth(const std::map> &id_to_ide_map, const std::vector &mc_particles, const geo::GeometryCore& geometry, const geo::WireReadoutGeom& wireReadout, @@ -163,10 +162,10 @@ namespace caf { // Load the hits // match on the channel, which is unique - const std::vector> &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID); + const std::vector &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID); std::map> chan_2_ides; - for (auto const &ide_pair: match_ides) { - chan_2_ides[wireReadout.PlaneWireToChannel(ide_pair.first)].push_back(ide_pair.second); + for (auto const &ide_p: match_ides) { + chan_2_ides[wireReadout.PlaneWireToChannel(ide_p.wire)].push_back(ide_p.ide); } // pre-compute partial ranges @@ -634,15 +633,15 @@ namespace caf { void FillTrueG4Particle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, const std::vector> &neutrinos, caf::SRTrueParticle &srparticle) { - std::vector> empty; - const std::vector> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + std::vector empty; + const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; std::vector> emptyHits; const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; @@ -662,9 +661,9 @@ namespace caf { } } - for (auto const &ide_pair: particle_ides) { - const geo::WireID &w = ide_pair.first; - const sim::IDE *ide = ide_pair.second; + for (auto const &ide_p: particle_ides) { + const geo::WireID &w = ide_p.wire; + const sim::IDE *ide = ide_p.ide; if(w.Plane >= 0 && w.Plane < 3 && w.Cryostat < 2){ srparticle.plane[w.Cryostat][w.Plane].visE += ide->energy / 1000. /* MeV -> GeV*/; @@ -701,7 +700,7 @@ namespace caf { } // get the wall if (entry_point > 0) { - srparticle.wallin = GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); + srparticle.wallin = sbn::GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); } int exit_point = -1; @@ -781,7 +780,7 @@ namespace caf { exit_point++; // to avoid exactly the same start and end positions when single index is inside the active volumne } if (exit_point >= 0 && ((unsigned)exit_point) < particle.NumberTrajectoryPoints() - 1) { - srparticle.wallout = GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); + srparticle.wallout = sbn::GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); } // other truth information @@ -802,8 +801,8 @@ namespace caf { srparticle.endp = (exit_point >= 0) ? particle.Momentum(exit_point).Vect() : TVector3(-9999, -9999, -9999); srparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; - srparticle.start_process = GetG4ProcessID(particle.Process()); - srparticle.end_process = GetG4ProcessID(particle.EndProcess()); + srparticle.start_process = sbn::GetG4ProcessID(particle.Process()); + srparticle.end_process = sbn::GetG4ProcessID(particle.EndProcess()); srparticle.G4ID = particle.TrackId(); srparticle.parent = particle.Mother(); @@ -874,37 +873,6 @@ namespace caf { return ret; } - std::map>> PrepTrueHits(const std::vector> &allHits, - const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker) { - std::map>> ret; - for (const art::Ptr h: allHits) { - for (int ID: backtracker.HitToTrackIds(clockData, *h)) { - ret[abs(ID)].push_back(h); - } - } - return ret; - } - - std::map>> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { - std::map>> ret; - - for (const art::Ptr sc : simchannels) { - // Lookup the wire of this channel - raw::ChannelID_t channel = sc->Channel(); - std::vector maybewire = wireReadout.ChannelToWire(channel); - geo::WireID thisWire; // Default constructor makes invalid wire - if (maybewire.size()) thisWire = maybewire[0]; - - for (const auto &item : sc->TDCIDEMap()) { - for (const sim::IDE &ide: item.second) { - // indexing initializes empty vector - ret[abs(ide.trackID)].push_back({thisWire, &ide}); - } - } - } - return ret; - } - } // end namespace @@ -1256,135 +1224,6 @@ bool FRFillNueCC(const simb::MCTruth &mctruth, return true; } -caf::Wall_t caf::GetWallCross(const geo::BoxBoundedGeo &volume, const TVector3 p0, const TVector3 p1) { - TVector3 direction = (p1 - p0) * ( 1. / (p1 - p0).Mag()); - std::vector intersections = volume.GetIntersections(p0, direction); - - /* - * There are either two, or one, or zero intersection points. - * As per larcorealg/Geometry/BoxBoundedGeo.h documentation: - * If the return std::vector is empty the trajectory does not intersect with the box. - * Normally the return value should have one (if the trajectory originates in the box) or two (else) entries. - * If the return value has two entries the first represents the entry point and the second the exit point - */ - - if( intersections.empty() ) return caf::kWallNone; - - // get the intersection point closer to p0 - TVector3 closestIntersection; - double minDistance2 = std::numeric_limits::max(); - - for( TVector3 const & point : intersections ) { - const double d2 = (point - p0).Mag2(); - if( d2 > minDistance2 ) continue; - minDistance2 = d2; - closestIntersection = point; - } - - double eps = 1e-3; - if (abs(closestIntersection.X() - volume.MinX()) < eps) { - return caf::kWallLeft; - } - else if (abs(closestIntersection.X() - volume.MaxX()) < eps) { - return caf::kWallRight; - } - else if (abs(closestIntersection.Y() - volume.MinY()) < eps) { - return caf::kWallBottom; - } - else if (abs(closestIntersection.Y() - volume.MaxY()) < eps) { - return caf::kWallTop; - } - else if (abs(closestIntersection.Z() - volume.MinZ()) < eps) { - return caf::kWallFront; - } - else if (abs(closestIntersection.Z() - volume.MaxZ()) < eps) { - return caf::kWallBack; - } - else assert(false); - - return caf::kWallNone; -}//GetWallCross - -//------------------------------------------ - -caf::g4_process_ caf::GetG4ProcessID(const std::string &process_name) { -#define MATCH_PROCESS(name) if (process_name == #name) {return caf::kG4 ## name;} -#define MATCH_PROCESS_NAMED(strname, id) if (process_name == #strname) {return caf::kG4 ## id;} - MATCH_PROCESS(primary) - MATCH_PROCESS(CoupledTransportation) - MATCH_PROCESS(FastScintillation) - MATCH_PROCESS(Decay) - MATCH_PROCESS(anti_neutronInelastic) - MATCH_PROCESS(neutronInelastic) - MATCH_PROCESS(anti_protonInelastic) - MATCH_PROCESS(protonInelastic) - MATCH_PROCESS(hadInelastic) - MATCH_PROCESS_NAMED(kaon+Inelastic, kaonpInelastic) - MATCH_PROCESS_NAMED(kaon-Inelastic, kaonmInelastic) - MATCH_PROCESS_NAMED(kaon+Inelastic, kaonpInelastic) - MATCH_PROCESS_NAMED(kaon-Inelastic, kaonmInelastic) - MATCH_PROCESS_NAMED(sigma+Inelastic, sigmapInelastic) - MATCH_PROCESS_NAMED(sigma-Inelastic, sigmamInelastic) - MATCH_PROCESS_NAMED(pi+Inelastic, pipInelastic) - MATCH_PROCESS_NAMED(pi-Inelastic, pimInelastic) - MATCH_PROCESS_NAMED(xi+Inelastic, xipInelastic) - MATCH_PROCESS_NAMED(xi-Inelastic, ximInelastic) - MATCH_PROCESS(kaon0LInelastic) - MATCH_PROCESS(kaon0SInelastic) - MATCH_PROCESS(lambdaInelastic) - MATCH_PROCESS_NAMED(anti-lambdaInelastic, anti_lambdaInelastic) - MATCH_PROCESS(He3Inelastic) - MATCH_PROCESS(ionInelastic) - MATCH_PROCESS(xi0Inelastic) - MATCH_PROCESS(alphaInelastic) - MATCH_PROCESS(tInelastic) - MATCH_PROCESS(dInelastic) - MATCH_PROCESS(anti_neutronElastic) - MATCH_PROCESS(neutronElastic) - MATCH_PROCESS(anti_protonElastic) - MATCH_PROCESS(protonElastic) - MATCH_PROCESS(hadElastic) - MATCH_PROCESS_NAMED(kaon+Elastic, kaonpElastic) - MATCH_PROCESS_NAMED(kaon-Elastic, kaonmElastic) - MATCH_PROCESS_NAMED(pi+Elastic, pipElastic) - MATCH_PROCESS_NAMED(pi-Elastic, pimElastic) - MATCH_PROCESS(conv) - MATCH_PROCESS(phot) - MATCH_PROCESS(annihil) - MATCH_PROCESS(nCapture) - MATCH_PROCESS(nKiller) - MATCH_PROCESS(muMinusCaptureAtRest) - MATCH_PROCESS(muIoni) - MATCH_PROCESS(eBrem) - MATCH_PROCESS(CoulombScat) - MATCH_PROCESS(hBertiniCaptureAtRest) - MATCH_PROCESS(hFritiofCaptureAtRest) - MATCH_PROCESS(photonNuclear) - MATCH_PROCESS(muonNuclear) - MATCH_PROCESS(electronNuclear) - MATCH_PROCESS(positronNuclear) - MATCH_PROCESS(compt) - MATCH_PROCESS(eIoni) - MATCH_PROCESS(muBrems) - MATCH_PROCESS(hIoni) - MATCH_PROCESS(ionIoni) - MATCH_PROCESS(hBrems) - MATCH_PROCESS(muPairProd) - MATCH_PROCESS(hPairProd) - MATCH_PROCESS(LArVoxelReadoutScoringProcess) - MATCH_PROCESS(Transportation) - MATCH_PROCESS(msc) - MATCH_PROCESS(StepLimiter) - MATCH_PROCESS(RadioactiveDecayBase) - std::cerr << "Error: Process name with no match (" << process_name << ")\n"; - assert(false); - return caf::kG4UNKNOWN; // unreachable in debug mode -#undef MATCH_PROCESS -#undef MATCH_PROCESS_NAMED - -}//GetG4ProcessID -//------------------------------------------- - float ContainedLength(const TVector3 &v0, const TVector3 &v1, const std::vector &boxes) { static const geoalgo::GeoAlgo algo; diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index 419ed6ee9..d4134d55c 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -33,6 +33,8 @@ #include "sbnanaobj/StandardRecord/StandardRecord.h" #include "sbnanaobj/StandardRecord/SRMeVPrtl.h" +#include "RecoUtils/RecoUtils.h" + namespace caf { struct HitsEnergy { @@ -40,13 +42,6 @@ namespace caf float totE; }; - // Helpers - caf::Wall_t GetWallCross( const geo::BoxBoundedGeo &volume, - const TVector3 p0, - const TVector3 p1); - - caf::g4_process_ GetG4ProcessID(const std::string &name); - void FillSRGlobal(const sbn::evwgh::EventWeightParameterSet& pset, caf::SRGlobal& srglobal, std::map& weightPSetIndex); @@ -73,7 +68,7 @@ namespace caf void FillTrueG4Particle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, @@ -103,7 +98,7 @@ namespace caf caf::SRTrack& srtrack, bool allowEmpty = false); - void FillTrackCaloTruth(const std::map>> &id_to_ide_map, + void FillTrackCaloTruth(const std::map> &id_to_ide_map, const std::vector &mc_particles, const geo::GeometryCore & geometry, const geo::WireReadoutGeom& wireReadout, @@ -133,9 +128,6 @@ namespace caf CLHEP::HepRandomEngine &rand, std::vector &srfakereco); - std::map>> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); - std::map>> PrepTrueHits(const std::vector> &allHits, - const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); std::map SetupIDHitEnergyMap(const std::vector> &allHits, const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); diff --git a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc index d8e9f210f..a2cdc8922 100644 --- a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc +++ b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc @@ -107,3 +107,151 @@ int CAFRecoUtils::GetShowerPrimary(const int g4ID) return primary_id; } + +std::map> sbn::PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { + std::map> ret; + + for (const art::Ptr sc : simchannels) { + // Lookup the wire of this channel + raw::ChannelID_t channel = sc->Channel(); + std::vector maybewire = wireReadout.ChannelToWire(channel); + geo::WireID thisWire; // Default constructor makes invalid wire + if (maybewire.size()) thisWire = maybewire[0]; + + for (const auto &item : sc->TDCIDEMap()) { + for (const sim::IDE &ide: item.second) { + // indexing initializes empty vector + ret[abs(ide.trackID)].push_back({thisWire, item.first, &ide}); + } + } + } + return ret; +} + +std::map>> sbn::PrepTrueHits(const std::vector> &allHits, + const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker) { + std::map>> ret; + for (const art::Ptr h: allHits) { + for (int ID: backtracker.HitToTrackIds(clockData, *h)) { + ret[abs(ID)].push_back(h); + } + } + return ret; +} + +caf::Wall_t sbn::GetWallCross(const geo::BoxBoundedGeo &volume, const TVector3 p0, const TVector3 p1) { + TVector3 direction = (p1 - p0) * ( 1. / (p1 - p0).Mag()); + std::vector intersections = volume.GetIntersections(p0, direction); + + assert(intersections.size() == 2); + + // get the intersection point closer to p0 + int intersection_i = ((intersections[0] - p0).Mag() < (intersections[1] - p0).Mag()) ? 0 : 1; + + double eps = 1e-3; + if (abs(intersections[intersection_i].X() - volume.MinX()) < eps) { + //std::cout << "Left\n"; + return caf::kWallLeft; + } + else if (abs(intersections[intersection_i].X() - volume.MaxX()) < eps) { + //std::cout << "Right\n"; + return caf::kWallRight; + } + else if (abs(intersections[intersection_i].Y() - volume.MinY()) < eps) { + //std::cout << "Bottom\n"; + return caf::kWallBottom; + } + else if (abs(intersections[intersection_i].Y() - volume.MaxY()) < eps) { + //std::cout << "Top\n"; + return caf::kWallTop; + } + else if (abs(intersections[intersection_i].Z() - volume.MinZ()) < eps) { + //std::cout << "Front\n"; + return caf::kWallFront; + } + else if (abs(intersections[intersection_i].Z() - volume.MaxZ()) < eps) { + //std::cout << "Back\n"; + return caf::kWallBack; + } + else assert(false); + //std::cout << "None\n"; + + return caf::kWallNone; +}//GetWallCross + +//------------------------------------------ +caf::g4_process_ sbn::GetG4ProcessID(const std::string &process_name) { +#define MATCH_PROCESS(name) if (process_name == #name) {return caf::kG4 ## name;} +#define MATCH_PROCESS_NAMED(strname, id) if (process_name == #strname) {return caf::kG4 ## id;} + MATCH_PROCESS(primary) + MATCH_PROCESS(CoupledTransportation) + MATCH_PROCESS(FastScintillation) + MATCH_PROCESS(Decay) + MATCH_PROCESS(anti_neutronInelastic) + MATCH_PROCESS(neutronInelastic) + MATCH_PROCESS(anti_protonInelastic) + MATCH_PROCESS(protonInelastic) + MATCH_PROCESS(hadInelastic) + MATCH_PROCESS_NAMED(kaon+Inelastic, kaonpInelastic) + MATCH_PROCESS_NAMED(kaon-Inelastic, kaonmInelastic) + MATCH_PROCESS_NAMED(kaon+Inelastic, kaonpInelastic) + MATCH_PROCESS_NAMED(kaon-Inelastic, kaonmInelastic) + MATCH_PROCESS_NAMED(sigma+Inelastic, sigmapInelastic) + MATCH_PROCESS_NAMED(sigma-Inelastic, sigmamInelastic) + MATCH_PROCESS_NAMED(pi+Inelastic, pipInelastic) + MATCH_PROCESS_NAMED(pi-Inelastic, pimInelastic) + MATCH_PROCESS_NAMED(xi+Inelastic, xipInelastic) + MATCH_PROCESS_NAMED(xi-Inelastic, ximInelastic) + MATCH_PROCESS(kaon0LInelastic) + MATCH_PROCESS(kaon0SInelastic) + MATCH_PROCESS(lambdaInelastic) + MATCH_PROCESS_NAMED(anti-lambdaInelastic, anti_lambdaInelastic) + MATCH_PROCESS(He3Inelastic) + MATCH_PROCESS(ionInelastic) + MATCH_PROCESS(xi0Inelastic) + MATCH_PROCESS(alphaInelastic) + MATCH_PROCESS(tInelastic) + MATCH_PROCESS(dInelastic) + MATCH_PROCESS(anti_neutronElastic) + MATCH_PROCESS(neutronElastic) + MATCH_PROCESS(anti_protonElastic) + MATCH_PROCESS(protonElastic) + MATCH_PROCESS(hadElastic) + MATCH_PROCESS_NAMED(kaon+Elastic, kaonpElastic) + MATCH_PROCESS_NAMED(kaon-Elastic, kaonmElastic) + MATCH_PROCESS_NAMED(pi+Elastic, pipElastic) + MATCH_PROCESS_NAMED(pi-Elastic, pimElastic) + MATCH_PROCESS(conv) + MATCH_PROCESS(phot) + MATCH_PROCESS(annihil) + MATCH_PROCESS(nCapture) + MATCH_PROCESS(nKiller) + MATCH_PROCESS(muMinusCaptureAtRest) + MATCH_PROCESS(muIoni) + MATCH_PROCESS(eBrem) + MATCH_PROCESS(CoulombScat) + MATCH_PROCESS(hBertiniCaptureAtRest) + MATCH_PROCESS(hFritiofCaptureAtRest) + MATCH_PROCESS(photonNuclear) + MATCH_PROCESS(muonNuclear) + MATCH_PROCESS(electronNuclear) + MATCH_PROCESS(positronNuclear) + MATCH_PROCESS(compt) + MATCH_PROCESS(eIoni) + MATCH_PROCESS(muBrems) + MATCH_PROCESS(hIoni) + MATCH_PROCESS(ionIoni) + MATCH_PROCESS(hBrems) + MATCH_PROCESS(muPairProd) + MATCH_PROCESS(hPairProd) + MATCH_PROCESS(LArVoxelReadoutScoringProcess) + MATCH_PROCESS(Transportation) + MATCH_PROCESS(msc) + MATCH_PROCESS(StepLimiter) + MATCH_PROCESS(RadioactiveDecayBase) + std::cerr << "Error: Process name with no match (" << process_name << ")\n"; + assert(false); + return caf::kG4UNKNOWN; // unreachable in debug mode +#undef MATCH_PROCESS +#undef MATCH_PROCESS_NAMED +}//GetG4ProcessID diff --git a/sbncode/CAFMaker/RecoUtils/RecoUtils.h b/sbncode/CAFMaker/RecoUtils/RecoUtils.h index a99121d8e..0b43760a4 100644 --- a/sbncode/CAFMaker/RecoUtils/RecoUtils.h +++ b/sbncode/CAFMaker/RecoUtils/RecoUtils.h @@ -24,11 +24,13 @@ #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/Track.h" #include "lardataobj/RecoBase/Shower.h" +#include "sbnanaobj/StandardRecord/SREnums.h" //#include "lardataobj/AnalysisBase/MVAPIDResult.h" //#include "lardataobj/AnalysisBase/ParticleID.h" #include "larsim/MCCheater/BackTrackerService.h" #include "larsim/MCCheater/ParticleInventoryService.h" #include "larcore/Geometry/Geometry.h" +#include "larcorealg/Geometry/WireReadoutGeom.h" // c++ #include @@ -37,6 +39,25 @@ // ROOT #include "TTree.h" +namespace sim { + class IDE; +} + +namespace sbn { + struct ReadoutIDE { + geo::WireID wire; ///< Wire on a given plane closest to the drift path of the charge. + unsigned short tick = 0; ///< Time tick at which the charge passes closest to the wire. + const sim::IDE *ide = nullptr; ///< Deposited charge information. + }; + + // Helpers + std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); + std::map>> PrepTrueHits(const std::vector> &allHits, + const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); + caf::Wall_t GetWallCross( const geo::BoxBoundedGeo &volume, const TVector3 p0, const TVector3 p1); + caf::g4_process_ GetG4ProcessID(const std::string &name); +} + namespace CAFRecoUtils{ std::vector> AllTrueParticleIDEnergyMatches(const detinfo::DetectorClocksData &clockData, const std::vector >& hits, bool rollup_unsaved_ids=1); @@ -47,4 +68,7 @@ namespace CAFRecoUtils{ int GetShowerPrimary(const int g4ID); } + + + #endif diff --git a/sbncode/CMakeLists.txt b/sbncode/CMakeLists.txt index daed2590a..14a8e0265 100644 --- a/sbncode/CMakeLists.txt +++ b/sbncode/CMakeLists.txt @@ -20,6 +20,7 @@ add_subdirectory(GeometryTools) add_subdirectory(CosmicID) add_subdirectory(DetSim) add_subdirectory(Cluster3D) +add_subdirectory(WireMod) add_subdirectory(HitFinder) # Supera diff --git a/sbncode/Calibration/CMakeLists.txt b/sbncode/Calibration/CMakeLists.txt index 64439a086..a101a3069 100644 --- a/sbncode/Calibration/CMakeLists.txt +++ b/sbncode/Calibration/CMakeLists.txt @@ -25,6 +25,7 @@ cet_build_plugin( TrackCaloSkimmer art::module sbncode_CAFMaker sbnobj::Common_Calibration_dict larevt::SpaceCharge + caf_RecoUtils ) cet_build_plugin(TrackCaloSkimmerSelectStoppingTrack art::tool diff --git a/sbncode/Calibration/TrackCaloSkimmer.h b/sbncode/Calibration/TrackCaloSkimmer.h index 497959518..549bd8cfd 100644 --- a/sbncode/Calibration/TrackCaloSkimmer.h +++ b/sbncode/Calibration/TrackCaloSkimmer.h @@ -51,8 +51,8 @@ #include "larevt/SpaceCharge/SpaceCharge.h" #include "larevt/SpaceChargeServices/SpaceChargeService.h" -#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" #include "larcorealg/Geometry/fwd.h" @@ -68,6 +68,9 @@ #include "ITCSSelectionTool.h" +// Useful functions +#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h" // sbn::ReadoutIDE, sbn::PrepTrueHits()... + namespace sbn { class TrackCaloSkimmer; enum EDet {kNOTDEFINED, kSBND, kICARUS}; @@ -117,7 +120,6 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { int ID; }; - // Represents a "Snippet" of ADCs shared by a set of hits on a wire struct Snippet { geo::WireID wire; @@ -166,7 +168,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { const std::vector> &mcparticles, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> id_to_ide_map, + const std::map> id_to_ide_map, const std::map>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc index 00a02b7fc..8a639fc17 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -11,13 +11,11 @@ //////////////////////////////////////////////////////////////////////// #include "TrackCaloSkimmer.h" +#include "sbnobj/SBND/CRT/CRTTrack.hh" +#include "sbnanaobj/StandardRecord/SREnums.h" #include "art/Utilities/make_tool.h" -// Useful functions -#include "sbncode/CAFMaker/FillTrue.h" -#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h" - #include "larcore/Geometry/WireReadout.h" #include "larcore/Geometry/Geometry.h" #include "larcore/CoreUtils/ServiceUtil.h" @@ -291,15 +289,15 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) // Prep truth-to-reco-matching info // - // Use helper functions from CAFMaker/FillTrue - std::map>> id_to_ide_map; + // Use helper functions from CAFMaker/RecoUtils + std::map> id_to_ide_map; std::map>> id_to_truehit_map; const cheat::BackTrackerService *bt = NULL; if (simchannels.size()) { art::ServiceHandle bt_serv; - id_to_ide_map = caf::PrepSimChannels(simchannels, *wireReadout); - id_to_truehit_map = caf::PrepTrueHits(allHits, clock_data, *bt_serv.get()); + id_to_ide_map = sbn::PrepSimChannels(simchannels, *wireReadout); + id_to_truehit_map = sbn::PrepTrueHits(allHits, clock_data, *bt_serv.get()); bt = bt_serv.get(); } @@ -582,14 +580,14 @@ geo::Point_t WireToTrajectoryPosition(const geo::Point_t &loc, const geo::TPCID sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, const geo::WireReadoutGeom *wireReadout) { - std::vector> empty; - const std::vector> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + std::vector empty; + const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; std::vector> emptyHits; const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; @@ -606,9 +604,9 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, trueparticle.plane0nhit = 0; trueparticle.plane1nhit = 0; trueparticle.plane2nhit = 0; - for (auto const &ide_pair: particle_ides) { - const geo::WireID &w = ide_pair.first; - const sim::IDE *ide = ide_pair.second; + for (auto const &ide_p: particle_ides) { + const geo::WireID &w = ide_p.wire; + const sim::IDE *ide = ide_p.ide; if (w.Plane == 0) { trueparticle.plane0VisE += ide->energy / 1000. /* MeV -> GeV*/; @@ -659,7 +657,7 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, // get the wall if (entry_point > 0) { - trueparticle.wallin = (int)caf::GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); + trueparticle.wallin = (int)sbn::GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); } int exit_point = -1; @@ -727,7 +725,7 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, exit_point = particle.NumberTrajectoryPoints() - 1; } if (exit_point >= 0 && ((unsigned)exit_point) < particle.NumberTrajectoryPoints() - 1) { - trueparticle.wallout = (int)caf::GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); + trueparticle.wallout = (int)sbn::GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); } // other truth information @@ -748,8 +746,8 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, trueparticle.endp = ConvertTVector((exit_point >= 0) ? particle.Momentum(exit_point).Vect() : TVector3(-9999, -9999, -9999)); trueparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; - trueparticle.start_process = (int)caf::GetG4ProcessID(particle.Process()); - trueparticle.end_process = (int)caf::GetG4ProcessID(particle.EndProcess()); + trueparticle.start_process = (int)sbn::GetG4ProcessID(particle.Process()); + trueparticle.end_process = (int)sbn::GetG4ProcessID(particle.EndProcess()); trueparticle.G4ID = particle.TrackId(); trueparticle.parent = particle.Mother(); @@ -757,10 +755,11 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, // Organize deposition info into per-wire true "Hits" -- key is the Channel Number std::map truehits; - for (auto const &ide_pair: particle_ides) { - const geo::WireID &w = ide_pair.first; + for (auto const &ide_part: particle_ides) { + const geo::WireID &w = ide_part.wire; unsigned c = wireReadout->PlaneWireToChannel(w); - const sim::IDE *ide = ide_pair.second; + const sim::IDE *ide = ide_part.ide; + unsigned short tick = ide_part.tick; // Set stuff truehits[c].cryo = w.Cryostat; @@ -776,6 +775,8 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, truehits[c].p.y = (truehits[c].p.y*old_elec + ide->y*ide->numElectrons) / new_elec; truehits[c].p.z = (truehits[c].p.z*old_elec + ide->z*ide->numElectrons) / new_elec; + truehits[c].time = (truehits[c].time*old_elec + tick*ide->numElectrons) / new_elec; + // Also get the position with space charge un-done geo::Point_t ide_p(ide->x, ide->y, ide->z); geo::Point_t ide_p_scecorr = WireToTrajectoryPosition(ide_p, w); @@ -791,10 +792,10 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, } // Compute widths - for (auto const &ide_pair: particle_ides) { - const geo::WireID &w = ide_pair.first; + for (auto const &ide_part: particle_ides) { + const geo::WireID &w = ide_part.wire; unsigned c = wireReadout->PlaneWireToChannel(w); - const sim::IDE *ide = ide_pair.second; + const sim::IDE *ide = ide_part.ide; geo::Point_t ide_p(ide->x, ide->y, ide->z); geo::Point_t ide_p_scecorr = WireToTrajectoryPosition(ide_p, w); @@ -819,8 +820,6 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, // Compute the time of each hit for (sbn::TrueHit &h: truehits_v) { - h.time = dprop.ConvertXToTicks(h.p.x, h.plane, h.tpc, h.cryo); - double xdrift = abs(h.p.x - wireReadout->Plane(geo::PlaneID(h.cryo, h.tpc, 0)).GetCenter().X()); h.tdrift = xdrift / dprop.DriftVelocity(); } @@ -1047,7 +1046,7 @@ void sbn::TrackCaloSkimmer::FillTrackTruth(const detinfo::DetectorClocksData &cl const std::vector> &mcparticles, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> id_to_ide_map, + const std::map> id_to_ide_map, const std::map>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, diff --git a/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc b/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc new file mode 100644 index 000000000..54ee42f59 --- /dev/null +++ b/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc @@ -0,0 +1,212 @@ +//////////////////////////////////////////////////////////////////////// +// Class: CCDeltaRadiative +// Plugin Type: filter (art v2_05_00) +// File: CCDeltaRadiative_module.cc +// +// Selects CC events with a photon from a Delta radiative decay. +// Mirrors NCDeltaRadiative_module.cc but requires CCNC == 0 (charged +// current) and checks all four Delta charge states (1114, 2114, 2214, 2224). +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDFilter.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 "art_root_io/TFileService.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileDirectory.h" + +#include + +#include "nusimdata/SimulationBase/MCTruth.h" + +#include "TTree.h" + +class CCDeltaRadiative : public art::EDFilter { + + TTree * ftree; + + int frun; + int fsubrun; + int fevent; + int fnu_pdg; + int fccnc; + int fmode; + int finteraction_type; + int fis_cc_delta_radiative; + int fparent_status_code; + int fparent_pdg; + +public: + explicit CCDeltaRadiative(fhicl::ParameterSet const & p); + + CCDeltaRadiative(CCDeltaRadiative const &) = delete; + CCDeltaRadiative(CCDeltaRadiative &&) = delete; + CCDeltaRadiative & operator = (CCDeltaRadiative const &) = delete; + CCDeltaRadiative & operator = (CCDeltaRadiative &&) = delete; + + void cout_stuff(art::Event & e, bool passed); + void FillTree(art::Event & e, + size_t const mct_index, + size_t const parent_index, + bool const is_cc_delta_radiative); + void Reset(); + bool filter(art::Event & e) override; + +}; + + +CCDeltaRadiative::CCDeltaRadiative(fhicl::ParameterSet const & p) : + art::EDFilter(p), + ftree(nullptr) { + + art::ServiceHandle tfs; + ftree = tfs->make("CCDeltaRadiativeFilter", ""); + + ftree->Branch("run", &frun, "run/I"); + ftree->Branch("subrun", &fsubrun, "subrun/I"); + ftree->Branch("event", &fevent, "event/I"); + ftree->Branch("nu_pdg", &fnu_pdg, "nu_pdg/I"); + ftree->Branch("ccnc", &fccnc, "ccnc/I"); + ftree->Branch("mode", &fmode, "mode/I"); + ftree->Branch("is_cc_delta_radiative", &fis_cc_delta_radiative, "is_cc_delta_radiative/I"); + ftree->Branch("parent_status_code", &fparent_status_code, "parent_status_code/I"); + ftree->Branch("parent_pdg", &fparent_pdg, "parent_pdg/I"); + +} + + +void CCDeltaRadiative::cout_stuff(art::Event & e, bool passed = false) { + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + + std::cout << passed << "\n" + << "==========================\n"; + for(simb::MCTruth const & mct : *ev_mct) { + std::cout << "----------------------------\n"; + for(int i = 0; i < mct.NParticles(); ++i) { + simb::MCParticle const & mcp = mct.GetParticle(i); + std::cout << mcp.TrackId() << " " << mcp.PdgCode() << " " << mcp.Mother() << " " << mcp.StatusCode() << "\n"; + } + } + +} + + +void CCDeltaRadiative::Reset() { + + frun = -1; + fsubrun = -1; + fevent = -1; + fnu_pdg = 0; + fccnc = -1; + fmode = -2; + finteraction_type = -2; + fis_cc_delta_radiative = -1; + fparent_status_code = -1; + fparent_pdg = 0; + +} + + +void CCDeltaRadiative::FillTree(art::Event & e, + size_t const mct_index, + size_t const parent_index, + bool const is_cc_delta_radiative) { + + Reset(); + + frun = e.id().run(); + fsubrun = e.id().subRun(); + fevent = e.id().event(); + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + simb::MCNeutrino const & mcn = ev_mct->at(mct_index).GetNeutrino(); + + fnu_pdg = mcn.Nu().PdgCode(); + fccnc = mcn.CCNC(); + fmode = mcn.Mode(); + finteraction_type = mcn.InteractionType(); + + fis_cc_delta_radiative = is_cc_delta_radiative; + if(parent_index != SIZE_MAX) { + fparent_status_code = ev_mct->at(mct_index).GetParticle(parent_index).StatusCode(); + fparent_pdg = ev_mct->at(mct_index).GetParticle(parent_index).PdgCode(); + } + + ftree->Fill(); + +} + + +// Returns true if |pdg| matches any Delta charge state (Δ-, Δ0, Δ+, Δ++). +static bool isDelta(int pdg) { + int apdg = std::abs(pdg); + return apdg == 1114 || apdg == 2114 || apdg == 2214 || apdg == 2224; +} + + +bool CCDeltaRadiative::filter(art::Event & e) { + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + + for(size_t i = 0; i < ev_mct->size(); ++i) { + + simb::MCTruth const & mct = ev_mct->at(i); + // Require charged-current interaction (CCNC == 0). + if(mct.GetNeutrino().CCNC() != 0) continue; + + std::vector exiting_photon_parents; + for(int j = 0; j < mct.NParticles(); ++j) { + simb::MCParticle const & mcp = mct.GetParticle(j); + if(mcp.TrackId() != j) { + std::cout << "ERROR: " << __LINE__ << " " << __PRETTY_FUNCTION__ << "\nTrackId does not match index\n"; + return false; + } + if(!(mcp.StatusCode() == 1 && mcp.PdgCode() == 22)) continue; + exiting_photon_parents.push_back(mcp.Mother()); + } + + std::vector in_nucleus_photons; + for(size_t const s : exiting_photon_parents) { + simb::MCParticle const & mcp = mct.GetParticle(s); + // Require parent vertex inside TPC. + if(std::abs(mcp.Vx()) > 210 || std::abs(mcp.Vy()) > 210 || mcp.Vz() > 510 || mcp.Vz() < -1) { + std::cout << "OUTSIDE TPC x y z =" << mcp.Vx() << " " << mcp.Vy() << " " << mcp.Vz() << std::endl; + return false; + } + if(isDelta(mcp.PdgCode())) { + if(ftree) FillTree(e, i, s, true); + return true; + } + else if(mcp.PdgCode() == 22) { + in_nucleus_photons.push_back(mcp.Mother()); + } + } + + for(size_t const s : in_nucleus_photons) { + simb::MCParticle const & mcp = mct.GetParticle(s); + if(isDelta(mcp.PdgCode())) { + if(ftree) FillTree(e, i, s, true); + return true; + } + } + + } + + if(ftree) FillTree(e, 0, SIZE_MAX, false); + return false; + +} + + +DEFINE_ART_MODULE(CCDeltaRadiative) diff --git a/sbncode/SinglePhotonAnalysis/CMakeLists.txt b/sbncode/SinglePhotonAnalysis/CMakeLists.txt index 35f24c2bb..4d144d72c 100644 --- a/sbncode/SinglePhotonAnalysis/CMakeLists.txt +++ b/sbncode/SinglePhotonAnalysis/CMakeLists.txt @@ -23,6 +23,13 @@ art_make_library( BASENAME_ONLY ROOT::Core ) +cet_build_plugin(CCDeltaRadiative art::module BASENAME_ONLY LIBRARIES + sbncode_SinglePhotonAnalysis_Libraries + sbncode_SinglePhotonAnalysis_SEAview + art_root_io::TFileService_service + ROOT::Tree + ) + cet_build_plugin(NCDeltaRadiative art::module BASENAME_ONLY LIBRARIES sbncode_SinglePhotonAnalysis_Libraries sbncode_SinglePhotonAnalysis_SEAview diff --git a/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl b/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl new file mode 100644 index 000000000..71178e589 --- /dev/null +++ b/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl @@ -0,0 +1,106 @@ +# Simulates GENIE neutrino interactions from the BNB beam +# forcing one interaction per event, inside the TPC volume +# (with 10 cm padding on each side), +# selecting only CC events with photons coming out from the Delta + + +#include "prodgenie_nu_singleinteraction_tpc_sbnd.fcl" + +#------ this could be a separated file +# +# services +# + +# +# modules +# + + + +process_name: GenieGenFiltered + + +# +# services +# +services: { + + TFileService: { fileName: "hists_genie.root" } + IFDH: {} # required by GENIEGen + @table::sbnd_basic_services # from simulationservices_sbnd.fcl + @table::sbnd_random_services # from simulationservices_sbnd.fcl + FileCatalogMetadata: @local::sbnd_file_catalog_mc # from sam_sbnd.fcl + + # since this is a configuration expected to be run for production, + # we set up message configuration accordingly: + message: @local::sbnd_message_services_prod + +} # services + + +# +# input +# + + +# +# processing +# +physics: +{ + + producers: + { + generator: @local::sbnd_genie_simple + rns: { module_type: "RandomNumberSaver" } + } + + filters: + { + CCDeltaRadFilter: + { + module_type: "CCDeltaRadiative" + } + + } + + + #define the producer and filter modules for this path, order matters, + simulate: [ rns, generator, CCDeltaRadFilter] + + #define the output stream, there could be more than one if using filters + stream1: [ out1 ] + + #trigger_paths is a keyword and contains the paths that modify the art::event, + #ie filters and producers + trigger_paths: [simulate] + + #end_paths is a keyword and contains the paths that do not modify the art::Event, + #ie analyzers and output streams. these all run simultaneously + end_paths: [stream1] + +} # physics + + +# +# output +# +outputs: +{ + out1: + { + module_type: RootOutput + fileName: "prodgenie_bnb_nu_filtered_CCDeltaRad_sbnd_%p-%tc.root" # default file name, can override from command line with -o or --output + SelectEvents: [simulate] + dataTier: "generated" + compressionLevel: 1 + } +} # outputs + + +# +# override +# THIS DOES NOT WORK, CHECK! +physics.producers.generator.TopVolume: "volTPCActive" +physics.producers.generator.BeamName: "booster" +physics.producers.generator.EventGeneratorList: "CCRES" diff --git a/sbncode/WireMod/CMakeLists.txt b/sbncode/WireMod/CMakeLists.txt new file mode 100644 index 000000000..9919ce928 --- /dev/null +++ b/sbncode/WireMod/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(Utility) diff --git a/sbncode/WireMod/Utility/CMakeLists.txt b/sbncode/WireMod/Utility/CMakeLists.txt new file mode 100644 index 000000000..8955de261 --- /dev/null +++ b/sbncode/WireMod/Utility/CMakeLists.txt @@ -0,0 +1,45 @@ +cet_make_library( + LIBRARY_NAME + sbncode_WireMod_Utility + + SOURCE + WireModUtility.cc + + LIBRARIES + ${ART_FRAMEWORK_CORE} + ${MF_MESSAGELOGGER} + ${Boost_SYSTEM_LIBRARY} + ${ROOT_BASIC_LIB_LIST} + ${ROOT_HIST} + art_root_io::TFileService_service + larcorealg::Geometry + larcore::Geometry_Geometry_service + lardataobj::AnalysisBase + lardataobj::RawData + lardataobj::RecoBase + lardataalg::DetectorInfo + lardata::RecoObjects + lardata::Utilities + larreco::RecoAlg + ${LARRECO_LIB} + ${LARDATA_LIB} + ${ART_FRAMEWORK_CORE} + ${ART_FRAMEWORK_PRINCIPAL} + ${ART_FRAMEWORK_SERVICES_REGISTRY} + ${ART_FRAMEWORK_SERVICES_OPTIONAL} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_RANDOMNUMBERGENERATOR_SERVICE} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_TFILESERVICE_SERVICE} + ${MF_MESSAGELOGGER} + ${FHICLCPP} + ${CLHEP} + ${ROOT_GEOM} + ${ROOT_XMLIO} + ${ROOT_GDML} + ${ROOT_BASIC_LIB_LIST} + ${ROOT_HIST} + BASENAME_ONLY +) + +install_headers() +install_source() + diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc new file mode 100644 index 000000000..126044753 --- /dev/null +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -0,0 +1,692 @@ +#include "sbncode/WireMod/Utility/WireModUtility.hh" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" + +//--- CalcROIProperties --- +sys::WireModUtility::ROIProperties_t sys::WireModUtility::CalcROIProperties(recob::Wire const& wire, size_t const& roi_idx) +{ + // get the ROI + recob::Wire::RegionsOfInterest_t::datarange_t const& roi = wire.SignalROI().get_ranges()[roi_idx]; + + // initialize the return value + ROIProperties_t roi_vals; + roi_vals.channel = wire.Channel(); + roi_vals.view = wire.View(); + roi_vals.begin = roi.begin_index(); + roi_vals.end = roi.end_index(); + roi_vals.center = 0; + roi_vals.total_q = 0; + roi_vals.sigma = 0; + + // loop over the roi and find the charge-weighted center and total charge + auto const& roi_data = roi.data(); + for (size_t i_t = 0; i_t < roi_data.size(); ++i_t) + { + roi_vals.center += roi_data[i_t]*(i_t+roi_vals.begin); + roi_vals.total_q += roi_data[i_t]; + } + roi_vals.center = roi_vals.center/roi_vals.total_q; + + // get the width (again charge-weighted) + // if the ROI is only one tick set the cent to the middle of the tick and width to 0.5 + for (size_t i_t = 0; i_t> sys::WireModUtility::GetTargetROIs(sim::SimEnergyDeposit const& shifted_edep, double offset) +{ + // initialize return value + // the pairs are + std::vector> target_roi_vec; + + // if the SimEnergyDeposit isn't inside the TPC then it's not going to match a wire + geo::TPCGeo const* curTPCGeomPtr = geometry->PositionToTPCptr(shifted_edep.MidPoint()); + if (curTPCGeomPtr == nullptr) + return target_roi_vec; + + // iterate over planes + for (auto const& plane : wireReadout->Iterate(curTPCGeomPtr->ID())) { + // check the wire exists in this plane + int wireNumber = int(0.5 + plane.WireCoordinate(shifted_edep.MidPoint())); + if ((wireNumber < 0) || (wireNumber >= (int)plane.Nwires())) + continue; + + // reconstruct the wireID from the position + geo::WireID edep_wireID = plane.NearestWireID(shifted_edep.MidPoint()); + + // if the deposition is inside the range where it would leave a signal on the wire, look for the ROI there + if (planeXInWindow(shifted_edep.X(), plane, *curTPCGeomPtr, offset + tickOffset)) + target_roi_vec.emplace_back(wireReadout->PlaneWireToChannel(edep_wireID), std::round(planeXToTick(shifted_edep.X(), plane, *curTPCGeomPtr, offset + tickOffset))); + } + /* + for (size_t i_p = 0; i_p < curTPCGeomPtr->Nplanes(); ++i_p) + { + // check the wire exists in this plane + int wireNumber = int(0.5 + curTPCGeomPtr->Plane(i_p).WireCoordinate(shifted_edep.MidPoint())); + if ((wireNumber < 0) || (wireNumber >= (int) curTPCGeomPtr->Plane(i_p).Nwires())) + continue; + + // reconstruct the wireID from the position + geo::WireID edep_wireID = curTPCGeomPtr->Plane(i_p).NearestWireID(shifted_edep.MidPoint()); + + // if the deposition is inside the range where it would leave a signal on the wire, look for the ROI there + if (planeXInWindow(shifted_edep.X(), i_p, *curTPCGeomPtr, offset + tickOffset)) + target_roi_vec.emplace_back(geometry->PlaneWireToChannel(edep_wireID), std::round(planeXToTick(shifted_edep.X(), i_p, *curTPCGeomPtr, offset + tickOffset))); + } + */ + + return target_roi_vec; +} + +//--- GetHitTargetROIs --- +std::vector> sys::WireModUtility::GetHitTargetROIs(recob::Hit const& hit) +{ + // initialize return value + // the pairs are + std::vector> target_roi_vec; + + int hit_wire = hit.Channel(); + int hit_tick = int(round(hit.PeakTime())); + + if (hit_tick < tickOffset || hit_tick >= readoutWindowTicks + tickOffset) + return target_roi_vec; + + + target_roi_vec.emplace_back((unsigned int) hit_wire, (unsigned int) hit_tick); + + return target_roi_vec; +} + +//--- FillROIMatchedEdepMap --- +void sys::WireModUtility::FillROIMatchedEdepMap(std::vector const& edepVec, std::vector const& wireVec, double offset) +{ + // clear the map in case it was already set + ROIMatchedEdepMap.clear(); + + // get the channel from each wire and set up a map between them + std::unordered_map wireChannelMap; + for (size_t i_w = 0; i_w < wireVec.size(); ++i_w) + wireChannelMap[wireVec[i_w].Channel()] = i_w; + + // loop over the energy deposits + for (size_t i_e = 0; i_e < edepVec.size(); ++i_e) + { + // get the ROIs + // + std::vector> target_rois = GetTargetROIs(edepVec[i_e], offset); + + // loop over ROI and match the energy deposits with wires + for (auto const& target_roi : target_rois) + { + // if we can't find the wire, skip it + if (wireChannelMap.find(target_roi.first) == wireChannelMap.end()) + continue; + + // get the wire + auto const& target_wire = wireVec.at(wireChannelMap[target_roi.first]); + + // if there are no ticks in in the wire skip it + // likewise if there's nothing in the region of interst + if (not target_wire.SignalROI().is_valid() || + target_wire.SignalROI().empty() || + target_wire.SignalROI().n_ranges() == 0 || + target_wire.SignalROI().size() <= target_roi.second || + target_wire.SignalROI().is_void(target_roi.second) ) + continue; + + // how far into the range is the ROI? + auto range_number = target_wire.SignalROI().find_range_iterator(target_roi.second) - target_wire.SignalROI().begin_range(); + + // popluate the map + ROIMatchedEdepMap[std::make_pair(target_wire.Channel(),range_number)].push_back(i_e); + } + } +} + +//--- FillROIMatchedHitMap --- +void sys::WireModUtility::FillROIMatchedHitMap(std::vector const& hitVec, std::vector const& wireVec) +{ + // clear the map in case it was already set + ROIMatchedHitMap.clear(); + + // get the channel from each wire and set up a map between them + std::unordered_map wireChannelMap; + for (size_t i_w = 0; i_w < wireVec.size(); ++i_w) + wireChannelMap[wireVec[i_w].Channel()] = i_w; + + // loop over hits + for (size_t i_h = 0; i_h < hitVec.size(); ++i_h) + { + // get the ROIs + // + std::vector> target_rois = GetHitTargetROIs(hitVec[i_h]); + + // loop over ROI and match the energy deposits with wires + for (auto const& target_roi : target_rois) + { + // if we can't find the wire, skip it + if (wireChannelMap.find(target_roi.first) == wireChannelMap.end()) + continue; + + auto const& target_wire = wireVec.at(wireChannelMap[target_roi.first]); + + // if there are no ticks in in the wire skip it + // likewise if there's nothing in the region of interst + if (not target_wire.SignalROI().is_valid() || + target_wire.SignalROI().empty() || + target_wire.SignalROI().n_ranges() == 0 || + target_wire.SignalROI().size() <= target_roi.second || + target_wire.SignalROI().is_void(target_roi.second) ) + continue; + + // which range is it? + auto range_number = target_wire.SignalROI().find_range_iterator(target_roi.second) - target_wire.SignalROI().begin_range(); + + // pupluate the map + ROIMatchedHitMap[std::make_pair(target_wire.Channel(),range_number)].push_back(i_h); + } + } +} + +//--- CalcSubROIProperties --- +std::vector sys::WireModUtility::CalcSubROIProperties(sys::WireModUtility::ROIProperties_t const& roi_properties, std::vector const& hitPtrVec) +{ + std::vector subroi_properties_vec; + sys::WireModUtility::SubROIProperties_t subroi_properties; + subroi_properties.channel = roi_properties.channel; + subroi_properties.view = roi_properties.view; + + // if this ROI doesn't contain any hits, define subROI based on ROI properities + // otherwise, define subROIs based on hits + if (hitPtrVec.size() == 0) + { + subroi_properties.key = std::make_pair(roi_properties.key, 0); + subroi_properties.total_q = roi_properties.total_q; + subroi_properties.center = roi_properties.center; + subroi_properties.sigma = roi_properties.sigma; + subroi_properties_vec.push_back(subroi_properties); + } else + { + for (unsigned int i_h=0; i_h < hitPtrVec.size(); ++i_h) + { + auto hit_ptr = hitPtrVec[i_h]; + subroi_properties.key = std::make_pair(roi_properties.key, i_h); + subroi_properties.total_q = hit_ptr->Integral(); + subroi_properties.center = hit_ptr->PeakTime(); + subroi_properties.sigma = hit_ptr->RMS(); + subroi_properties_vec.push_back(subroi_properties); + } + } + + return subroi_properties_vec; +} + +//--- MatchEdepsToSubROIs --- +std::map> sys::WireModUtility::MatchEdepsToSubROIs(std::vector const& subROIPropVec, + std::vector const& edepPtrVec, double offset) +{ + // for each TrackID, which EDeps are associated with it? keys are TrackIDs + std::map> TrackIDMatchedEDepMap; + + // total energy of EDeps matched to the ROI (not strictly necessary, but useful for understanding/development + double total_energy = 0.0; + + // loop over edeps, fill TrackIDMatchedEDepMap and calculate total energy + for (auto const& edep_ptr : edepPtrVec) + { + TrackIDMatchedEDepMap[edep_ptr->TrackID()].push_back(edep_ptr); + total_energy += edep_ptr->E(); + } + + // calculate EDep properties by TrackID + std::map TrackIDMatchedPropertyMap; + for (auto const& track_edeps : TrackIDMatchedEDepMap) + TrackIDMatchedPropertyMap[track_edeps.first] = CalcPropertiesFromEdeps(track_edeps.second, offset); + + // map it all out + std::map> EDepMatchedSubROIMap; // keys are indexes of edepPtrVec, values are vectors of indexes of subROIPropVec + std::map> TrackIDMatchedSubROIMap; // keys are TrackIDs, values are sets of indexes of subROIPropVec + std::map> SubROIMatchedEDepMap; // keys are indexes of subROIPropVec, values are vectors of indexes of edepPtrVec + std::map> SubROIMatchedTrackEnergyMap; // keys are indexes of subROIPropVec, values are maps of TrackIDs to matched energy (in MeV) + + // loop over EDeps + for (unsigned int i_e = 0; i_e < edepPtrVec.size(); ++i_e) + { + // get EDep properties + auto edep_ptr = edepPtrVec[i_e]; + const geo::TPCGeo& curTPCGeom = geometry->PositionToTPC(edep_ptr->MidPoint()); + const auto plane0 = wireReadout->FirstPlane(curTPCGeom.ID()); + double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now + double zeroTick = detPropData.ConvertXToTicks(0, plane0.ID()); + auto edep_tick = ticksPercm * edep_ptr->X() + (zeroTick + offset) + tickOffset; + edep_tick = detPropData.ConvertXToTicks(edep_ptr->X(), plane0.ID()) + offset + tickOffset; + + // loop over subROIs + unsigned int closest_hit = std::numeric_limits::max(); + float min_dist = std::numeric_limits::max(); + for (unsigned int i_h = 0; i_h < subROIPropVec.size(); ++i_h) + { + auto subroi_prop = subROIPropVec[i_h]; + if (edep_tick > subroi_prop.center-subroi_prop.sigma && edep_tick < subroi_prop.center+subroi_prop.sigma) + { + EDepMatchedSubROIMap[i_e].push_back(i_h); + TrackIDMatchedSubROIMap[edep_ptr->TrackID()].emplace(i_h); + } + float hit_dist = std::abs(edep_tick - subroi_prop.center) / subroi_prop.sigma; + if (hit_dist < min_dist) + { + closest_hit = i_h; + min_dist = hit_dist; + } + } + + // if EDep is less than 2.5 units away from closest subROI, assign it to that subROI + if (min_dist < 5) // try 5 for testing purposes + { + auto i_h = closest_hit; + SubROIMatchedEDepMap[i_h].push_back(i_e); + SubROIMatchedTrackEnergyMap[i_h][edep_ptr->TrackID()] += edep_ptr->E(); + } + } + + // convert to desired format (possibly a better way to do this...?) + std::map> ReturnMap; + for (auto it_h = SubROIMatchedEDepMap.begin(); it_h != SubROIMatchedEDepMap.end(); ++it_h) + { + auto key = subROIPropVec[it_h->first].key; + for (auto const& i_e : it_h->second) + { + ReturnMap[key].push_back(edepPtrVec[i_e]); + } + } + + return ReturnMap; +} + +//--- CalcPropertiesFromEdeps --- +sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEdeps(std::vector const& edepPtrVec, double offset) +{ + //split the edeps by TrackID + std::map< int, std::vector > edepptrs_by_trkid; + std::map< int, double > energy_per_trkid; + for(auto const& edep_ptr : edepPtrVec) + { + edepptrs_by_trkid[edep_ptr->TrackID()].push_back(edep_ptr); + energy_per_trkid[edep_ptr->TrackID()]+=edep_ptr->E(); + } + + int trkid_max = std::numeric_limits::min(); + double energy_max = std::numeric_limits::min(); + for(auto const& e_p_id : energy_per_trkid) + { + if(e_p_id.second > energy_max) + { + trkid_max = e_p_id.first; + energy_max = e_p_id.second; + } + } + + auto edepPtrVecMaxE = edepptrs_by_trkid[trkid_max]; + + //first, let's loop over all edeps and get an average weight scale... + sys::WireModUtility::TruthProperties_t edep_props; + double total_energy_all = 0.0; + sys::WireModUtility::ScaleValues_t scales_e_weighted[3]; + for(size_t i_p = 0; i_p < 3; ++i_p) + { + scales_e_weighted[i_p].r_Q = 0.0; + scales_e_weighted[i_p].r_sigma = 0.0; + } + + for(auto const edep_ptr : edepPtrVec) + { + if (edep_ptr->StepLength() == 0) + continue; + + edep_props.x = edep_ptr->X(); + edep_props.y = edep_ptr->Y(); + edep_props.z = edep_ptr->Z(); + + edep_props.dxdr = (edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); + edep_props.dydr = (edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); + edep_props.dzdr = (edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); + + edep_props.dedr = edep_ptr->E() / edep_ptr->StepLength(); + + total_energy_all += edep_ptr->E(); + + const geo::TPCGeo& curTPCGeom = geometry->PositionToTPC(edep_ptr->MidPoint()); + for (auto const& plane : wireReadout->Iterate(curTPCGeom.ID())) { + int i_p = plane.ID().Plane; + auto scales = GetViewScaleValues(edep_props, plane.View()); + scales_e_weighted[i_p].r_Q += edep_ptr->E()*scales.r_Q; + scales_e_weighted[i_p].r_sigma += edep_ptr->E()*scales.r_sigma; + } + } + + for(size_t i_p = 0; i_p < 3; ++i_p) + { + if (total_energy_all > 0) + { + scales_e_weighted[i_p].r_Q = scales_e_weighted[i_p].r_Q / total_energy_all; + scales_e_weighted[i_p].r_sigma = scales_e_weighted[i_p].r_sigma / total_energy_all; + } + if (scales_e_weighted[i_p].r_Q == 0) + scales_e_weighted[i_p].r_Q = 1; + if (scales_e_weighted[i_p].r_sigma == 0) + scales_e_weighted[i_p].r_sigma = 1; + } + + TruthProperties_t edep_col_properties; + + //copy in the scales that were calculated before + for(size_t i_p = 0; i_p < 3; ++i_p) + { + edep_col_properties.scales_avg[i_p].r_Q = scales_e_weighted[i_p].r_Q; + edep_col_properties.scales_avg[i_p].r_sigma = scales_e_weighted[i_p].r_sigma; + } + + // calculations happen here + edep_col_properties.x = 0.; + edep_col_properties.x_rms = 0.; + edep_col_properties.x_rms_noWeight = 0.; + edep_col_properties.x_min = std::numeric_limits::max(); + edep_col_properties.x_max = std::numeric_limits::min(); + + edep_col_properties.y = 0.; + edep_col_properties.z = 0.; + edep_col_properties.dxdr = 0.; + edep_col_properties.dydr = 0.; + edep_col_properties.dzdr = 0.; + + edep_col_properties.dedr = 0.; + + double total_energy = 0.0; + for (auto const& edep_ptr : edepPtrVecMaxE) + { + edep_col_properties.x += edep_ptr->X()*edep_ptr->E(); + edep_col_properties.x_min = (edep_ptr->X() < edep_col_properties.x_min) ? edep_ptr->X() : edep_col_properties.x_min; + edep_col_properties.x_max = (edep_ptr->X() > edep_col_properties.x_max) ? edep_ptr->X() : edep_col_properties.x_max; + total_energy += edep_ptr->E(); + edep_col_properties.y += edep_ptr->Y()*edep_ptr->E(); + edep_col_properties.z += edep_ptr->Z()*edep_ptr->E(); + + if (edep_ptr->StepLength() == 0) + continue; + + edep_col_properties.dxdr += edep_ptr->E()*(edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); + edep_col_properties.dydr += edep_ptr->E()*(edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); + edep_col_properties.dzdr += edep_ptr->E()*(edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); + + edep_col_properties.dedr += edep_ptr->E()*edep_ptr->E() / edep_ptr->StepLength(); + } + + if (total_energy > 0) + { + edep_col_properties.x = edep_col_properties.x / total_energy; + edep_col_properties.y = edep_col_properties.y / total_energy; + edep_col_properties.z = edep_col_properties.z / total_energy; + edep_col_properties.dxdr = edep_col_properties.dxdr / total_energy; + edep_col_properties.dydr = edep_col_properties.dydr / total_energy; + edep_col_properties.dzdr = edep_col_properties.dzdr / total_energy; + + edep_col_properties.dedr = edep_col_properties.dedr / total_energy; + } + + for (auto const& edep_ptr : edepPtrVecMaxE) + { + edep_col_properties.x_rms += (edep_ptr->X()-edep_col_properties.x)*(edep_ptr->X()-edep_col_properties.x)*edep_ptr->E(); + edep_col_properties.x_rms_noWeight += (edep_ptr->X()-edep_col_properties.x)*(edep_ptr->X()-edep_col_properties.x); + } + edep_col_properties.x_rms_noWeight = std::sqrt(edep_col_properties.x_rms_noWeight); + + if (total_energy > 0) + edep_col_properties.x_rms = std::sqrt(edep_col_properties.x_rms/total_energy); + + const geo::TPCGeo& tpcGeom = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}); + const auto plane0 = wireReadout->FirstPlane(tpcGeom.ID()); + double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now + edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_rms = ticksPercm*edep_col_properties.x_rms; + edep_col_properties.tick_rms_noWeight = ticksPercm*edep_col_properties.x_rms_noWeight; + edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, plane0.ID()) + offset + tickOffset; + edep_col_properties.total_energy = total_energy; + + + return edep_col_properties; +} + +//--- GetScaleValues --- +sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, sys::WireModUtility::ROIProperties_t const& roi_vals) +{ + sys::WireModUtility::ScaleValues_t scales; + sys::WireModUtility::ScaleValues_t channelScales = GetChannelScaleValues(truth_props, roi_vals.channel); + sys::WireModUtility::ScaleValues_t viewScales = GetViewScaleValues(truth_props, roi_vals.view); + scales.r_Q = channelScales.r_Q * viewScales.r_Q; + scales.r_sigma = channelScales.r_sigma * viewScales.r_sigma; + return scales; +} + +//--- GetChannelScaleValues --- +sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetChannelScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, raw::ChannelID_t const& channel) +{ + // initialize return + sys::WireModUtility::ScaleValues_t scales; + scales.r_Q = 1.0; + scales.r_sigma = 1.0; + + // try to get geo + // if not in a TPC return default values + double const truth_coords[3] = {truth_props.x, truth_props.y, truth_props.z}; + geo::TPCGeo const* curTPCGeomPtr = geometry->PositionToTPCptr(geo::vect::makePointFromCoords(truth_coords)); + if (curTPCGeomPtr == nullptr) + return scales; + + if (applyChannelScale) + { + if (spline_Charge_Channel == nullptr || + spline_Sigma_Channel == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply channel scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= spline_Charge_Channel->Eval(channel); + scales.r_sigma *= spline_Sigma_Channel ->Eval(channel); + } + + return scales; +} + +//--- GetViewScaleValues --- +sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, geo::View_t const& view) +{ + // initialize return + sys::WireModUtility::ScaleValues_t scales; + scales.r_Q = 1.0; + scales.r_sigma = 1.0; + + double temp_scale=1.0; + + // try to get geo + // if not in a TPC return default values + double const truth_coords[3] = {truth_props.x, truth_props.y, truth_props.z}; + geo::TPCGeo const* curTPCGeomPtr = geometry->PositionToTPCptr(geo::vect::makePointFromCoords(truth_coords)); + if (curTPCGeomPtr == nullptr) + return scales; + + // get the plane number by the view + auto const& plane_obj = wireReadout->Plane(curTPCGeomPtr->ID(), view); + unsigned int plane = plane_obj.ID().Plane; + + if (applyXScale) + { + if (splines_Charge_X[plane] == nullptr || + splines_Sigma_X [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply X scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= splines_Charge_X[plane]->Eval(truth_props.x); + scales.r_sigma *= splines_Sigma_X [plane]->Eval(truth_props.x); + } + + if (applyYZScale) + { + if (graph2Ds_Charge_YZ[plane] == nullptr || + graph2Ds_Sigma_YZ [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply YZ scale factor, but could not find graphs. Check that you have set those in the utility."; + temp_scale = graph2Ds_Charge_YZ[plane]->Interpolate(truth_props.z, truth_props.y); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_YZ [plane]->Interpolate(truth_props.z, truth_props.y); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXZAngleScale) + { + if (splines_Charge_XZAngle[plane] == nullptr || + splines_Sigma_XZAngle [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply XZ-angle scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= splines_Charge_XZAngle[plane]->Eval(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + scales.r_sigma *= splines_Sigma_XZAngle [plane]->Eval(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + } + if (applyYZAngleScale) + { + if (splines_Charge_YZAngle[plane] == nullptr || + splines_Sigma_YZAngle [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply YZ-angle scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= splines_Charge_YZAngle[plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + scales.r_sigma *= splines_Sigma_YZAngle [plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + } + + if (applydEdXScale) + { + if (splines_Charge_dEdX[plane] == nullptr || + splines_Sigma_dEdX [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply dEdX scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= splines_Charge_dEdX[plane]->Eval(truth_props.dedr); + scales.r_sigma *= splines_Sigma_dEdX [plane]->Eval(truth_props.dedr); + } + + if (applyXXZAngleScale) + { + if (graph2Ds_Charge_XXZAngle[plane] == nullptr || + graph2Ds_Sigma_XXZAngle [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply XXZAngle scale factor, but could not find graphs. Check that you have set those in the utility."; + temp_scale = graph2Ds_Charge_XXZAngle[plane]->Interpolate(truth_props.x, ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_XXZAngle [plane]->Interpolate(truth_props.x, ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXdQdXScale) + { + if (graph2Ds_Charge_XdQdX[plane] == nullptr || + graph2Ds_Sigma_XdQdX [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply XdQdX scale factor, but could not find graphs. Check that you have set those in the utility."; + temp_scale = graph2Ds_Charge_XdQdX[plane]->Interpolate(truth_props.x, truth_props.dqdr); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_XdQdX [plane]->Interpolate(truth_props.x, truth_props.dqdr); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXZAngledQdXScale) + { + if (graph2Ds_Charge_XZAngledQdX[plane] == nullptr || + graph2Ds_Sigma_XZAngledQdX [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply XZAngledQdX scale factor, but could not find graphs. Check that you have set those in the utility."; + temp_scale = graph2Ds_Charge_XZAngledQdX[plane]->Interpolate(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), truth_props.dqdr); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_XZAngledQdX [plane]->Interpolate(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), truth_props.dqdr); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + return scales; +} + +//--- ModifyROI --- +void sys::WireModUtility::ModifyROI(std::vector & roi_data, + sys::WireModUtility::ROIProperties_t const& roi_prop, + std::vector const& subROIPropVec, + std::map const& subROIScaleMap) +{ + // do you want a bunch of messages? + bool verbose = false; + + // initialize some values + double q_orig = 0.0; + double q_mod = 0.0; + double scale_ratio = 1.0; + + // loop over the ticks + for(size_t i_t = 0; i_t < roi_data.size(); ++i_t) + { + // reset your values + q_orig = 0.0; + q_mod = 0.0; + scale_ratio = 1.0; + + // loop over the subs + for (auto const& subroi_prop : subROIPropVec) + { + // get your scale vals + auto scale_vals = subROIScaleMap.find(subroi_prop.key)->second; + + q_orig += gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); + q_mod += gausFunc(i_t + roi_prop.begin, subroi_prop.center, scale_vals.r_sigma * subroi_prop.sigma, scale_vals.r_Q * subroi_prop.total_q); + + if (verbose) + std::cout << " Incrementing q_orig by gausFunc(" << i_t + roi_prop.begin << ", " << subroi_prop.center << ", " << subroi_prop.sigma << ", " << subroi_prop.total_q << ")" << '\n' + << " Incrementing q_mod by gausFunc(" << i_t + roi_prop.begin << ", " << subroi_prop.center << ", " << scale_vals.r_sigma * subroi_prop.sigma << ", " << scale_vals.r_Q * subroi_prop.total_q << ")" << std::endl; + } + + // do some sanity checks + if (isnan(q_orig)) + { + if (verbose) + std::cout << "WARNING: obtained q_orig = NaN... setting scale to 1" << std::endl; + scale_ratio = 1.0; + } else if (isnan(q_mod)) { + if (verbose) + std::cout << "WARNING: obtained q_mod = NaN... setting scale to 0" << std::endl; + scale_ratio = 0.0; + } else { + scale_ratio = q_mod / q_orig; + } + if(isnan(scale_ratio) || isinf(scale_ratio)) + { + if (verbose) + std::cout << "WARNING: obtained scale_ratio = " << q_mod << " / " << q_orig << " = NAN/Inf... setting to 1" << std::endl; + scale_ratio = 1.0; + } + + roi_data[i_t] = scale_ratio * roi_data[i_t]; + if (verbose) + std::cout << "\t tick " << i_t << ":" + << " data=" << roi_data[i_t] + << ", q_orig=" << q_orig + << ", q_mod=" << q_mod + << ", ratio=" << scale_ratio << std::endl; + } + + // we're done now + return; +} diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh new file mode 100644 index 000000000..1ae74c152 --- /dev/null +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -0,0 +1,243 @@ +#ifndef sbncode_WireMod_Utility_WireModUtility_hh +#define sbncode_WireMod_Utility_WireModUtility_hh + +//std includes +#include + +//ROOT includes +#include "TFile.h" +#include "TSpline.h" +#include "TGraph2D.h" +#include "TNtuple.h" + +//Framework includes +#include "larcorealg/Geometry/GeometryCore.h" +#include "larcorealg/Geometry/WireReadoutGeom.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" +#include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/TrackHitMeta.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/Simulation/SimEnergyDeposit.h" +#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" +#include "nusimdata/SimulationBase/MCParticle.h" +#include "nusimdata/SimulationBase/MCTruth.h" + +#include +#include + +namespace sys { + class WireModUtility{ + // for now make everything public, though it's probably a good idea to think about what doesn't need to be + public: + // detector constants, should be set by geometry service + // the notes at the end refer to their old names in the MircoBooNE code which preceded this + // TODO: how best to initialize the splines/graphs? + const geo::GeometryCore* geometry; // save the TPC geometry + const geo::WireReadoutGeom* wireReadout; // new for LarSoft v10 + const detinfo::DetectorPropertiesData& detPropData; // save the detector property data + bool applyChannelScale; // do we scale with channel? + bool applyXScale; // do we scale with X? + bool applyYZScale; // do we scale with YZ? + bool applyXZAngleScale; // do we scale with XZ angle? + bool applyYZAngleScale; // do we scale with YZ angle? + bool applydEdXScale; // do we scale with dEdx? + bool applyXXZAngleScale; // do we scale with X vs XZ angle? + bool applyXdQdXScale; // do we scale with X vs dQ/dX? + bool applyXZAngledQdXScale; // do we scale with XZ angle vs dQ/dX? + double readoutWindowTicks; // how many ticks are in the readout window? + double tickOffset; // do we want an offset in the ticks? + + TSpline3* spline_Charge_Channel; // the spline for the charge correction by channel + TSpline3* spline_Sigma_Channel; // the spline for the width correction by channel + std::vector splines_Charge_X; // the splines for the charge correction in X + std::vector splines_Sigma_X; // the splines for the width correction in X + std::vector splines_Charge_XZAngle; // the splines for the charge correction in XZ angle + std::vector splines_Sigma_XZAngle; // the splines for the width correction in XZ angle + std::vector splines_Charge_YZAngle; // the splines for the charge correction in YZ angle + std::vector splines_Sigma_YZAngle; // the splines for the width correction in YZ angle + std::vector splines_Charge_dEdX; // the splines for the charge correction in dEdX + std::vector splines_Sigma_dEdX; // the splines for the width correction in dEdX + std::vector graph2Ds_Charge_YZ; // the graphs for the charge correction in YZ + std::vector graph2Ds_Sigma_YZ; // the graphs for the width correction in YZ + + std::vector graph2Ds_Charge_XXZAngle; // the graphs for the charge correction in X vs XZ angle + std::vector graph2Ds_Sigma_XXZAngle; // the graphs for the width correction in X vs XZ angle + std::vector graph2Ds_Charge_XdQdX; // the graphs for charge correction in X vs dQ/dX + std::vector graph2Ds_Sigma_XdQdX; // the graphs for width correction in X vs dQ/dX + std::vector graph2Ds_Charge_XZAngledQdX; // the graphs for charge correction in XZ angle vs dQ/dX + std::vector graph2Ds_Sigma_XZAngledQdX; // the graphs for width correction in XZ angle vs dQ/dX + + // lets try making a constructor here + // assume we can get a geometry service, a detector clcok, and a detector properties + // pass the CryoStat and TPC IDs because it's IDs all the way down + // set some optional args fpr the booleans, the readout window, and the offset + WireModUtility(const geo::GeometryCore* geom, + const geo::WireReadoutGeom* wireRead, + const detinfo::DetectorPropertiesData& detProp, + const bool& arg_ApplyChannelScale = false, + const bool& arg_ApplyXScale = true, + const bool& arg_ApplyYZScale = true, + const bool& arg_ApplyXZAngleScale = true, + const bool& arg_ApplyYZAngleScale = true, + const bool& arg_ApplydEdXScale = true, + const bool& arg_ApplyXXZAngleScale = false, + const bool& arg_ApplyXdQdXScale = false, + const bool& arg_ApplyXZAngledQdXScale = false, + const double& arg_TickOffset = 0) + : geometry(geom), + wireReadout(wireRead), + detPropData(detProp), + applyChannelScale(arg_ApplyChannelScale), + applyXScale(arg_ApplyXScale), + applyYZScale(arg_ApplyYZScale), + applyXZAngleScale(arg_ApplyXZAngleScale), + applyYZAngleScale(arg_ApplyYZAngleScale), + applydEdXScale(arg_ApplydEdXScale), + applyXXZAngleScale(arg_ApplyXXZAngleScale), + applyXdQdXScale(arg_ApplyXdQdXScale), + applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), + readoutWindowTicks(detProp.ReadOutWindowSize()), // the default A2795 (ICARUS TPC readout board) readout window is 4096 samples + tickOffset(arg_TickOffset) // tick offset is for MC truth, default to zero and set only as necessary + { + } + + // typedefs + typedef std::pair ROI_Key_t; + typedef std::pair SubROI_Key_t; + + typedef struct ROIProperties + { + ROI_Key_t key; + raw::ChannelID_t channel; + geo::View_t view; + float begin; + float end; + float total_q; + float center; //charge weighted center of ROI + float sigma; //charge weighted RMS of ROI + } ROIProperties_t; + + typedef struct SubROIProperties + { + SubROI_Key_t key; + raw::ChannelID_t channel; + geo::View_t view; + float total_q; + float center; + float sigma; + } SubROIProperties_t; + + typedef struct ScaleValues + { + double r_Q; + double r_sigma; + } ScaleValues_t; + + typedef struct TruthProperties + { + float x; + float x_rms; + float x_rms_noWeight; + float tick; + float tick_rms; + float tick_rms_noWeight; + float total_energy; + float x_min; + float x_max; + float tick_min; + float tick_max; + float y; + float z; + double dxdr; + double dydr; + double dzdr; + double dqdr; + double dedr; + ScaleValues_t scales_avg[3]; + } TruthProperties_t; + + // internal containers + std::map< ROI_Key_t,std::vector > ROIMatchedEdepMap; + std::map< ROI_Key_t,std::vector > ROIMatchedHitMap; + + // some useful functions + // geometries + // TODO is this the most efficient for new v10 iterators? + double planeXToTick(double xPos, const geo::PlaneGeo& plane, const geo::TPCGeo& tpcGeom, double offset = 0) { + return detPropData.ConvertXToTicks(xPos, plane.ID()) + offset; + } + + bool planeXInWindow(double xPos, const geo::PlaneGeo& plane, const geo::TPCGeo& tpcGeom, double offset = 0) + { + double tick = planeXToTick(xPos, plane, tpcGeom, offset); + return (tick > 0 && tick <= detPropData.ReadOutWindowSize()); + } + + // for this function: in the future if we want to use non-gaussian functions make this take a vector of parameters + // the another wiremod utility could overwrite the ``fitFunc'' with some non-standard function + // would require a fair bit of remodling (ie q and sigma would need to be replace with, eg, funcVar[0] and funcVar[1] and probs a bunch of loops) + // so lets worry about that later + double gausFunc(double t, double mean, double sigma, double a = 1.0) + { + return (a / (sigma * std::sqrt(2 * util::pi()))) * std::exp(-0.5 * std::pow((t - mean)/sigma, 2)); + } + + double FoldAngle(double theta) + { + return (std::abs(theta) > 0.5 * util::pi()) ? util::pi() - std::abs(theta) : std::abs(theta); + } + + double ThetaXZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) + { + double planeAngleRad = planeAngle * (util::pi() / 180.0); + double sinPlaneAngle = std::sin(planeAngleRad); + double cosPlaneAngle = std::cos(planeAngleRad); + + //double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; // don't need to rotate Y for this angle + double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; + + double theta = std::atan2(dxdr, dzdrPlaneRel); + return FoldAngle(theta); + } + + double ThetaYZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) + { + double planeAngleRad = planeAngle * (util::pi() / 180.0); + double sinPlaneAngle = std::sin(planeAngleRad); + double cosPlaneAngle = std::cos(planeAngleRad); + + double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; + double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; + + double theta = std::atan2(dydrPlaneRel, dzdrPlaneRel); + return FoldAngle(theta); + } + + // theste are set in the .cc file + ROIProperties_t CalcROIProperties(recob::Wire const&, size_t const&); + + std::vector> GetTargetROIs(sim::SimEnergyDeposit const&, double offset); + std::vector> GetHitTargetROIs(recob::Hit const&); + + void FillROIMatchedEdepMap(std::vector const&, std::vector const&, double offset); + void FillROIMatchedHitMap(std::vector const&, std::vector const&); + + std::vector CalcSubROIProperties(ROIProperties_t const&, std::vector const&); + + std::map> MatchEdepsToSubROIs(std::vector const&, std::vector const&, double offset); + + TruthProperties_t CalcPropertiesFromEdeps(std::vector const&, double offset); + + ScaleValues_t GetScaleValues(TruthProperties_t const&, ROIProperties_t const&); + ScaleValues_t GetChannelScaleValues(TruthProperties_t const&, raw::ChannelID_t const&); + ScaleValues_t GetViewScaleValues(TruthProperties_t const&, geo::View_t const&); + + void ModifyROI(std::vector &, + ROIProperties_t const &, + std::vector const&, + std::map const&); + }; // end class +} // end namespace + +#endif // sbncode_WireMod_Utility_WireModUtility_hh