// Marjor Changes...   March 27, 2001      T. Abe
// For V3.0 of FastMC, April  2, 1999      J. Bogart
//
#ifndef LCDSWIMTRAJ_H
#define LCDSWIMTRAJ_H

#include "LCDDetectorVolume.h"
#include "TVector3.h"
#include <float.h>   // for DBL_EPSILON


// Provide methods to propagate a particle in a helix or straight line.
// Since it is expected that only one particle will need to be propagated
// at a time, it is possible to re-use a SwimTraj for a new particle.

// Model does *not* include energy loss. 

// Units used are Tesla, GeV/c and cm

class LCDSwimTraj : public TObject {

 public:
  // Return values for PropagateToCyl
  enum PropagateRet {INTERSECT_CYLIN     = -2, 
		     INTERSECT_ZMINUS    = -1,
		     INTERSECT_NONE      =  0,
		     INTERSECT_ZPLUS     =  1, 
		     INTERSECT_CYLOUT    =  2,
                     INTERSECT_XYMOM     =  3,
                     INTERSECT_UNCHARGED =  4,
                     INTERSECT_PATHOK    =  5,
                     INTERSECT_ERROR     =  6 };

  LCDSwimTraj();
  ~LCDSwimTraj() {}

  // Propagate along a suitable helix until it intersects a cylinder
  // with inner radius rCylIn, outer radius rCylOut and ends at 
  // z = +- zCyl.  Return code indicates
  // whether intersection was with -z end, +z end,  inner barrel surface,
  // outer barrel surface,
  // or never happened.  Return path length as well as final position 
  // and momentum vectors.
  // fld argument is in units of Tesla
  PropagateRet HelixToCyl(LCDDetectorVolume *pCur, Double_t *pS);

  // Set up for first call to a propagate routine by storing charge, 
  // init pos, mom, etc.
  void SetParticleParams(const TVector3& pos, const TVector3& mom,
			 Int_t charge);
  void Update();  // copies new pos, mom to slots for initial

  PropagateRet HelixPath(Double_t fld, Double_t s);
  PropagateRet HelixPath(Double_t s);

  // Propagate along a line.  Same call interface as HelixPath.
  PropagateRet LinePath(Double_t s);

  // Similar routine to propagate along straight line to cylinder
  PropagateRet LineToCyl(LCDDetectorVolume *pCur, Double_t *pS);

  TVector3  GetNewPos()    {return  xnew;}
  TVector3* GetNewPosPtr() {return &xnew;}
  TVector3  GetNewMom()    {return  pnew;}
  TVector3* GetNewMomPtr() {return &pnew;}
  TVector3  GetIniPos()    {return  xtrj;}
  TVector3* GetIniPosPtr() {return &xtrj;}
  TVector3  GetIniMom()    {return  ptrj;}
  TVector3* GetIniMomPtr() {return &ptrj;}
  Int_t    GetCharge() {return q;}
  Double_t GetS() {return totalS;}

 private:
  Int_t     q;     // charge

  TVector3  xtrj;   // position vector
  TVector3  ptrj;   // momentum vector
  TVector3  xnew;   // new position vector
  TVector3  pnew;   // new momentum vector

  // Now for a bunch of values used for helix
  Double_t  m_tkpar[5];//helix parameter vector...

  //Double_t  ptor;   // 

  Double_t  deltaS;    // delta path length for last piece of propagation
  Double_t  totalS;    // (geometrical) path length traversed by particle

  void CalcTrackParameters(Double_t fld);

  //  Bool_t checkInside(Double_t rIn, Double_t rOut, Double_t z);
  // This is only called immediately after sending the particle "a
  // little way" but not updating, so we want to use new position vector.
  Bool_t CheckInside(LCDDetectorVolume *pCur) {
    return pCur->isInside(xnew.X(), xnew.Y(), xnew.Z());
  }
  
  Int_t  CircleIntersect(Double_t rv,    // radius of cylinder
			 TVector3& p1,   // (x,y) for 1st solution
			 TVector3& p2);   // (x,y) for 2nd solution
  //check
  Double_t PhaseChange(TVector2 newP);

  // Return largest of roots to quadratic
  static Double_t QuadRoot(const Double_t a, const Double_t b, 
			   const Double_t c, Bool_t* status); 
  PropagateRet ShortestPath(Double_t sEnd, Double_t sIn, 
			    Double_t sOut, Double_t* pS);

 public:
  ClassDef(LCDSwimTraj,2)  // Used to propagate MC particle in calorimeter
};
#endif