MucRec2DRoad.cxx

Go to the documentation of this file.

//$id$
//
//$log$
 
/*
 *    2003/12/22   Zhengyun You     Peking University
 * 
 *    2005/02/24   Zhengyun You      Peking University
 *                 transplanted to Gaudi framework
 */
 
using namespace std;
 
#include <iostream>
#include <vector>
#include <CLHEP/Vector/ThreeVector.h>
#include <CLHEP/Geometry/Point3D.h>
 
#include "Identifier/MucID.h"
#include "MucGeomSvc/MucGeomSvc.h"
#include "MucGeomSvc/MucConstant.h"
#include "MucRecEvent/MucRec2DRoad.h"
#include "MucRecEvent/MucRecLineFit.h"
#include "MucRecEvent/MucRecQuadFit.h"
 
const double muckInfinity = 1e30;
// Constructor.
MucRec2DRoad::MucRec2DRoad(const int& part,
               const int& seg,
               const int& orient,
               const float& xVertex,
               const float& yVertex,
               const float& zVertex,
               const int& fittingMethod)
  : m_VertexPos(xVertex, yVertex, zVertex),
    m_VertexSigma(0.0, 0.0, 0.0),
    m_Part(part), m_Seg(seg), m_Orient(orient),
    m_Chi2(0.0), m_DOF(0),
    m_MaxHitsPerGap(0), m_LastGap(0), m_TotalHits(0),
    m_FitOK(false), m_QuadFitOK(false),
    m_HitDistance(5, float(muckInfinity)),
    m_pHits(0), m_fittingMethod(fittingMethod)
{ }
 
// Assignment constructor.
MucRec2DRoad&
MucRec2DRoad::operator=(const MucRec2DRoad& orig)
{
  // Assignment operator.
  if ( this != &orig ) {             // Watch out for self-assignment!
    m_VertexPos     = orig.m_VertexPos;
    m_VertexSigma   = orig.m_VertexSigma;
    m_Index         = orig.m_Index;
    m_Part          = orig.m_Part;
    m_Seg           = orig.m_Seg;
    m_Orient        = orig.m_Orient;
    m_Chi2          = orig.m_Chi2; 
    m_DOF           = orig.m_DOF;
    m_FitOK         = orig.m_FitOK;
    m_MaxHitsPerGap = orig.m_MaxHitsPerGap;
    m_LastGap       = orig.m_LastGap;
    m_TotalHits     = orig.m_TotalHits;
    m_HitDistance   = orig.m_HitDistance;
    m_pHits         = orig.m_pHits;
    m_fittingMethod = orig.m_fittingMethod;
  }
  
  return *this;
}
 
// Copy constructor.
MucRec2DRoad::MucRec2DRoad(const MucRec2DRoad& source)
  : m_VertexPos(source.m_VertexPos),
    m_VertexSigma(source.m_VertexSigma),
    m_Index(source.m_Index),
    m_Part(source.m_Part), m_Seg(source.m_Seg), m_Orient(source.m_Orient),
    m_Chi2(source.m_Chi2), m_DOF(source.m_DOF),
    m_MaxHitsPerGap(source.m_MaxHitsPerGap),
    m_LastGap(source.m_LastGap), m_TotalHits(source.m_TotalHits),
    m_FitOK(source.m_FitOK), 
    m_HitDistance(source.m_HitDistance),
    m_pHits(source.m_pHits),
    m_fittingMethod(source.m_fittingMethod)
{ }
 
// Destructor.
MucRec2DRoad::~MucRec2DRoad()
{ }
 
// Set the index for this road.
void
MucRec2DRoad::SetIndex(const int& index)
{
  if (index >= 0) m_Index = index;
}
 
// Attach the given hit this road.
// Assume that this hit has been verified to be consistent with the road.
void
MucRec2DRoad::AttachHit(MucRecHit* hit)
{
  //  cout << "MucRec2DRoad::AttachHit-I0  hit = " << hit << endl;
  if (!hit) {
    cout << "MucRec2DRoad::AttachHit-E1  null hit pointer!" << endl;
    return ;
  }
  
  int gap = hit->Gap();
  if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
    // The gap number of the hit is out of range.
    cout << "MucRec2DRoad::AttachHit-W3  bad gap number = " << gap
     << endl;
    return;
  }
 
  // Attach the hit to the road.
  //bool hitExist = false;
  for (int i = 0; i < (int)m_pHits.size(); i++) {
    if (m_pHits[i] == hit) return;
  }
  m_pHits.push_back(hit);
  //cout << "before Count " << m_pHits.size() << endl;
   //  m_HitDistance[gap] = dX;
 
  // Now recalculate the total number of hits and the max. number of
  // hits per gap.
  CountHits();
  //cout << "after Count " << m_pHits.size() << endl;
 
  // Redo the "simple" least-squares fit to the positions ...
  float a, sa, b, sb, chi;
  int n;
  int status = SimpleFit(a, b, sa, sb, chi, n);
  if (status < 0) {
    //cout << "MucRec2DRoad::AttachHit-E5  simple fit fail status = " << status << endl;
  }
 
  //cout << "Hit center = " << hit->GetCenterPos() << endl;
  //cout << "After attach hit, slope = " << a << "  intercept = " << b << endl;
}
 
// Attach the given hit this road.
// Assume that this hit has been verified to be consistent with the road.
void
MucRec2DRoad::AttachHitNoFit(MucRecHit* hit)
{
  //  cout << "MucRec2DRoad::AttachHit-I0  hit = " << hit << endl;
  if (!hit) {
    cout << "MucRec2DRoad::AttachHit-E1  null hit pointer!" << endl;
    return ;
  }
  
  int gap = hit->Gap();
  if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
    // The gap number of the hit is out of range.
    cout << "MucRec2DRoad::AttachHit-W3  bad gap number = " << gap
     << endl;
    return;
  }
 
  // Attach the hit to the road.
  //bool hitExist = false;
  for (int i = 0; i < (int)m_pHits.size(); i++) {
    if (m_pHits[i] == hit) return;
  }
  m_pHits.push_back(hit);
  //cout << "before Count " << m_pHits.size() << endl;
   //  m_HitDistance[gap] = dX;
 
  // Now recalculate the total number of hits and the max. number of
  // hits per gap.
  CountHits();
  //cout << "after Count " << m_pHits.size() << endl;
 
}
 
