//
//  TofGeometry.cc _ Definitions for TOF detector geometry
//
//  27 June 2002 _ m jones
//

#include <iostream.h>
#include <assert.h>

#include "TofGeometry.hh"

using std::cout;
using std::endl;

// ------------------------------------------------------------------------

TofHitBar::TofHitBar(int ibar) {
  _barNumber = ibar;
}

void TofHitBar::setInOut(const TVector3 &in,const TVector3 &out) {
  _in = in;
  _out = out;
}

void TofHitBar::setDirInOut(const TVector3 &dirin,
                            const TVector3 &dirout) {
  _dirin = dirin;
  _dirout = dirout;
}

void TofHitBar::setArcLengthInOut(double ss_in, double ss_out) {
  _ss_in = ss_in;
  _ss_out = ss_out;
}

void TofHitBar::setTimeInOut(double t_in, double t_out) {
  _t_in = t_in;
  _t_out = t_out;
}

void TofHitBar::setSurfaceInOut(int went_in, int went_out) {
  _went_in = went_in;
  _went_out = went_out;
}

void TofHitBar::setUVInOut(double u_in, double v_in,
                        double u_out, double v_out) {
  _u_in  = u_in;
  _v_in  = v_in;
  _u_out = u_out;
  _v_out = v_out;
}

void TofHitBar::setUErrorInOut(double err_u_in, double err_u_out) {
  _err_u_in  = err_u_in;
  _err_u_out = err_u_out;
}

TofBar::TofBar() {
  double width_top, width_bot, height, length, angle;

  width_top = 3.92;
  width_bot = 4.04;
  height = 4.00;
  length = 279.5;
  angle = atan((width_bot-width_top)/(2*height));

  //  Top surface
 
  _normal[0] = TVector3(0,1,0);
  _point[0] = TVector3(0,height/2,0);
  _axis[0] = TVector3(0,0,1);
  _length[0] = length;
  _width[0] = width_top;
  _dwdl[0] = 0;
 
  //  Bottom surface
 
  _normal[1] = TVector3(0,-1,0);
  _point[1] = TVector3(0,-height/2,0);
  _axis[1] = TVector3(0,0,1);
  _length[1] = length;
  _width[1] = width_bot;
  _dwdl[1] = 0;

  //  Left side

  _normal[2] = TVector3(-cos(angle),sin(angle),0);
  _point[2] = TVector3(-(width_top+width_bot)/4,0,0);
  _axis[2] = TVector3(0,0,1);
  _length[2] = length;
  _width[2] = height/cos(angle);
  _dwdl[2] = 0;

  //  Right side

  _normal[3] = TVector3(cos(angle),sin(angle),0);
  _point[3] = TVector3((width_top+width_bot)/4,0,0);
  _axis[3] = TVector3(0,0,1);
  _length[3] = length;
  _width[3] = height/cos(angle);
  _dwdl[3] = 0;

  //  Front end

  _normal[4] = TVector3(0,0,-1);
  _point[4] = TVector3(0,0,-length/2);
  _axis[4] = TVector3(0,1,0);
  _length[4] = height;
  _width[4] = (width_top+width_bot)/2;
  _dwdl[4] = (width_top-width_bot)/height;

  //  Back end

  _normal[5] = TVector3(0,0,1);
  _point[5] = TVector3(0,0,length/2);
  _axis[5] = TVector3(0,1,0);
  _length[5] = height;
  _width[5] = (width_top+width_bot)/2;
  _dwdl[5] = (width_top-width_bot)/height;

  _box = NULL;
  _hit = NULL;
}

void TofBar::transform(const TVector3 &d,const TRotation &r) {
  for ( int i=0; i<6; i++ ) {
    _normal[i] = r*_normal[i];
    _point[i] = d+r*_point[i];
    _axis[i] = r*_axis[i];
  }
}

