Tools/InclinedShowers/TankResponse/Simple/TankResponse.cc

Go to the documentation of this file.

#include "TankResponse.h"

#include <tls/TankResponseUtilities.h>

#include <utl/Reader.h>
#include <utl/Math.h>
#include <utl/NormalDistribution.h>

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


TankResponse::TankResponse(const Branch branch) :
  fVerticalTankTrackLength(TankMeanTrackLength(0))
{
  branch.GetChild("ParameterMean").GetData(fParMean);
  branch.GetChild("ParameterSigma").GetData(fParSigma);
}


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


double
TankResponse::CDF(const double threshold,
                                    const double theta,
                                    const double /* r */,
                                    const ulong muons)
  const
{
  double m, s;
  CalculateAverageAndSigma(m, s, theta, muons);
  return NormalCDF(threshold, m, s);
}


double
TankResponse::PDF(const double signal,
                          const double theta,
                          const double /* r */,
                          const ulong muons)
  const
{
  double m, s;
  CalculateAverageAndSigma(m, s, theta, muons);
  return NormalPDF(signal, m, s);
}


double
TankResponse::Mean(const double theta,
                   const double /* r */,
                   const ulong muons)
  const
{
  double m, s;
  CalculateAverageAndSigma(m, s, theta, muons);
  return m;
}


double
TankResponse::StDev(const double theta,
                    const double /* r */,
                    const ulong muons)
  const
{
  double m, s;
  CalculateAverageAndSigma(m, s, theta, muons);
  return s;
}


void
TankResponse::CalculateAverageAndSigma(double& mean, double& sigma,
                                       const double theta, const ulong muons)
  const
{
  // signal of vertical muon := 1 VEM
  // signal of inclined muon = signalVertical * averageTrackLength(theta)/averageTrackLength(0)
  mean = TankMeanTrackLength(theta)/fVerticalTankTrackLength * muons;
  const double nmuWithCut = muons < 10 ? muons : 10;
  const double meanCorrection =
    fParMean[0] + theta*(fParMean[1] + theta*(fParMean[3] + theta*fParMean[5])) +
    nmuWithCut*(fParMean[2] + nmuWithCut*fParMean[4]);
  mean *= (1 + meanCorrection);

  const double th = theta/degree;
  const double trackMoment2 =
    fParSigma[0] + th*(fParSigma[1] + th*(fParSigma[2] + th*(fParSigma[3] + th*fParSigma[4])));

  const double trackSigma2 = trackMoment2 - pow(TankMeanTrackLength(theta), 2);
  const double energyMu = 3; // GeV
  const double b = (2.33726e-01 + 2.26614e-02*log(energyMu))*trackSigma2;
  sigma = b*sqrt(muons);
}