RdSdCheck/Validation.cc

Go to the documentation of this file.

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

#include <sevt/SEvent.h>
#include <sevt/SortCriteria.h>

#include <utl/ErrorLogger.h>
#include <utl/Branch.h>
#include <utl/AugerUnits.h>

#include <tst/Validatrix.h>

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>

#include "Validation.h"

using namespace std;
using namespace evt;
using namespace sevt;
using namespace fevt;
using namespace fwk;
using namespace utl;
using namespace tst::Validatrix;
using namespace ValidationNS;


VModule::ResultFlag
Validation::Init()
{
  const Branch topB = CentralConfig::GetInstance()->GetTopBranch("Validation");
  topB.GetChild("RefFilename").GetData(fRefFileName);

  fResults = new ofstream(fRefFileName + ".new");

  return eSuccess;
}


VModule::ResultFlag
Validation::Run(evt::Event& event)
{
  if (!event.HasSEvent()) {
    ERROR("Event should contain an SEvent.");
    return eFailure;
  }

  event.GetSEvent().SortStations(ByIncreasingId());
  const SEvent& sEvent = event.GetSEvent();

  const string eventId = event.GetHeader().GetId();

  *fResults << BeLabel("SD_rec");

  const double sigEps = 0.05;

  for (const auto& station : sEvent.StationsRange()) {
    if (station.HasRecData()) {
      ostringstream label;
      label << "E" << eventId << "_S" << station.GetId();
      *fResults << BeLabel(label.str());
      *fResults << BeCloseRel(station.GetRecData().GetTotalSignal(), sigEps);
    }
  }

  *fResults << BeLabel("Shower_rec");

  // Shower-level rec data
  if (!event.HasRecShower())
    ERROR("Event does not contain a RecShower");

  const ShowerRecData& showerRec = event.GetRecShower();

  // SD rec (shower level)

  if (!showerRec.HasSRecShower())
    ERROR("RecShower does not have SRecShower");

  *fResults << BeLabel("Shower_S_rec");

  const ShowerSRecData& showerSRec = showerRec.GetSRecShower();

  const Point& corePos = showerSRec.GetCorePosition();
  const Point& baryPos = showerSRec.GetBarycenter();

  const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(corePos);

  const CoordinateSystemPtr baryCS = LocalCoordinateSystem::Create(baryPos);

  const double relEps = 1e-3;

  ostringstream label;

  label << eventId << "_ShowerSize";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetShowerSize(), relEps);
  label.str("");

  label << eventId << "_ShowerSizeError";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetShowerSizeError(), relEps);
  label.str("");

  label << eventId << "_LDFChi2";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetLDFChi2(), relEps);
  label.str("");

  label << eventId << "_LDFNdof";
  *fResults << BeLabel(label.str());
  *fResults << BeEqual(showerSRec.GetLDFNdof());
  label.str("");

  label << eventId << "_Theta";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetAxis().GetTheta(coreCS), relEps);
  label.str("");

  label << eventId << "_Phi";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetAxis().GetPhi(coreCS), relEps);
  label.str("");

  label << eventId << "_ThetaError";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetThetaError(), relEps);
  label.str("");

  label << eventId << "_PhiError";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetPhiError(), relEps);
  label.str("");

  const double posEps = 0.1*meter;

  label << eventId << "_CoreX";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseAbs(corePos.GetX(baryCS), posEps);
  label.str("");

  label << eventId << "_CoreY";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseAbs(corePos.GetY(baryCS), posEps);
  label.str("");

  label << eventId << "_CoreZ";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseAbs(corePos.GetZ(baryCS), posEps);
  label.str("");

  /* comment out because for events that only have a barycenter reconstruction, comparison of 0 with 0 fails
  label << eventId << "_CorrelationXY";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetCorrelationXY(), relEps);
  label.str("");*/

  label << eventId << "_CoreXVariance";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetParameterCovariance(sevt::eShowerAxisX, sevt::eShowerAxisX), relEps);
  label.str("");

  label << eventId << "_CoreYVariance";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetParameterCovariance(sevt::eShowerAxisY, sevt::eShowerAxisY), relEps);
  label.str("");

  label << eventId << "_CoreXYCovariance";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetParameterCovariance(sevt::eShowerAxisX, sevt::eShowerAxisY), relEps);
  label.str("");

  label << eventId << "_CoreTSecond";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetCoreTime().GetGPSSecond(), relEps);
  label.str("");

  label << eventId << "_CoreTNanoSecond";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetCoreTime().GetGPSNanoSecond(), relEps);
  label.str("");

  label << eventId << "_CoreXError";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseAbs(showerSRec.GetCoreError().GetX(baryCS), posEps);
  label.str("");

  label << eventId << "_CoreYError";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseAbs(showerSRec.GetCoreError().GetY(baryCS), posEps);
  label.str("");

  label << eventId << "_CoreZError";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseAbs(showerSRec.GetCoreError().GetZ(baryCS), posEps);
  label.str("");

  label << eventId << "_CoreTError";
  *fResults << BeLabel(label.str());
  *fResults << BeCloseRel(showerSRec.GetCoreTimeError().GetInterval(), relEps);
  label.str("");

  return eSuccess;
}


VModule::ResultFlag
Validation::Finish()
{
  fResults->close();

  ifstream newFile(fRefFileName + ".new");
  ifstream oldFile(fRefFileName);

  const bool comp = Compare(oldFile, newFile, /*failOnFirst=*/true);

  if (comp) {
    INFO("\n"
         "#--------------------------------#\n"
         ":all SD tests were successful:\n"
         "#--------------------------------#");
  } else {
    INFO("\n"
         "#------------------------------#\n"
         ":one or more SD tests failed:\n"
         "#------------------------------#");
  }

  return comp ? eSuccess : eFailure;
}