StMuDSTMaker/COMMON/StMuTrack.cxx

/***************************************************************************
 *
 * $Id: StMuTrack.cxx,v 1.48 2012/05/07 14:47:06 fisyak Exp $
 *
 * Author: Frank Laue, BNL, laue@bnl.gov
 ***************************************************************************/

#include "StMuTrack.h"
#include "StMuDebug.h"
#include "StMuException.hh"
#include "StEvent/StEventTypes.h"
#include "StEvent/StTrackGeometry.h"
#include "StEvent/StPrimaryVertex.h"
#include "StEvent/StDcaGeometry.h"
#include "StEvent/StBTofHit.h"
#include "StEvent/StBTofPidTraits.h"
#include "StarClassLibrary/SystemOfUnits.h"
#include "StEvent/StTpcDedxPidAlgorithm.h"
#include "StarClassLibrary/StThreeVectorD.hh"
#include "StarClassLibrary/StPhysicalHelixD.hh"
#include "StarClassLibrary/StParticleTypes.hh"
#include "StEventUtilities/StuProbabilityPidAlgorithm.h"
#include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
#include "StEmcUtil/projection/StEmcPosition.h"
#include "StEmcUtil/geometry/StEmcGeom.h"
#include "THelixTrack.h"
#include "TMath.h"
#include "TString.h"
namespace {
        const StThreeVectorF gDummy(-999,-999,-999);
}

