HistoView.cc

Go to the documentation of this file.

#include "AddText.h"
#include "HistoView.h"

#include <utl/Histogram.h>
#include <utl/Math.h>
#include <sdet/Station.h>
#include <sevt/Station.h>
#include <sevt/PMT.h>
#include <sevt/PMTCalibData.h>
#include <utl/QuadraticFitData.h>
#include <utl/ExponentialFitData.h>

#include <TCanvas.h>
#include <TH1D.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TPaveText.h>
#include <TArrow.h>

#include <vector>
#include <iomanip>

using namespace utl;
using namespace sevt;
using namespace std;
using namespace SdCalibPlotterOG;


void
MuonPeakView::Draw(const sdet::Station& dStation, const sevt::Station& sStation, const PMT& pmt, const int pmtId)
{
  if (HasNoCalib(pmt))
    return;
  ostringstream os;
  os << "MuPeak PMT " << pmtId;
  const auto& peakBinning = dStation.GetMuonPeakHistogramBinning<double>(pmt.GetType(), sStation.GetCalibData().GetVersion());
  auto& peak = Add(new TH1D(os.str().c_str(), os.str().c_str(), peakBinning.size()-1, &peakBinning[0]));
  const auto& pmtCalib = pmt.GetCalibData();
  SetContent(peak, peakBinning, pmtCalib.GetMuonPeakHisto(), "peak");
  peak.SetLineColor(kRed);
  peak.SetStats(false);
  peak.Draw();
  const int offset = pmtCalib.GetMuonPeakHistoOffset();
  const double lsValue = pmtCalib.GetOnlinePeak() + pmtCalib.GetBaseline() - offset;
  if (pmt.HasRecData())
    DrawQuadraticFit(pmt.GetRecData().GetMuonPeakFitData(), peak, lsValue);
  else
    DrawQuadraticFit(peak, lsValue);
  const double peakMax = peak.GetMaximum();
  DrawBaseline(pmtCalib, peakMax);
  DrawOffset(offset, peakMax);
  IsNotOK(pmt);
}


void
MuonChargeView::Draw(const sdet::Station& dStation, const sevt::Station& sStation, const PMT& pmt, const int pmtId)
{
  if (HasNoCalib(pmt))
    return;
  ostringstream os;
  os << "MuCharge PMT" << pmtId;
  const auto& chargeBinning = dStation.GetMuonChargeHistogramBinning<double>(pmt.GetType(), sStation.GetCalibData().GetVersion());
  const auto& pmtCalib = pmt.GetCalibData();
  auto& charge = Add(new TH1D(os.str().c_str(), os.str().c_str(), chargeBinning.size()-1, &chargeBinning[0]));
  SetContent(charge, chargeBinning, pmtCalib.GetMuonChargeHisto(), "charge");
  charge.SetLineColor(kRed);
  charge.SetStats(false);
  charge.Draw();
  const int offset = pmtCalib.GetMuonChargeHistoOffset();
  const int baseFact = pmtCalib.GetMuonShapeHisto().size();
  const double lsValue = pmtCalib.GetOnlineCharge()*1.045 + pmtCalib.GetBaseline() * baseFact - offset;
  //                                                ^^^^^ hump magic
  if (pmt.HasRecData()) {
    DrawQuadraticFit(pmt.GetRecData().GetMuonChargeFitData(), charge, lsValue);
    DrawExponentialFit(charge, pmt.GetRecData().GetMuonChargeSlopeFitData());
  } else
    DrawQuadraticFit(charge, lsValue);
  const double chargeMax = charge.GetMaximum();
  DrawBaseline(pmtCalib, chargeMax);
  DrawOffset(double(offset) / baseFact, chargeMax);
  IsNotOK(pmt);
}