// Max number of consecutive gaps allowed with no hits attached.
// This parameter affects the calculation of the last gap.
void
MucRec2DRoad::SetMaxNSkippedGaps(const int& numGaps)
{
  m_MaxNSkippedGaps = numGaps;
  CountHits();    // recalculate the last gap and the hit counts.
}
 
// Calculate the best-fit straight line with "simple" weights.
int
MucRec2DRoad::SimpleFit(float& slope, 
            float& intercept,
            float& sigmaSlope,
            float& sigmaIntercept,
            float& chisq,
            int&   ndof)
{
//hits in one track can not more than 100.
  if(m_pHits.size()>100){
    cout<<"MucRec2DRoad: too many hits in this track!"<<endl;
    return -1;
    }
  // Assign to temporary arrays to be used in fit.
  float px[100];
  float py[100];
  float pw[100];
  int npoints = 0;
 
  vector<MucRecHit*>::const_iterator iHit;
 
  float weight[100];
 
//    for (int i = 0; i < m_pHits.size(); i++) {
//  cout<<"info: "<<m_pHits[i]->Seg()<<" "<<m_pHits[i]->Gap()<<" "<< m_pHits[i]->Strip()<<endl;
//    }
  for (int i = 0; i < m_pHits.size(); i++) {
     weight[i] = 1;
      
     for(int j = 0; j < m_pHits.size(); j++){
 
    if(j == i) continue;
        if(m_pHits[i]->Part() == m_pHits[j]->Part() &&
       m_pHits[i]->Seg()  == m_pHits[j]->Seg() &&
       m_pHits[i]->Gap()  == m_pHits[j]->Gap() )
     {
       int deltaStrip = fabs(m_pHits[i]->Strip()- m_pHits[j]->Strip());
 
       //cout<<i<<"  "<<m_pHits[i]->Seg()<<"  "<<m_pHits[i]->Gap()<<"  "<<m_pHits[i]->Strip()<<" - "<<m_pHits[j]->Strip()<<endl;
       if(deltaStrip == 0){
             cout<<"deltaStrip == 0 ? check it"<<endl;
       } 
       else{
         weight[i] *=  (deltaStrip+1) * (deltaStrip+1);
       }
 
     }       
     }
  }
  
  for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if (*iHit) {  // Check for a null pointer.
      
      /*
    float localx, localy, localz;
      (*iHit)->GetStrip()->GetCenterPos(localx, localy, localz);
      if ( m_Orient == 0) {
    px[npoints] = localy;
    py[npoints] = localz;
      }
      else {
    px[npoints] = localx;
    py[npoints] = localz;
      }
      npoints++;
      }
      }
      */
      
      
      Hep3Vector center = (*iHit)->GetCenterPos();
      Hep3Vector sigma  = (*iHit)->GetCenterSigma();
      //Hep3Vector sigma(1.0, 1.0, 1.0);
 
      if (m_Part == 1) {
        if ( m_Orient == 0) {
          px[npoints] = center.z();
          py[npoints] = sqrt(center.x()*center.x() 
              + center.y()*center.y());
          if(m_Seg==2) py[npoints] = center.y();  //deal with seg2 seperately! because there is a hole here. 2006.11.9
 
          pw[npoints] = 4.0 * 1.0/(sigma.y()*sigma.y()) / weight[npoints];
        } 
        else {
          px[npoints] = center.x();
          py[npoints] = center.y();
          pw[npoints] = 4.0 * 1.0/(sigma.x()*sigma.x()) / weight[npoints];
        }
      }
      else {
        if ( m_Orient == 0) {
          px[npoints] = center.z();
          py[npoints] = center.y();
          pw[npoints] = 4.0 * 1.0/(sigma.y()*sigma.y()) / weight[npoints];
        }
        else {
          px[npoints] = center.z();
          py[npoints] = center.x();
          pw[npoints] = 4.0 * 1.0/(sigma.x()*sigma.x()) / weight[npoints];
        }
      }  
 
      //cout << " in MucRec2DRoad  ID: " <<(*iHit)->Part()<<" "<<(*iHit)->Seg()<<" "<<(*iHit)->Gap()<<" "<< (*iHit)->Strip()<<" "<<px[npoints] << "  " << py[npoints] << "  " << pw[npoints] << endl;
      npoints++;
 
      if(npoints > 99) 
      {  cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
        return -1;
      }   
  }
}
/*
float px2[100];
float py2[100];
for (int i = 0; i < m_pHits.size(); i++) {
    int hitsInGap = 1;
    px2[i] = -999; py2[i] = -999;
    for(int j = 0; j < m_pHits.size(); j++){
 
        if(j == i) continue;
        if(m_pHits[i]->Part() == m_pHits[j]->Part() &&
                m_pHits[i]->Seg()  == m_pHits[j]->Seg() &&
                m_pHits[i]->Gap()  == m_pHits[j]->Gap() )
        {
            hitsInGap++;
            px2[i] = (px[i]*(hitsInGap-1) + px[j])/hitsInGap;
            py2[i] = (py[i]*(hitsInGap-1) + py[j])/hitsInGap;
            cout<<hitsInGap<<"  "<<px[i]<<"  "<<px[j]<<"  "<<px2[i]<<endl;
        }
    }
}
 
for(int i = 0; i < m_pHits.size(); i++){
    if(px2[i] != -999&&py2[i] != -999){
            px[i] = px2[i]; py[i] = py2[i];
    }
}
*/
 
ndof = npoints - 2;
if (ndof < 0) return -1;
 
if (npoints == 1) {
  if (m_Part == 1) {
    if ( m_Orient == 0) {
      px[npoints] = m_VertexPos.z();
      py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x() 
          + m_VertexPos.y()*m_VertexPos.y());
      if(m_Seg==2) py[npoints] = m_VertexPos.y();
    } 
    else {
      px[npoints] = m_VertexPos.x();
      py[npoints] = m_VertexPos.y();
    }
  }
  else {
    if ( m_Orient == 0) {
      px[npoints] = m_VertexPos.z();
      py[npoints] = m_VertexPos.y();
    }
    else {
      px[npoints] = m_VertexPos.z();
      py[npoints] = m_VertexPos.x();
    }
  }  
  pw[npoints] = 1.0;
  npoints++;
} 
else {
  if (npoints == 0 ) {
    return -1;
  }
}
 
