CgemHitOnTrack.cxx

Go to the documentation of this file.

// $Id: CgemHitOnTrack.cxx,v 1.4 2018/01/04 00:11:49 huangzhen Exp $
//
#include "CgemData/CgemHitOnTrack.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/Matrix.h"
#include "TROOT.h"
#include "CLHEP/Geometry/Point3D.h"
 
#include "MdcGeom/Constants.h"
#include "MdcGeom/BesAngle.h"
// MdcGeom needed to verify if hit is inside of chamber...
// and to find the trajectory describing the hit, i.e. wire
#include "MdcGeom/MdcLayer.h"
#include "MdcGeom/MdcSWire.h"
#include "TrkBase/TrkPoca.h"
#include "TrkBase/TrkDifTraj.h"
#include "TrkBase/TrkRep.h"
#include "TrkBase/TrkSimpTraj.h"
#include "TrkBase/TrkRecoTrk.h"
#include "MdcGeom/EntranceAngle.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/Service.h"
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrackUtil/Helix.h"
using std::cout;
using std::endl;
 
const CgemGeomSvc* CgemHitOnTrack::_cgemGeomSvc = NULL;
const CgemCalibFunSvc* CgemHitOnTrack::_cgemCalibFunSvc = NULL;
//-------------
// Constructors
//-------------
CgemHitOnTrack::CgemHitOnTrack(//const TrkFundHit& fundHit,
            const RecCgemCluster& baseHit,
            double t0)
: TrkHitOnTrk(NULL,10.e-4),
    //: TrkHitOnTrk(&fundHit,10.e-4),
    //_hitTraj(baseHit.hitTraj()), //yzhang FIXME add hitTraj of this Cluster ?
    //_fitTime(baseHit.rawTime()-t0*1e-9),
    _fitTime(t0*1e-9),//yzhang FIXME
    _dHit(&baseHit)
{
    ICgemGeomSvc* icgemGeomSvc = NULL;
    StatusCode sc = Gaudi::svcLocator()->service("CgemGeomSvc", icgemGeomSvc);
    CgemGeomSvc* cgemGeomSvc = dynamic_cast<CgemGeomSvc*> (icgemGeomSvc);
    if(sc==StatusCode::SUCCESS && cgemGeomSvc!=NULL){
        _clusterr = cgemGeomSvc->getCgemLayer(baseHit.getlayerid())->getMiddleROfGapD()/10.;//r[baseHit.getlayerid()];
        //std::cout<<__FILE__<<"   "<<__LINE__<<" _clusterr  "<<_clusterr<<std::endl;
    }else{
        _clusterr = -999;
        //std::cout<<__FILE__<<"   "<<__LINE__<<" _clusterr Could not find CgemGeomSvc "<<_clusterr<<std::endl;
    }
    //hitTraj not implemented
    // need to flag somehow that that we haven't computed things yet...
    //   now, we know that _drift[0] is for ambig <0, and _drift[1] is ambig >0
    //   and _drift is signed according to ambig. Thus _drift[0] should be -tive
    //   and _drift[1] should be +tive. To indicate this are not yet initialized,
    //   put 'impossible' values here...
    //_drift[0] = +9999;
    //_drift[1] = -9999;
    setHitResid(-21212121.0);
    setHitRms( 0.02 );
    //setHitLen(0.5 * baseHit.layer()->zLength());//yzhang FIXME
    setFltLen(0.);
    //layer 
    double length = cgemGeomSvc->getCgemLayer(baseHit.getlayerid())->getLengthOfCgemLayer();
    //_startLen = hitTraj()->lowRange() - 5.;//lowRange?
    //_endLen   = hitTraj()->hiRange()  + 5.;
    _startLen = -1.*length/2. - 5.;
    _endLen   = length/2.  + 5.;
 
    //double r[3] = {7.9838 ,12.5104 ,16.7604};//cm 
}
 
TrkHitOnTrk*
CgemHitOnTrack::clone(TrkRep *rep, const TrkDifTraj *trkTraj) const
{
    return new CgemHitOnTrack(*this,rep,trkTraj);
}
//CgemHitOnTrack::CgemHitOnTrack(const TrkFundHit *fundHit, int ambig, double fittime,
//                             int layer, int wire)
//  : TrkHitOnTrk(fundHit,10.e-4),
//    _ambig(ambig), _fitTime(fittime)
//{
//  cout << "using old CgemHitOnTrack constructor..." << endl;
 
//  IService* ser;
//  StatusCode sc = Gaudi::svcLocator()->getService("MdcGeomSvc",ser);
//  if (!sc.isSuccess())
//    cout <<"CgemHitOnTrack::Could not open Geometry Service"<<endl;
//  MdcGeomSvc *m_gmSvc = dynamic_cast<MdcGeomSvc*> (ser);
//  if(!m_gmSvc) cout <<"CgemHitOnTrack::Could not open Geometry Service"<<endl;
 
 
 
