#ifndef __FdRecShower_H
#define __FdRecShower_H

#include <FdRecLevel.h>
#include <RecShower.h>

#include <vector>
#include <map>
#include <iostream>

//
//  Shower data definition
//

class FdRecShower : public RecShower {

public:
  FdRecShower();
  virtual ~FdRecShower() { }

  // ================  getters =================

  // -------- geometrical variables useful for cuts ---------

  /// get angular track length [rad]
  Double_t GetAngTrackObs() const { return fAngTrackObs; }
  /// get time length of track [100 ns]
  Double_t GetTimTrackObs() const { return fTimTrackObs; }
  /// get minimum angle eye-track [rad]
  Double_t GetMinAngle() const { return fMinAngle; }
  /// get maximum angle eye-track [rad]
  Double_t GetMaxAngle() const { return fMaxAngle; }
  /// get mean angle eye-track [rad]
  Double_t GetMeanAngle() const { return fMeanAngle; }
  /// get angular velocity [rad/(100 ns)]
  Double_t GetAngularVelocity() const
  { return fTimTrackObs > 0 ? fAngTrackObs/fTimTrackObs : 6.2831; }
  /// get apparent depth velocity [(g/cm^2)/(100 ns)]
  Double_t GetDepthVelocity() const
  { return fTimTrackObs > 0 ? GetTrackLength()/fTimTrackObs : 1.e4; }

  // -------- shower profile variables ---------

  /// get geometrical field of view minimum boundary [g/cm^2]
  Double_t GetXFOVMin() const { return fXFOVMin; }
  /// get geometrical field of view maximum boundary [g/cm^2]
  Double_t GetXFOVMax() const { return fXFOVMax; }
  /// get (min, max) viewable depth given [g/cm^2] max Xmax uncertainty [g/cm^2] and min MVA [rad], returns min>=max for empty range
  std::pair<double, double> GetExpectedFOVRange(const double xiXmaxCut,
                                                const double mvaCut) const;
  /// get estimated shower maximum uncertainty at Xmax (a.k.a xi(Xmax))
  Double_t GetZXmax() const { return fZXmax; }
  /// get vector of estimated shower maximum uncertainties
  /// (xi(X)) equally spaced between XFOVMin and XFOVMax
  std::vector<double>& GetZVector() { return fZValVec; }
  /// get vector of estimated shower maximum uncertainties
  /// (xi(X)) equally spaced between XFOVMin and XFOVMax
  const std::vector<double>& GetZVector() const { return fZValVec; }
  /// get vector of estimated minimum viewing angles
  /// (xi(X)) equally spaced between XFOVMin and XFOVMax
  std::vector<double>& GetMinViewingAngles() { return fMinViewVec; }
  /// get vector of estimated minimum viewing angles
  /// (xi(X)) equally spaced between XFOVMin and XFOVMax
  const std::vector<double>& GetMinViewingAngles() const { return fMinViewVec; }
  /// get depth for a given z-value (not yet implemented)
  double GetDepthFromZ(const double z) const;
  /// get total depth in FOV [g/cm^2]
  Double_t GetTrackLengthFOV() const { return fXFOVMax - fXFOVMin; }
  /// get minimum slant depth observed [g/cm^2]
  Double_t GetXTrackMin() const { return fXTrackMin; }
  Double_t GetXTrackMinChi() const { return fXTrackMinChi; }
  /// get maximum slant depth observed [g/cm^2]
  Double_t GetXTrackMax() const { return fXTrackMax; }
  /// get total depth observed [g/cm^2]
  Double_t GetTrackLength() const { return fXTrackMax-fXTrackMin; }
  /// get energy deposit at maximum [PeV/(g/cm^2)]
  Double_t GetdEdXmax() const { return fdEdXmax; }

  /// get error of energy deposit at maximum [PeV/gh(g/cm^2)] (deprecated)
  Double_t
  GetdEdXError()
    const
  {
    std::cerr << " FdRecShower::GetdEdXError() is deprecated and will be "
                 " removed soon..." << std::endl;
    return fdEdXmaxError;
  }