// Do the fits here.
MucRecLineFit fit;
  int status = fit.LineFit(px, py, pw, npoints,
               &slope, &intercept, &chisq,
               &sigmaSlope, &sigmaIntercept);
 
 
  float tempslope, tempintercept,tempchisq, tempsigmaslope, sigmaPos;
  int status4 = fit.LineFit(px, py, pw,m_Part,m_Seg,m_Orient, npoints,
      &tempslope, &tempintercept, &tempchisq,
      &tempsigmaslope, &sigmaPos);
  
  MucRecQuadFit quadfit;
  float quad_a, quad_b, quad_c, sigmaquad_a, sigmaquad_b, sigmaquad_c, chisq_quad;
  int whichhalf, status2;
 
  if(m_fittingMethod == 2){
  status2 = quadfit.QuadFit(px, py, pw, npoints,
                &quad_a, &quad_b, &quad_c, &whichhalf, &chisq_quad,
                &sigmaquad_a, &sigmaquad_b, &sigmaquad_c);
 
  }
  //cout << " in MucRec2DRoad slope " << slope << " "<<intercept<<endl;
  
  if (slope > 1.0e10 || slope < -1.0e10) {
    if (m_Seg > 4) slope *= -1.0;  // to avoid wrong direction
  }
 
  m_SimpleSlope          = slope;
  m_SimpleSlopeSigma     = sigmaSlope;
  m_SimpleIntercept      = intercept;
  m_SimpleInterceptSigma = sigmaIntercept;
  m_SimplePosSigma       = sigmaPos;   //new 20071227
  m_SimpleChi2           = chisq;
  m_SimpleDOF            = ndof;
  m_SimpleFitOK          = (status == 0) && (chisq < 1000.0);
  m_QuadFitOK            = (status2 == 1);
 
  m_SimpleQuad_a         = quad_a;
  m_SimpleQuad_b         = quad_b;
  m_SimpleQuad_c         = quad_c;
  m_SimpleQuad_whichhalf = whichhalf;
  m_SimpleQuad_aSigma    = sigmaquad_a;
  m_SimpleQuad_bSigma    = sigmaquad_b;
  m_SimpleQuad_cSigma    = sigmaquad_c;
 
  return status;
}
 
 
/////////////////////////////////////////////////
// Calculate the best-fit straight line with "simple" weights. add not use current gap!!!
int   
MucRec2DRoad::SimpleFitNoCurrentGap(
      int currentgap,
      float& slope,
      float& intercept,
      float& sigmaSlope,
      float& sigmaIntercept,
      float& chisq,
      int&   ndof)
{
  // Assign to temporary arrays to be used in fit.
  float px[100];
  float py[100];
  float pw[100];
  int npoints = 0;
 
  int notused = 0;
 
  vector<MucRecHit*>::const_iterator iHit;
 
  float pw2[100];
  for(int i = 0; i< 100; i++) {
    pw2[i] = 1;
    //----if( m_pHits[i]->Gap()==currentgap ) pw2[i] = 9999;
  }
 
  for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if (*iHit) {  // Check for a null pointer.
      if( (*iHit)->Gap()==currentgap ) continue;
 
      Hep3Vector center = (*iHit)->GetCenterPos();
      Hep3Vector sigma  = (*iHit)->GetCenterSigma();
      //Hep3Vector sigma(1.0, 1.0, 1.0);
 
      if (m_Part == 1) {
        if ( m_Orient == 0) {
          px[npoints] = center.z();
          py[npoints] = sqrt(center.x()*center.x()
              + center.y()*center.y());
          if(m_Seg==2) py[npoints] = center.y();  //deal with seg2 seperately! because there is a hole here. 2006.11.9
 
          pw[npoints] = 1.0/(sigma.y()*sigma.y())/pw2[npoints]/pw2[npoints];
        }
        else {
          px[npoints] = center.x();
          py[npoints] = center.y();
          pw[npoints] = 1.0/(sigma.x()*sigma.x())/pw2[npoints]/pw2[npoints];
        }
      }
      else {
        if ( m_Orient == 0) {
          px[npoints] = center.z();
          py[npoints] = center.y();
          pw[npoints] = 1.0/(sigma.y()*sigma.y())/pw2[npoints]/pw2[npoints];
        }
        else {
          px[npoints] = center.z();
          py[npoints] = center.x();
          pw[npoints] = 1.0/(sigma.x()*sigma.x())/pw2[npoints]/pw2[npoints];
        }
      }
 
      //cout << " in MucRec2DRoad  ID: " <<(*iHit)->Part()<<" "<<(*iHit)->Seg()<<" "<<(*iHit)->Gap()<<" "<< (*iHit)->Strip()<<" "<<px[npoints] << "  " << py[npoints] << "  " << pw[npoints] << endl;
      //cout<<"process ngap: "<<currentgap<<" current: "<<(*iHit)->Gap()<<"  x: "<<px[npoints]<<"  y: "<<py[npoints]<<"  weight: "<<pw[npoints]<<endl;
      npoints++;
 
      if(npoints > 99)
      {  cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
        return -1;
      }
    }
  }
 
  ndof = npoints - 2;
  if (ndof < 0) return -1;
 
  if (npoints == 1) {
    if (m_Part == 1) {
      if ( m_Orient == 0) {
        px[npoints] = m_VertexPos.z();
        py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x()
            + m_VertexPos.y()*m_VertexPos.y());
        if(m_Seg==2) py[npoints] = m_VertexPos.y();
      }
      else {
        px[npoints] = m_VertexPos.x();
        py[npoints] = m_VertexPos.y();
      }
    }
    else {
      if ( m_Orient == 0) {
        px[npoints] = m_VertexPos.z();
        py[npoints] = m_VertexPos.y();
      }
      else {
        px[npoints] = m_VertexPos.z();
        py[npoints] = m_VertexPos.x();
      }
    }
    pw[npoints] = 1.0;
    npoints++;
  }
else {
  if (npoints == 0 ) {
    return -1;
  }
 
}
 
// Do the fits here.
MucRecLineFit fit;
int status = fit.LineFit(px, py, pw, npoints,
    &slope, &intercept, &chisq,
    &sigmaSlope, &sigmaIntercept);
 
MucRecQuadFit quadfit;
float quad_a, quad_b, quad_c, sigmaquad_a, sigmaquad_b, sigmaquad_c, chisq_quad;
int whichhalf, status2;
 
if(m_fittingMethod == 2){
  status2 = quadfit.QuadFit(px, py, pw, npoints,
      &quad_a, &quad_b, &quad_c, &whichhalf, &chisq_quad,
      &sigmaquad_a, &sigmaquad_b, &sigmaquad_c);
 
}
//cout << " in MucRec2DRoad slope " << slope << " "<<intercept<<endl;
 
