ParametricMuonProfile.cc

Go to the documentation of this file.

#include "ParametricMuonProfile.h"
#include "IndexProcessor.h"

#include <utl/MathConstants.h>
#include <utl/AugerUnits.h>
#include <utl/Branch.h>

#include <cmath>
#include <string>
#include <iostream>

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

void ParametricMuonProfile::ReadModelData(const Branch& modelBranch,
                                      ModelDataType& tensor,
                                      LocalDataType& tensorLocal) const
{
  std::vector<size_t> extent;
  modelBranch.GetChild("matrixExtension").GetData(extent);

  tensor.resize(extent);
  tensorLocal.resize(extent);

  std::vector<double> data;
  modelBranch.GetChild("data").GetData(data);

  IndexProcessor<6> ip(extent);
  size_t kjvmlu[6];
  for (size_t i=0,ii=0; i < tensor.num_elements(); ++i)
  {
    tensor.data()[i] = 0;
    ip.Delinearize(i, kjvmlu);
    // skip coefficients which are multiplied with zero in expansion
    if ( (kjvmlu[1]==0 && kjvmlu[2]==1) || (kjvmlu[4]==0 && kjvmlu[5]==1) )
      continue;
    tensor.data()[i] = data[ii];
    ++ii;
  }
}


ParametricMuonProfile::ParametricMuonProfile(Branch topBranch) :
  fTheta(0.), fPhi(0.), fEnergy(10.*EeV), fBPsi(0.)
{
  topBranch.GetChild("thetaRange").GetData(fThetaRange);
  topBranch.GetChild("radiusMax").GetData(fSqrtRadiusMax);
  fSqrtRadiusMax = sqrt(fSqrtRadiusMax);

  topBranch.GetChild("dModel").GetData(fDMaxData);

  topBranch.GetChild("nMuModel").GetData(fNMuData);

  Branch b;
  std::vector<double> data;

  b = topBranch.GetChild("nMuTankModel");
  b.GetChild("biasCorrectionFactor").GetData(fBiasCorrection[eNMuTank]);
  ReadModelData(b, fModelData[eNMuTank], fLocalData[eNMuTank]);

  b = topBranch.GetChild("cosThetaModel");
  b.GetChild("biasCorrectionFactor").GetData(fBiasCorrection[eCosTheta]);
  ReadModelData(b, fModelData[eCosTheta], fLocalData[eCosTheta]);

  b = topBranch.GetChild("logEnergyModel");
  b.GetChild("biasCorrectionFactor").GetData(fBiasCorrection[eLogEnergy]);
  ReadModelData(b, fModelData[eLogEnergy], fLocalData[eLogEnergy]);
}


double ParametricMuonProfile::D(double theta) const
{
  /* Original analytical function was expanded (approximated) in a series of
     Chebyshev polynomials. Evaluation uses Clenshaw's recurrence formula
     with following recursion for Chebyshev polynomials T[n](z):
       T[n+1](z) = 2 z T[n](z) - T[n-1](z)                                */

  const unsigned int n = fDMaxData.size();
  const double z = 2*theta/(0.5*kPi)-1.0;
  double d = 0.0;
  double dd = 0.0;
  for (unsigned int i=n-1;i!=0;--i)
  {
    const double sv = d;
    d = 2*z*d - dd + fDMaxData[i];
    dd = sv;
  }
  return z*d + 0.5*fDMaxData[0] - dd;
}


double ParametricMuonProfile::TotalMuonNumber(double theta, double energy) const
{
  const double lgE = log10(energy/EeV);
  const double lgD = log10(D(theta));
  return pow(10.,fNMuData[0]+fNMuData[1]*lgE+fNMuData[2]*lgD);
}


