Merge pull request #557 from SBNSoftware/feature/hlay_calib_ntuples_crt_t0

Calib NTuple CRT Timing Additions

Author: kjplows

sbncode/Calibration/TrackCaloSkimmer.h

@@ -63,6 +63,8 @@
 #include "sbnobj/Common/Calibration/TrackCaloSkimmerObj.h"
 #include "sbnobj/Common/CRT/CRTHitT0TaggingInfo.hh"
 #include "sbnobj/Common/CRT/CRTHitT0TaggingTruthInfo.hh"
+#include "sbnobj/SBND/CRT/CRTTrack.hh"
+#include "sbnobj/SBND/CRT/CRTSpacePoint.hh"
 
 #include "ITCSSelectionTool.h"
 
@@ -98,9 +100,12 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer {
     double  t0Pandora;
     double  t0CRTTrack;
     double  t0CRTHit;
+    double  t0CRTSpacePoint;
     bool    hasT0Pandora;
     bool    hasT0CRTTrack;
     bool    hasT0CRTHit;
+    bool    hasT0CRTSpacePoint;
+    double  crtMatchingScore;
   };
 
   // Internal data struct
@@ -178,7 +183,8 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer {
     const geo::WireReadoutGeom *wireReadout,
     const detinfo::DetectorClocksData &dclock,
     const cheat::BackTrackerService *bt_serv,
-    const detinfo::DetectorPropertiesData &dprop);
+    const detinfo::DetectorPropertiesData &dprop,
+    const float &xshift);
   void FillTrackCRTHitInfo(const std::vector<art::Ptr<sbn::crt::CRTHitT0TaggingInfo>> &tag);
 
   // helpers
@@ -192,6 +198,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer {
   art::InputTag fPFPT0producer;
   art::InputTag fCRTTrackT0producer;
   art::InputTag fCRTHitT0producer;
+  art::InputTag fCRTSpacePointT0producer;
   art::InputTag fCALOproducer;
   art::InputTag fTRKproducer;
   art::InputTag fTRKHMproducer;
@@ -216,6 +223,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer {
   double fTopCRTDistanceCutPassing;
   double fSideCRTDistanceCutStopping;
   double fSideCRTDistanceCutPassing;  
+  bool fAllowShowerLikePFPs;
 
   // tools
   std::vector<std::unique_ptr<sbn::ITCSSelectionTool>> fSelectionTools;

sbncode/Calibration/TrackCaloSkimmer_module.cc

@@ -11,7 +11,6 @@
 ////////////////////////////////////////////////////////////////////////
 
 #include "TrackCaloSkimmer.h"
-#include "sbnobj/SBND/CRT/CRTTrack.hh"
 
 #include "art/Utilities/make_tool.h"
 
@@ -70,6 +69,7 @@ sbn::TrackCaloSkimmer::TrackCaloSkimmer(fhicl::ParameterSet const& p)
   fPFPproducer  = p.get< art::InputTag > ("PFPproducer","pandoraGausCryo0");
   fPFPT0producer = p.get< art::InputTag > ("PFPT0producer", "pandoraGausCryo0");
   fCRTTrackT0producer = p.get< art::InputTag >("CRTTrackT0producer", "crttrackmatching");
+  fCRTSpacePointT0producer = p.get< art::InputTag >("CRTSpacePointT0producer", "crtspacepointmatching");
   fCRTHitT0producer = p.get< art::InputTag >("CRTHitT0producer", "CRTT0Tagging");
 
   fCALOproducer = p.get< art::InputTag > ("CALOproducer");
@@ -92,7 +92,8 @@ sbn::TrackCaloSkimmer::TrackCaloSkimmer(fhicl::ParameterSet const& p)
   fTopCRTDistanceCutPassing = p.get<double>("TopCRTDistanceCut_throughgoing", 100.);
   fSideCRTDistanceCutStopping = p.get<double>("SideCRTDistanceCut_stopping", 100.);
   fSideCRTDistanceCutPassing = p.get<double>("SideCRTDistanceCut_throughgoing", 100.);
-  
+  fAllowShowerLikePFPs = p.get<bool>("AllowShowerLikePFPs", false);
+
   if (fTailFitResidualRange > 10.) {
     std::cout << "sbn::TrackCaloSkimmer: Bad tail fit residual range config :(" << fTailFitResidualRange << "). Fits will not be meaningful.\n";
   }
@@ -231,6 +232,8 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e)
   //
   // Tracks (SBND style)
   art::FindManyP<sbnd::crt::CRTTrack, anab::T0> fmT0CRTTrack(tracks, e, fCRTTrackT0producer);
+  // SpacePoints (SBND style)
+  art::FindManyP<sbnd::crt::CRTSpacePoint, anab::T0> fmT0CRTSpacePoint(tracks, e, fCRTSpacePointT0producer);
 
   // Hits (ICARUS style)
   art::FindManyP<anab::T0> fmT0CRTHit(tracks, e, fCRTHitT0producer);
@@ -313,6 +316,9 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e)
   for (art::Ptr<recob::PFParticle> p_pfp: PFParticleList) {
     const recob::PFParticle &pfp = *p_pfp;
 
+    if(p_pfp->PdgCode() == 11 && !fAllowShowerLikePFPs)
+      continue;
+
     const std::vector<art::Ptr<recob::Track>> thisTrack = fmTracks.at(p_pfp.key());
     if (thisTrack.size() != 1)
       continue;
@@ -345,8 +351,9 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e)
     }
 
     // Collect T0s
