MultipleScatterer.cc

Go to the documentation of this file.

#include "MultipleScatterer.h"

#include <utl/Reader.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/AugerUnits.h>
#include <utl/RandomEngine.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>

#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>

#include <fevt/FdConstants.h>

#include <atm/Atmosphere.h>
#include <det/Detector.h>

#include <vector>

#include <CLHEP/Random/Randomize.h>
using CLHEP::RandFlat;

using namespace utl;

MultipleScatterer::MultipleScatterer()
{
  fwk::CentralConfig* cc = fwk::CentralConfig::GetInstance();
  utl::Branch topB = cc->GetTopBranch("MultipleScatterer");
  topB.GetChild("MaxTime").GetData(MaxTime);
  topB.GetChild("MaxZeta").GetData(MaxZeta);
  const string modelString = topB.GetChild("fAerosolModel").Get<string>();
  if (modelString == "eLongtin")
    fAerosolModel = eLongtin;
  else if (modelString == "eHG01")
    fAerosolModel = eHG01;
  else if (modelString == "eHG05")
    fAerosolModel = eHG05;
  else if (modelString == "eHG09")
    fAerosolModel = eHG09;
  else {
    ostringstream err;
    err << " Unknown aerosol model: " << modelString
        << ", switching to default eHG05";
    ERROR(err);
    fAerosolModel = eHG05;
  }
  fRandomEngine = &fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);
}


void
MultipleScatterer::AddPhotons(const utl::Point& PointOfOrigin,
                              const utl::Photon& PhotonDirect,
                              const utl::Point& PointOfDetection,
                              std::vector<utl::Photon>& photonVector)
{
  const double wavelength = PhotonDirect.GetWavelength();
  const atm::Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  const double mieAtt = atmo.EvaluateMieAttenuation(PointOfOrigin,PointOfDetection,wavelength);
  const double rayAtt = atmo.EvaluateRayleighAttenuation(PointOfOrigin,PointOfDetection,wavelength);

  const double Transmission = rayAtt*mieAtt;

  AddPhotons(PointOfOrigin, PhotonDirect, Transmission, photonVector);
}



void
MultipleScatterer::AddPhotons(const utl::Point& PointOfOrigin,
                              const utl::Photon& PhotonDirect,
                              const double Transmission,
                              std::vector<utl::Photon>& photonVector)
{

  if (PhotonDirect.GetSource() == fevt::FdConstants::eFluorTotal) {
    INFO(" Photon source eFluorTotal ");
    return;
  }
  if (PhotonDirect.GetSource() != fevt::FdConstants::eFluorDirect)
     return;

  const utl::ReferenceEllipsoid& wgs84 = utl::ReferenceEllipsoid::GetWGS84();
  double height1 = wgs84.PointToLatitudeLongitudeHeight(PointOfOrigin).get<2>();
  double height = height1 / kilometer;
  double OptDistance = - log (Transmission);

  if ((height > 0.) && (height < 15.) && (OptDistance < 5.)) {
                               //limits of application of parameterization

    double TimeBin = 100 * nanosecond;

    double zeta;
    double weightfactor;

    for (double Time = (0 * nanosecond); Time < MaxTime; Time+=TimeBin) {
       if (GetRandomZeta(zeta, weightfactor, Time, OptDistance, height, MaxZeta)) {


        Photon newPhoton = PhotonDirect;
        const Vector& photonDir = PhotonDirect.GetDirection();
        double timeScatt = PhotonDirect.GetTime().GetInterval();

        if (weightfactor > 1.) {
          double weight = PhotonDirect.GetWeight();
          newPhoton.SetWeight(weight * weightfactor);
        }

        const CoordinateSystemPtr& PhotonDirCS = photonDir.GetCoordinateSystem();

        // new coordinate system, photonDir along Z axis
        double theta = photonDir.GetTheta(PhotonDirCS);
        double phi = photonDir.GetPhi(PhotonDirCS);

        CoordinateSystemPtr newCS(
            CoordinateSystem::RotationY(theta,
            CoordinateSystem::RotationZ(phi,PhotonDirCS)));

        // random change of photon direction
        double thetaRandom = zeta;
        double phiRandom = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 2*kPi);

        Vector newphotonDir(cos(phiRandom)*sin(thetaRandom),
            sin(phiRandom)*sin(thetaRandom),cos(thetaRandom),newCS);

        newPhoton.SetDirection(newphotonDir);

        // random change of time within given time bin
        double timeRandom = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 100.);
        timeScatt += (Time + timeRandom * nanosecond);

        newPhoton.SetTime(TimeInterval(timeScatt));
        newPhoton.SetSource(fevt::FdConstants::eFluorMultScattered);

        // make photon
        photonVector.push_back(newPhoton);
      }

    }
  }
}

