SdRecPlotter.cc

Go to the documentation of this file.

#include <vector>
#include <sstream>
#include <iomanip>

#include <boost/tokenizer.hpp>

#include <TROOT.h>
#include <TStyle.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <TEllipse.h>
#include <TLatex.h>
#include <TLine.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TGraphBentErrors.h>
#include <TMultiGraph.h>
#include <TPaveText.h>
#include <TPolyLine.h>

#include <fwk/CentralConfig.h>

#include <det/Detector.h>

#include <sdet/SDetector.h>
#include <sdet/STimeVariance.h>

#include <evt/ShowerSimData.h>
#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <sevt/SEvent.h>
#include <sevt/Header.h>
#include <sevt/Station.h>
#include <sevt/StationRecData.h>
#include <sevt/StationTriggerData.h>

#include <utl/PhysicalConstants.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/TraceAlgorithm.h>
#include <utl/AxialVector.h>
#include <utl/UTMPoint.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/SaveCurrentTDirectory.h>

#include "SdRecPlotter.h"

using namespace SdRecPlotterOG;
using namespace sdet;
using namespace fwk;
using namespace evt;
using namespace sevt;
using namespace utl;
using namespace std;


namespace SdRecPlotterOG {

  // perpendicular distance from the axis
  inline
  double
  RPerp(const Vector& axis, const Vector& station)
  {
    const double scal = axis*station;
    return sqrt(station*station - scal*scal);
  }


  inline
  string
  GetT4TriggerName(const int t4)
  {
    ostringstream trig;
    if (t4) {
      if (t4 & ShowerSRecData::eT4_FD)
        trig << "FD";
      if (t4 & ShowerSRecData::eT4_3TOT)
        trig << (trig.str().length() ? "+" : "") << "3TOT";
      if (t4 & ShowerSRecData::eT4_4C1)
        trig << (trig.str().length() ? "+" : "") << "4C1";
    } else
      trig << "None";
    return trig.str();
  }


  inline
  string
  GetT5TriggerName(const int t5)
  {
    ostringstream trig;
    if (t5) {
      if (t5 & ShowerSRecData::eT5_Prior)
        trig << "Prior";
      if (t5 & ShowerSRecData::eT5_Posterior)
        trig << (trig.str().length() ? "+" : "") << "Posterior";
    } else
      trig << "None";
    return trig.str();
  }

  void
  AddText(TPaveText& pt, const string& str)
  {
    typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
    boost::char_separator<char> separator("\n");
    tokenizer tokens(str, separator);
    for (tokenizer::iterator it = tokens.begin();
         it != tokens.end(); ++it)
      pt.AddText(it->c_str());
  }

}


VModule::ResultFlag
SdRecPlotter::Init()
{
  CentralConfig* cc = CentralConfig::GetInstance();
  Branch topBranch = cc->GetTopBranch("SdRecPlotter");

  topBranch.GetChild("FileNamePrefix").GetData(fFileNamePrefix);

  fGenerateROOTFile = bool(topBranch.GetChild("GenerateROOTFile"));

  fGeneratePDFFile = bool(topBranch.GetChild("GeneratePDFFile"));

  fAppendOutput = bool(topBranch.GetChild("AppendOutput"));

  if (!fGenerateROOTFile && !fGeneratePDFFile)
    WARNING("Both ROOT and PS options switched off, no output will be generated.");

  if (fGenerateROOTFile) {
    const string fileName = fFileNamePrefix + ".root";
    fROOTFile = new TFile(fileName.c_str(), "recreate", "SdRecPlotter", 9);
    INFO(fileName + " file created.");
  }

  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1, 0);

  if (fGeneratePDFFile) {
    gStyle->SetPaperSize(TStyle::kA4);
    gStyle->SetTitlePS(fFileNamePrefix.c_str());
    if (fAppendOutput)
      fPDFFileName = fFileNamePrefix + ".pdf";
  }

  return eSuccess;
}


VModule::ResultFlag
SdRecPlotter::Run(Event& event)
{
  if (!fGenerateROOTFile && !fGeneratePDFFile)
    return eSuccess;
  if (!event.HasRecShower()) {
    WARNING("No reconstructed shower data found.");
    return eSuccess;
  }
  if (!event.GetRecShower().HasSRecShower()) {
    WARNING("No SD reconstructed shower data found.");
    return eSuccess;
  }

  if (fGeneratePDFFile && !fAppendOutput) {
    ostringstream os;
    os << fFileNamePrefix << '_' << event.GetSEvent().GetHeader().GetId() << ".pdf";
    fPDFFileName = os.str();
  }

  DrawEvent(event);

  return eSuccess;
}


VModule::ResultFlag
SdRecPlotter::Finish()
{
  if (fROOTFile) {
    fROOTFile->Close();
    delete fROOTFile;
  }
  if (fGeneratePDFFile && fAppendOutput) {
    TCanvas c;
    c.Print((fPDFFileName + ")").c_str(), "pdf");
  }
  return eSuccess;
}


//=================================================================================================

namespace SdRecPlotterOG {

  class CandidateInfo {
  public:
    CandidateInfo() { }
    ~CandidateInfo();
    void PushStation(const int id, const double x, const double y,
                     const double r, const double sig, const double tCurv,
                     const double t0, const double sigma_t,
                     const double t10, const double t50, const double t90,
                     const int groupSize, const int saturation, const bool tot);
    void SetSeed();
    void CalculateLimits(const double zoom = 1);
    double GetMinX() const { return fMinX; }
    double GetMaxX() const { return fMaxX; }
    double GetMinY() const { return fMinY; }
    double GetMaxY() const { return fMaxY; }
    double GetMinR() const { return fMinR; }
    double GetMaxR() const { return fMaxR; }
    double GetMaxSignal() const { return fMaxSignal; }
    double GetMinTime() const { return fMinTime; }
    double GetMaxTime() const { return fMaxTime; }
    void DrawArray();
    TGraph* GetLDFGraph();
    TGraphBentErrors* GetSaturatedLDFGraph();
    TGraph* GetStartTimeGraph();
    TGraphBentErrors* GetRiseFallTimeGraph();
    TGraph* GetCurvatureFitGraph();

