feat: set histogram boundaries according to binning scheme

macro/analysis_modify_histos.C

@@ -0,0 +1,102 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2022 Christopher Dilks
+
+R__LOAD_LIBRARY(EpicAnalysis)
+
+// analysis in bins of (x,Q2); tests modifying bin boundaries of Histos histograms
+// according to the binning scheme
+void analysis_modify_histos(
+    TString configFile,
+    TString outfilePrefix,
+    TString dataSource="epic",
+    TString reconMethod="Ele"
+    )
+{
+
+  // setup analysis ========================================
+  Analysis *A;
+  if   (dataSource.Contains("epic"))     A = new AnalysisEpic(   configFile, outfilePrefix );
+  else if(dataSource.Contains("athena")) A = new AnalysisAthena( configFile, outfilePrefix );
+  else if(dataSource.Contains("ecce"))   A = new AnalysisEcce(   configFile, outfilePrefix );
+#ifndef EXCLUDE_DELPHES
+  else A = new AnalysisDelphes( configFile, outfilePrefix );
+#endif
+
+  A->SetReconMethod(reconMethod); // set reconstruction method
+  A->AddFinalState("pipTrack"); // pion final state
+  // A->writeSimpleTree = true;
+
+  // define cuts ===========================================
+  A->AddBinScheme("w");  A->BinScheme("w")->BuildBin("Min",3.0); // W > 3 GeV
+  A->AddBinScheme("y");  A->BinScheme("y")->BuildBin("Range",0.01,0.95); // 0.01 < y < 0.95
+  A->AddBinScheme("z");  A->BinScheme("z")->BuildBin("Range",0.2,0.9); // 0.2 < z < 0.9
+  A->AddBinScheme("xF"); A->BinScheme("xF")->BuildBin("Min",0.0); // xF > 0
+  A->AddBinScheme("ptLab");  A->BinScheme("ptLab")->BuildBin("Min",0.1); // pT_lab > 0.1 GeV (tracking limit)
+
+  // set binning scheme ====================================
+  A->AddBinScheme("x");  A->BinScheme("x")->BuildBins(  6, 0.001, 1,    true );
+  A->AddBinScheme("q2"); A->BinScheme("q2")->BuildBins( 4, 1,     3000, true );
+
+
+  //////////////////////////////////////////////////////////
+  // modify histograms according to the binning scheme
+
+  // get the HistosDAG (adage)
+  A->Prepare(); // (to initialize the DAG with the above bins)
+  auto HD = A->GetHistosDAG();
+
+  // define a function to modify histogram named `hist_name` to
+  // have the bin boundaries of the bin named `var_name`; this will
+  // run on every possible multi-dimensional bin
+  auto mod_bounds = [&HD] (TString hist_name, TString var_name) {
+
+    // payload function, to be run on every multi-dim. bin
+    auto payload = [&hist_name,&var_name] (NodePath *NP, Histos *H ) {
+
+      // get the histogram and bin
+      auto hist_orig   = H->Hist(hist_name);
+      auto hist_config = H->GetHistConfig(hist_name);
+      auto var_bin     = NP->GetBinNode(var_name);
+      if(hist_orig==nullptr or hist_config==nullptr or var_bin==nullptr)
+        return;
+
+      // get the binning info // FIXME: only works for 1D
+      auto n_bins  = hist_orig->GetXaxis()->GetNbins();
+      auto var_min = var_bin->GetCut()->GetMin();
+      auto var_max = var_bin->GetCut()->GetMax();
+
+      // get histogram info
+      auto hist_orig_name  = hist_orig->GetName();
+      auto hist_orig_title = hist_orig->GetTitle();
+
+      // create new histogram, replacing the old one
+      decltype(hist_orig) hist_new;
+      switch(hist_orig->GetDimension()) {
+        case 1:
+          // FIXME: `Warning in <TFile::Append>: Replacing existing TH1`
+          hist_new = new TH1D(hist_orig_name, hist_orig_title, n_bins, var_min, var_max);
+          if(hist_config->logx) BinSet::BinLog(hist_new->GetXaxis()); // preserve equal bin-width in log scale setting
+          H->ReplaceHist(hist_name, hist_new);
+          break;
+        default:
+          fmt::print(stderr,"ERROR: {}-dimensional histogram not fully supported yet\n");
+      }
+    };
+
+    // run the payload on every multi-dim. bin
+    HD->Payload(payload);
+    HD->ExecuteAndClearOps();
+  };
+
+  // make 1D x distribution have a range that matches the x-bin boundaries
+  mod_bounds("x", "x");
+
+  // repeat for Q2; note the variable and histogram names differ
+  mod_bounds("Q2", "q2");
+
+  //////////////////////////////////////////////////////////
+
+
+  // perform the analysis ==================================
+  A->Execute();
+};

src/Analysis.cxx

@@ -22,6 +22,7 @@ Analysis::Analysis(
   , ionBeamEn(100.0)
   , crossingAngle(-25.0)
   , totalCrossSection(1e6)
+  , prepared(false)
 {
   // available variables for binning
   // - availableBinSchemes is a map from variable name to variable title
@@ -170,6 +171,9 @@ void Analysis::AddFileGroup(
 //------------------------------------
 void Analysis::Prepare() {
 
+  // run only once
+  if(prepared) return;
+
   // parse config file --------------------
   std::ifstream fin(configFileName);
   std::string line;
@@ -564,6 +568,9 @@ void Analysis::Prepare() {
   wInclusiveTotal = 0.;
   wTrackTotal     = 0.;
   wJetTotal       = 0.;
+
+  // run only once
+  prepared = true;
 };
 
 void Analysis::CalculateEventQ2Weights() {

src/Analysis.h

@@ -106,12 +106,13 @@ class Analysis : public TNamed
     // run the analysis
     virtual void Execute() = 0;
 
-  protected:
-
     // prepare to perform the analysis; in derived classes, define a method `Execute()`, which
     // will run the event loop; the first line of `Execute()` should call `Analysis::Prepare()`,
     // which set up common things like output files, `HistosDAG`, etc.
     void Prepare();
+    bool prepared; // true if Prepared() has been called
+
+  protected:
 
     // finish the analysis; call `Analysis::Finish()` at the end of derived `Execute()` methods
     void Finish();

src/Histos.cxx

@@ -193,6 +193,13 @@ void Histos::RegisterHist4(TString varname_, Hist4D *hist_, HistConfig *config_)
   hist4Map.insert({varname_,hist_});
 };
 
+// replace a histogram (keeps its HistConfig)
+void Histos::ReplaceHist(TString varname_, TH1 *hist_) {
+  auto hist_orig = this->Hist(varname_);
+  if(hist_orig==nullptr) return;
+  delete hist_orig;
+  histMap[varname_] = hist_;
+};
 
 // access histogram by name
 TH1 *Histos::Hist(TString histName, Bool_t silence) {

src/Histos.h

@@ -123,6 +123,9 @@ class Histos : public TNamed
         Bool_t logz = false
         );
 
+    // replace a histogram (keeps its HistConfig); useful for re-binning
+    void ReplaceHist(TString varname_, TH1 *hist_);
+
     // writers
     void WriteHists(TFile *ofile) {
       ofile->cd("/");