HdRecValidation.cc

Go to the documentation of this file.

#include <evt/Event.h>
#include <evt/ShowerFRecData.h>

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

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

#include <utl/ErrorLogger.h>
#include <utl/Branch.h>
#include <utl/Particle.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <tst/Validatrix.h>

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

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <det/Detector.h>

#include <iostream>

#include "HdRecValidation.h"

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


namespace HdRecValidationNS {

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

    fOutputName = fRefFileName + string(".new");
    fOutput = new ofstream(fOutputName.c_str());

    return eSuccess;
  }


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

    if (!event.HasFEvent()) {
      ERROR("Event should have FEvent.");
      return eFailure;
    }

    ostringstream eventId;
    eventId << event.GetHeader().GetId();
    *fOutput << BeLabel(eventId.str());

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

    for (const auto& s : sEvent.StationsRange()) {

      if (s.HasRecData()) {
        ostringstream label;
        label << "Station_" << s.GetId();
        *fOutput << BeLabel(label.str());
        *fOutput << BeLabel("signal");
        *fOutput << BeCloseRel(s.GetRecData().GetTotalSignal(), 1*percent);
      }

    }

    *fOutput << BeLabel("FD_part");

    const FEvent& fEvent = event.GetFEvent();
    for (FEvent::ConstEyeIterator eIt = fEvent.EyesBegin(ComponentSelector::eHasData);
         eIt != fEvent.EyesEnd(ComponentSelector::eHasData); ++eIt) {
      ostringstream label;

      const int eyeId = eIt->GetId();
      label << "Eye_" << eyeId << endl;
      *fOutput << BeLabel(label.str());

      const EyeRecData& eRec = eIt->GetRecData();
      const ShowerFRecData& sFRec = eRec.GetFRecShower();

      *fOutput << BeLabel("Rp");
      *fOutput << BeCloseRel(eRec.GetRp()/kilometer, 1*percent);

      const Vector sdp = eRec.GetSDP();
      const Point eyePos = det::Detector::GetInstance().GetFDetector().GetEye(eyeId).GetPosition();
      const CoordinateSystemPtr eyeCS = LocalCoordinateSystem::Create(eyePos);

      *fOutput << BeLabel("theta");
      *fOutput << BeCloseRel(sdp.GetTheta(eyeCS)/degree, 1*percent);

      *fOutput << BeLabel("phi");
      *fOutput << BeCloseRel(sdp.GetPhi(eyeCS)/degree, 1*percent);

      const double Xmax = (sFRec.HasGHParameters() ?
                           sFRec.GetGHParameters().GetXMax()/(g/cm2) : 0);
      if (Xmax > 0) {
        *fOutput << BeLabel("Xmax");
        *fOutput << BeCloseRel(Xmax, 1*percent);
      }

      const double energy = sFRec.GetTotalEnergy();
      if (energy > 0) {
        *fOutput << BeLabel("total_energy");
        *fOutput << BeCloseRel(energy/EeV, 1*percent);
      }

      const double energyEM = sFRec.GetEmEnergy();
      if (energyEM > 0) {
        *fOutput << BeLabel("EM_energy");
        *fOutput << BeCloseRel(energyEM/EeV, 1*percent);
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  HdRecValidation::Finish()
  {
    delete fOutput;
    fOutput = nullptr;

    ifstream oldFile(fRefFileName.c_str());
    ifstream newFile(fOutputName.c_str());

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

    return comp ? eSuccess : eFailure;
  }

}