SDenseStationListXMLManager.cc

Go to the documentation of this file.

#include <fwk/CentralConfig.h>
#include <det/VManager.h>
#include <sdet/SManagerRegister.h> // needed for registration macro

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/UTMPoint.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/Line.h>
#include <utl/Plane.h>
#include <utl/GeometryUtilities.h>

#include <fwk/CoordinateSystemRegistry.h>

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

#include <fwk/RunController.h>
#include <fwk/LocalCoordinateSystem.h>

#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <sdet/SDenseStationListXMLManager.h>

#include <boost/range/iterator_range.hpp>

using namespace std;
using namespace fwk;
using namespace utl;
using namespace det;
using namespace evt;


REGISTER_S_MANAGER("SDenseStationListXMLManager", sdet::SDenseStationListXMLManager);


namespace sdet {

  void
  SDenseStationListXMLManager::Init(const string& configLink)
  {
    // call shadowed VManager::Init() first
    VManager::Init(configLink);

    fDenseStationList.clear();

    const string stationString = "denseStation";

    for (auto cB = fBranch.GetFirstChild(); cB; cB = cB.GetNextSibling()) {
      if (cB.GetName() != stationString)
        continue;
      const auto atts = cB.GetAttributes();
      const auto id = FindComponent<int>("id", atts);
      const auto it = fDenseStationList.find(id);
      if (it == fDenseStationList.end())
        fDenseStationList.insert(fDenseStationList.end(), id);
      else {
        ostringstream warn;
        warn << "Found multiple instances of dense station " << id << ". "
                "Only the first one will be used.";
        WARNING(warn);
      }
    }
  }


  VManager::Status
  SDenseStationListXMLManager::GetData(int& returnData,
                                       const string& componentProperty,
                                       const string& componentName,
                                       const IndexMap& componentIndex)
    const
  {
    if (componentName != "stationList" || !IsDense(componentIndex))
      return eNotFound;

    if (componentProperty == "isUUB") {
      returnData = 0;
      auto isUUBB = fBranch.GetChild("isUUB");
      if (isUUBB)
        isUUBB.GetData(returnData);
    } else if (componentProperty == "hasSmallPMT") {
      returnData = 0;
      auto hasSmallPMTB = fBranch.GetChild("hasSmallPMT");
      if (hasSmallPMTB)
        hasSmallPMTB.GetData(returnData);
    } else if (componentProperty == "hasScintillator") {
      returnData = 0;
      auto hasScintillatorB = fBranch.GetChild("hasScintillator");
      if (hasScintillatorB)
        hasScintillatorB.GetData(returnData);
    } else if (componentProperty == "groupId")
      returnData = 0;
    else if (componentProperty == "zone")
      fBranch.GetChild("UTMInfo").GetChild("zone").GetData(returnData);
    else
      return eNotFound;

    // final test: throw if "stationId" is not present in componentIndex
    FindComponent<string>("stationId", componentIndex);
    return eFound;
  }


  VManager::Status
  SDenseStationListXMLManager::GetData(double& returnData,
                                       const string& componentProperty,
                                       const string& componentName,
                                       const IndexMap& componentIndex)
    const
  {
    if (componentName != "stationList" || !IsDense(componentIndex))
      return eNotFound;

    if (componentProperty == "northing" ||
        componentProperty == "easting" ||
        componentProperty == "altitude") {

      const auto& event = RunController::GetInstance().GetCurrentEvent();

      fDenseStations.Validate(event);

      const auto id = FindComponent<int>("stationId", componentIndex);

      if (!fDenseStations.HasStation(id)) {
        // Dense station positions specified in (2-D) polar coordinates in the
        // showerCS.  Here we project coordinates onto the corsika observation
        // level, and convert resulting coordinates into UTM so that that
        // sdet::Station will interpret the data as though it were a regular
        // (non-dense) station.

        double deltaPhi = 0;
        if (GetData(deltaPhi, "deltaPhi", "stationList", componentIndex) == eNotFound)
          return eNotFound;
        double r = 0;
        if (GetData(r, "r", "stationList", componentIndex) == eNotFound)
          return eNotFound;

        fDenseStations.AddStation(id, r, deltaPhi);
      }

      const auto& pos = fDenseStations.GetStation(id);
      if (componentProperty == "northing")
        returnData = pos.fNorthing;
      else if (componentProperty == "easting")
        returnData = pos.fEasting;
      else if (componentProperty == "altitude")
        returnData = pos.fHeight;
      else
        return eNotFound;

      return eFound;

    } else if (componentProperty == "axis1" || componentProperty == "axis2") {

      returnData = 0;
      return eFound;

    } else {

      const auto it = componentIndex.find("stationId");

      if (it != componentIndex.end()) {

        const IndexMap atts({{ "id", it->second }});

        auto stationB = fBranch.GetChild("denseStation", atts);
        if (!stationB)
          return eNotFound;

        auto componentB = stationB.GetChild(componentProperty);
        if (!componentB)
          return eNotFound;

        componentB.GetData(returnData);
        return eFound;

      } else {

        ostringstream info;
        info << "Could not interpret the request with "
                "component property = '" << componentProperty << "', "
                "component name = '" << componentName << "', index ";
        for (const auto& kv : componentIndex)
          info << " (" << kv.first << ", " << kv.second << ')';
        INFO(info);
        return eNotFound;

      }

    }
  }


