SdMonteCarloEventSelector.cc

Go to the documentation of this file.

#include <string>

#include <fwk/CentralConfig.h>
#include <sdet/SDetector.h>
#include <utl/Reader.h>
#include <utl/String.h>
#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>
#include "SdMonteCarloEventSelector.h"

using namespace std;
using namespace fwk;
using namespace evt;
using namespace sevt;
using namespace utl;


namespace SdMonteCarloEventSelectorOG {

  void
  MakeInfo(const vector<int>& ids, const string& message)
  {
    if (ids.empty())
      return;
    ostringstream info;
    info << "station" << String::Plural(ids) << ' ' << String::OfSortedIds(ids) << ' ' << message;
    INFO(info);
  }


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

    string arrayLayout;
    topB.GetChild("ArrayLayout").GetData(arrayLayout);
    if (arrayLayout == "KeepOnlyDenseArray")
      fArrayLayout = eKeepOnlyDenseArray;
    else if (arrayLayout == "KeepOnlyRegularArray")
      fArrayLayout = eKeepOnlyRegularArray;
    else
      fArrayLayout = eKeepAll;

    string innerRadiusTreatment;
    topB.GetChild("InnerRadiusTreatment").GetData(innerRadiusTreatment);
    if (innerRadiusTreatment == "KeepInnerRadiusStation")
      fInnerRadiusTreatment = eKeepInnerRadiusStation;
    else if (innerRadiusTreatment == "RejectInnerRadiusStation")
      fInnerRadiusTreatment = eRejectInnerRadiusStation;
    else
      fInnerRadiusTreatment = eSkipEvent;

    return eSuccess;
  }


  VModule::ResultFlag
  SdMonteCarloEventSelector::Run(Event& event)
  {
    const sdet::SDetector& sDetector = det::Detector::GetInstance().GetSDetector();

    sevt::SEvent& sEvent = event.GetSEvent();

    vector<int> svInnerRadiusTreatment;
    vector<int> svArrayLayout;

    for (SEvent::StationIterator sIt = sEvent.StationsBegin();
         sIt != sEvent.StationsEnd(); ++sIt) {

      // Stations with distance in shower plane coordinates smaller than the
      // inner radius used in simulations (CORSIKA/AIRES) are eith removed, kept
      // or the hole event is skipped

      if (sIt->HasSimData()) {
        const sevt::StationSimData& sSim = sIt->GetSimData();
        if (sSim.IsInsideMinRadius()) {
          svInnerRadiusTreatment.push_back(sIt->GetId());
          if (fInnerRadiusTreatment == eRejectInnerRadiusStation)
            sIt->SetRejected(StationConstants::eMCInnerRadiusCut);
        }
      }

      const int sId = sIt->GetId();
      const bool isDense = sDetector.GetStation(sId).IsDense();

      //  removal of dense array stations or regular array stations

      switch (fArrayLayout) {
      case eKeepAll:
        if (isDense && sIt->IsRejected() &&
            sIt->GetRejectionStatus() ==
              (StationConstants::eDenseArray | StationConstants::eOffGrid))
          sIt->SetCandidate();
        break;
      case eKeepOnlyDenseArray:
        if (isDense) {
          if (sIt->IsRejected() &&
              sIt->GetRejectionStatus() ==
              (StationConstants::eDenseArray | StationConstants::eOffGrid)) {
            if (sIt->HasRecData())
              sIt->SetCandidate();
            else
              sIt->SetSilent();
          }
        } else {
          svArrayLayout.push_back(sId);
          sIt->SetRejected(StationConstants::eNotDenseArray);
        }
        break;
      case eKeepOnlyRegularArray:
        if (isDense) {
          svArrayLayout.push_back(sId);
          sIt->SetRejected(StationConstants::eDenseArray);
        }
        break;
      }
    }

    if (!svInnerRadiusTreatment.empty()) {
      MakeInfo(svInnerRadiusTreatment, string("within inner radius cut") +
                                       (fInnerRadiusTreatment == eKeepInnerRadiusStation ?
                                        " kept: " :
                                        (fInnerRadiusTreatment == eRejectInnerRadiusStation ?
                                         " rejected: " : ", event skipped.")));
      if (fInnerRadiusTreatment == eSkipEvent)
        return eContinueLoop;
    }

    if (!svArrayLayout.empty()) {
      MakeInfo(svArrayLayout, string("not in ") +
                              (fArrayLayout == eKeepOnlyDenseArray ?
                               "dense " : "regular ") +
                              "array rejected.");
    } else if (fArrayLayout == eKeepAll)
      INFO("Keeping both dense and regular stations.");

    return eSuccess;
  }


  VModule::ResultFlag
  SdMonteCarloEventSelector::Finish()
  {
    return eSuccess;
  }

}