UnivCalibConstants.cc

Go to the documentation of this file.

#include "UnivCalibConstants.h"

using namespace UnivCalibConstantsNS;


UnivCalibConstants::UnivCalibConstants(int CalibOpt_i, int RecSys_i, int RecMixture_i)
{
  CalibOpt = CalibOpt_i;
  RecSys = RecSys_i;
  HadronicModel = (CalibOpt_i >= 0 ? CalibOpt : UNDEF);
  RecMixture = RecMixture_i;

  //printf("UnivCalibConstants: CalibOpt=%d  RecSys=%d HadronicModel=%d RecMixture=%d \n",CalibOpt,RecSys,HadronicModel,
  //   RecMixture);

  if (CalibOpt < -3 || CalibOpt > 3) {
    printf("UnivCalibConstants: Fatal error CalibOpt \n");
    exit(1);
  }

  if (CalibOpt < 0) {
    RecMixture = UNDEF;
    HadronicModel = UNDEF;
  } else if (RecMixture < 0 || RecMixture >= nRecMixtures || HadronicModel < 0 || HadronicModel > 3) {
    printf("UnivCalibConstants: Fatal error RecMixture \n");
    exit(1);
  }

  if (CalibOpt <= -2  && (RecSys < 0 || RecSys >= nSys)) {
    printf("UnivCalibConstants: Fatal error RecSys %d %d \n", CalibOpt, RecSys);
    exit(1);
  }
}


//-----------------------------------------------------------------------------------
//--------- Proton fraction
//-----------------------------------------------------------------------------------
double
UnivCalibConstants::GetProtonFraction(double logE)
{
  double fp = UNDEF;
  if (CalibOpt < 0)
    return fp;

  const double fProton19[nRecMixtures] = {1, 0, 0.5, 0.1, 0.9, 1};
  const double fProton20[nRecMixtures] = {1, 0, 0.5, 0.1, 0.9, 0};

  fp = fProton19[RecMixture] + (logE - 19) * (fProton20[RecMixture] - fProton19[RecMixture]);

  if (fp < 0)
    fp = 0;

  if (fp > 1)
    fp = 1;

  return fp;
}


//-----------------------------------------------------------------------------------
//---------<Xmax>
//-----------------------------------------------------------------------------------
double
UnivCalibConstants::GetMeanXmax(double logE)
{
  double meanXmax = UNDEF;

  // Data
  if (CalibOpt <= -2) {
    logE -= log10(fEnergySys_Auger[RecSys]);

    const double logE0 = 18.36, Xmax0 = 754, ER0 = 71, ER1 = 19;

    if (logE < logE0)
      meanXmax = Xmax0 + (logE - logE0) * ER0;
    else
      meanXmax = Xmax0 + (logE - logE0) * ER1;

    meanXmax += fXmaxSys_Auger[RecSys];

  } else if (CalibOpt == -1) { // Photon analysis
    const double Xmax_PhotonSearch = 1050;
    meanXmax = Xmax_PhotonSearch;

  } else { // MC
    const double Xmax0_Models[4][2] =  {{788, 703}, {819, 712}, {790, 695}, { 805, 709}}; //[HadronicModel][Primary]
    const double A_Xmax_Models[4][2] = {{43, 51}, {63, 62}, {55, 56}, {56, 59}};

    double meanXmax_i[2]; // Proton/Iron
    for (int ip = 0; ip < 2; ++ip)
      meanXmax_i[ip] = Xmax0_Models[HadronicModel][ip] + (logE - logE_Calib) * A_Xmax_Models[HadronicModel][ip];

    double fp = GetProtonFraction(logE);
    meanXmax = meanXmax_i[0] * fp + meanXmax_i[1] * (1 - fp);
  }

  return meanXmax;
}


//-----------------------------------------------------------------------------------
//---------<Nmu>
//-----------------------------------------------------------------------------------
double
UnivCalibConstants::GetMeanNmu(double logE, double theta, double Xmax)
{
  double MeanNmu = GetMeanNmu(logE, theta);

  // <Nmu>
  if (Xmax < 0 || CalibOpt == -1)
    return MeanNmu;

  // Using Nmu(Xmax) correlation and normalizing so Nmu(<Xmax>)=<Nmu>
  double MeanXmax_pFe = GetMeanXmax_pFe(logE);
  double dX0 = GetMeanXmax(logE) - MeanXmax_pFe;
  double dX = Xmax - MeanXmax_pFe;
  double Nmu = MeanNmu * (1. + 0.5 * atan(-dX / 40) / M_PI) / (1. + 0.5 * atan(-dX0 / 40) / M_PI) ;

  return Nmu;
}


double
UnivCalibConstants::GetMeanNmu(double logE, double theta)
{
  double meanNmu = UNDEF;

  // Calib Data
  if (CalibOpt == -3)
    return meanNmu;
  else if (CalibOpt == -2) {  // Data
    logE -= log10(fEnergySys_Auger[RecSys]);

    const double Nmu0_Auger = 1.99, A_Nmu_Auger = 0.21, B_Nmu_Auger = 0.17;
    meanNmu = Nmu0_Auger + A_Nmu_Auger * (logE - logE_Calib) + B_Nmu_Auger * (1 / cos(theta) - 2) ;

    meanNmu *= fNmuSys_Auger[RecSys];
  } else if (CalibOpt == -1) { // Photon Analysis
    const double Nmu_PhotonSearch = 0.15;
    meanNmu = Nmu_PhotonSearch;
  } else { // MC
    const double Nmu0_Models[4][2] = {{0.98, 1.38}, {1.19, 1.70}, {1.20, 1.69}, {1.28, 1.82}}; // Nmu at theta=60 [deg] log(E)=19.5 [eV]
    const double A_Nmu_Models[4][2] = {{ -0.04, 0.02 }, { -0.14, -0.12}, {0.07, 0.13}, {0.01, 0.10}}; // Zenith dependence
    const double B_Nmu_Models[4][2] = {{ -0.01, -0.09}, { -0.04, -0.11}, { -0.01, -0.10}, { -0.01, -0.11}}; // Energy evolution

    double meanNmu_i[2];
    for (int ip = 0; ip < 2; ++ip)
      meanNmu_i[ip] = Nmu0_Models[HadronicModel][ip] + A_Nmu_Models[HadronicModel][ip] * (1 / cos(theta) - 2) + B_Nmu_Models[HadronicModel][ip] * (logE - logE_Calib);

    const double fp = GetProtonFraction(logE);
    meanNmu = meanNmu_i[0] * fp + meanNmu_i[1] * (1 - fp);
  }

  return meanNmu;
}


