Subrun ntuple

README.md

@@ -24,14 +24,16 @@ A [list of branches is available](./doc/branches.md)
 
 The help understand what all the branches and leaves mean, we have an [```ntuplehelper```](doc/ntuplehelper.md) tool
 
+Alongside the EventNtuple itself we have a SubRunNtuple that contains subrun-level quantities. These quantites are also summed into the "totals histograms".
+
 ## How to Check the EventNtuple
 You can get some basic information about the EventNtuple file useing ```checkEventNtuple``` like so:
 
 ```
 checkEventNtuple file1.root file2.root
 ```
 
-This will print the version number (for versions after v6.3.0) and the trigger branches in the file (for versions after v6.8.0)
+This will print the version number (for versions after v6.3.0) and the trigger branches in the file (for versions after v6.8.0). It will also print the number of entries in the "totals histograms" if any exist (for versions after v6.13.0)
 
 
 ## How to Analyze an EventNtuple

bin/checkEventNtuple

@@ -38,6 +38,19 @@ def main():
                     print("For trigger definitions, please contact the Trigger group.")
                 else:
                     print("{} tree does not exist in {}...".format(treename, filename))
+
+                # Print totals histograms
+                print("\nTotals histograms:")
+                totals_histnames = ["n_gen_events", "n_proc_events", "cosmic_livetime"]
+                any_found = False
+                for totals_histname in totals_histnames:
+                    if (folder.GetListOfKeys().Contains(totals_histname)):
+                        h = f.Get(foldername+"/"+totals_histname)
+                        print(" - {} contains {:.02f} entries".format(totals_histname, h.GetBinContent(1)))
+                        any_found = True
+                if not any_found:
+                    print(" - no totals histograms (e.g. n_proc_events) found")
+
             else:
                 print("{}/ folder does not exist in {}".format(foldername, filename))
         except OSError:

fcl/from_dig-OnSpill.fcl

@@ -56,3 +56,7 @@ physics.analyzers.EventNtupleTTMCMpr.trk.mc.kalSeedMCAssns : "TTMprDeKSFMC"
 end_paths : [ EndPath ]
 physics.trigger_paths : [ "TrigPath" ]
 services.TFileService.fileName: "nts.owner.trkana-triggerMC.version.sequencer.root"
+physics.analyzers.EventNtupleTTMCApr.subruns.cosmicLivetime.include : false
+physics.analyzers.EventNtupleTTMCTpr.subruns.cosmicLivetime.include : false
+physics.analyzers.EventNtupleTTMCCpr.subruns.cosmicLivetime.include : false
+physics.analyzers.EventNtupleTTMCMpr.subruns.cosmicLivetime.include : false

fcl/from_dig-calo.fcl

@@ -54,3 +54,5 @@ physics.analyzers.EventNtuple.calo.mc.fillDigiSim : true
 physics.analyzers.EventNtuple.calo.mc.fillDigis : true
 
 services.TFileService.fileName: "nts.owner.description.version.sequencer.root"
+
+physics.analyzers.EventNtuple.subruns.cosmicLivetime.include : false

fcl/from_dig-mockdata.fcl

@@ -36,7 +36,7 @@ physics.EventNtuplePath : [
   @sequence::Reconstruction.OnSpillPath,
   @sequence::EventNtuple.Path
 ]
-physics.EventNtupleEndPath : [ @sequence::EventNtuple.EndPathNoMC ]
+physics.EventNtupleEndPath : [ @sequence::EventNtuple.EndPath ]
 physics.analyzers.EventNtuple.SelectEvents : [ "EventNtuplePath" ]
 
 physics.trigger_paths : [ EventNtuplePath ]
@@ -54,3 +54,5 @@ services.TFileService.fileName: "nts.owner.description.version.sequencer.root"
 #include "Production/JobConfig/reco/drop_trigger.fcl"
 
 #include "Production/Validation/database.fcl"
+
+physics.analyzers.EventNtuple.subruns.cosmicLivetime.include : false

fcl/from_mcs-DeCalib.fcl

@@ -24,3 +24,4 @@ physics.end_paths : [ EndPath ]
 physics.producers.TrkQualDe.KalSeedPtrCollection : "MergeKKDeCalib"
 services.TimeTracker.printSummary: true
 services.TFileService.fileName: "nts.owner.EventNtupleDeCalib.version.sequence.root"