//  _drift[0] = +9999;
//  _drift[1] = -9999;
//  _hitTraj = gblEnv->getMdc()->getMdcLayer(layer)->makeHitTrajInGlobalCoords(wire);
//    _hitTraj = m_gmSvc->Layer(layer)->makeHitTrajInGlobalCoords(wire);
//  setHitResid(-21212121.0);
//  setHitRms( 0.01 );
//  setHitLen(2.);
//  setFltLen(0.);
//  _startLen = hitTraj()->lowRange() - 5.;
//  _endLen   = hitTraj()->hiRange()  + 5.;
//}
 
CgemHitOnTrack::CgemHitOnTrack(const CgemHitOnTrack &hot,TrkRep *newRep,
            const TrkDifTraj* trkTraj, const RecCgemCluster *hb)
: TrkHitOnTrk(hot,newRep,trkTraj)
{
    _hitTraj  = hot._hitTraj;
    _fitTime  = hot._fitTime;
    _drift[0] = hot._drift[0];
    _drift[1] = hot._drift[1];
    _startLen = hot._startLen;
    _endLen   = hot._endLen;
    _dHit = (hb==0?hot._dHit:hb);
    _clusterr = hot._clusterr;
}
 
CgemHitOnTrack::~CgemHitOnTrack()
{ ; }
 
void
CgemHitOnTrack::setT0(double t0)
{ 
    //_fitTime= _dHit->rawTime()-t0*1e-9; 
    _fitTime = t0*1e-9;//FIXME 
}
 
double
CgemHitOnTrack::dcaToWire() const
{
    double dca = -9999.;
    if ( getParentRep() == 0 ) {
        //     cout << "no parent rep" << endl;
        return dca;
    }
    // WARNING: cannot use the internal _poca, as it lags one iteration
    //          behind the fit... therfore use _EXTERNAL_ residual
    //if (isActive())  {  // FIXME: currently can only use 'resid()' if isActive..
    //  dca = resid()+drift();
    //} else {
    //  TrkPoca poca(getParentRep()->traj(), fltLen(), *hitTraj(), hitLen(),
    //      _tolerance);
    //  if (poca.status().success()) dca = poca.doca();
    //}
    //yzhang FIXME
    return dca;
}
 
double
CgemHitOnTrack::entranceAngleHit() const
{
    static Hep3Vector dir;
    static HepPoint3D pos;
    if (getParentRep() == 0) return 0.;
    getParentRep()->traj().getInfo(fltLen(), pos, dir);
 
    return BesAngle(dir.phi() - pos.phi());
}
 
double
CgemHitOnTrack::entranceAngle() const
{
    static Hep3Vector dir;
    static HepPoint3D pos;
    if (getParentRep() == 0) return 0.;
    getParentRep()->traj().getInfo(fltLen(), pos, dir);
 
    return entranceAngle(pos, dir);
}
 
double
CgemHitOnTrack::entranceAngle(const HepPoint3D pos, const Hep3Vector dir) const
{
    double angle = EntranceAngle(dir.phi() - _dHit->getrecphi());
    return angle;
}
 
//unsigned
//CgemHitOnTrack::layerNumber() const
//{
// return layernumber();
//}
 
double
CgemHitOnTrack::dipAngle() const
{
    return getParentRep()==0?0:Constants::pi/2-getParentRep()->traj().direction(fltLen()).theta();
}
 
