CherenkovModel.cc

Go to the documentation of this file.

/*
  \file
  Implementation of CherenkovModel

  \autor Diego Melo.
  \date 18 Dic 2003
*/


#include <iostream>
#include <string>
#include <vector>

#include <fwk/CentralConfig.h>

#include <utl/Vector.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>

#include <det/Detector.h>

#include <atm/Atmosphere.h>
#include <atm/CherenkovModel.h>
#include <atm/AttenuationResult.h>

using namespace std;
using namespace fwk;
using namespace utl;
using namespace det;
using namespace atm;


CherenkovModel::CherenkovModel()
{ }


// Gaisser-Hillas
double
CherenkovModel::GaisserHillas(const double chi, const double chi0, const double chimax, const double Nmax, const double lambda)
  const
{
  if (chi < chi0)
    return 0;

  const double x = (chi - chi0) / (chimax - chi0);
  const double p = (chimax - chi0) / lambda;
  const double q = (chimax - chi) / lambda;
  return Nmax * std::pow(x, p) * std::exp(q);
}


// Energy(beta)
double
CherenkovModel::Energy(const double beta)
  const
{
  return utl::ElectronMass/std::sqrt(1 - beta*beta);
}


// Beta(energy)
double
CherenkovModel::Beta(const double energy)
  const
{
  const double gammaInv = utl::ElectronMass / energy;
  return std::sqrt(1 - gammaInv*gammaInv);
}


// Threshold de Energia Cherenkov.
double
CherenkovModel::CKVThreshold(const double nIndex)
  const
{
  return utl::ElectronMass / std::sqrt(1 - 1 / (nIndex * nIndex));
}


// Distance between two points (m)
double
CherenkovModel::Distance(const double x1[], const double x2[])
  const
{
  double d2 = 0;

  for (int k = 0; k < 3; ++k) {
    const double dx = x1[k] - x2[k];
    d2 += dx * dx;
  }

  return std::sqrt(d2);
}


// Distance between two points (g/cm2).
double
CherenkovModel::GDistance(const double x1[], const double x2[])
  const
{
  double gdistance = 0;
  const double chi1 = det::Detector::GetInstance().GetAtmosphere().EvaluateDepth(x1[2]);
  const double chi2 = det::Detector::GetInstance().GetAtmosphere().EvaluateDepth(x2[2]);

  if (x1[2] != x2[2])
    return Distance(x1, x2) * std::abs((chi1 - chi2) / (x1[2] - x2[2]));
  if (x1[2] <= 12*utl::km)
    return Distance(x1, x2) * std::abs(chi1 / (0.19*(44340 - x1[2])));
  if (x1[2] <= 45.470*utl::km)
    return Distance(x1, x2) * std::abs(chi1 / (6.340*utl::km));
  return 0;
}


// Hillas -J.Phys G: Nucl. Phys. 8, p1461, 1982-
double
CherenkovModel::TrackLengthH(const double energy, const double chi, const double chimax)
  const
{
  const double s = 3 / (1 + 2*chimax/chi);
  const double e0 = (s < 0.4) ? 26 : 44 - 17*std::pow(s - 1.46, 2);

  return std::pow((0.89*e0 - 1.20) / (e0 + energy), s) / std::pow(1 + 1e-4 * s * energy, 2);
}


// Nerling -Malargue Nov. 2003-
double
CherenkovModel::TrackLengthN(const double energy, const double chi, const double chimax)
  const
{
  const double s = 3 / (1 + 2*chimax/chi);

  const double a0 = 102.24;                //Estimacion por el momento !!!!
  const double a1 = 6.425 - 1.532 * s;
  const double a2 = 168.13 - 42.14 * s;

  return a0 / ((energy + a1)*std::pow(energy + a2, s));
}


// Fotones Cherenkov Producidos por una unica Particula de Energia E
double
CherenkovModel::CKVPhot(const double trl, const double nIndex, const double beta, const double dl, const double l1, const double l2)
  const
{
  const double cte = 2*utl::kPi * utl::Alpha;

  return cte * trl * std::abs(1/l1 - 1/l2) * (1 - 1/std::pow(beta*nIndex, 2)) * dl;
}


// Funcion Auxiliar Empleada  para Integrar la Luz Cherenkov
double
CherenkovModel::Fu(const double x, const double a)
  const
{
  return std::exp(-x/a) * std::sin(x);
}