StuProbabilityPidAlgorithm* StMuTrack::mProbabilityPidAlgorithm=0;
double StMuTrack::mProbabilityPidCentrality=0;
 
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
StMuTrack::StMuTrack(const StEvent* event, const StTrack* track, const StVertex *vertex, int index2Global, int index2RichSpectra, bool l3, TObjArray *vtxList) : 
  mId(0), mType(0), mFlag(0), mFlagExtension(0), mIndex2Global(index2Global), mIndex2RichSpectra(index2RichSpectra), mNHits(0), mNHitsPoss(0), mNHitsDedx(0),mNHitsFit(0), 
  mPidProbElectron(0), mPidProbPion(0),mPidProbKaon(0),mPidProbProton(0), 
  /* mNSigmaElectron(__NOVALUE__), mNSigmaPion(__NOVALUE__), mNSigmaKaon(__NOVALUE__), mNSigmaProton(__NOVALUE__) ,*/ 
  mdEdx(0.), mPt(0.), mEta(0.), mPhi(0.), mIndex2Cov(-1), mIdTruth(0), mQuality(0), mIdParentVx(0)  {

  const StGlobalTrack* globalTrack = dynamic_cast<const StGlobalTrack*>(track->node()->track(global));

  mId = track->key();
  mType = track->type();
  mFlag = track->flag();
  mFlagExtension = track->flagExtension();
  mTopologyMap = track->topologyMap();
  mIdTruth = track->idTruth();
  mQuality = track->qaTruth();
  mIdParentVx = track->idParentVx();
  // while getting the bestGuess, the pidAlgorithm (StTpcDedxPidAlgorithm) is set up.
  // pointers to track and pidTraits are set 
  // So, even though BestGuess will generate a "variable not used" warning, DO NOT DELETE THE NEXT LINES
  static StTpcDedxPidAlgorithm PidAlgorithm;
  static StElectron* Electron = StElectron::instance();
  static StPionPlus* Pion = StPionPlus::instance();
  static StKaonPlus* Kaon = StKaonPlus::instance();
  static StProton* Proton = StProton::instance();
  const StParticleDefinition* pd = track->pidTraits(PidAlgorithm);
  if (pd) {
    mNSigmaElectron = pack2Int( fabsMin(PidAlgorithm.numberOfSigma(Electron),__SIGMA_SCALE__), __SIGMA_SCALE__ );
    mNSigmaPion =     pack2Int( fabsMin(PidAlgorithm.numberOfSigma(Pion),__SIGMA_SCALE__),     __SIGMA_SCALE__ );
    mNSigmaKaon =     pack2Int( fabsMin(PidAlgorithm.numberOfSigma(Kaon),__SIGMA_SCALE__),     __SIGMA_SCALE__ );
    mNSigmaProton =   pack2Int( fabsMin(PidAlgorithm.numberOfSigma(Proton),__SIGMA_SCALE__),   __SIGMA_SCALE__ );
  } else {
    // BEWARE : Data members initialized to __NOVALUE__ was an hazardous
    // bunsiness as both NOVALUE and SIGMA_SCALE would have to be changed 
    // simultaneously (error prone). The above is equivalent while safer ...
    mNSigmaElectron = int(__NOVALUE__*__SIGMA_SCALE__); // pack2Int( fabsMin(__NOVALUE__*1.0,__SIGMA_SCALE__), __SIGMA_SCALE__ );
    mNSigmaPion =     int(__NOVALUE__*__SIGMA_SCALE__); // pack2Int( fabsMin(__NOVALUE__*1.0,__SIGMA_SCALE__), __SIGMA_SCALE__ );
    mNSigmaKaon =     int(__NOVALUE__*__SIGMA_SCALE__); // pack2Int( fabsMin(__NOVALUE__*1.0,__SIGMA_SCALE__), __SIGMA_SCALE__ );
    mNSigmaProton =   int(__NOVALUE__*__SIGMA_SCALE__); // pack2Int( fabsMin(__NOVALUE__*1.0,__SIGMA_SCALE__), __SIGMA_SCALE__ );
  }

  // if we have pid traits
  if ( PidAlgorithm.traits() ) {
    mdEdx = PidAlgorithm.traits()->mean();
    mNHitsDedx = PidAlgorithm.traits()->numberOfPoints();
  }

  // the following just point to particle definitions in StEvent
  if ( track->detectorInfo() ) {
    mFirstPoint = track->detectorInfo()->firstPoint();
    mLastPoint = track->detectorInfo()->lastPoint();
    mNHits = track->detectorInfo()->numberOfPoints();
  } 

  mNHitsPoss = track->numberOfPossiblePoints();

  int tpc_hits;
  mNHitsPossTpc=0;
  if (track->numberOfPossiblePoints(kTpcId)==track->numberOfPossiblePoints(kFtpcEastId)) {
    // backward compatibility mode, figure out which TPC points are in
    if (track->topologyMap().hasHitInDetector(kTpcId))
      mNHitsPossTpc=track->numberOfPossiblePoints(kTpcId);
    else if (track->topologyMap().hasHitInDetector(kFtpcEastId))
      mNHitsPossTpc=(1 << 6) + track->numberOfPossiblePoints(kFtpcEastId);
    else if (track->topologyMap().hasHitInDetector(kFtpcWestId))
      mNHitsPossTpc=(2 << 6) + track->numberOfPossiblePoints(kFtpcWestId);
  } else {
    // new Ittf chain distinguishes the three TPCs properly
    if ( (tpc_hits=track->numberOfPossiblePoints(kTpcId)) )
      mNHitsPossTpc=tpc_hits; 
    else if ( (tpc_hits=track->numberOfPossiblePoints(kFtpcEastId)) )
      mNHitsPossTpc=(1 << 6) + tpc_hits;
    else if ( (tpc_hits=track->numberOfPossiblePoints(kFtpcWestId)) ) 
      mNHitsPossTpc=(2 << 6) + tpc_hits;
  }

  mNHitsPossInner=track->numberOfPossiblePoints(kSvtId) & 0x7;
  if (! mNHitsPossInner) mNHitsPossInner=track->numberOfPossiblePoints(kIstId) & 0x7;
  mNHitsPossInner|=(track->numberOfPossiblePoints(kSsdId) & 0x3) << 3;
  mNHitsPossInner|=(track->numberOfPossiblePoints(kPxlId) & 0x7) << 5;
  
  mNHitsFit = track->fitTraits().numberOfFitPoints();
  mNHitsFitTpc=0;
  if (track->numberOfPossiblePoints(kTpcId)==track->numberOfPossiblePoints(kFtpcEastId)) {
    // backward compatibility mode, figure out which TPC points are in
    if (track->topologyMap().hasHitInDetector(kTpcId))
      mNHitsFitTpc=track->fitTraits().numberOfFitPoints(kTpcId);
    else if (track->topologyMap().hasHitInDetector(kFtpcEastId))
      mNHitsFitTpc=(1 << 6) + track->fitTraits().numberOfFitPoints(kFtpcEastId);
    else if (track->topologyMap().hasHitInDetector(kFtpcWestId))
      mNHitsFitTpc=(2 << 6) + track->fitTraits().numberOfFitPoints(kFtpcWestId);
  } else {
    // new Ittf chain distinguishes the three TPCs properly
    if ( (tpc_hits=track->fitTraits().numberOfFitPoints(kTpcId)) ) 
      mNHitsFitTpc=tpc_hits; 
    else if ( (tpc_hits=track->fitTraits().numberOfFitPoints(kFtpcEastId)) ) 
      mNHitsFitTpc=(1 << 6) + tpc_hits;
    else if ( (tpc_hits=track->fitTraits().numberOfFitPoints(kFtpcWestId)) ) 
      mNHitsFitTpc=(2 << 6) + tpc_hits;
  }
  mNHitsFitInner=track->fitTraits().numberOfFitPoints(kSvtId) & 0x7;
  if (! mNHitsFitInner) mNHitsFitInner=track->fitTraits().numberOfFitPoints(kIstId) & 0x7;
  mNHitsFitInner|=(track->fitTraits().numberOfFitPoints(kSsdId) & 0x3) << 3;
  mNHitsFitInner|=(track->fitTraits().numberOfFitPoints(kPxlId) & 0x7) << 5;

  mChiSqXY = track->fitTraits().chi2(0);
  mChiSqZ = track->fitTraits().chi2(1);

  if ( track->geometry() && track->geometry()->charge()) {
    mHelix = StMuHelix(track->geometry()->helix(),event->runInfo()->magneticField());
    //if (event->primaryVertex()) {
    if (vertex) {
      mP = momentumAtPrimaryVertex(event,track,vertex);
      Int_t vtx_id=vtxList->IndexOf(vertex);
      if (vtx_id >= 0) {
        mVertexIndex=vtx_id;
      }
      else {
        gMessMgr->Warning() << "Track does not point to a primary vertex" << endm;
        mVertexIndex = -1;
        mDCA = StThreeVectorF(-999,-999,-999);
        mDCAGlobal = StThreeVectorF(-999,-999,-999);
      }
      mDCA = dca(track, vertex);
      if (globalTrack) {
        if (globalTrack->dcaGeometry()) {
          const StDcaGeometry *dcaGeometry = globalTrack->dcaGeometry();
          const StThreeVectorF &pvert = vertex->position();
          Double_t        vtx[3] = {pvert[0],pvert[1],pvert[2]};
          THelixTrack     thelix =  dcaGeometry->thelix();
          thelix.Move(thelix.Path(vtx));
          const Double_t *pos = thelix.Pos();
          const Double_t *mom = thelix.Dir();
          mDCAGlobal = StThreeVectorF(pos[0],pos[1],pos[2]) - vertex->position();
          if (track->type() == global) {
            mDCA = mDCAGlobal;
            mP   = StThreeVectorF(mom[0],mom[1],mom[2]);
            mP  *= dcaGeometry->momentum().mag();
          }
        }
        else {  // No StDcaGeometry; use helix
          mDCAGlobal = dca(globalTrack, vertex);
        }
      }
      else {  // No global track available
        mDCAGlobal = StThreeVectorF(-999,-999,-999);
      }

      mPt = mP.perp();
      mPhi = mP.phi();
      mEta = mP.pseudoRapidity();

      if (!l3) { // L3TRACKS seem to break pid    
        int charge = track->geometry()->charge();
        //      StParticleDefinition* pc = (*mProbabilityPidAlgorithm)( *track, track->pidTraits() );
        mProbabilityPidAlgorithm->processPIDAsFunction(mProbabilityPidCentrality, mDCA.mag(), charge, mP.mag()/charge, mEta, mNHitsDedx, mdEdx);
        mPidProbElectron= pack2UnsignedShort( (charge>0) ? mProbabilityPidAlgorithm->beingPositronProb() : mProbabilityPidAlgorithm->beingElectronProb(), __PROB_SCALE__) ;
        mPidProbPion=     pack2UnsignedShort( (charge>0) ? mProbabilityPidAlgorithm->beingPionPlusProb() : mProbabilityPidAlgorithm->beingPionMinusProb(), __PROB_SCALE__);
        mPidProbKaon=     pack2UnsignedShort( (charge>0) ? mProbabilityPidAlgorithm->beingKaonPlusProb() : mProbabilityPidAlgorithm->beingKaonMinusProb(), __PROB_SCALE__); 
        mPidProbProton=   pack2UnsignedShort( (charge>0) ? mProbabilityPidAlgorithm->beingProtonProb()   : mProbabilityPidAlgorithm->beingAntiProtonProb(), __PROB_SCALE__); 
      }
    }
    else {  // vertex == 0
      mVertexIndex = -1;
      mDCA = StThreeVectorF(-999,-999,-999);
      mDCAGlobal = StThreeVectorF(-999,-999,-999);
      if (globalTrack) {
        if (globalTrack->dcaGeometry()) {
          const StDcaGeometry *dcaGeometry = globalTrack->dcaGeometry();
          THelixTrack     thelix =  dcaGeometry->thelix();
          const Double_t *mom = thelix.Dir();
          if (track->type() == global) {
            mP   = StThreeVectorF(mom[0],mom[1],mom[2]);
            mP  *= dcaGeometry->momentum().mag();
          }
        }
          }
          mPt = mP.perp();
          mPhi = mP.phi();
          mEta = mP.pseudoRapidity();
        }
  }
  fillMuProbPidTraits(event,track);
  mIndex2BTofHit = -1;   // filled later
  fillMuBTofPidTraits(track);

  if ( track->outerGeometry() ) 
    mOuterHelix = StMuHelix(track->outerGeometry()->helix(),event->runInfo()->magneticField());
}