TrkErrCode
CgemHitOnTrack::updateMeasurement(const TrkDifTraj* traj, bool maintainAmb)
{
    //std::cout<<__FILE__<<"   "<<__LINE__<<" WARNING CgemHitOnTrack::updateMeasurement not implemented  "<<std::endl;
    traj = (traj!=0?traj:&(getParentRep()->traj()));//track traj
    //TrkErrCode status=updatePoca(traj,maintainAmb);
    //if (status.failure()) {
    //  std::cout<<" ErrMsg(warning) " << "CgemHitOnTrack::updateMeasurement failed " << status << std::endl;
    //  return status;
    //}
    //assert (_poca!=0);
    double dca=-999;//FIXME _poca->doca();
    ////yzhang FIXME
    bool forceIteration = true;//(maintainAmb&&ambig()!=0)?false:updateAmbiguity(dca);
    //if( whatView()==0) {
    BesAngle phiCluster = BesAngle(_dHit->getrecphi());//-pi 2 pi
    double flt=0;
    const HepVector par = getParentRep()->traj().localTrajectory(0.,flt)->parameters()->parameter();
    //cout<<par<<endl;
    double d0 = par[0];
    double phi0 = par[1];
    double omega = par[2];
    double radius = fabs(omega)<0.00001 ? -9999: 1./fabs(omega);
    double l = d0 + radius;
    double r = _clusterr;
    double xc = l * cos(phi0 - Constants::halfPi);
    double yc = l * sin(phi0 - Constants::halfPi);
    double phi_center = atan2(yc,xc);
 
    double cosPhi = (radius * radius + l * l  - r * r) / (2 * radius * l);
    if(cosPhi < -1 || cosPhi > 1) return TrkErrCode(TrkErrCode::fail);
    double dPhi = phi_center - acos(cosPhi) - (phi0-Constants::halfPi);
    //std::cout<<__FILE__<<"   "<<__LINE__<<" radius "<<radius<<" l "<<l<<" r "<<r<<" xc "<<xc<<" yc "<<yc <<" cosPhi "<<cosPhi<<" phi center "<<phi_center<<" phi0 - pi/2 "<<phi0-Constants::halfPi<<" "<<std::endl; 
    BesAngle phi(dPhi);
    //if(dPhi < -M_PI) dPhi += 2 * M_PI;
 
    double x = radius*sin(dPhi) - (radius + d0)*sin(phi0);
    double y = -1.*radius*cos(dPhi) + (radius + d0)*sin(phi0);
 
    //double phiTrack = atan2(y,x);
    BesAngle phiTrack = BesAngle(atan2(y,x));
 
    //std::cout<<__FILE__<<"   "<<__LINE__<<" phiCluster "<<phiCluster<<" phiTrack "<<phiTrack<<" x "<<x<<" y "<<y<<std::endl; 
    dca = (phiCluster - phiTrack).rad();
    //}
    //return TrkErrCode(TrkErrCode::fail);//not used for fit
    return TrkErrCode(TrkErrCode::succeed);//not used for fit
    //TrkErrCode(TrkErrCode::succeed, 11, "Ambiguity flipped");
}
 
void
CgemHitOnTrack::updateCorrections()
{
    std::cout<<__FILE__<<"   "<<__LINE__<<" WARNING CgemHitOnTrack::updateCorrections not implemented  "<<std::endl;
    const TrkRep* tkRep = getParentRep();
    assert(tkRep != 0);
    static HepPoint3D pos; static Hep3Vector dir;
    _trkTraj->getInfo(fltLen(), pos, dir);
 
    // for bes3 cosmic test
    //int wireAmb = wireAmbig();
    double tof = tkRep->arrivalTime(fltLen());
    // at this point, since dcaToWire is computed, _ambig must be either -1 or +1
    //assert( ambig() == -1 || ambig() == 1 );
    //double eAngle = entranceAngle(pos, dir);//2011-11-22
    //double dAngle = Constants::pi/2 - dir.theta();
    //double z = pos.z();
    // note the special case for INCOMING tracks:
    //    wire left of track implies hit on the right side of wire
    //int wireAmb= fabs(eAngle)<Constants::pi/2?_ambig:-_ambig; //yzhang delete 2012-07-17 
 
    // provide the underlying hit with the *external* information
    // needed to compute the drift distance, i.e. those numbers that
    // the hit cannot figure out by itself...
    double dist = -999;//FIXME
    //_dHit->driftDist(tof, wireAmb, eAngle, dAngle, z); 
    //<<" dir "<<dir<<" pos "<<pos
    //assert(dist>0);//yzhang 2011-12-05 delete
 
    //FIXME//_fitTime = _dHit->driftTime(tof,z);
    //
    //std::cout<<" CgemHitOnTrack  ("<<layernumber()<<","<<wire()<<") "<<mdcHit()->isCosmicFit()<<" fltLen "<<fltLen() <<" eAngle "<<eAngle<<" ambig "<<ambig()<<" wireAmb "<<wireAmb<<" tof "<<tof<<" dd "<<dist<<" dt "<<_fitTime<<std::endl;
    //_drift[ambig()<0?0:1] = ambig() * dist;
    //assert( driftCurrent() );
 
    double newSig = -999;//_dHit->sigma(dist, wireAmb, eAngle, dAngle, z);
 
    //assert(newSig>0);//yzhang 2011-12-05 delete
    setHitRms(newSig);
}
 
bool
CgemHitOnTrack::timeResid(double &t, double &tErr) const
{
    //double v = driftVelocity();
    //if (v <= 0) return false;
    //t = (fabs(drift())-fabs(dcaToWire()))/v;
    t = -999;//FIXME
    //tErr= hitRms()/v;
    return true;
}
 
bool
CgemHitOnTrack::timeAbsolute(double &t, double &tErr) const
{
    double tresid(-1.0);
    if(timeResid(tresid,tErr)){
        // add back the track time
        t = tresid + getParentRep()->parentTrack()->trackT0();
        return true;
    } else
      return false;
}
 
TrkEnums::TrkViewInfo
CgemHitOnTrack::whatView() const
{
    //enum TrkViewInfo {noView=-1,xyView=0, zView, bothView };
    TrkEnums::TrkViewInfo view = TrkEnums::TrkViewInfo(-1);
    if(_dHit->getflag()==0) view = TrkEnums::TrkViewInfo(0);
    else if(_dHit->getflag()==1) view = TrkEnums::TrkViewInfo(1);
    else if(_dHit->getflag()==2) view = TrkEnums::TrkViewInfo(2);
 
    return view;
}
 
