CFMatrixCalculator.h

Go to the documentation of this file.

#ifndef _CFMatrixCalculator_h_
#define _CFMatrixCalculator_h_

#include "LinearAlgebra.h"
#include "TelescopeData.h"
#include "LateralLightCalculator.h"

#include <fevt/FdConstants.h>
#include <utl/Vector.h>

#include <vector>
#include <boost/iterator/filter_iterator.hpp>


namespace fevt {
  class Eye;
  class TelescopeRecData;
}

namespace utl {
  class TabulatedFunctionErrors;
  class TabulatedFunction;
  class Point;
}

namespace FdEnergyDepositFinderKG {

  class LateralLightCalculator;


  class CFMatrixCalculator {

  public:
    enum EMethod {
      eFast,
      ePrecise
    };

  public:
    CFMatrixCalculator() = default;
    ~CFMatrixCalculator();
    void Init();
    void BuildMatrix(const fevt::Eye& eye, const bool leavingAtmoIsError = false, const unsigned int addDense = 1);
    bool HasMatrix() const { return GetSize(); }
    unsigned int GetSize() const { return fOnlyDirect ? fDirectCFMatrix.GetSize() : fCFMatrix.GetSize(); }
    unsigned int GetFOVSize() const { return fFOVsize; }
    double GetXmax() const { return fXmax; }
    unsigned int GetEyeId() const { return fEyeId; }
    bool GetOnlyDirect() const { return fOnlyDirect; }

    const LowerTriangularMatrix& GetCFMatrix() const
    { assert(!fOnlyDirect); return fCFMatrix; }  //would have to break const or always needlessly create

    const DiagonalMatrix& GetDirectCFMatrix() const
    { assert(fOnlyDirect); return fDirectCFMatrix; }

    LowerTriangularMatrix
    GetInverseCFMatrix()
      const
    {
      if (!fOnlyDirect)
        return fCFMatrix.GetInverse();
      return LowerTriangularMatrix(fNumberOfDepthBins) + fDirectCFMatrix.GetInverse();
    }

    //TODO? avoid blowing up diagonal matrix in direct case
    //      by doing dedx = CFM^-1* lightFlux in one method
    //      -- would not look elegant though.
    //      alternatively create cases in CalculateProfile.
    //ColumnVector inverseCFMTimes(ColumnVector& lightFlux) const;

    const LowerTriangularMatrix& GetMieScatteredCherenkovMatrix() const
    { return fMieScatteredCherenkovMatrix; }
    const LowerTriangularMatrix& GetRayScatteredCherenkovMatrix() const
    { return fRayScatteredCherenkovMatrix; }
    const DiagonalMatrix& GetDirectFluorescenceMatrix() const
    { return fDirectFluorescenceMatrix;}
    const DiagonalMatrix& GetDirectCherenkovMatrix() const
    { return fDirectCherenkovMatrix;}

    typedef std::vector<TelescopeData>::const_iterator ConstTelDataIterator;
    ConstTelDataIterator AllTelDataBegin() const
    { return fTelescopeData.begin(); }
    ConstTelDataIterator AllTelDataEnd() const
    { return fTelescopeData.end(); }

    ConstTelDataIterator NoiseTelDataBegin() const
    { return fNoiseTelescopeData.begin(); }
    ConstTelDataIterator NoiseTelDataEnd() const
    { return fNoiseTelescopeData.end(); }

    typedef boost::filter_iterator<TelescopeFilter, ConstTelDataIterator> FilteredTelDataIterator;
    FilteredTelDataIterator TelDataBegin(const TelescopeData::ETelDataType type)
    { return FilteredTelDataIterator(TelescopeFilter(type), AllTelDataBegin(), AllTelDataEnd()); }
    FilteredTelDataIterator TelDataEnd(const TelescopeData::ETelDataType type)
    { return FilteredTelDataIterator(TelescopeFilter(type), AllTelDataEnd(), AllTelDataEnd()); }

