LaserLightSimulator.cc

Go to the documentation of this file.

#include "LaserLightSimulator.h"

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

#include <evt/Event.h>
#include <evt/ShowerSimData.h>

#include <evt/LaserData.h>

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>

#include <atm/AttenuationResult.h>
#include <atm/ProfileResult.h>
#include <atm/ScatteringResult.h>
#include <atm/Atmosphere.h>

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/UTMPoint.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/Photon.h>
#include <utl/RandomEngine.h>

#include <cmath>
#include <sstream>
#include <vector>

#include <CLHEP/Random/RandGauss.h>

using namespace LaserLightSimulatorNA;
using CLHEP::RandGauss;

using namespace fwk;
using namespace evt;
using namespace det;
using namespace fdet;
using namespace atm;
using namespace utl;

using namespace std;

LaserLightSimulator::LaserLightSimulator() :
  fMaxLaserHeight(0),
  fEnergy(0),
  fNPhotons(0),
  fDoPoissonFluctuate(false),
  fRandomEngine(0)
{
}

LaserLightSimulator::~LaserLightSimulator()
{
}

VModule::ResultFlag
LaserLightSimulator::Init()
{
  CentralConfig* cc = CentralConfig::GetInstance();

  Branch topBranch =
    cc->GetTopBranch("LaserLightSimulatorNA");

  topBranch.GetChild("maxLaserHeight").GetData(fMaxLaserHeight);
  topBranch.GetChild("poissonFluctuate").GetData(fDoPoissonFluctuate);
  topBranch.GetChild("gpsTimingFluctuation").GetData(fGpsTimingFluctuation);

  fRandomEngine =
    &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

  ostringstream info;
  info << "Version: " << GetVersionInfo(VModule::eRevisionNumber)
       << "\n Parameters:\n"
       << "                   Maximum laser height: "
       << fMaxLaserHeight / km << " km\n"
       << "     Poisson-fluctuate the photon count: "
       << (fDoPoissonFluctuate ? "yes\n" : "no\n")
       << "                 GPS timing fluctuation: "
       << fGpsTimingFluctuation / nanosecond << " ns\n";
  INFO(info);

  return eSuccess;
}

