ReReconstruction/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 <sevt/PMTSimData.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 <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;


namespace ValidationNS {

  VModule::ResultFlag
  Validation::Init()
  {
    INFO("Doing the validation !!!!");
    Branch topB = CentralConfig::GetInstance()->GetTopBranch("Validation");
    topB.GetChild("RefFilename").GetData(fRefFileName);
    const string mode = topB.GetChild("Mode").Get<string>();

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

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

    fResults = new stringstream;

    return eSuccess;
  }


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

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

    *fResults << BeLabel("SD_rec");

    for (SEvent::ConstStationIterator sIt = sEvent.StationsBegin();
         sIt != sEvent.StationsEnd(); ++sIt) {
      if (sIt->HasRecData()) {
        ostringstream label;
        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()) {
        ostringstream eLabel;
        eLabel << "Eye_" << eIt->GetId();
        *fResults << BeLabel(eLabel.str());

        if (eIt->GetRecData().HasFRecShower()) {
          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*/5e-2);

          *fResults << BeLabel("Xmax(g/cm2)");
          *fResults << BeCloseRel(showerRec.HasGHParameters()?
                                  showerRec.GetGHParameters().GetXMax()/gram*cm2:
                                  0., 1e-5);
          
          *fResults << BeLabel("Xmax2Error");
          *fResults << BeCloseRel(showerRec.HasGHParameters()?
                                  showerRec.GetGHParameters().GetXMaxError()/gram*cm2 : 0., 5e-2);
        }
      }
      
      for (Eye::ConstTelescopeIterator tIt = eIt->TelescopesBegin(ComponentSelector::eHasData);
           tIt != eIt->TelescopesEnd(ComponentSelector::eHasData); ++tIt) {
        if (tIt->HasRecData()) {
          ostringstream tLabel;
          tLabel << "Telescope_" << tIt->GetId();
          *fResults << BeLabel(tLabel.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);

    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
      string newRef = fRefFileName + ".after-rereconstruction";
      ofstream fout(newRef.c_str());
      fout << fResults->str();
      fout.close();
    }

    delete fResults;
    fResults = nullptr;

    return comp ? eSuccess : eFailure;
  }

}