    void SetVerbosity(const int v) { fVerbosity = v; }
    void SetMethod(const CFMatrixCalculator::EMethod m) { fMethod = m; }
    void SetLDFMethod(const LateralLightCalculator::ECalculationMethod m) { fLateralLightCalculator->SetMethod(m); }
    void SetOnlyDirect(const bool od) { fOnlyDirect = od;}

  private:
    void Clear();
    bool CalculateTelescopeData(const fevt::Eye& eye, const bool leavingAtmoIsError, const unsigned int addDense);

    bool InitCalculation(const fevt::Eye& eye);
    void AddBinsOutsideFOV(const fevt::Eye& eye);

    void CalculateDiagonalParameters();
    void CalculateFluorescenceMatrix();
    void CalculateDirectCherenkovMatrix();
    void CalculateMieAndRayScattCherenkovMatrix();

    double MultipleScatteringFraction(const TelescopeDataBin& telDataBin,
                                      const std::vector<double>& waveLengths,
                                      const utl::TabulatedFunction& yield,
                                      const double zeta,
                                      const utl::Point& telescopePosition,
                                      const utl::TabulatedFunction& efficiency) const;

    utl::TabulatedFunctionErrors*
    CalculateLCEff(const fevt::FdConstants::LightSource lightSource,
                   const fevt::TelescopeRecData& telRecData) const;

    void SetTelescopeParameters(const fevt::Eye& eye,
                                TelescopeData& telData) const;

    // telescope data 'random access'
    std::pair<const TelescopeData*, const TelescopeDataBin*> GetTelescopeDataBin(const unsigned int i) const;
    // check for telescope range overlaps
    bool IsOverlapBin(const int i, const int j) const;

    int fVerbosity = 0;
    unsigned int fEyeId = 0;
    unsigned int fNumberOfDepthBins = 0;

    LowerTriangularMatrix fMieScatteredCherenkovMatrix;
    LowerTriangularMatrix fRayScatteredCherenkovMatrix;
    DiagonalMatrix fDirectCherenkovMatrix;
    DiagonalMatrix fDirectFluorescenceMatrix;
    LowerTriangularMatrix fCFMatrix;
    DiagonalMatrix fDirectCFMatrix;

    std::vector<TelescopeData> fTelescopeData;
    // for likelihood calculation of "near border" showers (for PCGF)
    std::vector<TelescopeData> fNoiseTelescopeData;

    // wavelengths to be used for Fluorescence and Cherenkov
    std::vector<double> fCherWaveLength;
    std::vector<double> fFluoWaveLength;

    // transmission vectors
    std::vector<std::vector<double>> fCherTransmissionToTel;
    std::vector<std::vector<double>> fFluoTransmissionToTel;
    std::vector<std::vector<double>> fCherTransmissionShower;
    // scattering vectors
    std::vector<std::vector<double>> fRayScatToTel;
    std::vector<std::vector<double>> fMieScatToTel;
    // multiple scattering factors
    std::vector<double> fFluorescenceMultipleScattering;
    std::vector<double> fCherenkovMultipleScattering;
    // cherenkov at track
    std::vector<std::vector<double>> fCherenkovAtTrack;
    // geometrical factor
    std::vector<double> fGeometricalFactor;
    // distance to shower axis given zeta
    std::vector<double> fZetaDistance;
    // average energy deposit per electron
    std::vector<double> fMeandEdXPerElectron;
    // energy threshold for electrons
    double fElectronEnergyThreshold = 0;
    // shower maximum
    double fXmax = 0;
    // shower axis (physical direction)
    utl::Vector fShowerAxis;
    // cos(zenith) in core CS
    // double fCosTheta;
    // cos(zenith) local in air
    std::vector<double> fCosTheta;

    EMethod fMethod = ePrecise;
    bool fOnlyDirect = false;
    unsigned int fFOVsize = 0;

    // for lateral light distributions
    LateralLightCalculator* fLateralLightCalculator = nullptr;
    bool fDoMultipleScattering = false;

  };

}


#endif