WireMod for Develop

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<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map;
+  std::map<int, std::vector<sbn::ReadoutIDE>>  id_to_ide_map;
   std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map;
   std::map<int, caf::HitsEnergy> id_to_hit_energy_map;
 
   if ( !isRealData ) {
     art::ServiceHandle<cheat::BackTrackerService> 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<recob::Vertex>(fmPFPart, evt,
              fParams.PFParticleLabel() + slice_tag_suff);
 
-    art::FindManyP<recob::Hit> fmTrackHit =
-      FindManyPStrict<recob::Hit>(slcTracks, evt,
+    art::FindManyP<recob::Hit, recob::TrackHitMeta> fmTrackHit =
+      FindManyPDStrict<recob::Hit, recob::TrackHitMeta>(slcTracks, evt,
           fParams.RecoTrackLabel() + slice_tag_suff);
 
     art::FindManyP<recob::Hit> 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) {

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
             )

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<art::Ptr<recob::Hit>> &hits,
+        const std::vector<const recob::TrackHitMeta*>& 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<recob::Hit> &h: hits) {
+        for (unsigned i_hit = 0; i_hit < hits.size(); i_hit++) {
+          const art::Ptr<recob::Hit> &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<art::Ptr<anab::Calorimetry>> &calos,
+                     const recob::Track& track,
                      const std::vector<art::Ptr<recob::Hit>> &hits,
+                     const std::vector<const recob::TrackHitMeta*>& 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]);
       }
     }
 

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<art::Ptr<recob::Hit>> &hits,
+                     const std::vector<const recob::TrackHitMeta*>& 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<sbn::ScatterClosestApproach> closestapproach,
@@ -257,9 +263,12 @@ namespace caf
                         bool allowEmpty = false);
 
   void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
+                     const recob::Track& track,
                      const std::vector<art::Ptr<recob::Hit>> &hits,
+                     const std::vector<const recob::TrackHitMeta*>& 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);