SDetector/Station.cc

Go to the documentation of this file.

#include <sdet/Station.h>
#include <sdet/SDetector.h>
#include <mdet/MDetector.h>
#include <utl/GeometryUtilities.h>
#include <utl/UTMPoint.h>
#include <utl/Test.h>
#include <utl/UTCDateTime.h>
#include <sevt/PMT.h>

#include <boost/lexical_cast.hpp>
#include <boost/dynamic_bitset.hpp>

#include <iostream>
#include <sstream>
#include <bitset>

using namespace det;
using namespace std;
using namespace utl;


namespace sdet {

  const Station&
  Station::InternalPartnerStationFunctor::operator()(const int stationId)
    const
  {
    return Detector::GetInstance().GetSDetector().GetStation(stationId);
  }


  Station::Station(const int id) :
    fId(id)
  {
    /* Adding station id and hardware version to station index map which is
     * passed to the managers with each GetStationData request (or when
     * Handle is used). */
    fStationIndexMap["stationId"] = boost::lexical_cast<string>(id);
    fStationIndexMap["isUUB"] = boost::lexical_cast<string>(fIsUUB);

    // Constructing PMTs
    for (unsigned int i = kFirstWaterCherenkovLargePMTId;
         i <= kLastWaterCherenkovLargePMTId; ++i) {
      const auto pmt = new PMT(i, id, IsUUB(), PMTConstants::eWaterCherenkovLarge);
      fPMTVector.push_back(pmt);
    }
  }


  Station::~Station()
  {
    for (const auto p : fPMTVector)
      delete p;
  }


  bool
  Station::ExistsAssociatedCounter()
    const
  {
    return Detector::GetInstance().GetMDetector().ExistsAssociatedCounter(GetId());
  }


  int
  Station::GetAssociatedCounterId()
    const
  {
    return Detector::GetInstance().GetMDetector().GetAssociatedCounterId(GetId());
  }


  const PMT&
  Station::GetPMT(const int id)
    const
  {
    for (const auto& p : fPMTVector)
      if (p->GetId() == id)
        return *p;
    ostringstream err;
    err << "Invalid request for PMT with Id: " << id;
    throw NonExistentComponentException(err.str());
  }


  const PMT&
  Station::GetPMT(const sevt::PMT& pmt)
    const
  {
    return GetPMT(pmt.GetId());
  }


  const PMT&
  Station::GetSmallPMT()
    const
  {
    // Check to ensure that there is only one PMT of the given type
    if (GetNPMTs(sdet::PMTConstants::eWaterCherenkovSmall) > 1) {
      std::ostringstream err;
      err << "The station appears to contain more than one small PMT"
             "which is neither expected nor allowed.";
      ERROR(err);
      throw utl::DoesNotComputeException(err.str());
    }
    // Find and return PMT of requested type
    for (const auto p : fPMTVector)
      if (p->GetType() == sdet::PMTConstants::eWaterCherenkovSmall)
        return *p;
    std::ostringstream err;
    err << "Small PMT does not exist.";
    ERROR(err);
    throw utl::NonExistentComponentException(err.str());
  }


  const PMT&
  Station::GetScintillatorPMT()
    const
  {
    // Check to ensure that there is only one PMT of the given type
    if (GetNPMTs(sdet::PMTConstants::eScintillator) > 1) {
      std::ostringstream err;
      err << "The station appears to contain more than one scintillator PMT"
             "which is neither expected nor allowed.";
      ERROR(err);
      throw utl::DoesNotComputeException(err.str());
    }
    // Find and return PMT of requested type
    for (const auto p : fPMTVector)
      if (p->GetType() == sdet::PMTConstants::eScintillator)
        return *p;
    std::ostringstream err;
    err << "Scintillator PMT does not exist.";
    ERROR(err);
    throw utl::NonExistentComponentException(err.str());
  }