  /// get error of energy deposit at maximum [PeV/gh(g/cm^2)]
  Double_t GetdEdXmaxError() const { return fdEdXmaxError; }
  /// get slant depth of first interaction [g/cm^2]
  Double_t GetX1() const { return fX1; }
  /// get error of slant depth of first interaction [g/cm^2]
  Double_t GetX1Error() const { return fX1Error; }
  /// get slant depth of shower maximum [g/cm^2]
  Double_t GetXmax() const { return fXmax; }
  /// get error of slant depth of shower maximum [g/cm^2]
  Double_t GetXmaxError() const { return fXmaxError; }
  /// get chi-angle of first interaction [g/cm^2]
  Double_t GetX1Chi() const { return fX1Chi; }
  /// get error of chi-angle of first interaction [g/cm^2]
  Double_t GetX1ChiError() const { return fX1ChiError; }
  /// get chi-angle of shower maximum [g/cm^2]
  Double_t GetXmaxChi() const { return fXmaxChi; }
  /// get error of chi-angle of shower maximum [g/cm^2]
  Double_t GetXmaxChiError() const { return fXmaxChiError; }
  /// get Gaisser-Hillas shape parameter X0 [g/cm^2]
  Double_t GetX0() const { return fX0; }
  /// get error of Gaisser-Hillas shape parameter X0 [g/cm^2]
  Double_t GetX0Error() const { return fX0Error; }
  /// get Gaisser-Hillas shape parameter lambda [g/cm^2]
  Double_t GetLambda() const { return fLambda; }
  /// get error Gaisser-Hillas shape parameter lambda [g/cm^2]
  Double_t GetLambdaError() const { return fLambdaError; }

  /// check if correlation coefficients are present
  bool HasProfileCorrelations() const;
  /// get error Gaisser-Hillas parameter correlation coefficient
  Double_t GetdEdXmaxXMaxCorrelation() const { return fRhodEdXmaxXMax; }
  /// get error Gaisser-Hillas parameter correlation coefficient
  Double_t GetdEdXmaxLambdaCorrelation() const { return fRhodEdXmaxLambda; }
  /// get error Gaisser-Hillas parameter correlation coefficient
  Double_t GetdEdXmaxX0Correlation() const { return fRhodEdXmaxX0; }
  /// get error Gaisser-Hillas parameter correlation coefficient
  Double_t GetXMaxLambdaCorrelation() const { return fRhoXMaxLambda; }
  /// get error Gaisser-Hillas parameter correlation coefficient
  Double_t GetXMaxX0Correlation() const { return fRhoXMaxX0; }
  /// get error Gaisser-Hillas parameter correlation coefficient
  Double_t GetLambdaX0Correlation() const { return fRhoLambdaX0; }

  /// get chi2 of Gaisser-Hillas fit
  Double_t GetGHChi2() const { return fGHChi2; }
  /// get NdF of Gaisser-Hillas fit
  Int_t GetGHNdf() const { return fGHNdf; }
  /// get chi2 of a linear fit to profile
  Double_t GetLinearProfileFitChi2() const { return fLineFitChi2; }
  /// get first Xmax of double Gaisser-Hillas fit
  Double_t GetDGHXmax1() const { return fDGHXmax1; }
  /// get second Xmax of double Gaisser-Hillas fit
  Double_t GetDGHXmax2() const { return fDGHXmax2; }
  /// get difference in chi-square of single vs double Gaisser-Hillas fit
  Double_t GetDGHChi2Improvement() const { return fDGHChi2Improvement; }
  /// get calorimetric energy [eV]
  Double_t GetEcal() const { return fCalEnergy; }
  /// get error of calorimetric energy [eV]
  Double_t GetEcalError() const { return fCalEnergyError; }
  /// get reconstructed fraction of Cherenkov light [%]
  Double_t GetCherenkovFraction() const { return fChkovFrac; }  // [%]
  /// get height of shower maximum ( [m] a.s.l. )
  Double_t GetHeightXmax() const { return fHeightXmax; }
  /// get distance of shower maximum from eye [m]
  Double_t GetDistXmax() const { return fDistXmax; }
  /// get distance of shower maximum from from core [m]
  Double_t GetCoreXmaxDist() const { return fCoreXmaxDist; }
  /// get light attenuation from shower maximum to eye
  Double_t GetAttXmax() const { return fAttXmax; }
  /// get Mie light attenuation from shower maximum to

  Double_t GetMieAttXmax() const { return fMieAttXmax; }
  /// get Rayleigh light attenuation from shower maximum to eye
  Double_t GetRayleighAttXmax() const { return fAttXmax / fMieAttXmax; }
  /// get weighted integral of Mie light attenuation
  Double_t GetMieCorrection() const { return fMieCorrection; }
  /// get weighted integral of Rayleigh light attenuation
  Double_t GetRayleighCorrection() const { return fRayleighCorrection; }