+physics.analyzers.EventNtuple.subruns.cosmicLivetime.include : false

fcl/from_mcs-ceSimReco.fcl

@@ -4,4 +4,3 @@ physics.analyzers.EventNtuple.trk.fits : [ @local::De ]
 physics.analyzers.EventNtuple.trk.fits[0].branchname : "trk"
 physics.analyzers.EventNtuple.FillTriggerInfo : false
 physics.EventNtuplePath : [ MergeKKDe, PBIWeight, TrkQualDe, TrkPIDDe ]
-physics.EventNtupleEndPath : [ EventNtuple, genCountLogger ]

fcl/from_mcs-extracted.fcl

@@ -4,5 +4,7 @@ physics.EventNtuplePath : [ @sequence::EventNtuple.PathExt ] # path for extracte
 physics.analyzers.EventNtuple.trk.fits : [ @local::Ext ]
 physics.analyzers.EventNtuple.FitType : KinematicLine
 physics.analyzers.EventNtuple.trk.fillTrkQual : false
+physics.analyzers.EventNtuple.subruns.genEventCount.include : true
+physics.analyzers.EventNtuple.subruns.cosmicLivetime.include : true
 
 services.GeometryService.inputFile: "Production/JobConfig/cosmic/geom_cosmic_extracted.txt"

fcl/from_mcs-mockdata.fcl

@@ -16,11 +16,13 @@ physics :
 }
 
 physics.EventNtuplePath : [ @sequence::EventNtuple.Path ]
-physics.EventNtupleEndPath : [ @sequence::EventNtuple.EndPathNoMC ]
+physics.EventNtupleEndPath : [ @sequence::EventNtuple.EndPath ]
 
 physics.trigger_paths : [ EventNtuplePath ]
 physics.end_paths : [ EventNtupleEndPath ]
 
 #include "Production/JobConfig/common/epilog.fcl"
 services.TFileService.fileName: "nts.owner.description.version.sequencer.root"
 
+physics.analyzers.EventNtuple.subruns.genEventCount.include : false
+physics.analyzers.EventNtuple.subruns.cosmicLivetime.include : false

fcl/from_mcs-mockdata_noMC.fcl

@@ -4,4 +4,3 @@ physics.analyzers.EventNtuple.mc.fill : false
 physics.analyzers.EventNtuple.calo.mc.fill : false
 
 physics.EventNtuplePath : [ @sequence::EventNtuple.PathNoMC ]
-physics.EventNtupleEndPath : [ @sequence::EventNtuple.EndPathNoMC ]

fcl/from_mcs-primary.fcl

@@ -1,3 +1,3 @@
 #include "EventNtuple/fcl/from_mcs-mockdata.fcl"
 
-physics.EventNtupleEndPath : [ @sequence::EventNtuple.EndPath ] # adds back genCountLogger
+physics.analyzers.EventNtuple.subruns.genEventCount.include : true

fcl/prolog.fcl