    static const double kAspect;
    static const double kLDFMinR;
    static const double kStationMinRadius;
    static const double kStationMaxRadius;
    static const double kMinDX;
    static const double kTwinSplit;
    static const double kTwinIdSplit;
    static const double kLowGainSatFactor;

  private:
    vector<int> fId;
    vector<double> fPositionX;
    vector<double> fPositionY;
    double fMinX = 0;
    double fMaxX = 0;
    double fMinY = 0;
    double fMaxY = 0;
    vector<double> fR;
    double fMinR = 1;
    double fMaxR = 1;
    vector<double> fSignal;
    double fMaxSignal = 0;
    vector<double> fPlotRadius;
    vector<double> fCurvTime;
    vector<double> fStartTime;
    vector<double> fStartTimeErr;
    double fMinStartTime = 0;
    double fMaxStartTime = 0;
    vector<double> fTime10;
    vector<double> fTime50;
    vector<double> fTime90;
    double fMinTime = 0;
    double fMaxTime = 0;
    vector<int> fTwin;
    vector<int> fSaturation;
    vector<bool> fTOT;
    // deletables
    vector<TEllipse*> fCircles;
    TGraph* fLDFGraph = nullptr;
    TGraphBentErrors* fSaturatedLDFGraph = nullptr;
    TGraphErrors* fLDFNonsaturatedGraph = nullptr;
    TGraph* fStartTimeGraph = nullptr;
    TGraphBentErrors* fRiseFallTimeGraph = nullptr;
    TGraph* fCurvatureFitGraph = nullptr;
  };


  const double CandidateInfo::kAspect = sqrt(2.);
  const double CandidateInfo::kLDFMinR = 0.2; //km
  const double CandidateInfo::kStationMinRadius = 0.07; //km
  const double CandidateInfo::kStationMaxRadius = 0.7; //km
  const double CandidateInfo::kMinDX = 7; //km
  const double CandidateInfo::kTwinSplit = 0.07; //km
  const double CandidateInfo::kTwinIdSplit = 0.3; //km
  const double CandidateInfo::kLowGainSatFactor = 0.8;

}


CandidateInfo::~CandidateInfo()
{
  for (const auto c : fCircles)
    delete c;
  delete fLDFGraph;
  delete fSaturatedLDFGraph;
  delete fLDFNonsaturatedGraph;
  delete fStartTimeGraph;
  delete fRiseFallTimeGraph;
  delete fCurvatureFitGraph;
}


void
CandidateInfo::PushStation(const int id, const double x, const double y,
                           const double r, const double sig, const double tCurv,
                           const double t0, const double sigma_t,
                           const double t10, const double t50, const double t90,
                           const int groupSize, const int sat, const bool tot)
{
  fId.push_back(id);
  fPositionX.push_back(x);
  fPositionY.push_back(y);
  fR.push_back(r);
  if (r < fMinR)
    fMinR = r;
  if (r > fMaxR)
    fMaxR = r;
  fSignal.push_back(sig);
  if (sig > fMaxSignal)
    fMaxSignal = sig;
  fCurvTime.push_back(tCurv);
  fStartTime.push_back(t0);
  fStartTimeErr.push_back(sigma_t);
  fTime10.push_back(t10);
  fTime50.push_back(t50);
  fTime90.push_back(t90);
  fTwin.push_back(groupSize);
  fSaturation.push_back(sat);
  fTOT.push_back(tot);
}


void
CandidateInfo::CalculateLimits(const double zoom)
{
  // limits
  const unsigned int nStations = fPositionX.size();
  if (!nStations)
    return;
  fMinX = fMaxX = fPositionX[0];
  fMinY = fMaxY = fPositionY[0];
  fMinTime = 0; // should always see zero
  fMinStartTime = fMaxStartTime = fMaxTime = fStartTime[0];
  const double sqrtMaxSignal = sqrt(fMaxSignal);
  for (unsigned int i = 0; i < nStations; ++i) {
    const double sqrtS = sqrt(fSignal[i]);
    const double r = kStationMinRadius +
      (kStationMaxRadius - kStationMinRadius)*sqrtS/sqrtMaxSignal;
    fPlotRadius.push_back(r);
    const double x = fPositionX[i];
    const double xdown = x - r;
    if (fMinX > xdown)
      fMinX = xdown;
    const double xup = x + r;
    if (fMaxX < xup)
      fMaxX = xup;
    const double y = fPositionY[i];
    const double ydown = y - r;
    if (fMinY > ydown)
      fMinY = ydown;
    const double yup = y + r;
    if (fMaxY < yup)
      fMaxY = yup;
    const double tCurv = fCurvTime[i];
    if (fMinTime > tCurv)
      fMinTime = tCurv;
    if (fMaxTime < tCurv)
      fMaxTime = tCurv;
    const double t0 = fStartTime[i];
    if (fMinStartTime > t0)
      fMinStartTime = t0;
    if (fMaxStartTime < t0)
      fMaxStartTime = t0;
    if (fMinTime > t0)
      fMinTime = t0;
    if (fMaxTime < t0)
      fMaxTime = t0;
    const double t10 = fTime10[i];
    if (fMinTime > t10)
      fMinTime = t10;
    if (fMaxTime < t10)
      fMaxTime = t10;
    const double t50 = fTime50[i];
    if (fMinTime > t50)
      fMinTime = t50;
    if (fMaxTime < t50)
      fMaxTime = t50;
    /*const double t90 = fTime90[i];
    if (fMinTime > t90)
      fMinTime = t90;
    if (fMaxTime < t90)
      fMaxTime = t90;*/
  }
  // R
  if (fMinR < kLDFMinR)
    fMinR = kLDFMinR;
  if (fMaxR < fMinR)
    fMaxR = 1;
  fMinR /= 2;
  fMaxR *= zoom;
  // ascpect issue
  double dx = zoom * (fMaxX - fMinX);
  if (dx < kMinDX)
    dx = kMinDX;
  double dy = zoom * (fMaxY - fMinY);
  if (dx/kAspect > dy) {
    dy = dx/kAspect;
  } else {
    dx = kAspect*dy;
  }
  dx /= 2;
  dy /= 2;
  const double meanX = (fMinX + fMaxX)/2;
  const double meanY = (fMinY + fMaxY)/2;
  fMinX = meanX - dx;
  fMaxX = meanX + dx;
  fMinY = meanY - dy;
  fMaxY = meanY + dy;
  // time
  double dt = zoom * (fMaxTime - fMinTime) / 2;
  if (!dt)
    dt = 1;
  const double meanTime = (fMinTime + fMaxTime) / 2;
  fMinTime = meanTime - dt;
  fMaxTime = meanTime + dt;
  if (fMinStartTime == fMaxStartTime)
    ++fMaxStartTime;
}


