USCMuonProfile.cc

Go to the documentation of this file.

#include "USCMuonProfile.h"
#include <tls/TankResponseUtilities.h>

#include <vector>
#include <sstream>
#include <cmath>
#include <string>

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

#include <TH2D.h>
#include <TFile.h>

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

inline
double
InterpolateHistrogram(const TH2D& h, double x, double y)
{
  double x1(0),x2(0),y1(0),y2(0);

  int ix = h.GetXaxis()->FindBin(x);
  int iy = h.GetYaxis()->FindBin(y);
  if (ix == 1) ix += 1;
  if (ix == h.GetNbinsX()+1) ix -= 1;

  x1 = h.GetXaxis()->GetBinCenter(ix);
  if (x<x1) {
    --ix;
    x1 = h.GetXaxis()->GetBinCenter(ix);
  }
  x2 = h.GetXaxis()->GetBinCenter(ix+1);

  y1 = h.GetYaxis()->GetBinCenter(iy);
  if (y<y1) {
    --iy;
    y1 = h.GetYaxis()->GetBinCenter(iy);
  }
  y2 = h.GetYaxis()->GetBinCenter(iy+1);

  const double v00 = h.GetBinContent(ix,iy);
  const double v01 = h.GetBinContent(ix,iy+1);
  const double v10 = h.GetBinContent(ix+1,iy);
  const double v11 = h.GetBinContent(ix+1,iy+1);
  const double zx = (x-x1)/(x2-x1);
  const double zy = (y-y1)/(y2-y1);
  return (1.0-zx)*(1.0-zy)*v00 + (1.0-zx)*zy*v01 + zx*(1.0-zy)*v10 + zx*zy*v11;
}


inline
double
Extrapolate(const TH2D& h, double rmax, double r, double zeta)
{ // extrapolation with formula ln N_\mu = a + b sqrt{r}
  const double a = log( InterpolateHistrogram( h,rmax*cos(zeta),rmax*sin(zeta) ) );
  const double b = ( a - log( InterpolateHistrogram( h,0.9*rmax*cos(zeta),0.9*rmax*sin(zeta) ) ) )/(sqrt(rmax)-sqrt(0.9*rmax));
  return exp(a + b*(sqrt(r)-sqrt(rmax)));
}


inline
double
LocalMuonAngle(double x, double y, double theta, double phi)
{
  // Calculate approximate distance of muon production point
  // with Chebychev polynomial

  const double c[64] = { 1.323609e+05,  1.187303e+05,  9.395542e+04,  6.895397e+04,
                         4.840035e+04,  3.251341e+04,  2.091493e+04,  1.282508e+04,
                         7.433439e+03,  4.007066e+03,  1.945498e+03,  7.860169e+02,
                         1.895583e+02, -7.845296e+01, -1.691854e+02, -1.726411e+02,
                        -1.390602e+02, -9.562383e+01, -5.631717e+01, -2.703336e+01,
                        -8.487988e+00,  1.335395e+00,  5.137142e+00,  5.495047e+00,
                         4.413872e+00,  3.072145e+00,  1.835842e+00,  6.221296e-01,
                        -6.599210e-01, -1.829552e+00, -2.470651e+00, -2.467011e+00,
                        -1.996112e+00, -1.483970e+00, -1.217739e+00, -1.242982e+00,
                        -1.399658e+00, -1.550564e+00, -1.492873e+00, -1.091996e+00,
                        -4.217213e-01,  3.443991e-01,  1.022913e+00,  1.490880e+00,
                         1.786550e+00,  2.054111e+00,  2.299077e+00,  2.372312e+00,
                         2.122376e+00,  1.491762e+00,  5.303430e-01, -4.858504e-01,
                        -1.286566e+00, -1.638024e+00, -1.536687e+00, -1.157604e+00,
                        -8.427614e-01, -9.519352e-01, -1.534490e+00, -2.439957e+00,
                        -3.454757e+00, -4.376512e+00, -4.991014e+00, -5.234057e+00  };

  const unsigned int n = 64;
  const double z = 2*theta/(0.5*utl::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 + c[i];
    dd = sv;
  }
  d = (z*d + 0.5*c[0] - dd)*utl::meter;

  const double cosThetaMu = cos(theta)/sqrt( pow(cos(theta),2) + pow(d*cos(phi)*sin(theta)-x,2) + pow(d*sin(phi)*sin(theta)-y,2) );
  return acos(cosThetaMu);
}


