////////////////////////////////////////////////////////////////////
//
//     File:  TrackingVariables.cc
//     Authors:  Heather Gerberich
//     Description: High level analysis tracking variables defined as functions
//     Created: 10/30/2000
//
//     Revision: 0.2
//
////////////////////////////////////////////////////////////////////
//august 29, 2001  change return values to float from double

#include <cstdio>
//#include <math.h>
#include "HighLevelObjects/TrackingVariables.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "TrackingObjects/Tracks/SimpleReconstructedTrack.hh"
#include "CalorGeometry/CalConstants.hh"
#include "Trajectory/Angle.hh"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Geometry/Vector3D.h"
#include "CLHEP/Geometry/Point3D.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "linalg/CT_Helix.hh"
#include "TrackingObjects/Storable/CdfTrackColl.hh"
#include "TrackingObjects/CT_Track/CT_Residual.hh"
#include "TrackingObjects/CT_Track/CT_Hit.hh"
#include "TrackingObjects/CT_Track/CT_HitLink.hh"
#include "CotGeometry/CT_Configuration.hh"
#include "VertexObjects/Beamline.hh"
#include "TrackingSI/Utils/SiExpected.hh"
#include "TrackingObjects/SiData/SiHitSet.hh"
#include "TrackingObjects/SiData/SiHit.hh"
#include "TrackingCT/CT_ReFit.hh"


namespace TrackingVariables {


  //Returns pointer to 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){
  //int testInt = 0;
  CdfTrack_clnk childLnk = track->child();
  const CdfTrack * childPtr = childLnk.operator->();
  const CdfTrack * currentPtr =  track;
  //std::cout<<"Before Best Fit while loop "<<testInt;
  while(childPtr!=0)
    {
      //std::cout<<"Best Fit while loop "<<testInt;
      //testInt++;
      currentPtr = childPtr;
      childLnk = currentPtr->child();
      childPtr = childLnk.operator->();
    }
  //std::cout<<"After Best Fit while loop "<<testInt;
  return currentPtr;
}
  


  //Returns pointer to CdfTrack which contains only COT information
  //if no Cot Only Track exists, then returns original track
  //for now assumes incoming track is valid

const CdfTrack * trkCotOnly(const CdfTrack * track){

  //  int testInt = 0;

  const CdfTrack::TrackingAlgorithm algorithm=(*track).algorithm();
  if(algorithm.value()==CdfTrack::CotStandAloneAlg) return track;
  CdfTrack_clnk parentLnk = track->parent();
  const CdfTrack * parentPtr = parentLnk.operator->();
  const CdfTrack * currentPtr = track;
  //std::cout<<"before COT Only while loop "<<testInt<<std::endl;
  while(parentPtr!=0)
    {
      //std::cout<<"In CotOnly While Loop "<<testInt<<std::endl;
      //testInt++;
      currentPtr = parentPtr;
      parentLnk =currentPtr->parent();
      parentPtr = parentLnk.operator->();
    }
  //std::cout<<"after COT Only while loop "<<testInt<<std::endl;
    const CdfTrack::TrackingAlgorithm finalAlgorithm = (*currentPtr).algorithm();
    if(finalAlgorithm.value()==CdfTrack::CotStandAloneAlg) return currentPtr;
    errlog.setSubroutine("TrackingVariables::trkCotOnly");
    ERRLOG(ELwarning, "No CotStandAloneAlg") 
      << "Error:  No CotStandAloneAlg, using original track" << endmsg;
    return track; 
}



  //in following functions, if trkFlag = 0, use best fit
  //in following functions, if trkFlag = 1, use COT 


  //Return pointer to "Best Fit" track if trkFlag is 0
  //or pointer to "COT Track" if trkFlag is 1
  //Default is 2

const CdfTrack * trkGetDesiredTrack(const CdfTrack * track, int trkFlag){
  if(trkFlag == 0) return trkBestFit(track);
  else if (trkFlag == 1) return trkCotOnly(track);
  else return track;  //use default track if trkFlag is neither 0 or 1
}


  // Momentum 3-Vector

Hep3Vector trkMomentum(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->momentum();
}


  // Track Px

float trkPx(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  return (ourTrack->momentum()).x();
}


  // Track Py

float trkPy(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return (ourTrack->momentum()).y();
}


  // Track Pz

float trkPz(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  return (ourTrack->momentum()).z();
}


  // Track Pt

float trkPt(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->pt();
}


  // Track P

float trkP(const CdfTrack * track, int trkFlag){
  Hep3Vector ourMomentum = trkMomentum(track, trkFlag);
  return ourMomentum.mag();
}

  // Track Curvature (1/2r) in cm.  (+ for CounterClockwise Track; - for Clockwise Track)

float trkCurvature(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->curvature();
}

  // track direction at point of closest approach to the Z-axis

Angle trkPhi0(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->phi0();
}


  // track Lamda = Cot(Theta) where Theta is the Polar Angle of the track

float trkLambda(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->lambda();
}


  // track pseudorapidity

float trkEta0(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->pseudoRapidity();
}


  // track Charge

float trkCharge(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->charge();  
}


  // track d0 (impact parameter) in cm (+ sign for negative Lz, - sign for positive Lz) 

float trkD0(const CdfTrack * track, int trkFlag){
  //const CdfTrack::TrackingAlgorithm algorithm=(*track).algorithm();
  //algorithm.print(std::cout);std::cout<<" ";
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  //const CdfTrack::TrackingAlgorithm algorithmNew=(*ourTrack).algorithm();
  //algorithmNew.print(std::cout); std::cout<<" Best Fit"<<std::endl;
  return ourTrack->d0();
}


  // track X0

float trkX0(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return -1*trkD0(track, trkFlag)*sin(trkPhi0(track, trkFlag));
}
  

  // track Y0