void TofBar::DrawBar(int color) {
  if ( _box == NULL ) {
    float x[5], y[5];
    TVector3 vec;
    TVector3 edge;
    vec = _normal[0].Cross(_axis[0]);
    edge = _point[0] + _width[0]/2*vec;
    x[0] = edge.X();
    y[0] = edge.Y();
    edge = _point[0] - _width[0]/2*vec;
    x[1] = edge.X();
    y[1] = edge.Y();
    vec = _normal[1].Cross(_axis[1]);
    edge = _point[1] - _width[1]/2*vec;
    x[3] = edge.X();
    y[3] = edge.Y();
    edge = _point[1] + _width[1]/2*vec;
    x[2] = edge.X();
    y[2] = edge.Y();
    x[4] = x[0];
    y[4] = y[0];
    _box = new TPolyLine(5,x,y);
  }
  _box->SetLineColor(color!=0?color:11);
  _box->Draw();
}

void TofBar::DrawHit(TofHitBar *hit) {
  if ( _hit == NULL ) {
    _hit = new TMarker(0,0,20);
    _hit->SetMarkerColor(4);
    _hit->SetMarkerSize(0.5);
  }
  _hit->DrawMarker(hit->inPoint().X(),hit->inPoint().Y());
  _hit->DrawMarker(hit->outPoint().X(),hit->outPoint().Y());
}

void TofDetector::DrawHit(TofHitBar *hit) {
  int i = hit->BarNumber();
  bars[i]->DrawHit(hit);
}

TofDetector::TofDetector(int ilev) {
  const double TofSurvey[72][5] = {
#include "TofSurvey.icc"
  };
  const double UsyninAngles[72] = {
#include "UsyninAngles.icc"
  };

//
// Create the surface descriptions for all bars
//
/*
  HepTransform3D disp[3];
  disp[0] = HepTranslate3D(-4.06908,0.05842,0.0)*
            HepRotateZ3D(-1.66*M_PI/180.0);
  disp[1] = HepTransform3D(HepRotation(),Hep3Vector(0,0,0));
  disp[2] = HepTranslate3D(4.06908,0.05842,0.0)*
            HepRotateZ3D(1.66*M_PI/180.0);
*/
  TVector3 disp[3];
  TRotation rot[3];
  disp[0] = TVector3(-4.06908,0.05842,0.0);
  disp[2] = TVector3(4.06908,0.05842,0.0);
  rot[0].RotateZ(-1.66*M_PI/180.0);
  rot[2].RotateZ(1.66*M_PI/180.0);

  for ( int i=0; i<216; i++ ) {
    bars[i] = new TofBar();
    bars[i]->SetBarNumber(i);
    bars[i]->SetPhi(i*2*M_PI/216.0);
    bars[i]->transform(disp[i%3],rot[i%3]);
  }
  _inrad = -1;
  _outrad = -1;
  for ( int i=0; i<72; i++ ) {
/*
    HepTranslate3D displacement(TofSurvey[i][0],
                                TofSurvey[i][1],
                                TofSurvey[i][2]);
    HepRotateZ3D nominalrotation(atan2(TofSurvey[i][1],
                                       TofSurvey[i][0]));
    HepTransform3D misalignment = HepRotateX3D(TofSurvey[i][4])*
                                  HepRotateZ3D(TofSurvey[i][3]);
    HepTransform3D tiny_rotation = (ilev==SURVEY)?HepRotateZ3D():
                                   HepRotateZ3D(UsyninAngles[i]/1e3);
    HepTransform3D notorsion = tiny_rotation*
                               displacement*
                               nominalrotation*
                               HepRotateZ3D(TofSurvey[i][3])*
                               HepRotateZ3D(M_PI/2.);
*/
    TVector3 displacement(TofSurvey[i][0],
                          TofSurvey[i][1],
                          TofSurvey[i][2]);
    TRotation nominalrotation;
    nominalrotation.RotateZ(atan2(TofSurvey[i][1],TofSurvey[i][0]));
    TRotation misalignment;
    misalignment.RotateZ(TofSurvey[i][3]);
    misalignment.RotateX(TofSurvey[i][4]);
    TRotation tiny_rotation;
    if ( ilev == SURVEY ) tiny_rotation.RotateZ(UsyninAngles[i]/1e3);
    TRotation tippy;
    tippy.RotateZ(M_PI/2);
    tippy.RotateZ(TofSurvey[i][3]);
    TRotation notorsion = tiny_rotation*nominalrotation*tippy;
    TVector3 wherever = tiny_rotation*displacement;

    for ( int j=3*i; j<3*i+3; j++ ) {
      bars[j]->transform(wherever,notorsion);
      double r_in, r_out;
      bars[j]->GetEdges( r_in, r_out );
      if ( _inrad  < 0 || _inrad  > r_in  ) _inrad  = r_in;
      if ( _outrad < 0 || _outrad < r_out ) _outrad = r_out;
    }
  }
  for ( int i=1; i<216; i++ ) {
    _phibnd[i] = 0.5*(bars[i]->getMinPhi()+bars[i-1]->getMaxPhi());
  }
  double p;
  p = bars[215]->getMaxPhi();
  if ( p > M_PI ) p -= 2*M_PI;
  _phibnd[0] = 0.5*(p+bars[0]->getMinPhi());
  p = bars[0]->getMinPhi();
  if ( p < 0 ) p += 2*M_PI;
  _phibnd[216] = 0.5*(p+bars[215]->getMaxPhi());
}

