From 77102ebae970bc1637a2fcf17f74e0fb8cfe80bd Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Fri, 30 Jan 2026 12:22:51 -0600 Subject: [PATCH 01/10] Add branch for mcsteps that is independent of the track fits --- fcl/prolog.fcl | 1 + inc/InfoMCStructHelper.hh | 2 +- inc/MCStepInfo.hh | 4 +++ src/EventNtupleMaker_module.cc | 52 ++++++++++++++++++++++++++++++++-- src/InfoMCStructHelper.cc | 19 +++++++++++++ 5 files changed, 75 insertions(+), 3 deletions(-) diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index 868f28f1..b3b92d0d 100644 --- a/fcl/prolog.fcl +++ b/fcl/prolog.fcl @@ -319,6 +319,7 @@ EventNtupleMaker : { ExtraMCStepCollectionTags : [] SurfaceStepCollectionTag : "compressRecoMCs" FitType : LoopHelix + StepPointMCTags : [ "compressRecoMCs:virtualdetector" ] } # instance for processing trigger (digitization) output from simulation EventNtupleMakerTTMC: { diff --git a/inc/InfoMCStructHelper.hh b/inc/InfoMCStructHelper.hh index 63d2d238..ba555ae6 100644 --- a/inc/InfoMCStructHelper.hh +++ b/inc/InfoMCStructHelper.hh @@ -75,7 +75,7 @@ namespace mu2e { void fillExtraMCStepInfos(KalSeedMC const& kseedmc, StepPointMCCollection const& mcsteps, std::vector& mcsics, std::vector& mcssis); void fillSurfaceStepInfos(KalSeedMC const& kseedmc, SurfaceStepCollection const& surfsteps,std::vector& ssic); - + void fillStepPointMCInfo(StepPointMCCollection const& mcsteps, MCStepInfos& mcstepinfos); }; } diff --git a/inc/MCStepInfo.hh b/inc/MCStepInfo.hh index b10a420d..f0b9a7d5 100644 --- a/inc/MCStepInfo.hh +++ b/inc/MCStepInfo.hh @@ -20,6 +20,10 @@ namespace mu2e { XYZVectorF pos; // particle position at the start of this step void reset() {*this = MCStepInfo(); } bool valid() { return vid>=0; } + int pdg = -1; // true PDG code + int startCode = -1; // creation process code + int stopCode = -1; // stop process code + }; using MCStepInfos = std::vector; diff --git a/src/EventNtupleMaker_module.cc b/src/EventNtupleMaker_module.cc index e0562ec4..a9d4af71 100644 --- a/src/EventNtupleMaker_module.cc +++ b/src/EventNtupleMaker_module.cc @@ -182,9 +182,10 @@ namespace mu2e { fhicl::Atom primaryParticleTag{Name("PrimaryParticleTag"), Comment("Tag for PrimaryParticle"), art::InputTag()}; fhicl::Atom kalSeedMCTag{Name("KalSeedMCAssns"), Comment("Tag for KalSeedMCAssn"), art::InputTag()}; // extra MC - fhicl::OptionalSequence extraMCStepTags{Name("ExtraMCStepCollectionTags"), Comment("Input tags for any other StepPointMCCollections you want written out")}; + fhicl::OptionalSequence extraMCStepTags{Name("ExtraMCStepCollectionTags"), Comment("Input tags for any other StepPointMCCollections you want written out. This will only write out steps associated with the KalSeed")}; // passive elements and Virtual Detector MC information fhicl::OptionalAtom SurfaceStepsTag{Name("SurfaceStepCollectionTag"), Comment("SurfaceStep Collection")}; + fhicl::OptionalSequence stepPointMCTags{Name("StepPointMCTags"), Comment("Input tags for any other StepPointMCCollections you want written out. This will write out all steps in the collection")}; // Calo MC fhicl::Atom fillCaloMC{ Name("FillCaloMC"),Comment("Fill CaloMC information"), true}; fhicl::Atom caloClusterMCTag{Name("CaloClusterMCTag"), Comment("Tag for CaloClusterMCCollection") ,art::InputTag()}; @@ -249,7 +250,7 @@ namespace mu2e { std::map _tmap; // map between path and trigger ID. ID should come from trigger itself FIXME! // MC truth (fcl parameters) bool _fillmc; - // MC truth inputs + // MC steps (associated with KalSeed) std::vector _extraMCStepTags; std::vector> _extraMCStepCollections; std::map>> _extraMCStepInfos; @@ -258,6 +259,11 @@ namespace mu2e { art::InputTag _surfaceStepsTag; std::map>> _surfaceStepInfos; art::Handle _surfaceStepsHandle; + // MC steps (all) + std::vector _stepPointMCTags; + std::vector> _stepPointMCCollections; + std::map _stepPointMCInfos; + // art::Handle _pph; art::Handle _ksmcah; @@ -441,6 +447,11 @@ namespace mu2e { _surfaceStepInfos[i_branch] = std::vector>(); } } + if(_conf.stepPointMCTags(_stepPointMCTags)){ + for (StepCollIndex i_stepPointMCTag = 0; i_stepPointMCTag < _stepPointMCTags.size(); ++i_stepPointMCTag) { + _stepPointMCInfos[i_stepPointMCTag] = MCStepInfos(); + } + } } } @@ -566,6 +577,18 @@ namespace mu2e { } // helix info if(_conf.helices()) _ntuple->Branch("helices.",&_hinfos,_buffsize,_splitlevel); + + // all MC information + if (_fillmc) { + if(_conf.stepPointMCTags(_stepPointMCTags)){ + for(size_t icoll=0;icoll < _stepPointMCTags.size(); ++icoll) { + auto const& tag = _stepPointMCTags[icoll]; + auto inst = tag.instance(); + std::string branchname = "mcsteps_" + inst + "."; + _ntuple->Branch(branchname.c_str(),&_stepPointMCInfos[icoll],_buffsize,_splitlevel); + } + } + } } void EventNtupleMaker::beginSubRun(const art::SubRun & subrun ) { @@ -675,6 +698,16 @@ namespace mu2e { } event.getByLabel(_surfaceStepsTag,_surfaceStepsHandle); + // find extra MCStep collections + _stepPointMCCollections.clear(); + for(size_t icoll = 0; icoll < _stepPointMCTags.size(); icoll++){ + auto const& tag = _stepPointMCTags[icoll]; + art::Handle mcstepch; + event.getByLabel(tag,mcstepch); + _stepPointMCCollections.push_back(mcstepch); + } + + // loop through all track types unsigned ntrks(0); for (BranchIndex i_branch = 0; i_branch < _allBranches.size(); ++i_branch) { @@ -709,6 +742,10 @@ namespace mu2e { _extraMCStepSummaryInfos.at(i_branch).at(i_extraMCStepTag).clear(); } _surfaceStepInfos.at(i_branch).clear(); + for (StepCollIndex i_stepPointMCTag = 0; i_stepPointMCTag < _stepPointMCTags.size(); ++i_stepPointMCTag) { + _stepPointMCInfos.at(i_stepPointMCTag).clear(); + } + if(_fillcalomc) { _allMCTCHIs.at(i_branch).clear(); } @@ -875,6 +912,17 @@ namespace mu2e { } } + + // fill MC information unrelated to any reconstructed object + if (_fillmc) { + // fill MCStepCollection branch + for(size_t icoll = 0; icoll < _stepPointMCTags.size(); icoll++){ + auto const& mcsteps = *_stepPointMCCollections[icoll]; + auto& mcstepinfos = _stepPointMCInfos[icoll]; + _infoMCStructHelper.fillStepPointMCInfo(mcsteps, mcstepinfos); + } + } + // fill this row in the TTree bool fill = true; // default to fliling event if(_hastrks && ntrks == 0) { // if we require tracks (_hastrks) but we have none, don't write diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc index 3acfe722..d92ec7fa 100644 --- a/src/InfoMCStructHelper.cc +++ b/src/InfoMCStructHelper.cc @@ -434,4 +434,23 @@ namespace mu2e { // std::sort(ssic.begin(),ssic.end(),[](const auto& a, const auto& b){return a.time < b.time;}); } + void InfoMCStructHelper::fillStepPointMCInfo(StepPointMCCollection const& mcsteps, MCStepInfos& mcstepinfos) { + GeomHandle det; + MCStepInfo mcstepinfo; + for(auto const& mcstep : mcsteps) { + std::cout << mcstep << std::endl; + mcstepinfo.reset(); + mcstepinfo.vid = mcstep.volumeId(); + mcstepinfo.time = mcstep.time(); + mcstepinfo.de = mcstep.totalEDep(); + mcstepinfo.dp = mcstep.postMomentum().mag()- mcstep.momentum().mag(); + mcstepinfo.mom = mcstep.momentum(); + mcstepinfo.pos = XYZVectorF(det->toDetector(mcstep.position())); + const auto& simParticle = mcstep.simParticle(); + mcstepinfo.pdg = simParticle->pdgId(); + mcstepinfo.startCode = simParticle->creationCode(); + mcstepinfo.stopCode = simParticle->stoppingCode(); + mcstepinfos.emplace_back(mcstepinfo); + } + } } From d010f49621680876ee647313255cd08fb5dd9997 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Fri, 30 Jan 2026 12:28:34 -0600 Subject: [PATCH 02/10] Add mcsteps branch to ntuplehelper --- helper/ntuplehelper.py | 7 ++++--- inc/MCStepInfo.hh | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/helper/ntuplehelper.py b/helper/ntuplehelper.py index 4136d36b..dae02e9a 100644 --- a/helper/ntuplehelper.py +++ b/helper/ntuplehelper.py @@ -3,14 +3,14 @@ class nthelper: single_object_branches = ['evtinfo', 'evtinfomc', 'hitcount', 'tcnt', 'crvsummary', 'crvsummarymc'] - vector_object_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'caloclusters', 'calohits', 'calorecodigis', 'calodigis', 'crvcoincs', 'crvcoincsmc', 'crvcoincsmcplane', 'trkqual', 'trkpid'] + vector_object_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'caloclusters', 'calohits', 'calorecodigis', 'calodigis', 'crvcoincs', 'crvcoincsmc', 'crvcoincsmcplane', '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'] trksegs_branches = ['trksegs', 'trksegpars_lh', 'trksegpars_ch', 'trksegpars_kl', 'trksegsmc'] straw_branches = ['trkhits', 'trkmats', 'trkhitsmc', 'trkhitcalibs'] - mc_branches = ['trkmcsim'] + mc_branches = ['trkmcsim', 'mcsteps'] calo_branches = ['caloclusters', 'calohits', 'calorecodigis', 'calodigis'] crv_branches = ['crvsummary','crvsummarymc','crvcoincs','crvcoincsmc','crvcoincsmcplane'] deprecated_branches = ['trkmcsci','trkmcssi'] @@ -58,7 +58,8 @@ class nthelper: "trkqual" : "MVAResultInfo", "trkpid" : "MVAResultInfo", "helices" : "HelixInfo", - "trksegsmc" : "SurfaceStepInfo" + "trksegsmc" : "SurfaceStepInfo", + "mcsteps" : "MCStepInfo" } # diff --git a/inc/MCStepInfo.hh b/inc/MCStepInfo.hh index f0b9a7d5..dfaca703 100644 --- a/inc/MCStepInfo.hh +++ b/inc/MCStepInfo.hh @@ -14,8 +14,8 @@ namespace mu2e { float time = 0; // time of this step WRT MC primary proton (ns) float de = 0; // energy deposit through this step (MeV) float dp = 0; // momentum magnitude change throw this step (MeV/c) - bool early = false; - bool late = false; // flag if this is the earliest or latest step + bool early = false; // flag if this is the earliest step + bool late = false; // flag if this is the latest step XYZVectorF mom; // particle momentum at the start of this step XYZVectorF pos; // particle position at the start of this step void reset() {*this = MCStepInfo(); } From 87ed2cfdbc7871bea8f0a8c1fab4acf8c0d7acca Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Fri, 30 Jan 2026 12:31:47 -0600 Subject: [PATCH 03/10] In ntuplehelper differentiate between MC info that is related to a trk and general MC information --- helper/ntuplehelper.py | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/helper/ntuplehelper.py b/helper/ntuplehelper.py index dae02e9a..7e1e67d3 100644 --- a/helper/ntuplehelper.py +++ b/helper/ntuplehelper.py @@ -10,7 +10,8 @@ class nthelper: trk_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'trkqual', 'trkpid'] trksegs_branches = ['trksegs', 'trksegpars_lh', 'trksegpars_ch', 'trksegpars_kl', 'trksegsmc'] straw_branches = ['trkhits', 'trkmats', 'trkhitsmc', 'trkhitcalibs'] - mc_branches = ['trkmcsim', 'mcsteps'] + trk_mc_branches = [ 'trkmcsim' ] + general_mc_branches = [ 'mcsteps' ] calo_branches = ['caloclusters', 'calohits', 'calorecodigis', 'calodigis'] crv_branches = ['crvsummary','crvsummarymc','crvcoincs','crvcoincsmc','crvcoincsmcplane'] deprecated_branches = ['trkmcsci','trkmcssi'] @@ -171,16 +172,34 @@ def list_all_branches(self, export_to_md=False): print("| " + tokens[0] + " | " + tokens[1] + " | " + tokens[2] + "| [see " + struct_file + "](../inc/"+struct_file+")") if not export_to_md: - print("\nMonte Carlo Branches") + print("\nTrk Monte Carlo Branches") print("================") else: - print("## Monte Carlo Branches\n") + print("## Trk Monte Carlo Branches\n") print("These branches contain 4 elements per event corresponding to different Kalman fit hypotheses (see Track branches).\n") print("Within each Kalman fit element, there is a vector containing Monte Carlo truth information about the particle making the track and its parent particles.\n") print("The vector is sorted in reverse chronological order, such that the last element is the initial particle simulated in GEANT4, and each element before correspond to one of its daughter particles.\n") print("| branch | structure | explanation | leaf information |") print("|--------|-----------|-------------|------------------|") - for branch in self.mc_branches: + for branch in self.trk_mc_branches: + explanation = self.get_branch_explanation(branch) + struct = self.branch_struct_dict[branch] + struct_file = struct + ".hh"; + if not export_to_md: + print(explanation) + else: + tokens=explanation.split(":") + print("| " + tokens[0] + " | " + tokens[1] + " | " + tokens[2] + "| [see " + struct_file + "](../inc/"+struct_file+")") + + if not export_to_md: + print("\nGeneral Monte Carlo Branches") + print("================") + else: + print("## General Monte Carlo Branches\n") + print("These branches contain MC information in the event and are not required to be related to tracks") + print("| branch | structure | explanation | leaf information |") + print("|--------|-----------|-------------|------------------|") + for branch in self.general_mc_branches: explanation = self.get_branch_explanation(branch) struct = self.branch_struct_dict[branch] struct_file = struct + ".hh"; From 2f2282bdb27429fefabbec908bcfafcffca008f9 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Fri, 30 Jan 2026 12:40:04 -0600 Subject: [PATCH 04/10] Add new leaves for trkmcsci branch. Even though this is a deprecated branch, I'm adding this for completeness --- src/InfoMCStructHelper.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc index d92ec7fa..b6fd0946 100644 --- a/src/InfoMCStructHelper.cc +++ b/src/InfoMCStructHelper.cc @@ -409,6 +409,9 @@ namespace mu2e { mcsi.dp = mcstep.postMomentum().mag()- mcstep.momentum().mag(); mcsi.mom = mcstep.momentum(); mcsi.pos = spos; + mcsi.pdg = mcstep.simParticle()->pdgId(); + mcsi.startCode = mcstep.simParticle()->creationCode(); + mcsi.stopCode = mcstep.simParticle()->stoppingCode(); } } } @@ -438,7 +441,6 @@ namespace mu2e { GeomHandle det; MCStepInfo mcstepinfo; for(auto const& mcstep : mcsteps) { - std::cout << mcstep << std::endl; mcstepinfo.reset(); mcstepinfo.vid = mcstep.volumeId(); mcstepinfo.time = mcstep.time(); From 191fbebc91a713658ea1d02a8556351aece74463 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Sat, 31 Jan 2026 11:37:27 -0600 Subject: [PATCH 05/10] Add mcsteps_virtualdetector branch to RooUtil --- rooutil/README.md | 16 +++++------ rooutil/examples/PrintEvents.C | 5 ++++ rooutil/inc/Event.hh | 3 ++ rooutil/inc/RooUtil.hh | 2 ++ validation/create_val_file_rooutil.C | 43 ++++++++++++++++++++++++++++ 5 files changed, 61 insertions(+), 8 deletions(-) diff --git a/rooutil/README.md b/rooutil/README.md index e163e63a..7f600dd6 100644 --- a/rooutil/README.md +++ b/rooutil/README.md @@ -114,6 +114,7 @@ Some branches are not contained in any of the above classes: * ```crvcoincsmcplane``` * ```calohits```, ```calorecodigis```, ```calodigis``` * ```trig_``` branches +* ```mcsteps_virtualdetector``` Reach out to the developers on the #analysis-tools Slack channel if you need to have these added somewhere. @@ -224,23 +225,22 @@ Beware: there are a lot of debug messages. If a new branch is added to the EventNtuple, then the following needs to be done to so that RooUtil can access the branch: 1. In [Event.hh](inc/Event.hh) add the pointers at the bottom of the file -2. In [Event.hh](inc/Event.hh) constructor, set the branch address - - make sure to test that the branch exists -3. In [Event.hh](inc/Event.hh) add the #include to the underlying object +2. In [Event.hh](inc/Event.hh) constructor, set the branch address using the ```CheckForBranch()``` function +3. In [Event.hh](inc/Event.hh) add the #include to the underlying info struct 4. In [RooUtil.hh](inc/RooUtil.hh) add it to the ```CreateOutputEventNtuple()``` function 5. Add to validation places: - [PrintEvents.C](examples/PrintEvents.C): at least the first and last leaf in the struct - [create_val_file_rooutil.C](../../validation/create_val_file_rooutil.C) - copy contents of struct into place where histograms are defined - then copy and replace "type " with "TH1F* h_branchname_" (e.g. "float " with "TH1F* h_trkcalohit_") - - then copy in " = new TH1F("h_branchname_", "", 100,0,100); //" + - then copy in " = new TH1F("h_branchname_", "", 100,0,100); //" after the histogram name - delete line after "//" - add in leaf names to histname (and make separate x, y, z histograms for XYZVectorF leaves) - copy histogram lines into main loop - - search and replace "TH1F* " with "" - - search and replace " = new TH1F("h_" with "->Fill(" - - search and replace "", "", 100,0,100);" with ");" - - search and replace "(branchname_" with "(branchname." + - search and replace ```TH1F* ``` (note space) with nothing + - search and replace ``` = new TH1F("h_``` with ```->Fill(``` + - search and replace ```", "", 100,0,100);``` with ```);``` + - search and replace ```(branchname_``` with ```(branchname.``` in the ```Fill()``` commands - update histogram ranges 6. (Optional) If appropriate, add branches to other classes (e.g. Track.hh) and to ```Event::Update()``` - be sure to test it with an example script diff --git a/rooutil/examples/PrintEvents.C b/rooutil/examples/PrintEvents.C index a4aa22dd..e38cd9e9 100644 --- a/rooutil/examples/PrintEvents.C +++ b/rooutil/examples/PrintEvents.C @@ -263,5 +263,10 @@ void PrintEvents(std::string filename) { std::cout << pair.first << ": " << event.trigger.Fired(pair.second) << std::endl; } + // mcsteps_virtualdetector branch + if (event.mcsteps_virtualdetector != nullptr) { + for (const auto& vdstep : *(event.mcsteps_virtualdetector)) { + std::cout << "vdstep: " << vdstep.vid << ", " << vdstep.time << ", " << vdstep.mom << ", " << vdstep.pos << ", " << vdstep.pdg << ", " << vdstep.startCode << ", " << vdstep.stopCode << std::endl; + } } } diff --git a/rooutil/inc/Event.hh b/rooutil/inc/Event.hh index b400a987..349d113b 100644 --- a/rooutil/inc/Event.hh +++ b/rooutil/inc/Event.hh @@ -37,6 +37,7 @@ #include "EventNtuple/inc/MVAResultInfo.hh" #include "EventNtuple/inc/TrigInfo.hh" +#include "EventNtuple/inc/MCStepInfo.hh" #include "EventNtuple/rooutil/inc/Track.hh" #include "EventNtuple/rooutil/inc/CrvCoinc.hh" @@ -87,6 +88,7 @@ namespace rooutil { CheckForBranch(ntuple, "calorecodigis", &this->calorecodigis); CheckForBranch(ntuple, "calodigis", &this->calodigis); + CheckForBranch(ntuple, "mcsteps_virtualdetector", &this->mcsteps_virtualdetector); } // Check if a branch exists in the TChain, and optionally set its address @@ -342,6 +344,7 @@ namespace rooutil { std::vector* crvcoincsmcplane = nullptr; std::vector>* trkmcsim = nullptr; + std::vector* mcsteps_virtualdetector = nullptr; // TODO: EventNtuple could have other mcsteps_* branches but for the time being just hardcode for the virtualdetector ones }; } // namespace rooutil #endif diff --git a/rooutil/inc/RooUtil.hh b/rooutil/inc/RooUtil.hh index 2f6c9855..3d6d84af 100644 --- a/rooutil/inc/RooUtil.hh +++ b/rooutil/inc/RooUtil.hh @@ -148,6 +148,8 @@ namespace rooutil { output_ntuple->Branch(("trig_" + pair.first).c_str(), &event->triginfo._triggerArray[pair.second]); } + if (event->mcsteps_virtualdetector) { output_ntuple->Branch("mcsteps_virtualdetector", event->mcsteps_virtualdetector); } + // Write out histograms from input to output hVersion->Write(); } diff --git a/validation/create_val_file_rooutil.C b/validation/create_val_file_rooutil.C index c172bb26..9788d72e 100644 --- a/validation/create_val_file_rooutil.C +++ b/validation/create_val_file_rooutil.C @@ -505,6 +505,25 @@ void create_val_file_rooutil(std::string filename, std::string outfilename) { TH1F* h_calodigis_caloRecoDigiIdx_ = new TH1F("h_calodigis_caloRecoDigiIdx_", "", 100,0,100); TH1F* h_triginfo = new TH1F("h_triginfo", "", mu2e::TrigInfo::ntrig_,0,mu2e::TrigInfo::ntrig_); + + TH1F* h_mcsteps_virtualdetector_vid = new TH1F("h_mcsteps_virtualdetector_vid", "", 150,0,150); + TH1F* h_mcsteps_virtualdetector_sid = new TH1F("h_mcsteps_virtualdetector_sid", "", 100,0,100); + TH1F* h_mcsteps_virtualdetector_iinter = new TH1F("h_mcsteps_virtualdetector_iinter", "", 100,0,100); + TH1F* h_mcsteps_virtualdetector_time = new TH1F("h_mcsteps_virtualdetector_time", "", 100,0,2000); + TH1F* h_mcsteps_virtualdetector_de = new TH1F("h_mcsteps_virtualdetector_de", "", 100,0,200); + TH1F* h_mcsteps_virtualdetector_dp = new TH1F("h_mcsteps_virtualdetector_dp", "", 100,0,200); + TH1F* h_mcsteps_virtualdetector_early = new TH1F("h_mcsteps_virtualdetector_early", "", 2,0,2); + TH1F* h_mcsteps_virtualdetector_late = new TH1F("h_mcsteps_virtualdetector_late", "", 2,0,2); + TH1F* h_mcsteps_virtualdetector_mom_x = new TH1F("h_mcsteps_virtualdetector_mom_x", "", 100,-200,200); + TH1F* h_mcsteps_virtualdetector_mom_y = new TH1F("h_mcsteps_virtualdetector_mom_y", "", 100,-200,200); + TH1F* h_mcsteps_virtualdetector_mom_z = new TH1F("h_mcsteps_virtualdetector_mom_z", "", 320,-200,200); + TH1F* h_mcsteps_virtualdetector_pos_x = new TH1F("h_mcsteps_virtualdetector_pos_x", "", 100,-1000,1000); + TH1F* h_mcsteps_virtualdetector_pos_y = new TH1F("h_mcsteps_virtualdetector_pos_y", "", 100,-1000,1000); + TH1F* h_mcsteps_virtualdetector_pos_z = new TH1F("h_mcsteps_virtualdetector_pos_z", "", 320,-1600,1600); + TH1F* h_mcsteps_virtualdetector_pdg = new TH1F("h_mcsteps_virtualdetector_pdg", "", 100,0,100); + TH1F* h_mcsteps_virtualdetector_startCode = new TH1F("h_mcsteps_virtualdetector_startCode", "", 100,0,100); + TH1F* h_mcsteps_virtualdetector_stopCode = new TH1F("h_mcsteps_virtualdetector_stopCode", "", 100,0,100); + // Grab the first event to name the bins const auto& evt = util.GetEvent(0); for (const auto pair : evt.trigger.NameToIndexMap()) { @@ -1153,6 +1172,30 @@ void create_val_file_rooutil(std::string filename, std::string outfilename) { h_triginfo->AddBinContent(i_trig+1, event.triginfo._triggerArray[i_trig]); } + // Creating mcsteps_virtualdetector histogram + std::cout << "Creating mcsteps_virtualdetector histogram" << std::endl; + if (event.mcsteps_virtualdetector != nullptr) { + for (const auto& vdstep : *(event.mcsteps_virtualdetector)) { + h_mcsteps_virtualdetector_vid->Fill(mcsteps_virtualdetector.vid); + h_mcsteps_virtualdetector_sid->Fill(mcsteps_virtualdetector.sid); + h_mcsteps_virtualdetector_iinter->Fill(mcsteps_virtualdetector.iinter); + h_mcsteps_virtualdetector_time->Fill(mcsteps_virtualdetector.time); + h_mcsteps_virtualdetector_de->Fill(mcsteps_virtualdetector.de); + h_mcsteps_virtualdetector_dp->Fill(mcsteps_virtualdetector.dp); + h_mcsteps_virtualdetector_early->Fill(mcsteps_virtualdetector.early); + h_mcsteps_virtualdetector_late->Fill(mcsteps_virtualdetector.late); + h_mcsteps_virtualdetector_mom_x->Fill(mcsteps_virtualdetector.mom.x()); + h_mcsteps_virtualdetector_mom_y->Fill(mcsteps_virtualdetector.mom.y()); + h_mcsteps_virtualdetector_mom_z->Fill(mcsteps_virtualdetector.mom.z()); + h_mcsteps_virtualdetector_pos_x->Fill(mcsteps_virtualdetector.pos.x()); + h_mcsteps_virtualdetector_pos_y->Fill(mcsteps_virtualdetector.pos.y()); + h_mcsteps_virtualdetector_pos_z->Fill(mcsteps_virtualdetector.pos.z()); + h_mcsteps_virtualdetector_pdg->Fill(mcsteps_virtualdetector.pdg); + h_mcsteps_virtualdetector_startCode->Fill(mcsteps_virtualdetector.startCode); + h_mcsteps_virtualdetector_stopCode->Fill(mcsteps_virtualdetector.stopCode); + + } + } } file->Write(); From 4213c635b35b953f0bbf32191c925966c99110d9 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Sat, 31 Jan 2026 11:39:49 -0600 Subject: [PATCH 06/10] Update developer instructions --- doc/developers.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/developers.md b/doc/developers.md index 8ec1f94b..3ced93b2 100644 --- a/doc/developers.md +++ b/doc/developers.md @@ -62,8 +62,9 @@ mu2e -c EventNtuple/fcl/from_mcs-mockdata.fcl -s test-art-file.art -n 10 8. Add to ```ntuplehelper``` * instructions [here](./ntuplehelper.md#Adding-a-branch) -9. Add to validation, RooUtil, and PyUtil - * [instructions to be completed... for the time being contact Andy and Sophie for this] +9. Add to validation, RooUtil, and pyutils + * instructions for validation and RooUtil are [here](../rooutil/README.md#For-Developers). Contact Andy for more info if needed + * contact Sophie for pyutils 10. Run validation * both [test_fcls.sh](../validation/README.md#Validating-EventNtuple-Runs), and * and [valCompare](../validation/README.md#Validating-EventNtuple-Contents) steps From 24f67246446492adc10acd8acfb773166b07e1bc Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Sat, 31 Jan 2026 12:09:47 -0600 Subject: [PATCH 07/10] Fix validation macro --- validation/create_val_file_rooutil.C | 34 ++++++++++++++-------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/validation/create_val_file_rooutil.C b/validation/create_val_file_rooutil.C index 9788d72e..87a62a25 100644 --- a/validation/create_val_file_rooutil.C +++ b/validation/create_val_file_rooutil.C @@ -1176,23 +1176,23 @@ void create_val_file_rooutil(std::string filename, std::string outfilename) { std::cout << "Creating mcsteps_virtualdetector histogram" << std::endl; if (event.mcsteps_virtualdetector != nullptr) { for (const auto& vdstep : *(event.mcsteps_virtualdetector)) { - h_mcsteps_virtualdetector_vid->Fill(mcsteps_virtualdetector.vid); - h_mcsteps_virtualdetector_sid->Fill(mcsteps_virtualdetector.sid); - h_mcsteps_virtualdetector_iinter->Fill(mcsteps_virtualdetector.iinter); - h_mcsteps_virtualdetector_time->Fill(mcsteps_virtualdetector.time); - h_mcsteps_virtualdetector_de->Fill(mcsteps_virtualdetector.de); - h_mcsteps_virtualdetector_dp->Fill(mcsteps_virtualdetector.dp); - h_mcsteps_virtualdetector_early->Fill(mcsteps_virtualdetector.early); - h_mcsteps_virtualdetector_late->Fill(mcsteps_virtualdetector.late); - h_mcsteps_virtualdetector_mom_x->Fill(mcsteps_virtualdetector.mom.x()); - h_mcsteps_virtualdetector_mom_y->Fill(mcsteps_virtualdetector.mom.y()); - h_mcsteps_virtualdetector_mom_z->Fill(mcsteps_virtualdetector.mom.z()); - h_mcsteps_virtualdetector_pos_x->Fill(mcsteps_virtualdetector.pos.x()); - h_mcsteps_virtualdetector_pos_y->Fill(mcsteps_virtualdetector.pos.y()); - h_mcsteps_virtualdetector_pos_z->Fill(mcsteps_virtualdetector.pos.z()); - h_mcsteps_virtualdetector_pdg->Fill(mcsteps_virtualdetector.pdg); - h_mcsteps_virtualdetector_startCode->Fill(mcsteps_virtualdetector.startCode); - h_mcsteps_virtualdetector_stopCode->Fill(mcsteps_virtualdetector.stopCode); + h_mcsteps_virtualdetector_vid->Fill(vdstep.vid); + h_mcsteps_virtualdetector_sid->Fill(vdstep.sid); + h_mcsteps_virtualdetector_iinter->Fill(vdstep.iinter); + h_mcsteps_virtualdetector_time->Fill(vdstep.time); + h_mcsteps_virtualdetector_de->Fill(vdstep.de); + h_mcsteps_virtualdetector_dp->Fill(vdstep.dp); + h_mcsteps_virtualdetector_early->Fill(vdstep.early); + h_mcsteps_virtualdetector_late->Fill(vdstep.late); + h_mcsteps_virtualdetector_mom_x->Fill(vdstep.mom.x()); + h_mcsteps_virtualdetector_mom_y->Fill(vdstep.mom.y()); + h_mcsteps_virtualdetector_mom_z->Fill(vdstep.mom.z()); + h_mcsteps_virtualdetector_pos_x->Fill(vdstep.pos.x()); + h_mcsteps_virtualdetector_pos_y->Fill(vdstep.pos.y()); + h_mcsteps_virtualdetector_pos_z->Fill(vdstep.pos.z()); + h_mcsteps_virtualdetector_pdg->Fill(vdstep.pdg); + h_mcsteps_virtualdetector_startCode->Fill(vdstep.startCode); + h_mcsteps_virtualdetector_stopCode->Fill(vdstep.stopCode); } } From e5ec8eb3a49ab0be9fac244627ef2bc7178a04f0 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Sat, 31 Jan 2026 12:10:26 -0600 Subject: [PATCH 08/10] Add example for plotting VD steps with RooUtil --- rooutil/examples/PlotVDSteps.C | 52 ++++++++++++++++++++++++++++++++++ rooutil/inc/common_cuts.hh | 1 + 2 files changed, 53 insertions(+) create mode 100644 rooutil/examples/PlotVDSteps.C diff --git a/rooutil/examples/PlotVDSteps.C b/rooutil/examples/PlotVDSteps.C new file mode 100644 index 00000000..bd51978e --- /dev/null +++ b/rooutil/examples/PlotVDSteps.C @@ -0,0 +1,52 @@ +// +// An example of how to plot reco vs true +// This uses cut functions defined in common_cuts.hh +// + +#include "EventNtuple/rooutil/inc/RooUtil.hh" +#include "EventNtuple/rooutil/inc/common_cuts.hh" + +#include "TCanvas.h" +#include "TH1F.h" +#include "TH2F.h" + +using namespace rooutil; +void PlotVDSteps(std::string filename) { + + // Create the histogram you want to fill + TH2F* hVD13_XY = new TH2F("hVD13_XY", "XY position at VD13 (tracker entrance)", 100,-500,500, 100,-500,500); + hVD13_XY->SetXTitle("X (det. coords) [mm]"); + hVD13_XY->SetYTitle("Y (det. coords) [mm]"); + + TH1F* hVD13_Mom = new TH1F("hVD13_Mom", "Momentum at VD13 (tracker entrance)", 120,0,120); + hVD13_Mom->SetXTitle("Momentum [MeV/c]"); + + // 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_minus tracks from the event + if (event.mcsteps_virtualdetector != nullptr) { + for (const auto& vdstep : *(event.mcsteps_virtualdetector)) { + if (vdstep.vid == mu2e::VirtualDetectorId::TT_FrontHollow) { + // Fill the histogram + hVD13_XY->Fill(vdstep.pos.x(), vdstep.pos.y()); + + hVD13_Mom->Fill(vdstep.mom.R()); + } + } + } + } + + // Draw the histogram + TCanvas* c = new TCanvas(); + c->Divide(2); + c->cd(1); + hVD13_XY->Draw("COLZ"); + c->cd(2); + hVD13_Mom->Draw("HIST E"); +} diff --git a/rooutil/inc/common_cuts.hh b/rooutil/inc/common_cuts.hh index ef22d70b..5f3928b2 100644 --- a/rooutil/inc/common_cuts.hh +++ b/rooutil/inc/common_cuts.hh @@ -10,6 +10,7 @@ #include "Offline/MCDataProducts/inc/ProcessCode.hh" #include "Offline/MCDataProducts/inc/GenId.hh" #include "Offline/MCDataProducts/inc/MCRelationship.hh" +#include "Offline/DataProducts/inc/VirtualDetectorId.hh" #include "EventNtuple/rooutil/inc/Track.hh" #include "EventNtuple/rooutil/inc/TrackSegment.hh" From a125640374dfffe276cc64c19d0124a318723118 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Sat, 31 Jan 2026 12:46:31 -0600 Subject: [PATCH 09/10] Update RooUtil validation --- rooutil/examples/PrintEvents.C | 1 + validation/test_rooutil_examples.sh | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/rooutil/examples/PrintEvents.C b/rooutil/examples/PrintEvents.C index e38cd9e9..07d70321 100644 --- a/rooutil/examples/PrintEvents.C +++ b/rooutil/examples/PrintEvents.C @@ -268,5 +268,6 @@ void PrintEvents(std::string filename) { for (const auto& vdstep : *(event.mcsteps_virtualdetector)) { std::cout << "vdstep: " << vdstep.vid << ", " << vdstep.time << ", " << vdstep.mom << ", " << vdstep.pos << ", " << vdstep.pdg << ", " << vdstep.startCode << ", " << vdstep.stopCode << std::endl; } + } } } diff --git a/validation/test_rooutil_examples.sh b/validation/test_rooutil_examples.sh index 852f6ffb..eb74676c 100644 --- a/validation/test_rooutil_examples.sh +++ b/validation/test_rooutil_examples.sh @@ -1,10 +1,10 @@ -filename="\"/pnfs/mu2e/tape/phy-nts/nts/mu2e/CeEndpointOnSpillTriggered/MDC2020aq_best_v1_3_v06_03_00/root/ed/df/nts.mu2e.CeEndpointOnSpillTriggered.MDC2020aq_best_v1_3_v06_03_00.001210_00000000.root\"" +filename="\"nts.ntuple.after.root\"" scripts=( "PrintEvents.C" "CreateNtuple.C" "CreateTrackNtuple.C" "PlotCRVPEs.C" "PlotCRVPEsVsMCEDep.C" "PlotEntranceFitPars.C" "PlotEntranceMomentum.C" "PlotEntranceMomentumCRVCut.C" "PlotEntranceMomentumResolution.C" "PlotEntranceMomentumResolution_TrkQualCut.C" "PlotMCParentPosZ.C" "PlotMCParticleMom.C" "PlotMuonPosZ.C" "PlotStoppingTargetFoilSegment.C" "PlotTrackNHits_RecoVsTrue.C" "PlotTrkCaloHitEnergy.C" "PrintEventsNoMC.C" "TrackCounting.C" "PlotTrackHitTimes.C" "PlotTrackHitTimesMC.C" "PlotStrawMaterials.C" - "PlotGenCosmicMom.C" "PlotCRVTotalPEs.C" "PlotEntranceMomentum_UpstreamDownstream.C" ) + "PlotGenCosmicMom.C" "PlotCRVTotalPEs.C" "PlotEntranceMomentum_UpstreamDownstream.C" "PlotVDSteps.C" ) for script in "${scripts[@]}" do From b93fd85d4002cdefc9b9808006016ad8aff230d4 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Sat, 31 Jan 2026 13:31:18 -0600 Subject: [PATCH 10/10] Turn off vd steps by default but add a fcl file that demonstrates how to use it --- fcl/README.md | 1 + fcl/from_mcs-primary_addVDSteps.fcl | 3 +++ fcl/prolog.fcl | 2 +- validation/test_fcls.sh | 14 ++++++++++++++ 4 files changed, 19 insertions(+), 1 deletion(-) create mode 100644 fcl/from_mcs-primary_addVDSteps.fcl diff --git a/fcl/README.md b/fcl/README.md index dc1838b5..9fe79829 100644 --- a/fcl/README.md +++ b/fcl/README.md @@ -22,6 +22,7 @@ where ```tier``` is the data tier of the input dataset, ```type``` is the type o | from_mcs-mockdata_separateTrkBranches.fcl | mock datasets | example on how to separate the tracks into separate branches again| | from_mcs-mockdata_selectorExample.fcl | mock datasets | example on how to use a selector to select certain types of tracks before putting them into the EventNtuple | | from_mcs-mixed_trkQualCompare.fcl | reconstructed mixed (i.e. primary+background hits) datasets | shows how to output result of more than one TrkQual | +| from_mcs-primary_addVDSteps.fcl | reconstructed primary (i.e. no background hits) datasets | shows how to add the branch for virtual detector steps | | from_mcs-DeMCalib.fcl | reconstructed primary or mixed datasets | only writes one track per event | | from_mcs-OffSpill.fcl | off spill datasets | only contains ```CentralHelix``` tracks (i.e. field-on cosmics) | | from_dig-mockdata.fcl | mock datasets (digis) | runs reconstruction and creates EventNtuple in one job | diff --git a/fcl/from_mcs-primary_addVDSteps.fcl b/fcl/from_mcs-primary_addVDSteps.fcl new file mode 100644 index 00000000..85906f29 --- /dev/null +++ b/fcl/from_mcs-primary_addVDSteps.fcl @@ -0,0 +1,3 @@ +#include "EventNtuple/fcl/from_mcs-primary.fcl" + +physics.analyzers.EventNtuple.StepPointMCTags : [ "compressRecoMCs:virtualdetector" ] diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index b3b92d0d..b18595c8 100644 --- a/fcl/prolog.fcl +++ b/fcl/prolog.fcl @@ -319,7 +319,7 @@ EventNtupleMaker : { ExtraMCStepCollectionTags : [] SurfaceStepCollectionTag : "compressRecoMCs" FitType : LoopHelix - StepPointMCTags : [ "compressRecoMCs:virtualdetector" ] + StepPointMCTags : [ ] } # instance for processing trigger (digitization) output from simulation EventNtupleMakerTTMC: { diff --git a/validation/test_fcls.sh b/validation/test_fcls.sh index 6f05f89a..6792ad42 100644 --- a/validation/test_fcls.sh +++ b/validation/test_fcls.sh @@ -2,6 +2,10 @@ # test_fcls.sh - runs fcl files to make sure that they complete successfully # Note: requires relevant filelists in a directory above this one # +if [ $# -ne 1 ]; then + echo "Usage: . ./validation/test_fcls.sh MDC202X" + return +fi vMDC=$1 log_file="test_fcls.log" @@ -137,6 +141,16 @@ else echo "FAIL" fi +echo -n "from_mcs-primary_addVDSteps.fcl... " +echo "mu2e -c fcl/from_mcs-primary_addVDSteps.fcl -S ../filelists/$primary_dataset.list --TFileName nts.ntuple.primaryVDSteps.root -n 100" >> ${log_file} 2>&1 +mu2e -c fcl/from_mcs-primary_addVDSteps.fcl -S ../filelists/$primary_dataset.list --TFileName nts.ntuple.primaryVDSteps.root -n 100 >> ${log_file} 2>&1 +if [ $? == 0 ]; then + echo "OK" +else + echo "FAIL" +fi + + echo -n "from_dig-mockdata.fcl... " echo "mu2e -c fcl/from_dig-mockdata.fcl -S ../filelists/${digi_dataset}.list --TFileName nts.ntuple.dig.root -n 100" >> ${log_file} 2>&1 mu2e -c fcl/from_dig-mockdata.fcl -S ../filelists/${digi_dataset}.list --TFileName nts.ntuple.dig.root -n 100 >> ${log_file} 2>&1