void
MuonShapeView::Draw(const sdet::Station& dStation, const PMT& pmt, const int pmtId, const int calibVersion)
{
  if (HasNoCalib(pmt))
    return;
  gPad->SetLogy();
  ostringstream os;
  os << "MuShape PMT" << pmtId;
  const auto& shapeBinning = dStation.GetMuonShapeHistogramBinning(calibVersion);
  auto& shape = Add(new TH1D(os.str().c_str(), os.str().c_str(), shapeBinning.size()-1, &shapeBinning[0]));
  const auto& pmtCalib = pmt.GetCalibData();
  SetContent(shape, shapeBinning, pmtCalib.GetMuonShapeHisto(), "shape", 0);
  shape.SetLineColor(kRed);
  shape.SetStats(false);
  shape.Draw();
  if (pmt.HasRecData()) {
    const auto& pmtRec = pmt.GetRecData();
    DrawExponentialFit(shape, pmtRec.GetMuonShapeFitData());
    DrawDecay(pmtRec.GetMuonPulseDecayTime());
  }
  IsNotOK(pmt);
}


void
PMTBaselineView::Draw(const sdet::Station& dStation, const PMT& pmt, const int pmtId)
{
  if (HasNoCalib(pmt))
    return;
  gPad->SetLogy();
  ostringstream os;
  os << "MuBase PMT" << pmtId;
  const auto& baseBinning = dStation.GetBaselineHistogramBinning<double>();
  auto& base = Add(new TH1D(os.str().c_str(), os.str().c_str(), baseBinning.size()-1, &baseBinning[0]));
  const auto& pmtCalib = pmt.GetCalibData();
  SetContent(base, baseBinning, pmtCalib.GetMuonBaseHisto(), "base");
  base.SetLineColor(kRed);
  base.SetStats(false);
  base.Draw();
  DrawAverage(base);
  IsNotOK(pmt);
}


void
HistoView::DividePage(TCanvas& c, const double xmargin, const double ymargin)
{
  const auto savePad = (TPad*)gPad;
  c.cd();
  const string name = c.GetName();
  ostringstream os;
  const int nRows = 6;
  int n = 1;
  for (int row = 0; row < nRows; ++row) {
    const double y1 = 1 - (row+1)*(1 - ymargin)/nRows;
    const double y2 = 1 - ymargin - row*(1 - ymargin)/nRows;
    const int nColumns = (row ? (row < 4 ? 2 : 3) : 1);
    for (int column = 0; column < nColumns; ++column) {
      const double x1 = xmargin + column*(1 - xmargin)/nColumns;
      const double x2 = (column+1)*(1 - xmargin)/nColumns;
      os.str(""); os << name << '_' << row+1 << column+1;
      const auto pad = new TPad(os.str().c_str(), os.str().c_str(), x1, y1, x2, y2, 0);
      pad->SetNumber(n++);
      pad->Draw();
    }
  }
  savePad->cd();
}


bool
HistoView::HasNoCalib(const PMT& pmt)
{
  if (!pmt.HasCalibData()) {
    auto& pt = Add(new TPaveText(0.3,0.4,0.7,0.6, "ndc"));
    pt.SetBorderSize(2);
    pt.SetFillColor(kRed);
    pt.SetTextAlign(12);
    pt.AddText("No Calibration Data");
    pt.Draw();
    return true;
  }
  return false;
}


bool
HistoView::IsNotOK(const PMT& pmt)
{
  if (!pmt.GetCalibData().IsTubeOk()) {
    auto& pt = Add(new TPaveText(0.5,0.43,0.8,0.57, "ndc"));
    pt.SetBorderSize(2);
    pt.SetFillColor(kRed);
    pt.SetTextAlign(12);
    pt.AddText("Tube not OK");
    pt.Draw();
    return true;
  }
  return false;
}


void
HistoView::SetContent(TH1D& rh,
                      const vector<double>& bins, const vector<int>& hist,
                      const string& what,
                      const int subtractBin)
  const
{
  if (bins.size()-1 != hist.size()) {
    ostringstream warn;
    warn << "Mismatch in PMT " << what << " histogram: " << bins.size()-1 << " bin intervals "
            "but " << hist.size() << " values!";
    WARNING(warn);
    return;
  }

  const VariableBinHistogramWrap<double, int> h(bins, hist);
  const double sub = (subtractBin >= 0 ? h.GetBinAverage(subtractBin) : 0);
  const int n = h.GetNBins();
  for (int i = 0; i < n; ++i) {
    const double val = h.GetBinAverage(i) - sub;
    rh.SetBinContent(i+1, (val > 0 ? val : 0));
  }
}


