OpticalHalo.cc

Go to the documentation of this file.

#include <utl/AugerUnits.h>
#include <utl/TabularStream.h>
#include "OpticalHalo.h"

#include <cmath>
#include <iostream>

using namespace std;
using namespace utl;
using namespace FdProfileReconstructorKG;


OpticalHalo::OpticalHalo(EHaloType type) :
  fHaloType(type)
{ }


double
OpticalHalo::GetHaloFraction(const double zeta, const double age, const double dist)
  const
{
  switch (fHaloType) {
  case eNone:
    return 1;
  case eFlasher2005:
    return Flasher2005Fraction(zeta, age, dist);
  case eFlasher2008:
    return Flasher2008Fraction(zeta, age, dist);
  case eSpotGroup:
    return SpotGroupFraction(zeta, age, dist);
  case eSpotGroup2:
    return SpotGroupFraction2(zeta, age, dist);
  }

  return 1;
}


double
OpticalHalo::Flasher2005Fraction(const double zetaAngle, const double /*age*/, const double /*dist*/)
  const
{
  const double p       = 9.30207e-01;
  const double zeta_1  = 2.26745e-01;
  const double zeta_2  = 2.90549e+00;
  const double gamma   = 1.48972e+00;
  const double zeta    = zetaAngle/degree;

  const double f1 =      p  * (1 - exp(-zeta/zeta_1));
  const double f2 = (1 - p) * (1 - exp(-pow(zeta/zeta_2, gamma)));

  return f1 + f2;
}


double
OpticalHalo::Flasher2008Fraction(const double zetaAngle, const double age, const double dist)
  const
{
  // Flasher2008 analysis only valid above 1.4 degree
  const unsigned int n = 21;
  const double angles[n] = {
    1.4, 1.5, 1.6, 1.7, 1.8,
    1.9, 2,   3,   4,   5,
    6,   7,   8,   9,   12,
    15,  18,  21,  24,  27, 30
  };

  const double cumulative[n] = {
    0.859170, 0.861439, 0.863476,
    0.865431, 0.867303, 0.869093,
    0.870800, 0.883968, 0.892440,
    0.898886, 0.904640, 0.910163,
    0.915655, 0.921485, 0.937852,
    0.953002, 0.966389, 0.977545,
    0.986638, 0.994117, 1.000000
  };

  const double zeta = zetaAngle / degree;

  if (zeta < angles[0]) {
    const double norm = cumulative[0] / Flasher2005Fraction(angles[0], age, dist);
    return Flasher2005Fraction(zetaAngle, age, dist) * norm;
  } else if (zeta > angles[n-1])
    return 1;
  else {

    if (!fFlasher2008Table.GetNPoints()) {
      for (unsigned int i = 0; i < n; ++i)
        fFlasher2008Table.PushBack(angles[i], cumulative[i]);
    }

    return fFlasher2008Table.Y(zeta);

  }
}


double
OpticalHalo::SpotGroupFraction(const double zetaAngle,
                               const double /*age*/,
                               const double /*dist*/)
  const
{
  const double zeta = zetaAngle / degree;

  const double par[4] = { 0.738, 0.0308, 0.818, 0.643 };

  const double zeta0 = 0.7;
  const double norm1 = par[0] * (1 - exp(-zeta0/par[1]));
  const double norm2 = (1 - par[0]) * (1 - exp(-pow(zeta0/par[2], par[3])));
  const double norm = 1 - norm1 - norm2;

  const double f1 = par[0] * (1 - exp(-zeta/par[1]));
  const double f2 = (1 - par[0]) * (1 - exp(-pow(zeta/par[2], par[3])));
  const double f = 0.15 * (1 - f1 - f2);

  return 1 - f / norm;
}


double
OpticalHalo::SpotGroupFraction2(const double zetaAngle,
                                const double age,
                                const double dist)
  const
{
  const double zeta = zetaAngle / degree;
  const double maxDistance = 30 * km;
  const double d = fmin(dist, maxDistance) / km;

  // zeta --> ZetaOpt in degrees
  // d --> FD-shower distance in kms
  // age --> shower age

  const double x10 = d - 10;
  const double x10Squared = x10 * x10;

  const double a0 = 18.269707 - x10 * 0.539974 + x10Squared * 0.007218;
  const double a1 = (0.192196 + x10 * 0.007970 + x10Squared * 0.005581) *
    exp(-0.001762 * x10Squared);

  const double b0 =  5.251938 - x10 * 0.294560;
  const double b1 =  1.245041 - x10 * 0.129958;
  const double b2 = -0.666616 + x10 * 0.067478;

  const double z1 = zeta - 1;
  const double p0 = a0 * exp(-a1 * z1) / (1 + 2 * z1);
  const double p1 = (b0 + b1 * z1 + b2 * pow(z1, 2)) * exp(-z1);

  if (p1 < 0) {
    const double lightExcess = p0;
    return (lightExcess < 0) ? 1 : 1 - lightExcess / 100;
  } else {
    const double lightExcess = p0 + (age - 1) / 0.2 * p1;
    return (lightExcess < 0) ? 1 : 1 - lightExcess / 100;
  }

  /*
  // nOutside / nInside
  const double lightExcess = fmax(p1 < 0 ? p0 : p0 + (age - 1) / 0.2 * p1, 0) / 100;

  // nOutSide / (nOutside + nInside) = 1 / (1 + nOutside / nInside)
  return 1 / (1 + lightExcess);
  */

}

void
OpticalHalo::PrintFractions()
  const
{
  TabularStream tab("|.|.|.|.|.|");
  tab << hline
      << " zeta " << endc
      << " Flasher2005 " << endc
      << " Flasher2008 " << endc
      << " SpotGroup " << endc
      << " SpotGroup2 " << endr
      << hline;

  for (double dist = 1000; dist < 1.8; dist += 500)
    for (double age = 0.7; age < 1.8; age += 0.1)
      for (double zeta = 0.5*degree; zeta < 10*degree; zeta += 0.5*degree)
        tab << zeta / degree << endc
            << Flasher2005Fraction(zeta, age, dist) << endc
            << Flasher2008Fraction(zeta, age, dist) << endc
            << SpotGroupFraction(zeta, age, dist) << endc
            << SpotGroupFraction2(zeta, age, dist) << endr;
  tab << hline;
  cout << "\n" << tab << endl;
}