// Integral Angular Light Cherenkov
void
CherenkovModel::CKVIntegral(const double xlo, const double xhi, const double xb, double& integral)
  const
{
  double s2 = 0;

  double suma = Fu(xlo, xb);
  const double step = (xhi - xlo) / 10;

  for (int i = 1; i <= 5; ++i) {
    const double a1 = xlo + (2*i - 1)*step;
    const double a2 = a1 + step;

    const double s1 = 4 * Fu(a1, xb);
    s2 = 2 * Fu(a2, xb);

    suma += s1 + s2;
   }

  integral = (suma - s2/2) * step / 3;
}


// Cherenkov Directo
void
CherenkovModel::CKVDirect(const double bin, const double theta, const double height, const double rp, const double chimax, std::vector<double>& dirCh)
  const
{
  std::vector<double> nphot(kN_WAVE_BIN);
  std::vector<double> ophot(kN_WAVE_BIN);
  double h1;
  double h2;

  for (int w = 0; w < kN_WAVE_BIN; ++w) {
    //nphot[w] = 0;
    //ophot[w] = 0;
    dirCh[w] = 0;
   }

  double th1 = theta - bin/2;
  if (th1 <= 0)
    th1 = theta/2;

  const double th2 = theta + bin/2;

  const double dl = rp * std::abs(1 / std::tan(th1) - 1 / std::tan(th2));

  const double a = 0.83;
  const double b = -0.67;

  const double chi = det::Detector::GetInstance().GetAtmosphere().EvaluateDepth(height);
  const double nIndex = det::Detector::GetInstance().GetAtmosphere().EvaluateRefractionIndex(height);

  const double em = min(CKVThreshold(nIndex), 2e3);

  const double th0 = a * std::pow(em, b);

  // Integrate in Theta
  double integral = 0;
  CKVIntegral(th1, th2, th0, integral);

  // Normalizes Integral
  const double norm = integral / th0;

  // Integrate in Energy
  const double dE = 10;

  // From Threshold to Ultrarelativistic "E = 361.33 <---> beta = 1 - 1e-6"
  double jmax = 361.33;

  double lambda = 0;
  for (double j = em; j <= jmax; j += dE) {
    const double tl1 = TrackLengthN(j, chi, chimax);
    const double tl2 = TrackLengthN(j+dE, chi, chimax);

    const double aux1 = tl1 - tl2;
    const double aux2 = j + dE/2;
    const double aux3 = Beta(aux2);
    lambda = kWAVE_MIN - kWAVE_STEP/2;

    for (int w = 0; w < kN_WAVE_BIN; ++w) {
      lambda += kWAVE_STEP;
      const double aux4 = lambda - kWAVE_STEP/2;
      const double aux5 = lambda + kWAVE_STEP/2;

      nphot[w] = CKVPhot(aux1, nIndex, aux3, dl, aux4, aux5) + ophot[w];
      ophot[w] = nphot[w];
    }

    jmax = j;
  }

  // Add Ultrarelativistic Energy Contribution "beta = 1.0"
  lambda = kWAVE_MIN - kWAVE_STEP/2;
  const double aux1 = TrackLengthN(jmax, chi, chimax);
  const double aux3 = 1;

  for (int w = 0; w < kN_WAVE_BIN; ++w) {
    lambda += kWAVE_STEP;
    const double aux4 = lambda - kWAVE_STEP/2;
    const double aux5 = lambda + kWAVE_STEP/2;

    nphot[w] = CKVPhot(aux1, nIndex, aux3, dl, aux4, aux5) + ophot[w];
  }

  // Area of Annulus on Ground
  const double dA = 2*utl::kPi * rp*rp * bin;

  for (int w = 0; w < kN_WAVE_BIN; ++w)
    dirCh[w] = nphot[w] * norm / dA;
}


