diff --git a/README.md b/README.md index caef2b2c..3267241a 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,8 @@ 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: @@ -31,7 +33,7 @@ You can get some basic information about the EventNtuple file useing ```checkEve 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 diff --git a/bin/checkEventNtuple b/bin/checkEventNtuple index 5899ec45..c73e6816 100755 --- a/bin/checkEventNtuple +++ b/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: diff --git a/fcl/from_dig-OnSpill.fcl b/fcl/from_dig-OnSpill.fcl index f8b54a19..f2b513a4 100644 --- a/fcl/from_dig-OnSpill.fcl +++ b/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 diff --git a/fcl/from_dig-calo.fcl b/fcl/from_dig-calo.fcl index d3318d7d..b1dcd364 100644 --- a/fcl/from_dig-calo.fcl +++ b/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 diff --git a/fcl/from_dig-mockdata.fcl b/fcl/from_dig-mockdata.fcl index af375ddd..a34489b1 100644 --- a/fcl/from_dig-mockdata.fcl +++ b/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 diff --git a/fcl/from_mcs-DeCalib.fcl b/fcl/from_mcs-DeCalib.fcl index b2c528e8..4c1476df 100644 --- a/fcl/from_mcs-DeCalib.fcl +++ b/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 diff --git a/fcl/from_mcs-ceSimReco.fcl b/fcl/from_mcs-ceSimReco.fcl index c32c7550..67ae1310 100644 --- a/fcl/from_mcs-ceSimReco.fcl +++ b/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 ] diff --git a/fcl/from_mcs-extracted.fcl b/fcl/from_mcs-extracted.fcl index 330f0efd..b1f020bc 100644 --- a/fcl/from_mcs-extracted.fcl +++ b/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" diff --git a/fcl/from_mcs-mockdata.fcl b/fcl/from_mcs-mockdata.fcl index 04f965fb..60b32434 100644 --- a/fcl/from_mcs-mockdata.fcl +++ b/fcl/from_mcs-mockdata.fcl @@ -16,7 +16,7 @@ physics : } physics.EventNtuplePath : [ @sequence::EventNtuple.Path ] -physics.EventNtupleEndPath : [ @sequence::EventNtuple.EndPathNoMC ] +physics.EventNtupleEndPath : [ @sequence::EventNtuple.EndPath ] physics.trigger_paths : [ EventNtuplePath ] physics.end_paths : [ EventNtupleEndPath ] @@ -24,3 +24,5 @@ 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 diff --git a/fcl/from_mcs-mockdata_noMC.fcl b/fcl/from_mcs-mockdata_noMC.fcl index 72698c77..e4b08dd0 100644 --- a/fcl/from_mcs-mockdata_noMC.fcl +++ b/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 ] diff --git a/fcl/from_mcs-primary.fcl b/fcl/from_mcs-primary.fcl index e9650a82..21e5605e 100644 --- a/fcl/from_mcs-primary.fcl +++ b/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 diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index 4b8b68d8..f809b366 100644 --- a/fcl/prolog.fcl +++ b/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 ] } diff --git a/inc/SubRunInfo.hh b/inc/SubRunInfo.hh new file mode 100644 index 00000000..8e4a9b70 --- /dev/null +++ b/inc/SubRunInfo.hh @@ -0,0 +1,16 @@ +// +// SubRunInfo: subrun information +// Andy Edmonds (2026) +// +#ifndef SubRunInfo_HH +#define SubRunInfo_HH +#include +namespace mu2e +{ + struct SubRunInfo { + int run = 0; // run number + int subrun = 0; // subrun number + void reset() {*this = SubRunInfo(); } + }; +} +#endif diff --git a/rooutil/README.md b/rooutil/README.md index 54af35a7..056faa59 100644 --- a/rooutil/README.md +++ b/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```, 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: diff --git a/rooutil/examples/PlotCosmicTrackCrvZResidual_LivetimeNormalized.C b/rooutil/examples/PlotCosmicTrackCrvZResidual_LivetimeNormalized.C new file mode 100644 index 00000000..3238bfa7 --- /dev/null +++ b/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 + +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"); +} diff --git a/rooutil/examples/SelectSubRuns.C b/rooutil/examples/SelectSubRuns.C new file mode 100644 index 00000000..7195be64 --- /dev/null +++ b/rooutil/examples/SelectSubRuns.C @@ -0,0 +1,37 @@ +#include "EventNtuple/rooutil/inc/RooUtil.hh" + +#include + +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; +} diff --git a/rooutil/examples/SubRunCounting.C b/rooutil/examples/SubRunCounting.C new file mode 100644 index 00000000..41b34395 --- /dev/null +++ b/rooutil/examples/SubRunCounting.C @@ -0,0 +1,29 @@ +#include "EventNtuple/rooutil/inc/RooUtil.hh" + +#include + +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; + } +} diff --git a/rooutil/inc/BranchUtils.hh b/rooutil/inc/BranchUtils.hh new file mode 100644 index 00000000..3a382231 --- /dev/null +++ b/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 diff --git a/rooutil/inc/Event.hh b/rooutil/inc/Event.hh index e8789f3e..a7205bd6 100644 --- a/rooutil/inc/Event.hh +++ b/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 diff --git a/rooutil/inc/RooUtil.hh b/rooutil/inc/RooUtil.hh index 0ba2cdd4..993ab791 100644 --- a/rooutil/inc/RooUtil.hh +++ b/rooutil/inc/RooUtil.hh @@ -2,25 +2,30 @@ #define RooUtil_hh_ #include +#include +#include +#include +#include #include "TFile.h" #include "TTree.h" #include "TH1I.h" #include "EventNtuple/rooutil/inc/Event.hh" +#include "EventNtuple/rooutil/inc/SubRun.hh" namespace rooutil { class RooUtil { public: - RooUtil(std::string filename, bool debug = false, std::string treename = "EventNtuple/ntuple") : debug(debug), n_proc_events(0) { + RooUtil(std::string filename, bool debug = false, std::string treename = "EventNtuple/ntuple") : debug(debug) { ntuple = new TChain(treename.c_str()); + subrun_ntuple = new TChain("EventNtuple/subrunNtuple"); // Check if the given filename contains .root at the end std::string root_suffix = ".root"; - if (filename.compare(filename.size() - root_suffix.size(), root_suffix.size(), root_suffix) == 0) { - ntuple->Add(filename.c_str()); + if (filename.size() >= root_suffix.size() && filename.compare(filename.size() - root_suffix.size(), root_suffix.size(), root_suffix) == 0) { + AddFile(filename); SetVersionNumber(filename); - SetNProcessedEvents(filename); } else { // assume its a file list std::ifstream filelist(filename); @@ -29,13 +34,13 @@ namespace rooutil { std::string line; bool first_line = true; while (std::getline(filelist, line)) { - ntuple->Add(line.c_str()); + if (line.empty()) continue; + AddFile(line); if (first_line) { SetVersionNumber(line); first_line = false; } - SetNProcessedEvents(line); } filelist.close(); } else { @@ -44,6 +49,7 @@ namespace rooutil { } event = new Event(ntuple); + subrun = new SubRun(subrun_ntuple); } void Debug(bool dbg) { debug = dbg; } @@ -68,21 +74,63 @@ namespace rooutil { delete file; } - void SetNProcessedEvents(std::string filename) { - TFile* file = new TFile(filename.c_str(), "READ"); - TH1I* hProcEvents = (TH1I*) file->Get("EventNtuple/n_proc_events"); - if (!hProcEvents) { - std::cout << "Warning: this EventNtuple file does not contain the n_proc_events histogram. It is either v06_09_02 or older. This is just a warning..." << std::endl; + void AddFile(const std::string& filename) { + ntuple->Add(filename.c_str()); + TFile file(filename.c_str(), "READ"); + if (!file.IsZombie() && file.Get("EventNtuple/subrunNtuple") != nullptr) { + subrun_ntuple->Add(filename.c_str()); } - else { - n_proc_events += hProcEvents->GetBinContent(1); + SetTotals(filename); + } + + void SetTotals(const std::string& filename) { + TFile* file = new TFile(filename.c_str(), "READ"); + for (const auto& total_name : TotalNames()) { + TH1* histogram = dynamic_cast(file->Get(("EventNtuple/" + total_name).c_str())); + if (histogram == nullptr) { + incomplete_totals.insert(total_name); + continue; + } + totals[total_name] += histogram->GetBinContent(1); } file->Close(); delete file; } int GetNEvents() { return ntuple->GetEntries(); } - int GetNProcEvents() { return n_proc_events; } + int GetNProcEvents() { + const auto total = GetTotal("n_proc_events"); + return total ? static_cast(*total) : -1; + } + + std::optional GetTotal(const std::string& total_name) const { + if (incomplete_totals.count(total_name) != 0) return std::nullopt; + const auto total = totals.find(total_name); + if (total == totals.end()) return std::nullopt; + return total->second; + } + + bool HasTotal(const std::string& total_name) const { + return GetTotal(total_name).has_value(); + } + + std::vector GetAvailableTotals() const { + std::vector available_totals; + for (const auto& total_name : TotalNames()) { + if (HasTotal(total_name)) available_totals.push_back(total_name); + } + return available_totals; + } + + std::optional GetGeneratedEvents() const { return GetTotal("n_gen_events"); } + std::optional GetProcessedEvents() const { return GetTotal("n_proc_events"); } + std::optional GetCosmicLivetime() const { return GetTotal("cosmic_livetime"); } + + std::optional GetRate(double count) const { + const auto livetime = GetCosmicLivetime(); + if (!livetime || *livetime == 0.0) return std::nullopt; + return count / *livetime; + } Event& GetEvent(int i_event) { if (debug) { std::cout << "RooUtil::GetEvent(): Getting event " << i_event << std::endl; } @@ -95,6 +143,33 @@ namespace rooutil { return *event; } + int GetNSubRuns() { return subrun_ntuple->GetEntries(); } + + SubRun& GetSubRun(int i_subrun) { + if (i_subrun < 0 || i_subrun >= GetNSubRuns()) { + throw std::out_of_range("RooUtil::GetSubRun(): SubRun index is out of range"); + } + subrun_ntuple->GetEntry(i_subrun); + return *subrun; + } + + SubRun* FindSubRun(int run, int subrun_number) { + BuildSubRunIndex(); + const auto entry = subrun_indices.find(std::make_pair(run, subrun_number)); + if (entry == subrun_indices.end()) return nullptr; + return &GetSubRun(entry->second); + } + + using SubRunCut = bool (*)(const SubRun&); + + std::vector GetSubRunIndices(SubRunCut cut) { + std::vector selected_subruns; + for (int i_subrun = 0; i_subrun < GetNSubRuns(); ++i_subrun) { + if (cut(GetSubRun(i_subrun))) selected_subruns.push_back(i_subrun); + } + return selected_subruns; + } + void TurnOffBranch(const std::string& branchname) { ntuple->SetBranchStatus(branchname.c_str(), 0); } @@ -117,9 +192,22 @@ namespace rooutil { void TurnOnAllBranches() { TurnOnBranch("*"); } + void TurnOffSubRunBranch(const std::string& branchname) { + subrun_ntuple->SetBranchStatus(branchname.c_str(), 0); + } + void TurnOnSubRunBranch(const std::string& branchname) { + subrun_ntuple->SetBranchStatus(branchname.c_str(), 1); + } + void TurnOffAllSubRunBranches() { + TurnOffSubRunBranch("*"); + } + void TurnOnAllSubRunBranches() { + TurnOnSubRunBranch("*"); + } void CreateOutputEventNtuple(TFile* outfile) { - auto dir = outfile->mkdir("EventNtuple"); + auto dir = outfile->GetDirectory("EventNtuple"); + if (dir == nullptr) dir = outfile->mkdir("EventNtuple"); dir->cd(); output_ntuple = new TTree("ntuple", "reduced ntuple"); @@ -172,22 +260,64 @@ namespace rooutil { if (event->mcsteps_virtualdetector) { output_ntuple->Branch("mcsteps_virtualdetector", event->mcsteps_virtualdetector); } // Write out histograms from input to output - hVersion->Write(); + if (hVersion != nullptr) hVersion->Write(); } void FillOutputEventNtuple() { output_ntuple->Fill(); } + void CreateOutputSubRunNtuple(TFile* outfile) { + auto dir = outfile->GetDirectory("EventNtuple"); + if (dir == nullptr) dir = outfile->mkdir("EventNtuple"); + dir->cd(); + output_subrun_ntuple = new TTree("subrunNtuple", "reduced subrun ntuple"); + + if (subrun->has_srinfo) { output_subrun_ntuple->Branch("srinfo", subrun->srinfo); } + if (subrun->has_gen_event_count) { output_subrun_ntuple->Branch("genEventCount", &subrun->genEventCount); } + if (subrun->has_proc_event_count) { output_subrun_ntuple->Branch("procEventCount", &subrun->procEventCount); } + if (subrun->has_cosmic_livetime) { output_subrun_ntuple->Branch("cosmicLivetime", &subrun->cosmicLivetime); } + } + + void FillOutputSubRunNtuple() { + output_subrun_ntuple->Fill(); + } + private: + static const std::vector& TotalNames() { + static const std::vector total_names{ + "n_gen_events", "n_proc_events", "cosmic_livetime" + }; + return total_names; + } + + void BuildSubRunIndex() { + if (subrun_index_built) return; + if (!subrun->has_srinfo) { + subrun_index_built = true; + return; + } + for (int i_subrun = 0; i_subrun < GetNSubRuns(); ++i_subrun) { + const auto& current_subrun = GetSubRun(i_subrun); + subrun_indices[std::make_pair(current_subrun.srinfo->run, current_subrun.srinfo->subrun)] = i_subrun; + } + subrun_index_built = true; + } + TChain* ntuple; Event* event; // holds all the variables for SetBranchAddress + TChain* subrun_ntuple; + SubRun* subrun; bool debug; - TH1I* hVersion; - int n_proc_events; + TH1I* hVersion = nullptr; + std::map totals; + std::set incomplete_totals; + std::map, int> subrun_indices; + bool subrun_index_built = false; - TTree* output_ntuple; // for output + TTree* output_ntuple = nullptr; // for output + TTree* output_subrun_ntuple = nullptr; }; } // namespace rooutil #endif diff --git a/rooutil/inc/SubRun.hh b/rooutil/inc/SubRun.hh new file mode 100644 index 00000000..dfe087b7 --- /dev/null +++ b/rooutil/inc/SubRun.hh @@ -0,0 +1,27 @@ +#ifndef SubRun_hh_ +#define SubRun_hh_ + +#include "EventNtuple/inc/SubRunInfo.hh" +#include "EventNtuple/rooutil/inc/BranchUtils.hh" + +namespace rooutil { + struct SubRun { + explicit SubRun(TChain* ntuple) { + has_srinfo = CheckForBranch(ntuple, "srinfo", &srinfo); + has_gen_event_count = CheckForBranch(ntuple, "genEventCount", &genEventCount); + has_proc_event_count = CheckForBranch(ntuple, "procEventCount", &procEventCount); + has_cosmic_livetime = CheckForBranch(ntuple, "cosmicLivetime", &cosmicLivetime); + } + + mu2e::SubRunInfo* srinfo = nullptr; + long genEventCount = 0; + long procEventCount = 0; + float cosmicLivetime = 0.0; + + bool has_srinfo = false; + bool has_gen_event_count = false; + bool has_proc_event_count = false; + bool has_cosmic_livetime = false; + }; +} // namespace rooutil +#endif diff --git a/src/EventNtupleMaker_module.cc b/src/EventNtupleMaker_module.cc index 94991962..5f7bd378 100644 --- a/src/EventNtupleMaker_module.cc +++ b/src/EventNtupleMaker_module.cc @@ -10,6 +10,8 @@ #include "Offline/MCDataProducts/inc/CaloClusterMC.hh" #include "Offline/MCDataProducts/inc/CaloHitMC.hh" #include "Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh" +#include "Offline/MCDataProducts/inc/GenEventCount.hh" +#include "Offline/MCDataProducts/inc/CosmicLivetime.hh" #include "Offline/RecoDataProducts/inc/KalSeed.hh" #include "Offline/RecoDataProducts/inc/KalSeedAssns.hh" #include "Offline/RecoDataProducts/inc/CaloCluster.hh" @@ -85,6 +87,7 @@ #include "EventNtuple/inc/MCStepSummaryInfo.hh" #include "EventNtuple/inc/CaloDigiMCInfo.hh" #include "Offline/MCDataProducts/inc/CaloShowerSim.hh" +#include "EventNtuple/inc/SubRunInfo.hh" // C++ includes. #include @@ -263,6 +266,30 @@ namespace mu2e { fhicl::Table mc{Name("mc"), Comment("CRV MC filling options")}; }; + // Configuration for specific subrun quantities + struct SubRunBranchConfig { + using Name=fhicl::Name; + using Comment=fhicl::Comment; + + fhicl::Atom include{Name("include"), Comment("True/false to include this quantity or not (master switch)")}; + fhicl::Atom inputTag{Name("inputTag"), Comment("InputTag for subrun quantity")}; + fhicl::Atom fillNtuple{Name("fillNtuple"), Comment("True/false to add this quantity to the ntuple")}; + fhicl::Atom fillTotalsHistogram{Name("fillTotalsHistogram"), Comment("True/false to make a totals histogram for this quantity")}; + }; + + // Overal configuration for subrun ntuple and histograms + struct SubRunConfig { + using Name=fhicl::Name; + using Comment=fhicl::Comment; + + fhicl::Atom makeNtuple{Name("makeNtuple"), Comment("True/false of whether to make the subrun ntuple")}; + fhicl::Atom makeTotalsHistograms{Name("makeTotalsHistograms"), Comment("True/false of whether to make the totals histograms")}; + fhicl::Table genEventCount{Name("genEventCount"), Comment("Number of generated events (MC)")}; + fhicl::Table procEventCount{Name("procEventCount"), Comment("Number of processed events")}; + fhicl::Table cosmicLivetime{Name("cosmicLivetime"), Comment("Cosmic livetime [s]")}; + }; + + struct Config { using Name=fhicl::Name; using Comment=fhicl::Comment; @@ -293,6 +320,7 @@ namespace mu2e { fhicl::Table helices {Name("helices"), Comment("Helix seed branches config")}; fhicl::Table timeclusters{Name("timeclusters"),Comment("Time cluster branch config")}; fhicl::Table mcsteps {Name("mcsteps"), Comment("MC step collection branches config")}; + fhicl::Table subruns {Name("subruns"), Comment("SubRun ntuple/histogram config")}; }; typedef art::EDAnalyzer::Table Parameters; @@ -302,6 +330,7 @@ namespace mu2e { void beginJob() override; void beginSubRun(const art::SubRun & subrun ) override; void analyze(const art::Event& e) override; + void endSubRun(const art::SubRun & subrun ) override; private: @@ -310,7 +339,11 @@ namespace mu2e { // main TTree TTree* _ntuple; TH1I* _hVersion; - TH1I* _hProcEvents; + // subrun quantities + TTree* _subrunNtuple; + TH1I* _hGenEventCount; + TH1I* _hProcEventCount; + TH1F* _hCosmicLivetime; // general event info branch EventInfo _einfo; EventInfoMC _einfomc; @@ -440,6 +473,12 @@ namespace mu2e { // for trigger branch: bool firstEvent = true; + // For subrun ntuple + mu2e::SubRunInfo _srinfo; + long _genEventCount; + long _procEventCount; + float _cosmicLivetime; + // ── Gating helpers ───────────────────────────────────────────────────── // Event-level MC gate (simParticles, mcTrajectories, primaryParticle). // Also gates tracker MC and CRV MC. @@ -567,7 +606,6 @@ namespace mu2e { _hVersion->GetXaxis()->SetBinLabel(1, "major"); _hVersion->SetBinContent(1, 6); _hVersion->GetXaxis()->SetBinLabel(2, "minor"); _hVersion->SetBinContent(2, 12); _hVersion->GetXaxis()->SetBinLabel(3, "patch"); _hVersion->SetBinContent(3, 1); - _hProcEvents = tfs->make("n_proc_events", "number of processed events", 1,0,1); // event info branch _ntuple->Branch("evtinfo",&_einfo,_buffsize,_splitlevel); if (fillEventMC()) { @@ -710,10 +748,80 @@ namespace mu2e { _ntuple->Branch(("mcsteps_"+inst+".").c_str(),&_stepPointMCInfos[icoll],_buffsize,_splitlevel); } } + + // Subrun tree + if (_conf.subruns().makeNtuple()) { + _subrunNtuple = tfs->make("subrunNtuple","Mu2e SubRun Ntuple"); + _subrunNtuple->Branch("srinfo", &_srinfo, _buffsize, _splitlevel); + + if (_conf.subruns().genEventCount().include() && _conf.subruns().genEventCount().fillNtuple()) { + _subrunNtuple->Branch("genEventCount", &_genEventCount, _buffsize, _splitlevel); + } + + if (_conf.subruns().procEventCount().include() && _conf.subruns().procEventCount().fillNtuple()) { + _subrunNtuple->Branch("procEventCount", &_procEventCount, _buffsize, _splitlevel); + } + + if (_conf.subruns().cosmicLivetime().include() && _conf.subruns().cosmicLivetime().fillNtuple()) { + _subrunNtuple->Branch("cosmicLivetime", &_cosmicLivetime, _buffsize, _splitlevel); + } + + } + if (_conf.subruns().makeTotalsHistograms()) { + if (_conf.subruns().genEventCount().include() && _conf.subruns().genEventCount().fillTotalsHistogram()) { + _hGenEventCount = tfs->make("n_gen_events", "total number of generated events", 1,0,1); + } + if (_conf.subruns().procEventCount().include() && _conf.subruns().procEventCount().fillTotalsHistogram()) { + _hProcEventCount = tfs->make("n_proc_events", "total number of processed events", 1,0,1); + } + if (_conf.subruns().cosmicLivetime().include() && _conf.subruns().cosmicLivetime().fillTotalsHistogram()) { + _hCosmicLivetime = tfs->make("cosmic_livetime", "total cosmic livetime [sec]", 1,0,1); + } + + } } void EventNtupleMaker::beginSubRun(const art::SubRun & subrun ) { _infoStructHelper.updateSubRun(); + + _srinfo.run = subrun.run(); + _srinfo.subrun = subrun.subRun(); + + } + + void EventNtupleMaker::endSubRun(const art::SubRun & subrun ) { + + if (_conf.subruns().genEventCount().include()) { + art::Handle genEventCountHandle; + subrun.getByLabel(_conf.subruns().genEventCount().inputTag(), genEventCountHandle); + + _genEventCount = genEventCountHandle->count(); + if (_conf.subruns().makeTotalsHistograms() && _conf.subruns().genEventCount().fillTotalsHistogram()) { + _hGenEventCount->Fill(0.0, genEventCountHandle->count()); + } + } + + if (_conf.subruns().cosmicLivetime().include()) { + art::Handle cosmicLivetimeHandle; + subrun.getByLabel(_conf.subruns().cosmicLivetime().inputTag(), cosmicLivetimeHandle); + + _cosmicLivetime = cosmicLivetimeHandle->liveTime(); + if (_conf.subruns().makeTotalsHistograms() && _conf.subruns().cosmicLivetime().fillTotalsHistogram()) { + _hCosmicLivetime->Fill(0.0, cosmicLivetimeHandle->liveTime()); + } + } + + if (_conf.subruns().makeNtuple()) { + _subrunNtuple->Fill(); + } + + if (_conf.subruns().procEventCount().include()) { + if (_conf.subruns().makeTotalsHistograms() && _conf.subruns().procEventCount().fillTotalsHistogram()) { + _hProcEventCount->Fill(0.0, _procEventCount); + } + } + + _procEventCount = 0; // reset procEventCount every subrun } void EventNtupleMaker::resolveCrvPlanes() { @@ -1167,7 +1275,7 @@ namespace mu2e { if(fill) { _ntuple->Fill(); } - _hProcEvents->Fill(0); + _procEventCount++; } diff --git a/src/classes.h b/src/classes.h index 95f93a4f..356b6341 100644 --- a/src/classes.h +++ b/src/classes.h @@ -37,3 +37,4 @@ #include "EventNtuple/inc/MVAResultInfo.hh" #include "EventNtuple/inc/SurfaceStepInfo.hh" #include "EventNtuple/inc/TrigInfo.hh" +#include "EventNtuple/inc/SubRunInfo.hh" diff --git a/src/classes_def.xml b/src/classes_def.xml index c086ee96..b4177d73 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -78,4 +78,5 @@ + diff --git a/validation/test_rooutil_examples.sh b/validation/test_rooutil_examples.sh index 4fabbb25..6a1920bf 100644 --- a/validation/test_rooutil_examples.sh +++ b/validation/test_rooutil_examples.sh @@ -5,7 +5,7 @@ scripts=( "PrintEvents.C" "CreateNtuple.C" "CreateTrackNtuple.C" "PlotCRVPEs.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" "PlotVDSteps.C" "PlotCaloClusterEnergy_RecoVsTrue.C" - "PlotCaloClusterAndHits.C" "PlotCaloCluster_SimParticles.C" + "PlotCaloClusterAndHits.C" "PlotCaloCluster_SimParticles.C" "SubRunCounting.C" "SelectSubRuns.C" "PlotCosmicTrackCrvZResidual_LivetimeNormalized.C" ) for script in "${scripts[@]}"