Skip to content

Move to updated PID model, add track dt/dt fit info #386

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
michaelmackenzie wants to merge 5 commits into
base: main
Choose a base branch
from
+216 −104
Open

Move to updated PID model, add track dt/dt fit info #386

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
632acb0
Move to updated PID model/module
michaelmackenzie Jul 23, 2026
5f4ff9c
Add KalSeedDtDt info to the ntuple
michaelmackenzie Jul 23, 2026
6c1b2bb
Add track dt/dt info to rooutil
michaelmackenzie Jul 23, 2026
5a23f42
Add example plot, fix field naming
michaelmackenzie Jul 23, 2026
7131c56
Add helper/doc info for track dt/dt
michaelmackenzie Jul 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 39 additions & 27 deletions fcl/prolog.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,11 @@
#include "EventNtuple/fcl/prolog_mc_reweighters.fcl"

BEGIN_PROLOG

#################################
# Track quality
#################################

# Track quality module
TrkQual : {
module_type : TrackQuality
Expand Down Expand Up @@ -50,31 +55,36 @@ TrkQualProducers : {

TrkQualProducersPath : [ TrkQualDeM, TrkQualUeM, TrkQualDmuM, TrkQualDeP, TrkQualUeP, TrkQualDmuP, TrkQualUmuM, TrkQualUmuP ]

# track PID module
TrkPID : {
module_type : TrackPID
MaxDE : 500.0 # MeV
DeltaTOffset : -1.15 # specific to MDC2018h
datFilename : "ArtAnalysis/TrkDiag/data/TrackPID_v1.dat"
}
# this module only makes sense for downstream electron fits
TrkPIDDeM : @local::TrkPID
TrkPIDDeM.KalSeedPtrCollection : "MergeKKDeM"
TrkPIDDeP : @local::TrkPID
TrkPIDDeP.KalSeedPtrCollection : "MergeKKDeP"
TrkPIDDe : @local::TrkPID
TrkPIDDe.KalSeedPtrCollection : "MergeKKDe"
TrkPIDAll : @local::TrkPID
TrkPIDAll.KalSeedPtrCollection : "MergeKKAll"
#################################
# PID
#################################

TrkPIDProducers : {
TrkPIDDeM : @local::TrkPIDDeM
TrkPIDDeP : @local::TrkPIDDeP
TrkPIDDe : @local::TrkPIDDe
TrkPIDAll : @local::TrkPIDAll
# Reconstruct track dt_hit/dt_fit
TrkDtDt : {
module_type : TrackDtDt
doHistograms : false
diagLevel : 0
kalSeeds : [ "KKDe", "KKUe", "KKDmu", "KKUmu", "MergeKKAll" ]
}

# Track PID module
TrkPID : {
module_type : TrackPID
MaxDE : 500.0 # MeV
KalSeeds : [ "KKDe" , "KKUe" , "KKDmu" , "KKUmu" , "MergeKKAll" ]
KalSeedDtDts : ["TrkDtDt:KKDe", "TrkDtDt:KKUe", "TrkDtDt:KKDmu", "TrkDtDt:KKUmu", "TrkDtDt:MergeKKAll"]
MVAVersion : 1
DatFilename : "ArtAnalysis/TrkDiag/data/TrackPID_v1.dat"
DebugLevel : 0
}
}
TrkPIDProducersPath : [ TrkPIDDeM, TrkPIDDeP ]

# trigger
TrkPIDPath : [ TrkDtDt, TrkPID ]

#################################
# Trigger
#################################

TTProducers : {
MergeTTApr : {
Expand Down Expand Up @@ -263,7 +273,8 @@ All : { input : "MergeKKAll"
fill : true
options : { fillMC : true fillHits : true genealogyDepth : -1 matchDepth : -1 }
trkQualTags : ["TrkQualAll"]
trkPIDTags : ["TrkPIDAll"]
trkPIDTags : ["TrkPID:MergeKKAll"]
trkDtDtTag : "TrkDtDt:MergeKKAll"
}
DeCalib : { input : "MergeKKDeCalib"
branchname : "trk"
Expand All @@ -277,7 +288,7 @@ De : { input : "MergeKKDe"
fill : true
options : { fillMC : true fillHits : false genealogyDepth : -1 matchDepth : -1 }
trkQualTags : ["TrkQualDe"]
trkPIDTags : ["TrkPIDDe"]
trkPIDTags : ["TrkPID:KKDe"]
}

EventNtupleMaker : {
Expand Down Expand Up @@ -322,6 +333,7 @@ EventNtupleMaker : {
fillHitCalibs : false
fillTrkQual : true
fillTrkPID : false
fillTrkDtDt : true
fits : [ @local::All ]
mc : {
fill : true
Expand Down Expand Up @@ -428,12 +440,12 @@ EventNtuple : {
genCountLogger : @local::genCountLogger
}

# TrigSequence : [ PBIWeight, @sequence::TrkQualProducersPath, @sequence::TrkPIDProducersPath]
Path : [ @sequence::MergeKKProducersPath, PBIWeight, TrkQualAll, TrkPIDAll ]
# TrigSequence : [ PBIWeight, @sequence::TrkQualProducersPath, @sequence::TrkPIDPath]
Path : [ @sequence::MergeKKProducersPath, PBIWeight, TrkQualAll, @sequence::TrkPIDPath ]
PathExt : [ @sequence::MergeKKNoFieldPath ]
PathOff : [ @sequence::MergeKKOffSpillPath ]
PathSeparate : [ @sequence::MergeKKSeparatePath, PBIWeight, @sequence::TrkQualProducersPath]
PathNoMC : [ @sequence::MergeKKProducersPath, TrkQualAll, TrkPIDAll ]
PathNoMC : [ @sequence::MergeKKProducersPath, TrkQualAll, @sequence::TrkPIDPath ]
EndPath : [ EventNtuple, genCountLogger ]
EndPathNoMC : [ EventNtuple ]

Expand Down
3 changes: 3 additions & 0 deletions inc/InfoStructHelper.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#define EventNtuple_InfoStructHelper_hh
#include "Offline/RecoDataProducts/inc/StrawHitIndex.hh"
#include "Offline/RecoDataProducts/inc/KalSeed.hh"
#include "Offline/RecoDataProducts/inc/KalSeedDtDt.hh"
#include "Offline/RecoDataProducts/inc/StrawHitFlag.hh"
#include "Offline/RecoDataProducts/inc/RecoCount.hh"
#include "Offline/RecoDataProducts/inc/HelixSeed.hh"
Expand All @@ -24,6 +25,7 @@
#include "EventNtuple/inc/LoopHelixInfo.hh"
#include "EventNtuple/inc/CentralHelixInfo.hh"
#include "EventNtuple/inc/KinematicLineInfo.hh"
#include "EventNtuple/inc/TrkDtDtInfo.hh"
#include "EventNtuple/inc/TrkStrawHitInfo.hh"
#include "EventNtuple/inc/TrkStrawHitCalibInfo.hh"
#include "EventNtuple/inc/TrkStrawMatInfo.hh"
Expand Down Expand Up @@ -67,6 +69,7 @@ namespace mu2e {
void fillKinematicLineInfo(const KalSeed& kseed, std::vector<std::vector<KinematicLineInfo>>& klis);
void fillTrkQualInfo(const KalSeed& kseed, MVAResult mva, std::vector<MVAResultInfo>& all_mvas);
void fillTrkPIDInfo(const KalSeed& kseed, MVAResult mva, std::vector<MVAResultInfo>& all_mvas);
void fillTrkDtDtInfo(const KalSeedDtDt& dtdt,std::vector<TrkDtDtInfo>& trkinfo);
void fillHitInfo(const KalSeed& kseed, std::vector<std::vector<TrkStrawHitInfo>>& all_tshinfos, std::vector<std::vector<TrkStrawHitCalibInfo>>& all_tshcinfos, bool saveCalib );
void fillMatInfo(const KalSeed& kseed, std::vector<std::vector<TrkStrawMatInfo>>& tminfos );
void fillTrkCaloHitInfo(const KalSeed& kseed, std::vector<TrkCaloHitInfo>& tchinfo );
Expand Down
20 changes: 20 additions & 0 deletions inc/TrkDtDtInfo.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
//
// TrkDtDtInfo: Result of track dt_hit / dt_trk fit
// Michael MacKenzie
//
#ifndef TrkDtDtInfo_HH
#define TrkDtDtInfo_HH

namespace mu2e
{
struct TrkDtDtInfo {
float slope = 0.; // linear fit slope
float offset = 0.; // linear fit intercept
float slopeUnc = 0.; // uncertainty on the slope
float chisq = 0.; // linear fit chi-squared
int dof = 0; // degrees of freedom in linear fit (N(hits) - 2)

void reset() { *this = TrkDtDtInfo(); }
};
}
#endif
54 changes: 54 additions & 0 deletions rooutil/examples/PlotTrkDtDt.C
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
//
// An example of how to plot the electron vs. muon track dt/dt fit results
// This uses cut functions defined in common_cuts.hh
//

#include "EventNtuple/rooutil/inc/RooUtil.hh"
#include "EventNtuple/rooutil/inc/common_cuts.hh"

#include "TH1F.h"

using namespace rooutil;
void PlotTrkDtDt(std::string filename) {

// Create the histogram you want to fill
TH1F* hE = new TH1F("h_e_dtdt", "Track dt_{hit} / dt_{trk};dt_{hit}/dt_{trk};Entries", 50,-2.,4.);
TH1F* hM = new TH1F("h_m_dtdt", "Track dt_{hit} / dt_{trk};dt_{hit}/dt_{trk};Entries", 50,-2.,4.);

// Set up RooUtil
RooUtil util(filename);

// Loop through the events
for (int i_event = 0; i_event < util.GetNEvents(); ++i_event) {
// Get the next event
auto& event = util.GetEvent(i_event);

// Get the e- and mu- tracks from the event
auto e_minus_tracks = event.GetTracks(is_e_minus );
auto mu_minus_tracks = event.GetTracks(is_mu_minus);

// Loop through the e- tracks
for (auto& track : e_minus_tracks) {
// Fill the histogram
if (track.trkdtdt) hE->Fill(track.trkdtdt->slope);
}

// Loop through the mu- tracks
for (auto& track : mu_minus_tracks) {
// Fill the histogram
if (track.trkdtdt) hM->Fill(track.trkdtdt->slope);
}
}

// Draw the histogram
hE->Draw("HIST");
hM->SetLineColor(kRed);
hM->SetMarkerStyle(20);
hM->Draw("E1 SAMES");
hE->GetYaxis()->SetRangeUser(0.1, 1.2*max(hE->GetMaximum(), hM->GetMaximum()));

auto leg = new TLegend();
leg->AddEntry(hE, "Electron fit", "L");
leg->AddEntry(hM, "Muon fit", "PE");
leg->Draw();
}
4 changes: 4 additions & 0 deletions rooutil/inc/Event.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ namespace rooutil {
AddTriggerInfo(ntuple);

CheckForBranch(ntuple, "trk", &this->trk);
CheckForBranch(ntuple, "trkdtdt", &this->trkdtdt);
CheckForBranch(ntuple, "trksegs", &this->trksegs);
CheckForBranch(ntuple, "trkcalohit", &this->trkcalohit);
CheckForBranch(ntuple, "trkqual", &this->trkqual);
Expand Down Expand Up @@ -132,6 +133,7 @@ namespace rooutil {
if (debug) { std::cout << "Event::Update(): Creating Track " << i_track << "... " << std::endl; }
Track track(&(trk->at(i_track)), &(trksegs->at(i_track)), &(trkcalohit->at(i_track))); // passing the addresses of the underlying structs
UpdateObject(track.trkmc, trkmc, i_track, debug);
UpdateObject(track.trkdtdt, trkdtdt, i_track, debug);
UpdateObject(track.trksegsmc, trksegsmc, i_track, debug);
UpdateObject(track.trksegpars_lh, trksegpars_lh, i_track, debug);
UpdateObject(track.trksegpars_ch, trksegpars_ch, i_track, debug);
Expand Down Expand Up @@ -283,6 +285,7 @@ namespace rooutil {
for (int i_trk = trks_to_remove.size()-1; i_trk >= 0; --i_trk) {
trk->erase(trk->begin()+trks_to_remove[i_trk]);
if (trkmc) { trkmc->erase(trkmc->begin()+trks_to_remove[i_trk]); }
if (trkdtdt) { trkdtdt->erase(trkdtdt->begin()+trks_to_remove[i_trk]); }
if (trksegs) { trksegs->erase(trksegs->begin()+trks_to_remove[i_trk]); }
if (trksegsmc) { trksegsmc->erase(trksegsmc->begin()+trks_to_remove[i_trk]); }
if (trkcalohit) { trkcalohit->erase(trkcalohit->begin()+trks_to_remove[i_trk]); }
Expand Down Expand Up @@ -417,6 +420,7 @@ namespace rooutil {

std::vector<mu2e::TrkInfo>* trk = nullptr;
std::vector<mu2e::TrkInfoMC>* trkmc = nullptr;
std::vector<mu2e::TrkDtDtInfo>* trkdtdt = nullptr;
std::vector<mu2e::TrkCaloHitInfo>* trkcalohit = nullptr;
std::vector<mu2e::CaloClusterInfoMC>* trkcalohitmc = nullptr;
std::vector<mu2e::MVAResultInfo>* trkqual = nullptr;
Expand Down
32 changes: 17 additions & 15 deletions rooutil/inc/Track.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,26 +12,27 @@
#include "EventNtuple/inc/TrkStrawHitInfoMC.hh"
#include "EventNtuple/inc/TrkStrawMatInfo.hh"
#include "EventNtuple/inc/TrkSegInfo.hh"
#include "EventNtuple/inc/TrkDtDtInfo.hh"
#include "EventNtuple/inc/TrkCaloHitInfo.hh"
#include "EventNtuple/inc/SimInfo.hh"
#include "EventNtuple/inc/MVAResultInfo.hh"

#include "EventNtuple/rooutil/inc/TrackSegment.hh"
#include "EventNtuple/rooutil/inc/MCParticle.hh"
#include "EventNtuple/rooutil/inc/TrackHit.hh"

namespace rooutil {
struct Track {
Track(mu2e::TrkInfo* trk, std::vector<mu2e::TrkSegInfo>* trksegs, mu2e::TrkCaloHitInfo* trkcalohit)
: trk(trk), trksegs(trksegs), trkcalohit(trkcalohit) {

// Create the underlying track segments
for (int i_segment = 0; i_segment < nSegments(); ++i_segment) {
TrackSegment segment(&(trksegs->at(i_segment))); // passing the addresses of the underlying structs
segments.emplace_back(segment);
}
}

void Update(bool debug = false) {
if (trksegsmc != nullptr) { // if we have MC information
// search for corresponding SurfaceStepInfo (it will have the same sid and sindex)
Expand All @@ -44,14 +45,14 @@ namespace rooutil {
double max_dt = 5; // want to be at least this close
for (auto& segment : segments) {
if (segment.trkseg != nullptr) { // we might have added MC-only segments...

// Want to to match segments that are at the same surface (+ index) travelling in the same direction and closest in time
double dt = std::fabs(segment.trkseg->time - trksegsmc->at(i_segment_mc).time);
if (segment.trkseg->sid == trksegsmc->at(i_segment_mc).sid
&& segment.trkseg->sindex == trksegsmc->at(i_segment_mc).sindex
&& segment.trkseg->mom.Dot(trksegsmc->at(i_segment_mc).mom)>0
&& dt < max_dt) {

if (dt < min_dt) {
min_dt = dt;
segment.trksegmc = &(trksegsmc->at(i_segment_mc));
Expand All @@ -69,7 +70,7 @@ namespace rooutil {
}
}
}

if (trkmcsim != nullptr) {
if (debug) { std::cout << "Track::Update(): Updating trkmcsim..." << std::endl; }
// Create the underlying MCParticles if possible
Expand All @@ -79,7 +80,7 @@ namespace rooutil {
mc_particles.emplace_back(mc_particle);
}
}

if (trksegpars_lh != nullptr) { // if we LoopHelix info
if (debug) { std::cout << "Track::Update(): Updating trksegpars_lh..." << std::endl; }
for (int i_segment = 0; i_segment < nSegments(); ++i_segment) {
Expand All @@ -98,20 +99,20 @@ namespace rooutil {
segments[i_segment].trksegpars_kl = &(trksegpars_kl->at(i_segment));
}
}

if (trkhits != nullptr) {
if (debug) { std::cout << "Track::Update(): Updating trkhits..." << std::endl; }
// Create the underlying TrackHits
for (int i_trkhit = 0; i_trkhit < nHits(); ++i_trkhit) {
if (debug) { std::cout << "Track::Update(): trkhit " << i_trkhit+1 << " / " << nHits() << std::endl; }
TrackHit trkhit;
trkhit.reco = &(trkhits->at(i_trkhit)); // passing the addresses of the underlying structs

// the first trkhitsmc correspond 1:1 with the reco hits, the remaining trkhitsmc also go with the truth
if (trkhitsmc != nullptr) {
trkhit.mc = &(trkhitsmc->at(i_trkhit));
}

// trkhitcalibs have 1:1 correspondance with reco hits
if (trkhitcalibs != nullptr) {
trkhit.calib = &(trkhitcalibs->at(i_trkhit));
Expand All @@ -131,7 +132,7 @@ namespace rooutil {
}
}
}

int nSegments() const { return trksegs->size(); }
TrackSegments GetSegments() { return segments; }
TrackSegments GetSegments(TrackSegmentCut cut) {
Expand All @@ -144,7 +145,7 @@ namespace rooutil {
return select_segments;
}
TrackSegments segments;

int nMCParticles() const {
if (trkmcsim == nullptr) { return 0; }
else { return trkmcsim->size(); }
Expand All @@ -160,7 +161,7 @@ namespace rooutil {
return select_mc_particles;
}
MCParticles mc_particles;

int nHits() const { return trkhits->size(); }
TrackHits GetHits() { return hits; }
TrackHits GetHits(TrackHitCut cut) {
Expand All @@ -173,10 +174,11 @@ namespace rooutil {
return select_hits;
}
TrackHits hits;

// Pointers to the data
mu2e::TrkInfo* trk = nullptr;
mu2e::TrkInfoMC* trkmc = nullptr;
mu2e::TrkDtDtInfo* trkdtdt = nullptr;
std::vector<mu2e::TrkSegInfo>* trksegs = nullptr;
std::vector<mu2e::SurfaceStepInfo>* trksegsmc = nullptr;
std::vector<mu2e::LoopHelixInfo>* trksegpars_lh = nullptr;
Expand All @@ -192,7 +194,7 @@ namespace rooutil {
mu2e::MVAResultInfo* trkqual_alt = nullptr; // TODO: is there a way to allow for more than two...
mu2e::MVAResultInfo* trkpid = nullptr;
};

typedef std::function<bool(Track&)> TrackCut;
typedef std::vector<Track> Tracks;
} // namespace rooutil
Expand Down
Loading