HdSimWriteRec/Validation.cc

Go to the documentation of this file.

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

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

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/EyeRecData.h>
#include <fevt/TelescopeRecData.h>

#include <evt/ShowerFRecData.h>

#include <utl/ErrorLogger.h>
#include <utl/Branch.h>
#include <utl/Particle.h>
#include <utl/AugerUnits.h>
#include <utl/Point.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <tst/Validatrix.h>

#include <fwk/CentralConfig.h>

#include <iostream>

#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()
{
  INFO("Doing the validation !!!!");
  Branch topB = CentralConfig::GetInstance()->GetTopBranch("Validation");
  topB.GetChild("RefFilename").GetData(fRefFileName);
  string mode;
  topB.GetChild("Mode").GetData(mode);

  cout << "Validation mode = " << mode << endl;

  fMode = mode == "Create" ? Validation::eCreate : Validation::eCompare;

  fResults.str("");

  return eSuccess;
}


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

  SEvent& sEvent = event.GetSEvent();

  fResults << BeLabel("SD_rec");

  // NB you have to sort the stations before writing to ref
  // file, since there is no guarantee the
  // station list will come back in the same order after serialization

  sEvent.SortStations(ByIncreasingId());

  ostringstream label;
  for (SEvent::ConstStationIterator sIt = sEvent.StationsBegin();
       sIt != sEvent.StationsEnd(); ++sIt) {
    if (sIt->HasRecData()) {
      label.str("");
      label << "Station_" << sIt->GetId();
      fResults << BeLabel(label.str());
      fResults << BeCloseRel(sIt->GetRecData().GetTotalSignal(), /*tolerance*/0.05);
    }
  }

  // FD
  if (!event.HasFEvent()) {
    ERROR("Event should contain an FEvent.");
    return eFailure;
  }

  const FEvent& fEvent = event.GetFEvent();

  fResults << BeLabel("FD_rec");

  for (FEvent::ConstEyeIterator eIt = fEvent.EyesBegin(ComponentSelector::eHasData);
       eIt != fEvent.EyesEnd(ComponentSelector::eHasData); ++eIt) {

    if (eIt->HasRecData()) {
      label.str("");
      label << "Eye_" << eIt->GetId();
      fResults << BeLabel(label.str());

      const ShowerFRecData& showerRec = eIt->GetRecData().GetFRecShower();
      fResults << BeLabel("Stations_in_fit");
      fResults << BeEqual(showerRec.GetStationIds());

      fResults << BeLabel("TotalEnergy(EeV)");
      fResults << BeCloseRel(showerRec.GetTotalEnergy()/EeV, /*tolerance*/1e-5);

      fResults << BeLabel("TotalEnergyError(EeV)");
      fResults << BeCloseRel(showerRec.GetTotalEnergyError()/EeV, /*tolerance*/1e-1);

      fResults << BeLabel("Xmax(g/cm2)");
      fResults << BeCloseRel(showerRec.GetGHParameters().GetXMax()/gram*cm2, 1e-5);

      fResults << BeLabel("Xmax2Error");
      fResults << BeCloseRel(showerRec.GetGHParameters().GetXMaxError()/gram*cm2, 1e-1);
    }

    for (Eye::ConstTelescopeIterator tIt = eIt->TelescopesBegin(ComponentSelector::eHasData);
         tIt != eIt->TelescopesEnd(ComponentSelector::eHasData); ++tIt) {

      if (tIt->HasRecData()) {
        label.str("");
        label << "Telescope_" << tIt->GetId();
        fResults << BeLabel(label.str());

        const TelescopeRecData& tRec = tIt->GetRecData();
        //fResults << BeLabel("zeta");
        //fResults << BeCloseRel(tRec.GetZeta(), /*tolerance*/0.05);
        fResults << BeLabel("photon_trace_start_sec");
        fResults << BeEqual(tRec.GetPhotonsStartTime().GetGPSSecond());
        fResults << BeLabel("photon_trace_start_ns");
        fResults << BeCloseAbs(tRec.GetPhotonsStartTime().GetGPSNanoSecond(), 50*ns);
      }
    }
  }

  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();

  fResults << BeLabel("S1000");
  fResults << BeCloseRel(showerSRec.GetShowerSize(), 1e-5);

  // Confirm that the MC truth is stored and read back

  if (!event.HasSimShower())
    ERROR("Event does not have ShowerSimData");

  const ShowerSimData& sSim = event.GetSimShower();

  const Point& cP = sSim.GetPosition();
  // get reference coordinate system of detector (usually PampaAmarilla)
  const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

  fResults << BeLabel("core position (pampa CS)");
  fResults << BeCloseRel(cP.GetX(referenceCS)/km) << " "
           << BeCloseRel(cP.GetY(referenceCS)/km) << " "
           << BeCloseRel(cP.GetZ(referenceCS)/km);

  return eSuccess;
}


VModule::ResultFlag
Validation::Finish()
{
  bool comp = true;

  if (fMode == Validation::eCreate) {
    INFO("Writing reference file");
    ofstream fout(fRefFileName.c_str());
    fout << fResults.str();
    fout.close();
  } else {
    INFO("Comparing re-reconstruction results to reference file");
    ifstream savedFile(fRefFileName.c_str());
    comp = Compare(fResults, savedFile, /*failOnFirst=*/false);
    savedFile.close();

    // store the new results in a ref file for easier comparison in case of problems
    const string newRef = fRefFileName + ".after-rereconstruction";
    ofstream fout(newRef.c_str());
    fout << fResults.str();
    fout.close();
  }
  fResults.str("");

  return comp ? eSuccess : eFailure;
}