HybridSimValidation.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/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 "HybridSimValidation.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 HybridSimValidationNS;


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

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

  return eSuccess;
}


VModule::ResultFlag
HybridSimValidation::Run(evt::Event& event)
{
  ostringstream eventId;
  eventId << event.GetHeader().GetId();
  *fOutput << BeLabel(eventId.str());

  // Station/telescope level rec data

  // SD
  if (!event.HasSEvent()) {
    ERROR("Event does not contain an SEvent.");
    return eFailure;
  }

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

  const FEvent& fEvent = event.GetFEvent();

  *fOutput << 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();
      *fOutput << BeLabel(eLabel.str());

      if (eIt->HasRecData()) {
        if (eIt->GetRecData().HasFRecShower()) {
          const ShowerFRecData& showerRec = eIt->GetRecData().GetFRecShower();
          *fOutput << BeLabel("Stations_in_fit");
          *fOutput << BeEqual(showerRec.GetStationIds());

          *fOutput << BeLabel("TotalEnergy(EeV)");
          *fOutput << BeCloseRel(showerRec.GetTotalEnergy()/EeV, /*tolerance*/(eIt->GetId() == 4 ? 0.2 : 0.1));

          *fOutput << BeLabel("Xmax(g/cm2)");
          *fOutput << BeCloseAbs(showerRec.GetGHParameters().GetXMax()/(g/cm2), (eIt->GetId() == 4 ? 30 : 15));
        }
      }
    }

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

      if (tIt->HasRecData()) {
        ostringstream tLabel;
        tLabel << "Telescope_" << tIt->GetId();
        *fOutput << BeLabel(tLabel.str());

        const TelescopeRecData& tRec = tIt->GetRecData();
        *fOutput << BeLabel("Rp");
        *fOutput << BeCloseRel(tRec.GetRp(), /*tolerance*/0.05);
      }
    }
  }

  *fOutput << 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");

  *fOutput << BeLabel("Shower_S_rec");

  const ShowerSRecData& showerSRec = showerRec.GetSRecShower();

  *fOutput << BeLabel("S1000");
  *fOutput << BeCloseRel(showerSRec.GetShowerSize(), 0.1);

  return eSuccess;
}


VModule::ResultFlag
HybridSimValidation::Finish()
{
  fOutput->close();
  delete fOutput;

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

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

  oldFile.close();
  newFile.close();

  return comp ? eSuccess : eFailure;
}