Sd5T5Tester.cc

Go to the documentation of this file.

#include "Sd5T5Tester.h"
#include <utl/ErrorLogger.h>
#include <utl/AugerException.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/AugerUnits.h>
#include <utl/Branch.h>

#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <sstream>
#include <iostream>
#include <vector>

using namespace std;
using namespace utl;
using namespace fwk;

#define DEBUG(x) std::cout << #x << " " << x << std::endl


namespace Sd5T5Tester {

  static const int kNone = -1;


  VModule::ResultFlag
  Sd5T5Tester::Init()
  {
    // initialize internal variables
    fCounter = 0;
    for (int i = 0; i < 7; ++i)
      fCrownId[i] = kNone;

    // get constants to compute first crown from SdEventSelector
    CentralConfig* const cc = CentralConfig::GetInstance();
    Branch topB = cc->GetTopBranch("SdEventSelector");

    topB.GetChild("FirstCrownDistanceRange").GetData(fFirstCrownDistanceRange);

    return eSuccess;
  }


  VModule::ResultFlag
  Sd5T5Tester::Run(evt::Event& event)
  {
    if (!event.HasSEvent() ||
        !event.HasRecShower() ||
        !event.GetRecShower().HasSRecShower())
      return eContinueLoop;

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

    ostringstream msg;
    msg << "Counter = " << fCounter;
    INFO(msg);

    if (fCounter == 0) {
      // check whether we should use this event
      // and then reconstruct unmodified event

      const sdet::Station& sCenter = GetHottestStation(sEvent);

      const Crown crown = GetCrown(sEvent, sCenter);

      // skip events which already have missing stations in first or second crow
      if (!IsAccepted(sEvent, crown))
        return eBreakLoop;

      UpdateFirstCrown(sCenter, crown);
      // continue to reconstruct unmodified event

      ++fCounter;
      return eSuccess;
    }

    if (fCounter == 8) {
      // we now tried all possibilities for
      // current event, now exit loop and
      // start anew
      fCounter = 0;
      return eBreakLoop;
    }

    // modify event for fCounter 1 to 7

    // modify Auger ID to mark the modified events
    ostringstream augerId;
    augerId << event.GetHeader().GetId() << "-" << fCounter;
    event.GetHeader().SetId(augerId.str());

    // turn off next station, reactivate others
    for (int iStation = 1; iStation <= 7; ++iStation) {
      const int stationId = fCrownId[iStation - 1];
      if (stationId == kNone)
        continue;

      sevt::Station& station = event.GetSEvent().GetStation(stationId);

      if (iStation == fCounter) {
        ostringstream msg;
        msg << "Disabling station " << stationId;
        INFO(msg);
        station.SetRejected(sevt::StationConstants::eOffGrid);
      } else
        station.ClearRejectionStatus();
    }

    // count up for next iteration
    ++fCounter;
    return eSuccess;
  }


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


  const sdet::Station&
  Sd5T5Tester::GetHottestStation(const sevt::SEvent& sEvent)
    const
  {
    // find candidate station with largest signal

    const sevt::Station* maxStation = 0;

    double maxSignal = 0;
    for (sevt::SEvent::ConstCandidateStationIterator sIt = sEvent.CandidateStationsBegin(),
         sEnd = sEvent.CandidateStationsEnd(); sIt != sEnd; ++sIt) {
      const double signal = sIt->GetRecData().GetTotalSignal();
      if (maxSignal < signal) {
        maxSignal = signal;
        maxStation = &(*sIt);
      }
    }

    if (!maxStation)
      throw utl::DoesNotComputeException("Sd5T5Tester::GetHottestStation failed");

    return det::Detector::GetInstance().GetSDetector().GetStation(*maxStation);
  }


  Sd5T5Tester::Crown
  Sd5T5Tester::GetCrown(const sevt::SEvent& sEvent,
                        const sdet::Station& sCenter)
    const
  {
    // find direct neighbours that form a hexagon
    // aka "crown" around centered station,
    // following code in SdEventSelector

    Crown result;

    CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();

    // assign station to index according to its orientation
    const sdet::SDetector& sDet = det::Detector::GetInstance().GetSDetector();
    for (sevt::SEvent::ConstStationIterator sIt = sEvent.StationsBegin(),
         sEnd = sEvent.StationsEnd(); sIt != sEnd; ++sIt) {
      if (!(sIt->IsCandidate() || sIt->IsSilent()))
        continue;

      const sdet::Station& s = sDet.GetStation(*sIt);

      const Vector r = s.GetPosition() - sCenter.GetPosition();

      const double dist = r.GetMag();
      if (dist < fFirstCrownDistanceRange[0] ||
          fFirstCrownDistanceRange[1] < dist)
        continue;

      result.push_back(&s);
    }

    return result;
  }


  bool
  Sd5T5Tester::IsAccepted(const sevt::SEvent& sev, const Crown& crown)
    const
  {
    size_t nCandidates = 0;
    for (sevt::SEvent::ConstCandidateStationIterator it = sev.CandidateStationsBegin(),
         end = sev.CandidateStationsEnd(); it != end; ++it)
      ++nCandidates;

    if (nCandidates < 4)
      return false;

    if (crown.size() != 6)
      return false;

    for (Crown::const_iterator it = crown.begin();
         it != crown.end(); ++it) {
      const Crown c = GetCrown(sev, **it);
      if (c.size() != 6)
        return false;
    }

    return true;
  }


  void
  Sd5T5Tester::UpdateFirstCrown(const sdet::Station& sCenter, const Crown& crown)
  {
    // find direct neighbours that form a hexagon
    // aka "crown" around centered station,
    // following code in SdEventSelector

    // reset to list with dummy values
    for (int i = 0; i < 7; ++i)
      fCrownId[i] = kNone;

    fCrownId[6] = sCenter.GetId();

    CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();

    for (size_t iS = 0; iS < 6; ++iS) {
      const sdet::Station& sCrown = *crown[iS];

      const Vector r = sCrown.GetPosition() - sCenter.GetPosition();

      const int phi = r.GetPhi(siteCS) / degree;
      const size_t index = ((phi - 30 + 720) % 360) / 60;

      fCrownId[index] = sCrown.GetId();
    }
  }

}