USCMuonProfile::USCMuonProfile(Branch topBranch):
  fDensity(0), fThetaRange(),
  fTheta(0.), fPhi(0.), fEnergy(10.*EeV),
  fExponent(1.0), fNMuonScale(1.0), fNMuonScaleEnergy(1.0),
  fIndexThetaLow(0), fIndexThetaUp(0), fIndexPhiLow(0), fIndexPhiUp(0),
  fHistDen1(NULL),fHistDen2(NULL),fHistDen3(NULL),fHistDen4(NULL)
{
  Branch b;
  b = topBranch.GetChild("nMuModel");
  b.GetChild("exponent").GetData(fExponent);
  b.GetChild("scale").GetData(fNMuonScale);

  topBranch.GetChild("thetaRange").GetData(fThetaRange);

  string pathToDataFile;
  topBranch.GetChild("directory").GetData(pathToDataFile);

  fDensity = TFile::Open(pathToDataFile.c_str(),"READ");
  if (!fDensity) {
    ostringstream msg;
    msg << "Could not read histogram file " << pathToDataFile << "\n";
    throw NoDataForModelException(msg.str());
  }

  SetProfile(fTheta, fPhi, fEnergy);
}

void
USCMuonProfile::LoadHistograms(double theta, double phi)
{
  const double thetaDeg = theta/degree;
  const double phiDeg   = phi  /degree;

  fIndexThetaLow = Theta2IndexLower(thetaDeg);
  fIndexThetaUp  = Theta2IndexUpper(thetaDeg);

  const double phiMapDeg = PhiMap(phiDeg);
  fIndexPhiLow = Phi2IndexLower(phiMapDeg);
  fIndexPhiUp  = Phi2IndexUpper(phiMapDeg);

  unsigned int id1 = 60000+2000*(fIndexThetaLow)+fIndexPhiLow;
  unsigned int id2 = 60000+2000*(fIndexThetaLow)+fIndexPhiUp;
  unsigned int id3 = 60000+2000*(fIndexThetaUp)+fIndexPhiLow;
  unsigned int id4 = 60000+2000*(fIndexThetaUp)+fIndexPhiUp;

  ostringstream histoID1;
  histoID1 << "h" << id1 ;

  ostringstream histoID2;
  histoID2 << "h" << id2 ;

  ostringstream histoID3;
  histoID3 << "h" << id3 ;

  ostringstream histoID4;
  histoID4 << "h" << id4 ;

  fHistDen1 = static_cast<TH2D*>(fDensity->Get(histoID1.str().c_str()));
  fHistDen2 = static_cast<TH2D*>(fDensity->Get(histoID2.str().c_str()));
  fHistDen3 = static_cast<TH2D*>(fDensity->Get(histoID3.str().c_str()));
  fHistDen4 = static_cast<TH2D*>(fDensity->Get(histoID4.str().c_str()));

  if(!fHistDen1 || !fHistDen2 || !fHistDen3 || !fHistDen4 ) {
    ostringstream msg;
    msg << "the following histograms are missing in file "
        << fDensity->GetName() << ": ";
    if(!fHistDen1) msg << histoID1.str().c_str() << " ";
    if(!fHistDen2) msg << histoID2.str().c_str() << " ";
    if(!fHistDen3) msg << histoID3.str().c_str() << " ";
    if(!fHistDen4) msg << histoID4.str().c_str() << " ";
    msg  << endl << "\tTheta: " << thetaDeg
         << " Phi: " << phiDeg << endl;
    throw NoDataForModelException(msg.str());
  }
}

// Angle in AugerUnits
double
USCMuonProfile::NMuon(double xpos, double ypos, double theta, double phi, double energy)
{
  SetProfile(theta,phi,energy);

  const double rperp = rPerpendicular( xpos,  ypos, fTheta, fPhi);
  const double psi   = PsiPerpendicular( xpos,  ypos, fTheta, fPhi);

  const double zeta = psi-fBPsi+0.5*kPi;
  const double xperp = rperp*cos(zeta);
  const double yperp = rperp*sin(zeta);

  double dum1,dum2,dum3,dum4;
  if (fabs(xperp) < 3950. and fabs(yperp) < 3950.)
  {
    dum1 = InterpolateHistrogram(*fHistDen1,xperp,yperp);
    dum2 = InterpolateHistrogram(*fHistDen2,xperp,yperp);
    dum3 = InterpolateHistrogram(*fHistDen3,xperp,yperp);
    dum4 = InterpolateHistrogram(*fHistDen4,xperp,yperp);
  }
  else
  { // extrapolation with formula ln N_\mu = a + b sqrt{r}
    double rmax1 = 1e99;
    double rmax2 = 1e99;
    if (fabs(xperp)>3950.) rmax1 = fabs(3950./cos(zeta));
    if (fabs(yperp)>3950.) rmax2 = fabs(3950./sin(zeta));
    const double rmax = rmax1 < rmax2 ? rmax1 : rmax2;
    dum1 = Extrapolate(*fHistDen1,rmax,rperp,zeta);
    dum2 = Extrapolate(*fHistDen2,rmax,rperp,zeta);
    dum3 = Extrapolate(*fHistDen3,rmax,rperp,zeta);
    dum4 = Extrapolate(*fHistDen4,rmax,rperp,zeta);
  }

  double thetaMap = fTheta/utl::degree;
  double phiMap   = PhiMap(fPhi/utl::degree);

  double dummy;
  if (fIndexThetaLow==fIndexThetaUp)
  {
      if (fIndexPhiLow==fIndexPhiUp)
      {
        //do not interpolate. theta phi are exact
        dummy = dum1;
      }
      else
      {
        //interpolate only in phi
        double u1 = Index2Phi(fIndexPhiLow);
        double u2 = Index2Phi(fIndexPhiUp);

        double u = (phiMap-u1)/(u2 -u1);
        dummy = (1-u)*dum1+u*dum2;
      }
  }
  else
  {
    if (fIndexPhiLow==fIndexPhiUp)
    {
      //do not interpolate on phi. interpolate on theta.
      double t1 = Index2Theta(fIndexThetaLow);
      double t2 = Index2Theta(fIndexThetaUp);

      double t = (thetaMap-t1)/(t2-t1);

      dummy = (1-t)*dum1 +t*dum3;
    }
    else
    {
      //interpolate on theta and phi
      double t1 = Index2Theta(fIndexThetaLow);
      double t2 = Index2Theta(fIndexThetaUp);

      double t = (thetaMap-t1)/(t2-t1);

      double u1 = Index2Phi(fIndexPhiLow);
      double u2 = Index2Phi(fIndexPhiUp);

      double u = (phiMap-u1)/(u2 -u1);
      dummy = (1-t)*(1-u)*dum1 + u*(1-t)*dum2 + t*u*dum4 + t*(1-u)*dum3;
    }
  }

  const double localTheta = ThetaMuon(xpos,ypos,theta,phi,energy);
  const double scale = fNMuonScaleEnergy * cos(localTheta) * TankGroundArea(localTheta);
  return dummy*scale;
}