unsigned short StMuTrack::nHitsPoss() const {
  // Add 1 for primary tracks in old files
  if (mNHitsPossTpc==255 && type()==primary)
    return mNHitsPoss+1;
 
  return mNHitsPoss;
}

unsigned short StMuTrack::nHitsPoss(StDetectorId det) const {

  // Backward compatibility for old files
  if (mNHitsPossTpc==255) {
    if (det==kTpcId || det==kFtpcEastId || det==kFtpcWestId)
      return mTopologyMap.hasHitInDetector(det)*mNHitsPoss;
    else 
      return 0;
  }

  // New situation: decode point counts
  switch (det) {
  case kTpcId:
    return ((mNHitsPossTpc & 0xC0)==0)*mNHitsPossTpc;
    break;
  case kFtpcEastId:
    return ((mNHitsPossTpc & 0xC0)==0x40)*(mNHitsPossTpc & 0x3F);
    break;
  case kFtpcWestId:
    return ((mNHitsPossTpc & 0xC0)==0x80)*(mNHitsPossTpc & 0x3F);
    break;
  case kSvtId:
  case kIstId:
    return (mNHitsPossInner & 0x7);
    break;
  case kSsdId:
    return ((mNHitsPossInner >> 3) & 0x3);
    break;
  case kPxlId:
    return ((mNHitsPossInner >> 5) & 0x7);
    break;
  default:
    return 0;
  }

}