float trkY0(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return trkD0(track, trkFlag)*cos(trkPhi0(track, trkFlag));
} 


  // track Z0

float trkZ0(const CdfTrack * track, int trkFlag){
  //const CdfTrack::TrackingAlgorithm algorithm=(*track).algorithm();
  //algorithm.print(std::cout);std::cout<<" ";
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  //const CdfTrack::TrackingAlgorithm algorithmNew=(*ourTrack).algorithm();
  //algorithmNew.print(std::cout); std::cout<<" COT Only"<<std::endl;
  return ourTrack->z0();
}

  // corrected D0 to a Beamline (comes from CotBeam or SvxBeam)
  double trkCorrectedD0(const CdfTrack* track, Beamline& beamline) {
    if (! beamline.isValid()) {
      errlog.setSubroutine("TrackingVariables::trkCorrectedD0");
      ERRLOG(ELinfo, "Beamline not valid") << "Provided beamline not valid; defaulting to origin"
					      << endmsg;
      return track->d0();
    }
    // this is a slight cheat but the inaccuracy introduced is negligible
    HepPoint3D beamspot = beamline.position(track->z0());
    
    return trkCorrectedD0(track, beamspot.x(), beamspot.y());
  }
  

  // corrected D0 to an arbitrary point
double trkCorrectedD0(const CdfTrack* track, double vx, double vy) {
  double uncord0 = track->d0();
  // "curvature" is in fact half-curvature
  double radius = fabs(0.5/track->curvature());
  double phi0 = track->phi0();
  // (x1, y1): coordinates of closest approach to origin
  double x1 = uncord0*-sin(phi0);
  double y1 = uncord0*cos(phi0);
  int charge = (track->curvature() < 0) ? -1 : 1;
  // (xc, yc): coordinates of helix center
  double xc = (charge*uncord0+radius)/(charge*uncord0)*x1;
  double yc = (charge*uncord0+radius)/(charge*uncord0)*y1;

  // rho: distance from helix center to (vx,vy)
  double rho = sqrt(pow(xc-vx,2)+pow(yc-vy,2));
  // to match CDF convention, multiply by charge
  return charge*(rho-radius);
}


  // Chi-Squared

float trkChi2(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  // Stephen Levy; Sep 13, 2004
  // The Tracking webpage says chiSquared() is a temporary result of a
  // SimpleReconstructedTrack that should not be used.
  //return ourTrack->chiSquared();
  return ourTrack->chi2();
}


  // Number of Degrees of Freedom

int trkNumDOF(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->numDegreesOfFreedom();
}

  //Number of hits in specified superlayer

int trkNumHitsInSL(int superLayer, const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  return ourTrack->numCTHits(superLayer);
}

  
  // Number of COT Axial Superlayers hit by track

int trkNumCOT_ASL(const CdfTrack * track, int trkFlag, int minHits){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  int numASL = 0;
  for (int i = 1; i<8; i+=2){
    if (ourTrack->numCTHits(i) >= minHits) numASL++;
  }
  return numASL;
}


  //Number of COT Stereo Superlayers hit by track

int trkNumCOT_SSL(const CdfTrack * track, int trkFlag, int minHits){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  int numSSL = 0;
  for (int i = 0; i<8; i+=2){
    if (ourTrack->numCTHits(i) >= minHits) numSSL++;
  }
  return numSSL;
}


  //Number of Silicon Hits on Track

int trkNumSVX_Hits(const CdfTrack * track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  return ourTrack->numSIHits();
}

  // "Bitmap" for hit pattern
  int trkSLBitmap(const CdfTrack * track, int trkFlag) {
    int retval = 0;
    for (int i = 0; i < 8; i++) {
      retval += ((trkNumHitsInSL(i, track, trkFlag) & 0xf) << i*4);
    }
    return retval;
  }
  
  int trkId(const CdfTrack * track, int trkFlag){
    const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
    return (ourTrack->id()).value();
  }


