MCShowerChecker.cc

Go to the documentation of this file.

#include "MCShowerChecker.h"

#include <fwk/CentralConfig.h>

#include <det/Detector.h>

#include <evt/Header.h>
#include <evt/Event.h>
#include <evt/ShowerSimData.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/TimeStamp.h>
#include <utl/TabularStream.h>
#include <utl/TabulatedFunction.h>
#include <utl/AugerException.h>

#include <cstddef>
#include <cstdio>
#include <sstream>

using namespace MCShowerCheckerOG;
using namespace fwk;
using namespace evt;
using namespace utl;
using namespace std;


VModule::ResultFlag
MCShowerChecker::Init()
{
  Branch topBranch = CentralConfig::GetInstance()->GetTopBranch ("MCShowerChecker");

  topBranch.GetChild("verbose").GetData(fVerbose);

  const string input = topBranch.GetChild("ElectronNumberProfile").Get<string>();
  fElectronNumberProfile = ToEnum(input);

  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "               electron number profile: " << ToText(fElectronNumberProfile)
       << endl;

  return eSuccess;
}


std::string
MCShowerChecker::ToText(const EHandle flag)
  const
{
  switch (flag) {
  case eReject:
    return "reject";
    break;
  case eAssert:
    return "assert";
    break;
  case eCorrect:
    return "correct";
    break;
  }
  return "unknown";
}


MCShowerChecker::EHandle
MCShowerChecker::ToEnum(const string& text)
  const
{
  if (text == "reject")
    return eReject;
  if (text == "correct")
    return eCorrect;
  if (text == "assert")
    return eAssert;
  throw InvalidConfigurationException("MCShowerChecker:: encountered unknown configuration string: '" + text + "'");
}


VModule::ResultFlag
MCShowerChecker::Run(evt::Event& event)
{
  ++fNEventsProcessed;

  if (!event.HasSimShower())
    return eContinueLoop;
  evt::ShowerSimData& simShower = event.GetSimShower();

  if (!simShower.HasLongitudinalProfile())
    return eSuccess;

  TabulatedFunction& chargeProfile = simShower.GetLongitudinalProfile();
  bool foundNeProblem = false;
  for (unsigned int iPoint = 0; iPoint < chargeProfile.GetNPoints(); ++iPoint) {

    if (chargeProfile.GetY(iPoint) < 0) {

      foundNeProblem = true;

      switch (fElectronNumberProfile) {
      case eAssert:
        ERROR("Found Monte Carlo shower information with negative entries in electron number profile!");
        return eFailure;
      case eReject:
        ++fNEventsWithElectronProfileProblem;
        WARNING("Rejecting Monte Carlo shower: negative entries in electron number profile!");
        return eContinueLoop;
      case eCorrect:
        {
          ostringstream info;
          info << "Correcting negative entries in electron profile of Monte Carlo shower: bin="
               << iPoint
               << " X=" << chargeProfile.GetX(iPoint) / (g/cm2)
               << " Ne=" << chargeProfile.GetY(iPoint);
          WARNING(info);
          chargeProfile.GetY(iPoint) = 0;
        }
        break;
      }
    }
  }
  if (foundNeProblem)
    ++fNEventsWithElectronProfileProblem;

  ++fNEventsAccepted;

  return eSuccess;
}


VModule::ResultFlag
MCShowerChecker::Finish()
{
  if (fVerbose > 0) {

    if (fNEventsProcessed) {

      const int nDiscarded = fNEventsProcessed - fNEventsAccepted;

      TabularStream tab("r  r  .");

      tab << endc << 'N' << endc << '%' << endr
          << "Events processed:" << endc << fNEventsProcessed << endc << 100 << endr
          << "Events skipped:" << endc << nDiscarded << endc << int(nDiscarded*1000./fNEventsProcessed)/10. << endr
          << "Ne-profile:" << endc << fNEventsWithElectronProfileProblem << endc
                           << int(fNEventsWithElectronProfileProblem*1000./fNEventsProcessed)/10.;
      ostringstream info;
      info << "\n\n" << tab;
      INFO(info);

    } else
      INFO("No events processed.");

  }
  return eSuccess;
}