void
HistoView::DrawQuadraticFit(const QuadraticFitData& qf, TH1D& hist,
                            const double lsValue)
{
  const string parabolaName = string(hist.GetName()) + "_parabola";
  auto& par = Add(new TF1(parabolaName.c_str(),
                          [&qf](double* x, double*) { return qf(*x); }, qf.GetStart(), qf.GetStop(), 0));
  par.SetLineWidth(1);
  par.SetLineColor(kBlue);
  par.Draw("same");
  const double x = qf.GetExtremePosition();
  const double yext = qf(x);
  const double y = 1.15 * yext;
  const double dx = qf.GetExtremePositionError();
  // fit value with errorbar
  auto& gfit = Add(new TGraphErrors(1, &x, &y, &dx, 0));
  gfit.SetLineColor(kBlue);
  gfit.SetMarkerStyle(kFullCircle);
  gfit.SetMarkerColor(kBlue);
  gfit.SetMarkerSize(0.5);
  gfit.Draw("p same");
  const double yy = 1.1 * y;
  const double histMax = hist.GetMaximum();
  if (yy > histMax)
    hist.SetMaximum(yy);
  DrawQuadraticFit(hist, lsValue, y);
}


void
HistoView::DrawQuadraticFit(TH1D& hist, const double lsValue, const double fitValue)
{
  const double hMax = hist.GetMaximum();
  const double hValue = hist.GetBinContent(hist.FindBin(lsValue));
  const bool up = (hValue > 0.5*hMax);
  const double y2 = up ? 0.82*hValue : max(hValue + 0.15*hMax, 1.15*fitValue);
  const double y1 = up ? 0.5*y2 : y2 + 0.5*(hMax-y2);
  // ls reported value
  auto& als = Add(new TArrow(lsValue, y1, lsValue, y2, 0.01));
  als.SetAngle(20);
  als.SetLineColor(8);
  als.Draw();
}


void
HistoView::DrawExponentialFit(TH1D& hist, const ExponentialFitData& ef)
{
  const string exponentialName = string(hist.GetName()) + "_exponential";
  auto& ex = Add(new TF1(exponentialName.c_str(),
                         [&ef](double* x, double*) { return ef(*x); }, ef.GetStart(), ef.GetStop(), 0));
  ex.SetParameter(0, ef.GetAmplitude());
  ex.SetParameter(1, ef.GetSlope());
  ex.SetLineWidth(1);
  ex.SetLineColor(kBlue);
  ex.Draw("same");
}


void
HistoView::DrawBaseline(const PMTCalibData& pmtCalib, const double max)
{
  const double x = pmtCalib.GetBaseline();
  const double y = 0.19*max;
  const double dx = pmtCalib.GetBaselineRMS();
  auto& g = Add(new TGraphErrors(1, &x, &y, &dx, 0));
  g.SetLineColor(46);
  g.SetMarkerColor(46);
  g.Draw("same");
}


void
HistoView::DrawOffset(const double x, const double max)
{
  const double y = 0.15*max;
  auto& a = Add(new TArrow(x, y, x, 0, 0.01));
  a.SetAngle(20);
  a.SetLineColor(kBlack);
  a.Draw();
}


void
HistoView::DrawAverage(TH1D& h)
{
  const double x = h.GetMean();
  const double dx = sqrt(Sqr(h.GetMeanError()) + Sqr(h.GetBinWidth(1))/12);
  const double y = 0.01*h.GetMaximum();
  auto& g = Add(new TGraphErrors(1, &x, &y, &dx, &y));
  g.SetLineColor(kBlack);
  g.SetMarkerColor(kBlack);
  g.DrawClone("same");
  const double y2 = 0.1*y;
  const double dx2 = h.GetRMS();
  auto& g2 = Add(new TGraphErrors(1, &x, &y2, &dx2, 0));
  g2.SetLineColor(kBlack);
  g2.SetMarkerColor(kBlack);
  g2.Draw("same");
}


void
HistoView::DrawDecay(const double time)
{
  ostringstream os;
  os << "#tau = " << setprecision(4) << time/nanosecond << " ns";
  auto& pt = Add(new TPaveText(0.58,0.77,0.87,0.85, "ndc"));
  pt.SetBorderSize(1);
  pt.SetFillColor(kWhite);
  pt.SetTextAlign(22);
  pt.AddText(os.str().c_str());
  pt.Draw();
}