HoughHit.cxx

Go to the documentation of this file.

#include "HoughTransAlg/HoughHit.h"
#include "Identifier/MdcID.h"
#include "RawEvent/RawDataUtil.h"
#include "MdcGeom/Constants.h"
#include "TrkBase/TrkPoca.h"
#include "TrkBase/HelixTraj.h"
#include "TrkFitter/TrkCircleTraj.h"
#include "TrkBase/TrkExchangePar.h"
 
#include <iomanip>
#include <string>
 
using namespace std;
 
MdcGeomSvc*         HoughHit::m_mdcGeomSvc = NULL;
MdcCalibFunSvc*     HoughHit::m_mdcCalibFunSvc = NULL;
CgemGeomSvc*        HoughHit::m_cgemGeomSvc = NULL;
CgemCalibFunSvc*    HoughHit::m_cgemCalibFunSvc = NULL;
const MdcDetector*  HoughHit::m_mdcDetector = NULL;
HoughHit::HoughHit(){
    //m_det=Global::m_gm;
    m_cgemCluster=NULL;
    m_mdcDigi=NULL;
    m_cgemMcHit=NULL;
    m_mdcMcHit=NULL;
    //m_hitMap = NULL;
}
 
HoughHit::HoughHit(const MdcDigi* mdcDigi, double bunchTime, int hitID)
{
    m_hitID = hitID;
    m_hitType = MDCHIT;
    m_cgemCluster = NULL;
    m_mdcDigi = mdcDigi;
    m_cgemMcHit=NULL;
    m_mdcMcHit=NULL;
    m_layer = MdcID::layer(mdcDigi->identify());
    m_wire = MdcID::wire(mdcDigi->identify());
    m_flag = flag();
    m_use = 0;
    m_bunchTime = bunchTime;
    m_rawTime = RawDataUtil::MdcTime(mdcDigi->getTimeChannel());
    m_depositEnergy = mdcDigi->getChargeChannel();
    const MdcGeoWire* wire = m_mdcGeomSvc->Wire(m_layer,m_wire);
    m_hitPosition = (wire->Forward() + wire->Backward())/10/2;
    m_westPoint = wire->Forward();
    m_eastPoint = wire->Backward();
    m_driftDist = driftDistance();
    m_residual = 9999;
    //m_hitMap = NULL;
    //m_sLeft = 0;
    //m_zLeft = 0;
    //m_sRight = 0;
    //m_zRight = 0;
    m_trkID.clear();
    m_sz.clear();
    m_pairHit = NULL;
    m_halfCircle = 0;
    m_position.clear();
    m_mdcHit = NULL;
    //MdcHit mdcHit = *(new MdcHit(m_mdcDigi,m_mdcDetector));
    //mdcHit.setCalibSvc(m_mdcCalibFunSvc);
    //mdcHit.setCountPropTime(true);
}
 
HoughHit::HoughHit(const RecCgemCluster* cgemCluster, double bunchTime, int hitID)
{
    m_hitID = hitID;
    m_hitType = CGEMHIT;
    m_cgemCluster = cgemCluster;
    m_mdcDigi=NULL;
    m_cgemMcHit=NULL;
    m_mdcMcHit=NULL;
    m_layer = cgemCluster->getlayerid();
    m_wire = cgemCluster->getsheetid();
    //m_wire = cgemCluster->getclusterid();
    //m_wire = -1;
    m_flag = cgemCluster->getflag();
    m_use = 0;
    m_bunchTime = bunchTime;
    m_rawTime = 0; //FIXME
    m_depositEnergy = cgemCluster->getenergydeposit(); //FIXME
    m_hitPosition.setX((m_cgemGeomSvc->getCgemLayer(m_layer)->getMiddleROfGapD()/10.)*cos(cgemCluster->getrecphi()));
    m_hitPosition.setY((m_cgemGeomSvc->getCgemLayer(m_layer)->getMiddleROfGapD()/10.)*sin(cgemCluster->getrecphi()));
    m_hitPosition.setZ(cgemCluster->getRecZ()/10.);
    m_westPoint = HepPoint3D(0,0,0);
    m_eastPoint = HepPoint3D(0,0,0);
    m_driftDist = 0;
    m_residual = 9999;
    //m_hitMap = NULL;
    //m_sLeft = 0;
    //m_zLeft = 0;
    //m_sRight = 0;
    //m_zRight = 0;
    m_trkID.clear();
    m_sz.clear();
    m_pairHit = NULL;
    m_halfCircle = 0;
    m_position.clear();
    m_mdcHit = NULL;
}
 
