ShowerSizeFunction.cc

Go to the documentation of this file.

#include "ShowerSizeFunction.h"
#include "SdHorizontalReconstruction.h"
#include "Utilities.h"

#include <utl/AugerUnits.h>
#include <utl/Point.h>
#include <utl/EquationSolver.h>
#include <utl/Math.h>
#include <utl/AugerException.h>

#include <tls/MuonProfile.h>
#include <tls/VTankResponse.h>
#include <tls/EMComponent.h>

#include <limits>

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

ShowerSizeFunction::ShowerSizeFunction(const SdHorizontalReconstruction& config,
                                           const StationList& list,
                                           const SilentStationList& slist,
                                           const AxisData& ad,
                                           const ExternalGeometryData& gd)
  : fType(eFull), fConfig(config), fList(list), fSList(slist), fExt(gd)
{
  fUseAxisCovariance = 0;
  if (ad.fCov[eTheta] > 1e-8 && ad.fCov[ePhi] > 1e-8)
  {
    fUseAxisCovariance = 2;
    if (ad.fCov[eDistance] > 1e-8)
      fUseAxisCovariance = 3;
  }

  for (size_t i = 0; i < fUseAxisCovariance; ++i)
  {
    fAxisPar[i] = ad.fPar[i];
    for (size_t j = 0; j < fUseAxisCovariance; ++j)
      fAxisInvCov[i][j] = ad.fCov(i, j);
  }

  if (fUseAxisCovariance)
  {
    try {
      InvertMatrix(fUseAxisCovariance, fAxisInvCov);
      for (size_t i = 0; i < fUseAxisCovariance; ++i)
        if (fAxisInvCov[i][i] <= 0) throw utl::DoesNotComputeException();
    }
    catch(const utl::DoesNotComputeException& e)
    {
      fUseAxisCovariance = 0;
    }
  }
}