void
CandidateInfo::DrawArray()
{
  const unsigned int nStations = fPositionX.size();
  for (unsigned int i = 0; i < nStations; ++i) {
    const double r = fPlotRadius[i];
    const int sat = fSaturation[i];
    TEllipse* s = nullptr;
    TEllipse* s2 = nullptr;
    TEllipse* thresh = nullptr;
    switch (fTwin[i]) {
    case 0: // normal
      s = new TEllipse(fPositionX[i], fPositionY[i], (sat == 2 ? kLowGainSatFactor : 1)*r);
      if (sat)
        s2 = new TEllipse(fPositionX[i], fPositionY[i], r);
      if (!fTOT[i])
        thresh = new TEllipse(fPositionX[i], fPositionY[i], r/4);
      break;
    case 1: // lower twin
      s = new TEllipse(fPositionX[i], fPositionY[i] - kTwinSplit,
                       (sat == 2 ? kLowGainSatFactor : 1)*r);
      if (sat)
        s2 = new TEllipse(fPositionX[i], fPositionY[i] - kTwinSplit, r);
      if (!fTOT[i])
        thresh = new TEllipse(fPositionX[i], fPositionY[i] - kTwinSplit, r/4,
                              0, 180, 360, 0);
      break;
    default: // upper twin
      s = new TEllipse(fPositionX[i], fPositionY[i] + kTwinSplit,
                       (sat == 2 ? kLowGainSatFactor : 1)*r);
      if (sat)
        s2 = new TEllipse(fPositionX[i], fPositionY[i] + kTwinSplit, r);
      if (!fTOT[i])
        thresh = new TEllipse(fPositionX[i], fPositionY[i] + kTwinSplit,
                              r/4, 0, 0, 180, 0);
      break;
    }
    const int color =
      int(51 + (fStartTime[i] - fMinStartTime)/(fMaxStartTime - fMinStartTime)*49.);
    switch (sat) {
    case 0: // no saturation (no border)
      s->SetLineColor(color);
      if (thresh)
        thresh->SetLineColor(10);
      break;
    case 1: // high
    case 2: // low gain saturation (black border)
      s->SetLineWidth(2);
      s->SetLineColor(color);
      s2->SetLineWidth(2);
      s2->SetLineColor(kBlack);
      if (thresh)
        thresh->SetLineColor(kBlack);
      break;
    }
    s->SetFillColor(color);
    s->SetFillStyle(1001);
    s->Draw();
    fCircles.push_back(s);
    if (s2) {
      s2->Draw();
      fCircles.push_back(s2);
    }
    if (thresh) {
      thresh->SetFillColor(10);
      thresh->SetFillStyle(1001);
      thresh->Draw();
      fCircles.push_back(thresh);
    }
  }
  TLatex label;
  label.SetTextSize(0.03);
  for (unsigned int i = 0; i < nStations; ++i) {
    ostringstream id;
    id << fId[i];
    //twin text split
    const double textSplit = kTwinIdSplit * (fTwin[i] ? (fTwin[i] >= 1 ? -1 : 1) : 0);
    label.DrawLatex(fPositionX[i], fPositionY[i] + textSplit, id.str().c_str());
  }
}


TGraph*
CandidateInfo::GetLDFGraph()
{
  delete fLDFGraph;
  fLDFGraph = nullptr;
  vector<double> r;
  vector<double> signal;
  vector<double> signalErr;
  const unsigned int nR = fR.size();

  for (unsigned int i = 0; i < nR; ++i)
    if (fSaturation[i] != 2) {
      r.push_back(fR[i]);
      signal.push_back(fSignal[i]);
      signalErr.push_back(sqrt(fSignal[i]));
    }
  if (r.size()) {
    fLDFGraph = new TGraphErrors(r.size(), &r[0], &signal[0], 0, &signalErr[0]);
    fLDFGraph->SetMarkerSize(0.5);
    fLDFGraph->SetMarkerStyle(8);
  }
  return fLDFGraph;
}


TGraphBentErrors*
CandidateInfo::GetSaturatedLDFGraph()
{
  delete fSaturatedLDFGraph;
  fSaturatedLDFGraph = nullptr;
  vector<double> r;
  vector<double> signal;
  vector<double> signalErr;
  const unsigned int nR = fR.size();
  for (unsigned int i = 0; i < nR; ++i)
    if (fSaturation[i] == 2) {
      r.push_back(fR[i]);
      signal.push_back(fSignal[i]);
      signalErr.push_back(sqrt(fSignal[i]));
    }
  if (r.size()) {
    fSaturatedLDFGraph = new TGraphBentErrors(r.size(), &r[0], &signal[0], 0, 0, &signalErr[0]);
    fSaturatedLDFGraph->SetMarkerSize(0.5);
    fSaturatedLDFGraph->SetMarkerStyle(22);
    fSaturatedLDFGraph->SetMarkerColor(kBlue);
  }
  return fSaturatedLDFGraph;
}


TGraph*
CandidateInfo::GetStartTimeGraph()
{
  delete fStartTimeGraph;
  fStartTimeGraph = new TGraphErrors(fR.size(), &fR[0], &fStartTime[0], 0, &fStartTimeErr[0]);
  fStartTimeGraph->SetMarkerStyle(22);
  fStartTimeGraph->SetMarkerSize(0.5);
  return fStartTimeGraph;
}