-    bool hasT0 = false, hasPFPT0 = false, hasCRTTrackT0 = false, hasCRTHitT0 = false;
+    bool hasT0 = false, hasPFPT0 = false, hasCRTTrackT0 = false, hasCRTHitT0 = false, hasCRTSpacePointT0 = false;
     int whicht0 = -1;
+    float crtMatchingScore = std::numeric_limits<float>::signaling_NaN();
 
     double t0PFP = std::numeric_limits<float>::signaling_NaN();
     if (fmT0PFP.isValid() && fmT0PFP.at(p_pfp.key()).size()) {
@@ -356,12 +363,23 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e)
     }
 
     double t0CRTTrack = std::numeric_limits<float>::signaling_NaN();
+    double t0CRTTrackScore = std::numeric_limits<float>::signaling_NaN();
     if (fmT0CRTTrack.isValid() && fmT0CRTTrack.at(trkPtr.key()).size()) {
       t0CRTTrack = fmT0CRTTrack.data(trkPtr.key()).at(0)->Time();
+      t0CRTTrackScore = fmT0CRTTrack.data(trkPtr.key()).at(0)->TriggerConfidence();
       hasCRTTrackT0 = true;
     }
 
+    double t0CRTSpacePoint = std::numeric_limits<float>::signaling_NaN();
+    double t0CRTSpacePointScore = std::numeric_limits<float>::signaling_NaN();
+    if (fmT0CRTSpacePoint.isValid() && fmT0CRTSpacePoint.at(trkPtr.key()).size()) {
+      t0CRTSpacePoint = fmT0CRTSpacePoint.data(trkPtr.key()).at(0)->Time();
+      t0CRTSpacePointScore = fmT0CRTSpacePoint.data(trkPtr.key()).at(0)->TriggerConfidence();
+      hasCRTSpacePointT0 = true;
+    }
+
     double t0CRTHit = std::numeric_limits<float>::signaling_NaN();
+    double t0CRTHitScore = std::numeric_limits<float>::signaling_NaN();
     if (fIncludeCRTHitTagging && fmT0CRTHit.isValid() && fmT0CRTHit.at(trkPtr.key()).size()) {
       const sbn::crt::CRTHitT0TaggingInfo &tag = *fmCRTHitT0TaggingInfo.at(trkPtr.key()).at(0);
       double time = fmT0CRTHit.at(trkPtr.key()).at(0)->Time();
@@ -384,23 +402,42 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e)
 
       if (!crtHitSysRejected && !crtHitDistanceRejected) {
         t0CRTHit = time;
+        t0CRTHitScore = tag.Distance;
         hasCRTHitT0 = true;
       }
     }
 
-    T0TimingInfo thisTrackTimingInfo = {t0PFP, t0CRTTrack, t0CRTHit, hasPFPT0, hasCRTTrackT0, hasCRTHitT0};
-    hasT0 = hasPFPT0 || hasCRTTrackT0 || hasCRTHitT0;
+    hasT0 = hasPFPT0 || hasCRTTrackT0 || hasCRTHitT0 || hasCRTSpacePointT0;
 
     // "whicht0" should reflect the T0 used for the reconstruction of the drift coordinate.
-    if(!hasT0) whicht0 = -1 ;
+    if(!hasT0) whicht0 = -1;
     // In this way, if a track is T0 tagged from PFP and CRT tagged, which T0 reflects the PFP Tag.
