//*-- Authors: A.Ivashkin && A.Sadovsky (04.11.2004)
#ifndef HRUNGEKUTTA_H
#define HRUNGEKUTTA_H

#include "TString.h"
#include "TObject.h"
#include "TNamed.h"
#include "TMatrixD.h"
#include "hgeomvector.h"
#include "hgeomrotation.h"
#include "hgeomtransform.h"
#include "hmdctrackgfield.h"

class HMdcSizesCells;

class HRungeKutta {
 private:
   HMdcTrackGField* fieldMap;              // pointer to field map
   Double_t fieldScalFact;                 // factor of the magnetic field strength
   Bool_t mdcInstalled[4][6];              // remembers which MDCs have known geometry
   Double_t mdcPos[4][6][3];               // real position of mdc's in sector coordinate system
   Double_t normVecMdc[4][6][3];           // normal vector on each Mdc module in sector coordinate system
   Double_t pointForVertexRec[6][3];       // point on virtual plane in front of MDC1 used for vertex reconstruction
   Double_t dzfacMdc[2][6];                // factor to calculate dz from dy in the first two MDC planes
   Double_t p0Guess;                       // initial momentum

   //--------------------------------- RK-tracking ----------------------------
   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  minPrecision;                  // minimum required precision for a single step
   Float_t  maxPrecision;                  // maximum required precision for a single step
   Int_t    maxNumSteps;                   // maximum number of steps
   Double_t maxDist;                       // maximum distance for straight line approximation
   Float_t  pfit;                          // momentum after RK-fit
   Float_t  chi2;                          // normalized chi^2 after RK-fit
   Float_t  trackLength;                   // track length
   Float_t  dpar[5];                       // variations (momentum, xdir, ydir, xpos, ypos)
   Float_t  resolution[8];                 // resolution
   Float_t  sig[8];                        // sigmas for hit points
   Double_t hit[4][3],fit[4][3];           // hits and fitted points at mdc's
   Double_t fitdp[4][3],fitdh[4][3],fitdv[4][3],fit1h[4][3],fit1v[4][3]; // fitted positions for derivatives
   Int_t    mdcModules[4];                 // module number in hit/fit
   Int_t    mdcLookup[4];                  // index in mdcModules for all 4 modules
   Double_t dydpar[5][8];                  // derivatives
   Double_t ydata[8];                      // difference between fitted and measured hit positions
   Double_t fitpar[5];                     // fit parameters
   Int_t    ndf;                           // number of degrees of freedom
   Int_t    nMaxMod;                       // maximum hit module number (3 or 4)
   Int_t    sector;                        // sector index
   Double_t dirAtFirstMDC[3];              // fitted direction at first MDC
   Double_t posNearLastMDC[3];             // fitted direction at last MDC
   Double_t dirAtLastMDC[3];               // fitted direction at last MDC
   Float_t  stepSizeAtLastMDC;             // RK step size at last MDC
   Float_t  trackLengthAtLastMDC;          // track length from target to last MDC
   Double_t trackPosOnMETA[3];             // fitted trajectory point on META
   Int_t    jstep;                         // step number
   Double_t zSeg1,rSeg1,thetaSeg1,phiSeg1; // fitted data for inner segment
   Double_t zSeg2,rSeg2,thetaSeg2,phiSeg2; // fitted data for outer segment
   Float_t  polbending;                    // polarity of bending (field scaling factor * pfit)
   Double_t mTol;                          //! tolerance for invertion of matrix (default : DBL_EPSILON)
public:
   HRungeKutta();
   virtual ~HRungeKutta() {}
   void clear();
   void setField(HMdcTrackGField*,Float_t);
   void setFieldFactor(Float_t fpol) {fieldScalFact = fpol;}

   void setMdcPosition(Int_t,Int_t,HGeomTransform&);
   Bool_t fit4Hits(Double_t*,Double_t*,Double_t*,Double_t*,Float_t*,Int_t,Double_t pGuess=999999999.);
   Bool_t fit3Hits(Double_t*,Double_t*,Double_t*,Double_t*,Float_t*,Int_t,Double_t pGuess=999999999.);
   void traceToVertex(HMdcSizesCells*);
   void traceToMETA(HGeomVector&,HGeomVector&,HGeomVector* point = NULL,HGeomVector* norm = NULL);

   Float_t getPfit() {return pfit;}
   Float_t getChi2() {return chi2/Float_t(ndf);}
   Int_t getNMaxMod() {return nMaxMod;}

   // fitted points on MDC in Segment Coord System, Mdc={1,2,3,4}
   Float_t getXfit(Int_t iMdc) {
     return (mdcLookup[iMdc-1]!=-1) ? fit[(mdcLookup[iMdc-1])][0] : -999.;
   }
   Float_t getYfit(Int_t iMdc) {
     return (mdcLookup[iMdc-1]!=-1) ? fit[(mdcLookup[iMdc-1])][1] : -999.;
   }
   Float_t getZfit(Int_t iMdc) {
     return (mdcLookup[iMdc-1]!=-1) ? fit[(mdcLookup[iMdc-1])][2] : -999.;
   }