if (slope > 1.0e10 || slope < -1.0e10) {
  if (m_Seg > 4) slope *= -1.0;  // to avoid wrong direction
}
 
////////2009-03-12  
  m_SimpleSlope          = slope;
  m_SimpleSlopeSigma     = sigmaSlope;
  m_SimpleIntercept      = intercept;
  m_SimpleInterceptSigma = sigmaIntercept;
  //m_SimplePosSigma       = sigmaPos;   //new 20071227
  m_SimpleChi2           = chisq;
  m_SimpleDOF            = ndof;
  m_SimpleFitOK          = (status == 0) && (chisq < 1000.0);
  m_QuadFitOK            = (status2 == 1);
 
  m_SimpleQuad_a         = quad_a;
  m_SimpleQuad_b         = quad_b;
  m_SimpleQuad_c         = quad_c;
  m_SimpleQuad_whichhalf = whichhalf;
  m_SimpleQuad_aSigma    = sigmaquad_a;
  m_SimpleQuad_bSigma    = sigmaquad_b;
  m_SimpleQuad_cSigma    = sigmaquad_c;
 
 
return status;
}
 
 
 
 
// Fit (with weights) the hit centers to a straight line.
// This is a helper function for the Project methods.
// Give an estimated position (x,y,z) of the point to which we will
// project.
  int
MucRec2DRoad::Fit(const float& x,
    const float& y,
    const float& z,
    float& slope,      float& intercept,
    float& sigmaSlope, float& sigmaIntercept,
    float& chisq,      int&   ndof)
{
  int status;
 
  // If the "simple" fit has not been done yet, do it now.
  if (!m_SimpleFitOK) {
    // Least-squares fit to the positions ...
    float a, sa, b, sb, chi;
    int n;
    status = SimpleFit(a, b, sa, sb, chi, n);
    if (status < 0) {
      //cout << "MucRec2DRoad::Fit-E2  simple fit fail status = "
      //  << status << endl;
      return status;
    }
  }
 
  // Assign to temporary arrays to be used in fit.
  float px[100];
  float py[100];
  float pw[100];
  int npoints = 0;
  float dx, dy, dr;
 
 
  // include vertex point when do the fancy fitting
  //if (m_Orient == kHORIZ) {
  //  px[npoints] = m_VertexPos.y();
  //  dx = px[npoints] - y;
  //} else {
  //  px[npoints] = m_VertexPos.x();
  //  dx = px[npoints] - x;
  //}
  //pz[npoints] = m_VertexPos.z();
  //dz = pz[npoints] - z;
  //dr = sqrt(dx*dx + dz*dz);
  //pw[npoints] = WeightFunc(m_SimpleChi2,dr);  
  //npoints++;
 
  vector<MucRecHit*>::const_iterator iHit;
 
  // Determine the weights based on the chi-square of the fit
  // (the scatter of the points should be roughly related to the
  // momentum) and the distance from the center to the estimated
  // location.
 
  //cout << m_pHits.size() << endl;
  for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if (*iHit) {  // Check for a null pointer.
 
      Hep3Vector center = (*iHit)->GetCenterPos();
      Hep3Vector sigma  = (*iHit)->GetCenterSigma();
 
      if (m_Part == 1) {    // change from 0 to 1 at 2006.11.30
        if ( m_Orient == 0) {
          px[npoints] = center.z();
          dx = px[npoints] - z;
          py[npoints] = sqrt(center.x()*center.x() 
              + center.y()*center.y());
          if(m_Seg==2) py[npoints] = center.y();  //deal with seg2 seperately! because there is a hole here. 2006.11.9
          dy = py[npoints] - sqrt(x*x + y*y);
        } 
        else {
          px[npoints] = center.x();
          dx = px[npoints] - x;
          py[npoints] = center.y();
          dy = py[npoints] - y;
        }
      }
      else {
        if ( m_Orient == 0) {
          px[npoints] = center.z();
          dx = px[npoints] - z;
          py[npoints] = center.y();
          dy = py[npoints] - y;
        }
        else {
          px[npoints] = center.z();
          dx = px[npoints] - z;
          py[npoints] = center.x();
          dy = py[npoints] - x;
        }
      }  
 
      dr = sqrt(dx*dx + dy*dy);
      pw[npoints] = WeightFunc(m_SimpleChi2, dr);
      //pw[npoints] = WeightFunc(m_SimpleChi2,dr);
 
      //    cout << "       " << npoints  << "     "
      // << px[npoints] << "  " << py[npoints] << "  " << pw[npoints]
      // << "         " << dx << "  " << dy << "  " << dr
      // << endl;
 
      npoints++;
 
      if(npoints > 99)
      {  cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
        return -1;
      }   
 
    }
  }
 
  // Refit ...
  ndof = npoints - 2;
  if (npoints == 1) {
    if (m_Part == 1) {     // change from 0 to 1 at 2006.11.30
      if ( m_Orient == 0) {
        px[npoints] = m_VertexPos.z();
        py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x() 
            + m_VertexPos.y()*m_VertexPos.y());
        if(m_Seg==2) py[npoints] = m_VertexPos.y();
      } 
      else {
        px[npoints] = m_VertexPos.x();
        py[npoints] = m_VertexPos.y();
      }
    }
    else {
      if ( m_Orient == 0) {
        px[npoints] = m_VertexPos.z();
        py[npoints] = m_VertexPos.y();
      }
      else {
        px[npoints] = m_VertexPos.z();
        py[npoints] = m_VertexPos.x();
      }
    }  
    pw[npoints] = 1.0;
    npoints++;
  } 
  else {
    if (npoints == 0) {
      return -1;
    }
  }
 
  MucRecLineFit fit;
  //cout << "npoints " << npoints << endl;
  status = fit.LineFit(px, py, pw, npoints,
      &slope, &intercept, &chisq,
      &sigmaSlope, &sigmaIntercept);
  m_DOF  = ndof;
  m_Chi2 = chisq;
 
  if (status < 0) {
    //cout << "MucRec2DRoad::Fit-E3  fit fail status = " << status << endl;
  }
 
  return status;
}
 
