diff --git a/ParticleID/CMakeLists.txt b/ParticleID/CMakeLists.txt index 09e0c0f1f5..01d72cb9dd 100644 --- a/ParticleID/CMakeLists.txt +++ b/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 ) diff --git a/ParticleID/src/SConscript b/ParticleID/src/SConscript index 6297a3b3f2..7a45d2954c 100644 --- a/ParticleID/src/SConscript +++ b/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', diff --git a/ParticleID/src/TrackDtDt_module.cc b/ParticleID/src/TrackDtDt_module.cc new file mode 100644 index 0000000000..a06e917fb2 --- /dev/null +++ b/ParticleID/src/TrackDtDt_module.cc @@ -0,0 +1,116 @@ +// Fit the track dt_{hit} / dt_{track fit poca time} distribution +// +// Michael MacKenzie, 2026 + +#include +#include +#include + +#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 kalSeeds { Name("kalSeeds") , Comment("KalSeed (Ptr) collection names") }; + fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diag Level") ,0 }; + }; + + std::vector kalSeeds_; + int diagLevel_; + + public: + explicit TrackDtDt(const art::EDProducer::Table& config); + virtual void produce(art::Event& event); + KalSeedDtDt fitTrackDtDt(const KalSeed& seed); + }; + + //================================================================ + TrackDtDt::TrackDtDt(const art::EDProducer::Table& config) + : art::EDProducer{config} + , kalSeeds_(config().kalSeeds()) + , diagLevel_(config().diagLevel()) + { + for(const auto& name : kalSeeds_) { + produces(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 seedHandle; + event.getByLabel(name, seedHandle); + art::Handle 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 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) diff --git a/RecoDataProducts/inc/KalSeedDtDt.hh b/RecoDataProducts/inc/KalSeedDtDt.hh new file mode 100644 index 0000000000..d2d3aa17ce --- /dev/null +++ b/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 +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; +} +#endif diff --git a/RecoDataProducts/src/classes.h b/RecoDataProducts/src/classes.h index 1eda0cc472..10f747d571 100644 --- a/RecoDataProducts/src/classes.h +++ b/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" diff --git a/RecoDataProducts/src/classes_def.xml b/RecoDataProducts/src/classes_def.xml index 00202da02d..89d341a131 100644 --- a/RecoDataProducts/src/classes_def.xml +++ b/RecoDataProducts/src/classes_def.xml @@ -264,6 +264,11 @@ + + + + + @@ -452,7 +457,7 @@ - +