TGraphBentErrors*
CandidateInfo::GetRiseFallTimeGraph()
{
  vector<double> low;
  vector<double> high;
  const unsigned int nR = fR.size();
  for (unsigned int i = 0; i < nR; ++i) {
    const double t50 = fTime50[i];
    low.push_back(t50 - fTime10[i]);
    high.push_back(fTime90[i] - t50);
  }
  delete fRiseFallTimeGraph;
  fRiseFallTimeGraph = new TGraphBentErrors(fR.size(), &fR[0], &fTime50[0], 0, 0, &low[0], &high[0]);
  fRiseFallTimeGraph->SetMarkerStyle(21);
  fRiseFallTimeGraph->SetMarkerSize(0.3);
  fRiseFallTimeGraph->SetMarkerColor(38);
  fRiseFallTimeGraph->SetLineColor(38);
  return fRiseFallTimeGraph;
}


TGraph*
CandidateInfo::GetCurvatureFitGraph()
{
  delete fCurvatureFitGraph;
  fCurvatureFitGraph = new TGraph(fR.size(), &fR[0], &fCurvTime[0]);
  fCurvatureFitGraph->SetMarkerStyle(4);
  fCurvatureFitGraph->SetMarkerSize(0.9);
  fCurvatureFitGraph->SetMarkerColor(kGreen);
  return fCurvatureFitGraph;
}


//=================================================================================================

namespace SdRecPlotterOG {

  class AccidentalInfo {
  public:
    AccidentalInfo() { }
    ~AccidentalInfo();
    void SetLimits(const double minX, const double maxX,
                   const double minY, const double maxY,
                   const double maxR, const double maxSignal);
    void PushStation(const int id, const double x, const double y, const double r,
                     const double sig, const double t, const int groupSize);
    void DrawArray();
    TGraphErrors* GetSignalGraph();
    TGraph* GetTimeGraph();

  private:
    vector<int> fId;
    vector<double> fPositionX;
    vector<double> fPositionY;
    double fMinX = 0;
    double fMaxX = 0;
    double fMinY = 0;
    double fMaxY = 0;
    vector<double> fR;
    double fMaxR = 0;
    vector<double> fSignal;
    double fMaxSignal = 0;
    vector<double> fTime;
    vector<int> fTwin;
    // deletables
    vector<TLine*> fLines;
    vector<TEllipse*> fCircles;
    TGraphErrors* fSignalGraph = nullptr;
    TGraph* fTimeGraph = nullptr;
  };

}


AccidentalInfo::~AccidentalInfo()
{
  for (const auto line : fLines)
    delete line;
  for (const auto circ : fCircles)
    delete circ;
  delete fSignalGraph;
  delete fTimeGraph;
}


void
AccidentalInfo::SetLimits(const double minX, const double maxX,
                          const double minY, const double maxY,
                          const double maxR, const double maxSignal)
{
  const double d = 0.1;
  fMinX = minX + d;
  fMaxX = maxX - d;
  fMinY = minY + d;
  fMaxY = maxY - d;
  fMaxR = maxR;
  fMaxSignal = maxSignal;
}


void
AccidentalInfo::PushStation(const int id, const double x, const double y,
                            const double r, const double sig, const double t,
                            const int groupSize)
{
  fId.push_back(id);
  fPositionX.push_back(x);
  fPositionY.push_back(y);
  fR.push_back(r);
  fSignal.push_back(sig);
  fTime.push_back(t);
  fTwin.push_back(groupSize);
}


void
AccidentalInfo::DrawArray()
{
  const double sqrtMaxSignal = sqrt(fMaxSignal);
  const unsigned int nStations = fId.size();
  vector<bool> isDrawn(nStations, false);
  for (unsigned int i = 0; i < nStations; ++i) {
    const double x = fPositionX[i];
    const double y = fPositionY[i];
    const double sqrtS = sqrt(fSignal[i]);
    const double r = CandidateInfo::kStationMinRadius +
      (CandidateInfo::kStationMaxRadius - CandidateInfo::kStationMinRadius)*sqrtS/sqrtMaxSignal;
    const double size = r/sqrt(2.);
    if (fMinX + size < x && x < fMaxX - size &&
        fMinY + size < y && y < fMaxY - size) {
      TLine* line1 = nullptr;
      TLine* line2 = nullptr;
      TEllipse* c = nullptr;
      switch (fTwin[i]) {
      case 0:  // normal
        line1 = new TLine(x - size, y - size, x + size, y + size);
        line2 = new TLine(x - size, y + size, x + size, y - size);
        c = new TEllipse(x, y, r);
        break;
      case 1:  // lower twin
        {
          const double yy = y - CandidateInfo::kTwinSplit;
          line1 = new TLine(x - size, yy - size, x, yy);
          line2 = new TLine(x, yy, x + size, yy - size);
          c = new TEllipse(x, yy, r, 0, 180, 360, 0);
          break;
        }
      default:  // upper twin
        {
          const double yy = y + CandidateInfo::kTwinSplit;
          line1 = new TLine(x - size, yy + size, x, yy);
          line2 = new TLine(x, yy, x + size, yy + size);
          c = new TEllipse(x, yy, r, 0, 0, 180, 0);
          break;
        }
      }
      line1->SetLineColor(kRed);
      line2->SetLineColor(kRed);
      c->SetLineColor(kRed);
      line1->Draw();
      line2->Draw();
      c->Draw();
      fLines.push_back(line1);
      fLines.push_back(line2);
      fCircles.push_back(c);
      isDrawn[i] = true;
    }
  }
  TLatex label;
  label.SetTextSize(0.03);
  for (unsigned int i = 0; i < nStations; ++i)
    if (isDrawn[i]) {
      ostringstream id;
      id << fId[i];
      // twin text split
      const double textSplit =
        CandidateInfo::kTwinIdSplit * (fTwin[i] ? (fTwin[i] >= 1 ? -1 : 1) : 0);
      label.DrawLatex(fPositionX[i], fPositionY[i] + textSplit, id.str().c_str());
    }
}


