#ifndef __SdRecShower_H
#define __SdRecShower_H

#include <SdRecLevel.h>
#include <RecShower.h>
#include <LDF.h> // for the enums
#include <SdRiseTime.h>
#include <SdFootprintData.h>

#include <TVector3.h>
#include <TH1D.h>

#include <vector>
#include <iostream>

//
//  Shower data definition
//


class SdRecShower : public RecShower {

public:

  SdRecShower();
  virtual ~SdRecShower() { }

  ///> get the chi2 from the plane fit
  double GetAngleChi2() const { return fAngleChi2; }
  ///> get the no of degrees of freedom from the plane fit
  double GetAngleNDoF() const { return fAngleNdof; }

  ///> get curvature, unit: 1/m
  double GetCurvature() const { return fCurvature; }
  ///> get curvature error, unit: 1/m
  double GetCurvatureError() const { return fCurvatureError; }

  ///> get curvature starting point, unit: 1/m
  double GetCurvatureOffset() const {return fCt0;}
  ///> get curvature starting point error, unit: 1/m
  double GetCurvatureOffsetError() const {return fCt0Error;}
  
  ///> get radius of curvature, unit: m
  double GetRadiusOfCurvature() const { return fCurvature ? 1 / fCurvature : 0; }
  ///> get radius of curvature error, unit: m
  double GetRadiusOfCurvatureError() const
  { return fCurvatureError ? fCurvatureError / pow(fCurvature, 2) : 0; }

  ///> get the SD radius point (virtual shower center from curvature fit) in siteCS
  const TVector3& GetRadiusPoint() const { return fRadiusPoint; }
  TVector3& GetRadiusPoint() { return fRadiusPoint; }

  const SdRiseTime& GetRiseTimeResults() const { return fR1000; }
  SdRiseTime& GetRiseTimeResults() { return fR1000; }

  const SdFootprintData& GetFootprint() const { return fFootprint; }
  SdFootprintData& GetFootprint() { return fFootprint; }

  double GetROpt() const { return fRopt; }
  double GetROptError() const { return fRoptError; }

  /// get the mean of the time residual
  double GetTimeResidualMean() const { return fTimeResidualMean; }
  /// get the spread of the time residual
  double GetTimeResidualSpread() const { return fTimeResidualSpread; }

  /// get LDF parameters
  const LDF& GetLDF()const { return fLdf; }
  LDF& GetLDF() { return fLdf; }

  /// get label of shower size estimator
  const char* GetShowerSizeLabel() const { return fLdf.fShowerSizeLabel.c_str();}
  /// get shower size
  double GetShowerSize() const { return fLdf.fS1000; }
  /// get error of shower size
  double GetShowerSizeError() const { return fLdf.fS1000Error; }
  /// get systematic error of shower size
  double GetShowerSizeSys() const { return fLdf.fS1000BetaSys; }

  /// get S1000 unit: VEM (deprecated)
  double GetS1000() const { return fLdf.fS1000; }
  /// get S1000  error unit: VEM (deprecated)
  double GetS1000Error() const { return fLdf.fS1000Error; }
  /// get S1000 sys uncertainty from beta unit: VEM (deprecated)
  double GetS1000BetaSys() const { return fLdf.fS1000BetaSys; }
  /// get total S1000 uncertainty  unit: VEM (deprecated)
  double GetS1000TotalError() const
  { return sqrt(pow(fLdf.fS1000BetaSys, 2) + pow(fLdf.fS1000Error, 2)); }
  /// get beta
  double GetBeta() const { return fLdf.fBeta; }
  /// get beta error
  double GetBetaError() const { return fLdf.fBetaError; }
  /// get gamma
  double GetGamma() const { return fLdf.fGamma; }
  /// get gamma  error
  double GetGammaError() const { return fLdf.fGammaError; }

  /// get the chi2 from the LDF fit
  double GetLDFChi2() const { return fLdf.fLDFChi2; }
  /// get the NDoF from the LDF fit
  double GetLDFNdof() const { return fLdf.fLDFNdof; }
  /// get the maximum likelihood fit estimate
  double GetLDFLikelihood() const { return fLdf.fLDFLikelihood; }
  /// get the reconstruction status
  double GetLDFStatus() const { return fLdf.fLDFStatus; }

  double SignalUncertainty() const { return  0.32 + 0.42 / GetAxisCoreCS().Z(); }

  /// get muon production depth histogram (see GAP-2006-029 and 2010-100)
  bool HasMPDHistogram() const { return !fMPDData.empty(); }
  TH1D GetMPDHistogram() const;
  const std::vector<Double_t>& GetMPDData() const { return fMPDData; }
  double GetMPDMin() const { return fMPDMin; }
  double GetMPDMax() const { return fMPDMax; }

  // ----------- setters --------------------------------------

  void SetAngleChi2(const double chi2, const double ndof)
  { fAngleChi2 = chi2; fAngleNdof = ndof; }

  void SetTimeResidualMean(const double mean) { fTimeResidualMean = mean; }
  void SetTimeResidualSpread(const double spread) { fTimeResidualSpread = spread; }

  void SetCurvatureOffset(const double ct0, const double ct0Error){
    fCt0 = ct0; fCt0Error = ct0Error; 
  }
  
  void SetCurvature(const double curvature, const double error)
  { fCurvature = curvature; fCurvatureError = error; }

  /// set the SD radius point (shower front curvature fit) in siteCS
  void SetRadiusPoint(const TVector3& p) { fRadiusPoint = p; }

  void SetRiseTime1000(const double rt1000, const double rt1000er)
  { fRiseTime1000 = rt1000; fRiseTime1000Error = rt1000er; }

  void SetROpt(const double ropt, const double error)
  { fRopt = ropt; fRoptError = error; }

  void SetMPDMin(const double c) { fMPDMin = c; }
  void SetMPDMax(const double c) { fMPDMax = c; }
  void SetMPDData(const std::vector<double>& mpd) { fMPDData = mpd; }

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

private:
  double fCurvature;
  double fCurvatureError;
  double fCt0;
  double fCt0Error;

  TVector3 fRadiusPoint; //< the SD radius point (shower front curvature fit) in siteCS

  double fRiseTime1000;
  double fRiseTime1000Error;

  double fAngleChi2;
  double fAngleNdof;

  double fTimeResidualMean;
  double fTimeResidualSpread;

  double fRopt;
  double fRoptError;

  LDF fLdf;

  SdRiseTime fR1000;
  SdFootprintData fFootprint;

  std::vector<Double_t> fMPDData;
  double fMPDMin;
  double fMPDMax;

  ClassDef(SdRecShower, 18);
};


#endif