CachedShowerRegeneratorOG/StationPositionMatrix.h

Go to the documentation of this file.

#ifndef _CachedShowerRegeneratorOG_StationPositionMatrix_h_
#define _CachedShowerRegeneratorOG_StationPositionMatrix_h_

#include <utl/MathConstants.h>
#include <sdet/Station.h>
#include <sevt/Station.h>


namespace CachedShowerRegeneratorOG {

  class StationInfo {
  public:
    StationInfo(const sdet::Station& dStation, sevt::Station& eStation,
                const double phi, const double dphi,
                const double r, const double lnSqrR1, const double lnSqrR2,
                const double area) :
      fDetStation(&dStation),
      fEvtStation(&eStation),
      fR(r),
      fLnSqrR1(lnSqrR1),
      fLnSqrR2(lnSqrR2),
      fResamplingArea(area)
    {
      SetPhiRange(phi, dphi);
    }

    bool IsInArea(const double phi, const double lnSqrR) const
    { return IsPhiInRange(phi) && IsLnSqrRInRange(lnSqrR); }

    bool
    IsInScaledArea(const double scale, const double phi, const double lnSqrR)
      const
    {
      const double s = std::sqrt(scale);
      return IsPhiInScaledRange(s, phi) && IsLnSqrRInScaledRange(s, lnSqrR);
    }

    const sdet::Station& GetDetStation() const { return *fDetStation; }

    sevt::Station& GetEvtStation() const { return *fEvtStation; }

    double GetResamplingArea() const { return fResamplingArea; }

    double GetPhi1() const { return fPhi1; }

    double GetPhi2() const { return fPhi2; }

    double GetR() const { return fR; }

    double GetLnSqrR1() const { return fLnSqrR1; }

    double GetLnSqrR2() const { return fLnSqrR2; }

    void Dump() const;

  private:
    void
    SetPhiRange(const double phi, const double dphi)
    {
      fPhi1 = phi - dphi;
      fPhi2 = phi + dphi;
      if (fPhi1 < -utl::kPi) {
        if (fPhi2 > utl::kPi) {
          fPhi1 = -utl::kPi;
          fPhi2 = utl::kPi;
        } else
          fPhi1 += 2*utl::kPi;
      } else if (fPhi2 > utl::kPi)
        fPhi2 -= 2*utl::kPi;
    }

    static constexpr bool IsInInterval(const double x, const double a, const double b)
    { return a <= x && x <= b; }

    bool IsPhiInRange(const double phi) const
    { return (fPhi1 < fPhi2) ? IsInInterval(phi, fPhi1, fPhi2) : !IsInInterval(phi, fPhi2, fPhi1); }

    static constexpr double AbsShift(const double eps, const double a, const double b)
    { return 0.5*((1 - eps)*a + (1 + eps)*b); }

    bool
    IsPhiInScaledRange(const double scale, const double phi)
      const
    {
      if (!IsPhiInRange(phi))
        return false;
      double p1 = fPhi1;
      double p2 = fPhi2;
      if (fPhi1 > fPhi2) {
        if (phi < fPhi2)
          p1 -= utl::kTwoPi;
        else
          p2 += utl::kTwoPi;
      }
      return IsInInterval(phi, AbsShift(-scale, p1, p2), AbsShift(scale, p1, p2));
    }

    bool IsLnSqrRInRange(const double lnSqrR) const
    { return IsInInterval(lnSqrR, fLnSqrR1, fLnSqrR2); }

    static
    std::pair<double, double>
    RelShift(const double eps, const double a, const double b)
    {
      const double ea = std::exp(0.5*a);
      const double eb = std::exp(0.5*b);
      const double em = 0.5*(1 - eps);
      const double ep = 0.5*(1 + eps);
      return std::make_pair(2*std::log(ep*ea + em*eb), 2*std::log(em*ea + ep*eb));
    }

    bool
    IsLnSqrRInScaledRange(const double scale, const double lnSqrR)
      const
    {
      if (!IsLnSqrRInRange(lnSqrR))
        return false;
      const auto lnSqrR12 = RelShift(scale, fLnSqrR1, fLnSqrR2);
      return IsInInterval(lnSqrR, lnSqrR12.first, lnSqrR12.second);
    }

    const sdet::Station* const fDetStation;
    sevt::Station* const fEvtStation;
    double fPhi1 = 0;
    double fPhi2 = 0;
    double fR = 0;
    double fLnSqrR1 = 0;
    double fLnSqrR2 = 0;
    double fResamplingArea = 0;
  };


  using namespace std;
  using namespace utl;

  class StationPositionMatrix {
  public:
    typedef vector<StationInfo> StationInfoList;
    typedef vector<const StationInfo*> StationInfoPtrList;

    StationPositionMatrix(const double phiGranularity, const double lnSqrRGranularity) :
      fPhiGranularity(phiGranularity),
      fLnSqrRGranularity(lnSqrRGranularity)
    { }

    void PushBack(const sdet::Station& dStation, sevt::Station& eStation,
                  const double phi, const double dphi,
                  const double r, const double lnSqrR1, const double lnSqrR2,
                  const double area)
    { fStations.emplace_back(dStation, eStation, phi, dphi, r, lnSqrR1, lnSqrR2, area); }

    void CreateMatrix(const bool useSpatialStationMatrix = true);

    const StationInfoList& GetStationList() const
    { return fStations; }

    double GetMinLnSqrR() const { return fLnSqrR1; }

    const StationInfoPtrList&
    GetInjectionCandidates(const double phi, const double lnSqrR)
      const
    {
      const int iR = GetLnSqrRIndex(lnSqrR);
      if (iR < 0 || iR >= fNR)
        return fEmptyStationInfoPtrList;
      const int iPhi = GetPhiIndex(phi);
      return fStationMatrix[iPhi][iR];
    }

    StationInfoPtrList
    GetInjectionStations(const double phi, const double lnSqrR)
      const
    {
      StationInfoPtrList stations;
      for (const auto candidate : GetInjectionCandidates(phi, lnSqrR))
        if (candidate->IsInArea(phi, lnSqrR))
          stations.push_back(candidate);
      return stations;
    }

    void Clear();

    void DumpStats() const;

  private:
    int
    GetPhiIndex(const double phi)
      const
    {
      const int i = int((phi + utl::kPi) / fPhiStep);
      return (i < fNPhi) ? i : fNPhi - 1;
    }

    int GetLnSqrRIndex(const double lnSqrR) const
    { return int((lnSqrR - fLnSqrR1) / fLnSqrRStep); }

    void Resize(const int nPhi, const int nR);

    const double fPhiGranularity;
    const double fLnSqrRGranularity;

    int fNPhi = 0;
    double fPhiStep = 0;
    int fNR = 0;
    double fLnSqrR1 = 0;
    double fLnSqrRStep = 0;

    StationInfoList fStations;

    typedef vector<vector<StationInfoPtrList>> StationMatrix;
    StationMatrix fStationMatrix;

    StationInfoPtrList fEmptyStationInfoPtrList; // can contain all stations in case the matrix is not used
  };

}


#endif