VEMTraceView.cc

Go to the documentation of this file.

#include <algorithm>
#include <sstream>

#include <boost/format.hpp>

#include <TCanvas.h>
#include <TH1D.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TArrow.h>
#include <TBox.h>

#include <sevt/PMT.h>
#include <sevt/Station.h>

#include "AddText.h"
#include "VEMTraceView.h"

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


void
VEMTraceView::DividePage(TCanvas& c)
{
  const double xmargin = 0.003;
  const double ymargin = 0.003;
  const double vertDivision = 0.2;

  TPad* const savePad = (TPad*)gPad;
  c.cd();
  const string name = c.GetName();
  ostringstream os;
  int n = 1;
  os.str(""); os << name << "_1";
  const double x1 = xmargin;
  const double x2 = vertDivision - 0.5*xmargin;
  const double x3 = vertDivision + 0.5*xmargin;
  const double x4 = 1 - xmargin;
  const double y11 = ymargin;
  const double y12 = 1 - ymargin;
  TPad* const pad1 = new TPad(os.str().c_str(), os.str().c_str(), x1, y11, x2, y12, 0);
  pad1->SetNumber(n++);
  pad1->Draw();
  const int nRows = 4;
  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;
    os.str(""); os << name << '_' << n;
    TPad* const pad = new TPad(os.str().c_str(), os.str().c_str(), x3, y1, x4, y2, 0);
    pad->SetNumber(n++);
    pad->SetLeftMargin(0.5*pad->GetLeftMargin());
    pad->SetRightMargin(0);
    pad->Draw();
  }
  savePad->cd();
}