/*
   int
   CgemHitOnTrack::layernumber() const
   {
   return _dHit->layernumber();
   }
 
   const MdcLayer*
   CgemHitOnTrack::layer() const
   {
   return _dHit->getlayerid();
   }
 
   int
   CgemHitOnTrack::wire() const
   {
   return _dHit->wirenumber();
   }
 
   int
   CgemHitOnTrack::whichView() const
   {
   return _dHit->whichView();
   }
 
   double
   CgemHitOnTrack::rawTime() const
   {
   return _dHit->rawTime();
   }
*/
 
   double
   CgemHitOnTrack::charge() const
   {
   return _dHit->getenergydeposit();
   }
 
const Trajectory*
CgemHitOnTrack::hitTraj() const
{
    //return _hitTraj;
    std::cout<<__FILE__<<"   "<<__LINE__<<" CgemHitOnTrack::hitTraj() not implemented, return NULL "<<std::endl;
    return NULL;
}
 
// Replace underlying hit pointer (base class doesn't own it)
 
void 
CgemHitOnTrack::changeBase(RecCgemCluster* newBase)
{
    _dHit = newBase;
}
 
/*
   CgemHitOnTrack::isBeyondEndflange() const { 
   std::cout<<__FILE__<<" "<<__LINE__
   <<" hitLen "<<hitLen()
   <<" startLen "<<_startLen
   <<" endLen "<<_endLen
   <<  std::endl;
   return (hitLen() < _startLen || hitLen() > _endLen); 
   }
   */
 
void CgemHitOnTrack::print(std::ostream&) const{
    std::cout<<" CgemHOT id "<<_dHit->getclusterid()<<" layer "<<_dHit->getlayerid()<<" recphi "<< _dHit->getrecphi()<<" recv "<< _dHit->getrecv()<<" recz "<<_dHit->getRecZ() << std::endl;
}
 
TrkErrCode
//CgemHitOnTrack::getFitStuff(HepVector &derivs, double &deltaChi ) const
CgemHitOnTrack::getFitStuff(HepVector &derivs, HepVector &derivs2, double & sigma1, double & sigma2, double &deltaChi1, double &deltaChi2 ) const
{
    //std::cout<<__FILE__<<"   "<<__LINE__<<"  CgemHitOnTrack::getFitStuff implemented"<<std::endl;
    if (_poca==0 || _poca->status().failure()) {
        return TrkErrCode(TrkErrCode::fail);
    }
    double vec_phi = _dHit->getrecphi() ;
    double vec_z = _dHit->getRecZ() ;
    double flt=0;
    HepVector trkpar = getParentRep()->traj().localTrajectory(0.,flt)->parameters()->parameter(); //get para??????
    IMdcUtilitySvc*   m_mdcUtilitySvc;
    Gaudi::svcLocator()->service("MdcUtilitySvc",m_mdcUtilitySvc);
    if (m_mdcUtilitySvc == NULL) cout<< " Can not load MdcUtilitySvc"<<endl;
    HepVector trkparbes = m_mdcUtilitySvc->patPar2BesPar(trkpar);
    trkparbes[2]=-trkparbes[2]; // modification for the reversed sign in TrackUtil/Helix
 
    HepPoint3D pivot(0., 0., 0.);
    Helix cgem_helix_0(pivot, trkparbes);
    HepMatrix J_dmdx(2,3,0), J_dxda(3,5,0), J(2,5,0) ;
    HepSymMatrix V(2,0), W(2,0) ;
    HepVector a_new(5), dm(2);
    J_dmdx(1,1)= 1. ;
    J_dmdx(2,3)= 1. ;
    int layer = _dHit->getlayerid() ;
    double rl, R, x_reso, v_reso, a_stero;
 
    if(layer==0) { rl=79.838/10. ; R=87.5026/10.; x_reso=0.1372/10.; v_reso=0.1273/10.;a_stero=45.94*3.1415926/180;}
    else if (layer==1) { rl=125.104/10. ; R=132.7686/10.; x_reso=0.1476/10.;v_reso=0.1326/10.;a_stero=31.10*3.1415926/180;}
    else { rl=167.604/10.; R=175.2686/10.; x_reso=0.1412/10.;v_reso=0.1378/10.;a_stero=32.99*3.1415926/180;} //cm 
    double d_phi = cgem_helix_0.IntersectCylinder(rl); //const can only call const 
    HepPoint3D pos = cgem_helix_0.x(d_phi);
    double x = pos.x();
    double y = pos.y();
    double z = pos.z();
    double r2 = x*x+y*y;
    double tmp_phi = atan2(y, x); //[-pi,pi), return azimuthal angle
    double dphi = vec_phi - tmp_phi;
    if(dphi < -M_PI) dphi += 2*M_PI;
    if(dphi >= M_PI) dphi -= 2*M_PI;
    dm[0]=-dphi*rl*sin(vec_phi);
    dm[1]=vec_z-z;
    J_dxda=cgem_helix_0.delXDelA(d_phi); //3X5; derivative of x,y,z to 5 paras
    J=J_dmdx*J_dxda;
    for (int i=0;i<5;i++){
        derivs2[i]=J(1,i) ;
        derivs2[i+5]=J(2,i) ;
    }
    V(1,1)=rl*rl/R/R*sin(vec_phi)*sin(vec_phi)*x_reso*x_reso;
    V(2,1)=rl/R*sin(vec_phi)/tan(a_stero)*x_reso*x_reso;
    V(2,2)=v_reso*v_reso/sin(a_stero)/sin(a_stero)+x_reso*x_reso/tan(a_stero)/tan(a_stero);
    W=V;
 
    int ifail=99;
    V.invert(ifail);
    if (ifail) {
        return TrkErrCode(TrkErrCode::fail);
    }
    J=V*J; 
 
    for (int i=0;i<5;i++){
        derivs[i]=J(1,i) ;
        derivs[i+5]=J(2,i) ;
    }
    assert (_trkTraj == &(getParentRep()->traj()));
    sigma1=sqrt(W(1,1));
    sigma2=sqrt(W(2,2));
    deltaChi1=dm[0]; //V(1,2)=V(2,1) 
    deltaChi2=dm[1];
 
    return TrkErrCode(TrkErrCode::succeed);
    // FIXME: I wish I could tell poca to NOT iterate
    //        and ONLY compute the distance & derivatives...
    //  TrkDifPoca poca(*_trkTraj,fltLen(),*hitTraj(), hitLen(),_tolerance);
    /*
       if (poca.status().failure()) {
       return TrkErrCode(TrkErrCode::fail);
       }
       if (derivs.num_row() != 0) {
       poca.fetchDerivs(derivs);
       } else {
       derivs = poca.derivs();
       }
       double sigInv = 1. / hitRms();
       deltaChi = _resid * sigInv; // NOTE: use _INTERNAL_ residual
       derivs *= sigInv;
       */
    //return TrkErrCode(TrkErrCode::fail);//yzhang FIXME
}
 