unsigned short StMuTrack::nHitsFit(StDetectorId det) const {
  // Backward compatibility for old files
  if (mNHitsFitTpc==255) {
    if (det==kTpcId || det==kFtpcEastId || det==kFtpcWestId)
      return mTopologyMap.hasHitInDetector(det)*mNHitsFit;
    else 
      return 0;
  }

  // New situation: decode point counts
  switch (det) {
  case kTpcId:
    return ((mNHitsFitTpc & 0xC0)==0)*mNHitsFitTpc;
    break;
  case kFtpcEastId:
    return ((mNHitsFitTpc & 0xC0)==0x40)*(mNHitsFitTpc & 0x3F);
    break;
  case kFtpcWestId:
    return ((mNHitsFitTpc & 0xC0)==0x80)*(mNHitsFitTpc & 0x3F);
    break;
  case kSvtId:
  case kIstId:
    return (mNHitsFitInner & 0x7);
    break;
  case kSsdId:
    return ((mNHitsFitInner >> 3) & 0x3);
    break;
  case kPxlId:
    return ((mNHitsFitInner >> 5) & 0x7);
    break;
  default:
    return 0;
  }
}

Float_t StMuTrack::dcaD(Int_t vtx_id) const {
  if ((vtx_id == -1 && mVertexIndex == StMuDst::currentVertexIndex()) ||
       vtx_id == mVertexIndex) {
    StThreeVectorF dir;
    if (mType == global)
      dir = mP;
    else if (globalTrack())
      dir = globalTrack()->p();
    else
      return -999; 
    if (dir.mag() == 0)
      return -999;
    dir /= dir.mag();  
    Float_t cosl = dir.perp();
    return -dir[1]/cosl * mDCAGlobal[0] + dir[0]/cosl * mDCAGlobal[1]; 
  }
  else return -999;
}

Float_t StMuTrack::dcaZ(Int_t vtx_id) const {
  if ((vtx_id == -1 && mVertexIndex == StMuDst::currentVertexIndex()) ||
       vtx_id == mVertexIndex) {
    return mDCAGlobal.z();
  }
  else return -999;
}

StThreeVectorF StMuTrack::dca(Int_t vtx_id) const {
  if (vtx_id == -1)
    vtx_id = StMuDst::currentVertexIndex();
  if (vtx_id==mVertexIndex)
    return mDCA;
  else if (mVertexIndex == -1)  // should not happen
        return gDummy;
    if(((StMuPrimaryVertex*)StMuDst::array(muPrimaryVertex)->UncheckedAt(vtx_id))) 
                return dca(((StMuPrimaryVertex*)StMuDst::array(muPrimaryVertex)->UncheckedAt(vtx_id))->position());
        else return gDummy;
}

StThreeVectorF StMuTrack::dcaGlobal(Int_t vtx_id) const {
  if (vtx_id == -1)
    vtx_id = StMuDst::currentVertexIndex();
  if (vtx_id==mVertexIndex) {
    return mDCAGlobal;
  }
  else if (mVertexIndex == -1)  { // should not happen
    return gDummy;
  }

        const StMuTrack *gTrack = (mType == global) ? this : globalTrack();

        if (gTrack && ((StMuPrimaryVertex*)StMuDst::array(muPrimaryVertex)->UncheckedAt(vtx_id))){
                return gTrack->dca(((StMuPrimaryVertex*)StMuDst::array(muPrimaryVertex)->UncheckedAt(vtx_id))->position());
        }
        else {
    return gDummy;
  }
}

StThreeVectorF StMuTrack::dca(const StThreeVectorF &pos) const {
  const StPhysicalHelixD &hlx = helix();
  double pathlength = hlx.pathLength(pos, false); // do not scan periods
  return hlx.at(pathlength)-pos;
}

StThreeVectorD StMuTrack::dca(const StTrack* track, const StVertex *vertex) const{
  double pathlength = track->geometry()->helix().pathLength( vertex->position(), false ); // do not scan periods
  return track->geometry()->helix().at(pathlength)-vertex->position();
}

StThreeVectorD StMuTrack::momentumAtPrimaryVertex(const StEvent* event, const StTrack* track, const StVertex *vertex) const{
  double pathlength = track->geometry()->helix().pathLength( vertex->position() );
  return track->geometry()->helix().momentumAt(pathlength,event->runInfo()->magneticField()*kilogauss);
}

StPhysicalHelixD StMuTrack::helix() const 
{
  return StPhysicalHelixD(mHelix.p(),mHelix.origin(), mHelix.b()*kilogauss, mHelix.q());
}
  
