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..83bd0f0d00 --- /dev/null +++ b/ParticleID/src/TrackDtDt_module.cc @@ -0,0 +1,187 @@ +// 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" +#include "TGraph.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 doHistograms{ Name("doHistograms"), Comment("Make histograms") , false }; + fhicl::Atom diagLevel { Name("diagLevel") , Comment("Diag Level") ,0 }; + }; + + std::vector kalSeeds_; + bool doHistograms_; + int diagLevel_; + + // Histograms + TH1* h_slope_ = nullptr; // fit results + TH1* h_offset_ = nullptr; + TH1* h_slopeUnc_ = nullptr; + TH1* h_chisq_ = nullptr; + TH1* h_dof_ = nullptr; + TH1* h_chisqDof_ = nullptr; + TH1* h_trk_chisq_ = nullptr; // track parameters + TH1* h_trk_fitcon_ = nullptr; + TGraph* g_t_hit_vs_z = nullptr; + TGraph* g_t_poca_vs_z = nullptr; + TGraph* g_t_hit_vs_t_poca = nullptr; + + public: + explicit TrackDtDt(const art::EDProducer::Table& config); + virtual void produce(art::Event& event); + KalSeedDtDt fitTrackDtDt(const KalSeed& seed); + + void bookHistograms(); + void fillHistograms(const KalSeedDtDt& dtDt, const KalSeed& seed); + }; + + //================================================================ + TrackDtDt::TrackDtDt(const art::EDProducer::Table& config) + : art::EDProducer{config} + , kalSeeds_(config().kalSeeds()) + , doHistograms_(config().doHistograms()) + , diagLevel_(config().diagLevel()) + { + for(const auto& name : kalSeeds_) { + produces(name); + } + if(doHistograms_) bookHistograms(); + } + + //================================================================ + void TrackDtDt::bookHistograms() { + art::ServiceHandle tfs; + h_slope_ = tfs->make("h_slope", "Slope of dt_{hit} vs dt_{track fit poca time};Slope;Entries", 100, -2., 4.); + h_offset_ = tfs->make("h_offset", "Offset of dt_{hit} vs dt_{track fit poca time};Offset [ns];Entries", 100, -100., 100.); + h_slopeUnc_ = tfs->make("h_slopeUnc", "Uncertainty of slope;Slope Uncertainty;Entries", 100, 0., 0.1); + h_chisq_ = tfs->make("h_chisq", "Chi2 of fit;Chi2;Entries", 100, 0., 100.); + h_dof_ = tfs->make("h_dof", "Degrees of freedom of fit;DOF;Entries", 100, 0., 100.); + h_chisqDof_ = tfs->make("h_chisqDof", "Chi2/DOF of fit;Chi2/DOF;Entries", 100, 0., 5.); + h_trk_chisq_ = tfs->make("h_trk_chisq", "Track chi2;Chi2;Entries", 100, 0., 100.); + h_trk_fitcon_ = tfs->make("h_trk_fitcon", "Track fit convergence;Convergence;Entries", 100, 0., 1.); + g_t_hit_vs_z = tfs->makeAndRegister("t_hit_vs_z", "t_{hit} vs. z;z [mm];t_{hit} [ns]"); + g_t_poca_vs_z = tfs->makeAndRegister("t_poca_vs_z", "t_{POCA} vs. z;z [mm];t_{POCA} [ns]"); + g_t_hit_vs_t_poca = tfs->makeAndRegister("t_hit_vs_t_poca", "t_{hit} vs. t_{POCA};t_{POCA} [ns];t_{hit} [ns]"); + } + + //================================================================ + void TrackDtDt::fillHistograms(const KalSeedDtDt& dtDt, const KalSeed& seed) { + if(!doHistograms_) return; + h_slope_->Fill(dtDt.slope_); + h_offset_->Fill(dtDt.offset_); + h_slopeUnc_->Fill(dtDt.slopeUnc_); + h_chisq_->Fill(dtDt.chisq_); + h_dof_->Fill(dtDt.dof_); + h_chisqDof_->Fill((dtDt.dof_ > 0) ? dtDt.chisq_/dtDt.dof_ : 0.); + h_trk_chisq_->Fill(seed.chisquared()); + h_trk_fitcon_->Fill(seed.fitConsistency()); + + // Fill one example event + if(g_t_hit_vs_z->GetN() == 0) { + for(const auto& hit : seed.hits()) { + const double t_trk = hit.particleToca(); // estimated time of the particle + const double t_hit = t_trk + hit.fitDt(); // add dt to get the hit time + const double z_trk = hit._upoca.z(); + g_t_hit_vs_z ->AddPoint(z_trk, t_hit); + g_t_poca_vs_z ->AddPoint(z_trk, t_trk); + g_t_hit_vs_t_poca->AddPoint(t_trk, t_hit); + } + } + } + + //================================================================ + KalSeedDtDt TrackDtDt::fitTrackDtDt(const KalSeed& seed) { + LsqSums2 fitter; + for(const auto& hit : seed.hits()) { + const double t_trk = hit.particleToca(); // estimated time of the particle + const double t_hit = t_trk + hit.fitDt(); // add dt to get the hit time + const double t_var = hit.fitTocaVar(); // variance of the track poca time and the hit time + const double weight = 1./t_var; // inverse variance weighting + fitter.addPoint(t_trk, t_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_ = (result.dof_ > 0) ? fitter.chi2Dof()*result.dof_ : 0.f; // chi2Dof() not defined for dof <= 0 + + if(diagLevel_ > 1) { + std::cout << "[TrackDtDt::" << __func__ << "] Fit results for seed with " << seed.hits().size() << " hits:\n" + << " slope = " << result.slope_ << " +/- " << result.slopeUnc_ << "\n" + << " offset = " << result.offset_ << "\n" + << " chi2 = " << result.chisq_ << "\n" + << " dof = " << result.dof_ << "\n"; + } + if(doHistograms_) fillHistograms(result, seed); + 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/ParticleID/test/trkdtdt.fcl b/ParticleID/test/trkdtdt.fcl new file mode 100644 index 0000000000..9dcdbec8e5 --- /dev/null +++ b/ParticleID/test/trkdtdt.fcl @@ -0,0 +1,33 @@ +# Test TrackDtDt_module.cc + +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" + +process_name : TestTrkDtDt + +services : @local::Services.Reco +services.GeometryService.bFieldFile : "Offline/Mu2eG4/geom/bfgeom_reco_v01.txt" +services.scheduler.wantSummary: true + +source: { + module_type: RootInput + inputCommands: [ "keep *_*_*_*" + ] +} + +physics : { + producers : { + TrackDtDt : { + module_type : TrackDtDt + kalSeeds : [ "KKDe" ] + diagLevel : 0 + doHistograms: true + } + } + RecoPath : [ TrackDtDt ] + end_paths : [ ] + trigger_paths : [ RecoPath ] +} + +#include "Production/Validation/database.fcl" 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 @@ - +