-    else if (hasPFPT0) whicht0 = 0 ;
-    else if (hasCRTTrackT0) whicht0 = 1 ;
-    else if (hasCRTHitT0) whicht0 = 2 ;
+    else if (hasPFPT0) whicht0 = 0;
+    else if (hasCRTTrackT0)
+      {
+        whicht0 = 1;
+        crtMatchingScore = t0CRTTrackScore;
+      }
+    else if (hasCRTHitT0)
+      {
+        whicht0 = 2;
+        crtMatchingScore = t0CRTHitScore;
+      }
+    else if (hasCRTSpacePointT0)
+      {
+        whicht0 = 3;
+        crtMatchingScore = t0CRTSpacePointScore;
+      }
+
+    T0TimingInfo thisTrackTimingInfo = {t0PFP, t0CRTTrack, t0CRTHit, t0CRTSpacePoint, hasPFPT0, hasCRTTrackT0, hasCRTHitT0, hasCRTSpacePointT0, crtMatchingScore};
 
     if (fRequireT0 && !hasT0) continue;
 
-    if (fVerbose) std::cout << "Processing new track! ID: " << trkPtr->ID() << " time: " << t0PFP << " timeCRTTrack: " << t0CRTTrack << " timeCRTHit "<<t0CRTHit<<std::endl;
+    if (fVerbose) std::cout << "Processing new track! ID: " << trkPtr->ID() << "\n"
+                            << "\tPandora time: " << t0PFP << "\n"
+                            << "\tCRTTrack (SBND) time: " << t0CRTTrack << "\n"
+                            << "\tCRTHit (ICARUS) time: " << t0CRTHit << "\n"
+                            << "\tCRTSpacePoint (SBND) time: " << t0CRTSpacePoint << std::endl;
 
     // Reset the track object
     *fTrack = sbn::TrackInfo();
@@ -1124,6 +1161,8 @@ void sbn::TrackCaloSkimmer::FillTrack(const recob::Track &track,
   fTrack->t0PFP = t0Info.t0Pandora;
   fTrack->t0CRTTrack = t0Info.t0CRTTrack;
   fTrack->t0CRTHit = t0Info.t0CRTHit;
+  fTrack->t0CRTSpacePoint = t0Info.t0CRTSpacePoint;
+  fTrack->crtMatchingScore = t0Info.crtMatchingScore;
   fTrack->id = track.ID();
   fTrack->clear_cosmic_muon = pfp.Parent() == recob::PFParticle::kPFParticlePrimary;
 
@@ -1135,30 +1174,48 @@ void sbn::TrackCaloSkimmer::FillTrack(const recob::Track &track,
   fTrack->dir.x = track.StartDirection().X();
   fTrack->dir.y = track.StartDirection().Y();
   fTrack->dir.z = track.StartDirection().Z(); 
-  //If track is only CRTt0 tagged, undo the assumed trigger correction
-  if (t0Info.hasT0Pandora) {
-    fTrack->start.x = track.Start().X();
-    fTrack->end.x = track.End().X();
-  } else if (t0Info.hasT0CRTTrack) {
-    int driftDir = geo->TPC(hits[0]->WireID()).DriftDir().X();
-    const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature()));
-    fTrack->start.x = track.Start().X() + driftDir*driftv*t0Info.t0CRTTrack*1e-3;
-    fTrack->end.x = track.End().X() + driftDir*driftv*t0Info.t0CRTTrack*1e-3;
-  } else if (t0Info.hasT0CRTHit){ 
-    // If the track does not have a a Pandora T0, the tracks will always be either on the left or (ex Or) right of the cathode. 
-    int driftDir = geo->TPC(hits[0]->WireID()).DriftDir().X();
-    const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature()));
-    fTrack->start.x = track.Start().X() + driftDir*driftv*t0Info.t0CRTHit*1e-3;
-    fTrack->end.x = track.End().X() + driftDir*driftv*t0Info.t0CRTHit*1e-3;
-  }
+
+  if(t0Info.hasT0Pandora)
+    {
+      fTrack->start.x = track.Start().X();
+      fTrack->end.x = track.End().X();
+    }
+  else
+    {
+      // If CRT T0 tagged we can shift the x positions appropriately
+      int driftDir = geo->TPC(hits[0]->WireID()).DriftDir().X();
+      const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature()));
+
+      if(t0Info.hasT0CRTTrack)
+        {
+          const double xshift = driftDir*driftv*t0Info.t0CRTTrack*1e-3;
+          fTrack->start.x = track.Start().X() + xshift;
+          fTrack->end.x = track.End().X() + xshift;
+          fTrack->xShiftCRT = xshift;
+        }
+      else if(t0Info.hasT0CRTHit)
+        {
+          const double xshift = driftDir*driftv*t0Info.t0CRTHit*1e-3;
+          fTrack->start.x = track.Start().X() + xshift;
+          fTrack->end.x = track.End().X() + xshift;
+          fTrack->xShiftCRT = xshift;
+        }
+      else if(t0Info.hasT0CRTSpacePoint)
+        {
+          const double xshift = driftDir*driftv*t0Info.t0CRTSpacePoint*1e-3;
+          fTrack->start.x = track.Start().X() + xshift;
+          fTrack->end.x = track.End().X() + xshift;
+          fTrack->xShiftCRT = xshift;
+        }
+    }
 
   if (hits.size() > 0) {
     fTrack->cryostat = hits[0]->WireID().Cryostat;
   }
 
   // Fill each hit
   for (unsigned i_hit = 0; i_hit < hits.size(); i_hit++) {
-    sbn::TrackHitInfo hinfo = MakeHit(*hits[i_hit], hits[i_hit].key(), *thms[i_hit], track, t0Info, sps[i_hit], calo, geo, wireReadout, clock_data, bt_serv, dprop);
+    sbn::TrackHitInfo hinfo = MakeHit(*hits[i_hit], hits[i_hit].key(), *thms[i_hit], track, t0Info, sps[i_hit], calo, geo, wireReadout, clock_data, bt_serv, dprop, fTrack->xShiftCRT);
     if (hinfo.h.plane == 0) {
       fTrack->hits0.push_back(hinfo);
     }
@@ -1325,7 +1382,8 @@ sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit(const recob::Hit &hit,
     const geo::WireReadoutGeom *wireReadout,
     const detinfo::DetectorClocksData &dclock,
     const cheat::BackTrackerService *bt_serv,
-    const detinfo::DetectorPropertiesData &dprop) {
+    const detinfo::DetectorPropertiesData &dprop,
+    const float &xshift) {
 
   // TrackHitInfo to save
   sbn::TrackHitInfo hinfo;
@@ -1400,13 +1458,6 @@ sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit(const recob::Hit &hit,
     }
   }
 
-  // This is needed to reconstrut drift coordinate using a different Time
-  int driftDir = geo->TPC(hit.WireID()).DriftDir().X();
-  const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature()));
-  double anodeDistance = std::numeric_limits<float>::signaling_NaN();
-  if(t0Info.hasT0CRTHit) anodeDistance = (hit.PeakTime()-dclock.Time2Tick(dclock.TriggerTime())-t0Info.t0CRTHit*1e-3/dclock.TPCClock().TickPeriod())*dclock.TPCClock().TickPeriod()*driftv;
-  double wirePlaneX = wireReadout->Plane(geo::PlaneID(hit.WireID().Cryostat, hit.WireID().TPC, hit.WireID().Plane)).GetCenter().X();
-  double recoX = wirePlaneX - driftDir*anodeDistance;
   // Information from the TrackHitMeta
   bool badhit = (thm.Index() == std::numeric_limits<unsigned int>::max()) ||
                     (!trk.HasValidPoint(thm.Index()));
