From 632acb01547f26b55d48c90e7d777e1ac3bc51b3 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Thu, 23 Jul 2026 10:35:09 -0500 Subject: [PATCH 1/5] Move to updated PID model/module --- fcl/prolog.fcl | 64 +++++++++++++++++++++++++++++--------------------- 1 file changed, 37 insertions(+), 27 deletions(-) diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index 4b8b68d..a98ca74 100644 --- a/fcl/prolog.fcl +++ b/fcl/prolog.fcl @@ -6,6 +6,11 @@ #include "EventNtuple/fcl/prolog_mc_reweighters.fcl" BEGIN_PROLOG + +################################# +# Track quality +################################# + # Track quality module TrkQual : { module_type : TrackQuality @@ -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 : { @@ -263,7 +273,7 @@ All : { input : "MergeKKAll" fill : true options : { fillMC : true fillHits : true genealogyDepth : -1 matchDepth : -1 } trkQualTags : ["TrkQualAll"] - trkPIDTags : ["TrkPIDAll"] + trkPIDTags : ["TrkPID:MergeKKAll"] } DeCalib : { input : "MergeKKDeCalib" branchname : "trk" @@ -277,7 +287,7 @@ De : { input : "MergeKKDe" fill : true options : { fillMC : true fillHits : false genealogyDepth : -1 matchDepth : -1 } trkQualTags : ["TrkQualDe"] - trkPIDTags : ["TrkPIDDe"] + trkPIDTags : ["TrkPID:KKDe"] } EventNtupleMaker : { @@ -428,12 +438,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 ] From 5f4ff9c0043563d75e6632d708398708de5044d5 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Thu, 23 Jul 2026 14:04:53 -0500 Subject: [PATCH 2/5] Add KalSeedDtDt info to the ntuple --- fcl/prolog.fcl | 2 ++ inc/InfoStructHelper.hh | 3 +++ inc/TrkDtDtInfo.hh | 20 ++++++++++++++++++++ src/EventNtupleMaker_module.cc | 29 +++++++++++++++++++++++++++++ src/InfoStructHelper.cc | 10 ++++++++++ src/classes.h | 1 + src/classes_def.xml | 2 ++ 7 files changed, 67 insertions(+) create mode 100644 inc/TrkDtDtInfo.hh diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index a98ca74..d81c076 100644 --- a/fcl/prolog.fcl +++ b/fcl/prolog.fcl @@ -274,6 +274,7 @@ All : { input : "MergeKKAll" options : { fillMC : true fillHits : true genealogyDepth : -1 matchDepth : -1 } trkQualTags : ["TrkQualAll"] trkPIDTags : ["TrkPID:MergeKKAll"] + trkDtDtTag : "TrkDtDt:MergeKKAll" } DeCalib : { input : "MergeKKDeCalib" branchname : "trk" @@ -332,6 +333,7 @@ EventNtupleMaker : { fillHitCalibs : false fillTrkQual : true fillTrkPID : false + fillTrkDtDt : true fits : [ @local::All ] mc : { fill : true diff --git a/inc/InfoStructHelper.hh b/inc/InfoStructHelper.hh index 3d6707b..1295b1a 100644 --- a/inc/InfoStructHelper.hh +++ b/inc/InfoStructHelper.hh @@ -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" @@ -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" @@ -67,6 +69,7 @@ namespace mu2e { void fillKinematicLineInfo(const KalSeed& kseed, std::vector>& klis); void fillTrkQualInfo(const KalSeed& kseed, MVAResult mva, std::vector& all_mvas); void fillTrkPIDInfo(const KalSeed& kseed, MVAResult mva, std::vector& all_mvas); + void fillTrkDtDtInfo(const KalSeedDtDt& dtdt,std::vector& trkinfo); void fillHitInfo(const KalSeed& kseed, std::vector>& all_tshinfos, std::vector>& all_tshcinfos, bool saveCalib ); void fillMatInfo(const KalSeed& kseed, std::vector>& tminfos ); void fillTrkCaloHitInfo(const KalSeed& kseed, std::vector& tchinfo ); diff --git a/inc/TrkDtDtInfo.hh b/inc/TrkDtDtInfo.hh new file mode 100644 index 0000000..3f60672 --- /dev/null +++ b/inc/TrkDtDtInfo.hh @@ -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 diff --git a/src/EventNtupleMaker_module.cc b/src/EventNtupleMaker_module.cc index 9499196..7552ee6 100644 --- a/src/EventNtupleMaker_module.cc +++ b/src/EventNtupleMaker_module.cc @@ -12,6 +12,7 @@ #include "Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh" #include "Offline/RecoDataProducts/inc/KalSeed.hh" #include "Offline/RecoDataProducts/inc/KalSeedAssns.hh" +#include "Offline/RecoDataProducts/inc/KalSeedDtDt.hh" #include "Offline/RecoDataProducts/inc/CaloCluster.hh" #include "Offline/RecoDataProducts/inc/CaloHit.hh" #include "Offline/RecoDataProducts/inc/CaloRecoDigi.hh" @@ -61,6 +62,7 @@ #include "EventNtuple/inc/LoopHelixInfo.hh" #include "EventNtuple/inc/CentralHelixInfo.hh" #include "EventNtuple/inc/KinematicLineInfo.hh" +#include "EventNtuple/inc/TrkDtDtInfo.hh" #include "EventNtuple/inc/SimInfo.hh" #include "EventNtuple/inc/EventWeightInfo.hh" #include "EventNtuple/inc/TrkStrawHitInfo.hh" @@ -120,6 +122,7 @@ namespace mu2e { using Comment=fhicl::Comment; fhicl::Atom input{Name("input"), Comment("KalSeedCollection input tag")}; fhicl::Atom branchname{Name("branchname"), Comment("Name of output branch")}; + fhicl::Atom trkDtDtTag{Name("trkDtDtTag"), Comment("Input tag for KalSeedDtDtCollection")}; fhicl::Atom fill{Name("fill"), Comment("Set false to skip this branch entirely (no collection reads, no output branches)")}; fhicl::Sequence trkQualTags{Name("trkQualTags"), Comment("Input tags for MVAResultCollection to use for TrkQuals")}; fhicl::Sequence trkPIDTags{Name("trkPIDTags"), Comment("Input tags for MVAResultCollection to use for TrkPID")}; @@ -151,6 +154,7 @@ namespace mu2e { fhicl::Atom fillHitCalibs{Name("fillHitCalibs"), Comment("Fill hit calibration branches")}; fhicl::Atom fillTrkQual{Name("fillTrkQual"), Comment("Fill TrkQual MVA branches")}; fhicl::Atom fillTrkPID{Name("fillTrkPID"), Comment("Fill TrkPID MVA branches")}; + fhicl::Atom fillTrkDtDt{Name("fillTrkDtDt"), Comment("Fill TrkDtDt branches")}; // per-branch configurations fhicl::Sequence> fits{Name("fits"), Comment("KalSeed collections to write into a single track branch")}; // tracker MC sub-config @@ -325,6 +329,7 @@ namespace mu2e { TrkCount _tcnt; // track branches (inputs) std::vector > _allKSPCHs; + std::vector> _allKSDtDtHs; // track branches (outputs) std::map> _allTIs; std::map>> _allTSIs; @@ -332,6 +337,7 @@ namespace mu2e { std::map>> _allCHIs; std::map>> _allKLIs; std::map> _allTCHIs; + std::map> _allTDtDtIs; // quality branches (inputs) std::vector > > _allRQCHs; std::vector >> _allTrkQualCHs; @@ -526,6 +532,7 @@ namespace mu2e { _allLHIs[i_trk_fit_branch] = {}; _allCHIs[i_trk_fit_branch] = {}; _allKLIs[i_trk_fit_branch] = {}; + _allTDtDtIs[i_trk_fit_branch] = {}; _allMCVDInfos[i_trk_fit_branch] = {}; _allTCHIs[i_trk_fit_branch] = {}; _allMCTIs[i_trk_fit_branch] = {}; @@ -594,6 +601,10 @@ namespace mu2e { if(_ftype == LoopHelix) _ntuple->Branch((branch+"segpars_lh.").c_str(),&_allLHIs.at(i_trk_fit_branch),_buffsize,_splitlevel); if(_ftype == CentralHelix) _ntuple->Branch((branch+"segpars_ch.").c_str(),&_allCHIs.at(i_trk_fit_branch),_buffsize,_splitlevel); if(_ftype == KinematicLine) _ntuple->Branch((branch+"segpars_kl.").c_str(),&_allKLIs.at(i_trk_fit_branch),_buffsize,_splitlevel); + if(_conf.trk().fillTrkDtDt()) { + _ntuple->Branch((branch+"dtdt.").c_str(),&_allTDtDtIs.at(i_trk_fit_branch),_buffsize,_splitlevel); + } + // TrkCaloHit: currently only 1 _ntuple->Branch((branch+"calohit.").c_str(),&_allTCHIs.at(i_trk_fit_branch)); for (size_t i_trkQualTag = 0; i_trkQualTag < i_trkFitConfig.trkQualTags().size(); ++i_trkQualTag) { @@ -793,6 +804,7 @@ namespace mu2e { _allBestCrvAssns.clear(); _allTrkQualCHs.clear(); _allTrkPIDCHs.clear(); + _allKSDtDtHs.clear(); art::Handle khaH; if(_conf.helices().fill()){ @@ -810,6 +822,7 @@ namespace mu2e { _allTrkQualCHs.emplace_back(); _allRQCHs.push_back({}); _allTrkPIDCHs.emplace_back(); + _allKSDtDtHs.push_back(art::Handle()); continue; } art::Handle kalSeedPtrCollHandle; @@ -817,6 +830,13 @@ namespace mu2e { event.getByLabel(kalSeedPtrInputTag,kalSeedPtrCollHandle); _allKSPCHs.push_back(kalSeedPtrCollHandle); + art::Handle kalSeedDtDtCollHandle; + if(_conf.trk().fillTrkDtDt()) { // if not filling, pass invalid handle + art::InputTag kalSeedDtDtInputTag = i_trkFitConfig.trkDtDtTag(); + event.getByLabel(kalSeedDtDtInputTag,kalSeedDtDtCollHandle); + } + _allKSDtDtHs.push_back(kalSeedDtDtCollHandle); + std::vector> trkQualCollHandles; for (const auto& i_trkQualTag : i_trkFitConfig.trkQualTags()) { art::Handle trkQualCollHandle; @@ -905,6 +925,7 @@ namespace mu2e { _allCHIs.at(i_trk_fit_branch).clear(); _allKLIs.at(i_trk_fit_branch).clear(); _allTCHIs.at(i_trk_fit_branch).clear(); + _allTDtDtIs.at(i_trk_fit_branch).clear(); _allTSHIs.at(i_trk_fit_branch).clear(); _allTSHCIs.at(i_trk_fit_branch).clear(); @@ -1244,6 +1265,14 @@ namespace mu2e { _infoStructHelper.fillMatInfo(kseed, _allTSMIs.at(i_trk_fit_branch)); } + if(_conf.trk().fillTrkDtDt()) { + auto handle = _allKSDtDtHs.at(i_trk_fit_branch); + if(handle.isValid()) { + const auto& dtdt = handle->at(i_kseedptr); + _infoStructHelper.fillTrkDtDtInfo(dtdt, _allTDtDtIs.at(i_trk_fit_branch)); + } + } + _infoStructHelper.fillTrkCaloHitInfo(kseed, _allTCHIs.at(i_trk_fit_branch)); if(kseed.hasCaloCluster()) { _tcnt._ndec = 1; diff --git a/src/InfoStructHelper.cc b/src/InfoStructHelper.cc index 14e8de5..6169e34 100644 --- a/src/InfoStructHelper.cc +++ b/src/InfoStructHelper.cc @@ -224,6 +224,16 @@ namespace mu2e { all_mvas.push_back(result_TrkPID); } + void InfoStructHelper::fillTrkDtDtInfo(const KalSeedDtDt& dtdt,std::vector& trkinfo) { + TrkDtDtInfo info; + info.slope_ = dtdt.slope_ ; + info.offset_ = dtdt.offset_ ; + info.slopeUnc_ = dtdt.slopeUnc_; + info.chisq_ = dtdt.chisq_ ; + info.dof_ = dtdt.dof_ ; + trkinfo.push_back(info); + } + void InfoStructHelper::fillTrkInfoHits(const KalSeed& kseed, TrkInfo& trkinfo) { static StrawHitFlag active(StrawHitFlag::active); std::set planes; diff --git a/src/classes.h b/src/classes.h index 95f93a4..bf8bd56 100644 --- a/src/classes.h +++ b/src/classes.h @@ -13,6 +13,7 @@ #include "EventNtuple/inc/LoopHelixInfo.hh" #include "EventNtuple/inc/CentralHelixInfo.hh" #include "EventNtuple/inc/KinematicLineInfo.hh" +#include "EventNtuple/inc/TrkDtDtInfo.hh" #include "EventNtuple/inc/TimeClusterInfo.hh" #include "EventNtuple/inc/SimInfo.hh" #include "EventNtuple/inc/MCStepInfo.hh" diff --git a/src/classes_def.xml b/src/classes_def.xml index c086ee9..472dfb3 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -22,6 +22,8 @@ + + From 6c1b2bb06fa69f0646d39e9ae580dc148e861b03 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Thu, 23 Jul 2026 14:23:30 -0500 Subject: [PATCH 3/5] Add track dt/dt info to rooutil --- rooutil/inc/Event.hh | 4 ++++ rooutil/inc/Track.hh | 32 +++++++++++++++++--------------- 2 files changed, 21 insertions(+), 15 deletions(-) diff --git a/rooutil/inc/Event.hh b/rooutil/inc/Event.hh index e8789f3..2fd6942 100644 --- a/rooutil/inc/Event.hh +++ b/rooutil/inc/Event.hh @@ -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); @@ -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); @@ -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]); } @@ -417,6 +420,7 @@ namespace rooutil { std::vector* trk = nullptr; std::vector* trkmc = nullptr; + std::vector* trkdtdt = nullptr; std::vector* trkcalohit = nullptr; std::vector* trkcalohitmc = nullptr; std::vector* trkqual = nullptr; diff --git a/rooutil/inc/Track.hh b/rooutil/inc/Track.hh index b98763e..b44567f 100644 --- a/rooutil/inc/Track.hh +++ b/rooutil/inc/Track.hh @@ -12,6 +12,7 @@ #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" @@ -19,19 +20,19 @@ #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* 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) @@ -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)); @@ -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 @@ -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) { @@ -98,7 +99,7 @@ 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 @@ -106,12 +107,12 @@ namespace rooutil { 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)); @@ -131,7 +132,7 @@ namespace rooutil { } } } - + int nSegments() const { return trksegs->size(); } TrackSegments GetSegments() { return segments; } TrackSegments GetSegments(TrackSegmentCut cut) { @@ -144,7 +145,7 @@ namespace rooutil { return select_segments; } TrackSegments segments; - + int nMCParticles() const { if (trkmcsim == nullptr) { return 0; } else { return trkmcsim->size(); } @@ -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) { @@ -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* trksegs = nullptr; std::vector* trksegsmc = nullptr; std::vector* trksegpars_lh = nullptr; @@ -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 TrackCut; typedef std::vector Tracks; } // namespace rooutil From 5a23f42ec22c4913d32dbcaa818cf30f61583fb2 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Thu, 23 Jul 2026 14:47:31 -0500 Subject: [PATCH 4/5] Add example plot, fix field naming --- inc/TrkDtDtInfo.hh | 10 +++---- rooutil/examples/PlotTrkDtDt.C | 54 ++++++++++++++++++++++++++++++++++ src/InfoStructHelper.cc | 10 +++---- 3 files changed, 64 insertions(+), 10 deletions(-) create mode 100644 rooutil/examples/PlotTrkDtDt.C diff --git a/inc/TrkDtDtInfo.hh b/inc/TrkDtDtInfo.hh index 3f60672..51c8b7e 100644 --- a/inc/TrkDtDtInfo.hh +++ b/inc/TrkDtDtInfo.hh @@ -8,11 +8,11 @@ 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) + 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(); } }; diff --git a/rooutil/examples/PlotTrkDtDt.C b/rooutil/examples/PlotTrkDtDt.C new file mode 100644 index 0000000..de44b42 --- /dev/null +++ b/rooutil/examples/PlotTrkDtDt.C @@ -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(); +} diff --git a/src/InfoStructHelper.cc b/src/InfoStructHelper.cc index 6169e34..1bd16eb 100644 --- a/src/InfoStructHelper.cc +++ b/src/InfoStructHelper.cc @@ -226,11 +226,11 @@ namespace mu2e { void InfoStructHelper::fillTrkDtDtInfo(const KalSeedDtDt& dtdt,std::vector& trkinfo) { TrkDtDtInfo info; - info.slope_ = dtdt.slope_ ; - info.offset_ = dtdt.offset_ ; - info.slopeUnc_ = dtdt.slopeUnc_; - info.chisq_ = dtdt.chisq_ ; - info.dof_ = dtdt.dof_ ; + info.slope = dtdt.slope_ ; + info.offset = dtdt.offset_ ; + info.slopeUnc = dtdt.slopeUnc_; + info.chisq = dtdt.chisq_ ; + info.dof = dtdt.dof_ ; trkinfo.push_back(info); } From 7131c56bfd1054c19eb45187e41709b33c4d06d7 Mon Sep 17 00:00:00 2001 From: michaelmackenzie Date: Thu, 23 Jul 2026 15:07:37 -0500 Subject: [PATCH 5/5] Add helper/doc info for track dt/dt --- doc/branches.md | 94 +++++++++++++++--------------------------- helper/ntuplehelper.py | 5 ++- 2 files changed, 37 insertions(+), 62 deletions(-) diff --git a/doc/branches.md b/doc/branches.md index 6ea5a77..0dd94be 100644 --- a/doc/branches.md +++ b/doc/branches.md @@ -32,10 +32,11 @@ If a Kalman fit fails or there are multiple downstream tracks to fit, the number |--------|-----------|-------------|------------------| | trk | Vector branch | information about the reconstructed track| [see TrkInfo.hh](../inc/TrkInfo.hh) | trkmc | Vector branch | MC-truth information about the track| [see TrkInfoMC.hh](../inc/TrkInfoMC.hh) -| trkcalohit | Vector branch | calorimeter cluster associated with track with POCA and track-calo matching info (`active`, `did`, `poca`, `mom`, `cdepth`, `trkdepth`, `doca`, `dt`, `tresid`, `ctime`, `ctimeerr`, `edep`, `edeperr`)| [see TrkCaloHitInfo.hh](../inc/TrkCaloHitInfo.hh) -| trkcalohitmc | Vector branch | MC-truth information for calorimeter clusters matched to track| [see CaloClusterInfoMC.hh](../inc/CaloClusterInfoMC.hh) -| trkqual | Vector branch | quality output from multi-variate analysis (MVA)| [see MVAResultInfo.hh](../inc/MVAResultInfo.hh) -| trkpid | Vector branch | particle ID output from multi-variate analysis (MVA)| [see MVAResultInfo.hh](../inc/MVAResultInfo.hh) +| trkcalohit | Vector branch | the calorimeter cluster assigned to a track| [see TrkCaloHitInfo.hh](../inc/TrkCaloHitInfo.hh) +| trkcalohitmc | Vector branch | MC-truth infromation for calorimeter clusters| [see CaloClusterInfoMC.hh](../inc/CaloClusterInfoMC.hh) +| trkdtdt | Vector branch | Result of track dt_hit / dt_trk fit| [see TrkDtDtInfo.hh](../inc/TrkDtDtInfo.hh) +| trkqual | Vector branch | the output of a multi-variate analysis (MVA)| [see MVAResultInfo.hh](../inc/MVAResultInfo.hh) +| trkpid | Vector branch | the output of a multi-variate analysis (MVA)| [see MVAResultInfo.hh](../inc/MVAResultInfo.hh) ## Track segments Branches These branches contain 4 elements per event corresponding to different Kalman fit hypotheses (see Track branches). @@ -88,49 +89,33 @@ The branch is empty if there are no time clusters during the event. | timeclusters | Vector branch | Information in a reconstructed time cluster| [see TimeClusterInfo.hh](../inc/TimeClusterInfo.hh) ## Calorimeter Branches -These branches are vectors of calorimeter clusters/hits/recodigis/digis that occurred during the event. - -Each branch can be read independently, but elements contain indexes to other branches to maintain parentage links: -- Each cluster contains a vector `hits_` with indices of hits in this cluster -- Each hit contains indices of its recodigis (vector `recoDigis_`) and parent cluster index (`clusterIdx_`) -- Each recodi contains index of its raw digi (`caloDigiIdx_`) and parent hit index (`caloHitIdx_`) -- Each raw digi contains index of parent recodi if any (`caloRecoDigiIdx_`) - -**Hierarchy:** digis --> recodigis --> hits --> clusters - -Example: cluster 3 has `hits_ = {12, 13, 14, 15}`. Each of those hits will have `clusterIdx_ = 3`. -Calorimeter hits store crystal position in (x,y,z) with frame origin at tracker center. - -| branch | structure | explanation | key fields | leaf information | -|--------|-----------|-------------|-----------|------------------| -| caloclusters | Vector branch | calorimeter clusters with indices of hits| `diskID_`, `time_`, `timeErr_`, `energyDep_`, `energyDepErr_`, `cog_` (centroid), `hits_` (vector of hit indices), `size_`, `isSplit_` | [see CaloClusterInfo.hh](../inc/CaloClusterInfo.hh) -| calohits | Vector branch | calorimeter hits with crystal positions and indices to recodigis/cluster| `crystalId_`, `nSiPMs_`, `time_`, `timeErr_`, `eDep_`, `eDepErr_`, `crystalPos_`, `recoDigis_` (vector of recodi indices), `clusterIdx_` | [see CaloHitInfo.hh](../inc/CaloHitInfo.hh) -| calorecodigis | Vector branch | fitted calorimeter recodigis with indices to raw digi and parent hit| `eDep_`, `eDepErr_`, `time_`, `timeErr_`, `chi2_`, `ndf_`, `pileUp_`, `caloDigiIdx_`, `caloHitIdx_` | [see CaloRecoDigiInfo.hh](../inc/CaloRecoDigiInfo.hh) -| calodigis | Vector branch | raw calorimeter digis with waveform data and parent recodi index| `SiPMID_`, `diskID_`, `t0_` (first sample time), `waveform_`, `peakpos_` (peak index), `posX_`, `posY_`, `caloRecoDigiIdx_` | [see CaloDigiInfo.hh](../inc/CaloDigiInfo.hh) +These branches are vectors of calorimeter clusters/hits/recodigis/digis that happened during the event. +While each branch can be read independently, i.e. all the hits of the event, each element contains indexes to the other branches for parentage link. +The cluster element contains the vector 'hits_' containing the indexes of the hits branch belonging to this cluster. +Similarly, each hit contains the indexes of its two recodigis (left and right channels) and the index of its parent cluster. +Example: cluster 3 has hits_ = {12, 13, 14, 15}. Each of those hits will have 'clusterIdx_' = 3. +Calo digis and hits have the crystal (x,y,z) position saved. Frame origin is center of tracker. +| branch | structure | explanation | leaf information | +|--------|-----------|-------------|------------------| +| caloclusters | Vector branch | calorimeter clusters with indices of hits| [see CaloClusterInfo.hh](../inc/CaloClusterInfo.hh) +| calohits | Vector branch | calorimeter hits with indices of recodigis and of parent cluster| [see CaloHitInfo.hh](../inc/CaloHitInfo.hh) +| calorecodigis | Vector branch | calorimeter recodigis with index of raw digi and of parent hit| [see CaloRecoDigiInfo.hh](../inc/CaloRecoDigiInfo.hh) +| calodigis | Vector branch | calorimeter raw digis with index of parent recodigi if any| [see CaloDigiInfo.hh](../inc/CaloDigiInfo.hh) ## Calorimeter MC Branches -The calorimeter MC branches provide MC-truth information at different levels of the calorimeter reconstruction hierarchy. - -**Alignment with Reconstruction Branches:** -Most calorimeter MC branches are index-aligned with their reco counterparts (same size and element index). One important exception: - -> **Warning — `calohitsmc` is not aligned with `calohits`.** Reco hits come from `CaloHitMaker`; MC hits come from `compressRecoMCs`, which repacks `CaloHitMC` objects when copying cluster MC truth. The two branches can differ in size and order. Offline pairs reco and MC hits via `CaloHitMCTruthAssn` at reconstruction time, but that association is not kept in production art files read by EventNtuple. Each `calohitsmc` element therefore carries `caloHitIdx_` (index into `calohits`, or `-1`). RooUtil resolves MC hits through `caloHitIdx_`; do **not** assume `calohitsmc[i]` corresponds to `calohits[i]`. +The calorimeter mc clusters branch is a vector of MC clusters associated with the reconstructed clusters. +The vector is aligned with the reconstructed clusters (same size & indexes). -- `caloclustersmc` ← aligned with `caloclusters` (same size & indices) -- `calohitsmc` ← **not** index-aligned with `calohits` (see warning above); use `caloHitIdx_` to find the reco hit -- `calodigismc` ← aligned with `calodigis` when both branches are filled (same size & indices; raw digi MC truth) -- `calodigisim` ← unique SimParticles from raw digis -- `calomcsim` ← unique SimParticles from MC clusters - -Each MC element contains `simParticleIds` to index into the corresponding `*mcsim` branch for genealogy information. - -| branch | structure | explanation | key fields | leaf information | -|--------|-----------|-------------|-----------|------------------| -| caloclustersmc | Vector branch | MC-truth information for calorimeter clusters (aligned with caloclusters)| `nsim` (# of sim particles), `etot` (total true energy), `tavg` (avg time), `eprimary` (primary particle energy), `tprimary` (primary time), `simParticleIds`, `simRels` (MCRelationship), `hits_`, `prel` (primary to event primary relationship) | [see CaloClusterInfoMC.hh](../inc/CaloClusterInfoMC.hh) -| calohitsmc | Vector branch | MC-truth for calorimeter hits (**not** index-aligned with `calohits`; use `caloHitIdx_`)| `nsim`, `eDep`, `eDepG4`, `eprimary`, `tprimary`, `eDeps`, `tDeps`, `momentumIns`, `simParticleIds`, `simRels`, `clusterIdx_`, `caloHitIdx_`, etc. | [see CaloHitInfoMC.hh](../inc/CaloHitInfoMC.hh) -| calodigismc | Vector branch | MC-truth information for raw calorimeter digis (aligned with calodigis)| `nsim`, `eDep`, `eDepG4`, `eprimary`, `tprimary`, `eDeps`, `tDeps`, `momentumIns`, `simParticleIds`, `simRels`, `caloHitIdx_`, `crystalID_`, `diskID_`, `energyCorr_`, `timeCorr_`, `posX_`, `posY_` | [see CaloDigiMCInfo.hh](../inc/CaloDigiMCInfo.hh) -| calodigisim | Vector branch | unique SimParticles in all raw digis (genealogy information)| SimParticle genealogy info | [see SimInfo.hh](../inc/SimInfo.hh) -| calomcsim | Vector branch | unique SimParticles in all MC clusters (genealogy information)| SimParticle genealogy info | [see SimInfo.hh](../inc/SimInfo.hh) +The calorimeter sim info branch is a vector of SimParticles belonging to any MC cluster. +This branch is filled according to the MC clusters order, but no sim particles are repeated (if they belong to multiple clusters) +The correct SimParticle can be retrieved from the MC cluster via the simid number +| branch | structure | explanation | leaf information | +|--------|-----------|-------------|------------------| +| caloclustersmc | Vector branch | MC-truth infromation for calorimeter clusters| [see CaloClusterInfoMC.hh](../inc/CaloClusterInfoMC.hh) +| calohitsmc | Vector branch | MC-truth infromation for calorimeter Hits| [see CaloHitInfoMC.hh](../inc/CaloHitInfoMC.hh) +| calodigismc | Vector branch | MC-truth information for calorimeter raw digis| [see CaloDigiMCInfo.hh](../inc/CaloDigiMCInfo.hh) +| calomcsim | Vector branch | information about SimParticles in genealogy| [see SimInfo.hh](../inc/SimInfo.hh) +| calodigisim | Vector branch | information about SimParticles in genealogy| [see SimInfo.hh](../inc/SimInfo.hh) ## CRV Branches These branches contain a vector where each element is a CRV hit that happened during the event. @@ -143,25 +128,14 @@ The branch is empty if there are no CRV hit during the event. | crvcoincs | Vector branch | information about a cluster of CRV coincidence triplets| [see CrvHitInfoReco.hh](../inc/CrvHitInfoReco.hh) | crvcoincsmc | Vector branch | information about the MC track which most likely caused the CRV coincidence triplets| [see CrvHitInfoMC.hh](../inc/CrvHitInfoMC.hh) | crvcoincsmcplane | Vector branch | information about the point where the MC trajectory crosses the xz plane of CRV-T| [see CrvPlaneInfoMC.hh](../inc/CrvPlaneInfoMC.hh) -| crvpulses | Vector branch | information about CRV reco pulses, including ROC/FEB/FEB channel IDs. `crvHitIndex` gives the index in `crvcoincs`, or -1 if the pulse was not clustered into a CRV hit| [see CrvPulseInfoReco.hh](../inc/CrvPulseInfoReco.hh) -| crvpulsesmc | Vector branch | MC-truth information corresponding entry-by-entry to `crvpulses`| [see CrvHitInfoMC.hh](../inc/CrvHitInfoMC.hh) +| crvpulses | Vector branch | | [see CrvPulseInfoReco.hh](../inc/CrvPulseInfoReco.hh) +| crvdigis | Vector branch | | [see CrvWaveformInfo.hh](../inc/CrvWaveformInfo.hh) +| crvpulsesmc | Vector branch | | [see CrvPulseInfoReco.hh](../inc/CrvPulseInfoReco.hh) ## Trigger Branches -Trigger branches store trigger decision information for each event and track hypothesis. Each trigger branch corresponds to a different trigger path or decision criterion. - -Trigger branches are named with the prefix `trig_` followed by the trigger name. The structure of trigger branches is a single object per event containing an array of boolean decision bits, one for each trigger name. - -**TrigInfo Structure:** -- **`_triggerArray[ntrig_]`**: A boolean array with up to 500 trigger decision bits, where `ntrig_ = 500` -- Each element is a boolean that indicates whether the trigger passed (true) or failed (false) for that event -- The index corresponds to a specific trigger path or decision algorithm - -To discover the names and meanings of all trigger branches available in your file, run: -``` -checkEventNtuple filename.root -``` +Each element in these branch corresponds to a different Kalman fit hypotheses to reconstruct the track: -This will display the mapping of indices to trigger names, allowing you to correctly interpret the trigger array. For more details about the trigger structure, see [TrigInfo.hh](../inc/TrigInfo.hh). +Trigger branches are prefixed with ```trig_```. To find the names of the trigger branches run ```checkEventNtuple filename.root``` ## Deprecated Branches These branches are remnants from trkana and are deprecated. diff --git a/helper/ntuplehelper.py b/helper/ntuplehelper.py index cfd81d7..7fdc609 100644 --- a/helper/ntuplehelper.py +++ b/helper/ntuplehelper.py @@ -3,11 +3,11 @@ class nthelper: single_object_branches = ['evtinfo', 'evtinfomc', 'hitcount', 'tcnt', 'crvsummary', 'crvsummarymc'] - vector_object_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'timeclusters', 'caloclustersmc', 'calohitsmc', 'calodigismc', 'calomcsim', 'calodigisim', 'caloclusters', 'calohits', 'calorecodigis', 'calodigis', 'crvcoincs', 'crvcoincsmc', 'crvcoincsmcplane', 'crvpulses', 'crvdigis', 'crvpulsesmc', 'trkqual', 'trkpid', 'mcsteps'] + vector_object_branches = ['trk', 'trkmc', 'trkdtdt', 'trkcalohit', 'trkcalohitmc', 'timeclusters', 'caloclustersmc', 'calohitsmc', 'calodigismc', 'calomcsim', 'calodigisim', 'caloclusters', 'calohits', 'calorecodigis', 'calodigis', 'crvcoincs', 'crvcoincsmc', 'crvcoincsmcplane', 'crvpulses', 'crvdigis', 'crvpulsesmc', 'trkqual', 'trkpid', 'mcsteps'] vector_vector_object_branches = ['trksegs', 'trksegpars_lh', 'trksegpars_ch', 'trksegpars_kl', 'trkmcsim', 'trkhits', 'trkhitsmc', 'trkmats', 'trkhitcalibs', 'trkmcsci', 'trkmcssi', 'trksegsmc' ] evt_branches = ['evtinfo','evtinfomc','hitcount','tcnt'] - trk_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'trkqual', 'trkpid'] + trk_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'trkdtdt', 'trkqual', 'trkpid'] trksegs_branches = ['trksegs', 'trksegpars_lh', 'trksegpars_ch', 'trksegpars_kl', 'trksegsmc'] straw_branches = ['trkhits', 'trkmats', 'trkhitsmc', 'trkhitcalibs'] trk_mc_branches = [ 'trkmcsim' ] @@ -67,6 +67,7 @@ class nthelper: "crvcoincsmcplane" : "CrvPlaneInfoMC", "trkqual" : "MVAResultInfo", "trkpid" : "MVAResultInfo", + "trkdtdt" : "TrkDtDtInfo", "helices" : "HelixInfo", "trksegsmc" : "SurfaceStepInfo", "mcsteps" : "MCStepInfo",