#ifndef HMDCKICKPLANE_H
#define HMDCKICKPLANE_H

#include "TObject.h"
#include "hgeomvector.h"

class HMdcKickPlane : public TObject {
protected:
  //Plane equation: y=fDx*x + fDz[i]*z + fC[i]
  Double_t fDzLimit[2];
  Double_t fDx;
  Double_t fDz[3];
  Double_t fC[3];
  
  void setDefaultParam(void);
public:
  HMdcKickPlane(void) {setDefaultParam();}
  ~HMdcKickPlane(void) {}
  inline void calcIntersection(const HGeomVector &r,const HGeomVector &alpha,
			       HGeomVector &out) const;
  void calcIntersection(Double_t x1,Double_t y1,Double_t z1,
                        Double_t x2, Double_t y2, Double_t z2,
                        Double_t& x, Double_t& y, Double_t& z) const;
  
  ClassDef(HMdcKickPlane,0) //Defines a plane surface
};

inline void HMdcKickPlane::calcIntersection(const HGeomVector &p,const HGeomVector &dir,
				            HGeomVector &out) const {
  // Calcul. a cross of the line with plane       y=fDx*x+fDz*z+fC
  for(Int_t i=0;i<3;i++) {
    Double_t c2  = fDz[i]*dir.getZ();
    Double_t c3  = p.getY()-fC[i];
    
    Double_t cDx = fDx;
    Double_t c1  = cDx*dir.getX()-dir.getY();
    Double_t div = 1/(c1+c2);
    Double_t cx  = dir.getX()*(c3-fDz[i]*p.getZ())+ p.getX()*(c2-dir.getY());
    Double_t x   = cx*div;
    if(x<0) {
      cDx = -fDx;
      c1  = cDx*dir.getX()-dir.getY();
      div = 1/(c1+c2);
      x   = cx*div;
    }   
    Double_t z = (dir.getZ()*(c3-cDx*p.getX()) + p.getZ()*c1)*div;
    if(i<2 && z<fDzLimit[i]) continue;
    out.setXYZ(x,cDx*x+fDz[i]*z+fC[i],z);
    break;
  }
}

#endif   /*!HMDCKICKPLANE_H*/