diff --git a/macro/analysis_modify_histos.C b/macro/analysis_modify_histos.C new file mode 100644 index 00000000..fac4a629 --- /dev/null +++ b/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 : 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(); +}; diff --git a/src/Analysis.cxx b/src/Analysis.cxx index 7ddb545c..0407a856 100644 --- a/src/Analysis.cxx +++ b/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() { diff --git a/src/Analysis.h b/src/Analysis.h index e98b1832..cb650c52 100644 --- a/src/Analysis.h +++ b/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(); diff --git a/src/Histos.cxx b/src/Histos.cxx index 6e01ed9e..fb3da1e9 100644 --- a/src/Histos.cxx +++ b/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) { diff --git a/src/Histos.h b/src/Histos.h index bd03fbae..923933a0 100644 --- a/src/Histos.h +++ b/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("/");