TrkErrCode
//CgemHitOnTrack::getFitStuff(double &deltaChi)  const
CgemHitOnTrack::getFitStuff(double & sigma1, double & sigma2, double &deltaChi1, double &deltaChi2)  const
{
    assert (_trkTraj == &(getParentRep()->traj()));
    double vec_phi = _dHit->getrecphi() ;
    double vec_z = _dHit->getRecZ() ;
    double flt=0;
    HepVector trkpar = getParentRep()->traj().localTrajectory(0.,flt)->parameters()->parameter(); //get para??????
    IMdcUtilitySvc*   m_mdcUtilitySvc;
    Gaudi::svcLocator()->service("MdcUtilitySvc",m_mdcUtilitySvc);
    if (m_mdcUtilitySvc == NULL) cout<< " Can not load MdcUtilitySvc"<<endl;
    HepVector trkparbes = m_mdcUtilitySvc->patPar2BesPar(trkpar);
    trkparbes[2]=-trkparbes[2]; // modification for the reversed sign in TrackUtil/Helix
 
    HepPoint3D pivot(0., 0., 0.);
    Helix cgem_helix_0(pivot, trkparbes);
    HepSymMatrix V(2,0), W(2,0) ;
    HepVector dm(2);
    int layer = _dHit->getlayerid() ;
    double rl, R, x_reso, v_reso, a_stero;
    if(layer==0) { rl=79.838/10. ; R=87.5026/10.; x_reso=0.1372/10.; v_reso=0.1273/10.;a_stero=45.94*3.1415926/180;}
    else if (layer==1) { rl=125.104/10. ; R=132.7686/10.; x_reso=0.1476/10.;v_reso=0.1326/10.;a_stero=31.10*3.1415926/180;}
    else { rl=167.604/10.; R=175.2686/10.; x_reso=0.1412/10.;v_reso=0.1378/10.;a_stero=32.99*3.1415926/180;} //cm 
    double d_phi = cgem_helix_0.IntersectCylinder(rl); //const can only call const 
    HepPoint3D pos = cgem_helix_0.x(d_phi);
 
    double x = pos.x();
    double y = pos.y();
    double z = pos.z();
    double r2 = x*x+y*y;
    double tmp_phi = atan2(y, x); //[-pi,pi), return azimuthal angle
    double dphi = vec_phi - tmp_phi;
    if(dphi < -M_PI) dphi += 2*M_PI;
    if(dphi >= M_PI) dphi -= 2*M_PI;
 
    dm[0]=dphi*rl*sin(vec_phi);
    dm[1]=-(vec_z-z);
    V(1,1)=rl*rl/R/R*sin(vec_phi)*sin(vec_phi)*x_reso*x_reso;
    V(2,1)=rl/R*sin(vec_phi)/tan(a_stero)*x_reso*x_reso;
    V(2,2)=v_reso*v_reso/sin(a_stero)/sin(a_stero)+x_reso*x_reso/tan(a_stero)/tan(a_stero);
    W=V;
    int ifail=99;
    V.invert(ifail);
    if (ifail) {
        return TrkErrCode(TrkErrCode::fail);
    }
 
    sigma1=sqrt(W(1,1));
    sigma2=sqrt(W(2,2));
    deltaChi1=dm[0]; //V(1,2)=V(2,1) 
    deltaChi2=dm[1];
 
    return TrkErrCode(TrkErrCode::succeed);
}
 