// Where does the trajectory of this road intersect a specific gap?
void
MucRec2DRoad::Project(const int& gap, 
             float& x, float& y, float& z, float& x2, float& y2, float& z2)
{
  float sigmaX, sigmaY, sigmaZ;
 
  x = 0.0; sigmaX = 0.0; x2 = 0.0;
  y = 0.0; sigmaY = 0.0; y2 = 0.0;
  z = 0.0; sigmaZ = 0.0; z2 = 0.0;
 
  if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
    cout << "MucRec2DRoad::Project-E1  invalid gap = " << gap << endl;
    return;
  }
  
  // Determine the projection point of the "simple" fit to the desired gap.
  float x0, y0, z0, sigmaX0, sigmaY0, sigmaZ0;
  float phi = m_Seg*0.25*kPi;
  MucGeoGeneral *geom = MucGeoGeneral::Instance();
 
  if (m_Part == 1) {
    if (m_Orient == 0) {
      geom->FindIntersection(m_Part, m_Seg, gap,
                 m_SimpleSlope*cos(phi), m_SimpleSlope*sin(phi), 1.0,
                 m_SimpleIntercept*cos(phi),m_SimpleIntercept*sin(phi), 0.0,
                 m_SimpleSlopeSigma,     0.0, 0.0,
                 m_SimpleInterceptSigma, 0.0, 0.0, 
                 x0, y0, z0,
           sigmaX0, sigmaY0, sigmaZ0);
 
      if(m_SimpleSlope>10000){  //the line is in right center of this segment, we can not get intersection in y coordinate, in this case, m_SimpleIntercept is in z coordinate.                                               
        geom->FindIntersection(m_Part, m_Seg, gap,
        m_SimpleSlope*cos(phi), m_SimpleSlope*sin(phi), 1.0,
        0.0, 0.0, m_SimpleIntercept,
        m_SimpleSlopeSigma,     0.0, 0.0,
        m_SimpleInterceptSigma, 0.0, 0.0,
        x0, y0, z0,
        sigmaX0, sigmaY0, sigmaZ0);
 
      }
 
    }
    else {
      geom->FindIntersection(m_Part, m_Seg, gap,
                 1.0, m_SimpleSlope,          0.0,
                 0.0, m_SimpleIntercept,      0.0,
                 0.0, m_SimpleSlopeSigma,     0.0,
                 0.0, m_SimpleInterceptSigma, 0.0, 
                 x0, y0, z0,
                 sigmaX0, sigmaY0, sigmaZ0);
      //cout<<"in MucRec2DRoad line Project xyz0 = "<<x0<<" "<<y0<<" "<<z0<<endl;
 
    }
  }
  else {
    if (m_Orient == 0) {
      geom->FindIntersection(m_Part, m_Seg, gap,
                 m_SimpleSlope, m_SimpleSlope, 1.0,
                 0.0, m_SimpleIntercept,      0.0,
                 0.0, m_SimpleSlopeSigma,     0.0,
                 0.0, m_SimpleInterceptSigma, 0.0, 
                 x0, y0, z0,
                 sigmaX0, sigmaY0, sigmaZ0);
    }
    else {
      geom->FindIntersection(m_Part, m_Seg, gap,
                 m_SimpleSlope, m_SimpleSlope, 1.0,
                 m_SimpleIntercept,      0.0, 0.0,
                 m_SimpleSlopeSigma,     0.0, 0.0,
                 m_SimpleInterceptSigma, 0.0, 0.0,
                 x0, y0, z0,
                 sigmaX0, sigmaY0, sigmaZ0);
    }
    //cout << " part " << m_Part 
    // << " seg "  << m_Seg 
    // << " gap "  << gap
    // << " orient "  << m_Orient
    // << " slope = " << m_SimpleSlope
    // << endl;
  }
    
  //cout << "In find intersection x0 = " << x0 << " y0 = " << y0 << " z0 = " << z0 << endl; 
 
  float a,b,sa,sb,chi2;
  int ndof;
 
  int status = Fit(x0, y0, z0, a, b, sa, sb, chi2, ndof);
 
//   m_FitOK = (status == 0) && (chi2<1000.0);
//   if (!fFitOK) {
//     cout << "MucRec2DRoad::Project-E2  fit fail status = "
//   << status << "  npoints = " << npoints << "  chi-sq = "
//   << chi2 << endl;
//   }
 
//   // Assign to fields of TMuiRoad object.
//   m_DOF  = npoints - 2;
//   m_Chi2 = chi2;
 
  if (m_Part == 1) {    //change from 0 to 1 at 2006.11.30
    if (m_Orient == 0) {
      geom->FindIntersection(m_Part, m_Seg, gap,
                 a*cos(phi), a*sin(phi), 1.0,
                 //a,  0.0, 1.0,
                 b*cos(phi), b*sin(phi), 0.0,
                 sa, 0.0, 0.0,
                 sb, 0.0, 0.0, 
                 x,  y,   z,
                 sigmaX, sigmaY, sigmaZ);
 
      if(fabs(a)>10000){ /// 
        geom->FindIntersection(m_Part, m_Seg, gap,
            a*cos(phi), a*sin(phi), 1.0,
            0.0, 0.0, b,
            //a,  0.0, 1.0,
            //b,  0.0, 0.0,
            sa, 0.0, 0.0,
            sb, 0.0, 0.0,
            x,  y,   z,
            sigmaX, sigmaY, sigmaZ);
      
    }
    }
    else {
      geom->FindIntersection(m_Part, m_Seg, gap,
                 1.0, a,  0.0,
                 0.0, b,  0.0,
                 0.0, sa, 0.0,
                 0.0, sb, 0.0, 
                 x,   y,  z,
                 sigmaX, sigmaY, sigmaZ);
 
      if(m_fittingMethod == 2 && m_QuadFitOK ){
    float a, b, c;
    float sa, sb, sc, chi2x; int ydof; int whichhalf;
    
    GetSimpleFitParams(a, b, c, whichhalf, sa, sb, sc, chi2x, ydof);
    geom->FindIntersection(m_Part, m_Seg, gap,
                   10E30, 0.0, m_SimpleIntercept,  0.0,   //vy = infinite
                   a, b, c,  //y = a*x*x + b*x +c
                   whichhalf,
                   0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                   x, y, z, x2, y2, z2,
                   sigmaX, sigmaY, sigmaZ);
 
 
      }
 
    }
  }
  else {
    if (m_Orient == 0) {
      geom->FindIntersection(m_Part, m_Seg, gap,
                 a,   a,  1.0,
                 0.0, b,  0.0,
                 0.0, sa, 0.0,
                 0.0, sb, 0.0, 
                 x,  y,   z,
                 sigmaX, sigmaY, sigmaZ);
    }
    else {
      geom->FindIntersection(m_Part, m_Seg, gap,
                 a,  a, 1.0,
                 b,  0.0, 0.0,
                 sa, 0.0, 0.0,
                 sb, 0.0, 0.0,
                 x,  y,   z,
                 sigmaX, sigmaY, sigmaZ);
    }
  }
 
  return;
}
 
