Merge branch 'develop' into feature/kplows_flat_genie_record

Author: kjplows

CMakeLists.txt

@@ -16,7 +16,7 @@
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
 find_package(cetmodules 3.20.00 REQUIRED)
-project(sbncode VERSION 09.93.01 LANGUAGES CXX)
+project(sbncode VERSION 10.03.01 LANGUAGES CXX)
 
 message(STATUS "\n\n  ==========================   ${PROJECT_NAME}   ==========================")
 

sbncode/CAFMaker/CAFMaker_module.cc

@@ -87,6 +87,8 @@
 #include "messagefacility/MessageLogger/MessageLogger.h"
 
 // LArSoft includes
+#include "larcore/Geometry/WireReadout.h"
+#include "larcore/Geometry/Geometry.h"
 #include "lardataobj/RecoBase/PFParticle.h"
 #include "lardataobj/RecoBase/Slice.h"
 #include "lardataobj/RecoBase/Track.h"
@@ -99,8 +101,6 @@
 #include "nusimdata/SimulationBase/MCNeutrino.h"
 #include "nusimdata/SimulationBase/GTruth.h"
 
-#include "fhiclcpp/ParameterSetRegistry.h"
-
 #include "sbnobj/Common/EventGen/MeVPrtl/MeVPrtlTruth.h"
 #include "sbnobj/Common/Reco/RangeP.h"
 #include "sbnobj/Common/SBNEventWeight/EventWeightMap.h"
@@ -1351,7 +1351,9 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
   auto const dprop =
     art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
-  const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
+  const geo::GeometryCore* geom = lar::providerFrom<geo::Geometry>();
+  const geo::WireReadoutGeom &wireReadout =
+    art::ServiceHandle<geo::WireReadout>()->Get();
 
   auto const *sce = lar::providerFrom<spacecharge::SpaceChargeService>();
 
@@ -1379,7 +1381,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   if ( !isRealData ) {
     art::ServiceHandle<cheat::BackTrackerService> bt_serv;
 
-    id_to_ide_map = PrepSimChannels(simchannels, *geometry);
+    id_to_ide_map = PrepSimChannels(simchannels, wireReadout);
     id_to_truehit_map = PrepTrueHits(hits, clock_data, *bt_serv);
     id_to_hit_energy_map = SetupIDHitEnergyMap(hits, clock_data, *bt_serv);
   }
@@ -1638,22 +1640,17 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       }
     }
 
-  // Get all of the CRTPMT Matches ..
+  // Get all of the CRTPMT Matches
   std::vector<caf::SRCRTPMTMatch> srcrtpmtmatches;
-  std::cout << "srcrtpmtmatches.size = " << srcrtpmtmatches.size() << "\n";
   art::Handle<std::vector<sbn::crt::CRTPMTMatching>> crtpmtmatch_handle;
   GetByLabelStrict(evt, fParams.CRTPMTLabel(), crtpmtmatch_handle);
   if(crtpmtmatch_handle.isValid()){
-    std::cout << "valid handle! label: " << fParams.CRTPMTLabel() << "\n";
     const std::vector<sbn::crt::CRTPMTMatching> &crtpmtmatches = *crtpmtmatch_handle;
     for (unsigned i = 0; i < crtpmtmatches.size(); i++) {
       srcrtpmtmatches.emplace_back();
-      FillCRTPMTMatch(crtpmtmatches[i],srcrtpmtmatches.back());
+      FillCRTPMTMatch(crtpmtmatches[i], srcrtpmtmatches.back());
     }
   }
-  else{
-    std::cout << "crtpmtmatch_handle.isNOTValid!\n";
-  }
 
   // Get all of the OpFlashes
   std::vector<caf::SROpFlash> srflashes;
