Tools/InclinedShowers/TankResponse/Tabular/TankResponse.cc

Go to the documentation of this file.

#include "TankResponse.h"

#include <utl/AugerException.h>
#include <utl/Reader.h>
#include <utl/Branch.h>
#include <utl/AugerUnits.h>

#include <utl/Math.h>

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

#ifdef _OPENMP
  #include <omp.h>
#endif

#include <boost/lexical_cast.hpp>

using namespace TabularTankResponseNS;
using namespace utl;
using namespace std;
using namespace boost;
using namespace tls;

// not available in some compiler collections
inline double sqrt(const size_t x) { return sqrt(double(x)); }

inline
void
Trafo(size_t& i, double&z, const double x, const size_t n, const double x0, const double dx)
{
  if (std::isnan(x))
  {
    i = 0;
    z = x;
  }
  else if (x <= x0)
  {
    i = 0;
    z = 0;
  }
  else if (x >= x0+dx)
  {
    i = n-1;
    z = 0;
  }
  else
  {
    z = (x-x0)/dx*(n-1);
    i = static_cast<size_t>(z);
    z -= i;
  }
}

struct PdfData {
  double fX;
  double fY;
  vector<double> fPDF;
};

TankResponse&
TankResponse::GetInstance(const Branch branch)
{
  static TankResponse singleton(branch);
  return singleton;
}

