VTankResponse.cc

Go to the documentation of this file.

#include <tls/VTankResponse.h>
#include <utl/PoissonDistribution.h>
#include <utl/Math.h>
#include <utl/AugerException.h>

using namespace std;
using namespace utl;
using namespace tls;

typedef unsigned short ushort;

static const double precision = 1e-10;
static const ushort kmax = 1000;

double
VTankResponse::PoissonConvolvedPDF(const double signal,
                                   const double theta,
                                   const double r,
                                   const double muons)
  const
{
  double psum = 0;
  ushort k = 0;

  ulong nMuDown = ulong(muons);
  ulong nMuUp   = nMuDown+1;
  for (double p = 0; k < kmax && (k < 5 || p >= precision*psum); ++k)
  {
    p = PoissonPDF(nMuUp, muons) * PDF(signal, theta, r, nMuUp);
    ++nMuUp;

    if (nMuDown > 0)
    {
      p += PoissonPDF(nMuDown, muons) * PDF(signal, theta, r, nMuDown);
      --nMuDown;
    }

    psum += p;
  }

  // ~if (k == kmax)
  // ~{
    // ~ostringstream msg;
    // ~msg << "PoissonConvolvedPDF: iteration overflow, signal " << signal << " muons " << muons;
    // ~WARNING(msg);
  // ~}

  return psum;
}

double
VTankResponse::PoissonConvolvedCDF(const double sThreshold,
                                   const double theta,
                                   const double r,
                                   const double muons,
                                   const bool complement)
  const
{
  double psum = 0.0;
  ushort k = 0;

  if (complement) // saturation case
  {
    ulong nMuDown = ulong(muons);
    ulong nMuUp   = nMuDown+1;

    for (double p = 0; k < kmax && (k < 5 || p >= precision*psum); ++k)
    {
      p = PoissonPDF(nMuUp, muons) *
          (1.0 - CDF(sThreshold, theta, r, nMuUp));
      ++nMuUp;

      if (nMuDown > 0)
      {
        p += PoissonPDF(nMuDown, muons) *
             (1.0 - CDF(sThreshold, theta, r, nMuDown));
        --nMuDown;
      }

      psum += p;
    }
  }
  else // non-triggering case
  {
    // no muons: CDF(sThreshold) = 1 for sThreshold > 0
    psum = PoissonPDF(0, muons);
    k = 1;
    for (double p = 0; k < kmax && (k < 5 || p >= precision*psum); ++k)
    {
      p = PoissonPDF(k, muons) * CDF(sThreshold, theta, r, k);
      psum += p;
    }
  }

  // ~if (k == kmax)
  // ~{
    // ~ostringstream msg;
    // ~msg << "PoissonConvolvedCDF: iteration overflow, sThr " << sThreshold << " muons " << muons;
    // ~WARNING(msg);
  // ~}

  return psum;
}

void
VTankResponse::PoissonConvolvedMeanAndStDev(double& mean,
                                            double& stDev,
                                            const double theta,
                                            const double r,
                                            const double muons)
  const
{
  double cmom1 = 0;
  double cmom2 = 0;
  ushort k = 0;

  ulong nMuDown = ulong(muons);
  ulong nMuUp   = nMuDown+1;
  for (double mom1 = 0; k < kmax && (k < 5 || mom1 >= precision*cmom1); ++k)
  {
    double mom2 = 0;
    const double p = PoissonPDF(nMuUp,muons);
    const double mean = Mean(theta,r,nMuUp);
    const double sigma = StDev(theta,r,nMuUp);
    mom1 = p * mean;
    mom2 = p * (sigma*sigma + mean*mean);
    ++nMuUp;

    if (nMuDown > 0) {
      const double p = PoissonPDF(nMuDown,muons);
      const double mean = Mean(theta,r,nMuDown);
      const double sigma = StDev(theta,r,nMuDown);
      mom1 += p * mean;
      mom2 += p * (sigma*sigma + mean*mean);
      --nMuDown;
    }

    cmom1 += mom1;
    cmom2 += mom2;
  }

  // ~if (k == kmax)
  // ~{
    // ~ostringstream msg;
    // ~msg << "PoissonConvolvedMeanAndStDev: iteration overflow";
    // ~WARNING(msg);
  // ~}

  mean = cmom1;
  stDev = sqrt(cmom2 - cmom1*cmom1);
}