TTrackMC.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TTrackMC.cxx,v 1.5 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TTrackMC.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to have MC information of TTrack.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
/* for copysign */
#if defined( __sparc )
#  if defined( __EXTENSIONS__ )
#    include <cmath>
#  else
#    define __EXTENSIONS__
#    include <cmath>
#    undef __EXTENSIONS__
#  endif
#elif defined( __GNUC__ )
#  if defined( _XOPEN_SOURCE )
#    include <cmath>
#  else
#    define _XOPEN_SOURCE
#    include <cmath>
#    undef _XOPEN_SOURCE
#  endif
#endif
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TTrackMC.h"
// #include "TrkReco/TrkReco.h"
#include <cfloat>
// #include "tables/cdc.h"
#include "MdcTables/MdcTables.h"
 
#if defined( __alpha )
#  define DBL_MIN 2.2250738585072014E-308
#  define FLT_MIN 1.175494351E-38F
#endif
 
TTrackMC::TTrackMC( const TTrack& t )
    : _t( t )
    , _hep( 0 )
    , _hepID( -1 )
    , _wireFraction( -999. )
    , _wireFractionHEP( -999. )
    , _charge( false )
    , _ptFraction( -999. )
    , _pzFraction( -999. )
    , _ptResidual( -999. )
    , _pzResidual( -999. )
    , _ptPull( -999. )
    , _pzPull( -999. )
    , _state( 0 )
    , _quality( 0 ) {}
 
TTrackMC::~TTrackMC() {}
 
void TTrackMC::update( void ) {
  _state   = 0;
  _quality = 0;
  /*
      //...Prepare counters...
      unsigned nHep = TTrackHEP::list().length();
      unsigned nTrk = TrkReco::getTrkReco()->tracks().length();
      unsigned * N1 = (unsigned *) malloc(nHep * sizeof(unsigned));
      float * F1 = (float *) malloc(nHep * sizeof(float));
      unsigned N2 = 0;
      for (unsigned i = 0; i < nHep; i++) {
          N1[i] = 0;
          F1[i] = 0.;
      }
 
      //...Prepare for fraction F1...
  //  const AList<TMLink> & cores = _t.cores();
      const AList<TMLink> & cores = _t.finalHits();
      unsigned nCores = cores.length();
      for (unsigned i = 0; i < nCores; i++) {
          TMLink * t = cores[i];
          int hepID = t->hit()->mc()->hep()->id();
          ++N1[hepID];
      }
 
      //...Calculate fraction F1 and find the best HEP...
      int bestHep = -1;
      TTrackHEP * hep = 0;
      float bestF1 = 0.;
      for (unsigned i = 0; i < nHep; i++) {
          if (nCores) F1[i] = (float) N1[i] / (float) nCores;
          if (F1[i] > bestF1) {
              bestHep = i;
              bestF1 = F1[i];
          }
      }
 
      //...Check HEP...
      float F2 = 0.;
      if (bestHep != -1) {
          hep = TTrackHEP::list()[bestHep];
          unsigned nAll = 0;
          for (unsigned i = 0; i < hep->hits().length(); i++) {
              const TMDCWireHit * hit = hep->hits()[i]->hit();
              if (! hit) continue;
              if (hit->state() & WireHitInvalidForFit) continue;
 
              ++nAll;
              if (hit->track() == & _t) ++N2;
          }
 
          //...Calculate fraction F2...
          if (nAll) F2 = (float) N2 / (float) nAll;
      }
 
      //...Store results...
      _hepID = bestHep;
      _hep = hep;
      _wireFraction = bestF1;
      _wireFractionHEP = F2;
 
      //...Compare charge and momentum...
      compare();
 
      //...Classification...
      classify();
 
      //...Termination...
      free(N1);
      free(F1);
  */
}
 
void TTrackMC::dump( const std::string& msg, const std::string& pre ) const {
  std::cout << msg;
  std::cout << _t.name() << ":";
  std::cout << "state=" << _state << ":";
  if ( _quality & TTrackGood ) std::cout << "good   :";
  else if ( _quality & TTrackGhost ) std::cout << "ghost  :";
  else if ( _quality & TTrackBad ) std::cout << "bad    :";
  else if ( _quality & TTrackCharge ) std::cout << "bad    :";
  else if ( _quality & TTrackGarbage ) std::cout << "garbage:";
  else std::cout << "classification error:";
  bitDisplay( _quality, 23, 15 );
  std::cout << ":";
  std::cout << _hepID << ":";
  std::cout << _wireFraction << "," << _wireFractionHEP << ":";
  std::cout << _ptFraction << "," << _pzFraction;
  std::cout << std::endl;
}
 