TankResponse::TankResponse(const Branch branch)
{
  size_t nmuons;
  branch.GetChild("NMuons").GetData(nmuons);

  string filename;
  branch.GetChild("DataFile").GetData(filename);

  branch.GetChild("PDFSmoothing").GetData(fSmooth);

  const Branch b = Reader(filename).GetTopBranch();

  const double nan = numeric_limits<double>::quiet_NaN();
  fNX = 0;
  fNY = 0;
  double xrange[2] = { nan, nan };
  double yrange[2] = { nan, nan };
  fDS = 0;

  // discover grid properties and read pdfs
  vector<PdfData> pdfs;
  for (Branch bpdf = b.GetFirstChild(); bpdf; bpdf = bpdf.GetNextSibling())
  {
    AttributeMap am = bpdf.GetAttributes();
    const double x = lexical_cast<double>(am["x"])*degree;
    const double y = lexical_cast<double>(am["y"])+log10(meter);
    fDS = lexical_cast<double>(am["ds"]);

    if (std::isnan(xrange[0]) || x < xrange[0] || xrange[1] < x)
    {
      ++fNX;
      if (x < xrange[0] || std::isnan(xrange[0])) xrange[0] = x;
      if (x > xrange[1] || std::isnan(xrange[1])) xrange[1] = x;
    }

    if (std::isnan(yrange[0]) || y < yrange[0] || yrange[1] < y)
    {
      ++fNY;
      if (y < yrange[0] || std::isnan(yrange[0])) yrange[0] = y;
      if (y > yrange[1] || std::isnan(yrange[1])) yrange[1] = y;
    }
    pdfs.push_back(PdfData());
    pdfs.back().fX = x;
    pdfs.back().fY = y;
    bpdf.GetData(pdfs.back().fPDF);
  }

  if (yrange[0]==yrange[1])
  {
    // no y-dependency, set accepted range to something arbitrary,
    // constant extrapolation is done automatically
    yrange[0] = 0;
    yrange[1] = 1;
  }

  fX = xrange[0];
  fDX = xrange[1]-xrange[0];
  fY = yrange[0];
  fDY = yrange[1]-yrange[0];
  fMean.resize(extents[fNX][fNY]);
  fSigma.resize(extents[fNX][fNY]);
  fCdf.resize(extents[fNX][fNY][nmuons]);

  cout << "TabularTankResponse convolving";
  #ifdef _OPENMP
    #pragma omp parallel for
  #endif
  for (size_t i = 0; i < pdfs.size(); ++i)
  {
    cout << "." << flush;
    const vector<double>& pdf = pdfs[i].fPDF;
    const double x = pdfs[i].fX;
    const double y = pdfs[i].fY;

    const size_t ix = round((x-xrange[0])/(xrange[1]-xrange[0])*(fNX-1));
    const size_t iy = round((y-yrange[0])/(yrange[1]-yrange[0])*(fNY-1));

    const size_t n = pdf.size();

    // numerical auto-convolution step:
    // without compiler optimization assignment operators of vector
    // and multi_array are a serious bottlenecks here,
    // so I operate on memory directly
    const double* input = &pdf.front();
    for (size_t m = 0; m < nmuons; ++m)
    {
      fCdf[ix][iy][m].resize(1+(m+1)*n);
      double* data = &fCdf[ix][iy][m].front();
      data[0] = 0;

      if (m == 0)
      {
        double sum = 0;
        double vsum = 0;
        double vvsum = 0;
        /* Expectation of piecewise constant p.d.f.:
         * int dx x f(x) = sum_i f[i] 1/2 (x[i+1]^2-x[i]^2)
         *               = sum_i 1/2 f[i] ds^2 (2i+1)
         *               = sum_i f[i] ds^2 (i+1/2)
         * Second moment analog:
         * int dx x^2 f(x) = sum_i f[i] 1/3 (x[i+1]^3-x[i]^3)
         *                 = sum_i f[i] ds^3 (i^2 + i + 1/3)
         */
        for (size_t i = 0; i < n; ++i)
        {
          data[1+i] = input[i];
          sum += input[i];
          vsum += input[i]*(i+0.5);
          vvsum += input[i]*(i*i+i+1.0/3.0);
        }
        vsum *= fDS; // one ds cancels when dividing by sum term
        vvsum *= fDS*fDS; // one ds cancels when dividing by sum term
        fMean[ix][iy] = vsum/sum;
        fSigma[ix][iy] = sqrt(vvsum/sum - fMean[ix][iy]*fMean[ix][iy]);
      }
      else
      {
        const double* zero = &fCdf[ix][iy][0].front();
        const double* prev = &fCdf[ix][iy][m-1].front();
        for (size_t i = 0; i < (m+1)*n; ++i)
        {
          data[1+i] = 0;
          for (size_t j = i >= n ? i-n+1 : 0; j <= min(i,m*n-1); ++j)
            data[1+i] += zero[1+i-j]*prev[1+j];
        }
      }
    }

    for (size_t m = 0; m < nmuons; ++m)
    {
      // compute CDF from PDF
      double* data = &fCdf[ix][iy][m].front();
      size_t i = 0;
      for (; i < (m+1)*n && data[i+1] > 0; ++i)
        data[1+i] = data[i] + data[1+i];

      // Normalize CDF
      const double mean = fMean[ix][iy]*(m+1);
      const double sigma = fSigma[ix][iy]*sqrt(m+1);
      const double cTrans = NormalCDF((i+0.5)*fDS, mean, sigma);
      for (size_t j = 0; j < i; ++j)
        data[j] *= cTrans / data[i];
      for (size_t j = i; j < (m+1)*n; ++j)
        data[j] = NormalCDF((j+0.5)*fDS, mean, sigma);
    }
  } // loop over pdfs
  cout << endl;
}

double
TankResponse::PDF(const double signal,
                  const double theta,
                  const double r,
                  const ulong nmuons)
  const
{
  // no need to check ranges, constant extrapolation is done automatically
  double result = 0;

  if (nmuons <= fCdf.shape()[2])
  {
    const double h = fSmooth*fDS;
    double sm = signal - h; if (sm < 0) sm = 0;
    const double sp = signal + h;

    const double cm = CDF(sm, theta, r, nmuons);
    const double cp = CDF(sp, theta, r, nmuons);
    result = (cp-cm)/(sp-sm);
  }
  else
  {
    const double m = Mean(theta, r, nmuons);
    const double s = StDev(theta, r, nmuons);
    result = NormalPDF(signal, m, s);
  }

  return result;
}

