diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index f3c381ebd..770f546af 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -64,7 +64,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" @@ -92,6 +92,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" @@ -1386,15 +1387,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); } @@ -1984,8 +1985,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 = @@ -2239,10 +2240,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 76c3903d4..2326f8dee 100644 --- a/sbncode/CAFMaker/CMakeLists.txt +++ b/sbncode/CAFMaker/CMakeLists.txt @@ -29,6 +29,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 @@ -65,6 +66,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 d81cb5890..5089f8631 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -794,10 +794,30 @@ 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 @@ -832,7 +852,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; @@ -844,6 +865,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; + } } } @@ -897,9 +931,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) { @@ -912,7 +949,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 6e8ecf292..38c5d5f6b 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" @@ -204,9 +207,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, @@ -222,9 +228,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 491c813c0..aeeb0646d 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; @@ -778,7 +777,7 @@ namespace caf { exit_point = particle.NumberTrajectoryPoints() - 1; } 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 @@ -799,8 +798,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(); @@ -871,37 +870,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 @@ -1253,126 +1221,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); - - 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_ 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/Calibration/CMakeLists.txt b/sbncode/Calibration/CMakeLists.txt index 98059edc5..759c5c9f7 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 f9fb1930e..1096c056d 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" @@ -66,6 +66,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}; @@ -112,7 +115,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; @@ -161,7 +163,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 5fc33ab4d..ff73432be 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -12,13 +12,10 @@ #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" @@ -288,15 +285,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(); } @@ -545,14 +542,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; @@ -569,9 +566,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*/; @@ -622,7 +619,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; @@ -690,7 +687,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 @@ -711,8 +708,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(); @@ -720,10 +717,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; @@ -739,6 +737,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); @@ -754,10 +754,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); @@ -782,8 +782,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(); } @@ -1010,7 +1008,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,