void TTrackMC::compare( void ) {
  if ( !_hep ) return;
 
  //...Get charge of HEP particle(This part should be done by LUCHARGE)...
  int id  = _hep->pType();
  int aId = abs( id );
  if ( aId == 11 || aId == 13 || aId == 15 ) id *= -1;
 
  //...Compare charge...
  if ( (int)_t.charge() * id > 0 ) _charge = true;
 
  //...Get hep mom. at the inner-most hit...
  //  AList<TMLink> list = _t.cores();
  AList<TMLink> list  = _t.finalHits();
  unsigned      n     = list.length();
  bool          found = false;
  Vector3       pHep;
  Vector3       vHep;
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink* inner = InnerMost( list );
    if ( inner->hit()->mc()->hep() == _hep )
    {
      pHep  = inner->hit()->mc()->momentum();
      vHep  = inner->hit()->mc()->entrance();
      found = true;
      break;
    }
    list.remove( inner );
  }
  Helix hHep = Helix( vHep, pHep, copysign( 1., id ) );
  hHep.pivot( _t.helix().pivot() );
  pHep = hHep.momentum();
 
  //...For debug...
  if ( !found )
  {
    std::cout << "TTrackMC::compare !!! something wrong with mc hits" << std::endl;
 
    //...For debug...
    //      list = _t.cores();
    //      for (unsigned i = 0; i < list.length(); i++) {
    //          TMLink * inner = innerMost(list);
    //          std::cout << i << " " << inner << ":" << inner->hit()->mc()->hep();
    //          std::cout << " " << _hep << std::endl;
    //          if (inner->hit()->mc()->hep() == _hep) {
    //          pHep = inner->hit()->mc()->momentum();
    //          break;
    //          }
    //          list.remove(inner);
    //      }
    //      TMLink * t = 0;
    //      t->hit();
    //...For debug end...
 
    return;
  }
 
  //...Fill caches...
  _residual = _t.p() - pHep;
  _cosOpen  = pHep.unit().dot( _t.p().unit() );
 
  //...Compare pt...
  double pt      = _t.pt();
  double ptHep   = sqrt( pHep.x() * pHep.x() + pHep.y() * pHep.y() );
  _ptResidual    = pt - ptHep;
  const Helix& h = _t.helix();
  Vector       dPt( 5, 0 );
  dPt[2]          = -1. / ( h.kappa() * h.kappa() );
  double ptError2 = h.Ea().similarity( dPt );
  if ( ptError2 < 0.0 )
  { std::cout << h.kappa() << " " << h.Ea() << " dPt=" << dPt << std::endl; }
  double ptError = ( ptError2 > 0. ) ? sqrt( ptError2 ) : ( DBL_MIN );
  _ptPull        = ( ptError2 > 0. ) ? ( _ptResidual ) / ptError : ( FLT_MAX );
  _ptFraction    = ( fabs( ptHep ) > ( FLT_MIN ) ) ? _ptResidual / ptHep : 0.0;
 
  //...Compare pz...
  double pz    = _t.pz();
  double pzHep = pHep.z();
  _pzResidual  = pz - pzHep;
  Vector dPz( 5, 0 );
  dPz[2]          = -h.tanl() / ( h.kappa() * h.kappa() );
  dPz[4]          = 1. / h.kappa();
  double pzError2 = h.Ea().similarity( dPz );
  if ( pzError2 < 0.0 )
  { std::cout << h.kappa() << " " << h.Ea() << " dPz=" << dPz << std::endl; }
  double pzError = ( pzError2 > 0. ) ? sqrt( pzError2 ) : ( DBL_MIN );
  _pzPull        = ( pzError2 > 0. ) ? ( _pzResidual ) / pzError : ( FLT_MAX );
  _pzFraction    = ( abs( pzHep ) > FLT_MIN ) ? ( _pzResidual / pzHep ) : ( FLT_MAX );
}
 
void TTrackMC::classify( void ) {
  _state |= TTrackClassified;
  _quality = TTrackGarbage;
 
  //...HEP matching...
  if ( !_hep ) return;
  _quality |= TTrackHep;
  if ( fabs( _ptFraction ) < .1 ) _quality |= TTrackPt;
  if ( fabs( _pzFraction ) < .1 ) _quality |= TTrackPz;
  float momResidual = sqrt( _ptResidual * _ptResidual + _pzResidual * _pzResidual );
 
  if ( ( momResidual < 0.100 ) && ( _cosOpen > 0.99 ) ) _quality |= TTrackMatchingLoose;
  if ( ( momResidual < 0.020 ) && ( _cosOpen > 0.998 ) ) _quality |= TTrackMatchingTight;
 
  //...Wire fraction...
  if ( _wireFraction < 0.8 ) return;
  _quality |= TTrackWire;
  _quality |= TTrackCharge;
 
  //...Charge matching...
  if ( !_charge ) return;
  _quality |= TTrackBad;
 
  //...Momentum matching...
  if ( _quality & TTrackMatchingLoose ) _quality |= TTrackGhost;
 
  //...TTrackGood is set by TrkReco after checking uniqueness...
}
 
std::string TTrackMC::qualityString( void ) const { return TrackMCQualityString( _quality ); }
 
std::string TrackMCStatus( unsigned quality ) {
  //...This is a local function to hide from user...
 
  std::string matching;
  if ( quality & TTrackHep )
  {
    if ( quality & TTrackMatchingTight ) matching += "tight";
    else if ( quality & TTrackMatchingLoose ) matching += "loose";
    else matching = "bad";
  }
  return TrackMCQualityString( quality ) + " " + matching;
}
 
std::string TrackMCStatus( const TTrackMC& m ) { return TrackMCStatus( m.quality() ); }
 
std::string TrackMCStatus( const MdcRec_mctrk& m ) { return TrackMCStatus( m.quality ); }
 
std::string TrackMCQualityString( unsigned quality ) {
  if ( quality & TTrackGood ) return std::string( "Good" );
  else if ( quality & TTrackGhost ) return std::string( "Ghost" );
  else if ( quality & TTrackBad ) return std::string( "Bad" );
  else if ( quality & TTrackCharge ) return std::string( "Charge" );
  else if ( quality & TTrackGarbage ) return std::string( "Garbage" );
  return std::string( "Unknown" );
}