bool
MultipleScatterer::GetRandomZeta(double& zeta, double& weightfactor,
                                 double Time, double OptDistance,
                                 double height, double MaxZeta)
{
  double A(0), B(0), C(0), D(0), E(0), F(0), G(0), H(0);
  Time = Time / utl::ns;

  if (Time < 200){
    switch(fAerosolModel) {
    case eLongtin:
      A = (0.00152 + 0.0000165 * Time) * pow (OptDistance, (0.0826+0.00525*Time))
          * exp(-0.644*height) + 0.0000569 + 0.0000011 * Time;
      B = -0.0504 * height * OptDistance + (-0.17+0.000766*Time)*OptDistance
          - 0.0323 * height - 1.78 + 0.00814 * Time;
      C = 3.59e-6 * OptDistance + (2.03e-7 + 1.65e-9*Time)*height + 2.44e-6 - 6.2e-8 * Time;
      break;
    case eHG01:
      A = (0.000272 + 7.88e-6 * Time) * pow (OptDistance, (-0.248+0.0061*Time))
          * exp(-0.456*height) + 4.04e-5 + 1.29e-6 * Time;
      B = -0.0372 * height * OptDistance + (-0.272+0.000413*Time)*OptDistance
          - 0.00927 * height - 2.04 + 0.009 * Time;
      C = 1.61e-6 * OptDistance + (-1.56e-8 + 2.82e-9*Time)*height + 1.64e-6 - 5.09e-8 * Time;
      break;
    case eHG05:
      A = (0.00123 + 1.6e-5 * Time) * pow (OptDistance, (0.0346+0.00544*Time))
          * exp(-0.637*height) + 5.49e-5 + 1.15e-6 * Time;
      B = -0.0495 * height * OptDistance + (-0.186+0.000749*Time)*OptDistance
          - 0.0292 * height - 1.82 + 0.00831 * Time;
      C = 2.98e-6 * OptDistance + (1.41e-7 + 2.03e-9*Time)*height + 2.3e-6 - 5.62e-8 * Time;
      break;
    case eHG09:
      A = (0.0128 - 4.4e-5 * Time) * pow (OptDistance, (0.607+0.0026*Time))
          * exp(-0.43*height) + 0.000474 - 5.66e-6 * Time;
      B = (0.0343 - 0.000747 * Time) * height * OptDistance + (-0.178+0.0023*Time)*OptDistance
          + (-0.0132 - 0.000584 * Time) * height - 1.33 + 0.00665 * Time;
      C = (-4.34e-6 - 1.35e-6 * Time) * OptDistance + (4.73e-6 + 2.61e-7 * Time) * height
          - 1.18-5 - 8.11e-7 * Time;
      break;
    default:
      ERROR("No valid multiple satter model specified");
      break;
    }

    double ScattLight = IntegralA(A,B,C,MaxZeta);
    double ScattAngle = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 1.0);

    if (ScattLight>1) {
      weightfactor = ScattLight;
      ScattAngle *= weightfactor;
    }
    else
      weightfactor = 1.;


    if (ScattAngle > ScattLight)
      return false;
    else {
      double k = 0.*degree;
      double l = MaxZeta;
      double Accuracy = 0.1*degree;
      zeta = -1.*degree;
      int NMax = 20;

      // find angle by bisection of cumulative
      while (NMax>0) {
        double m = (k+l)/2;
        double val = IntegralA(A,B,C,m);
        if (val == ScattAngle || (l-k) < Accuracy) {
          zeta = m;
          return true;
        }
        if (val > ScattAngle)
          l = m;
        else
          k = m;
        --NMax;
      }
      WARNING(" maxIter reached!!! ");
      return false;
    }

  }
  else {
    switch(fAerosolModel) {
    case eLongtin:
      D = 128 * OptDistance * exp(-0.329 * height) * pow (Time, -1.76);
      E = (-81.3 * OptDistance + 21.8) * height * pow (Time, -0.857);
      F = 2004 * OptDistance * exp (-0.103*height) * pow (Time, -2.22);
      G = (-132 * height * OptDistance + 78.4 * height - 600) * pow (Time, -0.833);
      H = 0.0207 * exp(-0.367*OptDistance) * exp(-0.371*height) * pow(Time, (0.0812*height-1.54));
      break;
    case eHG01:
      D = 43 * OptDistance * exp(-0.098 * height) * pow (Time, -1.72);
      E = (-55.4 * OptDistance + 18) * height * pow (Time, -0.829);
      F = 2275 * OptDistance * exp (-0.0586*height) * pow (Time, -2.26);
      G = (-151 * height * OptDistance + 78.6 * height - 457) * pow (Time, -0.814);
      H = 22.1 * exp(-9.99*OptDistance) * exp(-2.99*height) * pow(Time, (0.0258*height-0.881));
      break;
    case eHG05:
      D = 90.7 * OptDistance * exp(-0.215 * height) * pow (Time, -1.77);
      E = (-74.6 * OptDistance + 19.8) * height * pow (Time, -0.863);
      F = 4260 * OptDistance * exp (-0.297*height) * pow (Time, -2.24);
      G = (-130 * height * OptDistance + 77.2 * height - 613) * pow (Time, -0.829);
      H = 0.0141 * exp(-0.286*OptDistance) * exp(-0.415*height) * pow(Time, (0.0748*height-1.44));
      break;
    case eHG09:
      D = 161 * OptDistance * exp(-0.068 * height) * pow (Time, -1.45);
      E = (-60.2 * OptDistance + 38.5) * height * pow (Time, -0.244);
      F = 934 * OptDistance * exp (-0.0437*height) * pow (Time, -2.16);
      G = (-49.5 * height * OptDistance + 7164 * height - 1340) * pow (Time, -0.824);
      H = -0.0854 * exp(0.263*OptDistance) * exp(-0.767*height) * pow(Time, (-1.11));
      break;
    default: break;
    }

    double ScattLight = IntegralB(D,E,F,G,H,MaxZeta);
    double ScattAngle = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 1.0);

    if (ScattLight>1) {
      weightfactor = ScattLight;
      ScattAngle *= weightfactor;
    }
    else
      weightfactor = 1.;


    if (ScattAngle > ScattLight)
      return false;
    else {
      double k = 0.*degree;
      double l = MaxZeta;
      double Accuracy = 0.1*degree;
      zeta = -1.*degree;
      int NMax = 20;

      // find angle by bisection of cumulative
      while (NMax>0) {
        double m = (k+l)/2;
        double val = IntegralB(D,E,F,G,H,m);
        if (val == ScattAngle || (l-k) < Accuracy) {
          zeta = m;
          return true;
        }
        if (val > ScattAngle)
          l = m;
        else
          k = m;
        --NMax;
      }
      WARNING(" maxIter reached!!! ");
      return false;
    }
  }
}