@@ -2136,7 +2133,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
         FillTrackRangeP(*thisTrack[0], rangePs, trk);
 
         if (fmChi2PID.isValid()) {
-           FillTrackChi2PID(fmChi2PID.at(iPart), lar::providerFrom<geo::Geometry>(), trk);
+           FillTrackChi2PID(fmChi2PID.at(iPart), trk);
         }
         if (fmScatterClosestApproach.isValid() && fmScatterClosestApproach.at(iPart).size()==1) {
            FillTrackScatterClosestApproach(fmScatterClosestApproach.at(iPart).front(), trk);
@@ -2151,7 +2148,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
           FillTrackCalo(fmCalo.at(iPart), fmTrackHit.at(iPart),
               (fParams.FillHitsNeutrinoSlices() && NeutrinoSlice) || fParams.FillHitsAllSlices(),
               fParams.TrackHitFillRRStartCut(), fParams.TrackHitFillRREndCut(),
-              lar::providerFrom<geo::Geometry>(), dprop, trk);
+              dprop, trk);
         }
         if (fmCRTHitMatch.isValid() && fDet == kICARUS) {
           art::FindManyP<sbn::crt::CRTHit> CRTT02Hit = FindManyPStrict<sbn::crt::CRTHit>
@@ -2189,7 +2186,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
             FillTrackTruth(fmTrackHit.at(iPart), id_to_hit_energy_map, true_particles, clock_data, trk);
             // Hit truth information corresponding to Calo-Points
             // Assumes truth matching and calo-points are filled
-            if (mc_particles.isValid() && fParams.FillTrackCaloTruth()) FillTrackCaloTruth(id_to_ide_map, *mc_particles, geometry, clock_data, sce, trk);
+            if (mc_particles.isValid() && fParams.FillTrackCaloTruth()) FillTrackCaloTruth(id_to_ide_map, *mc_particles, *geom, wireReadout, clock_data, sce, trk);
           }
         }
       } // thisTrack exists
@@ -2198,7 +2195,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
         assert(thisShower.size() == 1);
 
         SRShower& shw = pfp.shw;
-        FillShowerVars(*thisShower[0], vertex, fmShowerHit.at(iPart), lar::providerFrom<geo::Geometry>(), producer, shw);
+        FillShowerVars(*thisShower[0], vertex, fmShowerHit.at(iPart), wireReadout, producer, shw);
 
         // We may have many residuals per shower depending on how many showers ar in the slice
         if (fmShowerRazzle.isValid() && fmShowerRazzle.at(iPart).size()==1) {

sbncode/CAFMaker/FillReco.cxx

@@ -4,6 +4,9 @@
 // \author  $Author: psihas@fnal.gov
 //////////////////////////////////////////////////////////////////////
 
+#include "larcorealg/Geometry/GeometryCore.h"
+#include "larcorealg/Geometry/WireReadoutGeom.h"
+
 #include "FillReco.h"
 #include "RecoUtils/RecoUtils.h"
 
@@ -154,17 +157,25 @@ namespace caf
     srmatch.flashPosition = SRVector3D (match.flashPosition.X(), match.flashPosition.Y(), match.flashPosition.Z());
     srmatch.flashYWidth = match.flashYWidth;
     srmatch.flashZWidth = match.flashZWidth;
-    srmatch.flashClassification = static_cast<int>(match.flashClassification);
+    unsigned int topen = 0, topex = 0, sideen = 0, sideex = 0;
     for(const auto& matchedCRTHit : match.matchedCRTHits){
-      std::cout << "CRTPMTTimeDiff = "<< matchedCRTHit.PMTTimeDiff << "\n";
       caf::SRMatchedCRT matchedCRT;
       matchedCRT.PMTTimeDiff = matchedCRTHit.PMTTimeDiff; 
       matchedCRT.time = matchedCRTHit.time;
       matchedCRT.sys = matchedCRTHit.sys;
       matchedCRT.region = matchedCRTHit.region;
       matchedCRT.position = SRVector3D(matchedCRTHit.position.X(), matchedCRTHit.position.Y(), matchedCRTHit.position.Z());
+      if(matchedCRTHit.PMTTimeDiff < 0){
+        if(matchedCRTHit.sys == 0) topen++;
+        else if(matchedCRTHit.sys == 1) sideen++;
+      }
+      else if(matchedCRTHit.PMTTimeDiff >= 0){
+        if(matchedCRTHit.sys == 0) topex++;
+        else if(matchedCRTHit.sys == 1) sideex++;
+      }
       srmatch.matchedCRTHits.push_back(matchedCRT);
     }
+    srmatch.flashClassification = static_cast<int>(sbn::crt::assignFlashClassification(topen, topex, sideen, sideex, 0, 0));
   }
 
 
@@ -228,7 +239,7 @@ namespace caf
   void FillShowerVars(const recob::Shower& shower,
                       const recob::Vertex* vertex,
                       const std::vector<art::Ptr<recob::Hit>> &hits,
-                      const geo::GeometryCore *geom,
+                      const geo::WireReadoutGeom& wireReadout,
                       unsigned producer,
                       caf::SRShower &srshower,
                       bool allowEmpty)
@@ -281,9 +292,9 @@ namespace caf
     for(int p = 0; p < 3; ++p) srshower.plane[p].nHits = 0;
     for (auto const& hit:hits) ++srshower.plane[hit->WireID().Plane].nHits;
 
