ExtMdcTrack.cxx

Go to the documentation of this file.

 
////////////////
// #define Helix Ext_Helix
// File: ExtMdcTrack.cxx
// Author: L.L.Wang
// Deccription: A class to get MDC reconstuction track,
//              and supply information to  track extrapolation.
//              It is a interface.
// History: 2005.7.4 created by L.L.Wang
//
 
#include "CLHEP/Matrix/Matrix.h"
// #include "CLHEP/config/CLHEP.h"  //For constant PI(3.14......).
 
/////////////////////////////change
#include "Helix.h"
// #include "TrackUtil/Helix.h"
/////////////////////////////////////////
 
// #include "MdcFastTrkAlg/helix/Helix.h"
 
#include "ExtMdcTrack.h"
#include <float.h>
 
using namespace std;
 
const double mass[5] = { 0.511, 105.658, 139.57, 493.677, 938.27203 };
 
ExtMdcTrack::ExtMdcTrack()
    : myMsgFlag( true )
    , myInputTrk( "Kal" )
    , myParID( 2 )
    , myHelixPar( NumHelixPar )
    , myMdcErr( NdimMdcErr, 0 )
    , myPhiTerm( 0. ) {}
 
ExtMdcTrack::~ExtMdcTrack() {}
 
bool ExtMdcTrack::SetMdcRecTrkCol( RecMdcTrackCol* aPointer ) {
  myMdcTrackCol   = aPointer;
  myMdcRecTrkIter = myMdcTrackCol->begin();
  myMdcRecTrkIter--;
  myInputTrk = "Mdc";
  return true;
}
 
bool ExtMdcTrack::SetMdcKalTrkCol( RecMdcKalTrackCol* aPointer ) {
  myMdcKalTrackCol = aPointer;
  myMdcKalTrkIter  = myMdcKalTrackCol->begin();
  myMdcKalTrkIter--;
  myInputTrk = "Kal";
  return true;
}
 
bool ExtMdcTrack::GetOneGoodTrk() {
  if ( myMsgFlag ) cout << "ExtMdcTrack::GetOneGoodTrk() begin" << endl;
  if ( myInputTrk == "Mdc" )
  {
    myMdcRecTrkIter++;
    if ( myMdcRecTrkIter != myMdcTrackCol->end() ) return true;
    else
    {
      if ( myMsgFlag ) cout << "No more track!" << endl;
      return false;
    }
  }
  else if ( myInputTrk == "Kal" ) // use KalTrack as input
  {
    myMdcKalTrkIter++;
    if ( myMdcKalTrkIter != myMdcKalTrackCol->end() ) return true;
    else
    {
      if ( myMsgFlag ) cout << "No more track!" << endl;
      return false;
    }
  }
 
  std::cerr << __FILE__ << ":" << __LINE__ << ": should not reach here" << std::endl;
  exit( 1 );
}
 