TrkErrCode CgemHitOnTrack::getFitStuff(HepVector &Xderivs, double &XdeltaChi, HepVector &Vderivs, double &VdeltaChi) const
{
    //---get infomation from cluster
    int layer = _dHit->getlayerid();
    int sheet = _dHit->getsheetid();
    double phi_cluster = _dHit->getrecphi();
    double z_cluster = _dHit->getRecZ()/10;
    double Q = _dHit->getenergydeposit();
    //cout<<"layer:"<<layer<<endl;
 
    //---get geometry parameters from CgemGeomSvc
    CgemGeoReadoutPlane* readoutPlane = _cgemGeomSvc->getReadoutPlane(layer,sheet);
    double r_X = readoutPlane->getRX()/10;
    double r_V = readoutPlane->getRV()/10;
    double r_M = readoutPlane->getMidRAtGap()/10;
    double a_stero = readoutPlane->getStereoAngle();
    
    //---calculate the intersection of cgem layer and helix track
    double flt=0;
    HepVector trkpar = getParentRep()->traj().localTrajectory(0.,flt)->parameters()->parameter(); //get para??????
    //cout<<"barbar:"<<trkpar[0]<<"  "<<trkpar[1]<<"  "<<trkpar[2]<<"  "<<trkpar[3]<<"  "<<trkpar[4]<<endl;
    HepPoint3D pivot(0,0,0);
    trkpar[0] = -trkpar[0];// BarBar d0 -> BESIII dr
    trkpar[1] = trkpar[1]-M_PI/2;// BarBar phi0  -> BESIII phi0
    trkpar[2] = 333.567*trkpar[2];// BarBar  omega -> BESIII kappa
    trkpar[2] = -trkpar[2]; // modification for the reversed sign in TrackUtil/Helix
    trkpar[3] = trkpar[3];
    trkpar[4] = trkpar[4];
    while(trkpar[1] < 0) trkpar[1] += 2*M_PI;
    while(trkpar[1] >= 2*M_PI) trkpar[1] -= 2*M_PI;
    //cout<<"besIII:"<<trkpar[0]<<"  "<<trkpar[1]<<"  "<<trkpar[2]<<"  "<<trkpar[3]<<"  "<<trkpar[4]<<endl;
 
    Helix cgem_helix_0(pivot, trkpar);
    double dPhi = cgem_helix_0.IntersectCylinder(r_M);
    HepPoint3D pos = cgem_helix_0.x(dPhi);
    double x_cross = pos.x();
    double y_cross = pos.y();
    double z_cross = pos.z();
    double phi_cross = pos.phi(); //[-pi,pi)
    while( phi_cross < 0) phi_cross += 2*M_PI;
    while( phi_cross >= 2*M_PI) phi_cross -= 2*M_PI;
    double r2 = pos.perp2();
    //double dphi = phi_cluster - phi_cross;
    double dphi = phi_cross - phi_cluster;
    while(dphi < -M_PI) dphi += 2*M_PI;
    while(dphi >= M_PI) dphi -= 2*M_PI;
    //cout<<"radius:"<<cgem_helix_0.radius()<<endl;
    //cout<<"distance:"<<cgem_helix_0.center().perp()<<endl;
    
    //---get sigma of X and V strip from CgemCalibFunSvc
    int iView(0), mode(1);//FIXME
    double T(100);//FIXME
    double Phi_tangent = cgem_helix_0.direction(dPhi).phi();
    double Phi_position = phi_cross;//FIXME
    double delta_phi=Phi_tangent - Phi_position;
    while(delta_phi>M_PI)  delta_phi-=CLHEP::twopi;
    while(delta_phi<-M_PI) delta_phi+=CLHEP::twopi;
    double sigma_X = _cgemCalibFunSvc->getSigma(layer,0,mode,delta_phi,Q,T)/10;//  mm->cm
    double sigma_V = _cgemCalibFunSvc->getSigma(layer,1,mode,delta_phi,Q,T)/10;//  mm->cm
    
    //---calculate partial derivative of X and V w.r.t x,y,z
    HepMatrix J_dmdx(2,3,0), J_dxda(3,5,0), J(2,5,0);
    J_dmdx(1,1)=-r_X*y_cross/r2;
    J_dmdx(1,2)= r_X*x_cross/r2;
    J_dmdx(1,3)= 0.0;
    J_dmdx(2,1)=-r_V*cos(a_stero)*y_cross/r2;
    J_dmdx(2,2)= r_V*cos(a_stero)*x_cross/r2;
    J_dmdx(2,3)= sin(a_stero);
    
    //---calculate partial derivative of x,y,z w.r.t helix parameters:dr, phi0, kappa, dz, tanl
    double alpha = cgem_helix_0.alpha();
    //cout<<"alpha:"<<alpha<<endl;
    alpha = -alpha; // modification for the reversed sign in TrackUtil/Helix
    double dr = cgem_helix_0.dr();
    double phi0 = cgem_helix_0.phi0();
    double kappa = cgem_helix_0.kappa();
    kappa = -kappa; // modification for the reversed sign in TrackUtil/Helix
    double dz = cgem_helix_0.dz();
    double tanl = cgem_helix_0.tanl();
    
    double cosPhi=cos(dPhi);
    double sinPhi=sin(dPhi);
 
    double dPhiDdr = -(kappa/alpha-cosPhi/(dr+alpha/kappa))/sinPhi;
    //double dPhiDdr2= -kappa*kappa*(2.0*alpha*dr+kappa*(dr*dr+r2))/(2*alpha*(alpha+dr*kappa)*(alpha+dr*kappa)*sinPhi);
    //cout<<"dPhiDdr = "<<dPhiDdr<<endl;
    //cout<<"dPhiDdr2= "<<dPhiDdr2<<endl;
    double dPhiDkappa = -kappa*(dr*dr-r2)*(2*alpha+dr*kappa)/(2*alpha*(alpha+dr*kappa)*(alpha+dr*kappa)*sinPhi);
 
    double dxDdr = cos(phi0)+alpha/kappa*sin(phi0+dPhi)*dPhiDdr;
    double dyDdr = sin(phi0)-alpha/kappa*cos(phi0+dPhi)*dPhiDdr;
    double dxDphi0 = -dr*sin(phi0)+alpha/kappa*(-sin(phi0)+sin(phi0+dPhi));
    double dyDphi0 = -dr*cos(phi0)+alpha/kappa*( cos(phi0)-cos(phi0+dPhi));
    double dxDkappa = -alpha/(kappa*kappa)*(cos(phi0)-cos(phi0+dPhi))+alpha/kappa*sin(phi0+dPhi)*dPhiDkappa;
    double dyDkappa = -alpha/(kappa*kappa)*(sin(phi0)-sin(phi0+dPhi))-alpha/kappa*cos(phi0+dPhi)*dPhiDkappa;
    double dzDdr = -alpha/kappa*tanl*dPhiDdr;
    double dzDkappa = alpha/(kappa*kappa)*tanl*dPhi-alpha/kappa*tanl*dPhiDkappa;
    double dzDtanl = -alpha/kappa*dPhi;
 
    //---return the partial residuals and derivatives of X and V w.r.t helix parameters:d0, phi0, omega, z0, tanl
    J_dxda(1,1) = -1*dxDdr;//BESIII dr -> BarBar d0 
    J_dxda(1,2) = dxDphi0;
    J_dxda(1,3) = -alpha*dxDkappa;// BESIII kappa -> BarBar omega 
    J_dxda(2,1) = -1*dyDdr;//BESIII dr -> BarBar d0
    J_dxda(2,2) = dyDphi0;
    J_dxda(2,3) = -alpha*dyDkappa;// BESIII kappa -> BarBar omega
    J_dxda(3,1) = -1*dzDdr;//BESIII dr -> BarBar d0
    J_dxda(3,3) = -alpha*dzDkappa;// BESIII kappa -> BarBar omega
    J_dxda(3,4) = 1.0;
    J_dxda(3,5) = dzDtanl;
    J=J_dmdx*J_dxda;
    int nPar = getParentRep()->traj().localTrajectory(0.,flt)->nPar();
    for(int i=1;i<=nPar;i++){
        Xderivs(i) = J(1,i)/sigma_X;
        Vderivs(i) = J(2,i)/sigma_V;
    }
    XdeltaChi = dphi*r_X/sigma_X;//delta X
    //VdeltaChi = ((z_cluster-z_cross)*sin(a_stero)+r_V*cos(a_stero)*dphi)/sigma_V;
    VdeltaChi = ((z_cross-z_cluster)*sin(a_stero)+r_V*cos(a_stero)*dphi)/sigma_V;
    //cout<<"X:"<<phi_cluster<<"  "<<phi_cross<<"  "<<sigma_V<<endl;
    //cout<<"V:"<<z_cluster<<"  "<<z_cross<<"  "<<sigma_V<<endl;
    if(nPar == 3){
        Vderivs = HepVector(nPar,0);
        VdeltaChi = 0;
    }
    return TrkErrCode(TrkErrCode::succeed);
}
 