  const string&
  Station::GetName()
    const
  {
    if (!fName.IsValid()) {
      GetStationData(fName.Get(), "name", "stationList", "station name");
      fName.SetValid();
    }
    return fName.Get();
  }


  unsigned int
  Station::GetFADCTraceLength()
    const
  {
    if (!fFADCTraceLength.IsValid()) {
      const det::ManagerRegister& m = det::Detector::GetInstance().GetSManagerRegister();
      m.GetDataOrThrow(fFADCTraceLength.Get(), "electronics", "traceBins", fStationIndexMap);
      fFADCTraceLength.SetValid();
    }
    return fFADCTraceLength.Get();
  }


  double
  Station::GetFADCBinSize()
    const
  {
    if (!fFADCBinSize.IsValid()) {
      double frequency = 0;
      GetStationData(frequency, "frequency", "electronics", "fadc electronics frequency");
      fFADCBinSize = 1 / frequency;
    }
    return fFADCBinSize.Get();
  }


  double
  Station::GetTimingUncertainty()
    const
  {
    if (!fTimingUncertainty.IsValid()) {
      GetStationData(fTimingUncertainty.Get(), "timingUncertainty", "electronics", "timing uncertainty");
      fTimingUncertainty.SetValid();
    }
    return fTimingUncertainty.Get();
  }


  unsigned int
  Station::GetSaturationValue()
    const
  {
    if (!fSaturationValue.IsValid()) {
      int bits = 0;
      GetStationData(bits, "bits", "electronics", "fadc bits");
      fSaturationValue = (1 << bits) - 1;
    }
    return fSaturationValue.Get();
  }


  int
  Station::GetLatchBin()
    const
  {
    if (!fLatchBin.IsValid()) {
      GetStationData(fLatchBin.Get(), "latchBin", "electronics", "latch bin of trace");
      fLatchBin.SetValid();
    }
    return fLatchBin.Get();
  }


  const TabulatedFunction&
  Station::GetElectronicsImpulseResponse()
    const
  {
    if (!fElectronicsImpulseResponse) {
      fElectronicsImpulseResponse = new TabulatedFunction;
      GetStationData(*fElectronicsImpulseResponse, "impulseResponse", "electronics", "electronics impulse response");
    }
    return *fElectronicsImpulseResponse;
  }


  double
  Station::GetRadius()
    const
  {
    if (!fRadius.IsValid()) {
      GetStationData(fRadius.Get(), "insideRadius", "tankCasing", "tank casing inside radius");
      fRadius.SetValid();
    }
    return fRadius.Get();
  }


  double
  Station::GetHeight()
    const
  {
    if (!fHeight.IsValid()) {
      GetStationData(fHeight.Get(), "insideHeight", "tankCasing", "tank casing inside height");
      fHeight.SetValid();
    }
    return fHeight.Get();
  }


  double
  Station::GetThickness()
    const
  {
    if (!fThickness.IsValid()) {
      GetStationData(fThickness.Get(), "casingThickness", "tankCasing", "tank casing wall thickness");
      fThickness.SetValid();
    }
    return fThickness.Get();
  }


  const TabulatedFunction&
  Station::GetWaterAbsorptionLength()
    const
  {
    if (!fWaterAbsLength) {
      fWaterAbsLength = new TabulatedFunction;
      GetStationData(*fWaterAbsLength, "waterAbsorptionLength", "water", "water absorption length");
    }
    return *fWaterAbsLength;
  }


  const TabulatedFunction&
  Station::GetWaterRefractionIndex()
    const
  {
    if (!fWaterRefractionIndex) {
      fWaterRefractionIndex = new TabulatedFunction;
      GetStationData(*fWaterRefractionIndex, "waterRefractionIndex", "water", "water refraction index");
    }
    return *fWaterRefractionIndex;
  }


  const TabulatedFunction&
  Station::GetLinerReflectivity()
    const
  {
    if (!fLinerReflectivity) {
      fLinerReflectivity = new TabulatedFunction;
      GetStationData(*fLinerReflectivity, "reflectivity", "liner", "liner reflectivity");
    }
    return *fLinerReflectivity;
  }