HoughHit::HoughHit(const MdcMcHit* mdcMcHit, double bunchTime, int hitID)
{
    m_hitID = hitID;
    m_hitType = MDCMCHIT;
    m_cgemCluster = NULL;
    m_mdcDigi = NULL;
    m_cgemMcHit = NULL;
    m_mdcMcHit = mdcMcHit;
    m_layer = MdcID::layer(mdcMcHit->identify());
    m_wire = MdcID::wire(mdcMcHit->identify());
    m_flag = flag();
    m_use = mdcMcHit->getPositionFlag();
    m_bunchTime = bunchTime;
    m_rawTime = 0;
    m_depositEnergy = mdcMcHit->getDepositEnergy();
    m_hitPosition.setX(mdcMcHit->getPositionX()/10.);
    m_hitPosition.setY(mdcMcHit->getPositionY()/10.);
    m_hitPosition.setZ(mdcMcHit->getPositionZ()/10.);
    m_westPoint = HepPoint3D(0,0,0);
    m_eastPoint = HepPoint3D(0,0,0);
    m_driftDist = mdcMcHit->getDriftDistance()/10.;
    m_residual = 9999;
    //m_hitMap = NULL;
    //m_sLeft = 0;
    //m_zLeft = 0;
    //m_sRight = 0;
    //m_zRight = 0;
    m_trkID.push_back(mdcMcHit->getTrackIndex());
    //cout<<m_layer<<"  "<<mdcMcHit->getTrackIndex()<<endl;
    m_sz.clear();
    m_pairHit = NULL;
    m_halfCircle = 0;
    m_position.clear();
    m_mdcHit = NULL;
}
 
HoughHit::HoughHit(const CgemMcHit* cgemMcHit, double bunchTime, int hitID)
{
    m_hitID = hitID;
    m_hitType = CGEMMCHIT;
    m_cgemCluster = NULL;
    m_mdcDigi = NULL;
    m_cgemMcHit = cgemMcHit;
    m_mdcMcHit = NULL;
    m_layer = cgemMcHit->GetLayerID();
    //m_wire = cgemMcHit->GetParentID();
    //m_flag = cgemMcHit->GetPDGCode();
    m_wire = -1;
    m_flag = 2;
    m_use = 0;
    m_bunchTime = bunchTime;
    m_rawTime = 0;
    m_depositEnergy = 0;
    m_hitPosition.setX((cgemMcHit->GetPositionXOfPrePoint()+cgemMcHit->GetPositionXOfPostPoint())/2/10.);
    m_hitPosition.setY((cgemMcHit->GetPositionYOfPrePoint()+cgemMcHit->GetPositionYOfPostPoint())/2/10.);
    m_hitPosition.setZ((cgemMcHit->GetPositionZOfPrePoint()+cgemMcHit->GetPositionZOfPostPoint())/2/10.);
    m_westPoint = HepPoint3D(0,0,0);
    m_eastPoint = HepPoint3D(0,0,0);
    m_driftDist = 0;
    m_residual = 9999;
    //m_hitMap = NULL;
    //m_sLeft = 0;
    //m_zLeft = 0;
    //m_sRight = 0;
    //m_zRight = 0;
    m_trkID.push_back(cgemMcHit->GetTrackID());
    //cout<<m_layer<<"  "<<cgemMcHit->GetTrackID()<<endl;
    m_sz.clear();
    m_pairHit = NULL;
    m_halfCircle = 0;
    m_position.clear();
    m_mdcHit = NULL;
}
 