TofDetector::~TofDetector() {
  for ( int i=0; i<216; i++ ) {
    delete bars[i];
  }
}

void TofBar::GetEdges(double &rmin,double &rmax) {
  rmin = -1;
  rmax = -1;
  TVector3 vec;
  TVector3 edge;

  double phi = _point[0].Phi();
  if ( phi < 0 ) phi += 2*M_PI;
  
  for ( int i=0; i<2; i++ ) {
    double r_edge, p0, p1, pmin_edge, pmax_edge;
    vec = _normal[i].Cross(_axis[i]);
    edge = _point[i] + _width[i]/2*vec;
    r_edge = edge.Perp();
    if ( rmin < 0 || r_edge < rmin ) rmin = r_edge;
    if ( rmax < 0 || r_edge > rmax ) rmax = r_edge;
    p0 = edge.Phi();
    if ( phi > 0.5*M_PI && p0 < 0 ) p0 += 2*M_PI;
    edge = _point[i] - _width[i]/2*vec;
    r_edge = edge.Perp();
    if ( rmin < 0 || r_edge < rmin ) rmin = r_edge;
    if ( rmax < 0 || r_edge > rmax ) rmax = r_edge;
    p1 = edge.Phi();
    if ( phi > 0.5*M_PI && p1 < 0 ) p1 += 2*M_PI;
    if ( p0 < p1 ) {
      pmin_edge = p0;
      pmax_edge = p1;
    }
    else {
      pmax_edge = p0;
      pmin_edge = p1;
    }
    if ( i == 0 || pmin_edge < _phimin ) _phimin = pmin_edge;
    if ( i == 0 || pmax_edge > _phimax ) _phimax = pmax_edge;
  }
}

int TofDetector::barNumberAtPhi(double phi) const {
  double p = phi;
  while ( p < 0 ) p += 2*M_PI;
  while ( p > 2*M_PI ) p -= 2*M_PI;
  int ibar = (int)(p*216/(2*M_PI));
  do {
    if ( _phibnd[ibar] > p ) {
      ibar -= 1;
    }
    else if ( _phibnd[ibar+1] <= p ) {
      ibar += 1;
    }
  }
  while ( _phibnd[ibar] > p || _phibnd[ibar+1] <= p );
  return ibar<216?ibar:0;
}

void TofDetector::Draw() {
  for ( int i=0; i<216; i++ ) {
    bars[i]->DrawBar();
  }
}

void TofDetector::DrawBar(int i,int color) {
  bars[i]->DrawBar(color);
}