  const TabulatedFunction&
  Station::GetLinerSpecularLobe()
    const
  {
    if (!fLinerSpecularLobe) {
      fLinerSpecularLobe = new TabulatedFunction;
      GetStationData(*fLinerSpecularLobe, "specularLobe", "liner", "liner specular lobe constant");
    }
    return *fLinerSpecularLobe;
  }


  const TabulatedFunction&
  Station::GetLinerSpecularSpike()
    const
  {
    if (!fLinerSpecularSpike) {
      fLinerSpecularSpike = new TabulatedFunction;
      GetStationData(*fLinerSpecularSpike, "specularSpike", "liner", "liner specular spike constant");
    }
    return *fLinerSpecularSpike;
  }


  double
  Station::GetLinerSigmaAlpha()
    const
  {
    if (!fLinerSigmaAlpha.IsValid()) {
      double sa;
      GetStationData(sa, "sigmaAlpha", "liner", "liner sigma_alpha");
      fLinerSigmaAlpha = sa;
    }
    return fLinerSigmaAlpha.Get();
  }


  bool
  Station::IsInAcquisition()
    const
  {
    if (!fIsInAcquisition.IsValid()) {
      int is;
      GetStationData(is, "station_in_acquisition", "is_in_acquisition", "station in acquisition flag");
      fIsInAcquisition = is;
    }
    return fIsInAcquisition.Get();
  }


  int
  Station::GetGroupId()
    const
  {
    if (!fGroupId.IsValid()) {
      GetStationData(fGroupId.Get(), "groupId", "stationList", "group id");
      fGroupId.SetValid();
    }
    return fGroupId.Get();
  }


  const TimeRange&
  Station::GetCommissionTimeRange()
    const
  {
    if (!fCommissionTimeRange.IsValid()) {
      string timeString;
      GetStationData(timeString, "commission", "stationList", "commission time");
      const auto start = boost::lexical_cast<UTCDateTime>(timeString).GetTimeStamp();
      GetStationData(timeString, "decommission", "stationList", "commission time");
      const auto stop = boost::lexical_cast<UTCDateTime>(timeString).GetTimeStamp();
      fCommissionTimeRange = TimeRange(start, stop);
    }
    return fCommissionTimeRange.Get();
  }


  const TimeStamp&
  Station::GetUUBCommissionTime()
    const
  {
    if (!fUUBCommissionTime.IsValid()) {
      string timeString;
      GetStationData(timeString, "commissionUUB", "stationList", "uub commissioning time");
      fUUBCommissionTime = boost::lexical_cast<UTCDateTime>(timeString).GetTimeStamp();
    }
    return fUUBCommissionTime.Get();
  }


  Point  // <-- we don't want to return a const reference here due to the creep bug!
  Station::GetPosition()
    const
  {
    if (!fPosition.IsValid()) {

      double northing = 0;
      GetStationData(northing, "northing", "stationList", "northing");
      double easting = 0;
      GetStationData(easting, "easting", "stationList", "easting");
      double altitude = 0;
      GetStationData(altitude, "altitude", "stationList", "altitude");

      string ellipsoid;
      GetStationData(ellipsoid, "ellipsoid", "stationList", "ellipsoid for UTM conversion");
      fEllipsoid = ReferenceEllipsoid::GetEllipsoidIDFromString(ellipsoid);
      string band;
      GetStationData(band, "band", "stationList", "band for UTM conversion");
      int zone = 0;
      GetStationData(zone, "zone", "stationList", "zone for UTM conversion");

      const UTMPoint utmPosition(northing, easting, altitude, zone, band[0], fEllipsoid);
      const CoordinateSystemPtr cs = Detector::GetInstance().GetSiteCoordinateSystem();
      fPosition = utmPosition.GetPoint(cs);

    }

    return fPosition.Get();
  }