HoughHit::HoughHit(const HoughHit& other):
    m_hitID(other.m_hitID),
    m_hitType(other.m_hitType),
    m_cgemCluster(other.m_cgemCluster),
    m_mdcDigi(other.m_mdcDigi),
    m_cgemMcHit(other.m_cgemMcHit),
    m_mdcMcHit(other.m_mdcMcHit),
    m_layer(other.m_layer),
    m_wire(other.m_wire),
    m_flag(other.m_flag),
    m_use(other.m_use),
    m_bunchTime(other.m_bunchTime),
    m_rawTime(other.m_rawTime),
    m_depositEnergy(other.m_depositEnergy),
    m_hitPosition(other.m_hitPosition),
    m_westPoint(other.m_westPoint),
    m_eastPoint(other.m_eastPoint),
    m_driftDist(other.m_driftDist),
    m_residual(other.m_residual),
    //m_hitMap(other.m_hitMap),
    //m_sLeft(other.m_sLeft),
    //m_zLeft(other.m_zLeft),
    //m_sRight(other.m_sRight),
    //m_zRight(other.m_zRight), 
    m_trkID(other.m_trkID),
    m_sz(other.m_sz),
    m_pairHit(other.m_pairHit),
    m_halfCircle(other.m_halfCircle),
    m_position(other.m_position),
    m_mdcHit(other.m_mdcHit)
{}
 
HoughHit& HoughHit::operator=(const HoughHit& other)
{
    m_hitID = other.m_hitID;
    m_hitType = other.m_hitType;
    m_cgemCluster = other.m_cgemCluster;
    m_mdcDigi = other.m_mdcDigi;
    m_cgemMcHit = other.m_cgemMcHit;
    m_mdcMcHit = other.m_mdcMcHit;
    m_layer = other.m_layer;
    m_wire = other.m_wire;
    m_flag = other.m_flag;
    m_use = other.m_use;
    m_bunchTime = other.m_bunchTime;
    m_rawTime = other.m_rawTime;
    m_depositEnergy = other.m_depositEnergy;
    m_hitPosition = other.m_hitPosition;
    m_westPoint = other.m_westPoint;
    m_eastPoint = other.m_eastPoint;
    m_driftDist = other.m_driftDist;
    m_residual = other.m_residual;
    //m_hitMap = other.m_hitMap;
    //m_sLeft = other.m_sLeft;
    //m_zLeft = other.m_zLeft;
    //m_sRight = other.m_sRight;
    //m_zRight = other.m_zRight;
    m_trkID = other.m_trkID;
    m_sz = other.m_sz;
    m_pairHit = other.m_pairHit;
    m_halfCircle = other.m_halfCircle;
    m_position = other.m_position;
    m_mdcHit =other.m_mdcHit;
}
 
/*
void HoughHit::buildMap(int x_bin, double x_min, double x_max, int y_bin, double y_min, double y_max, int nPoint, int charge)
{
    if(m_hitMap!=NULL)delete m_hitMap;
    stringstream ssname;
    ssname<<"id"<<m_hitID<<",layer"<<m_layer<<",wire(sheet)"<<m_wire;
    string sname = ssname.str();
    const char * name = sname.c_str();
    m_hitMap = new TH2D(name,name, x_bin, x_min, x_max, y_bin, y_min, y_max);
    int N = x_bin*nPoint;
    double delta_alpha = (x_max - x_min )/N;
    double alpha = x_min-delta_alpha/2;
    double u = 2*m_hitPosition.x()/(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y()-m_driftDist*m_driftDist);
    double v = 2*m_hitPosition.y()/(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y()-m_driftDist*m_driftDist);
    double w = 2*m_driftDist/(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y()-m_driftDist*m_driftDist);
    for(int i=0;i<N;i++){
        alpha += delta_alpha;
        if(alpha>M_PI)alpha -= M_PI;
        double rho1 = u*cos(alpha) + v*sin(alpha) + w;
        double rho2 = u*cos(alpha) + v*sin(alpha) - w;
        //cout<<alpha<<"  "<<rho1<<"  "<<rho2<<endl;
        double slantOfLine = v*cos(alpha) - u*sin(alpha);
        int chargeOfHitOnCir1 = (slantOfLine/fabs(slantOfLine))*(rho1/fabs(rho1));
        if(chargeOfHitOnCir1 != charge && y_min<(rho1) && (rho1)<y_max){
        //if( y_min<(rho1) && (rho1)<y_max){
            m_hitMap->Fill(alpha,rho1);
        }
        int chargeOfHitOnCir2 = (slantOfLine/fabs(slantOfLine))*(rho2/fabs(rho2));
        if(chargeOfHitOnCir2 != charge && y_min<(rho2) && (rho2)<y_max){
        //if( y_min<(rho2) && (rho2)<y_max){
            if(m_hitType == CGEMHIT)continue;
            //int bin = m_hitMap->FindBin(alpha,rho2);
            //int bin = m_hitMap->FindFixBin(alpha,rho2);
            //if(m_hitMap->GetBinContent(bin)>0)continue;
            m_hitMap->Fill(alpha,rho2);
        }
    }
}
 
void HoughHit::clearMap()
{
    if(m_hitMap!=NULL)delete m_hitMap;
    m_hitMap = NULL;
}
*/
 
