CachedShowerRegeneratorOG/StationPositionMatrix.cc

Go to the documentation of this file.

#include <utl/Accumulator.h>
#include "StationPositionMatrix.h"

using namespace std;
using namespace utl;


namespace CachedShowerRegeneratorOG {

  void
  StationInfo::Dump()
    const
  {
    std::cerr
      << "ds " << fDetStation->GetId() << " es " << fEvtStation->GetId() << " "
         "phi " << fPhi1 << ' ' << fPhi2 << " "
         "r " << fR << " "
         "lnSqrR " << fLnSqrR1 << ' ' << fLnSqrR2 << " "
         "a " << fResamplingArea << "\n";
  }


  void
  StationPositionMatrix::Clear()
  {
    // ensure negative index is returned for all possible R
    fLnSqrR1 = -1000;
    fLnSqrRStep = -1;
    fStations.clear();
    fStationMatrix.clear();
    fEmptyStationInfoPtrList.clear();
  }


  void
  StationPositionMatrix::CreateMatrix(const bool useSpatialStationMatrix)
  {
    const double tpi = 2*utl::kPi;
    if (!useSpatialStationMatrix || fStations.empty()) {
      Clear();
      // "empty" list of stations will be returned
      fEmptyStationInfoPtrList.clear();
      for (const auto& s : fStations)
        fEmptyStationInfoPtrList.push_back(&s);
      return;
    }
    // find range
    Accumulator::Min<double> lnSqrR1(fStations[0].GetLnSqrR1());
    Accumulator::Max<double> lnSqrR2(fStations[0].GetLnSqrR2());
    Accumulator::Average deltaPhi;
    Accumulator::Average deltaLnSqrR;
    for (const auto& s : fStations) {
      lnSqrR1(s.GetLnSqrR1());
      lnSqrR2(s.GetLnSqrR2());
      deltaLnSqrR(s.GetLnSqrR2() - s.GetLnSqrR1());
      const double dphi = s.GetPhi2() - s.GetPhi1();
      deltaPhi((dphi >= 0) ? dphi : tpi-dphi);
    }
    fNPhi = ceil(tpi / deltaPhi.GetAverage() * fPhiGranularity);
    fPhiStep = tpi / fNPhi;
    const double dLnSqrR = lnSqrR2.GetMax() - lnSqrR1.GetMin();
    fNR = ceil(dLnSqrR / deltaLnSqrR.GetAverage() * fLnSqrRGranularity);
    fLnSqrR1 = lnSqrR1.GetMin();
    fLnSqrRStep = dLnSqrR / fNR;

    Resize(fNPhi, fNR);

    // fill matrix with station info pointers
    for (const auto& s : fStations) {
      const int iPhi1 = GetPhiIndex(s.GetPhi1());
      const int iPhi2 = GetPhiIndex(s.GetPhi2());
      const int iR1 = GetLnSqrRIndex(s.GetLnSqrR1());
      int iR2 = GetLnSqrRIndex(s.GetLnSqrR2());
      if (iR2 == fNR)
        iR2 = fNR - 1;
      if (iPhi1 <= iPhi2) {
        for (int i = iPhi1; i <= iPhi2; ++i)
          for (int j = iR1; j <= iR2; ++j)
          fStationMatrix[i][j].push_back(&s);
      } else {
        for (int j = iR1; j <= iR2; ++j) {
          for (int i = 0; i <= iPhi2; ++i)
            fStationMatrix[i][j].push_back(&s);
          for (int i = iPhi1; i < fNPhi; ++i)
            fStationMatrix[i][j].push_back(&s);
        }
      }
    }
  }


  void
  StationPositionMatrix::DumpStats()
    const
  {
    Accumulator::MinMax<int> minmax(fStationMatrix[0][0].size());
    Accumulator::Average size;
    for (int i = 0; i < fNPhi; ++i) {
      for (int j = 0; j < fNR; ++j) {
        const int s = fStationMatrix[i][j].size();
        minmax(s);
        size(s);
      }
    }
    cerr
      << "\n"
         "nPhi = " << fNPhi << "\n"
         "lnSqrR = " << fLnSqrR1 << " : " << fLnSqrR1 + fNR * fLnSqrRStep << ", " << fNR << "\n"
         "matrix = " << minmax.GetMin() << ", " << size.GetAverage() << ", "
      << minmax.GetMax() << '\n';
    for (const auto& s : fStations)
      s.Dump();
  }


  void
  StationPositionMatrix::Resize(const int nPhi, const int nR)
  {
    fStationMatrix.clear();
    fStationMatrix.resize(nPhi);
    for (auto& s : fStationMatrix)
      s.resize(nR);
  }

}