TofHitBar *TofDetector::extrapolateTrack(int ibar,const float *par,
                                                  const float *fullcov) {

  double kappa = par[0];
  double phi0 = par[1];
  double d0 = par[2];
  double tanlambda = par[3];
  double z0 = par[4];
  double *cov = NULL, covrz[3];
  if ( fullcov != NULL ) {
    covrz[0] = fullcov[14];
    covrz[1] = fullcov[13];
    covrz[2] = fullcov[9];
    cov = covrz;
  }

  double alpha[2], beta[2], gamma[2], xi[2];
  const double expand = 0.01;

  const TVector3 *normal = bars[ibar]->getNormal();
  const TVector3  *point = bars[ibar]->getPoint();
  const TVector3 *axis = bars[ibar]->getAxis();

  const double *length = bars[ibar]->getLength();
  const double *width = bars[ibar]->getWidth();
  const double *dwdl = bars[ibar]->getDwdl();

  alpha[0] = (1+2*d0*kappa)*sin(phi0);
  beta[0] = -(1+2*d0*kappa)*cos(phi0);
  gamma[0] = d0*(1+d0*kappa);
  xi[0] = 1+2*kappa*d0;

  TofHitBar *hit = new TofHitBar(ibar);

  TVector3 in, out;
  TVector3 dir_in, dir_out;
  int went_in  = -1;
  int went_out = -1;
  double ss_in = 0, ss_out = 0;
  double u_in = 0, v_in = 0, u_out = 0, v_out = 0;
  
  for ( int i=0; ( went_in < 0 || went_out < 0 ) && i<6; i++ ) {
    
//  First extrapolate to planes with normal in the z-direction.

    if ( fabs(normal[i].Z()) > 1.e-6 &&
         _outrad*2*tanlambda > fabs(point[i].Z()-z0) ) {
      double ss = (point[i].Z()-z0)/tanlambda;
      if ( ss > 0 ) {
        double xint, yint;
        if ( fabs(kappa) < 2.e-6 ) {
//
//  Approximation for tracks with pT > 1 TeV/c...  oh yeah.
//
          xint = -d0*sin(phi0) + ss*cos(phi0);
          yint = d0*cos(phi0) + ss*sin(phi0);
          cout << "xint = " << xint << ", yint = " << yint << endl;
        }
        else {
          xint = (-alpha[0]+sin(phi0+2*kappa*ss))/(2*kappa);
          yint = (-beta[0]-cos(phi0+2*kappa*ss))/(2*kappa);
        }
        double zint = point[i].Z();
        TVector3 hit(xint,yint,zint);
        TVector3 vec = hit - point[i];
//
//  There was reportedly some unreliability near here...
//
        vec.Dump();
        axis[i].Dump();
	double perp2 = vec.Perp2(axis[i]);
        if ( fabs(perp2) < 1.e-11 ) perp2 = 0;
        cout << " perp2 = " << perp2 << endl;
        assert(perp2>=0);
        double perp = sqrt(perp2);
        double parallel = vec.Dot(axis[i]);
        if ( fabs(parallel) < length[i]/2 ) {
          if ( fabs(perp) < (width[i]+dwdl[i]*parallel)/2 ) {
            TVector3 dir(-(2*kappa*yint+beta[0]),
                            2*kappa*xint+alpha[0],
                            tanlambda);
            if ( dir.Dot(normal[i]) < 0 ) {
              went_in = i;
              ss_in = ss;
              in = hit;
              dir_in = dir;
              u_in = parallel;
              v_in = normal[i].Dot(axis[i].Cross(vec));
            }
            else {
              went_out = i;
              ss_out = ss;
              out = hit;
              dir_out = dir;
              u_out = parallel;
              v_out = normal[i].Dot(axis[i].Cross(vec));
            }
          }
        }
      }
    }
    else if ( fabs(kappa) > 1.e-6 ) {

//  Extrapolate to planes with normal perpendicular to z axis.

      double xint[2], yint[2], zint[2];
      double a[3], c[3];
   
      alpha[1] = -normal[i].X();
      beta[1] = -normal[i].Y();
      gamma[1] = normal[i].X()*point[i].X() + normal[i].Y()*point[i].Y();
      a[0] = alpha[0]*alpha[0]+beta[0]*beta[0];
      a[1] = 1;
      a[2] = alpha[0]*alpha[1]+beta[0]*beta[1];
      double det = alpha[0]*beta[1] - alpha[1]*beta[0];
      c[0] = kappa*kappa;
      c[1] = a[2]*a[2] - 0.5*(a[0]+a[1]) - 2*a[2]*gamma[1]*kappa;
      c[2] = a[1]*gamma[0]*gamma[0] - 2*a[2]*gamma[0]*gamma[1] +
             a[0]*gamma[1]*gamma[1];
      double root = c[1]*c[1] - 4*c[0]*c[2];
      if ( root > 0 ) {
        if ( c[0] > 0 ) {
          double cinv = 1/c[0];
          if ( c[1]*cinv < 0 && c[2]*cinv > 0 && fabs(det) > 1.e-10 ) {
            root = sqrt(root);
            double rad2[2];
            rad2[0] = 0.5*cinv*(-c[1]-root);
            rad2[1] = 0.5*cinv*(-c[1]+root);
            double idet = 1/det;
            double ab1 = -idet*(beta[1]*gamma[0] - beta[0]*gamma[1]);
            double ab2 = idet*(alpha[1]*gamma[0] - alpha[0]*gamma[1]);
            double ac1 = -idet*beta[1]*kappa;
            double ac2 = idet*alpha[1]*kappa;
            TVector3 dir[2];
            double ss[2];
            for ( int j=0; j<2; j++ ) {
              xint[j] = ab1 + ac1*rad2[j];
              yint[j] = ab2 + ac2*rad2[j];
              dir[j] = TVector3(-(2*kappa*yint[j]+beta[0]),
                                   2*kappa*xint[j]+alpha[0],
                                   tanlambda);
//
//  Calculate arc lengths...
//
              double dh1 = kappa*(xint[j]*xint[j]+yint[j]*yint[j])+
                           alpha[0]*xint[j]+beta[0]*yint[j]+gamma[0];
              double dh = dh1*(1-kappa*dh1);
              double kdh = 4*kappa*dh1;
              if ( fabs(dh*kappa) > expand && kdh > -1 ) {
                dh = 0.5/kappa*(sqrt(1+kdh)-1);
              }
              double xon = (xint[j]-dh*alpha[0])/(1+2*kappa*dh);
              double yon = (yint[j]-dh*beta[0])/(1+2*kappa*dh);
              double cross = alpha[0]*yon - beta[0]*xon;
              double a1 = fabs(2*kappa)*cross/xi[0];
              double a2 = 2*kappa*(alpha[0]*xon+beta[0]*yon)/xi[0] + xi[0];
              double ang;
              if ( a2 != 0 ) {
                ang = atan2(a1,a2);
              }
              else {
                ang = (a1>0)?M_PI/2:-M_PI/2;
              }
              if ( ang < 0 ) ang += 2*M_PI;
              ss[j] = ang/fabs(2*kappa);
              zint[j] = z0 + ss[j]*tanlambda;
              TVector3 hit(xint[j],yint[j],zint[j]);
              TVector3 vec = hit - point[i];
//
//  Agagin, some unreliability reported near here...
//
              double perp2 = vec.Perp2(axis[i]);
              if ( fabs(perp2) < 1.e-11 ) perp2 = 0;
              assert(perp2>=0);
              double perp = sqrt(perp2);
              double parallel = vec.Dot(axis[i]);
              if ( fabs(parallel) < length[i]/2 ) {
                if ( fabs(perp) < (width[i]+dwdl[i]*parallel)/2 ) {
                  if ( dir[j].Dot(normal[i]) < 0 ) {
                    went_in = i;
                    ss_in = ss[j];
                    in = hit;
                    dir_in = dir[j];
                    u_in = parallel;
                    v_in = normal[i].Dot(axis[i].Cross(vec));
                  }
                  else {
                    went_out = i;
                    ss_out = ss[j];
                    out = hit;
                    dir_out = dir[j];
                    u_out = parallel;
                    v_out = normal[i].Dot(axis[i].Cross(vec));
                  }
                }
              }
            }
          }
        }
      }
    }
    else {
//
//  Handle the case where the curvature vanishes...
//
      TVector3 x0(-d0*sin(phi0),d0*cos(phi0),z0);
      TVector3 dir(cos(phi0),sin(phi0),tanlambda);
      if ( fabs(normal[i].Dot(dir)) > 1.e-6 ) {
        double ss = (normal[i].Dot(point[i]-x0))/normal[i].Dot(dir);
        if ( ss > 0 ) {
          double xint = -d0*sin(phi0) + ss*cos(phi0);
          double yint = d0*cos(phi0) + ss*sin(phi0);
          double zint = z0 + ss*tanlambda;
          TVector3 hit(xint,yint,zint);
          TVector3 vec = hit - point[i];
//
//  There was reportedly some unreliability near here...
//
          double perp2 = vec.Perp2(axis[i]);
          if ( fabs(perp2) < 1.e-11 ) perp2 = 0;
          assert(perp2>=0);
          double perp = sqrt(perp2);
          double parallel = vec.Dot(axis[i]);
          if ( fabs(parallel) < length[i]/2 ) {
            if ( fabs(perp) < (width[i]+dwdl[i]*parallel)/2 ) {
              if ( dir.Dot(normal[i]) < 0 ) {
                went_in = i;
                ss_in = ss;
                in = hit;
                dir_in = dir;
                u_in = parallel;
                v_in = normal[i].Dot(axis[i].Cross(vec));
              }
              else {
                went_out = i;
                ss_out = ss;
                out = hit;
                dir_out = dir;
                u_out = parallel;
                v_out = normal[i].Dot(axis[i].Cross(vec));
              }
            }
          }
        }
      }
    }
  }
  if ( went_in >= 0 && went_out >= 0 ) {
    hit->setInOut(in,out);
    hit->setArcLengthInOut(ss_in,ss_out);
    hit->setDirInOut(dir_in,dir_out);
    hit->setSurfaceInOut(went_in,went_out);
    hit->setUVInOut(u_in,v_in,u_out,v_out);

    if ( cov != NULL ) {
//
//  Test that covariance matrix is positive definite.
//
      double det = cov[0]*cov[2] - cov[1]*cov[1];
      assert(det>0);
      double uerr_in = sqrt(cov[0]+2*ss_in*cov[1]+ss_in*ss_in*cov[2]);
      double uerr_out = sqrt(cov[0]+2*ss_out*cov[1]+ss_out*ss_out*cov[2]);
      hit->setUErrorInOut(uerr_in,uerr_out);
    }
    return hit;
  }

  delete hit;

  if ( went_in == -1 &&
       went_out == -1 ) return NULL;

  if ( went_in >= 0 && went_out < 0 ) {
    cerr << "Track went into side " << went_in << " of bar "
         << ibar << " but did not come out." << endl;
    return NULL;
  }
  else {
    cerr << "Track came out of side " << went_out << " of bar "
         << ibar << " but did not go in." << endl;
    return NULL;
  }
}