double
MultipleScatterer::IntegralA(double A, double B, double C, double zeta)
{
  double integr, integr0, integr1;
  double x = zeta/radian;
  double x0 = 0.1*degree/radian;
  A *= pow(57.295779513, B);

  if (x > x0) {
    // series expansion of the integral
    if (B == -2) {
      integr1 = log(x) - x*x/12 + pow(x, 4)/480;
      integr0 = log(x0) - x0*x0/12 + pow(x0, 4)/480;
    } else if (B == -4) {
      integr1 = -1/(2*x*x) - log(x)/6 + x*x/240;
      integr0 = -1/(2*x0*x0) - log(x0)/6 + x0*x0/240;
    } else if (B == -6) {
      integr1 = -1/(4*pow(x, 4)) + 1/(12*x*x) + log(x)/120;
      integr0 = -1/(4*pow(x0, 4)) + 1/(12*x0*x0) + log(x0)/120;
    } else {
      integr1 = pow(x, B)*(x*x/(2+B) - pow(x, 4)/(6*(4+B)) + pow(x, 6)/(120*(6+B)));
      integr0 = pow(x0, B)*(x0*x0/(2+B) - pow(x0, 4)/(6*(4+B)) + pow(x0, 6)/(120*(6+B)));
    }
    integr = (A * integr1 - C * cos(zeta)) - (A * integr0 - C * cos(x0*radian));
    integr += (A * pow(x0, B) + C) * 1.5230867e-6;      //normalized contribution from image center
    integr *= 20626.48062;                              //2*pi*180^2/pi^2
  } else {
    integr = (A * pow(x0, B) + C) * (1. - cos(zeta));
    integr *= 20626.48062;
  }

  return integr;
}

double
MultipleScatterer::IntegralB(double D, double E, double F,
                             double G, double H, double zeta)
{
  E *= 57.295779513;
  G *= 57.295779513;
  double integr;
  integr = -H*cos(zeta) + (D*exp(E * zeta/radian)*(-cos(zeta)+E*sin(zeta)))/(1+E*E)
        + (F*exp(G * zeta/radian)*(-cos(zeta)+G*sin(zeta)))/(1+G*G)
        - (-H - D/(1+E*E) - F/(1+G*G));
  integr *= 20626.48062;                                //2*pi*180^2/pi^2
  return integr;
}