  /// get uncertainty of total energy from atmospheric uncertainty (Note: this is uncert. is also included in GetEnergyError()!)
  Double_t GetEnergyAtmError() const { return fEnergyAtmError; }

  /// get uncertainty of Xmax from atmospheric uncertainty (Note: this is uncert. is also included in GetXmaxError()!)
  Double_t GetXmaxAtmError() const { return fXmaxAtmError; }

  /// returns true if Xlow+deltaX<Xmax<Xup-deltaX
  bool IsXmaxInFOV(const double deltaX = 0) const;

  /// returns dEdX from Gaisser-Hillas function evaluated at X
  double EvaluateGaisserHillas(const double X) const;

  // ---- the profile itself

  ///returns vector of number of charged particles  at shower track
  const std::vector<Double_t>& GetElectrons() const { return fElectrons; }

  ///returns vector of error on number of charged particles  at shower track
  const std::vector<Double_t>& GetElectronsError() const { return fElectrons_err; }

  ///returns vector of energy deposit at shower track
  const std::vector<Double_t>& GetEnergyDeposit() const { return fEnergyDeposit; }

  ///returns vector of energy deposit at shower track
  const std::vector<Double_t>& GetEnergyDepositError() const { return fEnergyDeposit_err; }

  ///returns vector of slant depth \f$[g/cm^2]\f$
  const std::vector<Double_t>& GetDepth() const { return fDepth; }

  ///returns vector of error on slant depth \f$[g/cm^2]\f$
  const std::vector<Double_t>& GetDepthError() const { return fDepth_err; }

  ///returns vector of heights corresponding to the slant depths from GetDepth \f$[m]\f$
  const std::vector<Double_t>& GetHeight() const { return fHeight; }

  ///returns vector of combined <depth,dEdX,dEdXErr> given slant depth binning \f$[g/cm^2]\f$
  std::vector<std::vector<Double_t> > GetCombineddEdXProfile(const double depthBin) const;

  ///returns vector of approximate shower ages \f$s = 3 / (1 + 2 X_{max}/X)\f$
  std::vector<Double_t> GetShowerAge() const;

  ///was the cloud fraction calculated?
  bool HasMaxCloudFractionBetweenEyeAndShower()
  { return fMaxCloudFracEyeShower >= 0; }

  ///returns the maximum cloud fraction between eye and shower
  double GetMaxCloudFractionBetweenEyeAndShower() const
  { return fMaxCloudFracEyeShower; }

  // ================  setters =================

  // -------- variables useful for cuts ---------

  /// set minimum viewable depth [g/cm^2]
  void SetXFOVMin(const Double_t arg) { fXFOVMin = arg; }
  /// set maximum viewable depth [g/cm^2]
  void SetXFOVMax(const Double_t arg) { fXFOVMax = arg; }
  /// set z = -lg(r_diaphragm^2*AttXmax/DistXmax^2)
  void SetZXmax(const Double_t arg) { fZXmax = arg; }
  /// set observed angular track length [rad]
  void SetAngTrackObs(const Double_t arg) { fAngTrackObs = arg; }
  /// set total depth observed [g/cm^2]
  void SetXTrackObs(const Double_t arg) { fXTrackObs = arg; }
  /// set time length of track [ns]
  void SetTimTrackObs(const Double_t arg) { fTimTrackObs = arg; }
  /// set minimum angle eye-track [rad]
  void SetMinAngle(const Double_t arg) { fMinAngle = arg; }
  /// set maximum angle eye-track [rad]
  void SetMaxAngle(const Double_t arg) { fMaxAngle = arg; }
  /// set mean angle eye-track [rad]
  void SetMeanAngle(const Double_t arg) { fMeanAngle = arg; }

  // -------- shower profile variables ---------