TofHitBar *TofDetector::extrapolateTrack(int ibar,
                                         TVector3 &p,TVector3 &v) {
  double kappa = 0;
  TVector3 t = v*(1/v.Perp());
  double d0 = -p.X()*t.Y() + p.Y()*t.X();
  double phi0 = atan2(t.Y(),t.X());
  if ( phi0 < 0 ) phi0 += 2*M_PI;
  double tanlambda = t.Z();
  double s0 = p.X()*t.X() + p.Y()*t.Y();
  double z0 = p.Z() - t.Z()*s0;
  float par[5] = { kappa, phi0, d0, tanlambda, z0 };

  return extrapolateTrack(ibar,par);
}

int TofDetector::linearizeTrack(double r,const float *par,
                                TVector3 &pnt,TVector3 &dir) {
  double k = par[0];
  double p = par[1];
  double d = par[2];

  double arg = (2*k*k*d*d+2*k*d+1-2*r*r*k*k)/(2*k*d+1);
  if ( fabs(arg) > 1 ) return 0;
  double s = fabs(acos(arg)/(2*k));
  double x = -(d+1/(2*k))*sin(p) + sin(p+2*k*s)/(2*k);
  double y = (d+1/(2*k))*cos(p) - cos(p+2*k*s)/(2*k);
  double z = par[4] + s*par[3];
  TVector3 pp(x,y,z);
  pnt = pp;
  double dx = cos(p+2*k*s);
  double dy = sin(p+2*k*s);
  double dz = par[3];
  TVector3 vv(dx,dy,dz);
  dir = vv;
  return 1;
}

