AtmCorr.cc

Go to the documentation of this file.

#include "AtmCorr.h"
#include "Atmosphere.h"

using namespace AtmCorrNS;


AtmCorr::AtmCorr(int DetectorType)
{
  if (DetectorType < 0 || DetectorType > 2) {
    printf("AtmCorr: Incorrect Detector type %d \n", DetectorType);
    exit(1);
  }
  fDetectorType = DetectorType;
};


//----------------------------------------------------------------------------------------------------------------------------
//             Type I atmospheric corection:
//                       Multiplying factor to convert signal for a shower with (r,rhoGround) to the reference case defined as:
//                       Density at ground 1.043e-3 g/cm^3

double
AtmCorr::GetAlphaRho(double r, int icomp) // Alpha rho was obtained from simulations with hground=1452.e2 for Water Cherenkov
{
  double r_km = r / 1e5;
  double alphaRho = parAtmTypeI[0][icomp] * pow(r_km, parAtmTypeI[1][icomp]);
  return alphaRho;
}


double
AtmCorr::GetCorrTypeI(double r, int icomp, double rhoGround, double hGround)
{
  if (fDetectorType == 2)
    return 1;

  AtmosphereNS::Atmosphere fatm;
  const double rhoGroundRef = fatm.GetDensity(hGround);

  //Moliere correction
  if (icomp > 0 && icomp < 4) {
    double alphaRho = GetAlphaRho(r, icomp);
    double fcorr = (1 - alphaRho * (rhoGround - rhoGroundRef) * 1e3 / 0.03);
    return 1 / fcorr;
  } else
    return 1;
}


//----------------------------------------------------------------------------------------------------------------------------
//             Type II atmospheric corection:
//                       Multiplying factor to convert signal for a shower with theta to the reference zenith angle 0 [deg]


double
AtmCorr::GetCorrTypeII(double r, double theta, int icomp, double hground)
{
  if (fDetectorType == 2)
    return 1;

  AtmosphereNS::Atmosphere fatm;
  const double rhoAboveGroundRef = fatm.GetDensityAboveGround(hground, 0., 2); //theta=0

  if (icomp > 0 && icomp < 4) {
    // double r_km = r / 1.e5;
    double alphaRho = GetAlphaRho(r, icomp);
    double rhoAboveGround = fatm.GetDensityAboveGround(hground, theta, 2);

    double fcorr = (1. - alphaRho * (rhoAboveGround - rhoAboveGroundRef) * 1.e3 / 0.03);
    return 1 / fcorr;
  } else
    return 1;
}


//----------------------------------------------------------------------------------------------------------------------------
//             Type III atmospheric corection + Other possible effects
//                       Multiplying factor to convert signal for a shower with psi to the reference psi angle 90 [deg]

double
AtmCorr::GetParCorrTypeIII(double r,  int icomp, int ipar)
{
  if (fDetectorType == 2 && icomp > 0)
    return 0;
  if (ipar < 0 || ipar > 1)
    return 0;

  const double* par = 0;
  if (ipar == 0 && fDetectorType == 0)
    par = f_mod1_par_WCD[icomp];
  else if (ipar == 1 && fDetectorType == 0)
    par = f_mod2_par_WCD[icomp];
  else if (ipar == 0 && fDetectorType == 1)
    par = f_mod1_par_Scin[icomp];
  else if (ipar == 1 && fDetectorType == 1)
    par = f_mod2_par_Scin[icomp];
  else if (ipar == 0 && fDetectorType == 2)
    par = f_mod1_par_MD[icomp];
  else if (ipar == 1 && fDetectorType == 2)
    par = f_mod2_par_MD[icomp];
  else {
    printf("unknown case\n");
    exit(1);
  }

  return par[0] + par[1] * r / 1.e5 + par[2] * r * r / 1e10;
}


double
AtmCorr::GetCorrTypeIII(double r, double theta, double psi, int icomp)
{
  if (fDetectorType == 2 && icomp > 0)
    return 1;

  //f_mod1/f_mod2
  const double f_mod1 = GetParCorrTypeIII(r, icomp, 0);
  const double f_mod2 = GetParCorrTypeIII(r, icomp, 1);

  //----

  const double theta_deg = theta * 180 / M_PI;
  const double f_mod = f_mod1 * theta_deg / 60 + f_mod2 * pow(theta_deg / 60, 2);
  const double fcorr = (1 + cos(psi) * f_mod);

  return 1 / fcorr;
}