// Cherenkov Beam
void
CherenkovModel::CKVBeam(const double h1, const double h2, const double theta, const double chimax, std::vector<double>& beamCh)
  const
{
  std::vector<double> nphot(kN_WAVE_BIN);
  std::vector<double> ophot(kN_WAVE_BIN);
  double chimed;
  double em;
  double lambda;

  for (int w = 0; w < kN_WAVE_BIN; ++w) {
    //nphot[w] = 0;
    //ophot[w] = 0;
    beamCh[w] = 0;
  }

  const double chi1 = det::Detector::GetInstance().GetAtmosphere().EvaluateDepth(h1);
  const double chi2 = det::Detector::GetInstance().GetAtmosphere().EvaluateDepth(h2);

  const double dl = std::abs(h2 - h1) / std::cos(theta);

  const double nIndex = det::Detector::GetInstance().GetAtmosphere().EvaluateRefractionIndex((h1 + h2) / 2);

  const double em = min(CKVThreshold(n_index), 2e3);

  // Integrate in Energy
  const double dE = 10;

  // From Threshold to Ultrarelativistic "E = 361.33 <---> beta = 1 - 1e-6"
  double jmax = 361.33;

  const double aux1 = (chi1 + chi2) / 2;

  for (double j = em; j <= jmax; j += dE) {
    const double tl1 = TrackLengthN(j, aux1, chimax);
    const double tl2 = TrackLengthN(j+dE, aux1, chimax);

    const double aux2 = tl1 - tl2;
    const double aux3 = j + dE/2;
    const double aux4 = Beta(aux3);
    lambda = kWAVE_MIN - kWAVE_STEP/2;

    for (int w = 0; w < kN_WAVE_BIN; ++w) {
      lambda += kWAVE_STEP;
      const double aux5 = lambda - kWAVE_STEP/2;
      const double aux6 = lambda + kWAVE_STEP/2;

      nphot[w] = CKVPhot(aux2, nIndex, aux4, dl, aux5, aux6) + ophot[w];
      ophot[w] = nphot[w];
    }

    jmax = j;
  }

  // Add Ultrarelativistic Energy Contribution "beta = 1.0"
  lambda = kWAVE_MIN - kWAVE_STEP/2;
  const double aux2 = TrackLengthN(jmax, aux1, chimax);
  const double aux4 = 1;

  for (int w = 0; w < kN_WAVE_BIN; ++w) {
    lambda += kWAVE_STEP;
    const double aux4 = lambda - kWAVE_STEP/2;
    const double aux5 = lambda + kWAVE_STEP/2;

    nphot[w] = CKVPhot(aux2, nIndex, aux4, dl, aux5, aux6) + ophot[w];
  }

  for (int w = 0; w < kN_WAVE_BIN; ++w)
    beamCh[w] = nphot[w];
}


// Cherenkov Rayleigh
void
CherenkovModel::CKVRayleigh(const double h1, const double h2, const double theta, const double rr, std::vector<double>& rayCh)
  const
{
  std::vector<double> lambda(kN_WAVE_BIN);

  const double aux = kWAVE_MIN - kWAVE_STEP/2;

  for (int w = 0; w < kN_WAVE_BIN; ++w)
    lambda[w] += aux + kWAVE_STEP;

  const CoordinateSystemPtr cs = det::Detector::GetInstance().GetSiteCoordinateSystem();

  const Point p1(0, 0, h1, cs);
  const Point p2(0, 0, h2, cs);

  AttenuationResult AttRes =
    det::Detector::GetInstance().GetAtmosphere().EvaluateRayleighAttenuation(p1, p2, lambda);

  const vector<double> attenuationVector = AttRes.GetTransmissionFactor();

  const double anfc = 1 + std::pow(std::cos(theta), 2);
  const double factor = 3 / (16*utl::kPi) * anfc / (rr*rr);
  for (int wl = 0; wl < kN_WAVE_BIN; ++wl) {
    const double probScat = 1 - attenuationVector[wl];
    rayCh[wl] = factor * probScat;
  }
}


// Cherenkov Aerosol
void
CherenkovModel::CKVAerosol(const double h1, const double h2, const double theta, const double rr, std::vector<double>& aerCh)
  const
{
  const double a = 0.8;
  const double thM = 26.7*utl::degree;

  std::vector<double> lambda(kN_WAVE_BIN);

  const double aux = kWAVE_MIN - kWAVE_STEP/2;

  for (int w = 0; w < kN_WAVE_BIN; ++w)
    lambda[w] += aux + kWAVE_STEP;

  const CoordinateSystemPtr cs = det::Detector::GetInstance().GetSiteCoordinateSystem();

  const Point p1(0, 0, h1, cs);
  const Point p2(0, 0, h2, cs);

  AttenuationResult AttRes =
    det::Detector::GetInstance().GetAtmosphere().EvaluateMieAttenuation(p1, p2, lambda);

  const vector<double> attenuationVector = AttRes.GetTransmissionFactor();

  const double factor = a * std::exp(-theta/thM) / (rr*rr);
  for (int wl = 0; wl < kN_WAVE_BIN; ++wl) {
    const double prob_scat = 1 - attenuationVector[wl];
    aerCh[wl] = factor * prob_scat;
  }
}