MdcTrack.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: MdcTrack.cxx,v 1.3 2017/12/27 06:29:02 huangzhen Exp $
//
// Description:
//     
//
// Environment:
//      Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author(s): 
//  Steve Schaffner
//  Zhang Yao([email protected])  Migrate to BESIII
//
//------------------------------------------------------------------------
#include "MdcTrkRecon/MdcTrack.h"
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include <map>
#include "MdcGeom/MdcLayer.h"
#include "MdcData/MdcHitOnTrack.h"
#include "MdcGeom/BesAngle.h"
#include "MdcGeom/Constants.h"
#include "CLHEP/Alist/AList.h"
#include "TrkBase/TrkRep.h"
#include "TrkFitter/TrkHelixRep.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrkBase/TrkRecoTrk.h"
#include "TrkFitter/TrkCircleMaker.h"
#include "TrkBase/TrkContext.h"
 
#include "MdcData/MdcRecoHitOnTrack.h"
#include "BField/BField.h"
#include "MdcData/MdcHit.h"
#include "MdcData/MdcHitOnTrack.h"
#include "MdcRawEvent/MdcDigi.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
//#include "MdcData/CgemHitOnTrack.h"
#include "CgemData/CgemHitOnTrack.h"
 
using namespace std;
//************************************************************************/
MdcTrack::MdcTrack(TrkRecoTrk* aTrack){
//************************************************************************/
  _theTrack = aTrack;
  _firstLayer = _lastLayer = 0;
  _haveCurled = 0;
}
 
//************************************************************************/
MdcTrack::MdcTrack(int nsupers, const TrkExchangePar& par, double chisq, 
             TrkContext& context, double trackT0) {
//************************************************************************/
    TrkCircleMaker maker;
    _theTrack = maker.makeTrack(par, chisq, context, trackT0);
    _firstLayer = _lastLayer = 0;
    _haveCurled = 0;
}
 
//************************************************************************/
MdcTrack::~MdcTrack() {
//************************************************************************/
  delete _theTrack;
}
//************************************************************************/
int MdcTrack::projectToR(double radius, BesAngle &phiIntersect, int lCurl) 
  const {
//************************************************************************/
  // note that these are currently circle-only routines
 
  //return -1 if no intersection
  
  const TrkFit* tkFit = track().fitResult();
  if (tkFit == 0) return -1;
  TrkExchangePar par = tkFit->helix(0.0);
  double d0 = par.d0();
  double omega = par.omega();
  double phi0 = par.phi0();
  double r2d2 = radius * radius - d0 * d0;
  if (r2d2 < 0) return -1;
  double rinv = 1./radius;
  double k2dinv = 1./(1 + omega * d0);
  if (!lCurl) {
    double arg = d0 * rinv + 0.5*omega * rinv * r2d2 * k2dinv;
    if (fabs(arg) > 1.0) return -1;
    phiIntersect.setRad( phi0 + asin(arg) );
  }
  else {
    double arg = -d0 * rinv - 0.5*omega * rinv * r2d2 * k2dinv;
    if (fabs(arg) > 1.0) return -1;
    phiIntersect.setRad( phi0 + Constants::pi + asin(arg) );
  }
   
  return 0;
}
 
//************************************************************************/
int MdcTrack::projectToR(double radius, BesAngle &phiIntersect, 
                double &arcLength, int lCurl) const {
//************************************************************************/
 
  //return -1 if no intersection
  // Only valid for projecting track outwards from point of closest approach
  // Returns arclength from POCA
 
  const TrkFit* tkFit = track().fitResult();
  if (tkFit == 0) return -1;
  TrkExchangePar par = tkFit->helix(0.0);
  double d0 = par.d0();
  double omega = par.omega();
  double phi0 = par.phi0();
  double r2d2 = radius * radius - d0 * d0;
  if (r2d2 < 0) return -1;
  double rinv = 1./radius;
  double k2dinv = 1./(1 + omega * d0);
  if (k2dinv < 0.0) return -1;
  double arg;
 
  if (!lCurl) {
    arg = d0 * rinv + 0.5*omega * rinv * r2d2 * k2dinv;
    if (fabs(arg) > 1.0) return -1;
    phiIntersect.setRad( phi0 + asin(arg) );
  }
  else {
    arg = -d0 * rinv - 0.5*omega * rinv * r2d2 * k2dinv;
    if (fabs(arg) > 1.0) return -1;
    phiIntersect.setRad( phi0 + Constants::pi + asin(arg));
  }
 
  arg = 0.5*omega * sqrt(r2d2 * k2dinv);
  arcLength = 2. * asin(arg) / omega;
 
  return 0;
}
 