StPhysicalHelixD StMuTrack::outerHelix() const {
  return StPhysicalHelixD(mOuterHelix.p(),mOuterHelix.origin(), mOuterHelix.b()*kilogauss, mOuterHelix.q());
}  

double StMuTrack::length() const { 
  return fabs( helix().pathLength(StThreeVectorD(mLastPoint)) ); }
double StMuTrack::lengthMeasured() const { 
  return fabs( helix().pathLength(StThreeVectorD(mLastPoint)) - helix().pathLength(StThreeVectorD(mFirstPoint)) ); }

int StMuTrack::bad() const 
{
   if (mFlag <= 0  )            return 10;
   if (mHelix.bad())            return 20;
   if (mOuterHelix.bad())       return 30;
   return 0;
}
#include "StEvent/StProbPidTraits.h"
void StMuTrack::fillMuProbPidTraits(const StEvent* e, const StTrack* t) {
  // get vector of traits; 
  StPtrVecTrackPidTraits traits = t->pidTraits(kTpcId);
  // get the StDedxPidTraits
  StDedxPidTraits* dedxPidTraits =0;
  unsigned int size = traits.size();

  if (StMuDebug::level()>=3) {
      cout << " dedxPidTraits->method() ";
  }
  for (unsigned int i = 0; i < size; i++) {
      if ( !(dedxPidTraits=dynamic_cast<StDedxPidTraits*>(traits[i])) ) continue;
      if (StMuDebug::level()>=3) {
          cout << " " << dedxPidTraits->method();
      }
      if (dedxPidTraits->method() == kTruncatedMeanIdentifier)  {
          mProbPidTraits.setdEdxTruncated( dedxPidTraits->mean() ); 
          mProbPidTraits.setdEdxErrorTruncated( dedxPidTraits->errorOnMean() ); 
          mProbPidTraits.setLog2dX( dedxPidTraits->log2dX() );
      }
      if (dedxPidTraits->method() == kLikelihoodFitIdentifier)  {
          mProbPidTraits.setdEdxFit( dedxPidTraits->mean() ); 
          mProbPidTraits.setdEdxErrorFit( dedxPidTraits->errorOnMean() ); 
          mProbPidTraits.setdEdxTrackLength( dedxPidTraits->length() ); 
          mProbPidTraits.setLog2dX( dedxPidTraits->log2dX() );
      }
  }
  if (StMuDebug::level()>=3) {
      cout << endl;
  }

  // get the StProbPidTraits 
  StProbPidTraits* probPidTraits =0;
  size = traits.size();
  for (unsigned int i = 0; i < size; i++) {
    if ( (probPidTraits=dynamic_cast<StProbPidTraits*>(traits[i])) ) {
      for (int i=0; i<mProbPidTraits.numberOfParticles(); i++)  mProbPidTraits.setProbability(i,probPidTraits->GetProbability(i));
      mProbPidTraits.setNdf(probPidTraits->GetNDF());
    }
  } 
  
}

void StMuTrack::fillMuBTofPidTraits(const StTrack* t) {
  StBTofPidTraits* btofPidTraits = 0;
  StPtrVecTrackPidTraits traits = t->pidTraits(kTofId);
  unsigned int size = traits.size();
  for (unsigned int i = 0; i < size; i++) {
    if ( (btofPidTraits=dynamic_cast<StBTofPidTraits*>(traits[i])) ) {   
      mBTofPidTraits.setBTofPidTraits(btofPidTraits);
    }    
  }
}
#if 0
void StMuTrack::Print(Option_t *option) const {
  //
  // Print out some of the essential track info. 
  // Note: more is stored in the track; may expand printout
  // or specify options some time
  //


  if (mType == global)
    cout << "Global ";
  else if (mType == primary)
    cout << "Primary ";
  else    
    cout << "Other type ";
  cout << "track, id " << mId << ", flag " << mFlag << " (>0 is OK)" << endl;

  if (mVertexIndex != 0) 
    cout << "Not assigned to primary vertex ( vertex Index " << mVertexIndex << " )" << endl;

  cout << "momentum " << mP << endl;
  cout << "eta  " << mEta << ", phi " << mPhi << ", pt " << mPt << endl;
  cout << "DCA  " << mDCA << endl;
  cout << "\t radial " << dcaD()  << ", z " << dcaZ() << endl;
  cout << "global DCA " << mDCAGlobal << endl;
  cout << "Total hits: " << nHits() << ", fitted " << nHitsFit()
       << "\t ( TPC "
       << nHitsFit(kTpcId) << ", FTPC "
       << nHitsFit(kFtpcEastId) + nHitsFit(kFtpcWestId) << ", SVT "  
       << nHitsFit(kSvtId) << ", SSD "
       << nHitsFit(kSsdId) << " ) " << endl;

  cout << "Possible points: " << nHitsPoss() << " \t( TPC "
       << nHitsPoss(kTpcId) << ", FTPC "
       << nHitsPoss(kFtpcEastId) + nHitsPoss(kFtpcWestId) << ", SVT "
       << nHitsPoss(kSvtId) << ", SSD "
       << nHitsPoss(kSsdId) << " ) " << endl;

  cout << "\t first point " << mFirstPoint << endl;
  cout << "\t last point  " << mLastPoint << endl;
  //cout << "\t inner helix origin " << pr_track->helix().origin() << endl; 

}
#endif
ostream&  operator<<(ostream& os, const StMuTrack& v) {
  if      (v.type() == global)  os << "Gl ";
  else if (v.type() == primary) os << "Pr ";
  else                          os << "   ";
  os << Form("id:%5i fl:%5i vx:%3i p:%8.3f %8.3f %8.3f",v.id(),v.flag(),v.vertexIndex(), v.p().x(), v.p().y(), v.p().z());
  os << Form(" q:%2i eta:%6.3f phi:%6.3f pT: %6.3f",v.charge(),v.eta(),v.phi(),v.pt());
  os << Form(" DCA:%6.3f %6.3f %6.3f",v.dca().x(),v.dca().y(),v.dca().z());
  os << Form(" Total hits:%2i fitted:%2i poss:%2i",v.nHits(),v.nHitsFit(),v.nHitsPoss());
  os << Form(" Points F: %6.3f %6.3f %6.3f L: %6.3f %6.3f %6.3f", 
             v.firstPoint().x(),v.firstPoint().y(),v.firstPoint().z(),
             v.lastPoint().x(),v.lastPoint().y(),v.lastPoint().z());
  os << Form(" idT %4i qa %2i",v.idTruth(), v.qaTruth());
  return os;
}

