// $Header: SmearTraj.h

#ifndef NLD_SMEARTRAJ_H
#define NLD_SMEARTRAJ_H
#include "McPart.h"
#include "SmearVolume.h"

// 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 SmearTraj for a new particle.

// Model does *not* include energy loss. 

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

class SmearTraj : 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};
  

  //  SmearTraj(Double_t fld = 4.0) {fieldCnv = 0.003 * fld;}
  SmearTraj() {}

  // Set up for first call to a propagate routine by storing charge, 
  // init pos, mom, etc.
  void SetParticle(McPart *pPart, Int_t ix);

  void Update();  // copies new pos, mom to slots for initial

  PropagateRet HelixPath(Double_t fld, Double_t s);

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

  // 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(Double_t fld, Double_t rCylIn, Double_t rCylOut,
  //		  Double_t zCyl, Double_t *pS);
  PropagateRet HelixToCyl(SmearVolume *pCur, Double_t *pS);


  // Similar routine to propagate along straight line to cylinder
  //  PropagateRet LineToCyl(Double_t rCylIn, Double_t rCylOut, 
  //			 Double_t zCyl, 
  //			 Double_t *pS);
  PropagateRet LineToCyl(SmearVolume *pCur, Double_t *pS);

  Double_t* GetNewPos() {return &newP[0];}
  Double_t* GetNewMom() {return &newM[0];}
  Int_t     GetCharge() {return q;}
  Double_t  GetS() {return totalS;}
  McPart*   GetMc() {return pPart;}
  Int_t     GetMcIx() {return particleIx;}

private:
  // Must know the field to compute radius of circular projectio of helix.
  // Also will need to multiply it by a constant (approx. 0.003) to deal
  // with converion of units of momentum in GeV/c and distance in cm.

  //  Double_t fieldCnv;
  
  //  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(SmearVolume *pCur) 
    {return pCur->isInside(newP[0], newP[1], newP[2]);}
  
  void   makeHelixQuantities(Double_t fld);  // compute e.g. rad, momMag, etc.

  Int_t  circleIntersect(Double_t rv,    // radius of cylinder
			 Double_t p1[],   // (x,y) for 1st solution
			 Double_t p2[]);   // (x,y) for 2nd solution

  Double_t phaseChange(Double_t newX, Double_t newY);

  // 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); 

  // Numerically stable methods 
  static Double_t pythagPlus(const Double_t x, Double_t y);
  static Double_t pythagMinus(const Double_t x, Double_t y);

  McPart   *pPart;       // ptr to particle to be smeared
  Int_t     particleIx;  // index of part. to be smeared (in list of same)
  Int_t     q;        // charge
  Double_t  iP[3];    // position at starting point
  Double_t  iM[3];    // mom. at starting point
  Double_t  newP[3];  // position after calling a swim routine
  Double_t  newM[3];  // mom. after calling a swim routine

  // Now for a bunch of values used for helix
  Double_t  momMag;    // Magnitude of momentum
  Double_t  xyMomMag;  // Mag. of momentum projected to xy-plane
  Double_t  sinLambda;
  Double_t  cosLambda;

  Double_t  rad;      // radius of helix projected to xy-plane

  Double_t  radQ;     // ..take into account sign of q
  Double_t  xCenter, yCenter;  // ... of helix circle

  Double_t  fieldCnv;  // Parameter derived from mag. field
  Double_t  deltaS;    // delta path length for last piece of propagation
  Double_t  totalS;    // (geometrical) path length traversed by particle

  //  static const Double_t VERY_SMALL_PT = 0.010;
public:
  ClassDef(SmearTraj,1)  // Used to propagate MC particle in calorimeter
};
#endif