MdSdReReconstruction/Validation.cc

Go to the documentation of this file.

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

#include <sevt/SEvent.h>
#include <sevt/SortCriteria.h>
#include <sevt/PMTSimData.h>
#include <mevt/MEvent.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 mevt;
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) {
    /*********************
    Shower-Rec Validation
    **********************/
    *fResults << BeLabel("Shower_rec");

    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("S450");
    *fResults << BeEqual(showerSRec.GetShowerSize());
    *fResults << BeLabel("S450e");
    *fResults << BeCloseAbs(showerSRec.GetShowerSizeError(), /*tolerance*/ 1e-4);

    // MD rec (shower level)
    if (!showerRec.HasMRecShower())
      ERROR("RecShower does not have MRecShower");

    *fResults << "\n";
    *fResults << BeLabel("Shower_M_rec");

    const ShowerMRecData & showerMRec = showerRec.GetMRecShower();

    *fResults << BeLabel("rho450");
    *fResults << BeCloseAbs(showerMRec.GetNMuRef(), /*tolerance*/ 1e-5);
    *fResults << BeLabel("rho450e");
    *fResults << BeCloseAbs(showerMRec.GetNMuRef(), /*tolerance*/ 1e-4);
    *fResults << "\n";

    /**************************
    SD Station-Level Validation
    ***************************/
    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 << BeEqual(sIt -> GetRecData().GetTotalSignal());
      }
    }

    /**************************
    MD Counter-Level Validation
    ***************************/
    if (!event.HasMEvent()) {
      ERROR("Event should contain an MEvent.");
      return eFailure;
    }

    const MEvent & mEvent = event.GetMEvent();

    *fResults << "\n";
    *fResults << BeLabel("MD_rec");

    for (SEvent::ConstStationIterator sIt = sEvent.StationsBegin(); sIt != sEvent.StationsEnd(); ++sIt) {
      if (sIt -> HasRecData()) {
        // getting the counter 
        const mevt::Counter & counter = mEvent.GetCounter(1e5 + sIt -> GetId());
        // if at least one of the modules is a candidate
        if (!counter.IsRejected()) {

          string counterLabel;
          counterLabel = "Counter_" + to_string(counter.GetId());

          //Vectors to allocate the number of muons per module
          vector<double> muonNumber(3, 0);
          vector<double> muonNumberMC(3, 0);

          // Iteration over module mIt of station sdId
          for (mevt::Counter::ModuleConstIterator mIt = counter.ModulesBegin(); mIt != counter.ModulesEnd(); ++mIt) {

            int mId = mIt -> GetId();

            if (mIt -> IsRejected()) {
              muonNumberMC.at(mId - 1) = -1;
              muonNumber.at(mId - 1) = -1;
              continue;
            }

            // Iteration over scillantors of module mdId -> to get the injected muons, i.e. True Label
            for (mevt::Module::ScintillatorConstIterator sIt = mIt -> ScintillatorsBegin(); sIt != mIt -> ScintillatorsEnd(); ++sIt) {
              if (sIt -> HasSimData()) {
                const mevt::ScintillatorSimData & ssd = sIt -> GetSimData();
                muonNumberMC.at(mId - 1) += ssd.GetNumberOfInjectedMuons();
              }
            } // end of iteration over scillantors of module mdId 

            if (mIt -> HasRecData()) {

              // Reconstructed number of muons at module level
              const mevt::ModuleRecData & mRecData = mIt -> GetRecData();

              if (mRecData.IsSaturated())
                muonNumber.at(mId - 1) += mRecData.GetNumberOfMuonsLowLimit();
              else
                muonNumber.at(mId - 1) += mRecData.GetNumberOfEstimatedMuons();

            } else
              muonNumber.at(mId - 1) = -1;

          } // end of teration over module mIt of station sdId

        for (size_t i = 0; i < muonNumber.size(); ++i) { 
            *fResults << BeLabel(counterLabel + "_Module_" + to_string(i) + "_rec");
            *fResults << BeCloseAbs(muonNumber.at(i), /*tolerance*/ 1e-4); // Access the corresponding element in vector 'v' using iterator difference
            *fResults << BeLabel(counterLabel + "_Module_" + to_string(i) + "_MC"); 
            *fResults << BeEqual(muonNumberMC.at(i)); // Access the corresponding element in vector 'vMC' using iterator difference
          }
        }
      }
    }

    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;
  }
}