int HoughHit::flag()
{
    double nomShift = m_mdcGeomSvc->Layer(m_layer)->nomShift();
    if(nomShift>0) return 1;
    else if(nomShift<0) return -1;
    else return 0;
}
 
double HoughHit::driftTime()
{
    double ToF = 1.e9*sqrt(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y())/Constants::c;
    double Tprop = m_mdcCalibFunSvc->getTprop(m_layer,m_hitPosition.z()*10.);
    double Twalk = m_mdcCalibFunSvc->getTimeWalk(m_layer, m_depositEnergy);
    double T0 = m_mdcCalibFunSvc->getT0(m_layer,m_wire);
    double driftTime = m_rawTime - 1.e9*m_bunchTime - ToF - Tprop - Twalk -T0;
    //cout<<m_rawTime - Twalk - T0 - ToF - Tprop<<endl;
    //cout<<"("<<m_layer<<", "<<m_wire<<"): ";
    //cout<<"bunchT = "<<1.e9*m_bunchTime<<" , rawT "<<m_rawTime<<" , Tof "<<ToF<<" , Tprop "<<Tprop<<" , T0walk "<<Twalk+T0<<" , driftT "<<driftTime<<" , rmid "<<sqrt(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y())<<" , z "<<m_hitPosition.z()*10.0<<endl;
    //cout<<endl;
    return driftTime;
}
 
double HoughHit::driftDistance()
{
    //if(fabs(m_hitPosition.z())>150. || driftTime>1500.)return -1;
    double driftT = driftTime();
    if(fabs(driftT)>1500.||fabs(m_hitPosition.z()>150.)){
        m_use = -1;
        return 9999.;
    }
    double driftDist(-1);
    int lrCalib=2;
    //if (m_flag==1) lrCalib = 0;
    //else if (m_flag==-1) lrCalib = 1;
    double entranceAngle = 0;
    driftDist = 0.1 * m_mdcCalibFunSvc->driftTimeToDist(driftT,m_layer,m_wire,lrCalib,entranceAngle);
    return driftDist;
}
 