void
PMTVEMTraceView::Draw(const PMT& pmt, const int stationStart, const int stationStop)
{
  if (!pmt.HasRecData()) {
    Missing("No Rec Data");
    return;
  }

  const PMTRecData& pmtRec = pmt.GetRecData();

  if (!pmtRec.HasVEMTrace()) {
    Missing("No VEM Trace");
    return;
  }

  // vem trace
  const TraceD& trace = pmtRec.GetVEMTrace();
  const int traceLength = trace.GetSize();
  ostringstream os;
  os << "PMT" << pmt.GetId() << "_VEMTrace";
  TH1D& th =
    Add(new TH1D(os.str().c_str(), os.str().c_str(), traceLength, 0, traceLength));
  for (int i = 0; i < traceLength; ++i)
    th.SetBinContent(i+1, trace[i]);
  th.SetStats(false);
  th.SetLineColor(kRed);

  // peak amplitude
  const double peak = pmtRec.GetPeakAmplitude();
  const double xPeak[2] = { 0, double(traceLength) };
  const double yPeak[2] = { peak, peak };
  TGraph& pg = Add(new TGraph(2, xPeak, yPeak));
  pg.SetTitle(os.str().c_str());
  pg.SetLineStyle(kDashed);
  pg.SetLineColor(17); // gray
  pg.Draw("AL");
  pg.GetXaxis()->SetLimits(0, traceLength);
  pg.GetYaxis()->SetTitle("[VEM peak]");
  pg.GetYaxis()->SetTitleOffset(0.3);
  pg.GetYaxis()->SetTickLength(0.01);
  TH1F& pgh = *pg.GetHistogram();
  pgh.SetMinimum(th.GetMinimum());
  pgh.SetMaximum(th.GetMaximum());

  // station start/stop
  if (stationStart >= 0) {
    const double x1[2] = { double(stationStart), double(stationStart) };
    const double x2[2] = { double(stationStop), double(stationStop) };
    const double y[2] = { -peak, 3*peak };
    TGraph& ss1 = Add(new TGraph(2, x1, y));
    ss1.SetLineColor(kBlue);
    ss1.SetLineStyle(kDashed);
    ss1.Draw("L same");
    TGraph& ss2 = Add(new TGraph(2, x2, y));
    ss2.SetLineColor(kBlue);
    ss2.SetLineStyle(kDashed);
    ss2.Draw("L same");
  }

  // zero line
  const double yZero[2] = { 0, 0 };
  TGraph& zg = Add(new TGraph(2, xPeak, yZero));
  zg.SetLineWidth(1);
  zg.SetLineColor(17); // gray
  zg.Draw("L same");

  // signals
  const SignalSegmentCollection& signals = pmtRec.GetSignals();
  {
    for (unsigned int s = 0; s < signals.size(); ++s) {
      double min = trace[signals[s].fStart];
      double max = min;
      for (int i = signals[s].fStart+1; i < signals[s].fStop; ++i) {
        double val = trace[i];
        if (val < min)
          min = val;
        if (val > max)
          max = val;
      }
      const double eps = 0.1*(max - min);
      const double y1 = min - eps;
      const double y2 = max + eps;
      if (y1 < pgh.GetMinimum())
        pgh.SetMinimum(y1);
      if (y2 > pgh.GetMaximum())
        pgh.SetMaximum(y2);
    }
    const double plotSpan = pgh.GetMaximum() - pgh.GetMinimum();
    for (unsigned int s = 0; s < signals.size(); ++s) {
      double min = trace[signals[s].fStart];
      double max = min;
      for (int i = signals[s].fStart; i < signals[s].fStop; ++i) {
        double val = trace[i];
        if (val < min)
          min = val;
        if (val > max)
          max = val;
      }
      const double diff = max - min;
      const double eps = 0.05*diff;
      const double y1 = min - eps;
      const double y2 = max + eps;
      TBox& sb = Add(new TBox(signals[s].fStart, y1, signals[s].fStop, y2));
      sb.SetFillStyle(0);
      sb.SetLineColor(8); // dark green
      sb.Draw();
      const string& chargeStr = (format("%.2f") % signals[s].fCharge).str();
      if (diff < 0.35*plotSpan) {
        const double x = 0.5*(signals[s].fStart + signals[s].fStop);
        TArrow& a = Add(new TArrow(x, max + 0.3*plotSpan, x, max + 0.05*plotSpan, 0.01));
        a.SetAngle(20);
        a.SetLineColor(8); // dark green
        a.Draw();
        TText& t = Add(new TText(signals[s].fStart, max + 0.34*plotSpan, chargeStr.c_str()));
        t.Draw();
      } else {
        TText& t = Add(new TText(signals[s].fStop + 3, 0.3*min + 0.7*max, chargeStr.c_str()));
        t.Draw();
      }
    }
  }

  th.Draw("same");

  const sdet::PMTConstants::PMTGain gainUsed = pmtRec.GetGainUsed();

  // gain used
  {
    TPaveText& pt = Add(new TPaveText(0.88,0.77,0.99,0.87, "ndc"));
    pt.SetBorderSize(2);
    pt.SetTextAlign(12);
    if (gainUsed == sdet::PMTConstants::eHighGain) {
      pt.SetFillColor(kWhite);
      pt.AddText("From High Gain");
    } else {
      pt.SetFillColor(kCyan);
      pt.AddText("From Low Gain");
    }
    pt.Draw();
  }

  // sturated bins
  {
    const int satBins = pmtRec.GetFADCSaturatedBins(gainUsed);
    if (satBins) {
      TPaveText& pt = Add(new TPaveText(0.92,0.64,0.99,0.74, "ndc"));
      pt.SetBorderSize(2);
      pt.SetTextAlign(12);
      pt.SetFillColor(kCyan);
      pt.AddText((format("N_{sat} = %d ") % satBins).str().c_str());
      pt.Draw();
    }
  }

  // total charge
  {
    TPaveText& pt = Add(new TPaveText(0.9,0.90,0.99,1, "ndc"));
    pt.SetBorderSize(2);
    pt.SetTextAlign(12);
    pt.SetFillColor(kWhite);
    pt.AddText((format("q = %.2f VEM") % pmtRec.GetTotalCharge()).str().c_str());
    pt.Draw();
  }
}