  const vector<int>&
  Station::SetPartnerIds()
    const
  {
    if (!fPartnerIds.IsValid()) {
      fPartnerIds.Get().clear();
      GetStationData(fPartnerIds.Get(), "stationGroup", "stationList", "get partners");
      fPartnerIds.SetValid();
    }
    return fPartnerIds.Get();
  }


  int
  Station::GetNumberOfPartners()
    const
  {
    return SetPartnerIds().size();
  }


  bool
  Station::IsInGrid(const SDetectorConstants::GridIndex index)
    const
  {
    if (!fInGrid.IsValid()) {
      InGridType& inGrid = fInGrid.Get();
      inGrid.reset();
      vector<bool> data;
      GetStationData(data, "inGrid", "stationList", "get is-in-grid");
      for (size_t i = 0, nd = data.size(), ni = inGrid.size(); i < nd && i < ni; ++i)
        inGrid[i] = data[i];
      fInGrid.SetValid();
    }

    const InGridType& inGrid = fInGrid.Get();
    if (index == SDetectorConstants::eAny)
      return inGrid.any();
    else
      return inGrid[index];
  }


  const vector<double>&
  Station::GetAxes()
    const
  {
    if (!fAxes.IsValid()) {
      fAxes.Get().clear();
      double axis = 0;
      GetStationData(axis, "axis1", "stationList", "get axis1");
      fAxes.Get().push_back(axis);
      GetStationData(axis, "axis2", "stationList", "get axis2");
      fAxes.Get().push_back(axis);
      fAxes.SetValid();
    }
    return fAxes.Get();
  }


  const Station::StationIdCollection&
  Station::GetCrown(const int i)
    const
  {
    const auto levelIt = fCrowns.find(i);

    if (levelIt != fCrowns.end())

      return levelIt->second;

    else {

      auto indexMap = fStationIndexMap;
      indexMap["crown"] = boost::lexical_cast<string>(i);

      const auto& manager = Detector::GetInstance().GetSManagerRegister();

      StationIdCollection crown;
      const auto status = manager.GetData(crown, "crown", "stationList", indexMap);

      if (status == VManager::eNotFound) {
        ostringstream err;
        err << "Crown level " << i << " "
               "for station " << fId << " "
               "has not been found by any of the available managers.";
        NotFoundAndThrow(err.str());
      }

      const auto where = fCrowns.insert(make_pair(i, crown));
      return where.first->second;

    }
  }


  void
  Station::NotFoundAndThrow(const string& msg)
    const
  {
    ostringstream err;
    err << "Did not find requested component: '" << msg
        << "' for station = " << fStationIndexMap["stationId"];
    throw NonExistentComponentException(err.str());
  }


  utl::CoordinateSystemPtr
  Station::GetLocalCoordinateSystem()
    const
  {
    if (!fReferenceSystem)
      fReferenceSystem = fwk::LocalCoordinateSystem::Create(GetPosition());

    return fReferenceSystem;
  }


  void
  Station::Update()
    const
  {
    // Clean all the variables that might change with time and delete whatever was allocated in the Get*'s

    fIsInAcquisition.SetValid(false);
    fPartnerIds.SetValid(false);

    if (fDense) {
      fPosition.SetValid(false);
      fReferenceSystem = CoordinateSystemPtr();  // clear the CS ptr so positions can be Updated correctly
    } else {
      const auto& dTime = Detector::GetInstance().GetTime();
      const bool isUUB = (GetUUBCommissionTime() < dTime);
      if (isUUB != fIsUUB)
        UpdateElectronics(isUUB);
    }

    if (fScintillator)
      fScintillator->Update();

    for (const auto p : fPMTVector)
      p->Update();

    fCrowns.clear();
  }