void HoughHit::print()
{
    string hitType;
    switch(m_hitType){
        case CGEMHIT    : hitType = "CGEMHIT";      break;
        case MDCHIT     : hitType = "MDCHIT";       break;
        case CGEMMCHIT  : hitType = "CGEMMCHIT";    break;
        case MDCMCHIT   : hitType = "MDCMCHIT";     break;
        default         : hitType = "";
    } 
    cout<<"ID:"        <<setw(6)<<m_hitID;
    cout<<"("          <<setw(2)<<m_layer;
    cout<<","          <<setw(3)<<m_wire<<")  ";
    cout<<"  flag:"      <<setw(3)<<m_flag;
    if(m_hitType==CGEMHIT) {cout<<" clusterId "<<m_cgemCluster->getclusterid()<<" pos: "<<m_hitPosition;}
    if(m_hitType==CGEMHIT||m_hitType==MDCHIT){
        if(m_pairHit!=NULL){
            cout<<"  half:"<<setw(2)<<m_pairHit->getHalfCircle();
            cout<<"  trkID:"<<setw(3)<<(m_pairHit->getTrkID())[0];
            cout<<"  McHit:("<<setw(10)<<m_pairHit->getHitPosition().x()<<", "<<setw(10)<<m_pairHit->getHitPosition().y()<<", "<<setw(10)<<m_pairHit->getHitPosition().z()<<") ";
            if(m_hitType==MDCHIT){
                cout<<"  McLR:"<<setw(3)<<m_pairHit->getMdcMcHit()->getPositionFlag();
                cout<<" drift:"<<setw(10)<<m_pairHit->getMdcMcHit()->getDriftDistance();
            }
            //if(m_hitType==MDCHIT)
            //cout<<" trkID:"<<setw(3)<<m_pairHit->getMdcMcHit()->getTrackIndex();
        }
    }
    if(m_hitType==CGEMMCHIT||m_hitType==MDCMCHIT){
        cout<<"  half:"<<setw(2)<<getHalfCircle();
        cout<<"  trkID:"<<setw(3)<<m_trkID[0];
        cout<<"  McHit:("<<setw(10)<<m_hitPosition.x()<<", "<<setw(10)<<m_hitPosition.y()<<", "<<setw(10)<<m_hitPosition.z()<<") ";
        //if(m_hitType==MDCMCHIT)
        //cout<<" trkID:"<<setw(3)<<m_mdcMcHit->getTrackIndex();
        if(m_hitType==MDCMCHIT){
            cout<<"  McLR:"<<setw(3)<<m_mdcMcHit->getPositionFlag();
            cout<<" drift:"<<setw(10)<<m_mdcMcHit->getDriftDistance();
        }
    }
 
    //cout<<setw(7)<<"hitID:"     <<setw(6)<<m_hitID
        //<<setw(10)<<"hitType:"  <<setw(10)<<hitType
        //<<setw(8)<<"layer:"     <<setw(4)<<m_layer
        //<<setw(6)<<"wire:"      <<setw(6)<<m_wire
        //<<setw(10)<<"flag:"     <<setw(10)<<m_flag
        //<<setw(5)<<"use:"      <<setw(3)<<m_use
        //<<endl
        //<<setw(12)<<"bunchTime:"<<setw(12)<<m_bunchTime
        //<<setw(12)<<"rawTime:"  <<setw(12)<<m_rawTime
        //<<setw(12)<<"charge:"   <<setw(12)<<m_depositEnergy
        //<<setw(12)<<"driftTime:"<<setw(12)<<driftTime()
        //<<setw(11)<<"driftDist:"<<setw(12)<<m_driftDist
        //<<setw(12)<<"residual:" <<setw(12)<<m_residual
        //<<endl
        //<<setw(10)<<"sLeft:"    <<setw(10)<<m_sLeft
        //<<setw(10)<<"zLeft:"    <<setw(10)<<m_zLeft
        //<<setw(10)<<"sRight:"   <<setw(10)<<m_sRight
        //<<setw(10)<<"zRight:"   <<setw(10)<<m_zRight
        //<<endl
        //<<endl
        //cout<<m_trkID.size()<<" trkID:";
        //for(int i=0;i<m_trkID.size();i++)cout<<setw(4)<<m_trkID[i]<<endl;
        //;
        //if(m_layer>3)cout<<"TrkID:"<<setw(4)<<getDigi()->getTrackIndex(); 
        //cout<<endl;
    //;
    cout<<endl;
}
 
 
// erase trk id and residual from the Hough hit
void HoughHit::dropTrkID(int trkID)
{
    vector<int>::iterator it0 = m_trkID.begin();
    vector<int>::iterator result = find(m_trkID.begin(),m_trkID.end(),trkID); 
    if(result!=m_trkID.end()) {
        m_trkID.erase(result);
        vector<double>::iterator itRes = m_vecResid.begin()+(result-it0);
        if(itRes!=m_vecResid.end()) m_vecResid.erase(itRes);
    }
}
 