//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)  {
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  HepSymMatrix err=ourTrack->getHelixFit().getErrorMatrix();
   HepSymMatrix cov(5);
   for (int i=0; i < 5; i++) {
      for (int j=i; j < 5; j++) {
         if (err[i][i] <= 0 || err[j][j]<=0 ) {
            cov[i][j] = 999.;
         }
         else {
            cov[i][j] = err[i][j]/sqrt(err[i][i]*err[j][j]);
         }
      }
   }
  return cov;
}

  /*
HepSymMatrix trkCovMtrx(const CdfTrack * track, int trkFlag)  {
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  HepSymMatrix cov(5);  //Create 5X5 correlation matrix Pt=1, d0=2, phi0=3, z0=4, Lambda=5
  SimpleReconstructedTrack::TrackParameter Curve = SimpleReconstructedTrack::CURVATURE;  
  SimpleReconstructedTrack::TrackParameter D = SimpleReconstructedTrack::D0;
  SimpleReconstructedTrack::TrackParameter Phi = SimpleReconstructedTrack::PHI0;
  SimpleReconstructedTrack::TrackParameter Z = SimpleReconstructedTrack::Z0;  
  SimpleReconstructedTrack::TrackParameter Lambda = SimpleReconstructedTrack::LAMBDA;
  float k;
  if(ourTrack->curvature() > 0) k = 0.0021159;
  else k = -0.0021159;
  float sigmaCurvature = ourTrack->sigmaCurvature();
  float sigmaPt = sigmaCurvature*k/(ourTrack->curvature()*ourTrack->curvature());
  float sigmaD = ourTrack->sigmaD0();
  float sigmaPhi = ourTrack->sigmaPhi0();  
  float sigmaZ = ourTrack->sigmaZ0();
  float sigmaLambda = ourTrack->sigmaLambda();
  cov(1,1)=sigmaPt*sigmaPt; //sigmaPt^2
  cov(1,2)=-k*ourTrack->correlation(Curve, D)*sigmaCurvature*sigmaD/(ourTrack->curvature()*ourTrack->curvature()); //cov(Pt, D0) 
  cov(1,3)=-k*ourTrack->correlation(Curve, Phi)*sigmaCurvature*sigmaPhi/(ourTrack->curvature()*ourTrack->curvature()); //cov(Pt, Phi0)
  cov(1,4)=-k*ourTrack->correlation(Curve, Z)*sigmaCurvature*sigmaZ/(ourTrack->curvature()*ourTrack->curvature()); //cov(Pt, Z0)  
  cov(1,5)=-k*ourTrack->correlation(Curve, Lambda)*sigmaCurvature*sigmaLambda/(ourTrack->curvature()*ourTrack->curvature()); //cov(Pt, Lambda)
  cov(2,2)=sigmaD*sigmaD; //sigma(d0)^2
  cov(2,3)=ourTrack->correlation(D,Phi)*sigmaD*sigmaPhi; //cov(D,Phi)
  cov(2,4)=ourTrack->correlation(D,Z)*sigmaD*sigmaZ; //cov(D,Z)
  cov(2,5)=ourTrack->correlation(D,Lambda)*sigmaD*sigmaLambda; //cov(D,Lambda)
  cov(3,3)=sigmaPhi*sigmaPhi;  //Sigma(Phi)^2
  cov(3,4)=ourTrack->correlation(Phi,Z)*sigmaPhi*sigmaZ;  //cov(Phi, Z)
  cov(3,5)=ourTrack->correlation(Phi,Lambda)*sigmaPhi*sigmaLambda;  //Cov(Phi, Lambda)
  cov(4,4)=sigmaZ*sigmaZ;  //Sigma(Z)^2
  cov(4,5)=ourTrack->correlation(Z,Lambda)*sigmaZ*sigmaLambda;  //Cov(Z,Lambda)  
  cov(5,5)=sigmaLambda*sigmaLambda;  //Sigma(Lambda)^2
  return cov;
}
  */

  //Covariance Matrix element (Pt, Pt)

float trkCovMtrxPtPt(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(1,1);
}

  //Covariance Matrix element (Pt, D0)

float trkCovMtrxPtD0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(1,2);
}

  //Covariance Matrix element (Pt, Phi0)

float trkCovMtrxPtPhi0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(1,3);
}

  //Covariance Matrix element (Pt, Z0)

float trkCovMtrxPtZ0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(1,4);
}

  //Covariance Matrix element (Pt, Lambda)

float trkCovMtrxPtLambda(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(1,5);
}


  //Covariance Matrix element (D0, D0)

float trkCovMtrxD0D0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(2,2);
}

  //Covariance Matrix element (D0, Phi0)

float trkCovMtrxD0Phi0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(2,3);
}

  //Covariance Matrix element (D0, Z0)

float trkCovMtrxD0Z0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(2,4);
}

  //Covariance Matrix element (D0, Lambda)

float trkCovMtrxD0Lambda(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(2,5);
}

  //Covariance Matrix element (Phi0, Phi0)

float trkCovMtrxPhi0Phi0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(3,3);
}

  //Covariance Matrix element (Phi, Z0)

float trkCovMtrxPhi0Z0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(3,4);
}

  //Covariance Matrix element (Phi0, Lambda)

float trkCovMtrxPhi0Lambda(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(3,5);
}

  //Covariance Matrix element (Z0, Z0)

float trkCovMtrxZ0Z0(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(4,4);
}

  //Covariance Matrix element (Z0, Lambda)

float trkCovMtrxZ0Lambda(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(4,5);
}

  //Covariance Matrix element (Lambda, Lambda)

float trkCovMtrxLambdaLambda(const CdfTrack * track, int trkFlag){
  return trkCovMtrx(track, trkFlag)(5,5);
}

  //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){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  float c = ourTrack->curvature();
  float d = (ourTrack->d0());
  if( radius>fabs((1/c)+d)|| radius<fabs(d) ) return false; 
  else return true;
}

  //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. x, y, z

HepPoint3D trkExtPos(float radius, const CdfTrack * track, int trkFlag){
  if (!trkRIsValid(radius, track, trkFlag)) 
    {
      HepPoint3D DNE(-999999., -999999., -999999.);
      return DNE;
    }
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag); 
  const Helix & ourHelix = ourTrack->getTrajectory();  //trajectory of track
  float PathLength = ourHelix.getPathLengthAtRhoEquals(radius);  //length of path of track to given radius
  return ourHelix.getPosition(PathLength);  //return position of track when it has traveled distance PathLength 
}


  //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 x slope, y slope, z slope

HepVector3D trkExtSlope(float radius, const CdfTrack * track, int trkFlag){
  if (!trkRIsValid(radius, track, trkFlag))
    {
      HepVector3D DNE(-999999., -999999., -999999.);
      return DNE;
    }
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const Helix & ourHelix = ourTrack->getTrajectory();  //trajectory of track
  float PathLength = ourHelix.getPathLengthAtRhoEquals(radius);  //length of path of track to given radius
  return ourHelix.getDirection(PathLength);  //return position of track when it has traveled distance PathLength 
}


  //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.

float trkExtPathLength(float radius, const CdfTrack * track, int trkFlag){ 
  if (!trkRIsValid(radius, track, trkFlag)) return -999999.;
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const Helix & ourHelix = ourTrack->getTrajectory();
  return ourHelix.getPathLengthAtRhoEquals(radius);
}


  // 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){
  HepPoint3D Position = trkExtPos(radius, track, trkFlag);
  return Position.x();
}


  // 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 -999999999 for y