-    for (geo::PlaneGeo const& plane: geom->Iterate<geo::PlaneGeo>()) {
+    for (geo::PlaneGeo const& plane: wireReadout.Iterate<geo::PlaneGeo>()) {
 
-      const double angleToVert(geom->WireAngleToVertical(plane.View(), plane.ID()) - 0.5*M_PI);
+      const double angleToVert(wireReadout.WireAngleToVertical(plane.View(), plane.ID()) - 0.5*M_PI);
       const double cosgamma(std::abs(std::sin(angleToVert)*shower.Direction().Y()+std::cos(angleToVert)*shower.Direction().Z()));
 
       srshower.plane[plane.ID().Plane].wirePitch = plane.WirePitch()/cosgamma;
@@ -695,7 +706,6 @@ namespace caf
   }
 
   void FillTrackChi2PID(const std::vector<art::Ptr<anab::ParticleID>> particleIDs,
-                        const geo::GeometryCore *geom,
                         caf::SRTrack& srtrack,
                         bool allowEmpty)
   {
@@ -815,7 +825,7 @@ namespace caf
   void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
                      const std::vector<art::Ptr<recob::Hit>> &hits,
                      bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
-                     const geo::GeometryCore *geom, const detinfo::DetectorPropertiesData &dprop,
+                     const detinfo::DetectorPropertiesData &dprop,
                      caf::SRTrack& srtrack,
                      bool allowEmpty)
   {

sbncode/CAFMaker/FillReco.h

@@ -1,4 +1,3 @@
-
 #ifndef CAF_FILLRECO_H
 #define CAF_FILLRECO_H
 
@@ -9,8 +8,7 @@
 #include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h"
 
 // LArSoft includes
-#include "larcore/Geometry/Geometry.h"
-#include "larcorealg/Geometry/GeometryCore.h"
+#include "larcorealg/Geometry/fwd.h"
 
 #include "lardataobj/RecoBase/PFParticle.h"
 #include "lardataobj/RecoBase/Shower.h"
@@ -59,7 +57,7 @@ namespace caf
   void FillShowerVars(const recob::Shower& shower,
                       const recob::Vertex* vertex,
                       const std::vector<art::Ptr<recob::Hit>> &hits,
-                      const geo::GeometryCore *geom,
+                      const geo::WireReadoutGeom& wireReadout,
                       unsigned producer,
                       caf::SRShower& srshower,
                       bool allowEmpty = false);
@@ -165,7 +163,6 @@ namespace caf
 
   void FillPlaneChi2PID(const anab::ParticleID &particle_id, caf::SRTrkChi2PID &srpid);
   void FillTrackChi2PID(const std::vector<art::Ptr<anab::ParticleID>> particleIDs,
-                        const geo::GeometryCore *geom,
                         caf::SRTrack& srtrack,
                         bool allowEmpty = false);
 
@@ -190,7 +187,7 @@ namespace caf
   void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
                      const std::vector<art::Ptr<recob::Hit>> &hits,
                      bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
-                     const geo::GeometryCore *geom, const detinfo::DetectorPropertiesData &dprop,
+                     const detinfo::DetectorPropertiesData &dprop,
                      caf::SRTrack& srtrack,
                      bool allowEmpty = false);
 
@@ -225,8 +222,8 @@ namespace caf
                   caf::SROpFlash &srflash,
                   bool allowEmpty = false);
   void FillCRTPMTMatch(const sbn::crt::CRTPMTMatching &match,
-		  caf::SRCRTPMTMatch &srmatch,
-		  bool allowEmpty = false);
+		       caf::SRCRTPMTMatch &srmatch,
+		       bool allowEmpty = false);
 
   void FillTPCPMTBarycenterMatch(const sbn::TPCPMTBarycenterMatch *matchInfo,
                            caf::SRSlice& slice);

sbncode/CAFMaker/FillTrue.cxx

@@ -1,6 +1,8 @@
 #include "FillTrue.h"
 
 #include "larcorealg/GeoAlgo/GeoAlgo.h"
+#include "larcorealg/Geometry/GeometryCore.h"
+#include "larcorealg/Geometry/WireReadoutGeom.h"
 #include "larevt/SpaceCharge/SpaceCharge.h"
 #include "larcore/CoreUtils/ServiceUtil.h"
 #include "RecoUtils/RecoUtils.h"
