ProfileCalculator.cc

Go to the documentation of this file.

#include "ProfileCalculator.h"
#include "CFMatrixCalculator.h"
#include "LinearAlgebra.h"
#include "TelescopeData.h"


#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalFunctions.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeRecData.h>
#include <evt/ShowerFRecData.h>

#include <set>

using namespace fevt;
using namespace evt;
using namespace FdEnergyDepositFinderKG;
using namespace std;
using namespace utl;


bool
ProfileCalculator::CalculateProfile(Eye& eye,
                                    const CFMatrixCalculator& matrixBuilder,
                                    bool onlyEyedEdX)
{
  InitProfiles(eye);

  const unsigned int matrixSize =
    matrixBuilder.GetDirectFluorescenceMatrix().GetSize();
  ColumnVector lightFlux(matrixSize);
  vector<double> lightFluxVariance(matrixSize);
  vector<double> fluxFitVariance(matrixSize);
  vector<double> lightFluxTime(matrixSize);
  vector<double> lightFluxTimeBin(matrixSize);
  vector<double> gainVarVec(matrixSize);
  vector<double> bgVarVec(matrixSize);
  vector<double> fluxToPhotonVec(matrixSize);
  vector<double> depthVec(matrixSize);

  using namespace FdConstants;
  EyeRecData& eyeRecData = eye.GetRecData();
  ShowerFRecData& shower = eyeRecData.GetFRecShower();
  TabulatedFunctionErrors& dedxProfile = shower.GetEnergyDeposit();
  TabulatedFunctionErrors& sizeProfile = shower.GetLongitudinalProfile();
  TabulatedFunctionErrors& dirCher =
    eyeRecData.GetLightFlux(eCherDirect);
  TabulatedFunctionErrors& scatMCher =
    eyeRecData.GetLightFlux(eCherMieScattered);
  TabulatedFunctionErrors& scatRCher =
    eyeRecData.GetLightFlux(eCherRayleighScattered);
  TabulatedFunctionErrors& cherMS =
    eyeRecData.GetLightFlux(eCherMultScattered);
  TabulatedFunctionErrors& fluorTot =
    eyeRecData.GetLightFlux(eFluorTotal);
  TabulatedFunctionErrors& fluorDir =
    eyeRecData.GetLightFlux(eFluorDirect);
  TabulatedFunctionErrors& fluorMS =
    eyeRecData.GetLightFlux(eFluorMultScattered);
  // overwrite total light!
  TabulatedFunctionErrors& totalLight =
    eyeRecData.GetLightFlux(eTotal);
  totalLight.Clear();

  set<unsigned int> inFOVbins;
  unsigned int i = 0;
  for (CFMatrixCalculator::ConstTelDataIterator telIter = matrixBuilder.AllTelDataBegin();
        telIter != matrixBuilder.AllTelDataEnd(); ++telIter) {

    for (vector<TelescopeDataBin>::const_iterator binIter =
            telIter->TelDataBinsBegin();
          binIter !=  telIter->TelDataBinsEnd(); ++binIter) {

      if (telIter->GetType() == TelescopeData::eInsideFOV)
        inFOVbins.insert(i);

      lightFlux(i) = binIter->fSignal;
      lightFluxVariance[i] = pow(binIter->fSignalErr,2);
      lightFluxTime[i] = binIter->fTime;
      lightFluxTimeBin[i] = binIter->fTimeBinHalfWidth;
      depthVec[i] = binIter->GetMeanDepth();

      bgVarVec[i] = binIter->fBackgroundVariance;
      gainVarVec[i] = telIter->GetGainVariance();
      fluxToPhotonVec[i] = telIter->GetDiaphragmArea() *
                           telIter->GetPhotonToPhotoElectron();

      ++i;
    }
  }

  //bool onlyDirect=matrixBuilder.GetOnlyDirect();
  try {
    LowerTriangularMatrix inverseCFM = matrixBuilder.GetInverseCFMatrix();
    //ColumnVector dEdXVec = onlyDirect
    //  ? matrixBuilder.GetDirectCFMatrix().GetInverse()*lightFlux //these are two different operators
    //  : matrixBuilder.GetCFMatrix().GetInverse*lightFlux;
    ColumnVector dEdXVec = inverseCFM * lightFlux;

    // for dEdX/electron
    double xMax = shower.HasGHParameters() ?
      shower.GetGHParameters().GetXMax() :
      750.*g/cm2;

    const double electronEnergyThreshold = shower.GetEnergyCutoff();

    ColumnVector ghdEdXVec;
    // Calculate dedx vector from GH fit if available
    if (shower.HasGHParameters()) {
      ghdEdXVec.SetSize(matrixSize);
      for (unsigned int i = 0; i < matrixSize; ++i)
        ghdEdXVec(i) = shower.GetGHParameters().Eval(depthVec[i]);
    }

    // use GH-dEdXVec for the light fluxes if available
    ColumnVector& dEdXVecForLightFlux = (ghdEdXVec.GetSize() ? ghdEdXVec : dEdXVec);

    ColumnVector cherenkovFlux =
      matrixBuilder.GetDirectCherenkovMatrix()*dEdXVecForLightFlux;
    ColumnVector fluorescenceFlux =
      matrixBuilder.GetDirectFluorescenceMatrix()*dEdXVecForLightFlux;
    ColumnVector rayScatteredCherenkovFlux =
      matrixBuilder.GetRayScatteredCherenkovMatrix()*dEdXVecForLightFlux;
    ColumnVector mieScatteredCherenkovFlux =
      matrixBuilder.GetMieScatteredCherenkovMatrix()*dEdXVecForLightFlux;

    // recalculate dE/dX using forward folded GH outside FOV
    if (shower.HasGHParameters()) {
      for (unsigned int i = 0; i < matrixSize; ++i) {
        if (inFOVbins.find(i) == inFOVbins.end()) {
          const double lightSum = cherenkovFlux(i) +
                                  fluorescenceFlux(i) +
                                  rayScatteredCherenkovFlux(i) +
                                  mieScatteredCherenkovFlux(i);
          lightFlux(i) = lightSum;
        }
      }
      dEdXVec = inverseCFM*lightFlux;
    }

    // calculate lightflux errors based on fit result
    for (unsigned int i = 0; i < matrixSize; ++i) {
      if (inFOVbins.find(i) != inFOVbins.end()) {
        const double kFluxToPe = fluxToPhotonVec[i];
        const double kFluxToPe2 = kFluxToPe * kFluxToPe;
        const double gainVariance = gainVarVec[i];
        const double totalFitLight = cherenkovFlux(i) +
          fluorescenceFlux(i) +
          rayScatteredCherenkovFlux(i) +
          mieScatteredCherenkovFlux(i);
        const double bgVariance = bgVarVec[i];
        const double nPeBackground =
          bgVariance*kFluxToPe2 / (1 + gainVariance);
        // in first iteration (no GH fit) negative light fluxes
        // may occur
        const double nPeSignal = fmax(0, totalFitLight*kFluxToPe);
        const double nPeExp = nPeSignal + nPeBackground;
        const double fluxVariance = nPeExp*(1 + gainVariance) / kFluxToPe2;
        fluxFitVariance[i] = fluxVariance;
      } else
        fluxFitVariance[i] = 0;
    }

    // Fill eye data structures
    for (unsigned int i = 0; i < matrixSize; ++i) {

      if (inFOVbins.find(i) != inFOVbins.end()) {

        double varianceSum = 0.;
        for (unsigned int k = 0; k <= i; ++k)
          varianceSum += inverseCFM(i,k)*inverseCFM(i,k)*fluxFitVariance[k];

        dedxProfile.PushBack(depthVec[i], 0., dEdXVec(i), sqrt(varianceSum));
        if (onlyEyedEdX)
          continue;
        dirCher.PushBack(lightFluxTime[i], lightFluxTimeBin[i],
                         cherenkovFlux(i), 0.);
        scatMCher.PushBack(lightFluxTime[i], lightFluxTimeBin[i],
                           mieScatteredCherenkovFlux(i), 0.);
        scatRCher.PushBack(lightFluxTime[i], lightFluxTimeBin[i],
                           rayScatteredCherenkovFlux(i), 0.);
        fluorTot.PushBack(lightFluxTime[i], lightFluxTimeBin[i],
                          fluorescenceFlux(i), 0);
        fluorDir.PushBack(lightFluxTime[i], lightFluxTimeBin[i],
                          fluorescenceFlux(i), 0);
        totalLight.PushBack(lightFluxTime[i], lightFluxTimeBin[i],
                            lightFlux(i), sqrt(lightFluxVariance[i]));

        const double alpha = utl::EnergyDeposit(depthVec[i],
                                                xMax,
                                                electronEnergyThreshold);
        sizeProfile.PushBack(depthVec[i],0.,dEdXVec(i)/alpha,
                             sqrt(varianceSum)/alpha);

        // FIXME fill this stuff
        cherMS.PushBack(0,0,0,0);
        fluorMS.PushBack(0,0,0,0);
      }
    }

    // Fill telescope data structures
    int iMatrixBin = 0;

    for (CFMatrixCalculator::ConstTelDataIterator telIter = matrixBuilder.AllTelDataBegin();
         telIter != matrixBuilder.AllTelDataEnd(); ++telIter) {

      if (telIter->GetType() == TelescopeData::eInsideFOV) {

        Telescope& tel = eye.GetTelescope(telIter->GetId());
        if (!tel.HasRecData()) // shouldn't be necessary
          tel.MakeRecData();
        TelescopeRecData& telRecData = tel.GetRecData();

        // Fetch the telescope light fluxes
        TabulatedFunctionErrors& telDirCher =
          telRecData.GetLightFlux(eCherDirect);
        TabulatedFunctionErrors& telScatMCher =
          telRecData.GetLightFlux(eCherMieScattered);
        TabulatedFunctionErrors& telScatRCher =
          telRecData.GetLightFlux(eCherRayleighScattered);
        TabulatedFunctionErrors& telFluorTot =
          telRecData.GetLightFlux(eFluorTotal);
        TabulatedFunctionErrors& telFluorDir =
          telRecData.GetLightFlux(eFluorDirect);

        //TabulatedFunctionErrors& telCherMS =
        //  telRecData.GetLightFlux(eCherMultScattered);
        //TabulatedFunctionErrors& telFluorMS =
        //  telRecData.GetLightFlux(eFluorMultScattered);

        for (vector<TelescopeDataBin>::const_iterator binIter = telIter->TelDataBinsBegin();
              binIter !=  telIter->TelDataBinsEnd(); ++binIter) {
          if (inFOVbins.find(iMatrixBin) != inFOVbins.end()) {
            telDirCher.PushBack(lightFluxTime[iMatrixBin],
                                lightFluxTimeBin[iMatrixBin],
                                cherenkovFlux(iMatrixBin), 0.);
            telScatMCher.PushBack(lightFluxTime[iMatrixBin],
                                  lightFluxTimeBin[iMatrixBin],
                                  mieScatteredCherenkovFlux(iMatrixBin), 0.);
            telScatRCher.PushBack(lightFluxTime[iMatrixBin],
                                  lightFluxTimeBin[iMatrixBin],
                                  rayScatteredCherenkovFlux(iMatrixBin), 0.);
            telFluorTot.PushBack(lightFluxTime[iMatrixBin],
                                 lightFluxTimeBin[iMatrixBin],
                                 fluorescenceFlux(iMatrixBin), 0);
            telFluorDir.PushBack(lightFluxTime[iMatrixBin],
                                 lightFluxTimeBin[iMatrixBin],
                                 fluorescenceFlux(iMatrixBin), 0);
            //telTotalLight.PushBack(lightFluxTime[iMatrixBin],
            // lightFluxTimeBin[iMatrixBin],
            //                       lightFlux(iMatrixBin),
            // sqrt(lightFluxVariance[iMatrixBin]));
          }

          ++iMatrixBin;
        } // end foreach TelescopeDataBin
      } else
        iMatrixBin += telIter->GetTelescopeDataBins().size();

    } // end foreach TelescopeData

  } catch (SingularMatrixException& b) {
    ERROR(b.GetMessage());
    return false;
  }
  return true;
}