float trkExtY(float radius, const CdfTrack * track, int trkFlag){
  HepPoint3D Position = trkExtPos(radius, track, trkFlag);
  return Position.y();
}


  // 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){
  HepPoint3D Position = trkExtPos(radius, track, trkFlag);
  return Position.z();
}

  //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 x slope

float trkExtXSlope(float radius, const CdfTrack * track, int trkFlag){
  HepVector3D Slope = trkExtSlope(radius, track, trkFlag);
  return Slope.x();
}


  //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 y slope

float trkExtYSlope(float radius, const CdfTrack * track, int trkFlag){
  HepVector3D Slope = trkExtSlope(radius, track, trkFlag);
  return Slope.y();
}


  //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 z slope

float trkExtZSlope(float radius, const CdfTrack * track, int trkFlag){
  HepVector3D Slope = trkExtSlope(radius, track, trkFlag);
  return Slope.z();
}

  //returns a pointer to a CT_Residual object for the specified superLayer
 
const CT_Residual * trkResidual(int superLayer, const CdfTrack*track, int trkFlag){
    const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
    return ourTrack->beginCTResiduals(superLayer);
}

  //Delta D residual for a specified superlayer and wire
  //deltaD is hitD - |trackY|
float trkDeltaDRes(int superLayer, int wire, const CdfTrack*track, int trkFlag)
{
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const CT_Residual * resid = ourTrack->beginCTResiduals(superLayer);
  if (resid[wire].isValid()) return resid[wire].deltaD();
  else return -999999.;
}

  //Delta Y residual for a specified superlayer and wire 
  //DeltaY is hitD * (trackY / |trackY|) - trackY
float trkDeltaYRes(int superLayer, int wire, const CdfTrack*track, int trkFlag){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const CT_Residual * resid = ourTrack->beginCTResiduals(superLayer);
  if (resid[wire].isValid()) return resid[wire].deltaY();
  else return -999999.;
}

  //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){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const CT_Residual * resid = ourTrack->beginCTResiduals(superLayer);
  if (resid[wire].isValid()) return resid[wire].hitD();
  else return -999999.;
}

  //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){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const CT_Residual * resid = ourTrack->beginCTResiduals(superLayer);
  if (resid[wire].isValid()) return resid[wire].trackY();
  else return -999999.;
}

  //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){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const CT_Residual * resid = ourTrack->beginCTResiduals(superLayer);
  if (resid[wire].isValid()) return resid[wire].sinAspect();
  else return -999999.;
}

  //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){
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const CT_Residual * resid = ourTrack->beginCTResiduals(superLayer);
  if (resid[wire].isValid()) return resid[wire].hitT();
  else return -999999.;
}

//Hit width of a specified superlayer and wire with a valid residual
int trkHitWidth(int superLayer, int wire, const CdfTrack*track, int trkFlag)
{
   const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
   const CT_Residual * resid = ourTrack->beginCTResiduals(superLayer);
   const CT_HitLink * hits = ourTrack->beginCTHitBlock(superLayer);
   if (resid[wire].isValid()) return hits[wire].hit().width();
   else return -999999;
}

//Cell number a specified superlayer and wire with a valid residual
int trkHitCellNum(int superLayer, int wire, const CdfTrack*track, int trkFlag)
{    
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const CT_Residual * resid = ourTrack->beginCTResiduals(superLayer);
  const CT_HitLink * hits = ourTrack->beginCTHitBlock(superLayer);
  if (resid[wire].isValid()) return hits[wire].cell();
  else return -999999;
}