bool ExtMdcTrack::ReadTrk( int pid = 2 ) {
  if ( myMsgFlag ) cout << "ExtMdcTrack::ReadTrk() begin" << endl;
  myParID = pid;
  if ( myInputTrk == "Mdc" )
  {
    if ( myMdcRecTrkIter != myMdcTrackCol->end() )
    {
      RecMdcTrack* aMdcTrack = *myMdcRecTrkIter;
      //                if(!aMdcTrack->getStat())//If good track.Now status default 0
      {
        if ( myMsgFlag ) cout << "Read a good track!" << endl;
 
        // get MdcTrk data
        myTrackID  = aMdcTrack->trackId();
        myHelixPar = aMdcTrack->helix();
        myPivot    = aMdcTrack->getPivot();
        myMdcErr   = aMdcTrack->err();
        //              double aPhiTerm = aMdcTrack->getFiTerm();
        myPhiTerm   = aMdcTrack->getFiTerm();
        myTrkLength = GetTrackLength();
        double p    = ( GetMomentum() ).mag();
        double beta = p / sqrt( mass[myParID] * mass[myParID] + p * p );
        myTrkTof    = myTrkLength / ( beta * 299.792458 );
 
        /*              //Start caculation of myPhiTerm.(when pivot is
           (0,0,0).) double r = DBL_MAX; double rMdc = 81.0; double dro = myHelixPar(1);
                                        if(myHelixPar(3)!=0.0)
                                        r = 333.564095/myHelixPar(3);
                                        if( fabs(2.0*r+dro)>rMdc && fabs(dro)<rMdc )
                                        {
                                        double aValue =
           ((dro+r)*(dro+r)+r*r-rMdc*rMdc)/2.0/fabs(r)/fabs(dro+r); if(r<0.0) myPhiTerm =
           acos(aValue); else myPhiTerm = -1.0*acos(aValue);
                                        }
                                        else
                                        {
                                        if(myMsgFlag) cout<<"Get a track which can't go out of
           MDC!"<<endl; myMdcRecTrkIter++; continue;
                                        }
         */
        //              double
        // aTrackLength=fabs(r*myPhiTerm)*sqrt(myHelixPar(5)*myHelixPar(5)+1);
        if ( myMsgFlag )
        {
          cout << "MDC track:" << myTrackID << endl;
          cout << "Helix: " << myHelixPar;
          cout << "pivot: " << myPivot << "cm" << endl;
          cout << "Error: " << myMdcErr << endl;
 
          /*                Convert();
                                          cout<<"********after convert:"<<endl;
                                          cout<<"Helix: "<<myHelixPar;
                                          cout<<"pivot: "<<myPivot;
                                          cout<<"Error: "<<myMdcErr<<endl;
           */
          cout << "Pt: " << GetPt() << "GeV/c" << endl;
          cout << "PhiTerm: " << myPhiTerm << endl;
          //                    cout<<"recPhiTerm: "<<aPhiTerm<<endl;
          cout << "trackLength(mm): " << GetTrackLength() << endl;
          //                    cout<<"trackLengthTest(cm):
          //"<<aTrackLength<<endl;
        }
 
        return true;
      }
      //            else
      //            {
      //                myMdcRecTrkIter++;
      //                cout<<"Get a bad track!"<<endl;
      //            }
    }
    else
    {
      if ( myMsgFlag ) cout << "No more track!" << endl;
      return false;
    }
  }
  else if ( myInputTrk == "Kal" ) // use KalTrack as input
  {
    if ( myMdcKalTrkIter != myMdcKalTrackCol->end() )
    {
      RecMdcKalTrack* aMdcKalTrack = *myMdcKalTrkIter;
      {
        if ( myMsgFlag ) cout << "Get a good track!" << endl;
 
        // get MdcTrk data
        myTrackID  = aMdcKalTrack->getTrackId();
        myHelixPar = aMdcKalTrack->getLHelix( myParID );
        // myPivot = HepPoint3D(0,0,0);
        myPivot     = aMdcKalTrack->getLPivot( myParID );
        myMdcErr    = aMdcKalTrack->getLError( myParID );
        myPhiTerm   = aMdcKalTrack->getFiTerm( myParID );
        myPhiTerm   = 0.0;
        myTrkLength = aMdcKalTrack->getPathSM( myParID );
        myLPosition = aMdcKalTrack->getLPoint( myParID );
        myTrkTof    = aMdcKalTrack->getTof( myParID );
        // cout<<"myInitialTof = "<<myTrkTof<<endl;
 
        if ( myMsgFlag )
        {
          cout << "Kal track:" << myTrackID << endl;
          cout << "myParID:" << myParID << endl;
          cout << "Helix: " << myHelixPar;
          cout << "pivot: " << myPivot << "cm" << endl;
          cout << "Error: " << myMdcErr << endl;
          cout << "Pt: " << GetPt() << "GeV/c" << endl;
          cout << "PhiTerm: " << myPhiTerm << endl;
          //                    cout<<"trackLength(mm):
          //"<<GetTrackLength()<<endl;
          cout << "trackLength(cm): " << myTrkLength << endl;
          cout << "myTrkTof: " << myTrkTof << endl;
          cout << "Last point position: " << myLPosition << endl;
        }
        return true;
      }
    }
    else
    {
      if ( myMsgFlag ) cout << "No more track!" << endl;
      return false;
    }
  }
  std::cerr << __FILE__ << ":" << __LINE__ << ": should not reach here" << std::endl;
  exit( 1 );
}
 