@@ -1419,7 +1470,7 @@ sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit(const recob::Hit &hit,
 
     // The tp X coordinate is reconstructed only if it does not have a PandoraT0.
     hinfo.tp.x = loc.X();
-    if(t0Info.hasT0CRTHit && !t0Info.hasT0Pandora && !isnan(hinfo.tp.x)) hinfo.tp.x = recoX;
+    if((t0Info.hasT0CRTTrack || t0Info.hasT0CRTHit || t0Info.hasT0CRTSpacePoint) && !t0Info.hasT0Pandora && !isnan(hinfo.tp.x)) hinfo.tp.x = loc.X() + xshift;
     hinfo.tp.y = loc.Y();
     hinfo.tp.z = loc.Z();
 
@@ -1451,7 +1502,7 @@ sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit(const recob::Hit &hit,
   if (sp) {
     hinfo.h.sp.x = sp->position().x();
     // If the track has a Pandora T0, do not displace, track is already reconstructed at the correct position
-    if(t0Info.hasT0CRTHit && !t0Info.hasT0Pandora) hinfo.h.sp.x = sp->position().x() + driftDir*driftv*t0Info.t0CRTHit*1e-3;
+    if((t0Info.hasT0CRTTrack || t0Info.hasT0CRTHit || t0Info.hasT0CRTSpacePoint) && !t0Info.hasT0Pandora) hinfo.h.sp.x = sp->position().x() + xshift;
     hinfo.h.sp.y = sp->position().y();
     hinfo.h.sp.z = sp->position().z();
     hinfo.h.hasSP = true;

sbncode/Calibration/fcl/sbnd/sbnd_trackcalo_skimmer.fcl

@@ -83,6 +83,6 @@ a2c_selection: {
 }
 
 caloskim_nodigits_goldentracks: @local::caloskim_nodigits
-caloskim_nodigits_goldentracks.SelectionTools: [@local::stopping_selection, @local::a2c_selection, @local::throughgoing_selection]
+caloskim_nodigits_goldentracks.SelectionTools: [ @local::stopping_selection, @local::a2c_selection, @local::throughgoing_selection ]
 
 END_PROLOG