// Angle in AugerUnits
double
USCMuonProfile::ThetaMuon(double /*xpos*/, double /*ypos*/, double theta, double /*phi*/, double /*energy*/)
{
  // currently no own theta profile implemented in the USC case
  //   return LocalMuonAngle(xpos,ypos,theta,phi);
  return theta;
}

// Angle in AugerUnits
double
USCMuonProfile::EnergyMuon(double xpos, double ypos, double /*theta*/, double /*phi*/, double /*energy*/)
{
  // taken from HAS Model, paper Maximo et al.
  double r  = rPerpendicular(xpos,ypos,fTheta,fPhi);
  if (r>4000.*m) r = 4000*m;

  double th = fTheta/deg;

  double Ap0 = -1.92124e+01;
  double Ap1 =  8.75452e-01;
  double Ap2 = -1.19709e-02;
  double Ap3 =  5.74630e-05;
  double A = Ap0 + Ap1*th + Ap2*th*th + Ap3*th*th*th;


  double Bp0 =  3.65857e+00;
  double Bp1 = -1.61666e-01;
  double Bp2 =  1.97280e-03;
  double Bp3 = -8.13659e-06;
  double B = Bp0 + Bp1*th + Bp2*th*th + Bp3*th*th*th;

  double LogEMu = A + B*log10(r); //Eq.(8)
  double E = pow(10,LogEMu); //Energy in GeV

  //range (1GeV, 1TeV)
  if (std::isnan(E) || E<1.0) E = 1.0;
  else if (E>1000.) E=1000.;

  return E*GeV;
}

// Internal function, angle in degrees
unsigned int
USCMuonProfile::Theta2IndexLower(double theta)
{
  theta = approx(theta);
  return static_cast<unsigned int>(floor((theta - 60.0)/2.0));
}

// Internal function, angle in degrees
unsigned int
USCMuonProfile::Theta2IndexUpper(double theta)
{
  theta = approx(theta);
  return static_cast<unsigned int>(ceil((theta - 60.0)/2.0));
}

// Internal function, angle in degrees
unsigned int
USCMuonProfile::Phi2IndexLower(double phi)
{
  phi = approx(phi);
  return static_cast<unsigned int>(floor(phi/5.0))*5;
}

// Internal function, angle in degrees
unsigned int
USCMuonProfile::Phi2IndexUpper(double phi)
{
  phi = approx(phi);
  return static_cast<unsigned int>(ceil(phi/5.0))*5;
}

// Internal function, angle in degrees
double
USCMuonProfile::Index2Theta(unsigned int itheta)
{
    return 60.0+ static_cast<double>(2*(itheta));
}

// Internal function, angle in degrees
double
USCMuonProfile::Index2Phi(unsigned int iphi)
{
    return 5.0*static_cast<double>(iphi/5);
}

// Internal function, angle in degrees
double
USCMuonProfile::PhiMap(double phi)
{

  double phiAires = 90.0 - phi - tls::kDeclination/degree;

    if (phiAires<0) phiAires = 360+phiAires;

    if (phiAires>180) phiAires = 360 - phiAires;

    return phiAires;
}

double USCMuonProfile::GetThetaMin() { return fThetaRange[0]; };
double USCMuonProfile::GetThetaMax() { return fThetaRange[1]; };