//returns the inverse of beta (c/v) using residuals for calculation
//and the error on the inverse of beta
//both passed by reference (Kotwal)
void trkMinChi2BetaInv(float & MinChi2BetaInv, float & DeltaBetaInv,const CdfTrack * track, int trkFlag)
{
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  const float VDrift = 52.4; // drift velocity in microns / ns
  const float SpeedLight = 300000.0;  // speed of light in microns / ns 
  const float ResidueError = 0.02; // error on residue in cm
  const float WireSpacing = 0.762; // radial separation of wires in cm
  float BetaInv = 0.0; // inverse of track beta (velocity / speed of light)
  
  const int nChi2Points=600;
  const float BetaInvLow = -2.0;
  const float BetaInvHigh = 4.0;
  const float BetaInvStep = (BetaInvHigh-BetaInvLow)/float(nChi2Points);
  float Chi2Points[nChi2Points];
  
  for(int npts = 0; npts<nChi2Points; ++npts) {Chi2Points[npts] = 0.0;} 
  int NWires = 0; //count wires contributing to calculation
  
  float cottheta = ourTrack->lambda();
  float sintheta = 1.0/sqrt(1.0+cottheta*cottheta);  
  
  const CT_Configuration& config = CT_Configuration::reference();
  int numSL = config.numSuperLayers();  
   
  for ( int sl = 0 ; sl < numSL ; ++sl )
    {
      int numWires = config.numWires( sl ); //number of wires in superlayer sl
      for ( int wire = 0 ; wire < numWires ; ++wire )
	{
	  float cwire = float(wire);  //cell wire
	  float wire_radius = 46.0 + 12.0*float(sl) + (cwire-5.5)*WireSpacing;
	  float curv = ourTrack->curvature();
	  float sin_psi = fabs(wire_radius*curv);
	  if(sin_psi>1.0) {
	    MinChi2BetaInv=-9999.;  //looping track; not interesting
	    DeltaBetaInv=-9999.;  //looping track; not interesting
	    return;
	  }
	  float psi = asin(sin_psi);
	  float arc_length = wire_radius*psi/sin_psi;
	  float tof_light =  arc_length*VDrift/SpeedLight/sintheta;
	  float deltad = trkDeltaDRes(sl, wire, ourTrack, trkFlag);
            if(fabs(deltad)<9.9) 
	      {
		deltad = deltad * 10000.0;
	      }
	  float hitd = trkHitDRes(sl, wire, ourTrack,trkFlag);
	    if (fabs(deltad)<1000.0 && hitd>0.2 && hitd<0.7) 
	      {
		for(int npts = 0; npts<nChi2Points; ++npts) 
		  {
		    BetaInv = BetaInvLow + BetaInvStep/2.0 + BetaInvStep * npts;
		    float tof = tof_light*BetaInv; // tof is in cm 
		    float chi = hitd - tof;
		    float tracky = trkTrackYRes(sl, wire, ourTrack, trkFlag);
		    if (tracky<0.0) {chi = -chi;}
		    chi = (chi - tracky)/ResidueError; 
		    Chi2Points[npts] = Chi2Points[npts] +  chi*chi;
		  }
		NWires = NWires + 1; //num of wires contributing to calculation
	      }
	}
    } //end of all wires loop

  if( NWires < 7) {
    MinChi2BetaInv=-9999.;  
    DeltaBetaInv=-9999.;  
    return;
  }

  float MinChi2 = 9999.9;
  MinChi2BetaInv = 9999.9;
  bool foundMin_npts = false ;
  int min_npts;
  for(int npts = 0; npts<nChi2Points; ++npts) 
    {
      Chi2Points[npts] = Chi2Points[npts] / float(NWires - 6) ;
      BetaInv = BetaInvLow + BetaInvStep/2.0 + BetaInvStep * npts;
      if (Chi2Points[npts]<MinChi2) 
	{
	  MinChi2 =  Chi2Points[npts];
	  MinChi2BetaInv = BetaInv;
	  foundMin_npts = true;
	  min_npts = npts;
	}
    }
 
  if (!foundMin_npts) {
    DeltaBetaInv = -9999.;
    return;
  }
  float DeltaChi2 = MinChi2 / float(NWires - 6) ;
  int npts_hi;
  for(int npts = 0; npts<nChi2Points; ++npts) {
    if (Chi2Points[npts]-Chi2Points[min_npts]<DeltaChi2) {
      npts_hi = npts;
    }
  }

  int npts_lo;
  for(int npts = nChi2Points; npts>0; npts--) {
    if (Chi2Points[npts]-Chi2Points[min_npts]<DeltaChi2) {
      npts_lo = npts;      
    }
  }

  DeltaBetaInv = BetaInvStep * float(npts_hi-npts_lo+1) / 2.0 ;

  return; 

}



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

  float hitwidth_array[96];  

  const CT_Configuration& config = CT_Configuration::reference();
  int numSL = config.numSuperLayers();  
   
  int nHitWidths=0;  //hit widths counter
  for ( int sl = 0 ; sl < numSL ; ++sl )
    {
      int numWires = config.numWires( sl ); //number of wires in superlayer sl
      for ( int wire = 0 ; wire < numWires ; ++wire )
	{
	  int hitwidth = trkHitWidth(sl, wire, ourTrack, trkFlag);
	  if (hitwidth>15.0 && hitwidth<200.0) {
             hitwidth_array[nHitWidths] = hitwidth;
             nHitWidths = nHitWidths + 1;
           }
	}
    }//done looping over all wires

  // sort hit widths in ascending order
  for (Int_t k=0; k<nHitWidths-1; k++) {
    Int_t minindex = k;                // minimal element index
    for (Int_t j=k+1; j<nHitWidths; j++){
      if (hitwidth_array[j] < hitwidth_array[minindex]) {
	minindex = j;                  // new minimum, store index
      }
    }
    Float_t temp=hitwidth_array[k];
    hitwidth_array[k]=hitwidth_array[minindex];  // swap min, k elements
           hitwidth_array[minindex]=temp;
  }

  if(nHitWidths==0) return -9999.;

  // compute truncated mean
  nHitWidths = int(0.7*float(nHitWidths));
  float trunc_hitwidth = 0.0;
  for (Int_t k=0; k<nHitWidths; k++) {
    trunc_hitwidth = trunc_hitwidth + hitwidth_array[k];
  }
  trunc_hitwidth = trunc_hitwidth / float(nHitWidths);
  
  return trunc_hitwidth;
}


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

  float cotTheta=trkLambda(ourTrack, trkFlag);

  float hitwidth_array[96];  

  const CT_Configuration& config = CT_Configuration::reference();
  int numSL = config.numSuperLayers();  
   
  int nHitWidths=0;  //hit widths counter
  for ( int sl = 0 ; sl < numSL ; ++sl )
    {
     int numWires = config.numWires( sl ); //number of wires in superlayer sl
      for ( int wire = 0 ; wire < numWires ; ++wire )
        {
          float hitwidth = trkHitWidth(sl, wire, ourTrack, trkFlag);
          float sinAspectAngle = trkSinAspectRes(sl, wire, ourTrack,trkFlag);
          if (hitwidth>15.0 && hitwidth<200.0) {
            hitwidth = hitwidth - 9.0;
            if ( fabs(sinAspectAngle)<0.99 ) {
              float plc = sqrt ( cotTheta*cotTheta + 1.0 / (1.0 -
sinAspectAngle*sinAspectAngle) ) ; 
              if (plc<5.0) {hitwidth = hitwidth / plc ;}
            }
            hitwidth_array[nHitWidths] = hitwidth ;
            nHitWidths = nHitWidths + 1; 
          }
        }
    }//done looping over all wires

  // sort hit widths in ascending order
  for (Int_t k=0; k<nHitWidths-1; k++) {
    Int_t minindex = k;                // minimal element index
    for (Int_t j=k+1; j<nHitWidths; j++){
      if (hitwidth_array[j] < hitwidth_array[minindex]) {
        minindex = j;                  // new minimum, store index
      }
    }
    Float_t temp=hitwidth_array[k];
    hitwidth_array[k]=hitwidth_array[minindex];  // swap min, k elements
           hitwidth_array[minindex]=temp;
  }

  if(nHitWidths==0) return -9999.;

  // compute truncated mean
  nHitWidths = int(0.7*float(nHitWidths));
  float trunc_hitwidth = 0.0;
  for (Int_t k=0; k<nHitWidths; k++) {
    trunc_hitwidth = trunc_hitwidth + hitwidth_array[k];
  }
  trunc_hitwidth = trunc_hitwidth / float(nHitWidths);
  
  return trunc_hitwidth;
}