double ParametricMuonProfile::EvaluateProfile(ProfileType type, double x, double y,
                                          double theta, double phi, double energy)
{
  CachedSetProfile(theta,phi,energy);

  // translate ground to shower plane coordinates
  double sqr = sqrt( rPerpendicular( x,  y, theta, phi)/meter );
  const double zeta = PsiPerpendicular( x,  y, theta, phi ) - fBPsi + 0.5*kPi;

  double result = 0.0;

  switch(type)
  {
    case eNMuTank:
    if (sqr > fSqrtRadiusMax)
    { // extrapolation with formula ln N_\mu = a + b sqrt{r}
      const double a = EvaluateLocalProfile(eNMuTank, fSqrtRadiusMax, zeta);
      const double b = ( a - EvaluateLocalProfile(eNMuTank, fSqrtRadiusMax*0.95, zeta) )
                       /(0.05*fSqrtRadiusMax);
      result = fNMuonNumber
              *fBiasCorrection[eNMuTank]
              *pow(10., a + b*(sqr-fSqrtRadiusMax));
    }
    else
      result = fNMuonNumber
              *fBiasCorrection[eNMuTank]
              *pow(10., EvaluateLocalProfile(eNMuTank, sqr, zeta));
    break;

    case eCosTheta:
    if (sqr > fSqrtRadiusMax) sqr = fSqrtRadiusMax;
    result = fBiasCorrection[eCosTheta]
            *EvaluateLocalProfile(eCosTheta, sqr, zeta);
    if (result < 0) result = 0.0;
    if (result > 1) result = 1.0;
    result = acos(result);

    break;

    case eLogEnergy:
    if (sqr > fSqrtRadiusMax) sqr = fSqrtRadiusMax;
    result = fBiasCorrection[eLogEnergy]
            *pow(10.,EvaluateLocalProfile(eLogEnergy, sqr, zeta))*GeV;
    break;
  }

  return result;
}


double ParametricMuonProfile::EvaluateLocalProfile(ProfileType type, double sqr, double zeta) const
{
  double result = 0.0;

  for (size_t k=0; k < fLocalData[type].shape()[0]; ++k)
    for (size_t j=0; j < fLocalData[type].shape()[1]; ++j)
      for (size_t v=0; v < 2; ++v)
  {
    // skip coefficients which are multiplied with zero in expansion
    if (j==0 && v==1) continue;

    if (v==0)
      result += fLocalData[type][k][j][v]*pow(sqr, int(k))*cos(j*zeta);
    else
      result += fLocalData[type][k][j][v]*pow(sqr, int(k))*sin(j*zeta);
  }

  return result;
}


void ParametricMuonProfile::BuildLocalModel()
{
  const double dmax = D(90*degree);
  const double dmin = D(0);
  const double zA = 2*(D(fTheta)-dmin)/(dmax-dmin)-1;
  const double zB = 2*(fTheta/degree-60)/28-1;

  for (int iType=0;iType<3;++iType)
    for (size_t k=0; k<fModelData[iType].shape()[0]; ++k)
      for (size_t j=0; j<fModelData[iType].shape()[1]; ++j)
        for (size_t v=0; v<2; ++v)
  {
    fLocalData[iType][k][j][v] = 0;
    for (size_t m=0; m<fModelData[iType].shape()[3]; ++m)
      for (size_t l=0; l<fModelData[iType].shape()[4]; ++l)
        for (size_t u=0; u<2; ++u)
    {
      // skip coefficients which are multiplied with zero in expansion
      if ( (j==0 && v==1) || (l==0 && u==1) ) continue;
      const double z = iType == 0? zA : zB;
      if (u == 0)
        fLocalData[iType][k][j][v] += fModelData[iType][k][j][v][m][l][u]*pow(z, (double)m )*cos(l*fPhi);
      else
        fLocalData[iType][k][j][v] += fModelData[iType][k][j][v][m][l][u]*pow(z, (double)m )*sin(l*fPhi);
    }
  }
}


double ParametricMuonProfile::GetThetaMin() { return fThetaRange[0]; };


double ParametricMuonProfile::GetThetaMax() { return fThetaRange[1]; };