void StMuTrack::Print(Option_t *option) const {cout << *this << endl;}

const StMuTrack* StMuTrack::primaryTrack() const {

        if(mType==1) return this;
        const StMuTrack *prim = 0;

//      StMuDst *tMuDst = 
        
        //For old MuDsts where there was one vertex per event
        if(mIndex2Global==-1){
                int nVer = StMuDst::numberOfPrimaryVertices();
                StMuDst::fixTrackIndicesG(nVer);
        }
        if(mIndex2Global < 0 ) return prim;
        
        if(StMuDst::numberOfPrimaryVertices()==0) return StMuDst::primaryTracks(mIndex2Global);
        int curVer =  StMuDst::currentVertexIndex();
        StMuDst::setVertexIndex(mVertexIndex);
        prim =  StMuDst::primaryTracks(mIndex2Global);
        StMuDst::setVertexIndex(curVer);
        return prim;
}

int StMuTrack::vertexIndex() const {

        //For old MuDsts where there was one vertex per event
        if (StMuDst::numberOfPrimaryVertices()==0){
                if(!(fabs(StMuDst::event()->primaryVertexPosition().x()) < 1.e-5 && fabs(StMuDst::event()->primaryVertexPosition().y()) < 1.e-5 && fabs(StMuDst::event()->primaryVertexPosition().z()) < 1.e-5)){
                        if(primaryTrack()!=0) return 0;
                }
                else return -1;
        }
        if(this->type()==1) return mVertexIndex;
        if(primaryTrack()!=0) {
                int index = primaryTrack()->vertexIndex();
                return index;
        }
        else return -1;
}

TArrayI StMuTrack::getTower(bool useExitRadius,int det) const{ //1=BTOW, 3=BSMDE, 4=BSMDP... Returns TVector tower. tower[0] is module, tower[1] is eta, tower[2] is sub, and tower[3] is id
        
//      StMuTrack* track = this;
        TArrayI tower(4);
        tower[0] = -10;
        tower[1] = -10;
        tower[2] = -10;         
        tower[3] = -10;         

        StThreeVectorD momentum,position;
        Double_t radius;
        
        StEmcGeom* mEmcGeom = StEmcGeom::instance("bemc");
        StEmcGeom* mSmdEGeom= StEmcGeom::instance("bsmde");
        StEmcGeom* mSmdPGeom= StEmcGeom::instance("bsmdp");
        
        if(det==1) radius = mEmcGeom->Radius();
        if(det==2) radius = mEmcGeom->Radius();
        if(det==3) radius = mSmdEGeom->Radius();
        if(det==4) radius = mSmdPGeom->Radius();
                
        StEventSummary& evtSummary = StMuDst::event()->eventSummary();
        Double_t mField = evtSummary.magneticField()/10;
                
        //add 30 cm to radius to find out if track left same tower
        if(useExitRadius) radius += 30.0;
        
        StEmcPosition mEmcPosition;
        bool goodProjection;
        if(this) goodProjection = mEmcPosition.trackOnEmc(&position,&momentum,this,mField,radius);
        else return tower;
        if(goodProjection){
                int m,e,s,id=0;
                float eta=position.pseudoRapidity();
                float phi=position.phi();
                if(det==1){
                  mEmcGeom->getBin(phi,eta,m,e,s);
//                s = abs(s);
                  if(mEmcGeom->getId(m,e,s,id)==0){
                    tower[0] = m;
                    tower[1] = e;
                        tower[2] = s;           
                        tower[3] = id;          
                  }
                }
                else if(det==3){
                  int check=mSmdEGeom->getBin(phi,eta,m,e,s);
                  if(!check){
                    s = abs(s);
                    if(mSmdEGeom->getId(m,e,s,id)==0){
                    tower[0] = m;
                    tower[1] = e;
                        tower[2] = s;               
                        tower[3] = id;          }
                  }
                }
                else if(det==4){
                  int check=mSmdPGeom->getBin(phi,eta,m,e,s);
                  s = abs(s);
                  if(!check){
                    if(mSmdPGeom->getId(m,e,s,id)==0){
                    tower[0] = m;
                    tower[1] = e;
                        tower[2] = s;
                        tower[3] = s;   }
                  }
                }
        }
        return tower;
}

