#ifndef _TrackingVariables_HH_
#define _TrackingVariables_HH_ 1
////////////////////////////////////////////////////////////////////
//   File:  trackingVariables.hh
//   Authors: Heather Gerberich
//   Purpose: High level analysis tracking variables defined as functions
//   Created: 10/30/2000
//
//   Revision:  0.2
//
////////////////////////////////////////////////////////////////////

#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "Trajectory/Angle.hh"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "TrackingSI/Utils/SiExpected.hh"

class HepVector3D;
class HepPoint3D;
class Hep3Vector;
class HepSymMatrix;
class CdfTrack;
class CdfTrackView;
class CT_Residual;
class Vertex;
class ZVertex;
class TofMatches;
class TofMatch;
class Beamline;

namespace TrackingVariables {

//Returns pointer to the best fit CdfTrack which includes most information
//(is lowest in the parent-child track link list) 
//for now assumes incoming track is valid
const CdfTrack * trkBestFit(const CdfTrack * track);

//Returns pointer to CdfTrack which contains only COT information
//if no Cot Only Track exists return error message???
//for now assumes incoming track is valid
const CdfTrack * trkCotOnly(const CdfTrack * track);

//in following functions, if trkFlag = 0, use best fit
//in following functions, if trkFlag = 1, use COT
//otherwise use input track if trkFlag does not equal 0 or 1


//Return pointer to desired track ("COT Only", "Best Fit", or input track) 
const CdfTrack * trkGetDesiredTrack(const CdfTrack * track, int trkFlag=2);

// Momentum 3-Vector
Hep3Vector trkMomentum(const CdfTrack * track, int trkFlag=2);

// Px
float trkPx(const CdfTrack * track, int trkFlag=2);

// Momentum 3-Vector
float trkPy(const CdfTrack * track, int trkFlag=2);

// Momentum 3-Vector
float trkPz(const CdfTrack * track, int trkFlag=2);

// Track Pt
float trkPt(const CdfTrack * track, int trkFlag=2);

// Track P
float trkP(const CdfTrack * track, int trkFlag=2);

// Track Curvature
float trkCurvature(const CdfTrack * track, int trkFlag=2);

// track Phi at point of closest approach to the Z-axis
Angle trkPhi0(const CdfTrack * track, int trkFlag=2);

// Track Lamda = Cot(Theta) where Theta is the Polar Angle of the track
float trkLambda(const CdfTrack * track, int trkFlag=2);

// track Eta
float trkEta0(const CdfTrack * track, int trkFlag=2);

// track Charge
float trkCharge(const CdfTrack * track, int trkFlag=2);

// track d0 (impact parameter)
float trkD0(const CdfTrack * track, int trkFlag=2);

// track X0
float trkX0(const CdfTrack * track, int trkFlag=2);

// track Y0
float trkY0(const CdfTrack * track, int trkFlag=2);

// track Z0
float trkZ0(const CdfTrack * track, int trkFlag=2);

// corrected D0 to a Beamline (comes from CotBeam or SvxBeam)
double trkCorrectedD0(const CdfTrack* track, Beamline& beamline);

// corrected D0 to an arbitrary point
double trkCorrectedD0(const CdfTrack* track, double vx, double vy);

// Chi-Squared
float trkChi2(const CdfTrack * track, int trkFlag=2);

  // Track ID # in event
  int trkId(const CdfTrack * track, int trkFlag=2);

// Number of Degrees of Freedom
int trkNumDOF(const CdfTrack * track, int trkFlag=2);

//Number of hits in specified superlayer ***NOTE change in order of parameters****
int trkNumHitsInSL(int superLayer, const CdfTrack * track, int trkFlag=2);

// Number of COT Axial Superlayers hit by track
int trkNumCOT_ASL(const CdfTrack * track, int trkFlag=2, int minHits=1);

// Number of COT Stereo Superlayers hit by track
int trkNumCOT_SSL(const CdfTrack * track, int trkFlag=2, int minHits=1);

// Number of Silicon Hits on Track
int trkNumSVX_Hits(const CdfTrack * track, int trkFlag=2);