   // hit points on MDC (from segments) in Segment Coord System, Mdc={1,2,3,4}
   Float_t getXhit(Int_t iMdc) {
     return (mdcLookup[iMdc-1]!=-1) ? hit[(mdcLookup[iMdc-1])][0] : -999.;
   }
   Float_t getYhit(Int_t iMdc) {
     return (mdcLookup[iMdc-1]!=-1) ? hit[(mdcLookup[iMdc-1])][1] : -999.;
   }

   // difference between of fitted points and hit points
   Float_t getXfithit(Int_t iMdc) {
     return (mdcLookup[iMdc-1]!=-1)
       ? (fit[(mdcLookup[iMdc-1])][0] - hit[(mdcLookup[iMdc-1])][0]) : -999.;
   }
   Float_t getYfithit(Int_t iMdc) {
     return (mdcLookup[iMdc-1]!=-1)
       ? (fit[(mdcLookup[iMdc-1])][1] - hit[(mdcLookup[iMdc-1])][1]) : -999.;
   }

   // position and direction on first MDC after RungeKutta fit
   Float_t getXtrackFirstMDCFitPos()  { return fit[0][0]; }        // x
   Float_t getYtrackFirstMDCFitPos()  { return fit[0][1]; }        // y
   Float_t getZtrackFirstMDCFitPos()  { return fit[0][2]; }        // z
   Float_t getDXtrackFirstMDCFitPos() { return dirAtFirstMDC[0]; } // dx
   Float_t getDYtrackFirstMDCFitPos() { return dirAtFirstMDC[1]; } // dy
   Float_t getDZtrackFirstMDCFitPos() { return dirAtFirstMDC[2]; } // dz

   // position and direction on last MDC after RungeKutta fit
   Float_t getXtrackLastMDCFitPos()  { return fit[nMaxMod-1][0]; } // x
   Float_t getYtrackLastMDCFitPos()  { return fit[nMaxMod-1][1]; } // y
   Float_t getZtrackLastMDCFitPos()  { return fit[nMaxMod-1][2]; } // z
   Float_t getDXtrackLastMDCFitPos() { return dirAtLastMDC[0]; }   // dx
   Float_t getDYtrackLastMDCFitPos() { return dirAtLastMDC[1]; }   // dy
   Float_t getDZtrackLastMDCFitPos() { return dirAtLastMDC[2]; }   // dz

   // position and direction on META after RungeKutta fit
   Float_t getXtrackOnMETA()  { return trackPosOnMETA[0];} // x
   Float_t getYtrackOnMETA()  { return trackPosOnMETA[1];} // y
   Float_t getZtrackOnMETA()  { return trackPosOnMETA[2];} // z

   Float_t getTrackLength() { return trackLength; } // RK-track length

   Float_t getZSeg1() { return zSeg1; }         // z of inner segment after RK-fit
   Float_t getRSeg1() { return rSeg1; }         // r of inner segment after RK-fit
   Float_t getThetaSeg1() { return thetaSeg1; } // theta of inner segment after RK-fit
   Float_t getPhiSeg1() { return phiSeg1; }     // phi of inner segment after RK-fit

   Float_t getZSeg2() { return zSeg2; }         // z of outer segment after RK-fit
   Float_t getRSeg2() { return rSeg2; }         // r of outer segment after RK-fit
   Float_t getThetaSeg2() { return thetaSeg2; } // theta of outer segment after RK-fit
   Float_t getPhiSeg2() { return phiSeg2; }     // phi of outer segment after RK-fit

   void     setMTol(Double_t tol) { mTol=tol; }
   Double_t getMTol()             { return mTol; }

 private:
   Bool_t   fitMdc();
   void     findMdcIntersectionPoint(Double_t*,Double_t*,Int_t,Double_t*);
   Bool_t   findIntersectionPoint(Double_t*,Double_t*,Double_t*,Double_t*,Double_t*);
   void     initMom();
   Float_t  distance(Double_t*,Double_t*);
   void     parSetNew0(Float_t,Float_t,Int_t);
   void     parSetNewVar(Float_t,Float_t);
   void     cosines(Float_t,Float_t,Double_t*);
   void     gentrkNew(Float_t,Double_t*,Double_t*,Double_t*,Int_t);
   Float_t  rkgtrk(Float_t,Float_t,Double_t*,Double_t*,Int_t kind=0);
   void     rkeqfw(Double_t*,Float_t,Double_t*,Float_t*);
   void     derive(Double_t*,Int_t);
   Bool_t   linSys();
   Bool_t   decompose(TMatrixD &lu,Double_t &sign,Double_t tol,Int_t &nrZeros);
   ClassDef(HRungeKutta,0)
};


#endif