  void
  Station::UpdateElectronics(const bool isuub)
    const
  {
    // Do nothing if electronics is already in requested state
    if (fIsUUB == isuub)
      return;

    // Clearing variables which are different between hardware versions
    fFADCTraceLength.SetValid(false);
    fFADCBinSize.SetValid(false);
    fTimingUncertainty.SetValid(false);
    fLatchBin.SetValid(false);
    fSaturationValue.SetValid(false);
    fElectronicsImpulseResponse = nullptr;

    /* Because fIsUUB is passed to PMTs during
       construction and is necessary for PMTs to retrieve
       data through the managers, all PMTs are deleted
       and recreated here, with the new fIsUUB status
       being passed as an argument of the PMT constructor.
       This process could be replaced with a private
       SetIsUUB member function of the PMT class
       that is visible only to the station through
       friendship. Then again, the electronics-related
       properties of the PMT need to be updated at the same
       time anyway. Most of the same reasoning applies for
       fScintillator. */

    /* Store whether small PMT was present to know whether
       to recreate it with new electronics properties below. */
    const bool hadSmallPMT = HasSmallPMT();

    for (const auto p : fPMTVector)
      delete p;
    fPMTVector.clear();

    fIsUUB = isuub;
    fStationIndexMap["isUUB"] = boost::lexical_cast<string>(fIsUUB);

    /* The 3 PMTs of type eWaterCherenkovLarge must always be re-created
       if there is a change in electronics. See comment above for
       reasoning. */
    for (unsigned int i = kFirstWaterCherenkovLargePMTId;
         i <= kLastWaterCherenkovLargePMTId; ++i) {
      const auto pmt = new PMT(i, fId, fIsUUB, PMTConstants::eWaterCherenkovLarge);
      fPMTVector.push_back(pmt);
    }

    if (hadSmallPMT)
      MakeSmallPMT();

    if (fScintillator) {
      /* Must be set to 0 for re-creation with MakeScintillator() with correct
         This could be avoided with implementation inside a Scintillator Update
         method. */
      fScintillator = nullptr;
      MakeScintillator();
    }
  }


  void
  Station::MakeScintillator()
    const
  {
    // Do nothing if scintillator already exists
    if (fScintillator)
      return;

    fScintillator = new Scintillator(fId, fIsUUB);
    const auto pmt = new PMT(kScintillatorPMTId, fId, fIsUUB, sdet::PMTConstants::eScintillator);
    fPMTVector.push_back(pmt);
  }


  void
  Station::RemoveScintillator()
    const
  {
    // Do nothing if scintillator doesn't exist
    if (!fScintillator)
      return;
    fScintillator = nullptr;

    for (auto pmtIt = fPMTVector.begin(); pmtIt != fPMTVector.end(); ) {
      const auto pmtPtr = *pmtIt;
      if (pmtPtr->GetType() == sdet::PMTConstants::eScintillator) {
        delete pmtPtr;
        fPMTVector.erase(pmtIt);
        continue;
      }
      ++pmtIt;
    }
  }


  void
  Station::MakeSmallPMT()
    const
  {
    if (HasSmallPMT())
      return;

    const auto pmt = new PMT(kWaterCherenkovSmallPMTId, fId, fIsUUB, sdet::PMTConstants::eWaterCherenkovSmall);
    fPMTVector.push_back(pmt);
  }


  void
  Station::RemoveSmallPMT()
    const
  {
    // Do nothing if small PMT doesn't exist
    if (!HasSmallPMT())
      return;

    for (auto pmtIt = fPMTVector.begin(); pmtIt != fPMTVector.end(); ) {
      const auto pmtPtr = *pmtIt;
      if (pmtPtr->GetType() == sdet::PMTConstants::eScintillator) {
        delete pmtPtr;
        fPMTVector.erase(pmtIt);
        continue;
      }
      ++pmtIt;
    }
  }


  unsigned int
  Station::GetNPMTs(const sdet::PMTConstants::PMTType type)
    const
  {
    if (type == sdet::PMTConstants::eAnyType)
      return fPMTVector.size();

    unsigned int n = 0;
    for (const auto p : fPMTVector)
      if (p->GetType() == type)
        ++n;
    return n;
  }