//-----------------------------------------------------------------------------------
//---------  OffsetM_Mu
//-----------------------------------------------------------------------------------
double
UnivCalibConstants::GetOffsetM_Mu(double theta)
{
  return GetOffsetM_Mu(UNDEF, theta, UNDEF);
}


double
UnivCalibConstants::GetOffsetM_Mu(double /*logE*/, double theta, double /*Xmax*/)
{
  double offsetM_Mu[2] = {UNDEF, UNDEF};

  if (CalibOpt == -3) { // Calib Data
    offsetM_Mu[0] = OffsetM_Mu_Calib[0];
    offsetM_Mu[1] = OffsetM_Mu_Calib[1];
  } else if (CalibOpt == -2) { // Data
    offsetM_Mu[0] = OffsetM_Mu_Auger[0][RecSys];
    offsetM_Mu[1] = OffsetM_Mu_Auger[1][RecSys];
  } else if (CalibOpt == -1) { // Photon Analysis
    offsetM_Mu[0] = OffsetM_Mu_Photon[0];
    offsetM_Mu[1] = OffsetM_Mu_Photon[1];
  } else { // MC
    offsetM_Mu[0] = OffsetM_Mu_Models[0][HadronicModel];
    offsetM_Mu[1] = OffsetM_Mu_Models[1][HadronicModel];
  }

  const double offset = offsetM_Mu[0] + offsetM_Mu[1] * (1 / cos(theta) - 1);

  return offset;
}


//----------------------------------------------------------------------------------------------------
//--------- MC    ->  <Xmax>(logE) for a 50% p/Fe mixture
//--------- Data  ->  <Xmax>(logE) from Nmu(Xmax) fit at log(E)=19 [eV]
//----------------------------------------------------------------------------------------------------
double
UnivCalibConstants::GetMeanXmax_pFe(double logE)
{
  double MeanXmax_pFe, ER_pFe;
  if (CalibOpt == -3) {
    MeanXmax_pFe = MeanXmax_pFe_Calib;
    ER_pFe = ER_pFe_Calib;
  } else if (CalibOpt == -2) {
    logE -= log10(fEnergySys_Auger[RecSys]);
    MeanXmax_pFe = MeanXmax_pFe_Auger + Offset_MeanXmax_pFe_Auger[RecSys] + fXmaxSys_Auger[RecSys];
    ER_pFe = ER_pFe_Auger  + Offset_ER_pFe_Auger[RecSys];
  } else if (CalibOpt == -1) {
    MeanXmax_pFe = MeanXmax_pFe_Models[0];
    ER_pFe = ER_pFe_Models[0];
  } else {
    MeanXmax_pFe = MeanXmax_pFe_Models[HadronicModel];
    ER_pFe = ER_pFe_Models[HadronicModel];
  }

  return MeanXmax_pFe + (logE - logE_Calib) * ER_pFe;
}


//----------------------------------------------------------------------------------------------------
//---------- Xmax_SD / Xmax_FD /  Nmu_FD resolution /  RMS Xmax( RMS XmaxSD)
//----------------------------------------------------------------------------------------------------
double
UnivCalibConstants::GetSigmaXmax_SD(double logE, int ndev, bool IsInfill)
{
  double A_p[2] = {51.5, 25.5}, B_p [2] = {0.8, 0.5}; // Proton
  double A_fe[2] = {43.5, 21}, B_fe[2] = {0.8, 0.76}; // Iron

  unsigned int infill = IsInfill;

  double sigma_p = A_p[infill] * exp(-(logE - 19) * B_p[infill]);
  double sigma_fe = A_fe[infill] * exp(-(logE - 19) * B_fe[infill]);

  double sigma = 0.5 * (sigma_p + sigma_fe) + 0.5 * ndev * (sigma_p - sigma_fe);

  return sigma;
}


double
UnivCalibConstants::GetSigmaXmax_SD(double logE, int ndev)
{
  return GetSigmaXmax_SD(logE, ndev, false);
}


double
UnivCalibConstants::GetSigmaXmax_FD(double /*logE*/)
{
  return 20;
}


double
UnivCalibConstants::GetSigmaNmu_FD(double logE)
{
  return 0.2 - 0.05 * (logE - 19); // Nmu/Nmu(true)
}


double
UnivCalibConstants::GetTrueRMSXmax(double rms, double logE, bool IsInfill)
{
  double slope = 1.025 - 0.1 * (logE - 19);
  double sigma = (IsInfill ? 17.5 * exp(-(logE - 19) * 1) : 43 * exp(-(logE - 19) * 0.8));
  if (rms / slope < sigma)
    return UNDEF;
  return sqrt(rms * rms / slope / slope - sigma * sigma);
}