TrkErrCode CgemHitOnTrack::getFitStuff(double &XdeltaChi, double &VdeltaChi) const
{
    //---get infomation from cluster
    int layer = _dHit->getlayerid();
    int sheet = _dHit->getsheetid();
    double phi_cluster = _dHit->getrecphi();
    double z_cluster = _dHit->getRecZ()/10;
    double Q = _dHit->getenergydeposit();
 
    //---get geometry parameters from CgemGeomSvc
    CgemGeoReadoutPlane* readoutPlane = _cgemGeomSvc->getReadoutPlane(layer,sheet);
    double r_X = readoutPlane->getRX()/10;
    double r_V = readoutPlane->getRV()/10;
    double r_M = readoutPlane->getMidRAtGap()/10;
    double a_stero = readoutPlane->getStereoAngle();
 
    //---calculate the intersection of cgem layer and helix track
    double flt=0;
    HepVector trkpar = getParentRep()->traj().localTrajectory(0.,flt)->parameters()->parameter(); //get para??????
    HepPoint3D pivot(0,0,0);
    //cout<<"barbar parameters:"<<trkpar[0]<<"  "<<trkpar[1]<<"  "<<trkpar[2]<<"  "<<trkpar[3]<<"  "<<trkpar[4]<<endl;
    trkpar[0] = -trkpar[0];// BarBar d0 -> BESIII dr
    trkpar[1] = trkpar[1]-M_PI/2;// BarBar phi0  -> BESIII phi0
    trkpar[2] = 333.567*trkpar[2];// BarBar  omega -> BESIII kappa
    trkpar[2] = -trkpar[2]; // modification for the reversed sign in TrackUtil/Helix
    trkpar[3] = trkpar[3];
    trkpar[4] = trkpar[4];
    while(trkpar[1] < 0) trkpar[1] += 2*M_PI;
    while(trkpar[1] >= 2*M_PI) trkpar[1] -= 2*M_PI;
    //cout<<"besIII parameters:"<<trkpar[0]<<"  "<<trkpar[1]<<"  "<<trkpar[2]<<"  "<<trkpar[3]<<"  "<<trkpar[4]<<endl;
    Helix cgem_helix_0(pivot, trkpar);
    double dPhi = cgem_helix_0.IntersectCylinder(r_M);
    HepPoint3D pos = cgem_helix_0.x(dPhi);
    double x_cross = pos.x();
    double y_cross = pos.y();
    double z_cross = pos.z();
    double phi_cross = pos.phi(); //[-pi,pi)
    while( phi_cross < 0) phi_cross += 2*M_PI;
    while( phi_cross >= 2*M_PI) phi_cross -= 2*M_PI;
    double r2 = pos.perp2();
    //double dphi = phi_cluster - phi_cross;
    double dphi = phi_cross - phi_cluster;
    while(dphi < -M_PI) dphi += 2*M_PI;
    while(dphi >= M_PI) dphi -= 2*M_PI;
 
    //---get sigma of X and V strip from CgemCalibFunSvc
    int iView(0), mode(1);//FIXME
    double T(100);//FIXME
    double Phi_tangent = cgem_helix_0.direction(dPhi).phi();
    double Phi_position = phi_cross;//FIXME
    double delta_phi=Phi_tangent - Phi_position;
    while(delta_phi>M_PI)  delta_phi-=CLHEP::twopi;
    while(delta_phi<-M_PI) delta_phi+=CLHEP::twopi;
    double sigma_X = _cgemCalibFunSvc->getSigma(layer,0,mode,delta_phi,Q,T)/10;//  mm->cm
    double sigma_V = _cgemCalibFunSvc->getSigma(layer,1,mode,delta_phi,Q,T)/10;//  mm->cm
 
    //---return the partial derivatives and residuals of X and V
    int nPar = getParentRep()->traj().localTrajectory(0.,flt)->nPar();
    XdeltaChi = dphi*r_X/sigma_X;//delta X
    //VdeltaChi = ((z_cluster-z_cross)*sin(a_stero)+r_V*cos(a_stero)*dphi)/sigma_V;
    VdeltaChi = ((z_cross-z_cluster)*sin(a_stero)+r_V*cos(a_stero)*dphi)/sigma_V;
    //cout<<"X:"<<phi_cluster<<"  "<<phi_cross<<"  "<<sigma_V<<endl;
    //cout<<"V:"<<z_cluster<<"  "<<z_cross<<"  "<<sigma_V<<endl;
    if(nPar == 3){
        VdeltaChi = 0;
    }
    return TrkErrCode(TrkErrCode::succeed);
}
//cout<<__FILE__<<"   "<<__LINE__<<endl;