// A unique identifier for this road in the current event.
int
MucRec2DRoad::GetIndex() const
{
  return m_Index;
}
 
// In which part was this road found?
int
MucRec2DRoad::GetPart() const { return m_Part; }
 
// In which seg was this road found?
int
MucRec2DRoad::GetSeg() const { return m_Seg; }
 
// In which orientation was this road found?
int
MucRec2DRoad::GetOrient() const { return m_Orient; }
 
// Position of the vertex point.
void
MucRec2DRoad::GetVertexPos(float& x, float& y, float& z) const
{
  x = m_VertexPos.x();
  y = m_VertexPos.y();
  z = m_VertexPos.z();
 
  return;
}
 
// Which gap is the last one with hits attached to this road?
int
MucRec2DRoad::GetLastGap() const { return m_LastGap; }
 
// Return the number of hits in the gap with the most attached hits.
int
MucRec2DRoad::GetMaxHitsPerGap() const { return m_MaxHitsPerGap; }
 
// Return the total number of hits attached to this road.
int
MucRec2DRoad::GetTotalHits() const { return m_TotalHits; }
 
// How many gaps provide hits attached to this road?
int
MucRec2DRoad::GetNGapsWithHits() const
{
  const int ngap = (int)MucID::getGapMax();
  int gap, count = 0;
  vector<int> firedGap;
  for ( gap = 0; gap < ngap; gap++) {
    firedGap.push_back(0);
  }
 
  vector<MucRecHit*>::const_iterator iHit;
  for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if (*iHit) {  // Check for a null pointer.
      gap = (*iHit)->Gap();
      firedGap[gap] = 1;
    }
  }
  
  for ( gap = 0; gap < ngap; gap++) {
    count += firedGap[gap];
  }
 
  return count;
}
 
// How many hits per gap does this road contain?
int
MucRec2DRoad::GetHitsPerGap(const int& gap) const
{
  if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
    cout << "MucRec2DRoad::GetHitsPerGap-E1  invalid gap = " << gap << endl;
    return 0;
  } 
 
  vector<MucRecHit*>::const_iterator iHit;
  int hitsInGap = 0;
 
  for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    
    if ( !(*iHit) ) {
      cout << "MucRec2DRoad::GetHitsPerGap()-E2 null hit pointer !"  << endl;
      return 0;
    }
    else {
      if( gap == (*iHit)->Gap() ) hitsInGap += 1;
    }
  }
 
  return hitsInGap;
}
 
// Does this road contain any hits in the given seg in a gap?
bool
MucRec2DRoad::HasHitInGap(const int& gap) const
{
  if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
    cout << "MucRec2DRoad::HasHitInGap-E2  invalid gap = " << gap << endl;
    return false;
  }
 
  bool found = false;
  vector<MucRecHit*>::const_iterator iHit;
 
  for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if (*iHit) {  // Check for a null pointer.
      if ( (*iHit)->Gap() == gap ) {
    found = true;
      }
    }
  }
 
  return found;
}
 
// How many hits do two roads share?
int
MucRec2DRoad::GetNSharedHits(const MucRec2DRoad* road2) const
{
  if (!road2) {
    return 0;
  }
 
  int count = 0;
  vector<MucRecHit*>::const_iterator iHit1;
  vector<MucRecHit*>::const_iterator iHit2;
  MucRecHit *hit1, *hit2;
 
  for( iHit1 = m_pHits.begin(); iHit1 != m_pHits.end(); iHit1++){
    for( iHit2 = road2->m_pHits.begin(); 
     iHit2 != road2->m_pHits.end(); iHit2++){
      hit1 = (*iHit1);
      hit2 = (*iHit2);
 
      if ( (hit1 != 0) && (hit2 != 0) ) {
         if (hit1->GetID() == hit2->GetID()) {
       count++;
     }
      }
    }
  }
  
  return count;
}
 
// Slope of trajectory.
float 
MucRec2DRoad::GetSlope() const { return m_SimpleSlope; }
 
// Intercept of trajectory.
float 
MucRec2DRoad::GetIntercept() const { return m_SimpleIntercept; }
 
// Degrees of freedom of trajectory fit.
int
MucRec2DRoad::GetDegreesOfFreedom() const { return m_SimpleDOF; }
 
// Chi-square parameters (per degree of freedom) of trajectory fit.
float
MucRec2DRoad::GetReducedChiSquare() const
{
  if ( (!m_SimpleFitOK) || (m_SimpleDOF < 0) ) {
    return -1.0;
  }
  else if (m_SimpleDOF == 0) {
    return 0.0;
  }
  else {
    return (m_SimpleChi2 / m_SimpleDOF);
  }
}
 
// Get the parameters from the simple fit.
void
MucRec2DRoad::GetSimpleFitParams(float& slope,      float& intercept,
                    float& sigmaSlope, float& sigmaIntercept,
                    float& chisq,      int& ndof) const
{
  slope           = m_SimpleSlope;
  intercept       = m_SimpleIntercept;
  sigmaSlope      = m_SimpleSlopeSigma;
  sigmaIntercept  = m_SimpleInterceptSigma;
  chisq           = m_SimpleChi2;
  ndof            = m_SimpleDOF;
 
  return;
}
 
 
void
MucRec2DRoad::GetPosSigma(float &possigma)const
{
  possigma = m_SimplePosSigma;
 
}
 
 
// Get the parameters from the simple quad fit.
void
MucRec2DRoad::GetSimpleFitParams(float& a,      float& b, float& c, int& whichhalf,
                 float& sigmaa, float& sigmab,  float& sigmac,
                    float& chisq,      int& ndof) const
{
  a           = m_SimpleQuad_a;
  b           = m_SimpleQuad_b;
  c           = m_SimpleQuad_c;
  whichhalf   = m_SimpleQuad_whichhalf;
 
  sigmaa      = m_SimpleQuad_aSigma; 
  sigmab      = m_SimpleQuad_bSigma;
  sigmac      = m_SimpleQuad_cSigma;
 
  chisq           = m_SimpleChi2;
  ndof            = m_SimpleDOF;
 
  return;
}
 
 
 
