Skip to content

Commit 692e0df

Browse files
oksuzianoksuzian
authored
Merge pull request #1873 from michaelmackenzie/TrkDtDt
Add dt_hit / dt_ptoca linear fit information for hits on a track
2 parents f853866 + 3eff1ef commit 692e0df

7 files changed

Lines changed: 263 additions & 1 deletion

File tree

‎ParticleID/CMakeLists.txt‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,16 @@ cet_build_plugin(ParticleIDRead art::module
1515
LIBRARIES REG
1616
art_root_io::TFileService_service
1717
Offline::ParticleID
18-
18+
19+
Offline::RecoDataProducts
20+
)
21+
22+
cet_build_plugin(TrackDtDt art::module
23+
REG_SOURCE src/TrackDtDt_module.cc
24+
LIBRARIES REG
25+
art_root_io::TFileService_service
26+
Offline::Mu2eUtilities
27+
Offline::ParticleID
1928
Offline::RecoDataProducts
2029
)
2130

‎ParticleID/src/SConscript‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@ libs = [
2121
'mu2e_ConfigTools',
2222
'mu2e_DbTables',
2323
'mu2e_GeneralUtilities',
24+
'mu2e_Mu2eUtilities',
2425
'art_Framework_Core',
2526
'art_Framework_Principal',
2627
'art_Framework_Services_Registry',

‎ParticleID/src/TrackDtDt_module.cc‎

Lines changed: 187 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
// Fit the track dt_{hit} / dt_{track fit poca time} distribution
2+
//
3+
// Michael MacKenzie, 2026
4+
5+
#include <iostream>
6+
#include <string>
7+
#include <vector>
8+
9+
#include "cetlib_except/exception.h"
10+
11+
#include "fhiclcpp/types/Atom.h"
12+
#include "fhiclcpp/types/Sequence.h"
13+
#include "art/Framework/Core/EDProducer.h"
14+
#include "art/Framework/Principal/Event.h"
15+
#include "art/Framework/Principal/SubRun.h"
16+
#include "art/Framework/Principal/Run.h"
17+
#include "art/Framework/Principal/Handle.h"
18+
#include "art_root_io/TFileService.h"
19+
20+
#include "Offline/RecoDataProducts/inc/KalSeed.hh"
21+
#include "Offline/RecoDataProducts/inc/KalSeedDtDt.hh"
22+
#include "Offline/Mu2eUtilities/inc/LsqSums2.hh"
23+
24+
#include "TH1.h"
25+
#include "TGraph.h"
26+
27+
namespace mu2e {
28+
29+
//================================================================
30+
class TrackDtDt : public art::EDProducer {
31+
32+
struct Config
33+
{
34+
using Name = fhicl::Name;
35+
using Comment = fhicl::Comment;
36+
fhicl::Sequence<std::string> kalSeeds { Name("kalSeeds") , Comment("KalSeed (Ptr) collection names") };
37+
fhicl::Atom<bool> doHistograms{ Name("doHistograms"), Comment("Make histograms") , false };
38+
fhicl::Atom<int> diagLevel { Name("diagLevel") , Comment("Diag Level") ,0 };
39+
};
40+
41+
std::vector<std::string> kalSeeds_;
42+
bool doHistograms_;
43+
int diagLevel_;
44+
45+
// Histograms
46+
TH1* h_slope_ = nullptr; // fit results
47+
TH1* h_offset_ = nullptr;
48+
TH1* h_slopeUnc_ = nullptr;
49+
TH1* h_chisq_ = nullptr;
50+
TH1* h_dof_ = nullptr;
51+
TH1* h_chisqDof_ = nullptr;
52+
TH1* h_trk_chisq_ = nullptr; // track parameters
53+
TH1* h_trk_fitcon_ = nullptr;
54+
TGraph* g_t_hit_vs_z = nullptr;
55+
TGraph* g_t_poca_vs_z = nullptr;
56+
TGraph* g_t_hit_vs_t_poca = nullptr;
57+
58+
public:
59+
explicit TrackDtDt(const art::EDProducer::Table<Config>& config);
60+
virtual void produce(art::Event& event);
61+
KalSeedDtDt fitTrackDtDt(const KalSeed& seed);
62+
63+
void bookHistograms();
64+
void fillHistograms(const KalSeedDtDt& dtDt, const KalSeed& seed);
65+
};
66+
67+
//================================================================
68+
TrackDtDt::TrackDtDt(const art::EDProducer::Table<Config>& config)
69+
: art::EDProducer{config}
70+
, kalSeeds_(config().kalSeeds())
71+
, doHistograms_(config().doHistograms())
72+
, diagLevel_(config().diagLevel())
73+
{
74+
for(const auto& name : kalSeeds_) {
75+
produces<mu2e::KalSeedDtDtCollection>(name);
76+
}
77+
if(doHistograms_) bookHistograms();
78+
}
79+
80+
//================================================================
81+
void TrackDtDt::bookHistograms() {
82+
art::ServiceHandle<art::TFileService> tfs;
83+
h_slope_ = tfs->make<TH1D>("h_slope", "Slope of dt_{hit} vs dt_{track fit poca time};Slope;Entries", 100, -2., 4.);
84+
h_offset_ = tfs->make<TH1D>("h_offset", "Offset of dt_{hit} vs dt_{track fit poca time};Offset [ns];Entries", 100, -100., 100.);
85+
h_slopeUnc_ = tfs->make<TH1D>("h_slopeUnc", "Uncertainty of slope;Slope Uncertainty;Entries", 100, 0., 0.1);
86+
h_chisq_ = tfs->make<TH1D>("h_chisq", "Chi2 of fit;Chi2;Entries", 100, 0., 100.);
87+
h_dof_ = tfs->make<TH1D>("h_dof", "Degrees of freedom of fit;DOF;Entries", 100, 0., 100.);
88+
h_chisqDof_ = tfs->make<TH1D>("h_chisqDof", "Chi2/DOF of fit;Chi2/DOF;Entries", 100, 0., 5.);
89+
h_trk_chisq_ = tfs->make<TH1D>("h_trk_chisq", "Track chi2;Chi2;Entries", 100, 0., 100.);
90+
h_trk_fitcon_ = tfs->make<TH1D>("h_trk_fitcon", "Track fit convergence;Convergence;Entries", 100, 0., 1.);
91+
g_t_hit_vs_z = tfs->makeAndRegister<TGraph>("t_hit_vs_z", "t_{hit} vs. z;z [mm];t_{hit} [ns]");
92+
g_t_poca_vs_z = tfs->makeAndRegister<TGraph>("t_poca_vs_z", "t_{POCA} vs. z;z [mm];t_{POCA} [ns]");
93+
g_t_hit_vs_t_poca = tfs->makeAndRegister<TGraph>("t_hit_vs_t_poca", "t_{hit} vs. t_{POCA};t_{POCA} [ns];t_{hit} [ns]");
94+
}
95+
96+
//================================================================
97+
void TrackDtDt::fillHistograms(const KalSeedDtDt& dtDt, const KalSeed& seed) {
98+
if(!doHistograms_) return;
99+
h_slope_->Fill(dtDt.slope_);
100+
h_offset_->Fill(dtDt.offset_);
101+
h_slopeUnc_->Fill(dtDt.slopeUnc_);
102+
h_chisq_->Fill(dtDt.chisq_);
103+
h_dof_->Fill(dtDt.dof_);
104+
h_chisqDof_->Fill((dtDt.dof_ > 0) ? dtDt.chisq_/dtDt.dof_ : 0.);
105+
h_trk_chisq_->Fill(seed.chisquared());
106+
h_trk_fitcon_->Fill(seed.fitConsistency());
107+
108+
// Fill one example event
109+
if(g_t_hit_vs_z->GetN() == 0) {
110+
for(const auto& hit : seed.hits()) {
111+
const double t_trk = hit.particleToca(); // estimated time of the particle
112+
const double t_hit = t_trk + hit.fitDt(); // add dt to get the hit time
113+
const double z_trk = hit._upoca.z();
114+
g_t_hit_vs_z ->AddPoint(z_trk, t_hit);
115+
g_t_poca_vs_z ->AddPoint(z_trk, t_trk);
116+
g_t_hit_vs_t_poca->AddPoint(t_trk, t_hit);
117+
}
118+
}
119+
}
120+
121+
//================================================================
122+
KalSeedDtDt TrackDtDt::fitTrackDtDt(const KalSeed& seed) {
123+
LsqSums2 fitter;
124+
for(const auto& hit : seed.hits()) {
125+
const double t_trk = hit.particleToca(); // estimated time of the particle
126+
const double t_hit = t_trk + hit.fitDt(); // add dt to get the hit time
127+
const double t_var = hit.fitTocaVar(); // variance of the track poca time and the hit time
128+
const double weight = 1./t_var; // inverse variance weighting
129+
fitter.addPoint(t_trk, t_hit, weight);
130+
}
131+
KalSeedDtDt result;
132+
result.slope_ = fitter.dydx();
133+
result.offset_ = fitter.y0();
134+
result.slopeUnc_ = fitter.dydxErr();
135+
result.dof_ = fitter.qn() - 2; // two parameters fitted
136+
result.chisq_ = (result.dof_ > 0) ? fitter.chi2Dof()*result.dof_ : 0.f; // chi2Dof() not defined for dof <= 0
137+
138+
if(diagLevel_ > 1) {
139+
std::cout << "[TrackDtDt::" << __func__ << "] Fit results for seed with " << seed.hits().size() << " hits:\n"
140+
<< " slope = " << result.slope_ << " +/- " << result.slopeUnc_ << "\n"
141+
<< " offset = " << result.offset_ << "\n"
142+
<< " chi2 = " << result.chisq_ << "\n"
143+
<< " dof = " << result.dof_ << "\n";
144+
}
145+
if(doHistograms_) fillHistograms(result, seed);
146+
return result;
147+
}
148+
149+
//================================================================
150+
void TrackDtDt::produce(art::Event& event) {
151+
152+
// Loop over all KalSeed collections
153+
for(const auto& name : kalSeeds_) {
154+
155+
// Retrieve the KalSeed collection from the event, checking if it is a collection of KalSeed or KalSeedPtr
156+
art::Handle<KalSeedCollection> seedHandle;
157+
event.getByLabel(name, seedHandle);
158+
art::Handle<KalSeedPtrCollection> seedPtrHandle;
159+
const bool isSeedCollection = seedHandle.isValid();
160+
if(!isSeedCollection) {
161+
event.getByLabel(name, seedPtrHandle);
162+
if(!seedPtrHandle.isValid()) {
163+
throw cet::exception("RECO") << "TrackDtDt: No KalSeed or KalSeedPtr collection with label " << name << std::endl;
164+
}
165+
}
166+
167+
// Create the output collection
168+
std::unique_ptr<KalSeedDtDtCollection> dtDtCol(new KalSeedDtDtCollection());
169+
170+
// Loop over all seeds in the collection and fit the dt_{hit} / dt_{track fit poca time} distribution
171+
const auto nseeds = (isSeedCollection) ? seedHandle->size() : seedPtrHandle->size();
172+
if(diagLevel_ > 0) std::cout << "[TrackDtDt::" << __func__ << "] Processing " << nseeds << " seeds from collection " << name << std::endl;
173+
for(size_t iseed = 0; iseed < nseeds; ++iseed) {
174+
const auto& seed = (isSeedCollection) ? seedHandle->at(iseed) : *seedPtrHandle->at(iseed);
175+
if(diagLevel_ > 1) std::cout << "[TrackDtDt::" << __func__ << "] Processing seed " << iseed << " with " << seed.hits().size() << " hits" << std::endl;
176+
dtDtCol->emplace_back(fitTrackDtDt(seed));
177+
}
178+
179+
// Put the results into the event
180+
event.put(std::move(dtDtCol), name);
181+
}
182+
}
183+
184+
//================================================================
185+
} // namespace mu2e
186+
187+
DEFINE_ART_MODULE(mu2e::TrackDtDt)

‎ParticleID/test/trkdtdt.fcl‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
# Test TrackDtDt_module.cc
2+
3+
#include "Offline/fcl/minimalMessageService.fcl"
4+
#include "Offline/fcl/standardProducers.fcl"
5+
#include "Offline/fcl/standardServices.fcl"
6+
7+
process_name : TestTrkDtDt
8+
9+
services : @local::Services.Reco
10+
services.GeometryService.bFieldFile : "Offline/Mu2eG4/geom/bfgeom_reco_v01.txt"
11+
services.scheduler.wantSummary: true
12+
13+
source: {
14+
module_type: RootInput
15+
inputCommands: [ "keep *_*_*_*"
16+
]
17+
}
18+
19+
physics : {
20+
producers : {
21+
TrackDtDt : {
22+
module_type : TrackDtDt
23+
kalSeeds : [ "KKDe" ]
24+
diagLevel : 0
25+
doHistograms: true
26+
}
27+
}
28+
RecoPath : [ TrackDtDt ]
29+
end_paths : [ ]
30+
trigger_paths : [ RecoPath ]
31+
}
32+
33+
#include "Production/Validation/database.fcl"

‎RecoDataProducts/inc/KalSeedDtDt.hh‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
//
2+
// Class representing the fit information of a dt_{straw hit time} / dt_{track hit time of closest approach}
3+
// This is useful for PID, as the wrong hypothesis will typically have a slope
4+
// Original Author: Michael MacKenzie, 2026
5+
//
6+
#ifndef RecoDataProducts_KalSeedDtDt_HH
7+
#define RecoDataProducts_KalSeedDtDt_HH
8+
#include <Rtypes.h>
9+
namespace mu2e {
10+
struct KalSeedDtDt {
11+
Float_t slope_;
12+
Float_t offset_;
13+
Float_t slopeUnc_;
14+
Float_t chisq_;
15+
Int_t dof_;
16+
Float_t slope () const { return slope_ ; }
17+
Float_t offset () const { return offset_ ; }
18+
Float_t slopeUnc () const { return slopeUnc_ ; }
19+
Float_t chisq () const { return chisq_ ; }
20+
Int_t dof () const { return dof_ ; }
21+
KalSeedDtDt() : slope_(0.f), offset_(0.f), slopeUnc_(0.f), chisq_(0.f), dof_(0) {}
22+
};
23+
24+
using KalSeedDtDtCollection = std::vector<mu2e::KalSeedDtDt>;
25+
}
26+
#endif

‎RecoDataProducts/src/classes.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@
5252
#include "Offline/RecoDataProducts/inc/KalSeed.hh"
5353
#include "Offline/RecoDataProducts/inc/KalIntersection.hh"
5454
#include "Offline/RecoDataProducts/inc/KalSeedAssns.hh"
55+
#include "Offline/RecoDataProducts/inc/KalSeedDtDt.hh"
5556
#include "Offline/RecoDataProducts/inc/TrkCaloHitPID.hh"
5657
#include "Offline/RecoDataProducts/inc/TrkQual.hh"
5758
#include "Offline/RecoDataProducts/inc/RecoQual.hh"

‎RecoDataProducts/src/classes_def.xml‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -264,6 +264,11 @@
264264
<class name="mu2e::KKLineCollection" persistent="false"/>
265265
<class name="art::Wrapper<mu2e::KKLineCollection>"/>
266266

267+
<class name="mu2e::KalSeedDtDt"/>
268+
<class name="mu2e::KalSeedDtDtCollection"/>
269+
<class name="std::vector<mu2e::KalSeedDtDt>"/>
270+
<class name="art::Wrapper<mu2e::KalSeedDtDtCollection>"/>
271+
267272
<class name="mu2e::MVAStatus"/>
268273
<class name="mu2e::TrkQual"/>
269274
<class name="mu2e::TrkQualDetail"/>

0 commit comments

Comments
 (0)