// 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)
{
  const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);
  //const float VDrift = 52.4; // drift velocity in microns / ns
  //const float SpeedLight = 300000.0;  // speed of light in microns / ns 
  const float ResidueError = 0.02; // error on residue in cm
  const float WireSpacing = 0.762; // radial separation of wires in cm
  float t0 = 0.0; // track t0 (production time)
  
  const int nt0Points=100;
  const float t0Low = -0.05;
  const float t0High = 0.05;
  const float t0Step = (t0High-t0Low)/float(nt0Points);
  float Chi2Points[nt0Points];

  // initialize outputs  
  MinChi2t0 = -9999.0;  // t0 value at minimum chi2
  Deltat0 = -9999.0;    // error on t0
  Chi2t0 = -9999.0;     // value of chi2/dof at minimum
  NWirest0 = -9999;     // number of residuals used

  for(int nt0 = 0; nt0<nt0Points; ++nt0) {Chi2Points[nt0] = 0.0;} 
  int NWires = 0; //count wires contributing to calculation
  
  float cottheta = ourTrack->lambda();
  float sintheta = 1.0/sqrt(1.0+cottheta*cottheta);  
  
  const CT_Configuration& config = CT_Configuration::reference();
  int numSL = config.numSuperLayers();  
   
  for ( int sl = 0 ; sl < numSL ; ++sl )
    {
      int numWires = config.numWires( sl ); //number of wires in superlayer sl
      for ( int wire = 0 ; wire < numWires ; ++wire )
        {
          float cwire = float(wire);  //cell wire
          float wire_radius = 46.0 + 12.0*float(sl) + (cwire-5.5)*WireSpacing;
          float curv = ourTrack->curvature();
          float sin_psi = fabs(wire_radius*curv);
          if(sin_psi>1.0) return;  //looping track; not interesting
          float psi = asin(sin_psi);
          //float arc_length = wire_radius*psi/sin_psi;
          //float tof_light =  arc_length*VDrift/SpeedLight/sintheta;
          float deltad = trkDeltaDRes(sl, wire, ourTrack, trkFlag);
          if(fabs(deltad)<9.9) 
            {
              deltad = deltad * 10000.0;
            }
          float hitd = trkHitDRes(sl, wire, ourTrack,trkFlag);
            if (fabs(deltad)<1000.0 && hitd>0.2 && hitd<0.7) {
                for(int nt0 = 0; nt0<nt0Points; ++nt0) {
                    t0 = t0Low + t0Step * nt0;
                    float chi = hitd - t0;
                    float tracky = trkTrackYRes(sl, wire, ourTrack, trkFlag);
                    if (tracky<0.0) {chi = -chi;}
                    chi = (chi - tracky)/ResidueError; 
                    Chi2Points[nt0] += chi*chi;
                  }
                NWires = NWires + 1; //num of wires contributing to calculation
            }
        }
    } //end of all wires loop

  NWirest0 = NWires;
  if( NWires < 7) return;

  Chi2t0 = 9999.9;
  MinChi2t0 = 9999.9;
  int min_nt0;
  bool foundMin_nt0 = false ;
  for (int nt0 = 0; nt0<nt0Points; ++nt0) {
      Chi2Points[nt0] = Chi2Points[nt0] / float(NWires - 6) ;
      t0 = t0Low + t0Step * nt0;
      if (Chi2Points[nt0]<Chi2t0) {
          Chi2t0 =  Chi2Points[nt0];
          MinChi2t0 = t0;
          min_nt0 = nt0;
          foundMin_nt0 = true;
      }
    }

  if (!foundMin_nt0) return;

  float DeltaChi2 = Chi2t0 / float(NWires - 6) ;
  int nt0_hi;
  for(int nt0 = 0; nt0<nt0Points; ++nt0) {
    if (Chi2Points[nt0]-Chi2Points[min_nt0]<DeltaChi2) {
      nt0_hi = nt0;
    }
  }

  int nt0_lo;
  for(int nt0 = nt0Points; nt0>0; nt0--) {
    if (Chi2Points[nt0]-Chi2Points[min_nt0]<DeltaChi2) {
      nt0_lo = nt0;      
    }
  }

  Deltat0 = t0Step * float(nt0_hi-nt0_lo+1) / 2.0 ;

  return;
}

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

   const CdfTrack * ourTrack = trkGetDesiredTrack(track, trkFlag);

   const float VDrift = 52.4; // drift velocity in microns / ns
   const float SpeedLight = 300000.0;  // speed of light in microns / ns 
   const float ResidueError = 0.02; // error on residue in cm
   const float WireSpacing = 0.762; // radial separation of wires in cm
   float BetaInv = 0.0; // inverse of track beta (velocity / speed of light)
   float t0 = 0.0; // t0 for track (production time)

   const int max_loops = 3; // number of loops for increasing precision

   const int nChi2Points = 12;
   float BetaInvLow = -2.0;
   float BetaInvHigh = 4.0;
   float BetaInvStep = (BetaInvHigh-BetaInvLow)/float(nChi2Points);

   const int nt0Points = 8;
   float t0Low = -0.05;
   float t0High = 0.05;
   float t0Step = (t0High-t0Low)/float(nt0Points);

   float Chi2Points[nChi2Points][nt0Points];

   // initialize outputs  
   MinChi22Dt0 = -9999.0;        // t0 value at minimum chi2
   Delta2Dt0 = -9999.0;          // error on t0
   MinChi22DBetaInv = -9999.0;   // 1/Beta value at minimum chi2
   Delta2DBetaInv = -9999.0;     // error on 1/Beta
   Cov2D2DBetaInvt0 =  -9999.0;  // covariance between 1/Beta and t0
   Chi22DBetaInvt0 = -9999.0;    // value of chi2/dof at minimum
   NWiresBetaInvt0 = -9999;      // number of residuals used

  for(int nt0 = 0; nt0<nt0Points; ++nt0) { 
      for(int npts = 0; npts<nChi2Points; ++npts) {Chi2Points[npts][nt0] = 0.0;} }
  
  float cottheta = ourTrack->lambda();
  float sintheta = 1.0/sqrt(1.0+cottheta*cottheta);  
  
  const CT_Configuration& config = CT_Configuration::reference();
  int numSL = config.numSuperLayers();  
   
  float Sum;
  float t0Sum;
  float t0Sum2;
  float BetaInvSum;
  float BetaInvSum2;
  float t0BetaInvSum;
  float prob;
  float chi2;

  for ( int nloop = 0; nloop < max_loops; ++nloop ) {

  int NWires = 0; // count wires contributing to calculation

  for ( int sl = 0 ; sl < numSL ; ++sl ) {
      int numWires = config.numWires( sl ); //number of wires in superlayer sl
      for ( int wire = 0 ; wire < numWires ; ++wire )  {
        float cwire = float(wire);  //cell wire
        float wire_radius = 46.0 + 12.0*float(sl) + (cwire-5.5)*WireSpacing;
        float curv = ourTrack->curvature();
        float sin_psi = fabs(wire_radius*curv);
        if(sin_psi>1.0) return;  //looping track; not interesting
        float psi = asin(sin_psi);
        float arc_length = wire_radius*psi/sin_psi;
        float tof_light =  arc_length*VDrift/SpeedLight/sintheta;
        float deltad = trkDeltaDRes(sl, wire, ourTrack, trkFlag);
        if(fabs(deltad)<9.9)
          {
            deltad = deltad * 10000.0;
          }
        float hitd = trkHitDRes(sl, wire, ourTrack,trkFlag);
        if (fabs(deltad)<1000.0 && hitd>0.2 && hitd<0.7) {
             for(int nt0 = 0; nt0<nt0Points; ++nt0) { 
              t0 = t0Low + t0Step * nt0;   
              for(int npts = 0; npts<nChi2Points; ++npts) {
                BetaInv = BetaInvLow + BetaInvStep * npts;   
                float tof = tof_light*BetaInv; // tof is in cm 
                float chi = hitd - tof - t0;
                float tracky = trkTrackYRes(sl, wire, ourTrack, trkFlag);
                if (tracky<0.0) {chi = -chi;}
                chi = (chi - tracky)/ResidueError; 
                Chi2Points[npts][nt0] +=  chi*chi;
              }
             }
             NWires = NWires + 1;
        }
      }
  } // end wire and sl loop

  NWiresBetaInvt0 = NWires;
  if( NWires < 8) return;
  
  MinChi22Dt0 = 9999.0;  
  Delta2Dt0 = 9999.0;    
  MinChi22DBetaInv = 9999.0; 
  Delta2DBetaInv = 9999.0;   
  Cov2D2DBetaInvt0 =  9999.0; 
  Chi22DBetaInvt0 = 9999.0;   
  
  bool foundMin_nt0_npts = false ;
  
  for(int nt0 = 0; nt0<nt0Points; ++nt0) { 
    t0 = t0Low + t0Step * nt0;   
    for(int npts = 0; npts<nChi2Points; ++npts) {
      BetaInv = BetaInvLow + BetaInvStep * npts;
      Chi2Points[npts][nt0] = Chi2Points[npts][nt0] / float(NWires - 7) ;
      if (Chi2Points[npts][nt0]<Chi22DBetaInvt0) {
	Chi22DBetaInvt0 =  Chi2Points[npts][nt0];
	MinChi22DBetaInv = BetaInv;
	MinChi22Dt0 = t0;
	foundMin_nt0_npts = true;
      }
    }
  }
  if (!foundMin_nt0_npts) return;
  float DeltaChi2 = Chi22DBetaInvt0 / float(NWires - 7) ;
  
  Sum = 0.0;
  t0Sum = 0.0;
  t0Sum2 = 0.0;
  BetaInvSum = 0.0;
  BetaInvSum2 = 0.0;
  t0BetaInvSum = 0.0;
  int nfound = 0;
  
  for(int nt0 = 0; nt0<nt0Points; ++nt0) { 
    t0 = t0Low + t0Step * nt0;   
    for(int npts = 0; npts<nChi2Points; ++npts) {
      BetaInv = BetaInvLow + BetaInvStep * npts;
      chi2 = Chi2Points[npts][nt0] - Chi22DBetaInvt0 ; // subtract off minimum chi2
      if (chi2<3.0*DeltaChi2) {
	prob = exp ( - chi2 / 2.0 );    // 2-dof gamma distribution for chi2
	nfound++;
	Sum += prob; 
	t0Sum += t0 * prob;
	t0Sum2 += t0 * t0 * prob;
	BetaInvSum += BetaInv * prob;
	BetaInvSum2 += BetaInv * BetaInv * prob;
	t0BetaInvSum += t0 * BetaInv * prob;              
      }
    }
  }
  
  if(Sum>0. && nfound>1){
    MinChi22Dt0 = t0Sum / Sum;
    float meant02 = t0Sum2 / Sum;
    MinChi22DBetaInv = BetaInvSum / Sum;
    float meanBetaInv2 = BetaInvSum2 / Sum;
    float meant0BetaInv = t0BetaInvSum / Sum;      
    if ((meant02 -  MinChi22Dt0* MinChi22Dt0)>0. && 
        ( meanBetaInv2 - MinChi22DBetaInv*MinChi22DBetaInv )>0.) {
      Delta2Dt0 = sqrt ( meant02 -  MinChi22Dt0* MinChi22Dt0 ) ;
      Delta2DBetaInv = sqrt ( meanBetaInv2 - MinChi22DBetaInv*MinChi22DBetaInv ) ;
    }
    else {
      Delta2Dt0 = t0Step;
      Delta2DBetaInv = BetaInvStep;
    }
    Cov2D2DBetaInvt0 = meant0BetaInv -  MinChi22Dt0*MinChi22DBetaInv ;
  }
  else{
    Delta2Dt0 = t0Step;
    Delta2DBetaInv = BetaInvStep;
  }

  BetaInvLow = MinChi22DBetaInv - BetaInvStep ;
  BetaInvHigh = MinChi22DBetaInv + BetaInvStep ;
  BetaInvStep = 2.0 * BetaInvStep / float(nChi2Points) ;
  
  t0Low = MinChi22Dt0 - t0Step;
  t0High = MinChi22Dt0 + t0Step;
  t0Step =  2.0 * t0Step / float(nt0Points);
  
  }  

}

