ZombieGraveyard/LDFFinderOG/EnergyConversion.cc

Go to the documentation of this file.

#include <cmath>
#include <utl/Math.h>
#include "EnergyConversion.h"

using namespace std;
using namespace utl;
using namespace LDFFinderOG;


/*
  The energy estimate is rather independent of the LDF used
  (see the Auger ICRC 2005 contribution on LDFs).
  For now we use the energy conversion, which is given in the
  Auger ICRC 2007 paper on the energy spectrum (i.e. CIC).
*/
void
EnergyConversion::GetEnergy(const double cosTheta, const double s1000, const double s1000Err,
                            double& energy, double& energyErr)
  const
{
  const double cos2deg38 = 0.620960947799833907;  // Sqr(cos(38*degree))
  const double x = Sqr(cosTheta) - cos2deg38;
  const double cic = 1 + x * (fAttenuationPar1 + x * (fAttenuationPar2 + x * fAttenuationPar3));
  if (cic > 0) {
    const double s38 = s1000 / cic;
    energy = fEnergyS38Const * pow(s38, fEnergyS38Slope);
    energyErr = fEnergyS38Slope * energy * s1000Err / s1000;
  } else
    energy = energyErr = 0;
}


/*void
EnergyConversion::GetEnergyICRC2005(const double cosTheta, const double s1000, const double s1000Err,
                                    double& energy, double& energyErr)
  const
{
  
}*/


  /* ICRC05 parameters
  const double theta = fgShowerAxis.GetTheta(fgLocalCS)/degree;
  const double cic = 1.049 + 0.00974 * theta - 0.00029 * Sqr(theta);
  const double fS38 = fi.fS1000/cic;

  const double b0 = 0.16;
  const double b1 = 1.06;
  */

  /*const double cosTheta = fgShowerAxis.GetCosTheta(fgLocalCS);
  const double cos2deg38 = 0.620960947799833907;  // Sqr(cos(38*degree))
  const double normalizedAngle = Sqr(cosTheta) - cos2deg38;*/
  /* ICRC07 parameters
  const double cic = 1 + 0.92*normalizedAngle - 1.13*Sqr(normalizedAngle);
  const double b0 = 0.149;
  const double b1 = 1.078;
  */
  /*const double cic = 1 + fEnergyConversion.fAttenuationPar1 * normalizedAngle
    + fEnergyConversion.fAttenuationPar2 * Sqr(normalizedAngle)
    + fEnergyConversion.fAttenuationPar3 * pow(normalizedAngle, 3);

  const double fS38 = fi.fS1000/cic;
  if (cic > 0) {
    fi.fEnergy = fEnergyConversion.fEnergyS38Const * pow(fS38, fEnergyConversion.fEnergyS38Slope);
    fi.fEnergyErr = fEnergyConversion.fEnergyS38Slope * fi.fEnergy * fi.fS1000Err/fi.fS1000;
  } else
    fi.fEnergy = fi.fEnergyErr = 0;*/

  /* old parameterization: see P. Summers ICRC 2005

  const double secTheta = 1/fgShowerAxis.GetCosTheta(fgLocalCS);

  if (fgLDFType == ePL) {
    const double b0 = 1/7.8 * sqrt(1 + 11.8*Sqr(secTheta - 1));
    const double b1 = 1.052;
    const double e = b0 * fi.fS1000;
    fi.fEnergy = pow(e, b1) * EeV;
    fi.fEnergyErr = b1 * fi.fEnergy * fi.fS1000Err/fi.fS1000;
  } else {
    const double secTheta2 = Sqr(secTheta);
    const double b0 = 0.37 - 0.51*secTheta + 0.30*secTheta2;
    const double b1 = 1.27 - 0.27*secTheta + 0.08*secTheta2;
    fi.fEnergy = b0 * pow(fi.fS1000, b1) * EeV;
    fi.fEnergyErr = b1 * fi.fEnergy * fi.fS1000Err/fi.fS1000;
  }*/