VModule::ResultFlag
LaserLightSimulator::Run(evt::Event& theEvent)
{
  INFO("Calculating laser track");

  if (!theEvent.HasSimShower()){
    ERROR ("Event has no simulated shower.");
    return eFailure;
  }
  evt::ShowerSimData& theShower = theEvent.GetSimShower();

  if (!theShower.HasLaserData()){
    ERROR ("Event has no simulated laser.");
    return eFailure;
  }
  LaserData& theLaser = theShower.GetLaserData();

  fEnergy = theShower.GetEnergy();
  fLaserWavelength.push_back(theLaser.GetLaserWavelength());

  fNPhotons = (fEnergy/joule) / (kSpeedOfLightSI*kPlanckSI/fLaserWavelength[0]);
  ostringstream info;
  info << " Laser energy = " << fEnergy / (1e-3*joule) << " mJ;"
       << " wavelength = " << fLaserWavelength[0] / nanometer << " nm;"
       << " <N> = " << fNPhotons;

  if (fDoPoissonFluctuate) {
    const double N = RandGauss::shoot(&fRandomEngine->GetEngine(),
                                      fNPhotons,
                                      sqrt(fNPhotons));
    info << "; N - <N> = " << N - fNPhotons;
    fNPhotons = N;
  }

  double gpsTimingOffset = 0;
  if(fGpsTimingFluctuation != 0)
    gpsTimingOffset = RandGauss::shoot(&fRandomEngine->GetEngine(),
                                       gpsTimingOffset,
                                       kSpeedOfLight * fGpsTimingFluctuation);
  info << "; Laser offset = " << gpsTimingOffset / m;

  INFO(info);

  utl::Point posLsr = theShower.GetPosition();
  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
  CoordinateSystemPtr coreCS =
    AugerCoordinateSystem::Create(posLsr, wgs84, wgs84.GetECEF());

  Vector vecLsr = theShower.GetDirection();

  posLsr = posLsr + gpsTimingOffset * vecLsr;

  const CoordinateSystemPtr localCS = theEvent.GetSimShower().GetLocalCoordinateSystem();
  double zenith = (-theEvent.GetSimShower().GetDirection()).GetTheta(localCS);
  if (zenith > (0.5*kPi)) {
    zenith = kPi - zenith;
    INFO("UPGOING ->  Inverting Zenith Angle for CosTheta Calculation");
  }
  double cosTheta = cos(zenith);

  double zLaserGround = (wgs84.PointToLatitudeLongitudeHeight(posLsr)).get<2>();

  const double fluorLightBin = 50.*ns;
  double delta  = kSpeedOfLight * fluorLightBin;
  int nPoints   = int(abs((zLaserGround - fMaxLaserHeight)/cosTheta/delta));

  const int nWaves = 1;
  vector<double> distance;
  vector<double> laserPhotons[nWaves];
  double photons = fNPhotons;

  fLsrInitialPos = posLsr;
  fLsrFinalPos   = posLsr + (fMaxLaserHeight - zLaserGround)/cosTheta*vecLsr;

  // first point
  utl::Point LsrPoint = fLsrInitialPos + delta * 0.5 * vecLsr;
  photons = EvaluatePhotons(photons, fLsrInitialPos, LsrPoint);
  (laserPhotons[0]).push_back(photons);
  double height_point = (wgs84.PointToLatitudeLongitudeHeight(LsrPoint)).get<2>();
  double distance_point = (height_point - zLaserGround) / cosTheta;
  distance.push_back(distance_point);

  info.str("");
  info << " first point "
       << (wgs84.PointToLatitudeLongitudeHeight(fLsrInitialPos)).get<2>()
       << "; last point "
       << (wgs84.PointToLatitudeLongitudeHeight(fLsrFinalPos)).get<2>();
  INFO(info);

  for (int i=1; i<nPoints; i++){
    utl::Point LsrPoint2 = LsrPoint +  delta *vecLsr;
    height_point  =
      (wgs84.PointToLatitudeLongitudeHeight(LsrPoint2)).get<2>();

    distance_point = (height_point-zLaserGround)/cosTheta;
    distance.push_back(distance_point);

    photons = EvaluatePhotons(photons, LsrPoint, LsrPoint2);
    (laserPhotons[0]).push_back(photons);

    LsrPoint = LsrPoint2;
  }

  for (int iwl=0; iwl<nWaves; iwl++) {
    vector<double> laserLight = laserPhotons[iwl];

    TabulatedFunction f = TabulatedFunction(distance, laserLight);
    if (!theEvent.GetSimShower().HasFluorescencePhotons(iwl)) {
      theEvent.GetSimShower().AddFluorescencePhotons(f, iwl);
    }
  }

  return eSuccess;
}

VModule::ResultFlag
LaserLightSimulator::Finish()
{
  return eSuccess;
}

double
LaserLightSimulator::EvaluatePhotons(const double Nph,
  const utl::Point& LsrPos1, const utl::Point& LsrPos2)
  const
{
  const Atmosphere& atmos = Detector::GetInstance().GetAtmosphere();

  // Estimate how many photons in the laser beam survive molecular and
  // aerosol attenuation
  AttenuationResult rAttenLaserPoint =
    atmos.EvaluateRayleighAttenuation(LsrPos1, LsrPos2, fLaserWavelength);
  utl::TabulatedFunctionErrors rAtt = rAttenLaserPoint.GetTransmissionFactor();

  AttenuationResult mAttenLaserPoint =
    atmos.EvaluateMieAttenuation(LsrPos1, LsrPos2, fLaserWavelength);
  utl::TabulatedFunctionErrors mAtt = mAttenLaserPoint.GetTransmissionFactor();

  const double att = rAtt[0].Y() * mAtt[0].Y();

  return Nph * att;
}