void ExtMdcTrack::Convert() {
  myHelixPar( 1 ) *= 10.;
  myHelixPar( 3 ) *= 0.001;
  myHelixPar( 4 ) *= 10.;
  myPivot *= 10.;
  myMdcErr.fast( 1, 1 ) *= 100.;
  myMdcErr.fast( 2, 1 ) *= 10.;
  myMdcErr.fast( 3, 1 ) *= 0.01;
  myMdcErr.fast( 3, 2 ) *= 0.001;
  myMdcErr.fast( 3, 3 ) *= 0.000001;
  myMdcErr.fast( 4, 1 ) *= 100.;
  myMdcErr.fast( 4, 2 ) *= 10.;
  myMdcErr.fast( 4, 3 ) *= 0.01;
  myMdcErr.fast( 4, 4 ) *= 100.;
  myMdcErr.fast( 5, 1 ) *= 10.;
  myMdcErr.fast( 5, 3 ) *= 0.001;
  myMdcErr.fast( 5, 4 ) *= 10.;
}
 
RecMdcTrack* ExtMdcTrack::GetMdcRecTrkPtr() const { return ( *myMdcRecTrkIter ); }
 
const Hep3Vector ExtMdcTrack::GetPosition() const {
  //    if(myInputTrk=="Mdc") {
  Ext_Helix  aHelix( myPivot, myHelixPar, myMdcErr );
  Hep3Vector aPoint = aHelix.x( myPhiTerm );
  if ( myMsgFlag ) cout << "PhiTerm Position: " << aPoint << endl;
  aPoint *= 10.0;
  return ( aPoint );
  /*    }
          else if(myInputTrk=="Kal") {
                  return (myLPosition);
          }
  */
}
 
const Hep3Vector ExtMdcTrack::GetMomentum() const {
  Ext_Helix  aHelix( myPivot, myHelixPar, myMdcErr );
  Hep3Vector aMomentum = aHelix.momentum( myPhiTerm ); // In GeV/c units.
  if ( myMsgFlag ) cout << "PhiTerm Momentum: " << aMomentum << endl;
  aMomentum *= 1000.0; // In MeV/c uints.
  return ( aMomentum );
}
 