vector<HepPoint3D>& HoughHit::VHitPosition(HoughTrack* track)
{
    if(getFlag()==0)return m_position;
    m_position.clear();
    //if(getPairHit()==NULL)return m_position;
    //if(getPairHit()->getHalfCircle()!=1)return m_position;
    double z(0);
    double xc = track->center().x();
    double yc = track->center().y();
    double rTrack = track->radius();// signed FIXME
    if(getHitType()==HoughHit::MDCHIT){
        resetSZ();
        double drift = getDriftDist();
        const MdcGeoWire* wire = getMdcGeomSvc()->Wire(getLayer(),getWire());
        HepPoint3D westPoint = wire->Forward();
        HepPoint3D eastPoint = wire->Backward();
        double xEast = eastPoint.x()/10.0;
        double xWest = westPoint.x()/10.0;
        double yEast = eastPoint.y()/10.0;
        double yWest = westPoint.y()/10.0;
        double zEast = eastPoint.z()/10.0;
        double zWest = westPoint.z()/10.0;
        //cout<<"wast:x,y,z: "<<xWest<<", "<<yWest<<", "<<zWest<<endl;
        //cout<<"east:x,y,z: "<<xEast<<", "<<yEast<<", "<<zEast<<endl;
        //cout<<"Xc, Yc, Rc: "<<xc<<", "<<yc<<", "<<rTrack<<endl;
        double west2east = sqrt((xEast-xWest)*(xEast-xWest)+(yEast-yWest)*(yEast-yWest));
 
        double slope = (yEast-yWest)/(xEast-xWest);  
        double intercept = (yWest-slope*xWest+yEast-slope*xEast)/2;
        double a = 1+slope*slope;
        double b = -2*(xc+slope*yc-slope*intercept);
        double c1 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack+drift)*(rTrack+drift);
        double c2 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack-drift)*(rTrack-drift);
        //double c1 = intercept*(2*yc-intercept)+drift*(drift+rTrack);
        //double c2 = intercept*(2*yc-intercept)+drift*(drift-rTrack);
        double delta1 = (b*b-4*a*c1);
        double delta2 = (b*b-4*a*c2);
        //cout<<"a,b,c1,c2: "<<a<<", "<<b<<", "<<c1<<", "<<c2<<endl;
        //cout<<"delta: "<<delta1<<", "<<delta2<<endl;
        if(delta1>=0){
            double x1 = (-b+sqrt(delta1))/(2*a);
            double x2 = (-b-sqrt(delta1))/(2*a);
            double y1 = slope*x1+intercept;
            double y2 = slope*x2+intercept;
            if((x1-xWest)*(x1-xEast)<0){
                double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
                z = zWest + l/west2east*fabs((zEast-zWest));
                HepPoint3D position(x1,y1,z);
                m_position.push_back(position);
            }
            if((x2-xWest)*(x2-xEast)<0){
                double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
                z = zWest + l/west2east*fabs((zEast-zWest));
                HepPoint3D position(x2,y2,z);
                m_position.push_back(position);
            }
        }
 
        if( delta2>=0){
            double x1 = (-b+sqrt(delta2))/(2*a);
            double x2 = (-b-sqrt(delta2))/(2*a);
            double y1 = slope*x1+intercept;
            double y2 = slope*x2+intercept;
            if((x1-xWest)*(x1-xEast)<0){
                double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
                z = zWest + l/west2east*fabs((zEast-zWest));
                HepPoint3D position(x1,y1,z);
                m_position.push_back(position);
            }
            if((x2-xWest)*(x2-xEast)<0){
                //double l = zWest + sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
                //z = l/west2east*fabs((zEast-zWest));
                double l =  sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
                z = zWest + l/west2east*fabs((zEast-zWest));
                HepPoint3D position(x2,y2,z);
                m_position.push_back(position);
            }
        }
    }
    return m_position;
}
///*
double HoughHit::residual(HoughTrack* track)
{
    double residual(9999.);
    if(m_flag==0){
        HepPoint3D hitPoint = getHitPosition();
        if(getHitType() == HoughHit::MDCHIT || getHitType() == HoughHit::MDCMCHIT){
            //HepPoint3D circleCenter = m_helix.center();
            HepPoint3D circleCenter = track->center();
            //double distance = circleCenter.perp(hitPoint);
            double distance = (circleCenter-hitPoint).perp();
            //double Rc = m_helix.radius();
            double Rc = fabs(track->radius());
            double driftDist(0);
            //if(getHitType() == HoughHit::MDCMCHIT)driftDist = 0;
            if(getHitType() == HoughHit::MDCHIT)driftDist = getDriftDist();
            //residual = fabs(distance - Rc) - driftDist;
            residual = driftDist-fabs(distance - Rc);
            //cout<<"driftDist, distance, Rc = "<<driftDist<<", "<<distance<<", "<<Rc<<endl;
            //double distance = (m_helix.center()).perp(hit.getHitPosition());
            //res = fabs(distance - m_helix.radius()) - hit.getDriftDist();
        }else if(getHitType() == HoughHit::CGEMHIT || getHitType() == HoughHit::CGEMMCHIT){
            double rCgem = hitPoint.perp();
            //double phiTrkFlt = m_helix.IntersectCylinder(rCgem);
            double phiTrkFlt = track->IntersectCylinder(rCgem);
            phiTrkFlt = track->judgeHalf(this)*phiTrkFlt;
            //HepPoint3D crossPoint = m_helix.x(phiTrkFlt);
            HepPoint3D crossPoint = track->x(phiTrkFlt);
            double phiCrossPoint = crossPoint.phi();
            double phiMeasure = hitPoint.phi();
            residual = phiMeasure - phiCrossPoint;
            while(residual<-M_PI)residual += 2*M_PI;
            while(residual> M_PI)residual -= 2*M_PI;
            residual = rCgem*residual;
        }
    }else{
        double zTruth = getPairHit()->getHitPosition().z();
        double minRes(9999);
        if(getHitType() == HoughHit::MDCHIT || getHitType() == HoughHit::MDCMCHIT){
            vector<HepPoint3D> position = VHitPosition(track);
            for(vector<HepPoint3D>::iterator posIt = position.begin();posIt!=position.end();posIt++){
                HepPoint3D point = *posIt;
                double s = track->flightArc(point);
                double zTrk = track->dz()+s*track->tanl();
                double res = point.z() - zTrk;
                //res = point.z() - zTruth;
                if(fabs(res)<minRes){
                    minRes=fabs(res);
                    residual = res;
                }
            }
        }else if(getHitType() == HoughHit::CGEMHIT || getHitType() == HoughHit::CGEMMCHIT){
            double s = track->flightArc(m_hitPosition);
            double zTrk = track->dz()+s*track->tanl();
            residual = m_hitPosition.z() - zTrk;
            //residual = m_hitPosition.z() - zTruth;
        }
    }
    return residual;
}
 