void
ProfileCalculator::InitProfiles(Eye& eye)
{
  evt::ShowerFRecData& showerFRecData=eye.GetRecData().GetFRecShower();

  // create longitudinal profile if necessary
  if (!showerFRecData.HasLongitudinalProfile())
    showerFRecData.MakeLongitudinalProfile();
  showerFRecData.GetLongitudinalProfile().Clear();

  // create energy deposit profile if necessary
  if (!showerFRecData.HasEnergyDeposit())
    showerFRecData.MakeEnergyDeposit();
  showerFRecData.GetEnergyDeposit().Clear();

  // create fluorescence photons at track if necessary
  if (!showerFRecData.HasFluorescencePhotons())
    showerFRecData.MakeFluorescencePhotons();
  showerFRecData.GetFluorescencePhotons().Clear();

  // --- eye light fluxes
  EyeRecData& eyeRecData = eye.GetRecData();
  using namespace FdConstants;

  // direct cherenkov light at aperture
  if(!eyeRecData.HasLightFlux(eCherDirect))
    eyeRecData.MakeLightFlux(eCherDirect);
  eyeRecData.GetLightFlux(eCherDirect).Clear();

  // Rayleigh scattered cherenkov light at aperture
  if(!eyeRecData.HasLightFlux(eCherRayleighScattered))
    eyeRecData.MakeLightFlux(eCherRayleighScattered);
  eyeRecData.GetLightFlux(eCherRayleighScattered).Clear();

  // Mie scattered cherenkov light at aperture
  if(!eyeRecData.HasLightFlux(eCherMieScattered))
    eyeRecData.MakeLightFlux(eCherMieScattered);
  eyeRecData.GetLightFlux(eCherMieScattered).Clear();

  // direct and multiple scattered fluorescence light at aperture
  if(!eyeRecData.HasLightFlux(eFluorMultScattered))
    eyeRecData.MakeLightFlux(eFluorMultScattered);
  eyeRecData.GetLightFlux(eFluorMultScattered).Clear();

  if(!eyeRecData.HasLightFlux(eFluorDirect))
    eyeRecData.MakeLightFlux(eFluorDirect);
  eyeRecData.GetLightFlux(eFluorDirect).Clear();

  if(!eyeRecData.HasLightFlux(eFluorTotal))
    eyeRecData.MakeLightFlux(eFluorTotal);
  eyeRecData.GetLightFlux(eFluorTotal).Clear();

  // multiple scattered cherenkov light at aperture
  if(!eyeRecData.HasLightFlux(eCherMultScattered))
    eyeRecData.MakeLightFlux(eCherMultScattered);
  eyeRecData.GetLightFlux(eCherMultScattered).Clear();


  // ---- Now, we also initialize the light fluxes for each of the telescopes

  // We skip all telescopes without TelescopeRecData because they must have
  // had a total light at aperture at the very least before we started.
  for (Eye::TelescopeIterator telIt = eye.TelescopesBegin(ComponentSelector::eHasData),
       end = eye.TelescopesEnd(ComponentSelector::eHasData); telIt != end; ++telIt)
  {
    if (!telIt->HasRecData())
      continue;
    TelescopeRecData& telRecData = telIt->GetRecData();

    // direct cherenkov light at aperture
    if(!telRecData.HasLightFlux(eCherDirect))
      telRecData.MakeLightFlux(eCherDirect);
    else
      telRecData.GetLightFlux(eCherDirect).Clear();

    // Rayleigh scattered cherenkov light at aperture
    if(!telRecData.HasLightFlux(eCherRayleighScattered))
      telRecData.MakeLightFlux(eCherRayleighScattered);
    else
      telRecData.GetLightFlux(eCherRayleighScattered).Clear();

    // Mie scattered cherenkov light at aperture
    if(!telRecData.HasLightFlux(eCherMieScattered))
      telRecData.MakeLightFlux(eCherMieScattered);
    else
      telRecData.GetLightFlux(eCherMieScattered).Clear();

    // direct and multiple scattered fluorescence light at aperture
    if(!telRecData.HasLightFlux(eFluorMultScattered))
      telRecData.MakeLightFlux(eFluorMultScattered);
    else
      telRecData.GetLightFlux(eFluorMultScattered).Clear();

    if(!telRecData.HasLightFlux(eFluorDirect))
      telRecData.MakeLightFlux(eFluorDirect);
    else
      telRecData.GetLightFlux(eFluorDirect).Clear();

    if(!telRecData.HasLightFlux(eFluorTotal))
      telRecData.MakeLightFlux(eFluorTotal);
    else
      telRecData.GetLightFlux(eFluorTotal).Clear();

    // multiple scattered cherenkov light at aperture
    if(!telRecData.HasLightFlux(eCherMultScattered))
      telRecData.MakeLightFlux(eCherMultScattered);
    else
      telRecData.GetLightFlux(eCherMultScattered).Clear();

  } // end foreach telescope

}