#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     HKalMdcMeasLayer;
class     HKalPlane;


class HKalRungeKutta : public TObject {

private:

    Bool_t     bConstField;      // Use a constant magnetic field instead of fieldMap.
    Bool_t     bDoDEDX;          // Calculate energy loss.
    Bool_t     bDoMS;            // Calculate coulomb multiple scattering.
    Bool_t     bRotBfield;       //! Do/Don't do coordinate transformation of B-field vectors.
    Bool_t     direction;        // Propagation direction: kIterForward/kIterBackward.

    Double_t   fieldScaleFact;   //! Factor of the magnetic field strength.
    TVector3   magField;         //! Use a constant magnetic field instead of fieldMap.
    TRotation *pRotMat;          //! Rotation matrix for possible coordinate transformations  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 Runge-Kutta steps
    Int_t      maxNumStepsPCA;   //! maximum number of attempts to find the point of closest approach of the track to a wire
    Double_t   minLengthCalcQ;   //! minimum step length where process noise is still calculated.

    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     bElossErr;        // Encountered errors during energy loss calculation. Reset in every propagateToPlane() call.
    Bool_t     bPrintElossErr;   // Avoid repeated prints of the same error message.
    Bool_t     bPrintErr;        // Show error messages.
    Bool_t     bPrintWarn;       // Show warning messages.
    Bool_t     trackOk;          // Encountered errors during track propagation.
    Double_t   maxTan;           // Maximum allowed values for tx and ty.
    Double_t   maxPos;           // Maximum allowed values for x and y coordinates.
    Double_t   maxPropMat;       // Maximmum allowed value for an element of the propagator matrix.
    Double_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, Bool_t bCalcJac=kTRUE);

public:

    HKalRungeKutta();

    virtual ~HKalRungeKutta();

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

    virtual Double_t calcDEDXBetheBloch    (Double_t beta, Double_t mass, Double_t ZoverA, Double_t density,
					    Double_t exciteEner, Double_t z=1.);

    virtual Double_t calcDEDXIonLepton     (Double_t qp, Double_t ZoverA, 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 Bool_t   calcMaterialProperties(Double_t &A, Double_t &Z, Double_t &density, Double_t &radLength, Double_t &exciteEner,
                                            const TVector3 &posPreStep, const TVector3 &posPostStep,
                                            const HKalMdcMeasLayer &layerFrom, const HKalMdcMeasLayer &layerTo) const;

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

    virtual void     clear                 ();

    virtual Bool_t   checkPropMat          (TMatrixD &fProp) const;

    virtual Bool_t   checkTrack            (TVectorD &stateVec) const;

    virtual Bool_t   propagateToPlane      (TMatrixD &fProp, TMatrixD &fProc,
                                            TVectorD &stateVecTo, const TVectorD &stateVecFrom,
                                            const HKalPlane &planeFrom, const HKalPlane &planeTo,
                                            const HKalMdcMeasLayer &measLayFrom, const HKalMdcMeasLayer &measLayTo,
                                            Int_t pid, Bool_t propDir, Bool_t bCalcJac=kTRUE);

    virtual void     propagateStraightLine (TVectorD &stateVec, TMatrixD &DF, Double_t &zPos, Double_t dz);

    virtual Bool_t   propagateStraightLine (TVectorD &stateVec, TMatrixD &DF, Double_t &zPos, const HKalPlane &planeTo, Bool_t propDir);

    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     setDoEnerLoss         (Bool_t dedx)             { bDoDEDX = dedx; }

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

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

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

    virtual void     setFieldVector        (const TVector3 &B)       { magField = B; }

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

    virtual void     setPrintErrors        (Bool_t print)            { bPrintErr = print; }

    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     setRungeKuttaParams   (Float_t initialStpSize, Float_t stpSizeDec, Float_t stpSizeInc, Float_t maxStpSize, Float_t minStpSize, Float_t minPrec, Float_t maxPrec, Int_t maxNumStps, Int_t maxNumStpsPCA, Float_t maxDst, Double_t minLngthCalcQ);

    virtual void     setUseConstField      (Bool_t constField)       { bConstField = constField; }

    ClassDef(HKalRungeKutta,0)
};

#endif // HKalRungeKutta_h