  /// set minimum slant depth observed [g/cm^2]
  void SetXTrackMin(const Double_t arg, const double arg2) { fXTrackMin = arg; fXTrackMinChi = arg2; }
  /// set maximum slant depth observed [g/cm^2]
  void SetXTrackMax(const Double_t arg) { fXTrackMax = arg; }
  /// set energy deposit at maximum [PeV/(g/cm^2)]
  void SetdEdXmax(const Double_t arg) { fdEdXmax = arg; }
  /// set error of energy deposit at maximum [PeV/(g/cm^2)]
  void SetdEdXError(const Double_t arg) { fdEdXmaxError = arg; }
  /// set slant depth of first interaction [g/cm^2]
  void SetX1(const Double_t arg, const Double_t arg2) { fX1 = arg; fX1Chi = arg2; }
  /// set error of slant depth of first interaction [g/cm^2]
  void SetX1Error(const Double_t arg, const double arg2) { fX1Error = arg; fX1ChiError = arg2; }
  /// set slant depth of shower maximum [g/cm^2]
  void SetXmax(const Double_t arg, const Double_t arg2 = 0) { fXmax = arg; fXmaxChi = arg2; }
  /// set error of slant depth of shower maximum [g/cm^2]
  void SetXmaxError(const Double_t arg, const double arg2 = 0) { fXmaxError = arg; fXmaxChiError = arg2; }
  /// set Gaisser-Hillas shape parameter X0 [g/cm^2]
  void SetX0(const Double_t arg) { fX0 = arg; }
  /// set error of Gaisser-Hillas shape parameter X0 [g/cm^2]
  void SetX0Error(const Double_t arg) { fX0Error = arg; }
  /// set Gaisser-Hillas shape parameter lambda [g/cm^2]
  void SetLambda(const Double_t arg) { fLambda = arg; }
  /// set error Gaisser-Hillas shape parameter lambda [g/cm^2]
  void SetLambdaError(const Double_t arg) { fLambdaError = arg; }

  /// set Gaisser-Hillas parameter correlation coefficient
  void SetdEdXmaxXMaxCorrelation(const double rho) { fRhodEdXmaxXMax = rho; }
  /// set Gaisser-Hillas parameter correlation coefficient
  void SetdEdXmaxLambdaCorrelation(const double rho) { fRhodEdXmaxLambda = rho;}
  /// set Gaisser-Hillas parameter correlation coefficient
  void SetdEdXmaxX0Correlation(const double rho) { fRhodEdXmaxX0 = rho; }
  /// set Gaisser-Hillas parameter correlation coefficient
  void SetXMaxLambdaCorrelation(const double rho) { fRhoXMaxLambda = rho; }
  /// set Gaisser-Hillas parameter correlation coefficient
  void SetXMaxX0Correlation(const double rho) { fRhoXMaxX0 = rho; }
  /// set Gaisser-Hillas parameter correlation coefficient
  void SetLambdaX0Correlation(const double rho) { fRhoLambdaX0 = rho; }


  /// set chi2 of Gaisser-Hillas fit
  void SetGHChi2(const Double_t arg) { fGHChi2 = arg; }
  /// set NdF of Gaisser-Hillas fit
  void SetGHNdf(const Int_t arg) { fGHNdf = arg; }
  /// set chi2 of a linear fit to profile
  void SetLinearProfileFitChi2(const Double_t arg) { fLineFitChi2 = arg; }
  /// set first Xmax of double Gaisser-Hillas fit
  void SetDGHXmax1(const Double_t arg) { fDGHXmax1 = arg; }
  /// set second Xmax of double Gaisser-Hillas fit
  void SetDGHXmax2(const Double_t arg) { fDGHXmax2 = arg; }
  /// set difference in chi-square of single vs double Gaisser-Hillas fit
  void SetDGHChi2Improvement(const Double_t arg) { fDGHChi2Improvement = arg; }
  /// set calorimetric energy [eV]
  void SetEcal(const Double_t arg) { fCalEnergy = arg; }
  /// set error of calorimetric energy [eV]
  void SetEcalError(const Double_t arg) { fCalEnergyError = arg; }
  /// set reconstructed fraction of Cherenkov light [%]
  void SetCherenkovFraction(const Double_t arg) { fChkovFrac = arg; }  // [%]
  /// set distance of shower maximum from eye [m]
  void SetDistXmax(const Double_t arg) { fDistXmax = arg; }
  /// set height of shower maximum [m]
  void SetHeightXmax(const Double_t arg) { fHeightXmax = arg; }
  /// set distance of shower maximum from core [m]
  void SetCoreXmaxDist(const Double_t arg) { fCoreXmaxDist = arg; }
  /// set light attenuation from shower maximum to eye
  void SetAttXmax(const Double_t arg) { fAttXmax = arg; }
  /// set Mie light attenuation from shower maximum to eye
  void SetMieAttXmax(const Double_t arg) { fMieAttXmax = arg; }
  /// set weighted integral of Mie light attenuation
  void SetMieCorrection(const Double_t arg) { fMieCorrection = arg; }
  /// set weighted integral of Rayleigh light attenuation
  void SetRayleighCorrection(const Double_t arg) { fRayleighCorrection = arg; }