int TofDetector::getHitBarNumbers(const float *par,int &bar_in,int &bar_out) {
  double k = par[0];
  double p = par[1];
  double d = par[2];

  double a = k/(1+2*k*d);
  double b = d - a*d*d;

  double r_in  = _inrad;
  double r_out = _outrad;

  double phi_in, phi_out;
  double arg_in = a*r_in+b/r_in;
  double arg_out = a*r_out+b/r_out;

  if ( fabs(arg_in) > 1 ) return 0;

  phi_in = p + asin(arg_in);
  phi_out = fabs(arg_out)>1?p+M_PI-asin(arg_in):p+asin(arg_out);
  bar_in  = barNumberAtPhi(phi_in);
  bar_out = barNumberAtPhi(phi_out);

  return 1;
}

int TofDetector::getHitBarNumbers(TVector3 &p,TVector3 &v,
                                  int &bar_in,int &bar_out) {
  double kappa = 0;
  TVector3 t = v*(1/v.Perp());
  double d0 = -p.X()*t.Y() + p.Y()*t.X();
  double phi0 = atan2(t.Y(),t.X());
  if ( phi0 < 0 ) phi0 += 2*M_PI;
  float par[3] = { kappa, phi0, d0 };
  int in, out;
  int hit = getHitBarNumbers(par,in,out);
  bar_in = in;
  bar_out = out;
  if ( hit ) {
    cout << "Range of bars for zero-curvature track: "
         << in << "," << out << endl;
  }
  return hit;
}

