/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/TrkReco/TrkReco-00-08-59-patch4-slc6tag/src/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 <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");
}

---