Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-09-27 07:03:26 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Duane Byer
0003 
0004 #include "Hist4D.h"
0005 
0006 ClassImp(Hist4D);
0007 
0008 struct DifferentNumberOfBins : public std::runtime_error {
0009     DifferentNumberOfBins() :
0010         std::runtime_error("histograms have different numbers of bins") { }
0011 };
0012 struct DifferentAxisLimits : public std::runtime_error {
0013     DifferentAxisLimits() :
0014         std::runtime_error("histograms have different axis limits") { }
0015 };
0016 struct DifferentBinLimits : public std::runtime_error {
0017     DifferentBinLimits() :
0018         std::runtime_error("histograms have different bin limits") { }
0019 };
0020 struct DifferentLabels : public std::runtime_error {
0021     DifferentLabels() :
0022         std::runtime_error("histograms have different bin limits") { }
0023 };
0024 
0025 Hist4D::Hist4D(
0026         char const* name, char const* title,
0027         Int_t nbins_w, Double_t const* bins_w,
0028         Int_t nbins_x, Double_t const* bins_x,
0029         Int_t nbins_y, Double_t const* bins_y,
0030         Int_t nbins_z, Double_t const* bins_z) :
0031         TNamed(name, title),
0032         _w_axis(new TAxis(nbins_w, bins_w)),
0033         _x_axis(new TAxis(nbins_x, bins_x)),
0034         _y_axis(new TAxis(nbins_y, bins_y)),
0035         _z_axis(new TAxis(nbins_z, bins_z)),
0036         _hists() {
0037     CreateAxes();
0038     CreateHists();
0039 }
0040 
0041 Hist4D::Hist4D(
0042         char const* name, char const* title,
0043         Int_t nbins_w, Double_t w_lower, Double_t w_upper,
0044         Int_t nbins_x, Double_t x_lower, Double_t x_upper,
0045         Int_t nbins_y, Double_t y_lower, Double_t y_upper,
0046         Int_t nbins_z, Double_t z_lower, Double_t z_upper) :
0047         TNamed(name, title),
0048         _w_axis(new TAxis(nbins_w, w_lower, w_upper)),
0049         _x_axis(new TAxis(nbins_x, x_lower, x_upper)),
0050         _y_axis(new TAxis(nbins_y, y_lower, y_upper)),
0051         _z_axis(new TAxis(nbins_z, z_lower, z_upper)),
0052         _hists() {
0053     CreateAxes();
0054     CreateHists();
0055 }
0056 
0057 bool Hist4D::CheckConsistency(Hist4D const* h1, Hist4D const* h2) {
0058     if (h1 == h2) {
0059         return true;
0060     }
0061     Int_t nbins_w = h1->_w_axis->GetNbins();
0062     Int_t nbins_x = h1->_x_axis->GetNbins();
0063     Int_t nbins_y = h1->_y_axis->GetNbins();
0064     Int_t nbins_z = h1->_z_axis->GetNbins();
0065     if (nbins_w != h2->_w_axis->GetNbins()
0066             || nbins_x != h2->_x_axis->GetNbins()
0067             || nbins_y != h2->_y_axis->GetNbins()
0068             || nbins_z != h2->_z_axis->GetNbins()) {
0069         throw DifferentNumberOfBins();
0070         return false;
0071     }
0072     if (h1->_hists.size() != h2->_hists.size()) {
0073         throw DifferentNumberOfBins();
0074         return false;
0075     }
0076     bool result = true;
0077     result &= CheckAxisLimits(h1->_w_axis, h2->_w_axis);
0078     result &= CheckAxisLimits(h1->_x_axis, h2->_x_axis);
0079     result &= CheckAxisLimits(h1->_y_axis, h2->_y_axis);
0080     result &= CheckAxisLimits(h1->_z_axis, h2->_z_axis);
0081     result &= CheckBinLimits(h1->_w_axis, h2->_w_axis);
0082     result &= CheckBinLimits(h1->_x_axis, h2->_x_axis);
0083     result &= CheckBinLimits(h1->_y_axis, h2->_y_axis);
0084     result &= CheckBinLimits(h1->_z_axis, h2->_z_axis);
0085     result &= CheckBinLabels(h1->_w_axis, h2->_w_axis);
0086     result &= CheckBinLabels(h1->_x_axis, h2->_x_axis);
0087     result &= CheckBinLabels(h1->_y_axis, h2->_y_axis);
0088     result &= CheckBinLabels(h1->_z_axis, h2->_z_axis);
0089     return result;
0090 }
0091 
0092 bool Hist4D::CheckAxisLimits(TAxis const* ax1, TAxis const* ax2) {
0093     Double_t first_bin = ax1->GetBinWidth(1);
0094     Double_t last_bin = ax1->GetBinWidth(ax1->GetNbins());
0095     if (!TMath::AreEqualAbs(ax1->GetXmin(), ax2->GetXmin(), first_bin * 1e-10)
0096             || !TMath::AreEqualAbs(ax1->GetXmax(), ax2->GetXmax(), last_bin * 1e-10)) {
0097         throw DifferentAxisLimits();
0098         return false;
0099     }
0100     return true;
0101 }
0102 
0103 bool Hist4D::CheckBinLimits(TAxis const* ax1, TAxis const* ax2) {
0104     TArrayD const* h1arr = ax1->GetXbins();
0105     TArrayD const* h2arr = ax2->GetXbins();
0106     Int_t count = h1arr->GetSize();
0107     if (count != 0) {
0108         if (h2arr->GetSize() != count) {
0109             throw DifferentBinLimits();
0110             return false;
0111         } else {
0112             for (Int_t idx = 0; idx < count; ++idx) {
0113                 Double_t bin_width = ax1->GetBinWidth(idx);
0114                 if (!TMath::AreEqualAbs(h1arr->GetAt(idx), h2arr->GetAt(idx), bin_width * 1e-10)) {
0115                     throw DifferentBinLimits();
0116                     return false;
0117                 }
0118             }
0119         }
0120     }
0121     return true;
0122 }
0123 
0124 bool Hist4D::CheckBinLabels(TAxis const* ax1, TAxis const* ax2) {
0125     THashList* l1 = ax1->GetLabels();
0126     THashList* l2 = ax2->GetLabels();
0127     if (!l1 && !l2) {
0128         return true;
0129     } else if (!l1 || !l2) {
0130         throw DifferentLabels();
0131         return false;
0132     } else {
0133         for (Int_t idx = 1; idx <= ax1->GetNbins(); ++idx) {
0134             TString label1 = ax1->GetBinLabel(idx);
0135             TString label2 = ax2->GetBinLabel(idx);
0136             if (label1 != label2) {
0137                 throw DifferentLabels();
0138                 return false;
0139             }
0140         }
0141     }
0142     return true;
0143 }
0144 
0145 void Hist4D::CreateHists() {
0146     // Include overflow and underflow hists.
0147     Int_t nbins_w = _w_axis->GetNbins();
0148     Int_t nbins_x = _x_axis->GetNbins();
0149     Int_t nbins_y = _y_axis->GetNbins();
0150     Int_t nbins_z = _z_axis->GetNbins();
0151     for (Int_t idx_x = 0; idx_x < nbins_x + 2; ++idx_x) {
0152         for (Int_t idx_w = 0; idx_w < nbins_w + 2; ++idx_w) {
0153             std::string hist_name = std::string(fName.Data())
0154                 + "/hist" + std::to_string(idx_x) + ":" + std::to_string(idx_w);
0155             TArrayD const* bins_y = _y_axis->GetXbins();
0156             TArrayD const* bins_z = _z_axis->GetXbins();
0157             bool linear = (bins_y->GetSize() == 0) && (bins_z->GetSize() == 0);
0158             TH2D* hist;
0159             if (linear) {
0160                 hist = new TH2D(
0161                     hist_name.c_str(), "",
0162                     nbins_y, _y_axis->GetXmin(), _y_axis->GetXmax(),
0163                     nbins_z, _z_axis->GetXmin(), _z_axis->GetXmax());
0164             } else {
0165                 hist = new TH2D(
0166                     hist_name.c_str(), "",
0167                     bins_y->GetSize() - 1, bins_y->GetArray(),
0168                     bins_z->GetSize() - 1, bins_z->GetArray());
0169             }
0170             hist->GetXaxis()->SetTickLength(0.);
0171             hist->GetXaxis()->SetLabelOffset(999.);
0172             hist->GetXaxis()->SetLabelSize(0.);
0173             hist->GetYaxis()->SetTickLength(0.);
0174             hist->GetYaxis()->SetLabelOffset(999.);
0175             hist->GetYaxis()->SetLabelSize(0.);
0176             hist->SetStats(0);
0177             Int_t idx = (nbins_w + 2) * idx_x + idx_w;
0178             _hists.push_back(hist);
0179         }
0180     }
0181 }
0182 
0183 void Hist4D::CreateAxes() {
0184     std::string axes_hist_name = fName.Data();
0185     axes_hist_name += "/axes_hist";
0186     TArrayD const* bins_w = _w_axis->GetXbins();
0187     TArrayD const* bins_x = _x_axis->GetXbins();
0188     bool linear = (bins_w->GetSize() == 0) && (bins_x->GetSize() == 0);
0189     if (linear) {
0190         _axes = new TH2D(
0191             axes_hist_name.c_str(), "",
0192             _w_axis->GetNbins(), _w_axis->GetXmin(), _w_axis->GetXmax(),
0193             _x_axis->GetNbins(), _x_axis->GetXmin(), _x_axis->GetXmax());
0194     } else {
0195         _axes = new TH2D(
0196             axes_hist_name.c_str(), "",
0197             bins_w->GetSize() - 1, bins_w->GetArray(),
0198             bins_x->GetSize() - 1, bins_x->GetArray());
0199     }
0200     _axes->SetStats(0);
0201     _w_axis = _axes->GetXaxis();
0202     _x_axis = _axes->GetYaxis();
0203     _axes->GetXaxis()->SetTickLength(0.);
0204     //_axes->GetXaxis()->SetLabelOffset(999.);
0205     //_axes->GetXaxis()->SetLabelSize(0.);
0206     _axes->GetYaxis()->SetTickLength(0.);
0207     //_axes->GetYaxis()->SetLabelOffset(999.);
0208     //_axes->GetYaxis()->SetLabelSize(0.);
0209 }
0210 
0211 TObject* Hist4D::Clone(char const* new_name) const {
0212     Hist4D* result = new Hist4D();
0213     result->SetName(new_name);
0214     result->SetTitle(fTitle);
0215     result->_w_axis = _w_axis;
0216     result->_x_axis = _x_axis;
0217     result->_y_axis = _y_axis;
0218     result->_z_axis = _z_axis;
0219     result->_axes = static_cast<TH2D*>(_axes->Clone());
0220     for (TH2D* hist : _hists) {
0221         result->_hists.push_back(static_cast<TH2D*>(hist->Clone()));
0222     }
0223     CheckConsistency(this, result);
0224     return result;
0225 }
0226 
0227 void Hist4D::Fill(Double_t w, Double_t x, Double_t y, Double_t z) {
0228     Int_t idx_w = _w_axis->FindBin(w);
0229     Int_t idx_x = _x_axis->FindBin(x);
0230     Int_t nbins_w = _w_axis->GetNbins();
0231     Int_t nbins_x = _x_axis->GetNbins();
0232     Int_t idx = (nbins_w + 2) * idx_x + idx_w;
0233     _hists[idx]->Fill(y, z);
0234 }
0235 
0236 void Hist4D::Fill(Double_t w, Double_t x, Double_t y, Double_t z, Double_t weight) {
0237     Int_t idx_w = _w_axis->FindBin(w);
0238     Int_t idx_x = _x_axis->FindBin(x);
0239     Int_t nbins_w = _w_axis->GetNbins();
0240     Int_t nbins_x = _x_axis->GetNbins();
0241     Int_t idx = (nbins_w + 2) * idx_x + idx_w;
0242     _hists[idx]->Fill(y, z, weight);
0243 }
0244 
0245 Double_t Hist4D::GetEntries() {
0246     Double_t result = 0;
0247     for (TH2D* hist : _hists) {
0248         result += hist->GetEntries();
0249     }
0250     return result;
0251 }
0252 
0253 void Hist4D::Divide(Hist4D* other) {
0254     CheckConsistency(this, other);
0255     for (std::size_t idx = 0; idx < _hists.size(); ++idx) {
0256         _hists[idx]->Divide(other->_hists[idx]);
0257     }
0258 }
0259 
0260 void Hist4D::Scale(Double_t c) {
0261     for (TH2D* hist : _hists) {
0262         hist->Scale(c);
0263     }
0264 }
0265 
0266 TH2D* Hist4D::ProjectionYZ(char const* pname) {
0267     TArrayD const* bins_y = _y_axis->GetXbins();
0268     TArrayD const* bins_z = _z_axis->GetXbins();
0269     bool linear = (bins_y->GetSize() == 0) && (bins_z->GetSize() == 0);
0270     TH2D* hist_proj;
0271     if (linear) {
0272         hist_proj = new TH2D(
0273             pname, "",
0274             _y_axis->GetNbins(), _y_axis->GetXmin(), _y_axis->GetXmax(),
0275             _z_axis->GetNbins(), _z_axis->GetXmin(), _z_axis->GetXmax());
0276     } else {
0277         hist_proj = new TH2D(
0278             pname, "",
0279             bins_y->GetSize() - 1, bins_y->GetArray(),
0280             bins_z->GetSize() - 1, bins_z->GetArray());
0281     }
0282     for (TH2D* hist : _hists) {
0283         hist_proj->Add(hist);
0284     }
0285     return hist_proj;
0286 }
0287 
0288 TH2D* Hist4D::ProjectionWX(
0289         char const* pname,
0290         Int_t firstbiny, Int_t lastbiny,
0291         Int_t firstbinz, Int_t lastbinz) {
0292     TArrayD const* bins_w = _w_axis->GetXbins();
0293     TArrayD const* bins_x = _x_axis->GetXbins();
0294     Int_t nbins_w = _w_axis->GetNbins();
0295     Int_t nbins_x = _x_axis->GetNbins();
0296     bool linear = (bins_w->GetSize() == 0) && (bins_x->GetSize() == 0);
0297     TH2D* hist_proj;
0298     if (linear) {
0299         hist_proj = new TH2D(
0300             pname, "",
0301             _w_axis->GetNbins(), _w_axis->GetXmin(), _w_axis->GetXmax(),
0302             _x_axis->GetNbins(), _x_axis->GetXmin(), _x_axis->GetXmax());
0303     } else {
0304         hist_proj = new TH2D(
0305             pname, "",
0306             bins_w->GetSize() - 1, bins_w->GetArray(),
0307             bins_x->GetSize() - 1, bins_x->GetArray());
0308     }
0309     hist_proj->Sumw2(kFALSE);
0310     for (Int_t idx_x = 0; idx_x < nbins_x + 2; ++idx_x) {
0311         for (Int_t idx_w = 0; idx_w < nbins_w + 2; ++idx_w) {
0312             Int_t idx = (nbins_w + 2) * idx_x + idx_w;
0313             TH2D* hist = _hists[idx];
0314             Double_t error;
0315             Double_t integral = hist->IntegralAndError(
0316                 firstbiny, lastbiny, firstbinz, lastbinz, error);
0317             Int_t bin = hist_proj->GetBin(idx_w, idx_x);
0318             hist_proj->SetBinContent(bin, integral);
0319             hist_proj->SetBinError(bin, error);
0320         }
0321     }
0322     hist_proj->Sumw2(kTRUE);
0323     return hist_proj;
0324 }
0325 
0326 void Hist4D::DrawPad(TVirtualPad* pad, char const* options) {
0327     // Update the axes on all histograms.
0328     for (TH2D* hist : _hists) {
0329         TArrayD const* bins_y = _y_axis->GetXbins();
0330         TArrayD const* bins_z = _z_axis->GetXbins();
0331         bool linear = (bins_y->GetSize() == 0) && (bins_z->GetSize() == 0);
0332         if (linear) {
0333             hist->SetBins(
0334                 _y_axis->GetNbins(), _y_axis->GetXmin(), _y_axis->GetXmax(),
0335                 _z_axis->GetNbins(), _z_axis->GetXmin(), _z_axis->GetXmax());
0336         } else {
0337             hist->SetBins(
0338                 bins_y->GetSize() - 1, bins_y->GetArray(),
0339                 bins_z->GetSize() - 1, bins_z->GetArray());
0340         }
0341     }
0342     // Find minimum and maximum.
0343     Double_t min = _minimum;
0344     Double_t max = _maximum;
0345     if (TMath::IsNaN(min)) {
0346         min = TMath::Infinity();
0347         for (TH2D* hist : _hists) {
0348             Double_t next_min = hist->GetBinContent(hist->GetMinimumBin());
0349             if (next_min < min) {
0350                 min = next_min;
0351             }
0352         }
0353     }
0354     if (TMath::IsNaN(max)) {
0355         max = -TMath::Infinity();
0356         for (TH2D* hist : _hists) {
0357             Double_t next_max = hist->GetBinContent(hist->GetMaximumBin());
0358             if (next_max > max) {
0359                 max = next_max;
0360             }
0361         }
0362     }
0363     // Draw axes.
0364     pad->cd();
0365     _axes->SetTitle(fTitle);
0366     _axes->SetMinimum(min);
0367     _axes->SetMaximum(max);
0368     std::string axes_options = "colz ";
0369     axes_options += options;
0370     _axes->Draw(axes_options.c_str());
0371     pad->Paint();
0372     // Create sub-pads for individual plots.
0373     Int_t nbins_w = _w_axis->GetNbins();
0374     Int_t nbins_x = _x_axis->GetNbins();
0375     Int_t nbins_y = _y_axis->GetNbins();
0376     Int_t nbins_z = _z_axis->GetNbins();
0377     for (Int_t idx_x = 1; idx_x < nbins_x + 1; ++idx_x) {
0378         for (Int_t idx_w = 1; idx_w < nbins_w + 1; ++idx_w) {
0379             Int_t idx = (nbins_w + 2) * idx_x + idx_w;
0380             std::string pad_name = pad->GetName();
0381             pad_name += "/subpad" + std::to_string(idx_x) + ":" + std::to_string(idx_w);
0382             pad->cd();
0383             TPad* subpad = dynamic_cast<TPad*>(
0384                 pad->GetPrimitive(pad_name.c_str()));
0385             if (subpad == nullptr) {
0386                 Double_t width = pad->GetWw() * pad->GetAbsWNDC();
0387                 Double_t height = pad->GetWh() * pad->GetAbsHNDC();
0388                 Double_t pad_x_lo = pad->XtoPixel(pad->XtoPad(_w_axis->GetBinLowEdge(idx_w))) / width;
0389                 Double_t pad_y_lo = pad->YtoPixel(pad->YtoPad(_x_axis->GetBinLowEdge(idx_x))) / height;
0390                 Double_t pad_x_hi = pad->XtoPixel(pad->XtoPad(_w_axis->GetBinUpEdge(idx_w))) / width;
0391                 Double_t pad_y_hi = pad->YtoPixel(pad->YtoPad(_x_axis->GetBinUpEdge(idx_x))) / height;
0392                 Double_t pad_x_min = pad->XtoPixel(pad->GetUxmin()) / width;
0393                 Double_t pad_y_min = pad->YtoPixel(pad->GetUymin()) / height;
0394                 Double_t pad_x_max = pad->XtoPixel(pad->GetUxmax()) / width;
0395                 Double_t pad_y_max = pad->YtoPixel(pad->GetUymax()) / height;
0396                 if (!(pad_x_min < pad_x_max)) {
0397                     Double_t temp = pad_x_min;
0398                     pad_x_min = pad_x_max;
0399                     pad_x_max = temp;
0400                 }
0401                 if (!(pad_y_min < pad_y_max)) {
0402                     Double_t temp = pad_y_min;
0403                     pad_y_min = pad_y_max;
0404                     pad_y_max = temp;
0405                 }
0406                 if (!(pad_x_lo < pad_x_hi)) {
0407                     Double_t temp = pad_x_lo;
0408                     pad_x_lo = pad_x_hi;
0409                     pad_x_hi = temp;
0410                 }
0411                 if (!(pad_y_lo < pad_y_hi)) {
0412                     Double_t temp = pad_y_lo;
0413                     pad_y_lo = pad_y_hi;
0414                     pad_y_hi = temp;
0415                 }
0416                 if (pad_x_lo >= pad_x_min
0417                         && pad_y_lo >= pad_y_min
0418                         && pad_x_hi <= pad_x_max
0419                         && pad_y_hi <= pad_y_max) {
0420                     subpad = new TPad(
0421                         pad_name.c_str(), "",
0422                         pad_x_lo, 1. - pad_y_hi, pad_x_hi, 1. - pad_y_lo);
0423                     subpad->SetLeftMargin(0.01);
0424                     subpad->SetRightMargin(0.01);
0425                     subpad->SetTopMargin(0.01);
0426                     subpad->SetBottomMargin(0.01);
0427                     subpad->SetFillStyle(4000);
0428                     subpad->SetFrameFillStyle(4000);
0429                     subpad->SetFrameLineColor(0);
0430                     subpad->SetFrameLineWidth(0);
0431                     subpad->SetBorderMode(0);
0432                     if (pad->GetLogz()) {
0433                         subpad->SetLogz();
0434                     }
0435                 }
0436             }
0437             if (subpad != nullptr) {
0438                 pad->cd();
0439                 subpad->Draw();
0440                 subpad->cd();
0441                 std::string hist_options = "cola ";
0442                 hist_options += options;
0443                 _hists[idx]->SetMinimum(min);
0444                 _hists[idx]->SetMaximum(max);
0445                 _hists[idx]->Draw(hist_options.c_str());
0446             }
0447         }
0448     }
0449     pad->cd();
0450 }
0451 
0452 Hist4D::~Hist4D() {
0453   if(_w_axis) delete _w_axis;
0454   if(_x_axis) delete _x_axis;
0455   if(_y_axis) delete _y_axis;
0456   if(_z_axis) delete _z_axis;
0457   if(_axes) delete _axes;
0458   for(auto e : _hists)
0459     if(e) delete e;
0460   _hists.clear();
0461 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!