/* ------------------------------------------------

list<TofHitBar> TofDetector::getHitBars(const float *par,const float *cov) {
  int bar_in, bar_out;
  list<TofHitBar> hitList;
  if ( ! getHitBarNumbers(par,bar_in,bar_out) ) return hitList;
  int di;
  if ( bar_out - bar_in > 108 ) {
    di = -1;
  }
  else if ( bar_in - bar_out > 108 ) {
    di = 1;
  }
  else if ( bar_in > bar_out ) {
    di = -1;
  }
  else {
    di = 1;
  }

  int i = bar_in;
  int done = 0;
  while ( ! done ) {
    TofHitBar *hit = extrapolateTrack(i,par,cov);
    if ( hit != NULL ) {
      hitList.push_back(*hit);
      delete hit;
    }
    if ( i != bar_out ) {
      i += di;
      if ( i < 0 ) i = 215;
      if ( i > 215 ) i = 0;
    }
    else {
      done = 1;
    }
  }
  return hitList;
}

list<TofHitBar> TofDetector::getHitBars(TVector3 &p,TVector3 &v) {
  int bar_in, bar_out;
  list<TofHitBar> hitList;
  if ( ! getHitBarNumbers(p,v,bar_in,bar_out) ) return hitList;
  int di;
  if ( bar_out - bar_in > 108 ) {
    di = -1;
  }
  else if ( bar_in - bar_out > 108 ) {
    di = 1;
  }
  else if ( bar_in > bar_out ) {
    di = -1;
  }
  else {
    di = 1;
  }
  int i = bar_in;
  int done = 0;
  while ( ! done ) {
    TofHitBar *hit = extrapolateTrack(i,p,v);
    if ( hit != NULL ) {
      hitList.push_back(*hit);
      delete hit;
    }
    if ( i != bar_out ) {
      i += di;
      if ( i < 0 ) i = 215;
      if ( i > 215 ) i = 0;
    }
    else {
      done = 1;
    }
  }
  return hitList;
}

------------------------------------- */

void TofDetector::getBarCenter(int i,double &x,double &y) {
  x = bars[i]->point(4).X();
  y = bars[i]->point(4).Y();
}