  // bitmap for number of hits in each superlayer: nibble j contains the
  // number of hits for superlayer j
  int trkSLBitmap(const CdfTrack * track, int trkFlag=2);

//COVARIANCE MATRIX

//Covariance Matrix for Pt, d0, phi0, z0, Cot(Theta) (symmetric Matrix)
//Cov Matrix = (Sigma Pt^2     Cov(Pt, D0)    Cov(Pt, Phi0)   Cov(Pt, Z0)    Cov(Pt, Lambda)   )
//             (               Sigma D0^2     Cov(D0, Phi0)   Cov(D0, Z0)    Cov(D0, Lambda)   )
//             (                              Sigma Phi0^2    Cov(Phi0,Z0)   Cov(Phi0, Lambda) )
//             (                                              Sigma Z0^2     Cov(Z0, Lambda)   )
//             (                                                             Sigma Lambda^2    )
//Sigma Pt = k*SigmaCurvature/Curvature^2  where k is Pt*Curvature, equal to +/- 0.0021159
//Cov(Pt,x) = -k*cov(Curvature,x)/Curvature^2
//Cov(x,y) = Correlation(x, y)*SigmaX*SigmaY
HepSymMatrix trkCovMtrx(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Pt, Pt)
float trkCovMtrxPtPt(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Pt, D0)
float trkCovMtrxPtD0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Pt, Phi0)
float trkCovMtrxPtPhi0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Pt, z0)
float trkCovMtrxPtZ0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Pt, Lambda)
float trkCovMtrxPtLambda(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (D0, D0)
float trkCovMtrxD0D0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (D0, Phi0)
float trkCovMtrxD0Phi0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (D0, z0)
float trkCovMtrxD0Z0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (D0, Lambda)
float trkCovMtrxD0Lambda(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Phi0, Phi0)
float trkCovMtrxPhi0Phi0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Phi0, z0)
float trkCovMtrxPhi0Z0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Phi0, Lambda)
float trkCovMtrxPhi0Lambda(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (z0, z0)
float trkCovMtrxZ0Z0(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (z0, Lambda)
float trkCovMtrxZ0Lambda(const CdfTrack * track, int trkFlag=2);

//Covariance Matrix element (Lambda, Lambda)
float trkCovMtrxLambdaLambda(const CdfTrack * track, int trkFlag=2);


//EXTRAPOLATOR 

//Used in Extrapolator functions to determine if track position/
//slope is valid at a given r
bool trkRIsValid(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated Position of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for x,y,z
HepPoint3D trkExtPos(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated Slope of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for xslope, yslope, zslope
HepVector3D trkExtSlope(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated Pathlength of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for pathlength
float trkExtPathLength(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated x of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for x
float trkExtX(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated y of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for y
float trkExtY(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated Z of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for z
float trkExtZ(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated X Slope of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for xslope
float trkExtXSlope(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated Y Slope of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for yslope
float trkExtYSlope(float radius, const CdfTrack * track, int trkFlag=2);

//Extrapolated Z Slope of Track at some Radius (less than radius of solenoid for now)
//calculated using information from Helix object returned by CdfTrack::getTrajectory()
//if track will not reach requested radius or impact paramater is greater than radius,
//return -999999. for zslope
float trkExtZSlope(float radius, const CdfTrack * track, int trkFlag=2);
 

//RESIDUALS

//returns a pointer to a CT_Residual object for the specified superLayer 
const CT_Residual * trkResidual(int superLayer, const CdfTrack*track, int trkFlag=2);

//Delta D for a specified superlayer and wire.  deltaD is hitD - |trackY|
float trkDeltaDRes(int superLayer, int wire, const CdfTrack*track, int trkFlag=2);

//Delta Y for a specified superlayer and wire.  
//DeltaY is hitD * (trackY / |trackY|) - trackY
float trkDeltaYRes(int superLayer, int wire, const CdfTrack*track, int trkFlag=2);

//Hit D residual for a specified superlayer and wire
//hitD is hitT * driftspeed + a drift model correction that adjusts for
//     nonlinearities in the time to distance relationship
float trkHitDRes(int superLayer, int wire, const CdfTrack*track, int trkFlag=2);

//Track Y residual for a specified superlayer and wire
//trackY is the distance from the wire to the fitted track along the drift
//       direction. It's positive (negative) if the track is
//       counterclockwise (clockwise) from the wire.
float trkTrackYRes(int superLayer, int wire, const CdfTrack*track, int trkFlag=2);

//Sin Aspect residual for a specified superlayer and wire
//sinAspect is the sin of the angle between the particle momentum and the
//          drift direction. It's positive for positive charged
//          tracks coming from near the origin.
float trkSinAspectRes(int superLayer, int wire, const CdfTrack*track, int trkFlag=2);

//Hit T residual for a specified superlayer and wire
//hitT is the tdc time minus the interaction time, time of flight, signal
//     propagation time, and channel to channel calibration correction,
//     essentially the drift time.
float trkHitTRes(int superLayer, int wire, const CdfTrack*track, int trkFlag=2);

//Hit width of a specified superlayer and wire with a valid residual
int trkHitWidth(int superLayer, int wire, const CdfTrack*track, int trkFlag=2);

//Cell number a specified superlayer and wire with a valid residual
int trkHitCellNum(int superLayer, int wire, const CdfTrack*track, int trkFlag=2); 

//returns the inverse of beta (c/v) using residuals for calculation
//returns the error on the inverse of beta (c/v) (Kotwal)
void trkMinChi2BetaInv(float & MinChi2BetaInv, float & DeltaBetaInv,const CdfTrack * track, int trkFlag=2);

//returns truncated hit width mean using residuals for calculation  (Kotwal)
float trkHitWidthMean(const CdfTrack * track, int trkFlag=2);

//returns truncated hit width mean using residuals for calculation, 
//path length corrected (Kotwal)
float trkHitWidthMeanPLC(const CdfTrack * track, int trkFlag=2);

// returns the track t0 using residuals for calculation
// units are cm of drift time, should be divided by COT drift velocity
// to get t0 in ns
// particles are assumed to move at speed of light
// also returns error on t0, chi2 at minimum and number of hits used
void trkMinChi2t0( float & MinChi2t0, float & Deltat0, float & Chi2t0,                           int & NWirest0, const CdfTrack * track, int trkFlag=2);

// uses track residuals for a 2D fit to 1/Beta deviation from light speed (c/v-1)
// and track production time (t0) 
// returns the two best fit values, errors and covariance
// also returns chi2 at minimum and number of residuals used
void trkMinChi22DBetaInvt0( float & MinChi22Dt0, float & Delta2Dt0,           
                        float & MinChi22DBetaInv, float & Delta2DBetaInv,    
                        float & Cov2D2DBetaInvt0, float & Chi22DBetaInvt0,
                        int & NWiresBetaInvt0,                                
                        const CdfTrack * track, int trkFlag=2);

  // count hit patterns for rejecting decays in flight
  // max_seq = the longest sequence of hits on one side of the track;
  // n_trans = the number of times the hits switch from one side
  //   of the track to the other
  // return value is (n_trans<<8) | max_seq
  int trkTransitionPattern(const CdfTrack& track);

//************  Tracking Isolation

  // The following are overloaded ways to call the tracking isolation
  // currently flags are not implemented
  // Tracking isolation, track-based arguments
  double trkTrackIsolation(double phi, double eta, double z,
                           double cone = 0.4, double dzMaxCut = 5.0, 
                           double ptMinCut = 0.0, int flags = 0);

  double trkTrackIsolation(double phi, double eta, double z,
			   const CdfTrackView& trview,
                           double cone = 0.4, double dzMaxCut = 5.0, 
                           double ptMinCut = 0.0, int flags = 0);

  // Tracking isolation, track-based arguments
  double trkTrackIsolation(const CdfTrack& seedTrack,
                           double cone = 0.4, double dzMaxCut = 5.0, 
                           double ptMinCut = 0.0, int flags = 0);

  // Tracking isolation, track-based arguments, specified track view
  double trkTrackIsolation(const CdfTrack& seedTrack,
			   const CdfTrackView& trview,
                           double cone = 0.4, double dzMaxCut = 5.0, 
                           double ptMinCut = 0.0, int flags = 0);

  //Old utility routine with almost all the arguments
  double trkTrackIsolationCore(double phi, double eta, double z,
                const CdfTrackView& trackList, const CdfTrack* seedTrack,
                double cone, double dzMaxCut, double ptMinCut, int flags);

  //Utility routine with all the arguments
  double trkTrackIsolationInnerCore(double phi, double eta, double z,
                const CdfTrackView& trackList, const CdfTrack* seedTrack,
                double cone, double dzMaxCut, double ptMinCut, int flags,
				     int& nTrkInCone, double& ptMax);


//************  Vertices
// from vxprim.cc -> VertexColl

  //Primary Vertex
  const Vertex *  trkPrimaryVertex();

  //Primary Vertex
  int trkNPrimaryVertex();

  //Z of Primary Vertex
  double trkPrimaryVertexZ();

  //Ht of Primary Vertex
  double trkPrimaryVertexHt();

  //Number of tracks associated with the primary Vertex
  int trkPrimaryVertexNTracks();


//************  ZVertices
// from ZVertexModule -> ZVertexColl


  //Primary Vertex
  const ZVertex *  trkZVertex();

  //Number of vertices
  int trkNZVertex();

  //Number of good (cal>=12) vertices
  int trkNGoodZVertex();

  //Z of Primary Vertex
  double trkZVertexZ();

  //Ht of Primary Vertex
  double trkZVertexHt();

  //Number of tracks associated with the primary Vertex
  int trkZVertexNTracks();

  //Silicon hit info -- disabled until output format agreed upon
  int trkSiHitPatterns(CdfTrackView::const_iterator& iter, SiExpected& 
siexphit, std::vector<SiExpected::Info>& results, int k);

  //******************* TOF information

  //Return associated TofMatches
  TofMatches* trkTofMatches(const CdfTrack* track);

  //Return associated TofMatch.  This is the old tof object but is still
  //present in old production
  const TofMatch* trkOldTofMatch(const CdfTrack* track);

}//tracking Variables

#endif