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    : [email protected]
//-----------------------------------------------------------------------------
// 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 <cfloat>
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TrkReco.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TTrackMC.h"
#include "TrkReco/TTrackHEP.h"
//#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");
}