// Get a pointer to the first found hit attached in a particular gap.
MucRecHit*
MucRec2DRoad::GetHit(const int& gap) const
{
  if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
    cout << "MucRec2DRoad::Hit-E1  invalid gap = " << gap << endl;
    return 0;
  }
  
  vector<MucRecHit*>::const_iterator iHit;
  
  for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if (*iHit) {  // Check for a null pointer.
      if ( (*iHit)->Gap() == gap) { 
    return (*iHit);
      }
    }
  }
  
  return 0L;
}
 
// Distance from a hit to the projection of the road to the gap
// containing the hit.
float
MucRec2DRoad::GetHitDistance(const MucRecHit* hit)
{
  // Use straight-line projection?
  if (!hit) {
    cout << "MucRec2DRoad::GetHitDistance-E1  null hit pointer!" << endl;
    return muckInfinity;
  }
 
  int gap = hit->Gap();
  if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
    cout << "MucRec2DRoad::HitDistance-W2  bad gap number = " << gap << endl;
    return muckInfinity;
  }
 
  if ( hit->GetGap()->Orient() != m_Orient) {
    // The orientation of the hit is different from that of the road.
    cout << "MucRec2DRoad::GetHitDistance-W3 "
     << " hit has wrong orientation = "
     << hit->GetGap()->Orient()
     << endl;
    return muckInfinity;
  }
 
  HepPoint3D r  = hit->GetCenterPos();
  HepPoint3D rl = hit->GetGap()->TransformToGap(r);
  // cout << "hit center " << r << endl;
  // cout << "hit center local " << rl << endl;
 
  // Determine the projection point of the "simple" fit to the desired gap.
  // FIXME(?):  need to use fit with fancy weights instead?
  float delta, delta1,delta2;
  float x0, y0, z0;
  float x2, y2, z2;
  //float sigmaX0, sigmaY0, sigmaZ0;
 
  //cout << "y:x = " << m_SimpleSlope << endl;
 
  Project(gap, x0, y0, z0, x2, y2, z2);
  // cout << " gap = " << gap << " x0 = " << x0
  //      << " y0 = "  << y0  << " z0 = " << z0
  //      << endl;
 
  if(x0==y0&&x0==z0&&x0==0) {x0+=-9999;y0+=-9999;z0+=-9999;}//cout<<"wrong intersection"<<endl;}
  if(x2==y2&&x2==z2&&x2==0) {x2+=-9999;y2+=-9999;z2+=-9999;}//cout<<"wrong intersection"<<endl;}   //wrong intersection!!!
  HepPoint3D s  = HepPoint3D(x0, y0, z0);
  HepPoint3D s_2  = HepPoint3D(x2, y2, z2);
  HepPoint3D sl = hit->GetGap()->TransformToGap(s);
  HepPoint3D s_2l = hit->GetGap()->TransformToGap(s_2);
  //cout << "project center       " << s  << endl;
  
//   cout << "project center local sl= " << sl << endl;
//   cout << "project center local sl= " << s_2l << endl;
//   cout << "project center local rl= " << rl << endl;
  if ( m_Part == 0 ) {
    if ( m_Orient == 0 ) {
      delta1 = fabs( sl.y() - rl.y() );
      delta2= fabs( s_2l.y() - rl.y() );
    } 
    else {
      delta1 = fabs( sl.x() - rl.x() );
      delta2= fabs( s_2l.x() - rl.x() );
    }
  }
  else {
    if ( m_Orient == 0 ) {
      delta1 = fabs( sl.y() - rl.y() );
      delta2= fabs( s_2l.y() - rl.y() );
    } 
    else {
      delta1 = fabs( sl.x() - rl.x() );
      delta2= fabs( s_2l.x() - rl.x() );
    }
  }
//   if(which==0) {
//     if(delta1<delta2)delta = delta1;
//     else delta = delta2;
//   }
  //cout<<"which = "<<which<<"  distance = "<<delta1<<" "<<delta2<<endl;
  
  if(delta1 < delta2)  return delta1;
  else return delta2;
}
 
// Determine the size of the search window in the given gap.
float
MucRec2DRoad::GetSearchWindowSize( const int& gap) const
{
  if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
    cout << "MucRec2DRoad::GetSearchWindowSize-E1  invalid gap = " << gap << endl;
    return 0.0;
  }
 
  // Determine the projection point of the "simple" fit to the last
  // gap and the desired gap.
  // FIXME?  the "delta-x" variable is calculated as the scalar
  // difference between the positions obtained by projecting to the
  // last gap and to the desired gap.
 
  MucGeoGeneral *geom = MucGeoGeneral::Instance();
  float x1, y1, z1, x2, y2, z2, dx, dy, dr;
  float sigmaX, sigmaY, sigmaZ;
 
  if (m_Part == 0) {
    if (m_Orient == 0) {
      geom->FindIntersection(m_Part, 0, m_LastGap,
                 1.0, m_SimpleSlope,          0.0,
                 0.0, m_SimpleIntercept,      0.0,
                 0.0, m_SimpleSlopeSigma,     0.0,
                 0.0, m_SimpleInterceptSigma, 0.0, 
                 x1, y1, z1,
                 sigmaX, sigmaY, sigmaZ);
      geom->FindIntersection(m_Part, 0, gap,
                 1.0, m_SimpleSlope,          0.0,
                 0.0, m_SimpleIntercept,      0.0,
                 0.0, m_SimpleSlopeSigma,     0.0,
                 0.0, m_SimpleInterceptSigma, 0.0, 
                 x2, y2, z2,
                 sigmaX, sigmaY, sigmaZ);
      dx = z2 - z1;
      dy = sqrt(x2*x2 + y2*y2) - sqrt(x1*x1 + y1*y1);
    }
    else {
      geom->FindIntersection(m_Part, m_Seg, m_LastGap,
                 m_SimpleSlope,          0.0, 1.0,
                 m_SimpleIntercept,      0.0, 0.0,
                 m_SimpleSlopeSigma,     0.0, 0.0,
                 m_SimpleInterceptSigma, 0.0, 0.0, 
                 x1, y1, z1,
                 sigmaX, sigmaY, sigmaZ);
      geom->FindIntersection(m_Part, m_Seg, gap,
                 m_SimpleSlope,          0.0, 1.0,
                 m_SimpleIntercept,      0.0, 0.0,
                 m_SimpleSlopeSigma,     0.0, 0.0,
                 m_SimpleInterceptSigma, 0.0, 0.0, 
                 x2, y2, z2,
                 sigmaX, sigmaY, sigmaZ);
      dx = x2 - x1;
      dy = y2 - y1;
    }
  }
  else {
    if (m_Orient == 0) {
      geom->FindIntersection(m_Part, m_Seg, m_LastGap,
                 0.0, m_SimpleSlope,          1.0,
                 0.0, m_SimpleIntercept,      0.0,
                 0.0, m_SimpleSlopeSigma,     0.0,
                 0.0, m_SimpleInterceptSigma, 0.0, 
                 x1, y1, z1,
                 sigmaX, sigmaY, sigmaZ);
      geom->FindIntersection(m_Part, m_Seg, gap,
                 0.0, m_SimpleSlope,          1.0,
                 0.0, m_SimpleIntercept,      0.0,
                 0.0, m_SimpleSlopeSigma,     0.0,
                 0.0, m_SimpleInterceptSigma, 0.0, 
                 x2, y2, z2,
                 sigmaX, sigmaY, sigmaZ);
      dx = z2 - z1;
      dy = y2 - y1;
    }
    else {
      geom->FindIntersection(m_Part, m_Seg, m_LastGap,
                 m_SimpleSlope,          0.0, 1.0,
                 m_SimpleIntercept,      0.0, 0.0,
                 m_SimpleSlopeSigma,     0.0, 0.0,
                 m_SimpleInterceptSigma, 0.0, 0.0,
                 x1, y1, z1,
                 sigmaX, sigmaY, sigmaZ);
      geom->FindIntersection(m_Part, m_Seg, gap,
                 m_SimpleSlope,          0.0, 1.0,
                 m_SimpleIntercept,      0.0, 0.0,
                 m_SimpleSlopeSigma,     0.0, 0.0,
                 m_SimpleInterceptSigma, 0.0, 0.0,
                 x2, y2, z2,
                 sigmaX, sigmaY, sigmaZ);
      dx = z2 - z1;
      dy = x2 - x1;
    }
  }
  
  dr = sqrt(dx*dx + dy*dy);
 
  return WindowFunc(m_SimpleChi2,dr);
}
 