double
ShowerSizeFunction::operator()(const std::vector<double>& pars)
  const
{
  for (size_t i = 0; i < 6; ++i)
    if (std::isnan(pars[i]))
      return numeric_limits<double>::max();

  const double n19 = pars[0];
  const double coreX = pars[1]*km;
  const double coreY = pars[2]*km;
  const double theta0 = pars[3];
  const double phi0 = pars[4];
  const double distance0 = pars[5]*km;

  Point core, origin;
  fConfig.GetShowerAxis(core, origin, coreX, coreY, theta0, phi0, distance0);

  double result = 0;

  const Vector axis = origin - core;

  const double theta = axis.GetTheta(fConfig.fBaryCS);
  const double phi = axis.GetPhi(fConfig.fBaryCS);

  // station did trigger
  for (StationList::const_iterator
         sIt = fList.begin(),
         sEnd = fList.end();
       sIt != sEnd; ++sIt)
  {
    if (sIt->fRejected) continue;

    if (sIt->fSaturated)
    {
      if (sIt->fRecoveryErr > 0)
      { // recovered station: assume Gaussian fluctuations and include recovery error
        double mean, sigma;
        Predict(mean, sigma, sIt->fPos, origin, n19, core, sIt->fRecoveryErr);
        const double z = (sIt->fSignal - mean)/sigma;
        result += 0.5*z*z;
      }
      else
      { // non-recovered saturated station
        if (fType == eFull)
        {
          // measured signal is lower limit to unknown actual signal
          double xProj, yProj, rho, sTheta;
          fConfig.LocalCoordinates(xProj, yProj, rho, sTheta, sIt->fPos, origin, core);
          const double nMu =
            n19 * fConfig.fMuonProfile->NMuon(xProj, yProj, theta, phi, 10*EeV);
          const double ratioSEmSMu =
            fConfig.fEMComponent->SignalRatio(xProj, yProj, n19, theta, phi);

          const double signalMu = sIt->fSignal/(1.0+ratioSEmSMu);
          const double prob =
            fConfig.fTankResponse->PoissonConvolvedCDF(signalMu, sTheta, rho, nMu, true);
          result -= log(prob);
        }
        else
        {
          // use saturated signal as measured signal, assign huge error
          double mean, sigma;
          Predict(mean, sigma, sIt->fPos, origin, n19, core);
          const double z = (sIt->fSignal - mean)/mean;
          result += 0.5*z*z;
        }
      }
    } // is saturated
    else
    { // non-saturated regular station with signal
      if (fType == eFull)
      {
        if (fConfig.fDropStationsBelowThreshold &&
            (sIt->fSignal < fConfig.fSilentSignalThreshold))
        {
          // treat stations below assumed trigger threshold as non-triggering stations
          double xProj, yProj, rho, sTheta;
          fConfig.LocalCoordinates(xProj, yProj, rho, sTheta, sIt->fPos, origin, core);
          const double nMu =
            n19 * fConfig.fMuonProfile->NMuon(xProj, yProj, theta, phi, 10*EeV);
          const double ratioSEmSMu =
            fConfig.fEMComponent->SignalRatio(xProj, yProj, n19, theta, phi);
          const double sThrMu = fConfig.fSilentSignalThreshold/(1.0+ratioSEmSMu);
          const double prob =
            fConfig.fTankResponse->PoissonConvolvedCDF(sThrMu, sTheta, rho, nMu, false);
          result -= log(prob);
        }
        else
        {
          // probability of observing measured signal
          double xProj, yProj, rho, sTheta;
          fConfig.LocalCoordinates(xProj, yProj, rho, sTheta, sIt->fPos, origin, core);
          const double nMu =
            n19 * fConfig.fMuonProfile->NMuon(xProj, yProj, theta, phi, 10*EeV);
          const double ratioSEmSMu =
            fConfig.fEMComponent->SignalRatio(xProj, yProj, n19, theta, phi);

          const double signalMu = sIt->fSignal/(1.0+ratioSEmSMu);
          const double prob =
            fConfig.fTankResponse->PoissonConvolvedPDF(signalMu, sTheta, rho, nMu);
          result -= log(prob);
        }
      }
      else
      {
        // assume Gaussian fluctuations
        double mean, sigma;
        Predict(mean, sigma, sIt->fPos, origin, n19, core);
        const double z = (sIt->fSignal - mean)/sigma;
        result += 0.5*z*z;
      }
    } // regular station
  } // loop over candidate stations

  // station did not trigger
  for (SilentStationList::const_iterator
         sIt = fSList.begin(),
         sEnd = fSList.end();
       sIt != sEnd; ++sIt)
  {
    if (fType == eFull)
    {
      // probability of nMu to produce a signal below the T2 trigger threshold
      double xProj, yProj, rho, sTheta;
      fConfig.LocalCoordinates(xProj, yProj, rho, sTheta, sIt->fPos, origin, core);
      const double nMu =
        n19 * fConfig.fMuonProfile->NMuon(xProj, yProj, theta, phi, 10*EeV);
      const double ratioSEmSMu =
        fConfig.fEMComponent->SignalRatio(xProj, yProj, n19, theta, phi);
      const double sThrMu = fConfig.fSilentSignalThreshold/(1.0+ratioSEmSMu);
      // probability of nMu to produce a signal below the T2 trigger threshold
      const double prob =
        fConfig.fTankResponse->PoissonConvolvedCDF(sThrMu, sTheta, rho, nMu, false);
      result -= log(prob);
    }
    else
    {
      // assume Gaussian fluctuations, treat as if measured signal was zero
      double mean, sigma;
      Predict(mean, sigma, sIt->fPos, origin, n19, core);
      const double z = mean/sigma;
      result += 0.5*z*z;
    }
  } // loop over zero signal stations

  // axis may vary within its uncertainty
  if (fUseAxisCovariance)
  {
    const double delta[3] = {
      theta0 - fAxisPar[0],
      phi0 - fAxisPar[1],
      distance0 - fAxisPar[2]
    };

    for (size_t i = 0; i < fUseAxisCovariance; ++i)
      for (size_t j = 0; j < fUseAxisCovariance; ++j)
        result += 0.5 * delta[i] * fAxisInvCov[i][j] * delta[j];
  }

  // core and axis may vary within its external uncertainty
  {
    double delta[4];
    delta[eCoreXExt] = coreX - fExt.fPar[eCoreXExt];
    delta[eCoreYExt] = coreY - fExt.fPar[eCoreYExt];
    delta[eThetaExt] = theta0 - fExt.fPar[eThetaExt];
    delta[ePhiExt] = phi0 - fExt.fPar[ePhiExt];

    for (size_t i = 0; i < 4; ++i)
      for (size_t j = 0; j < 4; ++j)
        result += 0.5 * delta[i] * fExt.fInvCov[i][j] * delta[j];
  }

  return result;
}


void
ShowerSizeFunction::Predict(double& meanSignal, double& sigmaSignal,
                              const Point& station,
                              const Point& origin,
                              const double n19,
                              const Point& core,
                              const double recoveryErr)
  const
{
  double xProj, yProj, rho, sTheta;
  fConfig.LocalCoordinates(xProj, yProj, rho, sTheta, station, origin, core);

  const Vector axis = origin - core;

  const double theta = axis.GetTheta(fConfig.fBaryCS);
  const double phi = axis.GetPhi(fConfig.fBaryCS);

  const double nMu =
    n19 * fConfig.fMuonProfile->NMuon(xProj, yProj, theta, phi, 10*EeV);

  const double ratioSEmSMu =
    fConfig.fEMComponent->SignalRatio(xProj, yProj, n19, theta, phi);

  fConfig.fTankResponse->PoissonConvolvedMeanAndStDev(meanSignal, sigmaSignal,
                                                      sTheta, rho, nMu);

  meanSignal  *= (1 + ratioSEmSMu);
  sigmaSignal *= (1 + ratioSEmSMu);
  if (recoveryErr > 0)
    sigmaSignal = sqrt(Sqr(sigmaSignal) + Sqr(recoveryErr));
}