int trkSiHitPatterns(CdfTrackView::const_iterator& iter, SiExpected& 
siexphit,std::vector<SiExpected::Info>& results,int o) {

  // disable until format is agreed upon
  return 0;

  int layer = 999;
  int errlay1 = 0;
  int errlay2 = 0;
  int errlay3 = 0;
  int errlay4 = 0;
  int errlay5 = 0;
  int intlay1 = 0;
  int intlay2 = 0;
  int intlay3 = 0;
  int intlay4 = 0;
  int intlay5 = 0;
  int il1=0;
  int il2=0;
  int il3=0;
  int il4=0;
  int il5=0;

  if ( siexphit.findExpected( **iter, results ) ){
    for ( std::vector<SiExpected::Info>::iterator ii=results.begin(),iiend 
= results.end() ; ii != iiend ; ++ii ){     
      if ( (*ii).isPhi ){            //only count phi hits for now
	if ( (*ii). digiCode.getLayer()!=0 ){     //no layer 00
	  if ( (*ii). digiCode.getLayer()<=5 ){   //no ISL
	    if ( (*ii).hitByTrack ){
	      
	      layer = (*ii).digiCode.getLayer();
	      if ( (*ii).integrated ){
		  if (layer==1) intlay1 += 1;
		  if (layer==2) intlay2 += 1;
		  if (layer==3) intlay3 += 1;
		  if (layer==4) intlay4 += 1;
		  if (layer==5) intlay5 += 1;
		  if ( (*ii).readoutError ){
		      if (layer==1) errlay1 += 1;
		      if (layer==2) errlay2 += 1;
		      if (layer==3) errlay3 += 1;
		      if (layer==4) errlay4 += 1;
		      if (layer==5) errlay5 += 1;
		    }
	      }
	    }     // hit by track
	  }     //no ISL
	}      //forget layer 00 
      }      // only count phi side for now
      
      
    }  // loop over Info objects
    
  }

  for(CdfTrack::SiHitIterator siHit=(*iter)->beginSIHits();siHit !=(*iter)->endSIHits(); ++siHit){
     int il =(*siHit)->getHalfLad()->getDetectorCode().getLayer();
     if ((*siHit)->pSideHit()){
       if (il==1) il1 +=1;
       if (il==2) il2 +=1;
       if (il==3) il3 +=1;
       if (il==4) il4 +=1;
       if (il==5) il5 +=1;
     }
   }
   int pattern=il5*10000+il4*1000+il3*100+il2*10+il1;
   int errlay=errlay5*10000+errlay4*1000+errlay3*100+errlay2*10+errlay1;
   int intlay=intlay5*10000+intlay4*1000+intlay3*100+intlay2*10+intlay1;

   if (o==1) return pattern;
   else if (o==2) return errlay;
   else return intlay;
   
}


int trkTransitionPattern(const CdfTrack& track) {
  static CT_ReFit refit;
  refit.restore(track);
  int max_seq = 0;
  int n_trans = 0;
  
  if (!refit.ok()) return 0;
  
  // loop over COT layers
  int ilr = 0;
  int q = 0, seq  = 0;
  
  const CT_ReFit::LayerData* ldat = refit.layer_data();
  for (int sl = 0; sl < 8; ++sl) {
    for (int wr = 0; wr < 12; ++wr, ++ilr) {
      if (ldat[ilr].link.is_valid()) {
	
	// actual hit
	CT_Hit hit = ldat[ilr].link.hit();
	
	// residual
	float res = ldat[ilr].deltaY();
	if( q == 0 ) q = ( res >= 0 ) ? +1 : -1; // initialization
	
	if( res*q < 0 ){ // transition happend
	  n_trans++;
	  if(max_seq < seq) max_seq = seq;
	  seq = 0;
	  q = -q;
	}
	
      } // if layer data link is valid
      else {
      } // if layer data link is not valid
      seq++;
    } // next wire
  } // next sl
  

  if(max_seq < seq) max_seq = seq;

  int ret = (n_trans<<8) | max_seq;

  return ret;

} // end function

} // namespace TrackingVariables