// Calculate the hit counts and last gap quantities.
void
MucRec2DRoad::CountHits()
{
  m_MaxHitsPerGap = 0;
  m_TotalHits     = 0;
  m_LastGap       = 0;
 
  vector<MucRecHit*>::const_iterator iHit;
  const int ngap = (int)MucID::getGapNum(m_Part);
  vector<int> HitsPerGap;
  for (int gap = 0; gap < ngap; gap++) {
    HitsPerGap.push_back(0);
  }
  
  // Fill HitsPerGap
  for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if ( !(*iHit) ) {
      cout << "MucRec2DRoad::CountHits()-E1 null hit pointer !"  << endl;
      return;
    }
    else {
      int gap = (*iHit)->Gap(); 
      HitsPerGap[gap]++;
      //cout << "gap " << gap << endl;
    }
  } 
 
  // Find the gap from where we should stop count, 
  // namely in front of the gap there are more than 
  // m_MaxNSkippedGaps gaps containing no hit. 
 
  int StopGap = ngap;
   /*
  for (int gap = m_MaxNSkippedGaps; gap < ngap; gap++) { 
    int SubTotalHits = 0; 
    for (int igap = gap-m_MaxNSkippedGaps; igap <= gap; igap++) { 
      SubTotalHits += HitsPerGap[igap]; 
    } 
    if (SubTotalHits == 0 ) {
      StopGap = gap - m_MaxNSkippedGaps;
      break; 
    }
  } 
  */
 
  //cout << "StopGap " << StopGap << endl;
  // Now we might get correct counting on hits, last gap and MaxHitsPerGap 
  for (int gap = 0; gap < StopGap; gap++) { 
    m_TotalHits += HitsPerGap[gap]; 
    if(m_MaxHitsPerGap < HitsPerGap[gap]) m_MaxHitsPerGap = HitsPerGap[gap];  
    if(HitsPerGap[gap] > 0) m_LastGap = gap; 
  } 
} 
 
// Does the Hit exist in the road  .
bool 
MucRec2DRoad::HasHit(MucRecHit* hit) const
{
 
  vector<MucRecHit*>::const_iterator iHit;
  bool HitExist = false;
 
  // Also, if a track hits overlap region of two panels, we avoid
  // to double count hits in two panels
  
  Identifier id = hit->GetID();
 
  for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ ) {
    if ( *iHit ) {  // Check for a null pointer.
      if ( (*iHit)->GetID() == id ) HitExist = true;    
    }
    if (HitExist) break;
  }
  
  return HitExist;
}
 
// Get indices of all hits in the road.
vector<Identifier> 
MucRec2DRoad::GetHitsID() const
{
  vector<Identifier> idCon;
  
  vector<MucRecHit*>::const_iterator iHit;
  for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if (*iHit) {  // Check for a null pointer.
      Identifier id = (*iHit)->GetID();
      idCon.push_back(id);
      /*
      cout << " MucRec2DRoad::HitIndices  gap  orientation  twopack= "
       <<  (*iHit)->ChannelID().Plane() << "   "
       <<  (*iHit)->ChannelID().Orient()  << "   "
       <<  (*iHit)->ChannelID().TwoPack() << endl; 
      */
    }
  }
  return idCon;
}
 
vector<MucRecHit*>
MucRec2DRoad::GetHits() const
{
return m_pHits;
 
}
 
 
 
// Print Hits Infomation.
void
MucRec2DRoad::PrintHitsInfo() const
{
  vector<MucRecHit*>::const_iterator iHit;
  for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
    if (*iHit) {  // Check for a null pointer.
      float xl, yl, zl;
      (*iHit)->GetStrip()->GetCenterPos(xl, yl, zl);
      HepPoint3D vl(xl, yl, zl);
      HepPoint3D vg = (*iHit)->GetGap()->TransformToGlobal(vl);
 
      cout << " orient " << m_Orient
       << " part "   << (*iHit)->Part() 
       << " seg  "   << (*iHit)->Seg()
       << " gap  "   << (*iHit)->Gap()
       << " strip "  << (*iHit)->Strip()
       << " pos ("   << vg.x() << ", " << vg.y() << ", " << vg.z() << ")"  
       << endl;
    }
  }
 
}
 
// Function used to determine the search weight size
float 
MucRec2DRoad::WeightFunc(const float& chisq, const float& distance) const
{
  return 1.0;
}
 
// Function used to determine the search window size
float 
MucRec2DRoad::WindowFunc(const float& chisq, const float& distance) const
{
  return 1.0;
}