double StMuTrack::energyBEMC() const { //Return energy of negative 100 GeV is no BEMC hit is matched...

        double hitEnergy;
        TArrayI tower = getTower();
        unsigned int iMod = tower[0];
    unsigned int iEta = tower[1];
        int iSub = tower[2];
        if(iMod < 1 ||iMod > 120) return -100.0;

        if (StMuDst::emcCollection()) {
                StEmcDetector   *bemcDet = StMuDst::emcCollection()->detector(kBarrelEmcTowerId);
                StEmcModule     *mod = bemcDet->module(iMod);
                StSPtrVecEmcRawHit&  hits = mod->hits();
                for(unsigned int i=0; i<hits.size();i++){
                        if(hits[i]){
                                if((hits[i]->eta() == iEta) && (hits[i]->sub() == iSub)) {
                                        hitEnergy = hits[i]->energy();
                                        if(hitEnergy > 0) return hitEnergy;
                                        else return -100.0;
                                }
                        }
                }
        }
        return -100.0;
}

bool StMuTrack::matchBEMC() const {
        double mEmcThres = 0.15;
        if (energyBEMC() > mEmcThres) return true;
        return false;
}

ClassImp(StMuTrack)

/***************************************************************************
 *
 * $Log: StMuTrack.cxx,v $
 * Revision 1.48  2012/05/07 14:47:06  fisyak
 * Add handles for track to fast detector matching
 *
 * Revision 1.47  2011/10/17 00:19:14  fisyak
 * Active handing of IdTruth
 *
 * Revision 1.46  2011/04/08 01:25:51  fisyak
 * Add branches for MC track and vertex information, add IdTruth to  tracks and vertices, reserve a possiblity to remove Strange MuDst
 *
 * Revision 1.45  2010/06/02 18:53:13  tone421
 * Fixed a bug. Previously, only events where a vertex was found had global tracks with pt and eta filled properly...
 *
 * Revision 1.44  2010/04/27 20:47:17  tone421
 * Added extra functions for BEMC matching. See this post for more details:
 *
 * http://www.star.bnl.gov/HyperNews-star/get/starsofi/7816.html
 *
 * Revision 1.43  2010/03/08 19:06:51  tone421
 * Two things. Global tracks how are filled with an index to primary at birth. Added StMuDst::fixTrackIndicesG(), which is used for matching the primary track indices to global tracks. Previously, this was quite slow -  see this post:
 *
 * http://www.star.bnl.gov/HyperNews-star/protected/get/starsoft/8092/1/1/1.html
 *
 * for more details.
 *
 * Revision 1.42  2009/12/01 21:56:35  tone421
 * Implemented changes as per http://www.star.bnl.gov/rt2/Ticket/Display.html?id=1734
 *
 * Revision 1.38  2009/02/20 16:37:44  tone421
 * *** empty log message ***
 *
 * Revision 1.36  2009/01/07 20:50:59  tone421
 * Added more pointer protection in dcaGlobal(int) and dca(int)
 *
 * Revision 1.35  2008/05/01 20:53:59  mvl
 * Changed/fixed handling of DCA for global tracks without StDcaGeometry: now calculate DCA based on helix (previous behaviour was to default to (0,0,0); only affected globals without primaries).
 * Changed globalDca() function to also work for globals that are attached to a different primary vertex. Used to return -999, now use helix to calculate.
 *
 * Revision 1.34  2008/03/19 14:51:04  fisyak
 * Add two clone arrays for global and primary track covariance matrices, remove mSigmaDcaD and mSigmaDcaZ
 *
 * Revision 1.33  2007/10/18 03:44:23  mvl
 * Added Ist and Pixel hits to mNPossInner and mNFitInner
 *
 * Revision 1.32  2007/09/18 02:29:58  mvl
 * Added basic printing functionality. For convenience and to assist data consistency checks
 *
 * Revision 1.31  2007/07/12 19:46:20  fisyak
 * Add includes for ROOT 5.16
 *
 * Revision 1.30  2007/04/20 06:23:46  mvl
 * Removed Q-vector flag (will implement utility class)
 *
 * Revision 1.28  2007/01/29 18:34:44  mvl
 * Updates to use StDcaGeometry for global DCA and momentum.
 * Added helper functions for radial and Z component: dcaD and dcaZ.
 * Uncertainties on those are stored in sigmaDcaD and sigmaDcaZ.
 * dcaD and dcaZ only work for the primary vertex to which the track belongs (avoid long extrapolations).
 *
 * Revision 1.27  2007/01/05 20:19:44  jeromel
 * helix(helix()) cannot be interpreted by gcc 3.4.6
 *
 * Revision 1.26  2006/07/27 18:55:41  fisyak
 * Remove DCA hack used in SSD+SVT test production (P06id)
 *
 * Revision 1.24  2005/08/19 19:46:06  mvl
 * Further updates for multiple vertices. The main changes are:
 * 1) StMudst::primaryTracks() now returns a list (TObjArray*) of tracks
 *    belonging to the 'current' primary vertex. The index number of the
 *    'current' vertex can be set using StMuDst::setCurrentVertex().
 *    This also affects StMuDst::primaryTracks(int i) and
 *    StMuDst::numberOfprimaryTracks().
 * 2) refMult is now stored for all vertices, in StMuPrimaryVertex. The
 *    obvious way to access these numbers is from the StMuprimaryVertex structures,
 *    but for ebakcward compatibility a function is provided in StMuEvent as well
 *    (this is the only function taht works for existing MuDst)
 *
 * As an aside, I've also changes the internals of StMuDst::createStEvent and
 * StMuDst::fixTrackIndices() to be able to deal with a larger range of index numbers for tracks as generated by Ittf.
 *
 * BIG FAT WARNING: StMudst2StEventMaker and StMuDstFilterMaker
 * do not fully support the multiple vertex functionality yet.
 *
 * Revision 1.23  2005/07/20 04:08:47  mvl
 * Now also changed dca calculation in the other function (StMuTrack::dca(StThreeVectorF))
 * instead of only the helper function for StTrack, StVertex.
 *
 * Revision 1.22  2005/07/20 03:48:06  mvl
 * Changed dca calculation: no longer scan period of helix to get DCA.
 * This mainly affects tracks with large dcas in events with multiple vertices.
 *
 * Revision 1.21  2005/07/15 21:45:08  mvl
 * Added support for multiple primary vertices (StMuPrimaryVertex). Track Dcas are now calculated with repect to the first vertex in the list (highest rank), but another vertex number can be specified. Tarcks also store the index of the vertex they belong to (StMuTrack::vertexIndex())
 *
 * Revision 1.20  2005/03/17 21:55:00  mvl
 * Added StMuMomentumShiftMaker for applying a magnetic field scaling to the reconstructed MuDst. This class accesses StMuTrack, StMuEvent and StMuHelix and some Strangeness MuDst data members as 'friend'
 *
 * Revision 1.19  2005/02/05 01:08:12  perev
 * zero charge skip
 *
 * Revision 1.18  2004/12/15 21:35:06  mvl
 * Added initialisaion of eta, phi and pt. Affects events without main vertex.
 *
 * Revision 1.17  2004/08/14 00:53:41  mvl
 * Added 1 to possibel points for primary tracks, like in StEvent
 *
 * Revision 1.16  2004/08/12 01:34:25  mvl
 * Fixed typo when extracting number of fitted, possible points per tpc
 *
 * Revision 1.15  2004/08/11 00:51:54  mvl
 * Added topologyMap check for nHitsFit(StDetectorId) and nHitsPoss(StDetectorId) to better handle existing MuDst
 *
 * Revision 1.14  2004/08/10 22:38:30  mvl
 * Extended support for fitspoints in Svt and Tpc to work for old files + fixed some types in previous version.
 *
 * Revision 1.13  2004/08/07 02:44:05  mvl
 * Added support for fitted and possible points in different detectors, for ITTF
 *
 * Revision 1.12  2004/04/14 14:21:53  jeromel
 * Not sure why I made it - but better to preserve the sign
 *
 * Revision 1.11  2004/04/14 03:27:30  jeromel
 * Change init of mNSigma
 *
 * Revision 1.10  2004/04/08 16:21:57  subhasis
 * Initialization ,as suggested by Yuri
 *
 * Revision 1.9  2003/11/07 15:23:26  laue
 * added error on dEdx measurements to the StMuProbPidTraits
 *
 * Revision 1.8  2003/10/30 20:08:13  perev
 * Check of quality added
 *
 * Revision 1.7  2003/10/28 18:57:56  perev
 * BadData protection added
 *
 * Revision 1.6  2003/02/21 14:32:47  laue
 * Yuri's updates to the PID probabilities. dE/dx track length in TPC added
 *
 * Revision 1.5  2002/11/22 18:08:53  laue
 * Bug in logic fixed. Checking for null pointers was wrong.
 *
 * Revision 1.4  2002/11/18 14:29:31  laue
 * update for Yuri's new StProbPidTraits
 *
 * Revision 1.3  2002/09/19 21:54:01  laue
 * fix bug in length() method
 *
 * Revision 1.2  2002/08/20 19:55:49  laue
 * Doxygen comments added
 *
 * Revision 1.1  2002/03/08 17:04:18  laue
 * initial revision
 *
 *
 **************************************************************************/

---