  /// set uncertainty of total energy from atmospheric uncertainty (stat.)
  void SetEnergyAtmError(const Double_t arg) { fEnergyAtmError = arg; }
  /// set uncertainty of Xmax from atmospheric uncertainty (stat.)
  void SetXmaxAtmError(const Double_t arg) { fXmaxAtmError = arg; }

  // ---- the profile itself

  /// Set slant depth \f$[g/cm^2]\f$
  void SetDepth(const std::vector<Double_t>& dep, const std::vector<Double_t>& dep_err);
  /// Set number of charged particles  at shower track
  void SetElectrons(const std::vector<Double_t>& elec, const std::vector<Double_t>& elec_err);
  /// Set energy deposit at shower track
  void SetEnergyDeposit(const std::vector<Double_t>& edep, const std::vector<Double_t>& edep_err);

  /// sets vector of heights corresponding to the slant depths from GetDepth
  void SetHeight(const std::vector<Double_t>& height) { fHeight = height; }

  /// sets maximum cloud fraction between eye and shower
  void SetMaxCloudFractionBetweenEyeAndShower(const double value)
  { fMaxCloudFracEyeShower = value; }


  void DumpASCII(std::ostream& o = std::cout) const;

  // ================ static functions =================

  /// Calculates approximate shower age from depth and depth of maximum
  static double ShowerAge(const double x, const double xMax);

private:
  // -------- geometrical variables useful for cuts ---------

  Double_t fTimTrackObs;
  Double_t fMinAngle;
  Double_t fMaxAngle;
  Double_t fMeanAngle;

  // -------- shower profile variables ---------

  Double_t fXFOVMin;
  Double_t fXFOVMax;
  Double_t fZXmax;
  Double_t fAngTrackObs;
  Double_t fXTrackObs;
  Double_t fXTrackMinChi;
  Double_t fXTrackMin;
  Double_t fXTrackMax;
  Double_t fdEdXmax;
  Double_t fdEdXmaxError;
  Double_t fX1;
  Double_t fX1Error;
  Double_t fX1Chi;
  Double_t fX1ChiError;
  Double_t fXmax;
  Double_t fXmaxError;
  Double_t fXmaxChi;
  Double_t fXmaxChiError;
  Double_t fX0;
  Double_t fX0Error;
  Double_t fLambda;
  Double_t fLambdaError;
  Double_t fRhodEdXmaxXMax;
  Double_t fRhodEdXmaxLambda;
  Double_t fRhodEdXmaxX0;
  Double_t fRhoXMaxLambda;
  Double_t fRhoXMaxX0;
  Double_t fRhoLambdaX0;

  Double_t fGHChi2;
  Int_t fGHNdf;
  Double_t fLineFitChi2;
  Double_t fDGHXmax1;
  Double_t fDGHXmax2;
  Double_t fDGHChi2Improvement;
  Double_t fCalEnergy;
  Double_t fCalEnergyError;
  Double_t fChkovFrac;
  Double_t fDistXmax;
  Double_t fHeightXmax;
  Double_t fCoreXmaxDist;
  Double_t fAttXmax;
  Double_t fMieAttXmax;
  Double_t fMieCorrection;
  Double_t fRayleighCorrection;
  Double_t fEnergyAtmError;
  Double_t fXmaxAtmError;
  Double_t fMaxCloudFracEyeShower;

  ///slant depth \f$[g/cm^2]\f$
  std::vector<Double_t> fDepth;

  ///error on slant depth \f$[g/cm^2]\f$
  std::vector<Double_t> fDepth_err;

  ///Height \f$[m]\f$
  std::vector<Double_t> fHeight;

  ///number of charged particle at shower track
  std::vector<Double_t> fElectrons;

  ///error on number of particle at shower track
  std::vector<Double_t> fElectrons_err;

  ///number on energy deposit  at shower track
  std::vector<Double_t> fEnergyDeposit;

  ///error on energy deposit at shower track
  std::vector<Double_t> fEnergyDeposit_err;

  ///vector of xi(Xmax) along track
  ///(equally spaced between [fXFOVMin,fXFOVMax]
  std::vector<double> fZValVec;

  ///vector of MVAs along track
  ///(equally spaced between [fXFOVMin,fXFOVMax]
  std::vector<double> fMinViewVec;

  ClassDef(FdRecShower, 16);
};

#endif