double
TankResponse::CDF(const double signalmax,
                  const double theta,
                  const double r,
                  const ulong nmuons)
  const
{
  // no need to check ranges, constant extrapolation is done automatically
  if (nmuons == 0)
    return signalmax > 0;

  if (nmuons <= fCdf.shape()[2])
  {
    size_t ix, iy;
    double zx, zy;
    Trafo(ix,zx,theta,fNX,fX,fDX);
    Trafo(iy,zy,log10(r),fNY,fY,fDY);

    const vector<double>& v00 = fCdf[ix              ][iy              ][nmuons-1];
    const vector<double>& v01 = fCdf[ix              ][min(iy+1, fNY-1)][nmuons-1];
    const vector<double>& v10 = fCdf[min(ix+1, fNX-1)][iy              ][nmuons-1];
    const vector<double>& v11 = fCdf[min(ix+1, fNX-1)][min(iy+1, fNY-1)][nmuons-1];

    // interpolation idea from USC group (original in USCInterTankResponse)
    const double z = (signalmax-Mean(theta, r, nmuons))/StDev(theta, r, nmuons);
    const double sqrtn = sqrt(nmuons);
    const double m00 = fMean[ix][iy]*nmuons;
    const double m01 = fMean[ix][min(iy+1, fNY-1)]*nmuons;
    const double m10 = fMean[min(ix+1, fNX-1)][iy]*nmuons;
    const double m11 = fMean[min(ix+1, fNX-1)][min(iy+1, fNY-1)]*nmuons;
    const double s00 = fSigma[ix][iy]*sqrtn;
    const double s01 = fSigma[ix][min(iy+1, fNY-1)]*sqrtn;
    const double s10 = fSigma[min(ix+1, fNX-1)][iy]*sqrtn;
    const double s11 = fSigma[min(ix+1, fNX-1)][min(iy+1, fNY-1)]*sqrtn;
    const double p00 = EvaluateCdf(v00,s00*z+m00, m00, s00);
    const double p01 = EvaluateCdf(v01,s01*z+m01, m01, s01);
    const double p10 = EvaluateCdf(v10,s10*z+m10, m10, s10);
    const double p11 = EvaluateCdf(v11,s11*z+m11, m11, s11);
    return (1.0-zx)*(1.0-zy)*p00
          +(1.0-zx)*     zy *p01
          +     zx *(1.0-zy)*p10
          +     zx *     zy *p11;
  }
  else
  {
    const double m = Mean(theta, r, nmuons);
    const double s = StDev(theta, r, nmuons);
    return NormalCDF(signalmax, m, s);
  }
}

double
TankResponse::Mean(const double theta,
                   const double r,
                   const ulong nmuons)
  const
{
  // no need to check ranges, constant extrapolation is done automatically
  size_t ix, iy;
  double zx, zy;
  Trafo(ix,zx,theta,fNX,fX,fDX);
  Trafo(iy,zy,log10(r),fNY,fY,fDY);

  return ((1.0-zx)*(1.0-zy)*fMean[ix              ][iy              ]
         +(1.0-zx)*     zy *fMean[ix              ][min(iy+1, fNY-1)]
         +     zx *(1.0-zy)*fMean[min(ix+1, fNX-1)][iy              ]
         +     zx *     zy *fMean[min(ix+1, fNX-1)][min(iy+1, fNY-1)])*nmuons;
}

double
TankResponse::StDev(const double theta,
                    const double r,
                    const ulong nmuons)
  const
{
  // no need to check ranges, constant extrapolation is done automatically
  size_t ix, iy;
  double zx, zy;
  Trafo(ix,zx,theta,fNX,fX,fDX);
  Trafo(iy,zy,log10(r),fNY,fY,fDY);

  return ((1.0-zx)*(1.0-zy)*fSigma[ix              ][iy              ]
         +(1.0-zx)*     zy *fSigma[ix              ][min(iy+1, fNY-1)]
         +     zx *(1.0-zy)*fSigma[min(ix+1, fNX-1)][iy              ]
         +     zx *     zy *fSigma[min(ix+1, fNX-1)][min(iy+1, fNY-1)])*sqrt(nmuons);
}