  void
  Station::InitializeShape(vector<double>& histox, const int nBins, const double dt)
  {
    histox.clear();
    histox.reserve(nBins + 1);
    for (int i = 0; i <= nBins; ++i)
      histox.push_back(i * dt);
  }


  const vector<double>&
  Station::GetMuonShapeHistogramBinning(const int calibrationVersion)
    const
  {
    if (IsUUB()) {
      // different number of bins for UUBs
      static vector<double> muShapeUUB;
      if (muShapeUUB.empty())
        InitializeShape(muShapeUUB, 70, fFADCBinSize.Get());
      return muShapeUUB;
    }

    // UB

    if (calibrationVersion == 13) {
      // slightly different for LS calibration version 13
      static vector<double> muShape13;
      if (muShape13.empty())
        InitializeShape(muShape13, 19, fFADCBinSize.Get());
      return muShape13;
    }

    static vector<double> muShape;
    if (muShape.empty())
      InitializeShape(muShape, 20, fFADCBinSize.Get());
    return muShape;
  }


  bool
  Station::HasSmallPMT()
    const
  {
    for (const auto p : fPMTVector)
      if (p->GetType() == sdet::PMTConstants::eWaterCherenkovSmall)
        return true;
    return false;
  }


  bool
  Station::IsInsideStation(const utl::Point& pos)
    const
  {
    const auto& cs = GetLocalCoordinateSystem();
    const double stationRadius = GetRadius();
    const double stationHeight = GetHeight();
    return (pos.GetRho(cs) < stationRadius && pos.GetZ(cs) < stationHeight);
  }


  bool
  Station::IsHit(const Point& particlePoint, const Vector& particleDirection, Point& where)
    const
  {
    static const double farAway2 = Sqr(1000*kilometer);
    const Line particleLine(particlePoint, particleDirection);

    double minDistance2 = farAway2;

    /* top and bottom planes correspond to top and bottom of water volume
       (not top and bottom of tank including casing). */
    const double radius2 = Sqr(GetRadius());
    const auto& cs = GetLocalCoordinateSystem();
    Plane bottomPlane(Point(0, 0, GetThickness(), cs), Vector(0,0,1, cs));
    Plane topPlane(Point(0, 0, GetHeight() + GetThickness(), cs), Vector(0,0,1, cs));

    // try top
    where = Intersection(topPlane, particleLine);
    if ((where - topPlane.GetAnchor()).GetMag2() < radius2) {
      const double dist2 = (particlePoint - where).GetMag2();
      if (dist2 < minDistance2)
        minDistance2 = dist2;
    }

    // try bottom
    where = Intersection(bottomPlane, particleLine);
    if ((where - bottomPlane.GetAnchor()).GetMag2() < radius2) {
      const double dist2 = (particlePoint - where).GetMag2();
      if (dist2 < minDistance2)
        minDistance2 = dist2;
    }

    // try side
    const double ax = particleDirection.GetX(cs);
    const double ay = particleDirection.GetY(cs);
    const double px = particlePoint.GetX(cs);
    const double py = particlePoint.GetY(cs);
    const double a = Sqr(ax) + Sqr(ay);
    const double b = 2*(ax*px + ay*py);
    const double c = Sqr(px) + Sqr(py) - radius2;
    const double disc = Sqr(b) - 4*a*c;
    if (disc < 0)
      return false;
    const double ta = 2*a;
    const double s0 = -b/ta;
    const double delta = fabs(sqrt(disc)/ta);
    const double s = s0 - delta;  // this is the solution with shorter distance

    if (s > 0) {
      where = particleLine.Propagate(s);
      if (Sqr(s) < minDistance2 &&
          Distance(bottomPlane, where) >= 0 && Distance(topPlane, where) <= 0)
        minDistance2 = Sqr(s);
    }

    if (Test<Not<CloseTo>>(minDistance2, farAway2)) {
      where = particleLine.Propagate(sqrt(minDistance2));
      return true;
    }

    return false;
  }

}