TGraphErrors*
AccidentalInfo::GetSignalGraph()
{
  delete fSignalGraph;
  fSignalGraph = nullptr;
  vector<double> r;
  vector<double> signal;
  vector<double> signalErr;
  const unsigned int nR = fR.size();
  for (unsigned int i = 0; i < nR; ++i)
    if (fR[i] < fMaxR) {
      r.push_back(fR[i]);
      signal.push_back(fSignal[i]);
      signalErr.push_back(sqrt(fSignal[i]));
    }
  if (r.size()) {
    fSignalGraph = new TGraphErrors(r.size(), &r[0], &signal[0], 0, &signalErr[0]);
    fSignalGraph->SetLineColor(kRed);
    fSignalGraph->SetMarkerColor(kRed);
    fSignalGraph->SetMarkerStyle(8);
    fSignalGraph->SetMarkerSize(0.5);
  }
  return fSignalGraph;
}


TGraph*
AccidentalInfo::GetTimeGraph()
{
  delete fTimeGraph;
  fTimeGraph = nullptr;
  const unsigned int nR = fR.size();
  if (nR) {
    fTimeGraph = new TGraph(nR, &fR[0], &fTime[0]);
    fTimeGraph->SetMarkerStyle(22);
    fTimeGraph->SetMarkerSize(0.5);
    fTimeGraph->SetMarkerColor(kRed);
  }
  return fTimeGraph;
}


//=================================================================================================

namespace SdRecPlotterOG {

  class SilentInfo {
  public:
    SilentInfo() { }
    ~SilentInfo();
    void SetLimits(const double minX, const double maxX,
                   const double minY, const double maxY,
                   const double maxR);
    void PushStation(const double x, const double y,
                     const double r, const int groupSize);
    void DrawArray();
    TGraph* GetLimitsGraph();

    static double kSilentThreshold;

  private:
    vector<double> fPositionX;
    vector<double> fPositionY;
    double fMinX = 0;
    double fMaxX = 0;
    double fMinY = 0;
    double fMaxY = 0;
    vector<double> fR;
    double fMaxR = 0;
    vector<int> fTwin;
    // deletables
    vector<TEllipse*> fCircles;
    TGraph* fLimitsGraph = nullptr;
  };


  double SilentInfo::kSilentThreshold = 2.8;

}


SilentInfo::~SilentInfo()
{
  for (const auto circ : fCircles)
    delete circ;
  delete fLimitsGraph;
}


void
SilentInfo::SetLimits(const double minX, const double maxX,
                      const double minY, const double maxY,
                      const double maxR)
{
  const double d = 0.1;
  fMinX = minX + d;
  fMaxX = maxX - d;
  fMinY = minY + d;
  fMaxY = maxY - d;
  fMaxR = maxR;
}


void
SilentInfo::PushStation(const double x, const double y, const double r, const int groupSize)
{
  fPositionX.push_back(x);
  fPositionY.push_back(y);
  fR.push_back(r);
  fTwin.push_back(groupSize);
}


void
SilentInfo::DrawArray()
{
  const double size = 0.05; //km
  const double twinSplit = 0.03; //km
  const unsigned int nStations = fPositionX.size();
  for (unsigned int i = 0; i < nStations; ++i) {
    const double x = fPositionX[i];
    const double y = fPositionY[i];
    if (fMinX < x && x < fMaxX && fMinY < y && y < fMaxY) {
      TEllipse* c = nullptr;
      switch (fTwin[i]) {
      case 0: //normal
        c = new TEllipse(x, y, size);
        break;
      case 1: //lower twin
        c = new TEllipse(x, y - twinSplit, size, 0, 180, 360, 0);
        break;
      default: //upper twin
        c = new TEllipse(x, y + twinSplit, size, 0, 0, 180, 0);
        break;
      }
      c->Draw();
      fCircles.push_back(c);
    }
  }
}


TGraph*
SilentInfo::GetLimitsGraph()
{
  delete fLimitsGraph;
  fLimitsGraph = nullptr;
  vector<double> r;
  const unsigned int nR = fR.size();
  for (unsigned int i = 0; i < nR; ++i)
    if (fR[i] < fMaxR)
      r.push_back(fR[i]);
  const unsigned nLimits = r.size();
  if (nLimits) {
    vector<double> trigger(nLimits, kSilentThreshold);
    fLimitsGraph = new TGraph(nLimits, &r[0], &trigger[0]);
    fLimitsGraph->SetMarkerStyle(23);
    fLimitsGraph->SetMarkerColor(kBlue);
    fLimitsGraph->SetMarkerSize(0.5);
  }
  return fLimitsGraph;
}


//=================================================================================================