void
StationVEMTraceView::Draw(const Station& station)
{
  if (!station.HasVEMTrace()) {
    Missing("No VEM Trace");
    return;
  }

  const TraceD& trace = station.GetVEMTrace();
  const int traceLength = trace.GetSize();
  const string titleStr = (format("Station_%04d_VEMTrace") % station.GetId()).str();

  // zero line
  const double xZero[2] = { 0., double(traceLength) };
  const double yZero[2] = { 0, 0 };
  TGraph& zg = Add(new TGraph(2, xZero, yZero));
  zg.SetLineWidth(1);
  zg.SetLineColor(17); // gray
  zg.SetTitle(titleStr.c_str());
  zg.Draw("AL");
  zg.GetXaxis()->SetLimits(0, traceLength);
  zg.GetYaxis()->SetTitle("[VEM peak]");
  zg.GetYaxis()->SetTitleOffset(0.3);
  zg.GetYaxis()->SetTickLength(0.01);

  TH1D& th =
    Add(new TH1D(titleStr.c_str(), "", traceLength, 0, traceLength));
  for (int i = 0; i < traceLength; ++i)
    th.SetBinContent(i+1, trace[i]);
  th.SetStats(false);
  th.SetLineColor(kRed);

  TH1F& zgh = *zg.GetHistogram();
  zgh.SetMinimum(th.GetMinimum());
  zgh.SetMaximum(th.GetMaximum());

  // signals
  const SignalSegmentCollection& signals = station.GetSignals();
  {
    for (unsigned int s = 0; s < signals.size(); ++s) {
      double min = trace[signals[s].fStart];
      double max = min;
      for (int i = signals[s].fStart+1; i < signals[s].fStop; ++i) {
        double val = trace[i];
        if (val < min)
          min = val;
        if (val > max)
          max = val;
      }
      const double eps = 0.1*(max - min);
      const double y1 = min - eps;
      const double y2 = max + eps;
      if (y1 < zgh.GetMinimum())
        zgh.SetMinimum(y1);
      if (y2 > zgh.GetMaximum())
        zgh.SetMaximum(y2);
    }
    const double plotSpan = zgh.GetMaximum() - zgh.GetMinimum();
    for (unsigned int s = 0; s < signals.size(); ++s) {
      double min = trace[signals[s].fStart];
      double max = min;
      for (int i = signals[s].fStart; i < signals[s].fStop; ++i) {
        double val = trace[i];
        if (val < min)
          min = val;
        if (val > max)
          max = val;
      }
      const double diff = max - min;
      const double eps = 0.05*diff;
      const double y1 = min - eps;
      const double y2 = max + eps;
      TBox& sb = Add(new TBox(signals[s].fStart, y1, signals[s].fStop, y2));
      sb.SetFillStyle(0);
      sb.SetLineColor(8); // dark green
      sb.Draw();
      const string& chargeStr = (format("%.2f") % signals[s].fCharge).str();
      if (diff < 0.35*plotSpan) {
        const double x = 0.5*(signals[s].fStart + signals[s].fStop);
        TArrow& a = Add(new TArrow(x, max + 0.3*plotSpan, x, max + 0.05*plotSpan, 0.01));
        a.SetAngle(20);
        a.SetLineColor(8); // dark green
        a.Draw();
        TText& t = Add(new TText(signals[s].fStart, max + 0.34*plotSpan, chargeStr.c_str()));
        t.Draw();
      } else {
        TText& t = Add(new TText(signals[s].fStop + 3, 0.3*min + 0.7*max, chargeStr.c_str()));
        t.Draw();
      }
    }
  }

  if (!station.HasRecData()) {
    Missing("No Rec Data");
    th.Draw("same");
    return;
  }

  // final signal start/end time
  const StationRecData& stRec = station.GetRecData();
  const int start = stRec.GetSignalStartSlot();
  const int stop = stRec.GetSignalEndSlot();
  double min = trace[start];
  double max = min;
  for (int i = start; i < stop; ++i) {
    double val = trace[i];
    if (val < min)
      min = val;
    if (val > max)
      max = val;
  }
  const double diff = max - min;
  const double eps = 0.05*diff;
  const double y1 = min - eps;
  const double y2 = max + eps;
  TBox& sb = Add(new TBox(start, y1, stop, y2));
  sb.SetFillStyle(0);
  sb.SetLineColor(kBlue);
  sb.Draw();

  th.Draw("same");

  // total signal
  {
    TPaveText& pt = Add(new TPaveText(0.9,0.90,0.99,1, "ndc"));
    pt.SetBorderSize(2);
    pt.SetTextAlign(12);
    pt.SetFillColor(kWhite);
    pt.AddText((format("q = %.2f VEM") % stRec.GetTotalSignal()).str().c_str());
    pt.Draw();
  }
}


void
VEMTraceView::Missing(const string& what)
{
  TPaveText& pt = Add(new TPaveText(0.3,0.4,0.7,0.6, "ndc"));
  pt.SetBorderSize(2);
  pt.SetFillColor(kRed);
  pt.SetTextAlign(12);
  pt.AddText(what.c_str());
  pt.Draw();
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: