#ifndef _evt_GaisserHillas4Parameter_h_
#define _evt_GaisserHillas4Parameter_h_

static const char CvsId_evt_GaisserHillas4Parameter[] =
  "$Id$";

#include <evt/VGaisserHillasParameter.h>
#include <evt/GaisserHillasTypes.h>
#include <utl/TabulatedFunction.h>
#include <iostream>
#include <vector>
#include <utility>

namespace io {
  class GaisserHillasParameters_ROOT;
}

namespace evt {

  /**
     \class GaisserHillas4Parameter
     \brief Gaisser Hillas with 4 parameters

     To get some hints about proper usage of this class see
     the VGaisserHillasParameter documentation.

     \date Mo Jul 25 15:04:16 CEST 2005
     \version $Id$
     \ingroup shower
  */

  class GaisserHillas4Parameter : public VGaisserHillasParameter {

  public:

    typedef std::pair<gh::EShapeParameter, double> ShapeParameter;

    GaisserHillas4Parameter(const gh::EFunctionType functionType = gh::eClassic);
    GaisserHillas4Parameter(const double nmax, const double xmax,
                            const ShapeParameter par1, ShapeParameter par2,
                            const gh::EFunctionType functionType = gh::eClassic);
    virtual ~GaisserHillas4Parameter() { }

    GaisserHillas4Parameter* Create() const
    { return new GaisserHillas4Parameter; }

    GaisserHillas4Parameter* Clone() const
    { return new GaisserHillas4Parameter(*this); }

    gh::EFunctionType GetFunctionType() const { return fFunctionType; }

    /// access to all variants of shape parameters (see GaisserHillasTypes.h)
    double GetShapeParameter(const gh::EShapeParameter par) const;
    double GetShapeParameterError(const gh::EShapeParameter par) const;

    /** currently only the 'correct' internal correlations can be accessed, e.g.
        (eX0, eLambda) for type = eClassic, (eFWHM, eAsym) for type = eWidth etc.
        check with IsInternal(parameter) if in doubt.
    */
    double GetCorrelationXMaxShapeParameter(const gh::EShapeParameter par) const;
    double GetCorrelationNMaxShapeParameter(const gh::EShapeParameter par) const;
    double GetCorrelationShapeParameters() const { return fRhoShapePar1Par2; }

    /// check if parameter "par" is one of the internal shape parameters
    bool IsInternal(const gh::EShapeParameter par) const;

    /// evaluate function a X = depth
    double Eval(const double depth) const;
    /// dump the parameters
    void Dump(std::ostream& os = std::cout) const;
    /// calculate integral
    double GetIntegral() const;
    /// return relative error of integral
    double GetIntegralError() const;

    /// return depth left of XMax for which GH = h*NMax
    double InverseLeft(const double h) const;
    /// return depth right of XMax for which GH = h*NMax
    double InverseRight(const double h) const;

    /// Setters
    /** only the 'correct' internal parameters can be set, e.g.
        (eX0, eLambda) for type = eClassic, (eFWHM, eAsym) for type = eWidth etc.
        check with IsInternal(parameter) if in doubt.
    */
    void SetShapeParameter(const gh::EShapeParameter par,
                           const double value, const double error);
    void SetCorrelationXMaxShapeParameter(const gh::EShapeParameter par,
                                          const double rho);
    void SetCorrelationNMaxShapeParameter(const gh::EShapeParameter par,
                                          const double rho);
    void SetCorrelationShapeParameters(const double rho) { fRhoShapePar1Par2 = rho; }

    // \\deprecated
    double GetXZero() const;
    double GetXZeroError() const;
    double GetLambda() const;
    double GetLambdaError() const;
    double GetNMaxLambdaCorrelation() const;
    double GetNMaxX0Correlation() const;
    double GetXMaxLambdaCorrelation() const;
    double GetXMaxX0Correlation() const;
    double GetLambdaX0Correlation() const;

    friend class io::GaisserHillasParameters_ROOT;

  private:

    double Inverse(const double h, const int branch) const;
    enum EInternalShapeParameter {
      eFirst = 0,
      eLast
    };

    EInternalShapeParameter ExternalToInternal(const gh::EShapeParameter par) const;
    gh::EShapeParameter InternalToExternal(const EInternalShapeParameter) const;
    static gh::EShapeParameter InternalToExternal(const gh::EFunctionType,
                                                  const EInternalShapeParameter);

    double CalculateR(const double asym) const;

    gh::EFunctionType fFunctionType;

    double fShapePar[eLast+1];
    double fShapeParError[eLast+1];
    double fRhoNMaxShapePar[eLast+1];
    double fRhoXMaxShapePar[eLast+1];
    double fRhoShapePar1Par2;

    static utl::TabulatedFunction fRvsAsymTable;

  };

}


#endif