const HepSymMatrix ExtMdcTrack::GetErrorMatrix() const {
  Ext_Helix aHelix( myPivot, myHelixPar, myMdcErr );
  //    aHelix.pivot(aHelix.x(myPhiTerm));//Move the pivot to the PhiTerm position.
 
  Hep3Vector aPoint   = aHelix.x( myPhiTerm );
  double     phi0     = aHelix.phi0();
  double     kappaInv = 1.0 / aHelix.kappa();
  double     cosPhi0  = cos( phi0 );
  double     sinPhi0  = sin( phi0 );
  // 0 -> myPhiTerm, 2009-5-29, by wangll
  //    double px = aHelix.momentum(0).x();
  //    double py = aHelix.momentum(0).y();
  //    double pz = aHelix.momentum(0).z();
  double px = aHelix.momentum( myPhiTerm ).x();
  double py = aHelix.momentum( myPhiTerm ).y();
  double pz = aHelix.momentum( myPhiTerm ).z();
 
  // Calculate the Jacobian: d(xp)/d(helix).
  HepMatrix Helix2XpJcb( NdimExtErr, NdimMdcErr, 0 );
 
  Helix2XpJcb( 1, 1 ) = cosPhi0;
  Helix2XpJcb( 1, 2 ) = myPivot.y() - aPoint.y();
  Helix2XpJcb( 1, 3 ) = ( myPivot.x() - aPoint.x() + myHelixPar( 1 ) * cosPhi0 ) * kappaInv;
  Helix2XpJcb( 2, 1 ) = sinPhi0;
  Helix2XpJcb( 2, 2 ) = aPoint.x() - myPivot.x();
  Helix2XpJcb( 2, 2 ) = ( myPivot.y() - aPoint.y() + myHelixPar( 1 ) * sinPhi0 ) * kappaInv;
  Helix2XpJcb( 3, 3 ) = ( myPivot.z() - aPoint.z() + myHelixPar( 4 ) ) * kappaInv;
  Helix2XpJcb( 3, 4 ) = 1.0;
  Helix2XpJcb( 3, 5 ) = ( aPoint.z() - myPivot.z() - myHelixPar( 4 ) ) / myHelixPar( 5 );
  Helix2XpJcb( 4, 2 ) = -py;
  Helix2XpJcb( 4, 3 ) = -px * kappaInv;
  Helix2XpJcb( 5, 2 ) = px;
  Helix2XpJcb( 5, 3 ) = -py * kappaInv;
  Helix2XpJcb( 6, 3 ) = -pz * kappaInv;
  Helix2XpJcb( 6, 5 ) = fabs( kappaInv );
 
  // Calculate the Ext error matrix(6x6) in (x,p).
  //    if(myMsgFlag) cout<<"Error matrix at new pivot: "<<endl;
  //    aHelix.Ea();
  //    cout<<"haha"<<endl;
  HepSymMatrix aErrorMatrix = aHelix.Ea().similarity( Helix2XpJcb );
  if ( myMsgFlag ) cout << "PhiTerm Error matrix:" << aErrorMatrix << endl;
 
  // Uints Conversion.
  for ( int i = 1; i <= 6; i++ )
  {
    for ( int j = 1; j <= i; j++ )
    {
      if ( i <= 3 && j <= 3 ) aErrorMatrix.fast( i, j ) *= 100.0;
      if ( i >= 4 && j >= 4 ) aErrorMatrix.fast( i, j ) *= 1.0e+6;
      if ( i >= 4 && j <= 3 ) aErrorMatrix.fast( i, j ) *= 10000.0;
    }
  }
 
  return aErrorMatrix;
}
 
double ExtMdcTrack::GetTrackLength() const {
  if ( myInputTrk == "Mdc" )
  {
 
    /*  //old version
            Helix aHelix(myPivot,myHelixPar,myMdcErr);
 
            double phi0 = aHelix.phi0();
            HepPoint3D w = -(*this).GetParticleCharge() * aHelix.center();
            Hep3Vector center( w.x(),w.y(),w.z() );
            double centerX = center.x();
            double centerY = center.y();
            double phiCenter;
            if(fabs(centerX) > 1.0e-10)
            {
                    phiCenter = atan2(centerY,centerX);
                    if(phiCenter < 0.) phiCenter += 2.0*M_PI;
            }
            else
            {
                    phiCenter = (centerY>0)? M_PI_2:3.0*M_PI_2;
            }
            double dPhi = fabs(phi0 + myPhiTerm - phiCenter);
            if(dPhi >= 2.0*M_PI) dPhi -= 2.0*M_PI;
            double tanLambda = aHelix.tanl();
            double cosLambdaInv = sqrt(tanLambda*tanLambda + 1.0);
            return (10*fabs(aHelix.radius() * dPhi * cosLambdaInv));//10* for cm -> mm
    */
 
    Ext_Helix aHelix( myPivot, myHelixPar, myMdcErr );
    double    dPhi         = fabs( myPhiTerm );
    double    tanLambda    = aHelix.tanl();
    double    cosLambdaInv = sqrt( tanLambda * tanLambda + 1.0 );
    return ( 10 * fabs( aHelix.radius() * dPhi * cosLambdaInv ) ); // 10* for cm -> mm
  }
  else if ( myInputTrk == "Kal" )
  {
    return 10 * myTrkLength; // 10* for cm -> mm
  }
  std::cerr << __FILE__ << ":" << __LINE__ << ": should not reach here" << std::endl;
  exit( 1 );
}
 
double ExtMdcTrack::GetPt() const { return ( 1.0 / fabs( myHelixPar[2] ) ); }
 
double ExtMdcTrack::GetParticleCharge() const {
  return ( ( myHelixPar[2] > 0. ) ? 1.0 : -1.0 );
}