TPerfectFinder.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TPerfectFinder.cxx,v 1.9 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TPerfectFinder.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : [email protected]
//-----------------------------------------------------------------------------
// Description : A class to find tracks using MC info.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "TrkReco/TPerfectFinder.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TTrack.h"
//#include "helix/Helix.h"
//#include "TrkReco/Helix.h"
#include "TrackUtil/Helix.h"
 
TPerfectFinder::TPerfectFinder(int perfectFitting,
                   float maxSigma,
                   float maxSigmaStereo,
                   unsigned fittingFlag) 
: _perfectFitting(perfectFitting),
  _builder("conformal builder",
       maxSigma,
       maxSigmaStereo,
       0,
       0,
       0,
       fittingFlag),
  _fitter("helix fitter") {
}
 
TPerfectFinder::~TPerfectFinder() {
}
 
std::string
TPerfectFinder::version(void) const {
    return "2.04";
}
 
void
TPerfectFinder::dump(const std::string & msg, const std::string & pre) const {
    std::cout << pre;
    TFinderBase::dump(msg);
}
 
int
TPerfectFinder::doit(const AList<TMDCWireHit> & axialHits,
             const AList<TMDCWireHit> & stereoHits,
             AList<TTrack> & tracks,
             AList<TTrack> & tracks2D) {
 
    //...Preparations...
    static const HepPoint3D dummy(0, 0, 0);
    CAList<TMDCWireHit> hits;
    hits.append(axialHits);
    hits.append(stereoHits);
 
    //...Make a list of HEP track...
    CAList<TTrackHEP> heps;
    const unsigned nHits = hits.length();
    for (unsigned i = 0; i < nHits; i++) {
    const TMDCWireHitMC * const mc = hits[i]->mc();
    if (! mc) continue;
    const TTrackHEP * const hep = mc->hep();
 
    if (! hep) continue;
    heps.append(hep);
    }
    heps.purge();
 
    //...Create tracks...
    const AList<TMDCWireHitMC> & mcHits = TMDC::getTMDC()->hitsMC();
    const unsigned nHitsMC = mcHits.length();
    const unsigned nHeps = heps.length();
    for (unsigned i = 0; i < nHeps; i++) {
    const TTrackHEP & hep = * heps[i];
    const float chg = charge(hep.pType());
    AList<TMLink> hepLinks;
    HepPoint3D posIn, posOut;
    HepVector3D momIn, momOut;
    float rMin = 99999;
    float rMax = 0.;
    TMLink * linkIn;
    unsigned lastLayer = 0;
    for (unsigned j = 0; j < mcHits.length(); j++) {
        const TMDCWireHitMC * const mc = mcHits[j];
        if (! mc) continue;
        if (! mc->hit()) continue;
        if (& hep != mc->hep()) continue;
        if (! hits.hasMember(mc->hit())) continue;
 
        //...Remove hits by curl back...(assuming order of mc hits)
        if (mc->wire()->layerId() < lastLayer) break;
        lastLayer = mc->wire()->layerId();
        HepPoint3D ent = mc->entrance();
        HepVector3D dir = mc->direction();
        ent.setZ(0.);
        dir.setZ(0.);
        if ((ent.unit()).dot(dir.unit()) < 0.5) continue;
 
        //...Hit...
        TMLink * l = new TMLink(0, mc->hit(), dummy);
        l->leftRight(mc->leftRight());
        hepLinks.append(l);
        _links.append(l);
 
        //...Check r to get MC mom...
        const float r = ent.mag();
        if (r < rMin) {
        posIn = mc->entrance();
        momIn = mc->momentum();
        rMin = r;
        linkIn = l;
        }
        if (r > rMax) {
        posOut = mc->entrance();
        momOut = mc->momentum();
        rMax = r;
        }
    }
    if (_links.length() == 0) continue;
 
    //...Do perfect fitting...
    TTrack * t = 0;
//      Helix h(posIn, momIn, chg);
    Helix h(posOut, momOut, chg);
//      h.pivot(posIn);
    h.pivot(linkIn->wire()->xyPosition());
    t = new TTrack(h);
    t->append(hepLinks);
    t->assign(0);
 
//          Helix hX(posOut, momOut, chg);
    
//          std::cout << "    gen@cdc i = " << h.a() << std::endl;
//          std::cout << "       pivot  = " << h.pivot() << std::endl;
//          std::cout << "       mom    = " << h.momentum() << std::endl;
//          std::cout << "    gen@cdc o = " << hX.a() << std::endl;
//          std::cout << "       pivot  = " << hX.pivot() << std::endl;
//          std::cout << "       mom    = " << hX.momentum() << std::endl;
//          std::cout << "    momOut    = " << momOut << std::endl;
//          std::cout << "    posOut    = " << posOut << std::endl;
//          std::cout << "    ptOut     = " << momOut.perp() << std::endl;
//          std::cout << "    costhOut  = " << momOut.z() / momOut.mag();
//          std::cout << std::endl;
 
    if (_perfectFitting == 0) {
        std::cout << "special test in perfect finder" << std::endl;
        AList<TMLink> tmp;
        for (unsigned i = 0; i < t->nLinks(); i++) {
        bool rem = false;
        if (t->links()[i]->wire()->name() == "0-56")
            rem = true;
        else if (t->links()[i]->wire()->name() == "2-56")
            rem = true;
        else if (t->links()[i]->wire()->name() == "5-57")
            rem = true;
        else if (t->links()[i]->wire()->name() == "6=68")
            rem = true;
        else if (t->links()[i]->wire()->name() == "7=69")
            rem = true;
        else if (t->links()[i]->wire()->name() == "8=68")
            rem = true;
        else if (t->links()[i]->wire()->name() == "10-85")
            rem = true;
        else if (t->links()[i]->wire()->name() == "20-126")
            rem = true;
        else if (t->links()[i]->wire()->name() == "39-210")
            rem = true;
        else if (t->links()[i]->wire()->name() == "41=221")
            rem = true;
        else if (t->links()[i]->wire()->name() == "43=221")
            rem = true;
        else if (t->links()[i]->wire()->name() == "46-251")
            rem = true;
        else if (t->links()[i]->wire()->name() == "48-251")
            rem = true;
        if (rem) {
            std::cout << t->links()[i]->wire()->name() << " removed" << std::endl;
            tmp.append(* t->links()[i]);
        }
        }
        t->remove(tmp);
        t->dump("detail","@lastHit");
        Helix & h = (Helix &) t->helix();
        h.pivot(HepVector3D(0., 0., 0.));
        t->dump("detail","@origin ");
        _fitter.fit(* t);
        t->dump("detail","fitted  ");
        static const HepVector3D p0(1.226, -1.025, 0.120);
        std::cout << "Pdif mag=" << (t->p() - p0).mag() << std::endl;
    }
 
    tracks.append(t);
    }
    return 0;
}
 
void
TPerfectFinder::clear(void) {
    HepAListDeleteAll(_links);
}
 
float
TPerfectFinder::charge(int pType) const {
    float chg;
 
    if (pType == 11) chg = -1;
    else if (pType == -11) chg = 1;
    else if (pType == 13) chg = -1;
    else if (pType == -13) chg = 1;
    else if (pType == 211) chg = 1;
    else if (pType == -211) chg = -1;
    else if (pType == 321) chg = 1;
    else if (pType == -321) chg = -1;
    else if (pType == 2212) chg = 1;
    else if (pType == -2212) chg = -1;
    else {
    std::cout << "TPerfectFinder !!! charge of particle type=";
    std::cout << pType << " is unknown" << std::endl;
    }
 
    return chg;
}