void HoughHit::updateVHit(HoughTrack* track)
{
    if(m_flag!=0&&m_hitType==HoughHit::MDCHIT){
        double minRes(9999);
        HepPoint3D hitPoint(m_hitPosition);
        vector<HepPoint3D> position = VHitPosition(track);
        for(vector<HepPoint3D>::iterator posIt = position.begin();posIt!=position.end();posIt++){
            HepPoint3D point = *posIt;
            double s = track->flightArc(point);
            double zTrk = track->dz()+s*track->tanl();
            double res = point.z() - zTrk;
            if(fabs(res)<minRes){
                minRes=fabs(res);
                hitPoint = point;
            }
        }
        m_hitPosition = hitPoint;
        m_driftDist = driftDistance();
    }
}
//*/
double HoughHit::residual(HoughTrack* track, HepPoint3D& positionOntrack, HepPoint3D& positionOnDetector)
{
    double residual(9999.);
    HepPoint3D pivot(0,0,0);
    track->pivot(pivot);
    if(m_hitType==MDCHIT){
        const MdcDetector* mdcDetector(m_mdcDetector);
        MdcHit mdcHit(m_mdcDigi,mdcDetector); 
        const Trajectory* wire(mdcHit.hitTraj());
 
        double d0       = -track->dr();
        double phi0     = track->phi0()+M_PI/2;
        double omega    = track->kappa()/fabs(track->alpha());
        double z0       = track->dz();
        double tanl     = track->tanl();
 
        //HepVector a(5,0);
        //a(1) = d0;
        //a(2) = phi0;
        //a(3) = omega;
        //a(4) = z0;
        //a(5) = tanl;
        //HepSymMatrix Ea(5,0);
        //const HepVector pvec(a);
        //const HepSymMatrix pcov(Ea);
        //const HepPoint3D refpoint(track->pivot());
        //const HelixTraj trk(pvec, pcov, -99999., 99999., refpoint);
 
        double hitWireLength = m_hitPosition.z()-m_westPoint.z();
        const TrkExchangePar trkExchangePar(d0,phi0,omega,z0,tanl);
        double trkFlightLength = track->flightLength(m_hitPosition);
        const HelixTraj helixTraj(trkExchangePar);
        TrkPoca trkPoca(helixTraj,trkFlightLength,*wire,hitWireLength);
        trkFlightLength = trkPoca.flt1();
        hitWireLength = trkPoca.flt2();
        double doca = trkPoca.doca();
        //double trkFlightArc = track->flightArc(m_hitPosition);
        //const TrkCircleTraj circleTraj(trkExchangePar);
        //TrkPoca trkPoca(circleTraj,trkFlightArc,*wire,hitWireLength);
        //trkFlightArc = trkPoca.flt1();
        //hitWireLength = trkPoca.flt2();
        //double doca = trkPoca.doca();
 
        Hep3Vector wireDirection, trackDirection;
        //HepPoint3D hitOnWire, hitOntrack;
        wire->getInfo(hitWireLength, positionOnDetector, wireDirection);
        helixTraj.getInfo(trkFlightLength, positionOntrack, trackDirection);
        residual = m_driftDist - fabs(doca);
        //cout<<"("<<m_layer<<","<<m_wire<<") ";
        //cout<<setw(12)<<m_driftDist
            //<<setw(12)<<doca
            //<<setw(12)<<residual
            //<<endl;
 
        //cout<<"("<<m_layer<<","<<m_wire<<") "
            //<<setw(12)<<getPairHit()->getHitPosition().z()
            //<<setw(12)<<positionOntrack.z()
            //<<setw(12)<<positionOnDetector.z()
            //<<setw(12)<<tanl/sqrt(1+tanl*tanl)*trkFlightLength+z0
            //<<setw(12)<<tanl*(track->flightArc(m_hitPosition))+z0
            //<<endl;
    }
    else if(m_hitType==CGEMHIT){
        double rCgem = m_hitPosition.perp();
        //double phiTrkFlt = m_helix.IntersectCylinder(rCgem);
        double phiTrkFlt = track->IntersectCylinder(rCgem);
        phiTrkFlt = track->judgeHalf(this)*phiTrkFlt;
        //HepPoint3D crossPoint = m_helix.x(phiTrkFlt);
        positionOntrack = track->x(phiTrkFlt);
        positionOnDetector = track->x(phiTrkFlt);
        if(m_flag==0){
            double phiCrossPoint = positionOntrack.phi();
            double phiMeasure = m_hitPosition.phi();
            double dphi = phiMeasure - phiCrossPoint;
            while(dphi<-M_PI)dphi += 2*M_PI;
            while(dphi> M_PI)dphi -= 2*M_PI;
            residual = rCgem*dphi;
        }
        else{
            residual = (m_cgemCluster->getrecv() - m_cgemGeomSvc->getReadoutPlane(m_layer,m_wire)->getVFromPhiZ(positionOntrack.phi(),positionOntrack.z()*10.))/10.;
        }
    }
    return residual;
}