  VManager::Status
  SDenseStationListXMLManager::GetData(string& returnData,
                                       const string& componentProperty,
                                       const string& componentName,
                                       const IndexMap& componentIndex)
    const
  {
    if (componentName != "stationList" || !IsDense(componentIndex))
      return eNotFound;

    if (componentProperty == "name")
      returnData = "DenseStation" + FindComponent<string>("stationId", componentIndex);
    else if (componentProperty == "commission")
      returnData = "2000-01-01T00:00:00Z";
    else if (componentProperty == "decommission")
      returnData = "2030-01-01T00:00:00Z";
    else {
      auto componentB = fBranch.GetChild("UTMInfo").GetChild(componentProperty);
      if (componentB)
        componentB.GetData(returnData);
      else
        return eNotFound;
    }

    return eFound;
  }


  VManager::Status
  SDenseStationListXMLManager::GetData(vector<int>& returnData,
                                       const string& componentProperty,
                                       const string& componentName,
                                       const IndexMap& componentIndex)
    const
  {
    if (componentName != "stationList")
      return eNotFound;

    if (componentProperty == "denseStationList") {
      GetFullStationList(returnData);
      return eFound;
    } else if (componentProperty == "stationGroup") {
      if (IsDense(componentIndex)) {
        returnData.clear();
        return eFound;
      }
    } else if (componentProperty == "crown") {
      if (IsDense(componentIndex)) {
        returnData.clear();
        return eFound;
      }
    }

    return eNotFound;
  }


  VManager::Status
  SDenseStationListXMLManager::GetData(vector<bool>& returnData,
                                       const string& componentProperty,
                                       const string& componentName,
                                       const IndexMap& componentIndex)
    const
  {
    if (componentName != "stationList" || !IsDense(componentIndex))
      return eNotFound;

    if (componentProperty == "inGrid") {
      returnData.clear();
      return eFound;
    }

    return eNotFound;
  }


  void
  SDenseStationListXMLManager::GetFullStationList(vector<int>& returnList)
    const
  {
    returnList.assign(fDenseStationList.begin(), fDenseStationList.end());
  }


  bool
  SDenseStationListXMLManager::IsDense(const IndexMap& componentIndex)
    const
  {
    const auto id = FindComponent<int>("stationId", componentIndex);
    return fDenseStationList.find(id) != fDenseStationList.end();
  }


  //


  void
  SDenseStationListXMLManager::DenseStations::Validate(const Event& event)
  {
    Validate(event.GetHeader().GetId());
    if (fEventId.empty())
      Init(event);
  }


  void
  SDenseStationListXMLManager::DenseStations::Init(const Event& event)
  {
    // dense tank positions computed relative to sim shower. check
    // that the sim shower has all the necessary ingredients.

    if (!event.HasSimShower())
      throw MissingEventDataException(
        "Attempt to generate a dense array failed because "
        "no simulated shower is present in the event."
      );

    const auto& simShower = event.GetSimShower();

    if (!simShower.HasGeometry())
      throw MissingEventDataException(
        "Attempt to generate a dense array failed because "
        "the shower core position has not been set."
      );

    const auto& showerCS = simShower.GetShowerCoordinateSystem();
    const Point observationCore(0,0,0, showerCS);
    fShowerDirection = simShower.GetDirection();

    // The core has to be projected from the corsika observation-level
    // to the array ground. Since we do not know which array will be used
    // (ideal or real), we simply find a station which is closest to the shower
    // axis and take its height as the reference height.

    const Line showerAxis(observationCore, fShowerDirection);
    const auto& sDet = det::Detector::GetInstance().GetSDetector();
    const sdet::Station* closestStation = nullptr;
    {
      double minDistance = numeric_limits<double>::infinity();
      for (const auto& s : sDet.StationsRange()) {
        if (s.IsDense())
          continue;
        const double d = Distance(showerAxis, s.GetPosition());
        if (d < minDistance) {
          minDistance = d;
          closestStation = &s;
        }
      }
    }

    if (!closestStation)
      throw MissingEventDataException("Could not find closest station!");
    else {
      ostringstream info;
      info << "Reference height for dense stations taken from closest station "
              "id = " << closestStation->GetId();
      INFO(info);
    }

    // project shower axis (line) to the ground plane at height of the closest station

    const auto& refPos = closestStation->GetPosition();
    const auto& groundCS = LocalCoordinateSystem::Create(refPos);
    const Vector groundNormal(0,0,1, groundCS);
    fGroundPlane = Plane(refPos, groundNormal);
    fGroundCore = Intersection(fGroundPlane, showerAxis);
    fBaseX = Vector(1,0,0, showerCS);
    fBaseY = Vector(0,1,0, showerCS);

#if 0
    const auto& pampaCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    ostringstream info;
#define DUMP(_x_) #_x_ " (" << _x_.GetX(pampaCS) << ", " << _x_.GetY(pampaCS) << ", " << _x_.GetZ(pampaCS) << ")\n"
    info << "trigger calculation:\n"
         << DUMP(observationCore)
         << DUMP(fShowerDirection)
         << DUMP(refPos)
         << DUMP(groundNormal)
         << DUMP(fGroundCore)
         << DUMP(fBaseX)
         << DUMP(fBaseY)
      ;
#undef DUMP
    INFO(info);
#endif

    fEventId = event.GetHeader().GetId();
  }


  void
  SDenseStationListXMLManager::DenseStations::AddStation(const int id, const double r, const double deltaPhi)
  {
    // add the two unit vectors according to phi angle and stretch by r
    const double x = r * cos(deltaPhi);
    const double y = r * sin(deltaPhi);
    const Vector rho = x * fBaseX + y * fBaseY;
    const Point stationInShowerPlane = fGroundCore + rho;
    const Point stationInGroundPlane =
      Intersection(fGroundPlane, Line(stationInShowerPlane, fShowerDirection));
    const UTMPoint utm(stationInGroundPlane, ReferenceEllipsoid::eWGS84);
    const UTM neh{utm.GetNorthing(), utm.GetEasting(), utm.GetHeight()};
    fStations.insert({id, neh});
  }

}