#ifndef HKalRungeKutta_h
#define HKalRungeKutta_h

// from ROOT
#include "TMatrixD.h"
#include "TObjArray.h"
#include "TRotation.h"
class     TVector3;
#include "TVectorD.h"

// from hydra
#include "hmdctrackgfield.h"

#include "hkaldef.h"
class     HKalPlane;
class     HKalMdcMeasLayer;

using namespace Kalman;


class HKalRungeKutta : public TObject {

private:

    Bool_t     bDoDEDX;          //! Calculate energy loss.
    Bool_t     bDoMS;            //! Calculate coulomb multiple scattering.
    Bool_t     direction;        //! Propagation direction: kIterForward/kIterBackward.
    Double_t   fieldScaleFact;   //! Factor of the magnetic field strength.
    TRotation *pRotMat;          //! Rotation matrix for possible coordinate transformations  of B-field vectors.
    Bool_t     bRotBfield;       //! Do/Don't do coordinate transformation of B-field vectors.
    Double_t   trackPosAtZ;      //! Store z position of track during Runge-Kutta propagation.

    Double_t   energyLoss;       //! Energy loss during stepping in MeV/particle charge.
    Double_t   trackLength;      //! Track length in mm.
    Double_t   stepLength;       //! Length of last RK step in mm.
    Int_t      jstep;            //! Step number.
    Float_t    initialStepSize;  //! initial RK step size
    Float_t    stepSizeInc;      //! factor to increase stepsize
    Float_t    stepSizeDec;      //! factor to decrease stepsize
    Float_t    maxStepSize;      //! maximum stepsize
    Float_t    minStepSize;      //! minimum stepsize
    Float_t    maxPrecision;     //! maximum required precision for a single step
    Float_t    minPrecision;     //! minimum required precision for a single step
    Float_t    maxDist;          //! maximum distance for straight line approximation
    Int_t      maxNumSteps;      //! maximum number of steps
    Bool_t     bFillPointArrays; // == kTRUE if point arrays should be filled
    TObjArray  pointsTrack;      // Points along trajectory stored
    TObjArray  fieldTrack;       // Field at Points along trajectory stored

    Bool_t     bPrintWarn;       //! Show warning messages.
    Bool_t     tx2ty2Warn;       //! Warning about large tx and/or ty already shown for track.
    Double_t   maxTan;           //! Maximum allowed values for tx and ty.
    Double_t   maxPos;           //! Maximum allowed values for x and y coordinates.
    Float_t    minMom;           //! Minimum required value for momentum in MeV/c.

    HMdcTrackGField *fieldMap;   //! Pointer to B-field map.

protected:

    virtual Double_t rkSingleStep(TVectorD &stateVec, TMatrixD &fPropStep, Double_t stepSize);

public:

    HKalRungeKutta();

    virtual ~HKalRungeKutta();

    virtual Double_t calcEnergyLoss       (TMatrixD &fProc, TVectorD &stateVec, Double_t beta, Double_t length, Double_t mass, Double_t qp,
            Double_t A, Double_t Z, Double_t density, Double_t radLength, Double_t exciteEner, Int_t pid) const;

    virtual Double_t calcDEDXBetheBloch    (Double_t beta, Double_t mass, Double_t A, Double_t Z, Double_t density, Double_t exciteEner, Double_t z=1.) const;

    virtual Double_t calcDEDXIonLepton     (Double_t qp, Double_t A, Double_t Z, Double_t density, Double_t exciteEner, Int_t pid) const;

    virtual Double_t calcRadLoss           (TMatrixD &fProc, Double_t length, Double_t mass, Double_t qp, Double_t radLength) const;

    virtual TVector3 calcFieldIntegral     () const;

    virtual void     calcField_mm          (const TVector3& xv, TVector3& btos, Double_t fpol) const;

    virtual void     calcField_mm          (Double_t* xv, Double_t* btos,Double_t fpol) const;

    virtual void     calcMaterialProperties(Double_t &A, Double_t &Z, Double_t &density, Double_t &radLength, Double_t &exciteEner,
            const TVector3 &posPreStep, const TVector3 &posPostStep,
            const HKalMdcMeasLayer &planeFrom, const HKalMdcMeasLayer &planeTo) const;

    virtual void     calcProcNoise         (TMatrixD &fProc, const TVectorD &stateVec, Double_t length, Double_t radLength, Double_t beta) const;

    virtual void     clear                 ();

    virtual Bool_t   checkTrack            (TVectorD &statevec) const;

    virtual Double_t distance              (const TVector3 &vec1, const TVector3 &vec2) const;

    virtual Bool_t   propagateToPlane      (TMatrixD &fProp, TMatrixD &fProc, TVectorD &stateVecTo, const TVectorD &stateVecFrom, const HKalMdcMeasLayer &planeFrom, const HKalMdcMeasLayer &planeTo, Int_t pid, Bool_t propDir);

    virtual void     propagateStraightLine (TVectorD &stateVec, TMatrixD &DF, const HKalPlane &planeTo);

    virtual Double_t         getEnergyLoss      () const                  { return energyLoss; }

    virtual Int_t            getNrStep          () const                  { return jstep; }

    virtual Double_t         getTrackLength     () const                  { return trackLength; }

    virtual Double_t         getStepLength      () const                  { return stepLength; }

    virtual TObjArray const& getPointsTrack     () const                  { return pointsTrack; }

    virtual TObjArray const& getFieldTrack      () const                  { return fieldTrack; }

    virtual void     setDirection       (Bool_t dir)              { direction = dir; }

    virtual void     setFieldFact       (Double_t fpol)           { fieldScaleFact = fpol; }

    virtual void     setFieldMap        (HMdcTrackGField *fpField, Double_t fpol) { fieldMap = fpField; fieldScaleFact = fpol; }

    virtual void     setFillPointsArrays(Bool_t fill)             { bFillPointArrays = fill; }

    virtual void     setPrintWarnings   (Bool_t print)            { bPrintWarn = print; }

    virtual void     setRotationMatrix  (TRotation *pRM)          { pRotMat = pRM; }

    virtual void     setRotateBfieldVecs(Bool_t rotB)             { bRotBfield = rotB; }

    virtual void     switchMS           (Bool_t ms)               { bDoMS = ms; }

    virtual void     switchDEDX         (Bool_t dedx)             { bDoDEDX = dedx; }

    ClassDef(HKalRungeKutta,0)
};

#endif // HKalRungeKutta_h