void
SdRecPlotter::DrawEvent(const Event& event)
{
  const SEvent& sEvent = event.GetSEvent();
  const sevt::Header& sHeader = sEvent.GetHeader();
  //const int eventId = sHeader.GetId();
  const string globalId = event.GetHeader().GetId();
  const det::Detector& detector = det::Detector::GetInstance();
  const SDetector& sDetector = detector.GetSDetector();
  const STimeVariance& timeVariance = STimeVariance::GetInstance();
  const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();
  const ShowerSRecData& shRec = event.GetRecShower().GetSRecShower();
  const Point& core = shRec.GetCorePosition();
  const TimeStamp& coreTime = shRec.GetCoreTime();
  const Vector& axis = shRec.GetAxis();
  const double curvature = shRec.GetCurvature();
  const double curvR = curvature ? 1/curvature : 1000*kilometer;
  const Vector showerCenter(curvR*axis);
  const TimeStamp curvT0(coreTime - TimeInterval(curvR/kSpeedOfLight));

  CandidateInfo candidates;
  AccidentalInfo accidentals;
  SilentInfo silents;

  // curvature data
  double minHeight = 0;
  double maxHeight = 0;

  for (SEvent::ConstStationIterator sIt = sEvent.StationsBegin();
       sIt != sEvent.StationsEnd(); ++sIt) {

    const int sId = sIt->GetId();
    const sdet::Station& dStation = sDetector.GetStation(sId);
    const Point& sPosition = dStation.GetPosition();
    const double x = sPosition.GetX(siteCS)/kilometer;
    const double y = sPosition.GetY(siteCS)/kilometer;
    const Vector sRelPos(sPosition - core);
    const double r = RPerp(axis, sRelPos)/kilometer;
    const int groupSize = dStation.GetNumberOfPartners();

    if (sIt->IsCandidate()) {

      const double radius = (showerCenter - sRelPos).GetMag();
      const double height = axis*sRelPos;
      if (minHeight > height)
        minHeight = height;
      if (maxHeight < height)
        maxHeight = height;
      const StationRecData& sRec = sIt->GetRecData();
      const double signal = sRec.GetTotalSignal();
      const TimeStamp planeFrontTime(coreTime -
                                     TimeInterval(height/kSpeedOfLight));
      const double tCurv = (curvT0 + TimeInterval(radius/kSpeedOfLight) - planeFrontTime).GetInterval();
      const double t0 = (sRec.GetSignalStartTime() - planeFrontTime).GetInterval();
      const double traceTime = (sIt->GetTraceStartTime() - planeFrontTime).GetInterval();
      const int start = sRec.GetSignalStartSlot();
      const int end = sRec.GetSignalEndSlot();
      const TraceD& trace = sIt->GetVEMTrace();
      const double t10 = traceTime + TraceAlgorithm::TimeAtRelativeSignalX(trace, start, end, 10);
      const double t50 = traceTime + TraceAlgorithm::TimeAtRelativeSignalX(trace, start, end, 50);
      const double t90 = traceTime + TraceAlgorithm::TimeAtRelativeSignalX(trace, start, end, 90);
      const double sigma_t = sqrt(timeVariance.GetTimeSigma2(signal, t50 - t0, axis.GetCosTheta(siteCS), r));
      const int sat = sIt->IsLowGainSaturation() ? 2 : sIt->IsHighGainSaturation();
      const bool tot = sIt->GetTriggerData().GetAlgorithm() == StationTriggerData::eTimeOverThreshold;
      candidates.PushStation(sId, x, y, r, signal, curvature ? tCurv/nanosecond : 0,
                             t0/nanosecond, sigma_t/nanosecond, t10/nanosecond, t50/nanosecond, t90/nanosecond,
                             groupSize, sat, tot);

    } else if (sIt->IsRejected() && sIt->HasRecData()) {

      const StationRecData& sRec = sIt->GetRecData();
      const TimeStamp planeFrontTime(coreTime -
                                     TimeInterval((sPosition - core)*axis/kSpeedOfLight));
      const double t0 = (sRec.GetSignalStartTime() - planeFrontTime).GetInterval();
      accidentals.PushStation(sId, x, y, r, sIt->GetRecData().GetTotalSignal(), t0/nanosecond, groupSize);

    } else // Silent

      silents.PushStation(x, y, r, groupSize);

  }

  candidates.CalculateLimits(1.2);
  const double minX = candidates.GetMinX();
  const double maxX = candidates.GetMaxX();
  const double minY = candidates.GetMinY();
  const double maxY = candidates.GetMaxY();
  const double minR = candidates.GetMinR();
  const double maxR = candidates.GetMaxR();

  ostringstream cName;
  //  cName << "Event_" << eventId;
  cName << "Event_" << globalId;
  TCanvas canvas(cName.str().c_str(), cName.str().c_str(), 1);
  canvas.Divide(2, 2);

  // array view
  canvas.cd(1);
  TH2F array("array", "Array layout", 10, minX, maxX, 10, minY, maxY);
  array.GetXaxis()->SetTitle("x [km]");
  array.GetYaxis()->SetTitle("y [km]");
  array.SetStats(kFALSE);
  array.Draw();
  // accidentals
  accidentals.SetLimits(minX, maxX, minY, maxY, maxR, candidates.GetMaxSignal());
  accidentals.DrawArray();
  // silent
  silents.SetLimits(minX, maxX, minY, maxY, maxR);
  silents.DrawArray();
  // candidates
  candidates.DrawArray();
  // seed
  const vector<int>& seedStation = shRec.GetReconstructionSeed();
  if (seedStation.size()) {
    vector<double> sx;
    vector<double> sy;
    for (unsigned int i = 0; i < seedStation.size(); ++i) {
      const Point& sPos = sDetector.GetStation(seedStation[i]).GetPosition();
      sx.push_back(sPos.GetX(siteCS)/kilometer);
      sy.push_back(sPos.GetY(siteCS)/kilometer);
    }
    const Point& sPos = sDetector.GetStation(seedStation[0]).GetPosition();
    sx.push_back(sPos.GetX(siteCS)/kilometer);
    sy.push_back(sPos.GetY(siteCS)/kilometer);
    TPolyLine* const grSeed = new TPolyLine(sx.size(), &sx[0], &sy[0]);
    grSeed->SetBit(kCanDelete);
    grSeed->SetLineColor(kGreen);
    grSeed->Draw();
  }
  // coreMC
  const bool hasSimShower = event.HasSimShower();
  if (hasSimShower) {
    const ShowerSimData& simShower = event.GetSimShower();
    const Point& coreMC = simShower.GetPosition();
    const Vector& axisMC = -simShower.GetDirection();
    const double phiMC = axisMC.GetPhi(siteCS);
    const double coreMCX = coreMC.GetX(siteCS)/kilometer;
    const double coreMCY = coreMC.GetY(siteCS)/kilometer;
    TEllipse* const coreMCMark = new TEllipse(coreMCX, coreMCY, 0.25, 0, 31, 329, phiMC/degree);
    coreMCMark->SetBit(kCanDelete);
    coreMCMark->SetLineColor(kBlue);
    coreMCMark->SetLineWidth(2);
    coreMCMark->SetFillStyle(1001);
    coreMCMark->SetFillColor(10);
    coreMCMark->Draw();
    const double coreLength = 2;
    TLine* const coreMCLine =
      new TLine(coreMCX, coreMCY, coreMCX + coreLength*cos(phiMC), coreMCY + coreLength*sin(phiMC));
    coreMCLine->SetBit(kCanDelete);
    coreMCLine->SetLineColor(kBlue);
    coreMCLine->SetLineWidth(2);
    coreMCLine->Draw();
  }
  // core
  const double phi = axis.GetPhi(siteCS);
  const double coreX = core.GetX(siteCS)/kilometer;
  const double coreY = core.GetY(siteCS)/kilometer;
  TEllipse coreMark(coreX, coreY, 0.25, 0, 31, 329, phi/degree);
  coreMark.SetLineColor(kRed);
  coreMark.SetLineWidth(2);
  coreMark.SetFillStyle(1001);
  coreMark.SetFillColor(10);
  coreMark.Draw();
  const double coreLength = 2;
  TLine coreLine(coreX, coreY, coreX + coreLength*cos(phi), coreY + coreLength*sin(phi));
  coreLine.SetLineColor(kRed);
  coreLine.SetLineWidth(2);
  coreLine.Draw();

  // timing
  canvas.cd(2);
  double timeZeroX[] = { 0, maxR };
  double timeZeroY[] = { 0, 0 };
  TGraph* const grTimeZero = new TGraph(2, timeZeroX, timeZeroY);
  grTimeZero->SetBit(kCanDelete);
  grTimeZero->SetLineStyle(3);
  grTimeZero->SetTitle("Timing");
  grTimeZero->GetXaxis()->SetTitle("r [km]");
  grTimeZero->GetYaxis()->SetTitle("time to plane front [ns]");
  grTimeZero->GetXaxis()->SetLimits(0, maxR);
  grTimeZero->SetMinimum(candidates.GetMinTime());
  grTimeZero->SetMaximum(candidates.GetMaxTime());
  grTimeZero->Draw("al");
  TGraph* grStartTime = candidates.GetStartTimeGraph();
  if (curvature) {
    vector<double> timingR;
    vector<double> timingCurve;
    const double rc = 1/curvature;
    const int nPoints = 100;
    const double dr = maxR*kilometer/(nPoints - 1);
    const double rc1 = rc - maxHeight;
    for (int i = nPoints - 1; i >= 0; --i) {
      const double r = i*dr;
      timingR.push_back(r/kilometer);
      timingCurve.push_back((sqrt(r*r+rc1*rc1) - rc1)/kSpeedOfLight);
    }
    const double rc2 = rc - minHeight;
    for (int i = 0; i < nPoints; ++i) {
      const double r = i*dr;
      timingR.push_back(r/kilometer);
      timingCurve.push_back((sqrt(r*r+rc2*rc2) - rc2)/kSpeedOfLight);
    }
    TGraph* const grCurvature = new TGraph(timingR.size(), &timingR[0], &timingCurve[0]);
    grCurvature->SetBit(kCanDelete);
    grCurvature->SetLineColor(7);
    grCurvature->SetLineStyle(2);
    grCurvature->Draw("samel");
    TGraph* const grCurvatureFit = candidates.GetCurvatureFitGraph();
    grCurvatureFit->Draw("samep");
  }
  candidates.GetRiseFallTimeGraph()->Draw("samep");
  grStartTime->Draw("samep");
  TGraph* grAccidentalTime = accidentals.GetTimeGraph();
  if (grAccidentalTime)
    grAccidentalTime->Draw("samep");

  // LDF
  canvas.cd(3);
  gPad->SetLogy();
  // LDF signal
  TGraph* grSignal = candidates.GetLDFGraph();
  TGraphBentErrors* grSatSignal = candidates.GetSaturatedLDFGraph();

  if (shRec.HasLDF()) {
    const TabulatedFunctionErrors& ldf = shRec.GetLDF();
    // LDF error
    vector<double> signalR;
    vector<double> signal;
    vector<double> signalErr;
    const unsigned int nLDFPoints = ldf.GetNPoints();
    for (unsigned int i = 0; i < nLDFPoints; ++i) {
      signalR.push_back(ldf.GetX(i)/kilometer);
      const double sig = ldf.GetY(i);
      signal.push_back(sig);
      signalErr.push_back(sig + ldf.GetYErr(i));
    }
    for (int i = nLDFPoints - 1; i >= 0; --i) {
      signalR.push_back(ldf.GetX(i)/kilometer);
      const double low = ldf.GetY(i) - ldf.GetYErr(i);
      signalErr.push_back(low > 0.01 ? low : 0.01);
    }

    TGraph* const grLDFError = new TGraph(signalR.size(), &signalR[0], &signalErr[0]);
    grLDFError->SetBit(kCanDelete);
    grLDFError->SetTitle("LDF fit");
    grLDFError->GetYaxis()->SetTitle("Signal [VEM]");
    grLDFError->GetXaxis()->SetTitle("r [km]");
    grLDFError->SetFillStyle(1001);
    grLDFError->SetFillColor(18);
    grLDFError->GetXaxis()->SetLimits(minR, maxR);
    grLDFError->SetMinimum(min(SilentInfo::kSilentThreshold/2,
                               grSatSignal ? min(grSignal->GetYaxis()->GetXmin(), grSatSignal->GetYaxis()->GetXmin()) :
                               grSignal->GetYaxis()->GetXmin()));
    grLDFError->SetMaximum(1.2 * (grSatSignal ? max(grSignal->GetYaxis()->GetXmax(), grSatSignal->GetYaxis()->GetXmax()) :
                                  grSignal->GetYaxis()->GetXmax()));
    grLDFError->Draw("af");
    // LDF fit
    TGraph* grLDFFit = new TGraph(signal.size(), &signalR[0], &signal[0]);
    grLDFFit->SetBit(kCanDelete);
    grLDFFit->SetLineColor(kGreen);
    grLDFFit->SetLineWidth(2);
    grLDFFit->Draw("samel");
    double rOpt = shRec.GetLDFOptimumDistance();
    if (rOpt && ldf.GetX(0) <= rOpt && rOpt <= ldf.GetX(ldf.GetNPoints()-1)) {
      const double sROpt = ldf.Y(rOpt);
      rOpt /= kilometer;
      TGraph* const grLDFROpt = new TGraph(1, &rOpt, &sROpt);
      grLDFROpt->SetBit(kCanDelete);
      grLDFROpt->SetMarkerColor(kGreen);
      grLDFROpt->SetMarkerStyle(kStar);
      grLDFROpt->SetMarkerSize(1.3);
      grLDFROpt->Draw("samep");
    }
  } else {
    WARNING("No LDF tabulated function object found.");
    grSignal->SetTitle("LDF fit");
    grSignal->GetYaxis()->SetTitle("Signal [VEM]");
    grSignal->GetXaxis()->SetTitle("r [km]");
    grSignal->GetXaxis()->SetLimits(minR, maxR);
    grSignal->SetMinimum(min(SilentInfo::kSilentThreshold/2,
                             grSatSignal ? min(grSignal->GetYaxis()->GetXmin(), grSatSignal->GetYaxis()->GetXmin()) :
                             grSignal->GetYaxis()->GetXmin()));
    grSignal->SetMaximum(1.2 * (grSatSignal ? max(grSignal->GetYaxis()->GetXmax(), grSatSignal->GetYaxis()->GetXmax()) :
                                grSignal->GetYaxis()->GetXmax()));
    grSignal->Draw("ap");
  }
  // silent limits
  TGraph* grSilentLimits = silents.GetLimitsGraph();
  if (grSilentLimits)
    grSilentLimits->Draw("samep");
  // LDF accidental
  TGraphErrors* grAccidentalSignal = accidentals.GetSignalGraph();
  if (grAccidentalSignal)
    grAccidentalSignal->Draw("samep");
  // LDF candiates
  if (shRec.HasLDF())
    grSignal->Draw("samep");
  if (grSatSignal)
    grSatSignal->Draw("samep");
  gPad->RedrawAxis();

  // Data
  canvas.cd(4);
  gPad->Range(0,0,1,1);
  TPaveText pt(0.1,0.1,0.9,0.9, "br");
  pt.SetFillColor(kWhite);
  pt.SetTextAlign(11);
  ostringstream info;
  //info << "Event:  " << eventId;
  info << " \n"
          "Event:  " << globalId << '\n';
  Point coreMC;
  double eastingMC = 0;
  double northingMC = 0;
  double thetaMC = 0;
  double phiMC = 0;
  double energyMC = 0;
  if (hasSimShower) {
    const ShowerSimData& simShower = event.GetSimShower();
    coreMC = simShower.GetPosition();
    const Vector axisMC = -simShower.GetDirection();
    const UTMPoint utmMCCore(coreMC, ReferenceEllipsoid::GetWGS84());
    phiMC = axisMC.GetPhi(siteCS);
    thetaMC = axisMC.GetTheta(siteCS);
    eastingMC = utmMCCore.GetEasting();
    northingMC = utmMCCore.GetNorthing();
    energyMC = simShower.GetEnergy();
  }

  const TimeStamp& eventTime = sHeader.GetTime();
  info << "Time:  " << UTCDateTime(eventTime) << "\n"
          "GPS Time:  " << eventTime.GetGPSSecond() << " s, "
       << eventTime.GetGPSNanoSecond() << " ns\n"
          "T4:  " << GetT4TriggerName(shRec.GetT4Trigger())
       << "   T5:  " << GetT5TriggerName(shRec.GetT5Trigger());
  const int badId = shRec.GetBadPeriodId();
  if (badId)
    info << "   BadPeriod:  " << badId;
  info << "\nReconstruction stage:  " << shRec.GetLDFRecStage() << '\n';
  const UTMPoint utmCore(core, ReferenceEllipsoid::GetWGS84());
  const Vector& coreError = shRec.GetCoreError();
  info << "Easting:  " << setprecision(8) << utmCore.GetEasting()/meter << " #pm "
       << setprecision(3) << coreError.GetX(siteCS)/meter;
  if (hasSimShower)
    info << " (MC: " << setprecision(8) << eastingMC/meter << ')';
  info << " [m]\n"
          "Northing:  " << setprecision(9) << utmCore.GetNorthing()/meter << " #pm "
       << setprecision(3) << coreError.GetY(siteCS)/meter;
  if (hasSimShower)
    info << " (MC: " << setprecision(8) << northingMC/meter << ')';
  info << " [m]\n";
  if (hasSimShower) {
    info << "Distance to MC core:  " << setprecision(5) << (core - coreMC).GetMag()/meter
         << " [m]\n";
  }
  info << "#theta:  " << setprecision(4) << axis.GetTheta(siteCS)/degree << " #pm "
       << setprecision(2) << shRec.GetThetaError()/degree;
  if (hasSimShower)
    info << " (MC: " << setprecision(4) << thetaMC/degree << ')';
  info << " [#circ]      "
          "#phi:  " << setprecision(5) << axis.GetPhi(siteCS)/degree << " #pm "
       << setprecision(2) << shRec.GetPhiError()/degree;
  if (hasSimShower)
    info << " (MC: " << setprecision(5) << phiMC/degree << ')';
  info << " [#circ]\n";
  if (curvature) {
    const double rc = 1/curvature;
    info << "R_{c}:  " << setprecision(4) << rc/kilometer << " #pm "
         << setprecision(2) << shRec.GetCurvatureError()*rc*rc/kilometer << " [km]\n";
  }
  info << shRec.GetShowerSizeLabel() << ":  " << setprecision(4) << shRec.GetShowerSize()
       << " #pm " << setprecision(2) << shRec.GetShowerSizeError()
       << " (#pm" << shRec.GetShowerSizeSystematics() << " sys)";
  const double rOpt = shRec.GetLDFOptimumDistance();
  if (rOpt)
    info << "     r_{opt}:  " << setprecision(5) << rOpt << " [m]";
  info << "\n#beta:  " << setprecision(4) << shRec.GetBeta()
       << " #pm " << setprecision(2) << shRec.GetBetaError()
       << " (#pm" << shRec.GetBetaSystematics() << " sys)\n";
  info << "Energy:  " << setprecision(4) << shRec.GetEnergy()/EeV << " #pm "
       << setprecision(2) << shRec.GetEnergyError()/EeV;
  if (hasSimShower)
    info << " (MC: " << setprecision(4) << energyMC/EeV << ')';
  info << " [EeV]\n ";
  AddText(pt, info.str());
  pt.Draw();

  canvas.Update();
  if (fGeneratePDFFile) {
    static string multipage = "(";
    canvas.Print((fPDFFileName + multipage).c_str(), "pdf");
    multipage = "";
  }
  if (fROOTFile) {
    SaveCurrentTDirectory save;
    fROOTFile->cd();
    canvas.Write();
  }
}