//----------------------------------------------------------------
bool 
MdcTrack::operator==(const MdcTrack& tk) const {
//----------------------------------------------------------------
 
  return (track() == tk.track());
}
 
bool sortCluster(const RecCgemCluster* clusterA , const RecCgemCluster* clusterB){
    return clusterA->getlayerid()<clusterB->getlayerid();
}
 
//yzhang for store to TDS
void
MdcTrack::storeTrack(int trackId, RecMdcTrackCol* trackList, RecMdcHitCol* hitList,int tkStat){
    //check the fit succeeded
   // if (status().failure()) { return; }
    
    //get the results of the fit to this track
    const TrkFit* fitresult = track().fitResult();
    //check the fit worked
    if (fitresult == 0) return;
    
    //new a Rec. Track MdcTrack 
    RecMdcTrack* recMdcTrack = new RecMdcTrack();
    const TrkHitList* aList = track().hits();
    //some track info
    const BField& theField = track().bField();
    double Bz = theField.bFieldZ();
    //Fit info
    int hitId = 0;
    int nHits=0;
    int nSt=0;
    //int nAct=0; int nSeg=0; 
    //int maxlay = 0; int minlay = 43; int seg[11];/* 11 slayer */ int segme;
    //****************************
    //* active flag open this
    //****************************
    //* if (allHit>0){
    //*   nHits = aList->nHit();//add inActive hit also
    //* }else{
    //*   nHits = fitresult->nMdc();// = nActive()
    //* }
    //****************************
    //for 5.1.0 use all hits
    nHits = aList->nHit();
    //****************************
    // use active only
    // nHits = fitresult->nMdc();// = nActive()
    //****************************
 
    int    q        = fitresult->charge();
    double chisq    = fitresult->chisq();
    int    nDof     = fitresult->nDof();
    //track parameters at closest approach to beamline
    double fltLenPoca = 0.0;
    TrkExchangePar helix = fitresult->helix(fltLenPoca);
    //std::cout<< __FILE__ << " phi0  " << helix.phi0() << " omega  " <<helix.omega()<<std::endl;
    double phi0 = helix.phi0();
    double tanDip = helix.tanDip();
    //double z0 = helix.z0();
    double d0 = helix.d0();
    //momenta and position
    //CLHEP::Hep3Vector mom = fitresult->momentum(fltLenPoca);
    HepPoint3D poca = fitresult->position(fltLenPoca);
 
    double helixPar[5];
    //dro =-d0
    helixPar[0] = -d0;
    //phi0 = phi0 - pi/2   [0,2pi)
    double tphi0 = phi0 - Constants::pi/2.;
    if (tphi0 < 0. ) tphi0 += Constants::pi * 2.;
    helixPar[1] = tphi0;
    //kappa = q/pxy
    double pxy = fitresult->pt();
    if (pxy == 0.)  helixPar[2] = 9999.;
    else helixPar[2] = q/fabs(pxy);
    if(pxy>9999.) helixPar[2] = 0.00001;
    //dz = z0
    helixPar[3] = helix.z0();
    //tanl 
    helixPar[4] = tanDip;
    //error
    //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X)
    HepSymMatrix mS(helix.params().num_row(),0);
    mS[0][0]=-1.;//dr0=-d0
    mS[1][1]=1.;
    mS[2][2]=-333.567/Bz;//pxy -> cpa
    mS[3][3]=1.;
    mS[4][4]=1.;
    HepSymMatrix mVy = helix.covariance().similarity(mS);
    double errorMat[15];
    int k = 0;
    for (int ie = 0 ; ie < 5 ; ie ++){
      for (int je = ie ; je < 5 ; je ++) {
    errorMat[k] = mVy[ie][je];
    k++;
      }
    }
    double p,px,py,pz;
    px = pxy * (-sin(helixPar[1]));
    py = pxy * cos(helixPar[1]);
    pz = pxy * helixPar[4];
    p = sqrt(pxy*pxy + pz*pz);
    //theta, phi
    double theta = acos(pz/p);
    double phi = atan2(py,px);
    recMdcTrack->setTrackId(trackId);
    recMdcTrack->setHelix(helixPar); 
    recMdcTrack->setCharge(q);
    recMdcTrack->setPxy(pxy);
    recMdcTrack->setPx(px);
    recMdcTrack->setPy(py);
    recMdcTrack->setPz(pz);
    recMdcTrack->setP(p);
    recMdcTrack->setTheta(theta);
    recMdcTrack->setPhi(phi);
    recMdcTrack->setPoca(poca);
    recMdcTrack->setX(poca.x());//poca
    recMdcTrack->setY(poca.y());
    recMdcTrack->setZ(poca.z());
    recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
    HepPoint3D pivot(0.,0.,0.);
    recMdcTrack->setPivot(pivot);
    recMdcTrack->setVX0(0.);//pivot 
    recMdcTrack->setVY0(0.);
    recMdcTrack->setVZ0(0.);
    recMdcTrack->setError(errorMat);
    recMdcTrack->setError(mVy);
    recMdcTrack->setChi2(chisq);
    recMdcTrack->setStat(tkStat);
    recMdcTrack->setNhits(nHits);
    //-----------hit list-------------
    HitRefVec  hitRefVec;
    ClusterRefVec clusterRefVec;
    vector<int> vecHits;
    map<int,int> clusterFitStat;
    //for fiterm
    int maxLayerId = -1;
    int minLayerId = 43;
    double fiTerm = 999.;
    const MdcRecoHitOnTrack* fiTermHot = NULL;
    TrkHitList::hot_iterator hot = aList->begin();
    for (;hot!=aList->end();hot++){
      //cout<<__FILE__ <<"  "<<__LINE__ <<"  "<< typeid(*hot).name()<<endl;
      if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
        const MdcRecoHitOnTrack* recoHot;             
        recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
        //if (!recoHot->mdcHit()) return;
        RecMdcHit* recMdcHit = new RecMdcHit;
        //id
        recMdcHit->setId(hitId);
        //---------BESIII----------
        //phiWire - phiHit <0 flagLR=0 left driftleft<0 ;
        //phiWire - phiHit >0 flagLR=1 right driftright>0;
        //flag = 2 ambig;
        //-----Babar wireAmbig----
        //-1 -> right, 0 -> don't know, +1 -> left
        //+1 phiWire-phiHit<0 Left,
        //-1 phiWire-phiHit>0 Right, 
        //0  don't know   
        //ambig() w.r.t track direction
        //wireAmbig() w.r.t. wire location 
        double hotWireAmbig = recoHot->wireAmbig();
        double driftDist = fabs(recoHot->drift());
        double sigma = recoHot->hitRms();       
        double doca = fabs(recoHot->dcaToWire());
        int flagLR=2;
        if ( hotWireAmbig == 1){
          flagLR = 0; //left driftDist <0
          doca *= -1.;//2012-07-19 
        }else if( hotWireAmbig == -1){
          flagLR = 1; //right driftDist >0
        }else if( hotWireAmbig == 0){
          flagLR = 2; //don't know 
        }
        recMdcHit->setFlagLR(flagLR);
        recMdcHit->setDriftDistLeft((-1)*driftDist);
        recMdcHit->setErrDriftDistLeft(sigma);
        recMdcHit->setDriftDistRight(driftDist);
        recMdcHit->setErrDriftDistRight(sigma);
        //Mdc Id
        Identifier mdcId = recoHot->mdcHit()->digi()->identify();
        recMdcHit->setMdcId(mdcId);
        //corrected ADC
 
        //contribution to chisquare
        //fill chisq
        double res=999.,rese=999.;
        if (recoHot->resid(res,rese,false)){
        }else{}
        double deltaChi=0;
        recoHot->getFitStuff(deltaChi);//yzhang open 2010-09-20 
        recMdcHit->setChisqAdd( deltaChi * deltaChi );
        //set tracks containing this hit
        recMdcHit->setTrkId(trackId);
        //doca
        recMdcHit->setDoca(doca);//doca sign left<0
        //entrance angle
        
        recMdcHit->setEntra(recoHot->entranceAngle());
        //z of hit 
        HepPoint3D pos; Hep3Vector dir;
        double fltLen = recoHot->fltLen(); 
        recMdcHit->setFltLen(fltLen);
        recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
        recMdcHit->setZhit(pos.z());
        double tof = recoHot->getParentRep()->arrivalTime(fltLen);
        recMdcHit->setTdc(recoHot->mdcHit()->tdcIndex());
        recMdcHit->setAdc(recoHot->mdcHit()->adcIndex()); 
        //drift time      
        recMdcHit->setDriftT(recoHot->mdcHit()->driftTime(tof,pos.z()));
        //for fiterm
        int layerId = recoHot->mdcHit()->layernumber();
        int wireId = recoHot->mdcHit()->wirenumber();
        //std::cout<<" MdcTrack::store() ("<<layerId<<","<<wireId<<") fltLen "<<fltLen<<" wireAmbig  "<<hotWireAmbig<<"    y "<<pos.y()<<std::endl;
    //  <<recoHot->entranceAngle()<<std::endl;
        if (layerId >= maxLayerId){
          maxLayerId = layerId;
          fiTermHot = recoHot;
        } 
        if (layerId < minLayerId){
          minLayerId = layerId;
        }
        // status flag
        if (recoHot->isActive()) {
          recMdcHit->setStat(1);
        }else{
          recMdcHit->setStat(0);    
        }
        // for 5.1.0 use all hits
        if (recoHot->layer()->view()) {
          ++nSt;
        }
        hitList->push_back(recMdcHit);
        SmartRef<RecMdcHit> refHit(recMdcHit);
        hitRefVec.push_back(refHit);        
        vecHits.push_back(mdcId.get_value());
        ++hitId;          
      }else if(typeid(*hot)==typeid(CgemHitOnTrack)){
        const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));            
        const RecCgemCluster* recCgemCluster = recoHot->baseHit();
        int clusterid = recCgemCluster->getclusterid();
        int stat = recoHot->isActive();
        clusterFitStat[clusterid] = stat;
        //cout<<clusterid<<"  "<<stat<<endl;
        //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
        clusterRefVec.push_back(recCgemCluster);
      }
    }
    std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
    //fi terminal angle
    if (fiTermHot!=NULL){
      fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
    }
    recMdcTrack->setFiTerm(fiTerm);
    // number of  stereo hits contained
    recMdcTrack->setNdof(nDof);
    recMdcTrack->setNster(nSt);  
    //recMdcTrack->setFirstLayer(maxLayerId);  
    //recMdcTrack->setLastLayer(minLayerId);  
    recMdcTrack->setFirstLayer(minLayerId);  
    recMdcTrack->setLastLayer(maxLayerId);  
    recMdcTrack->setVecHits(hitRefVec);
    //recMdcTrack->setClusterFitStat(clusterFitStat);
    //recMdcTrack->setVecClusters(clusterRefVec);
    recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
    recMdcTrack->setNcluster(clusterRefVec.size());
    trackList->push_back(recMdcTrack);
 
    //cout<<__FILE__ <<"  "<<__LINE__ << endl;
    //cout<<"nDof "<<nDof<<endl;
    //cout<<"nSter "<<nSt<<endl;
    //cout<<"maxLayerId "<<maxLayerId<<endl;
    //cout<<"minLayerId "<<minLayerId<<endl;
    //cout<<"nHits "<<nHits<<endl;
    //cout<<"nHots "<<track().hots()->nHit()<<endl;
    //cout<<"nActs "<<track().hots()->nActive()<<endl;
    //cout<<"nCluster "<<clusterRefVec.size()<<endl;
    //cout<<"nActive "<<fitresult->nActive()<<endl;
    //cout<<endl;
 
    //TrkHotList* hotList = track().hots();
    //int nAct = 0;
    //cout<<hotList->end() - hotList->begin()<<endl;
    //for(TrkHotList::nc_hot_iterator i = hotList->begin();i!=hotList->end();++i){
        //cout<<i->layerNumber()<<"  "<<i->isActive()<<"  "<<i->whatView()<<endl;
    //}
    //hot = aList->begin();
    //for (;hot!=aList->end();hot++){
        //cout<<hot->layerNumber()<<"  "<<hot->isActive()<<"  "<<hot->whatView()<<endl;
    //}
       
    //cout<<__FILE__ <<"  "<<__LINE__ << endl;
    //for(int i=0;i<clusterRefVec.size();i++)
    //{
        //cout<<clusterRefVec[i]->getlayerid()<<endl;
    //}
    //hot = aList->begin();
    //for (;hot!=aList->end();hot++){
      //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
        //const MdcRecoHitOnTrack* recoHot;           
        //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
        //int layerId = recoHot->mdcHit()->layernumber();
        //int wireId = recoHot->mdcHit()->wirenumber();
        //cout<<"hit ("<<layerId<<","<<wireId<<")"<<endl;
      //}else if(typeid(*hot)==typeid(CgemHitOnTrack)){
        //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));          
        //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
        //cout<<"hit ("<<recCgemCluster->getlayerid()<<","<<recCgemCluster->getclusterid()<<")"<<endl;
        //clusterRefVec.push_back(recCgemCluster);
      //}
    //}
}//end of storeTrack