@@ -389,6 +389,15 @@ EventNtupleMaker : {
   mcsteps : {
     surfaceStepsTag : "compressRecoMCs"
   }
+
+  # -- SubRun quantities ---------------------------------------------------
+  subruns : {
+    makeNtuple : true
+    makeTotalsHistograms : true
+    genEventCount : { include : true inputTag : "genCounter" fillNtuple : true fillTotalsHistogram : true}
+    procEventCount : { include : true inputTag : "" fillNtuple : true fillTotalsHistogram : true }
+    cosmicLivetime : { include : true inputTag : "generate" fillNtuple : true fillTotalsHistogram : true }
+  }
 }
 # instance for processing trigger (digitization) output from simulation
 EventNtupleMakerTTMC : {
@@ -425,7 +434,6 @@ EventNtuple : {
 
   analyzers : {
     EventNtuple : { @table::EventNtupleMaker  }
-    genCountLogger : @local::genCountLogger
   }
 
   #  TrigSequence : [ PBIWeight, @sequence::TrkQualProducersPath, @sequence::TrkPIDProducersPath]
@@ -434,8 +442,7 @@ EventNtuple : {
   PathOff : [ @sequence::MergeKKOffSpillPath ]
   PathSeparate : [ @sequence::MergeKKSeparatePath, PBIWeight, @sequence::TrkQualProducersPath]
   PathNoMC : [ @sequence::MergeKKProducersPath, TrkQualAll, TrkPIDAll ]
-  EndPath : [ EventNtuple, genCountLogger ]
-  EndPathNoMC : [ EventNtuple ]
+  EndPath : [ EventNtuple ]
 
 }
 

inc/SubRunInfo.hh

@@ -0,0 +1,16 @@
+//
+// SubRunInfo: subrun information
+// Andy Edmonds (2026)
+//
+#ifndef SubRunInfo_HH
+#define SubRunInfo_HH
+#include <string>
+namespace mu2e
+{
+  struct SubRunInfo {
+    int run = 0; // run number
+    int subrun = 0; // subrun number
+    void reset() {*this = SubRunInfo(); }
+  };
+}
+#endif

rooutil/README.md

@@ -3,17 +3,18 @@
 ## Table of Contents
 1. [Introduction](#Introduction)
 2. [```RooUtil``` Class](#RooUtil-Class)
-3. [The ```Event``` Class](#The-Event-Class)
-4. [User-Friendly Classes](#User-Friendly-Classes)
-5. [Accessing User-Friendly Classes](#Accessing-User-Friendly-Classes)
-6. [Cut Functions](#Cut-Functions)
-7. [Common Cut Functions](#Common-Cut-Functions)
-8. [Combining Cut Function](#Combining-Cut-Functions)
-9. [Creating Ntuples From EventNtuple](#Creating-Ntuples-From-EventNtuple)
-10. [```roodask```](#roodask)
-11. [Speed Optimizations](#Speed-Optimizations)
-12. [Debugging](#Debugging)
-13. [For Developers](#For-Developers)
+3. [SubRun Data and Totals](#SubRun-Data-and-Totals)
+4. [The ```Event``` Class](#The-Event-Class)
+5. [User-Friendly Classes](#User-Friendly-Classes)
+6. [Accessing User-Friendly Classes](#Accessing-User-Friendly-Classes)
+7. [Cut Functions](#Cut-Functions)
+8. [Common Cut Functions](#Common-Cut-Functions)
+9. [Combining Cut Function](#Combining-Cut-Functions)
+10. [Creating Ntuples From EventNtuple](#Creating-Ntuples-From-EventNtuple)
+11. [```roodask```](#roodask)
+12. [Speed Optimizations](#Speed-Optimizations)
+13. [Debugging](#Debugging)
+14. [For Developers](#For-Developers)
 
 ## Introduction
 
@@ -29,6 +30,24 @@ The constructor takes two arguments:
 * ```filename``` can be the name of a single ROOT file (ending in ```.root```) or a list of ROOT files
 * (optional) ```treename``` the name of the tree
 
+## SubRun Data and Totals
+
+RooUtil loads the optional ```EventNtuple/subrunNtuple``` independently of the event tree. Each entry is represented by a ```SubRun``` object with the ```srinfo``` run/subrun identifier and any configured ```genEventCount```, ```procEventCount```, and ```cosmicLivetime``` quantities.
+
+```
+RooUtil util(filename);
+for (int i_subrun = 0; i_subrun < util.GetNSubRuns(); ++i_subrun) {
+  const auto& subrun = util.GetSubRun(i_subrun);
+  std::cout << subrun.srinfo->run << ":" << subrun.srinfo->subrun << std::endl;
+}
+```
+
+Use ```FindSubRun(run, subrun)``` to retrieve the SubRun corresponding to an event. It returns ```nullptr``` when the requested run/subrun is unavailable. ```GetSubRunIndices(cut)``` returns the indices of entries passing a function with signature ```bool cut(const SubRun&)```. ```CreateOutputSubRunNtuple()``` and ```FillOutputSubRunNtuple()``` create a reduced SubRun tree.
+
+The one-bin totals histograms are summed across every input file. ```GetTotal(name)``` returns ```std::optional<double>```, so a missing or incomplete total is distinct from a valid zero. The available names are ```n_gen_events```, ```n_proc_events```, and ```cosmic_livetime```. ```GetGeneratedEvents()```, ```GetProcessedEvents()```, and ```GetCosmicLivetime()``` are convenience accessors; ```GetRate(count)``` returns ```count / cosmic_livetime``` when a nonzero livetime is available. ```GetNProcEvents()``` remains available for compatibility and returns ```-1``` when the processed-event total is unavailable.
+
+Examples: [SubRunCounting.C](./examples/SubRunCounting.C), [SelectSubRuns.C](./examples/SelectSubRuns.C) for selecting SubRuns and retrieving the SubRun for an event, and [PlotCosmicTrackCrvZResidual_LivetimeNormalized.C](./examples/PlotCosmicTrackCrvZResidual_LivetimeNormalized.C) for normalizing a histogram with ```GetRate()```.
+
 ## The Event Class
 All branches and leaves can be accessed through the ```Event``` class like so:
 

rooutil/examples/PlotCosmicTrackCrvZResidual_LivetimeNormalized.C

@@ -0,0 +1,52 @@
+#include "EventNtuple/rooutil/inc/RooUtil.hh"
+#include "EventNtuple/rooutil/inc/common_cuts.hh"
+
+#include "TCanvas.h"
+#include "TH1F.h"
+
+#include <iostream>
+
+using namespace rooutil;
+
+bool crv_top_surface(TrackSegment& segment) {
+  if (segment.trkseg == nullptr) return false;
+
+  const int surface_id = segment.trkseg->sid;
+  return surface_id == mu2e::SurfaceIdDetail::TCRV; // note: this is a deprecated surface id, but it is what is in the file I'm testing
+}
+
+void PlotCosmicTrackCrvZResidual_LivetimeNormalized(std::string filename) {
+  RooUtil util(filename);
+  TH1F* crv_track_z_residual = new TH1F("crv_track_z_residual", "CRV coincidence-track z residual;z_{CRV coincidence} - z_{track} [mm];tracks / s", 100, -500, 500);
+
+  for (int i_event = 0; i_event < util.GetNEvents(); ++i_event) {
+    auto& event = util.GetEvent(i_event);
+    const auto tracks = event.GetTracks();
+    auto crv_coincs = event.GetCrvCoincs();
+    for (auto track : tracks) {
+      auto crv_top_segments = track.GetSegments([](TrackSegment& segment) {
+        return crv_top_surface(segment) && has_reco_step(segment);
+      });
+      for (auto& segment : crv_top_segments) {
+        for (auto& crv_coinc : crv_coincs) {
+          if (track_crv_coincidence(segment, crv_coinc)) {
+            crv_track_z_residual->Fill(crv_coinc.reco->pos.z() - segment.trkseg->pos.z());
+          }
+        }
+      }
+    }
+  }
+
+  const auto entries = crv_track_z_residual->Integral();
+  if (entries > 0.0) {
+    if (const auto rate = util.GetRate(entries)) {
+      crv_track_z_residual->Scale(*rate / entries);
+    } else {
+      std::cout << "Cannot normalize histogram: cosmic livetime is unavailable." << std::endl;
+    }
+  }
+
+  TCanvas* canvas = new TCanvas("canvas", "CRV coincidence-track z residual");
+  crv_track_z_residual->Draw("HIST");
+  canvas->SaveAs("crv_track_z_residual_rate.pdf");
+}

rooutil/examples/SelectSubRuns.C

@@ -0,0 +1,37 @@
+#include "EventNtuple/rooutil/inc/RooUtil.hh"
+
+#include <iostream>
+
+using namespace rooutil;
+
+bool has_more_than_100_processed_events(const SubRun& subrun) {
+  return subrun.has_proc_event_count && subrun.procEventCount > 100;
+}
+
+void SelectSubRuns(std::string filename) {
+  RooUtil util(filename);
+
+  const auto selected_subruns = util.GetSubRunIndices(has_more_than_100_processed_events);
+  std::cout << "SubRuns with more than 100 processed events:" << std::endl;
+  for (const auto i_subrun : selected_subruns) {
+    const auto& subrun = util.GetSubRun(i_subrun);
+    std::cout << subrun.srinfo->run << ":" << subrun.srinfo->subrun
+              << " has " << subrun.procEventCount << " processed events" << std::endl;
+  }
+
+  if (util.GetNEvents() == 0) return;
+
+  const auto& event = util.GetEvent(0);
+  const auto* event_subrun = util.FindSubRun(event.evtinfo->run, event.evtinfo->subrun);
+  if (event_subrun == nullptr) {
+    std::cout << "No SubRun information is available for the first event." << std::endl;
+    return;
+  }
+
+  std::cout << "First event belongs to " << event_subrun->srinfo->run << ":"
+            << event_subrun->srinfo->subrun;
+  if (event_subrun->has_cosmic_livetime) {
+    std::cout << " with " << event_subrun->cosmicLivetime << " s of cosmic livetime";
+  }
+  std::cout << std::endl;
+}

rooutil/examples/SubRunCounting.C

@@ -0,0 +1,29 @@
+#include "EventNtuple/rooutil/inc/RooUtil.hh"
+
+#include <iostream>
+
+using namespace rooutil;
+
+void SubRunCounting(std::string filename) {
+  RooUtil util(filename);
+
+  std::cout << filename << " contains " << util.GetNSubRuns() << " subruns" << std::endl;
+  for (int i_subrun = 0; i_subrun < util.GetNSubRuns(); ++i_subrun) {
+    const auto& subrun = util.GetSubRun(i_subrun);
+    std::cout << "run " << subrun.srinfo->run << ", subrun " << subrun.srinfo->subrun;
+    if (subrun.has_gen_event_count) { std::cout << ": " << subrun.genEventCount << " generated events"; }
+    if (subrun.has_proc_event_count) { std::cout << ", " << subrun.procEventCount << " processed events"; }
+    if (subrun.has_cosmic_livetime) { std::cout << ", " << subrun.cosmicLivetime << " s of cosmic livetime"; }
+    std::cout << std::endl;
+  }
+
+  if (const auto generated_events = util.GetGeneratedEvents()) {
+    std::cout << "Total generated events: " << *generated_events << std::endl;
+  }
+  if (const auto processed_events = util.GetProcessedEvents()) {
+    std::cout << "Total processed events: " << *processed_events << std::endl;
+  }
+  if (const auto cosmic_livetime = util.GetCosmicLivetime()) {
+    std::cout << "Total cosmic livetime: " << *cosmic_livetime << " s" << std::endl;
+  }
+}

rooutil/inc/BranchUtils.hh

@@ -0,0 +1,13 @@
+#ifndef BranchUtils_hh_
+#define BranchUtils_hh_
+
+#include "TChain.h"
+
+namespace rooutil {
+  inline bool CheckForBranch(TChain* ntuple, const char* branch_name, void* address = nullptr) {
+    if (ntuple->GetBranch(branch_name) == nullptr || ntuple->GetBranchStatus(branch_name) == 0) return false;
+    if (address != nullptr) ntuple->SetBranchAddress(branch_name, address);
+    return true;
+  }
+} // namespace rooutil
+#endif

rooutil/inc/Event.hh

@@ -47,8 +47,7 @@
 #include "EventNtuple/rooutil/inc/CaloCluster.hh"
 #include "EventNtuple/rooutil/inc/Trigger.hh"
 #include "EventNtuple/rooutil/inc/CaloHit.hh"
-
-#include "TChain.h"
+#include "EventNtuple/rooutil/inc/BranchUtils.hh"
 
 namespace rooutil {
   struct Event {
@@ -100,13 +99,6 @@ namespace rooutil {
       CheckForBranch(ntuple, "mcsteps_virtualdetector", &this->mcsteps_virtualdetector);
     }
 
-    // Check if a branch exists in the TChain, and optionally set its address
-    bool CheckForBranch(TChain* ntuple, const char* branch_name, void* address = nullptr) {
-      if(ntuple->GetBranch(branch_name) == nullptr || ntuple->GetBranchStatus(branch_name) == 0) return false;
-      if(address != nullptr) ntuple->SetBranchAddress(branch_name, address);
-      return true;
-    }
-
     // Add trigger branches and store the path name information
     void AddTriggerInfo(TChain* ntuple) {
       trigger.SetTrigInfo(&triginfo); // pointer to the underlying trigger data read in event-by-event