@@ -136,7 +138,8 @@ namespace caf {
   // Assumes truth matching and calo-points are filled
   void FillTrackCaloTruth(const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> &id_to_ide_map,
                           const std::vector<simb::MCParticle> &mc_particles,
-                          const geo::GeometryCore *geo,
+                          const geo::GeometryCore& geometry,
+                          const geo::WireReadoutGeom& wireReadout,
                           const detinfo::DetectorClocksData &clockData,
                           const spacecharge::SpaceCharge *sce,
                           caf::SRTrack& srtrack) {
@@ -163,7 +166,7 @@ namespace caf {
     const std::vector<std::pair<geo::WireID, const sim::IDE*>> &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID);
     std::map<unsigned, std::vector<const sim::IDE *>> chan_2_ides;
     for (auto const &ide_pair: match_ides) {
-      chan_2_ides[geo->PlaneWireToChannel(ide_pair.first)].push_back(ide_pair.second);
+      chan_2_ides[wireReadout.PlaneWireToChannel(ide_pair.first)].push_back(ide_pair.second);
     }
 
     // pre-compute partial ranges
@@ -267,18 +270,19 @@ namespace caf {
         // directly apply the true direction. We include the effect of space charge
         // on the true pitch here
         if (closest_dist >= 0. && direction.Mag() > 1e-4) {
-          geo::PlaneID plane(srtrack.producer, p.tpc, iplane);
-          float angletovert = geo->WireAngleToVertical(geo->View(plane), plane) - 0.5*::util::pi<>();
-          TVector3 loc_mdx_v = loc_nosce_v - direction * (geo->WirePitch(geo->View(plane)) / 2.);
-          TVector3 loc_pdx_v = loc_nosce_v + direction * (geo->WirePitch(geo->View(plane)) / 2.);
+          geo::PlaneID planeid(srtrack.producer, p.tpc, iplane);
+          geo::PlaneGeo const& plane = wireReadout.Plane(planeid);
+          float angletovert = wireReadout.WireAngleToVertical(plane.View(), planeid) - 0.5*::util::pi<>();
+          TVector3 loc_mdx_v = loc_nosce_v - direction * (plane.WirePitch() / 2.);
+          TVector3 loc_pdx_v = loc_nosce_v + direction * (plane.WirePitch() / 2.);
 
           // Convert types for helper functions
           geo::Point_t loc_mdx(loc_mdx_v.X(), loc_mdx_v.Y(), loc_mdx_v.Z());
           geo::Point_t loc_pdx(loc_pdx_v.X(), loc_pdx_v.Y(), loc_pdx_v.Z());
 
           // Map to wires
           if (sce && sce->EnableSimSpatialSCE()) { 
-            int corr = geo->TPC(plane).DriftDir().X();
+            int corr = geometry.TPC(plane.ID()).DriftDir().X();
 
             geo::Vector_t offset_m = sce->GetPosOffsets(loc_mdx);
             offset_m.SetX(offset_m.X()*corr); // convert from drift direction to detector direction
@@ -296,7 +300,7 @@ namespace caf {
           double cosgamma = std::abs(std::sin(angletovert)*dir.Y() + std::cos(angletovert)*dir.Z());
           double pitch;
           if (cosgamma) {
-            pitch = geo->WirePitch(geo->View(plane))/cosgamma;
+            pitch = plane.View()/cosgamma;
           }
           else {
             pitch = 0.;
@@ -306,7 +310,7 @@ namespace caf {
           geo::Point_t loc_atwires_alongpitch = loc_sce + dir*pitch; 
           geo::Point_t loc_alongpitch;
           if (sce && sce->EnableSimSpatialSCE()) { 
-            geo::Vector_t offset = sce->GetCalPosOffsets(loc_atwires_alongpitch, plane.TPC);
+            geo::Vector_t offset = sce->GetCalPosOffsets(loc_atwires_alongpitch, planeid.TPC);
             loc_alongpitch = loc_atwires_alongpitch + offset;
           }
           else loc_alongpitch = loc_atwires_alongpitch;
@@ -878,13 +882,13 @@ namespace caf {
     return ret;
   }
 
-  std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::GeometryCore &geo) {
+  std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout) {
     std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> ret;
 
     for (const art::Ptr<sim::SimChannel> sc : simchannels) {
       // Lookup the wire of this channel
       raw::ChannelID_t channel = sc->Channel();
-      std::vector<geo::WireID> maybewire = geo.ChannelToWire(channel);
+      std::vector<geo::WireID> maybewire = wireReadout.ChannelToWire(channel);
       geo::WireID thisWire; // Default constructor makes invalid wire
       if (maybewire.size()) thisWire = maybewire[0];