Add dt_hit / dt_ptoca linear fit information for hits on a track

ParticleID/CMakeLists.txt

@@ -15,7 +15,16 @@ cet_build_plugin(ParticleIDRead art::module
     LIBRARIES REG
       art_root_io::TFileService_service
       Offline::ParticleID
-
+
+      Offline::RecoDataProducts
+)
+
+cet_build_plugin(TrackDtDt art::module
+    REG_SOURCE src/TrackDtDt_module.cc
+    LIBRARIES REG
+      art_root_io::TFileService_service
+      Offline::Mu2eUtilities
+      Offline::ParticleID
       Offline::RecoDataProducts
 )
 

ParticleID/src/SConscript

@@ -21,6 +21,7 @@ libs =  [
     'mu2e_ConfigTools',
     'mu2e_DbTables',
     'mu2e_GeneralUtilities',
+    'mu2e_Mu2eUtilities',
     'art_Framework_Core',
     'art_Framework_Principal',
     'art_Framework_Services_Registry',

ParticleID/src/TrackDtDt_module.cc

@@ -0,0 +1,116 @@
+// Fit the track dt_{hit} / dt_{track fit poca time} distribution
+//
+// Michael MacKenzie, 2026
+
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "cetlib_except/exception.h"
+
+#include "fhiclcpp/types/Atom.h"
+#include "fhiclcpp/types/Sequence.h"
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/SubRun.h"
+#include "art/Framework/Principal/Run.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art_root_io/TFileService.h"
+
+#include "Offline/RecoDataProducts/inc/KalSeed.hh"
+#include "Offline/RecoDataProducts/inc/KalSeedDtDt.hh"
+#include "Offline/Mu2eUtilities/inc/LsqSums2.hh"
+
+#include "TH1.h"
+
+namespace mu2e {
+
+  //================================================================
+  class TrackDtDt : public art::EDProducer {
+
+    struct Config
+    {
+      using Name    = fhicl::Name;
+      using Comment = fhicl::Comment;
+      fhicl::Sequence<std::string> kalSeeds  { Name("kalSeeds") , Comment("KalSeed (Ptr) collection names") };
+      fhicl::Atom<int>             diagLevel { Name("diagLevel"), Comment("Diag Level")                    ,0 };
+    };
+
+    std::vector<std::string> kalSeeds_;
+    int diagLevel_;
+
+  public:
+    explicit TrackDtDt(const art::EDProducer::Table<Config>& config);
+    virtual void produce(art::Event& event);
+    KalSeedDtDt fitTrackDtDt(const KalSeed& seed);
+  };
+
+  //================================================================
+  TrackDtDt::TrackDtDt(const art::EDProducer::Table<Config>& config)
+    : art::EDProducer{config}
+    , kalSeeds_(config().kalSeeds())
+    , diagLevel_(config().diagLevel())
+  {
+    for(const auto& name : kalSeeds_) {
+      produces<mu2e::KalSeedDtDtCollection>(name);
+    }
+  }
+
+  //================================================================
+  KalSeedDtDt TrackDtDt::fitTrackDtDt(const KalSeed& seed) {
+    LsqSums2 fitter;
+    for(const auto& hit : seed.hits()) {
+      const double dt_hit = hit.time() - hit.TOTDriftTime();
+      const double dt_trk = hit.particleToca();
+      const double t_unc  = hit._udresid / hit._dvel;
+      const double weight = 1./t_unc;
+      fitter.addPoint(dt_trk, dt_hit, weight);
+    }
+    KalSeedDtDt result;
+    result.slope_     = fitter.dydx();
+    result.offset_    = fitter.y0();
+    result.slopeUnc_  = fitter.dydxErr();
+    result.dof_       = fitter.qn() - 2; // two parameters fitted
+    result.chisq_     = fitter.chi2Dof()*result.dof_;
+    return result;
+  }
+
+  //================================================================
+  void TrackDtDt::produce(art::Event& event) {
+
+    // Loop over all KalSeed collections
+    for(const auto& name : kalSeeds_) {
+
+      // Retrieve the KalSeed collection from the event, checking if it is a collection of KalSeed or KalSeedPtr
+      art::Handle<KalSeedCollection> seedHandle;
+      event.getByLabel(name, seedHandle);
+      art::Handle<KalSeedPtrCollection> seedPtrHandle;
+      const bool isSeedCollection = seedHandle.isValid();
+      if(!isSeedCollection) {
+        event.getByLabel(name, seedPtrHandle);
+        if(!seedPtrHandle.isValid()) {
+          throw cet::exception("RECO") << "TrackDtDt: No KalSeed or KalSeedPtr collection with label " << name << std::endl;
+        }
+      }
+
+      // Create the output collection
+      std::unique_ptr<KalSeedDtDtCollection> dtDtCol(new KalSeedDtDtCollection());
+
+      // Loop over all seeds in the collection and fit the dt_{hit} / dt_{track fit poca time} distribution
+      const auto nseeds = (isSeedCollection) ? seedHandle->size() : seedPtrHandle->size();
+      if(diagLevel_ > 0) std::cout << "[TrackDtDt::" << __func__ << "] Processing " << nseeds << " seeds from collection " << name << std::endl;
+      for(size_t iseed = 0; iseed < nseeds; ++iseed) {
+        const auto& seed = (isSeedCollection) ? seedHandle->at(iseed) : *seedPtrHandle->at(iseed);
+        if(diagLevel_ > 1) std::cout << "[TrackDtDt::" << __func__ << "] Processing seed " << iseed << " with " << seed.hits().size() << " hits" << std::endl;
+        dtDtCol->emplace_back(fitTrackDtDt(seed));
+      }
+
+      // Put the results into the event
+      event.put(std::move(dtDtCol), name);
+    }
+  }
+
+  //================================================================
+} // namespace mu2e
+
+DEFINE_ART_MODULE(mu2e::TrackDtDt)

RecoDataProducts/inc/KalSeedDtDt.hh

@@ -0,0 +1,26 @@
+//
+// Class representing the fit information of a dt_{straw hit time} / dt_{track hit time of closest approach}
+// This is useful for PID, as the wrong hypothesis will typically have a slope
+//  Original Author: Michael MacKenzie, 2026
+//
+#ifndef RecoDataProducts_KalSeedDtDt_HH
+#define RecoDataProducts_KalSeedDtDt_HH
+#include <Rtypes.h>
+namespace mu2e {
+  struct KalSeedDtDt {
+    Float_t slope_;
+    Float_t offset_;
+    Float_t slopeUnc_;
+    Float_t chisq_;
+    Int_t   dof_;
+    Float_t slope    () const { return slope_    ; }
+    Float_t offset   () const { return offset_   ; }
+    Float_t slopeUnc () const { return slopeUnc_ ; }
+    Float_t chisq    () const { return chisq_    ; }
+    Int_t   dof      () const { return dof_      ; }
+    KalSeedDtDt() : slope_(0.f), offset_(0.f), slopeUnc_(0.f), chisq_(0.f), dof_(0) {}
+  };
+
+  using KalSeedDtDtCollection = std::vector<mu2e::KalSeedDtDt>;
+}
+#endif

RecoDataProducts/src/classes.h

@@ -52,6 +52,7 @@
 #include "Offline/RecoDataProducts/inc/KalSeed.hh"
 #include "Offline/RecoDataProducts/inc/KalIntersection.hh"
 #include "Offline/RecoDataProducts/inc/KalSeedAssns.hh"
+#include "Offline/RecoDataProducts/inc/KalSeedDtDt.hh"
 #include "Offline/RecoDataProducts/inc/TrkCaloHitPID.hh"
 #include "Offline/RecoDataProducts/inc/TrkQual.hh"
 #include "Offline/RecoDataProducts/inc/RecoQual.hh"

RecoDataProducts/src/classes_def.xml

@@ -264,6 +264,11 @@
  <class name="mu2e::KKLineCollection" persistent="false"/>
  <class name="art::Wrapper<mu2e::KKLineCollection>"/>
 
+ <class name="mu2e::KalSeedDtDt"/>
+ <class name="mu2e::KalSeedDtDtCollection"/>
+ <class name="std::vector<mu2e::KalSeedDtDt>"/>
+ <class name="art::Wrapper<mu2e::KalSeedDtDtCollection>"/>
+
  <class name="mu2e::MVAStatus"/>
  <class name="mu2e::TrkQual"/>
  <class name="mu2e::TrkQualDetail"/>
@@ -452,7 +457,7 @@
  <class name="mu2e::STMFragmentSummary"/>
  <class name="mu2e::STMFragmentSummaryCollection"/>
  <class name="art::Wrapper<mu2e::STMFragmentSummaryCollection>"/>
- 
+
  <!--  ********* MSD ********* -->
   <class name="mu2e::MSDHit"/>
   <class name="mu2e::MSDHitCollection"/>