TBuilder0.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TBuilder0.cxx,v 1.9 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TBuilder0.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : [email protected]
//-----------------------------------------------------------------------------
// Description : A class to build a track.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#ifdef TRKRECO_DEBUG_DETAIL
#ifndef TRKRECO_DEBUG
#define TRKRECO_DEBUG
#endif
#endif
#include "TrkReco/TBuilder0.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TCircle.h"
#include "TrkReco/TLine0.h"
#include "TrkReco/TTrack.h"
////#include "TrkReco/TSegment0.h"
#ifdef TRKRECO_WINDOW
#include "TrkReco/TWindow.h"
#endif
 
TBuilder0::TBuilder0(const std::string & a)
: _name(a),
  _fitter("TBuilder0 Fitter"),
  _salvageLevel(30.),
  _stereoZ3(20.),
  _stereoZ4(20.),
  _stereoChisq3(15.),
  _stereoChisq4(9.),
  _stereoMaxSigma(30.) {
}
 
TBuilder0::TBuilder0(const std::string & a, float salvageLevel)
: _name(a),
  _fitter("TBuilder0 Fitter"),
  _salvageLevel(salvageLevel),
  _stereoZ3(20.),
  _stereoZ4(20.),
  _stereoChisq3(15.),
  _stereoChisq4(9.),
  _stereoMaxSigma(30.) {
}
 
TBuilder0::TBuilder0(const std::string & a,
           float stereoZ3,
           float stereoZ4,
           float stereoChisq3,
           float stereoChisq4,
           float stereoMaxSigma,
           unsigned corrections,
           float salvageLevel)
: _name(a),
  _fitter("TBuilder0 Fitter"),
  _stereoZ3(stereoZ3),
  _stereoZ4(stereoZ4),
  _stereoChisq3(stereoChisq3),
  _stereoChisq4(stereoChisq4),
  _stereoMaxSigma(stereoMaxSigma),
  _salvageLevel(salvageLevel) {
//    if (corrections & 1) _fitter.sag(true);
//    if (corrections & 2) _fitter.propagation(true);
//    if (corrections & 4) _fitter.tof(true);
//    if (corrections & 8) _fitter.freeT0(true);
}
 
TBuilder0::~TBuilder0() {
}
 
void
TBuilder0::dump(const std::string & msg, const std::string & pre) const {
}
 
TTrack *
TBuilder0::buildRphi(const AList<TMLink> & list) const {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << " ... building a rphi track" << std::endl;
    std::cout << _name << " ... selecting good hits" << std::endl;
//      TTrackBase tmp;
//      tmp.append(list);
//      rphiWindow->append(tmp, leda_red);
#endif
 
    //...Check # of links...
    if (list.length() < _circleSelector.nLinks()) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << " ... rejected by nLinks(";
    std::cout << list.length() << ") < ";
    std::cout << _circleSelector.nLinks() << std::endl;
#endif  
    }
 
    //...Select core hits...
    AList<TMLink> cores = list;
    selectHits(cores);
    if (cores.length() < 5) cores = list;
    //cout<<"list: "<<list.length()<<"  core: "<<cores.length()<<endl;
 
    //...Core check...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << " ... checking cores : cores=" << std::endl;
    Dump(cores, "detail");
#endif
    unsigned sLinks = SuperLayer(list);
    unsigned sUsed = SuperLayer(cores);
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    super layer ptn=" << sLinks;
    std::cout << ",super layer used=" << sUsed << std::endl;
#endif
    if (sLinks != sUsed) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "    ... appending hits to cores" << std::endl;
#endif
    unsigned diff = sLinks - sUsed;
    unsigned asl[5] = {2, 3, 4, 9, 10};
    for (unsigned j = 0; j < 5; j++) {
        //if (diff & (1 << ((5 - j) * 2))) {
        if (diff & (1 << asl[5 - j])) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "    super layer " << (5 - j) * 2 << "searching";
        std::cout << std::endl;
#endif
        AList<TMLink> links = SameSuperLayer(list, asl[5 - j]);
        TMLink * best = NULL;
        double bestD = 999.;
        for (unsigned i = 0; i < links.length(); i++) {
            double dist = links[i]->hit()->drift();
            if (dist < 0.02) {
#ifdef TRKRECO_DEBUG_DETAIL
            std::cout << "    " << links[i]->wire()->name();
            std::cout << " appended (small dist)" << std::endl;
#endif
            cores.append(links[i]);
            continue;
            }
            if (dist < bestD) {
            best = links[i];
            bestD = dist;
            }
        }
        if (best) {
            cores.append(best);
#ifdef TRKRECO_DEBUG_DETAIL
            std::cout << "    " << best->wire()->name();
            std::cout << " appended (best)" << std::endl;
#endif
        }
        }
    }
    }
 
    //...Check cores again...
    unsigned nCores = cores.length();
    AList<TMLink> realCores;
    for (unsigned i = 0; i < nCores; i++) {
    TMLink * l = cores[i];
    if (l->hit()->state() & WireHitFittingValid)
        realCores.append(l);
    }
    if (NSuperLayers(realCores) < 2) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... rejected by small number of super layers" << std::endl;
#endif  
    return NULL;
    }
    if (NLayers(realCores) < 5) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... rejected by small number of layers" << std::endl;
#endif  
    return NULL;
    }
 
    //...Make a circle...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << " ... making a circle : #cores=" << cores.length();
    std::cout << std::endl;
#endif
    AList<TMLink> hits = list;
    hits.remove(cores);
    TCircle c(cores);
 
    //...Test it...
    if (! _circleSelector.preSelect(c)) return NULL;
 
    //...Fitting...
    int err = c.fit();
    if (err < 0) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... rejected by failure of the 1st fit : ";
    std::cout << "err = " << err << std::endl;
#endif
    return NULL;
    }
 
    //...Test it...
    if (! _circleSelector.select(c)) return NULL;
 
    //...Make a track...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << " ... making a track" << std::endl;
#endif
    TTrack * t = new TTrack(c);
    if (! _trackSelector.preSelect(* t)) {
    delete t;
    return NULL;
    }
 
    //...Fitting...
    AList<TMLink> bad;
    //cout<<"track0: "<<t->nLinks()<<endl;
    err = fit(* t);
    if (err < 0) goto discard;
    t->refine(bad, _trackSelector.maxSigma() * 100.);
    //cout<<"track1: "<<t->nLinks()<<endl;
    err = fit(* t);
    if (err < 0) goto discard;
#ifdef TRKRECO_DEBUG_DETAIL
    t->dump("detail", "    1st> ");
#endif
 
    //...Test it...
    if (! _trackSelector.select(* t)) goto discard;
 
    //...Try to append non-core hits...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << " ... appending non-core hits" << std::endl;
#endif
    t->appendByApproach(hits, sqrt(_trackSelector.maxSigma()));
    err = fit(* t);
    if (err < 0) goto discard;
//    cout<<"  2D track1: "<<t->links().length()<<endl;
    t->refine(bad, _trackSelector.maxSigma() * 10.);
    err = fit(* t);
    if (err < 0) goto discard;
//    cout<<"  2D track2: "<<t->links().length()<<endl;
    t->refine(bad, _trackSelector.maxSigma());
    err = fit(* t);
    if (err < 0) goto discard;
//    cout<<"  2D track3: "<<t->links().length()<<endl;
#ifdef TRKRECO_DEBUG_DETAIL
    t->dump("detail", "    2nd> ");
#endif
 
    //...Test it...
    if (! _trackSelector.select(* t)) goto discard;
 
    //...OK...
#ifdef TRKRECO_DEBUG_DETAIL
//      rphiWindow->append(* t, leda_blue);
//      rphiWindow->draw();
//      rphiWindow->wait();
//      rphiWindow->remove(tmp);
//      rphiWindow->remove(* t);
#endif
 
    return t;
 
    //...Something happened...
discard:
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... rejected by fitting failure : ";
    std::cout << " err = " << err << std::endl;
//      rphiWindow->append(* t, leda_blue);
//      rphiWindow->draw();
//      rphiWindow->wait();
//      rphiWindow->remove(tmp);
//      rphiWindow->remove(* t);
#endif
    delete t;
    return NULL;
}
 
void
TBuilder0::selectHits(AList<TMLink> & list) const {
    AList<TMLink> bad;
    for (unsigned i = 0; i < list.length(); i++) {
    unsigned state = list[i]->hit()->state();
    if ((! (state & WireHitContinuous)) ||
        (! (state & WireHitIsolated))) {
        bad.append(list[i]);
        continue;
    }
/*  if ((! (state & WireHitPatternLeft)) &&
        (! (state & WireHitPatternRight))) {
        bad.append(list[i]);
        continue;
    }
*/    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(TBuilder0::selectHits) ... dump of candidates" << std::endl;
    Dump(list, "detail");
#endif
 
    list.remove(bad);
}
 
TTrack *
TBuilder0::buildStereo0(TTrack & track, const AList<TMLink> & list) const {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... dump of stereo candidate hits" << std::endl;
    Dump(list, "sort flag", "    ");
#endif
 
    //...Check # of links...
    if (list.length() < _lineSelector.nLinksStereo()) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... rejected by nLinks(";
    std::cout << list.length() << ") < ";
    std::cout << _lineSelector.nLinks() << std::endl;
#endif  
    return NULL;
    }
 
    //...Calculate s and z for every links...
    unsigned n = list.length();
    AList<TMLink> forLine;
    for (unsigned i = 0; i < n; i++) {
    TMLink * l = list[i];
    TMLink * t = new TMLink(* l);
 
    //...Assuming wire position...
    t->leftRight(2);
    int err = track.szPosition(* t);
    if (err) {
        delete t;
        continue;
    }
 
    //...Store the sz link...
    t->link(l);
    forLine.append(t);
    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... dump of sz links" << std::endl;
    Dump(forLine, "sort flag", "    ");
#endif    
 
    //...Check # of sz links...
    if (forLine.length() < _lineSelector.nLinksStereo()) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... rejected by sz nLinks(";
    std::cout << forLine.length() << ") < ";
    std::cout << _lineSelector.nLinks() << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    return NULL;
    }
 
    //...Make a line...
    unsigned nLine = forLine.length();
    TLine0 line(forLine);
    int err = line.fit2sp();
    if (err < 0) {
    err = line.fit2p();
    if (err < 0) err = line.fit();
    }
 
    //...Linear fit...
    if (err < 0) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... linear fit failure. nLinks(";
    std::cout << forLine.length() << ")" << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    return NULL;
    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... dump of left-right" << std::endl;
#endif
 
    //...Decide Left or Right...
    AList<TMLink> forNewLine;
    for (unsigned i = 0; i < nLine; i++) {
    TMLink * t = forLine[i];
    TMLink * tl = new TMLink(* t);
    TMLink * tr = new TMLink(* t);
 
    tl->leftRight(WireHitLeft);
    tr->leftRight(WireHitRight);
 
    int err = track.szPosition(* tl);
    if (err) {
        delete tl;
        tl = NULL;
    }
    err = track.szPosition(* tr);
    if (err) {
        delete tr;
        tr = NULL;
    }
    if ((tl == NULL) && (tr == NULL)) continue;
 
    TMLink * best;
    if (tl == NULL) best = tr;
    else if (tr == NULL) best = tl;
    else {
        if (line.distance(* tl) < line.distance(* tr)) {
        best = tl;
        delete tr;
        }
        else {
        best = tr;
        delete tl;
        }
    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ";
    std::cout << t->wire()->layerId() << "-";
    std::cout << t->wire()->localId();
    if (tl != NULL)
        std::cout << ",left " << tl->position() << "," << line.distance(* tl);
    if (tr != NULL)
        std::cout << ",right " << tr->position() << "," << line.distance(* tr);
    std::cout << std::endl;
#endif
 
    best->link(t->link());
    forNewLine.append(best);
    }
 
    //...Check # of sz links...
    if (forNewLine.length() < _lineSelector.nLinksStereo()) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... rejected by lr nLinks(";
    std::cout << forNewLine.length() << ") < ";
    std::cout << _lineSelector.nLinks() << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    HepAListDeleteAll(forNewLine);
    return NULL;
    }
 
    //...Create new line...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... creating a new line" << std::endl;
#endif  
    unsigned nNewLine = forNewLine.length();
    TLine0 newLine(forNewLine);
 
    //... Remove extremely bad points
    newLine.removeChits();
 
    //...Make a seed track again
    err = newLine.fit2sp();
    if (err < 0) {
    err = newLine.fit2p();
    if (err < 0) err = newLine.fit();
    }
 
    //...Linear fit...
    if (err < 0) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... 2nd linear fit failure. nLinks(";
    std::cout << forNewLine.length() << ")" << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    HepAListDeleteAll(forNewLine);
    return NULL;
    }
 
    //...Remove bad points...
    AList<TMLink> bad;
    newLine.refine(bad, 40.);
    err = newLine.fit();
    newLine.refine(bad, 30.);
    err = newLine.fit();
    newLine.refine(bad, 20.);
    err = newLine.fit();
    newLine.refine(bad, 10.);
    err = newLine.fit();
    float R = fabs(track.helix().curv());
    if (R > 80.) {
    newLine.refine(bad, 5.);
    err = newLine.fit();
    }
 
    //...Linear fit again...
    if (err < 0) {
    HepAListDeleteAll(forLine);
    HepAListDeleteAll(forNewLine);
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    appendStereo cut ... new line 2nd linear fit failure. ";
    std::cout << "# of links = " << n << "," << nLine;
    std::cout << "," << nNewLine << std::endl;
#endif  
    return NULL;
    }
 
    //...3D fit...
    const AList<TMLink> & good = newLine.links();
    unsigned nn = good.length();
    for (unsigned i = 0; i < nn; i++) {
    track.append(* good[i]->link());
    }
    Vector a(5);
    a = track.helix().a();
    a[4] = track.charge() * newLine.a();
    track._helix->a(a);
 
    //...Refine...
    err = fit(track);
    track.refine(bad, _trackSelector.maxSigma() * 100.);
    err = fit(track);
    track.refine(bad, _trackSelector.maxSigma() * 10.);
    err = fit(track);
    track.refine(bad, _trackSelector.maxSigma());
    err = fit(track);
 
    //...Test it...
    if (! _trackSelector.select(track)) {
    HepAListDeleteAll(forLine);
    HepAListDeleteAll(forNewLine);
    return NULL;
    }
 
    //...Termination...
    HepAListDeleteAll(forLine);
    HepAListDeleteAll(forNewLine);
    return & track;
}
 
TTrack *
TBuilder0::buildStereo(TTrack & track, const AList<TMLink> & list) const {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << _name << "(stereo) ... building in 3D" << std::endl;
    std::cout << "... dump of stereo candidate hits" << std::endl;
    Dump(list, "sort flag", "    ");
#endif
 
    //...Check # of links...
    if (list.length() < _lineSelector.nLinksStereo()) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... rejected by nLinks(";
    std::cout << list.length() << ") < ";
    std::cout << _lineSelector.nLinks() << std::endl;
#endif  
    return NULL;
    }
 
    //...Setup...
    int ichg;
    if (track.helix().curv() > 0.) ichg = 1;
    else                           ichg = -1;
    unsigned nlyr[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    unsigned llyr[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 
    //...Check 
    //... # of hits in a wire, 
    //... 2 contineus hits or not
    //... 
    //...and Calculate s and z for every links...
    unsigned n = list.length();
    AList<TMLink> forLine;
    unsigned iforLine = 0;
    for (unsigned i = 0; i < n; i++) {
    TMLink * l = list[i];
 
        if (i < n-1) {
        TMLink * lnext = list[i + 1];
 
        //...2 consective hits ?...
        // int LorR = check2CnHits(* l, * lnext, ichg);
        int LorR = consectiveHits(* l, * lnext, ichg);
//...For debug...
//          if (check2CnHits(* l, * lnext, ichg) !=
//          consectiveHits(* l, * lnext, ichg)) {
//          AList<TMLink> tmp;
//          tmp.append(l);
//          tmp.append(lnext);
//          std::cout << "!!! consective diff !!!" << std::endl;
//          std::cout << "    original = " << check2CnHits(* l, * lnext, ichg);
//          std::cout << ", iw = " << consectiveHits(* l, * lnext, ichg);
//          std::cout << std::endl;
//          Dump(tmp, "detail");
//          }
//...For debug end...
 
        //...Left/right solved by two consective hits...
        if (LorR == WireHitLeft || LorR == WireHitRight) {
 
        //... Set Left/Right
        if (LorR == WireHitLeft) {
            l->leftRight(WireHitLeft);
            lnext->leftRight(WireHitRight);
        }
        else {
            l->leftRight(WireHitRight);
            lnext->leftRight(WireHitLeft);
        }             
 
        //...Calculate z...
        int err1 = track.szPosition(* l);
        int err2 = track.szPosition(* lnext);
        if(err1 == 0 && err2 == 0){
            double deltaZ = fabs(l->position().y() -
                     lnext->position().y());
            if (deltaZ < 1.5) {
 
            //... O.K. l and lnext should be good 2 consectvie hits
            l->zStatus(20);
            lnext->zStatus(20);
            l->zPair(iforLine + 1);
            lnext->zPair(iforLine);
            }
#ifdef TRKRECO_DEBUG_DETAIL
            else {
            std::cout << "    ... rejected because delta z > 1.5";
            std::cout << std::endl;
            }
#endif
        } 
#ifdef TRKRECO_DEBUG_DETAIL
        else {
            if (err1) {
            std::cout << "    ... s-z calculation error with ";
            std::cout << l->wire()->name() << std::endl;
            }
            if (err2) {
            std::cout << "    ... s-z calculation error with ";
            std::cout << lnext->wire()->name() << std::endl;
            }
        }
#endif
        }
    }
 
    //... Calculate s and z...
    //... Aleady solved
        if(l->zStatus() == 20){
        TMLink * t = new TMLink(* l);
        t->zStatus(l->zStatus());
        t->zPair(l->zPair());
        int err = track.szPosition(* t);
        if (err) {
        delete t;
        continue;
        }
        else {
        //...Store the sz link...
        t->link(l);
        forLine.append(t);
        //... Check # of hits in a wire layer and Clustering them
        // nlyr[layertoStlayer(l->wire()->layerId())] += 1;
        nlyr[l->wire()->layer()->axialStereoLayerId()] += 1;
        // llyr[layertoStlayer(l->wire()->layerId())]  = iforLine;
        llyr[l->wire()->layer()->axialStereoLayerId()]  = iforLine;
        iforLine += 1;
        }
    }
    else {
        //...Assuming wire position...
        //... for left
        TMLink * tl = new TMLink(* l);
        tl->leftRight(WireHitLeft);
        tl->zStatus(-10);
        tl->zPair(0);
        int err = track.szPosition(* tl);
        if(err){
        delete tl;
        }
        else {
        //...Store the sz link...
        tl->link(l);
        forLine.append(tl);
        //... Check # of hits in a wire layer and Clustering them
        // nlyr[layertoStlayer(l->wire()->layerId())] += 1;
        nlyr[l->wire()->layer()->axialStereoLayerId()] += 1;
        // llyr[layertoStlayer(l->wire()->layerId())]  = iforLine;
        llyr[l->wire()->layer()->axialStereoLayerId()]  = iforLine;
        //... 
        iforLine += 1;
        }
 
        //... for right
        TMLink * tr = new TMLink(* l);
        tr->leftRight(WireHitRight);
        tr->zStatus(-10);
        tr->zPair(0);
        err = track.szPosition(* tr);
        if(err) {
        delete tr;
        }
        else{
        //...Store the sz link...
        tr->link(l);
        forLine.append(tr);
        //... Check # of hits in a wire layer and Clustering them
        // nlyr[layertoStlayer(l->wire()->layerId())] += 1;
        nlyr[l->wire()->layer()->axialStereoLayerId()] += 1;
        // llyr[layertoStlayer(l->wire()->layerId())]  = iforLine;
        llyr[l->wire()->layer()->axialStereoLayerId()]  = iforLine;
        iforLine += 1;
        }
    }
    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... dump of sz links" << std::endl;
    Dump(forLine, "sort flag", "    ");
    for (unsigned i = 0; i < 18; i++) {
    std::cout << i << " : " << nlyr[i] << ", " << llyr[i] << std::endl;
    }
#endif
 
    //...Check # of sz links...
    if (forLine.length() < _lineSelector.nLinksStereo()) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... rejected by sz nLinks(";
    std::cout << forLine.length() << ") < ";
    std::cout << _lineSelector.nLinks() << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    return NULL;
    }
 
    //... Select the best segment in each Superlayer
    AList<TMLink> gdSLine0;
    AList<TMLink> gdSLine1;
    AList<TMLink> gdSLine2;
    AList<TMLink> gdSLine3;
    AList<TMLink> gdSLine4;
    double min_chi2[5] = {9999.,9999.,9999.,9999.,9999.};
    double min_a[5] = {9999.,9999.,9999.,9999.,9999.};
 
#ifdef TRKRECO_DEBUG_DETAIL
//      TTrackBase base;
//      base.append(forLine);
//      TWindow baseWindow("s-z window");
//      baseWindow.append(base);
 
    bool display = false;
//      for (unsigned i = 0; i < forLine.length(); i++) {
//      if (forLine[i]->wire()->name() == "15=31")
//          display = true;
//      }
#endif
 
    //...For all stereo super layers...
    for (unsigned isl=0; isl < 5; isl++) {
    AList<TMLink> tmpLine;
    AList<TMLink> goodLine;
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... stereo supter layer " << isl << std::endl;
#endif
 
    //...3-stereo-layer case...
    if (isl < 2) {
 
        //...Initialize...
        unsigned ily1 = TMDC::getTMDC("1.0")->superLayer(isl * 2 + 1)
        ->first()->axialStereoLayerId();
        unsigned ily2 = ily1 + 1; 
        unsigned ily3 = ily2 + 1; 
 
        //...Loop for the first layer...
        unsigned ilst = llyr[ily1] + 1;
        unsigned ifst = ilst - nlyr[ily1];
        for (unsigned i =ifst; i<ilst; i++) {
        TMLink & l = * forLine[i];
        tmpLine.append(l);
        TMLink & m = *forLine[l.zPair()];           
        if (l.zStatus() == 20) {
            if (l.zPair() < i ) {
            tmpLine.append(m);
            }
            else {
            //...Already considered...
            tmpLine.remove(l);
            continue;
            }
        }
 
        //...Loop for the second layer...
        unsigned jlst = llyr[ily2] + 1;
        unsigned jfst = jlst - nlyr[ily2];
        for (unsigned j =jfst; j< jlst; j++) {
            TMLink & l2 = * forLine[j];
            tmpLine.append(l2);
            TMLink & m2 = *forLine[l2.zPair()];           
            if (l2.zStatus() == 20) {
            if (l2.zPair() < j) {
                tmpLine.append(m2);
            }
            else {
                //...Already considered.
                tmpLine.remove(l2);
                continue;
            }
            }
 
            //...Loop for the third layer...
            unsigned klst = llyr[ily3] + 1;
            unsigned kfst = klst - nlyr[ily3];
            for (unsigned k=kfst; k < klst; k++) {
            TMLink & l3 = * forLine[k];
            tmpLine.append(l3);
            TMLink & m3 = * forLine[l3.zPair()];           
            if (l3.zStatus() == 20) {
                if (l3.zPair() < k) {
                tmpLine.append(m3);
                }
                else {
                //...Already considered.
                tmpLine.remove(l3);
                continue;
                }
            }
 
            //... Check the hits in neighbor wirelayer
            //...
            //...    x|o|o|x   |
            //...      |o|     |
            //...              V IP ,  o means OK. 
            //...
            int relation12 = -1;
            int relation23 = -1;
 
            //... Check the hit in fist and it in the second layer
            relation12 = checkHits(l.hit()->wire()->localId(),
                           l2.hit()->wire()->localId(),
                           isl);
            //... 
            if (l.zStatus() == 20 && relation12 < 0)
                relation12 = checkHits(m.hit()->wire()->localId(),
                           l2.hit()->wire()->localId(),
                           isl);
 
            //...Check the hit in second  and it in the third layer
            relation23 = checkHits(l2.hit()->wire()->localId(),
                           l3.hit()->wire()->localId(),
                           isl);
            if (l.zStatus() == 20 && relation23 < 0)
                relation23 = checkHits(m2.hit()->wire()->localId(),
                           l3.hit()->wire()->localId(),
                           isl);
         
            //...Bad relation...
            if(relation12 || relation23 ) {
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "    ... bad relations";
                std::cout << " : segment rejected";
                Dump(tmpLine, "detail stereo", "    ");
                if (display) {
                TLine0 line(tmpLine);
                int err = line.fit();
//                  baseWindow.append(line);
//                  baseWindow.draw();
//                  baseWindow.wait();
//                  baseWindow.remove(line);
                }
#endif
 
                tmpLine.remove(l3);
                if (l3.zStatus() == 20) tmpLine.remove(m3);
                continue;
            }
 
            //...Make a segment
            unsigned ntmp = tmpLine.length();
            if(ntmp < 3) { //...Is this needed???...
                std::cout << "!!! is this possible !!!???" << std::endl;
 
                tmpLine.remove(l3);
                if(l3.zStatus() == 20) tmpLine.remove(m3);
                continue;
            }
 
            //...Do the linear fit 
            TLine0 line(tmpLine); 
            int err = line.fit();
            if (err < 0) {
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "    ... line fit error";
                std::cout << " : segment rejected" << std::endl;
                Dump(line.links(), "detail stereo", "    ");
                if (display) {
//                  baseWindow.append(line);
//                  baseWindow.draw();
//                  baseWindow.wait();
//                  baseWindow.remove(line);
                }
#endif
 
                tmpLine.remove(l3);
                if(l3.zStatus() == 20) tmpLine.remove(m3);
                continue;
            }
 
            //...Check quality for the fit
            //... if |z0| > 70 or chi2 > 0.7, remove it
            //              double chi2 = line.chi2()/(ntmp-2); 
            double chi2 = line.reducedChi2(); 
            if (fabs(line.b()) < _stereoZ3 &&
                chi2 < _stereoChisq3 &&
                chi2 < min_chi2[isl]) {
                goodLine = tmpLine;
                min_chi2[isl] = chi2;
                min_a[isl] = line.a();
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "    ... segment accepted : " << std::endl;
#endif
            }
#ifdef TRKRECO_DEBUG_DETAIL
            else {
                std::cout << "    ... bad quality : segment rejected :";
                std::cout << "chi2=" << chi2;
                std::cout << ", b=" << line.b() << std::endl;
            }
            Dump(line.links(), "detail stereo", "    ");
            if (display) {
//                  baseWindow.append(line);
//                  baseWindow.draw();
//                  baseWindow.wait();
//                  baseWindow.remove(line);
            }
#endif
 
            //...end of third loop
            if(l3.zStatus() ==20) tmpLine.remove(m3);
            tmpLine.remove(l3);
            }
            //... end of second loop
            if(l2.zStatus() ==20) tmpLine.remove(m2);
            tmpLine.remove(l2);
        }
        //... end of first loop
        if(l.zStatus() ==20) tmpLine.remove(m);
        tmpLine.remove(l);
        }
 
        //... Keep best segments
        if(isl==0) gdSLine0 = goodLine;
        if(isl==1) gdSLine1 = goodLine;
        goodLine.removeAll();
 
    }
    else {
 
        //... Nlayer == 4
 
        //...Initialize
        // unsigned ily1 = firstStlayer(isl);
        unsigned ily1 = TMDC::getTMDC("1.0")->superLayer(isl * 2 + 1)
        ->first()->axialStereoLayerId();
        unsigned ily2 = ily1 + 1; 
        unsigned ily3 = ily2 + 1; 
        unsigned ily4 = ily3 + 1; 
        if (nlyr[ily1] == 0 ||
        nlyr[ily2] == 0 ||
        nlyr[ily3] == 0 ||
        nlyr[ily4] == 0) continue;
 
        //...loop for the first layer
        unsigned ilst = llyr[ily1] + 1;
        unsigned ifst = ilst - nlyr[ily1];
        for(unsigned i =ifst; i<ilst;i++ ){
        TMLink & l = * forLine[i];
        tmpLine.append(l);
        TMLink & m = *forLine[l.zPair()];           
        if(l.zStatus() == 20 ) {
            if(l.zPair()< i) {
            tmpLine.append(m);
            } else {
            tmpLine.remove(l);
            continue;
            }
        }
 
        //...loop for the second layer
        unsigned jlst = llyr[ily2] + 1;
        unsigned jfst = jlst - nlyr[ily2];
        for(unsigned j =jfst; j< jlst; j++){
            TMLink & l2 = * forLine[j];
            tmpLine.append(l2);
            TMLink & m2 = *forLine[l2.zPair()];           
            if(l2.zStatus() == 20 ) {
            if(l2.zPair()< j) {
                tmpLine.append(m2);
            } else {
                tmpLine.remove(l2);
                continue;
            }
            }
            
            //...loop for the third layer
            unsigned klst = llyr[ily3] + 1;
            unsigned kfst = klst - nlyr[ily3];
            for( unsigned k=kfst; k < klst; k++ ){
            TMLink & l3 = * forLine[k];
            tmpLine.append(l3);
            TMLink & m3 = *forLine[l3.zPair()];           
            if(l3.zStatus() == 20 ) {
                if(l3.zPair()< k) {
                tmpLine.append(m3);
                } else {
                tmpLine.remove(l3);
                continue;
                }
            }
            
            //...loop for the 4th layer
            unsigned hlst = llyr[ily4] + 1;
            unsigned hfst = hlst - nlyr[ily4];
            for( unsigned h=hfst; h < hlst; h++ ){
                TMLink & l4 = * forLine[h];
                tmpLine.append(l4);
                TMLink & m4 = *forLine[l4.zPair()];           
                if(l4.zStatus() == 20 ) {
                if(l4.zPair()< h) {
                    tmpLine.append(m4);
                } else {
                    tmpLine.remove(l4);
                    continue;
                }
                }
 
                //...Check the relation between the hit
                //   in neighbor wirelayer
                //...
                //...    x|o|o|x      |
                //...      |o|        |
                //...                 V IP ,   o means OK. 
                //...
                int relation12 = -1;
                int relation23 = -1;
                int relation34 = -1;
 
//...For debug...
//              if (l.hit()->wire()->consective(* l2.hit()->wire())
//...For debug end...
 
                
                //... Check the hit in fist and it in the second
                relation12 = checkHits(l.hit()->wire()->localId(),
                           l2.hit()->wire()->localId(),
                           isl);
                if (l.zStatus() == 20 && relation12 < 0)
                relation12 =
                    checkHits(m.hit()->wire()->localId(),
                          l2.hit()->wire()->localId(),
                          isl);
                
                //... Check the hit in second  and it in the third
                relation23 = checkHits(l2.hit()->wire()->localId(),
                           l3.hit()->wire()->localId(),
                           isl);
                if (l.zStatus() == 20 && relation23 < 0)
                relation23 =
                    checkHits(m2.hit()->wire()->localId(),
                          l3.hit()->wire()->localId(),
                          isl);
                
                //... Check the hit in second  and it in the forth
                relation34 = checkHits(l3.hit()->wire()->localId(),
                           l4.hit()->wire()->localId(),
                           isl);
                if (l3.zStatus() == 20 && relation34 < 0)
                relation34 =
                    checkHits(m3.hit()->wire()->localId(),
                          l4.hit()->wire()->localId(),
                          isl);
                
                //...remove Bad segments
                if(relation12 || relation23 || relation34) {
 
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "    ... bad relations";
                std::cout << " : segment rejected";
                Dump(tmpLine, "detail stereo", "    ");
                if (display) {
                TLine0 line(tmpLine);
                int err = line.fit();
//                  baseWindow.append(line);
//                  baseWindow.draw();
//                  baseWindow.wait();
//                  baseWindow.remove(line);
                }
#endif
 
                tmpLine.remove(l4);
                if(l4.zStatus() == 20) tmpLine.remove(m4);
                continue;
                }
 
                //...Make a segment
                unsigned ntmp = tmpLine.length();
                if(ntmp < 4) {
                tmpLine.remove(l4);
                if(l4.zStatus() == 20) tmpLine.remove(m4);
                continue;
                }
 
                //...Do the linear fit
                TLine0 line(tmpLine); 
                int err = line.fit();
                if (err < 0) {
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "    ... line fit error";
                std::cout << " : segment rejected" << std::endl;
                Dump(line.links(), "detail stereo", "    ");
                if (display) {
//                  baseWindow.append(line);
//                  baseWindow.draw();
//                  baseWindow.wait();
//                  baseWindow.remove(line);
                }
#endif
 
                tmpLine.remove(l4);
                if(l4.zStatus() == 20) tmpLine.remove(m4);
                continue;
                }
 
                //...Check qualit of the fit
                //... if |z0| > 20 or chi2 > 0.5, remove it
                double chi2 = line.reducedChi2(); 
                if (fabs(line.b()) < _stereoZ4 &&
                chi2 < _stereoChisq4 &&
                chi2 < min_chi2[isl]) {
                    
                //...Keep good segments
                goodLine = tmpLine;
                min_chi2[isl] = chi2;
                min_a[isl] = line.a();
                    
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "    segment accepted : " << std::endl;
#endif
                }
#ifdef TRKRECO_DEBUG_DETAIL
            else {
                std::cout << "    ... bad quality : segment rejected :";
                std::cout << " chi2=" << chi2;
                std::cout << ", b=" << line.b() << std::endl;
            }
            Dump(line.links(), "detail stereo", "    ");
            if (display) {
//                  baseWindow.append(line);
//                  baseWindow.draw();
//                  baseWindow.wait();
//                  baseWindow.remove(line);
            }
#endif
 
                //...end of the 4th loop
                if(l4.zStatus() ==20) tmpLine.remove(m4);
                tmpLine.remove(l4);
            }
            //...end of the third loop 
            if(l3.zStatus() ==20) tmpLine.remove(m3);
            tmpLine.remove(l3);
            }
            //...end of the second loop
            if(l2.zStatus() ==20) tmpLine.remove(m2);
            tmpLine.remove(l2);
        }
        //... end of the first loop
        if(l.zStatus() ==20) tmpLine.remove(m);
        tmpLine.remove(l);
        }
 
        //...
        if(isl==2) gdSLine2 = goodLine;
        if(isl==3) gdSLine3 = goodLine;
        if(isl==4) gdSLine4 = goodLine;
        goodLine.removeAll();
    }
    } 
 
 
    //... Link segments
 
    //...Check how many segments are made
    unsigned Nsgmnts[5] = {0,0,0,0,0};
    Nsgmnts[0] = gdSLine0.length();
    Nsgmnts[1] = gdSLine1.length();
    Nsgmnts[2] = gdSLine2.length();
    Nsgmnts[3] = gdSLine3.length();
    Nsgmnts[4] = gdSLine4.length();
 
    unsigned NusedSgmnts = 0;
    for(unsigned jsl = 0; jsl < 5; jsl++) {
    if(Nsgmnts[jsl] > 0) NusedSgmnts +=1;
    }
 
    //...Require at least one Segment
    if (NusedSgmnts == 0) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... rejected because no segment found" << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    return NULL;
    }
 
    //... Make a line
    AList<TMLink> forNewLine;
    forNewLine.append(gdSLine0);
    forNewLine.append(gdSLine1);
    forNewLine.append(gdSLine2);
    forNewLine.append(gdSLine3);
    forNewLine.append(gdSLine4);
 
    //...
    unsigned nNewLine = forNewLine.length();
    float R = fabs(track.helix().curv());
    TLine0 newLine(forNewLine);
    
    //...Do linear fit
    int err = newLine.fit();
    //    double newLine_chi2 = newLine.chi2()/(nNewLine - 2);
    double newLine_chi2 = newLine.reducedChi2();
    double newLine_a = newLine.a();
 
    //...Check the quality of the line.
 
    //...If chi2 > 0.25 for R > 80.( > 0.8 for R < 80.), refine the line.
    //    if(((R > 80. && newLine_chi2 > 0.25) ||(R < 80. && newLine_chi2 > 0.8))&& NusedSgmnts > 1){
    if(((R > 80. && newLine_chi2 > 4.0) ||(R < 80. && newLine_chi2 > 13.0))&& NusedSgmnts > 1){
    //...Look at difference between the slope of the line and that of segment
    double max_diff_a = 0.;
    unsigned this_sly = 999;
    for(unsigned isl=0; isl < 5; isl++){
        if(Nsgmnts[isl]==0) continue;
        double diff_a = fabs((min_a[isl]-newLine_a)/newLine_a);
        if(diff_a > max_diff_a){
        max_diff_a = diff_a;
        this_sly = isl;
        }
    }
    
    //...If max slope diff. > 0.4 for R < 50.(> 0.3 for R < 50), remove it.
    if((R< 50. && max_diff_a > 0.4) || (R > 50. && max_diff_a > 0.3)) {
        
        //...clear # of entries in the segement
        Nsgmnts[this_sly]=0;
        
        //... remove the worst setment
        if(this_sly == 0){
        newLine.removeSLY(gdSLine0);
        } else if (this_sly == 1){
        newLine.removeSLY(gdSLine1);
        } else if (this_sly == 2){
        newLine.removeSLY(gdSLine2);
        } else if (this_sly == 3){
        newLine.removeSLY(gdSLine3);
        } else if (this_sly == 4){
        newLine.removeSLY(gdSLine4);
        } 
        
        //... fit again
        unsigned NnewLine_2 = newLine.links().length();
        if(NnewLine_2 < 3) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "... rejected because of few links for line" << std::endl;
#endif  
        HepAListDeleteAll(forLine);
        return NULL;
        }
        
        int err = newLine.fit();
        if (err < 0) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "... rejected because of line fit failure" << std::endl;
#endif  
        HepAListDeleteAll(forLine);
        return NULL;
        }
        double newLine_chi2_2 = newLine.chi2()/NnewLine_2;
    }
    }
    
    //...Remove bad points...
    AList<TMLink> bad;
    double maxSigma = 1.0;
    if(R < 80) maxSigma = 1.5;
    newLine.refine(bad, maxSigma);
    if(newLine.links().length()< 2) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... rejected because of few links for line after refine";
    std::cout << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    return NULL;
    }
    err = newLine.fit();
    if (err < 0) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... rejected because of line fit failure(2)" << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    return NULL;
    }
 
    //...Append hits, if the distance between the line and hits < 
    for(unsigned isl = 0; isl < 5; isl++){
      if(Nsgmnts[isl] == 0) {
        double maxdist = 0.5;
    if(R < 80) maxdist = 1.25;
        newLine.appendByszdistance(forLine, 2*isl+1, maxdist);
      }
    }
    err = newLine.fit();
    if (err < 0) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... rejected because of line fit failure(3)" << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    return NULL;
    }
 
    //...Remove bad points again...
    maxSigma = 1.0;
    newLine.refine(bad, maxSigma);
    if(newLine.links().length()< 2) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... rejected because of few links for line after 2nd refine";
    std::cout << std::endl;
#endif  
    HepAListDeleteAll(forLine);
    return NULL;
    }
 
    //...Linear fit again...
    err = newLine.fit();
    if (err < 0) {
    HepAListDeleteAll(forLine);
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    appendStereo cut ... new line 2nd linear fit failure. ";
    std::cout << "# of links = " << n << ",";
    std::cout << "," << nNewLine << std::endl;
#endif  
    return NULL;
    }
    
    //...
     unsigned NnewLine_3 = newLine.links().length();
     double newLine_chi2_3 = newLine.chi2()/NnewLine_3;
 
    //... Do we need quality cut?
    //     if(newLine_chi2_3 > 10.) return NULL;
 
    //...3D fit...
    const AList<TMLink> & good = newLine.links();
    unsigned nn = good.length();
    for (unsigned i = 0; i < nn; i++) {
    track.append(* good[i]->link());
    }
 
    //...Set initial values
    Vector a(5);
    a = track.helix().a();
    a[3] = newLine.b();
    a[4] = track.charge() * newLine.a();
    track._helix->a(a);
 
    //...Refine...
    err = fit(track);
    if (err < 0) goto discard;
    track.refine(bad, _stereoMaxSigma * 100.);
    err = fit(track);
    if (err < 0) goto discard;
    track.refine(bad, _stereoMaxSigma * 10.);
    err = fit(track);
    if (err < 0) goto discard;
    track.refine(bad, _stereoMaxSigma);
    err = fit(track);
    if (err < 0) goto discard;
 
    //...Test it...
    if (! _trackSelector.select(track)) goto discard;
 
    //...Termination...
    HepAListDeleteAll(forLine);
    return & track;
 
    //...Something happened...
discard:
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... rejected by fitting failure : ";
    std::cout << " err = " << err << std::endl;
#endif
    HepAListDeleteAll(forLine);
    return NULL;
}
 
const TMSelector &
TBuilder0::trackSelector(const TMSelector & a) {
    _circleSelector.nLinks(a.nLinks());
    _circleSelector.nSuperLayers(a.nSuperLayers());
    _circleSelector.minPt(a.minPt());
    _circleSelector.maxImpact(a.maxImpact());
 
    _trackSelector.nLinks(a.nLinks());
    _trackSelector.nSuperLayers(a.nSuperLayers());
    _trackSelector.maxSigma(a.maxSigma());
 
    _lineSelector.nLinksStereo(a.nLinksStereo());
    _lineSelector.maxDistance(a.maxDistance());
 
    return a;
}
 
void
TBuilder0::appendClusters(TTrack & track,
             const AList<TMLink> & list) const {
 
    AList<TMLink> tmp = list;
 
    //...Append them...
    track.appendByApproach(tmp, _trackSelector.maxSigma() * 2. / 3.);
 
    //...Refine it...
    AList<TMLink> bad;
    fit(track);
    track.refine(bad, _trackSelector.maxSigma());
}
 
int
TBuilder0::checkHits(unsigned i, unsigned j, unsigned isl) const {
 
     int nWr[5] = {80,128,160,208,256};
     int ilast = nWr[isl]-1;
     int ilocal = (int) i;
     int jlocal = (int) j;
     //...
     if(ilocal > 0 && ilocal < ilast){
       if(fabs(float(jlocal-ilocal)) > 1 ) {
         return -1;
       } else {
         return 0;
       }
     }else if(ilocal == 0){
       if(jlocal > 1 && jlocal < ilast){
         return -1;
       }else {
         if(jlocal == ilast){
         return 0;
       } else if(jlocal == 0){
         return 0;
       } else if(jlocal == 1){
         return 0;
       } else {
         return -1;
       }
       }
     }else if(ilocal == ilast){
       if(jlocal > 0 && jlocal < ilast-1){
         return -1;
       }else {
         if(jlocal == ilast-1){
         return 0;
       } else if(jlocal == ilast){
         return 0;
       } else if(jlocal == 0){
         return 0;
       } else {
         return -1;
       }
       }
     }
     //...
     return -1;
}
 
void
TBuilder0::salvage(TTrack & t, AList<TMLink> & hits) const {
 
//    if (t.type() == TrackTypeNormal) {
    salvageNormal(t, hits);   //Pt high enough; fromIP.
    return;
//    }
  //tmply for cosmic ray...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << name() << " ... salvaging" << std::endl;
    std::cout << "    # of given hits=" << hits.length() << std::endl;
#endif
 
    unsigned nHits = hits.length();
    if (nHits == 0) return;
 
    //...Hit loop...
    AList<TMLink> candidates;
    for (unsigned i = 0; i < nHits; i++) {
    TMLink & l = * hits[i];
    const TMDCWireHit & h = * l.hit();
 
    //...Already used?...
    if (h.state() & WireHitUsed) continue;
    candidates.append(l);
    }
 
    //...Try to append this hit...
//    t.appendByApproach(candidates, 10.);
    t.appendByApproach(candidates, _salvageLevel);
    fit(t);
    hits.remove(candidates);
    t.assign(WireHitConformalFinder);
    t.finder(TrackOldConformalFinder);
    // t.assign(WireHitConformalFinder, TrackOldConformalFinder);
}
 
void
TBuilder0::salvageNormal(TTrack & t, AList<TMLink> & hits) const {
    
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << name() << " ... normal salvaging : ";
    std::cout << " # of hits=" << hits.length() << std::endl;
#endif
 
    unsigned nHits = hits.length();
    if (nHits == 0) return;
 
    //...Determine direction to salvage...
    const HepPoint3D & center = t.helix().center();
    const HepVector3D a = ORIGIN - center;    
 
    //...Hit loop...
    AList<TMLink> candidates;
    for (unsigned i = 0; i < nHits; i++) {
    TMLink & l = * hits[i];
    const TMDCWireHit & h = * l.hit();
    if (a.cross(h.xyPosition()).z() * t.charge() > 0.) continue;
 
    //...Already used?...
    if (h.state() & WireHitUsed) continue;
    candidates.append(l);
    }
 
    //...Try to append this hit...
    t.appendByApproach(candidates, 30.);
    fit(t);
    hits.remove(candidates);
    t.assign(WireHitConformalFinder);
    t.finder(TrackOldConformalFinder);
    // t.assign(WireHitConformalFinder, TrackOldConformalFinder);
}
 
int
TBuilder0::check2CnHits(TMLink &l, TMLink & s, int ichg) const {
 
    //...Check same layer ?...
    if(l.hit()->wire()->layerId() != s.hit()->wire()->layerId()) return -1;
 
    //...Initialization...
    int nsl[11] = {64,80,96,128,144,160,192,208,240,256,288};
    float hcell[50] = {0.,0.,0.,0.,0.,0.,0.687499,0.747198,0.806896,0.,0.,0.,0.,0.,0.,0.782967,0.820598,0.858229,0.,0.,0.,0.,0.,0.878423,0.908646,0.939845,0.970068,0.,0.,0.,0.,0.,0.892908,0.916188,0.940219,0.963499,0.,0.,0.,0.,0.,0.901912,0.920841,0.940382,0.959312,0.,0.,0.,0.,0.};
 
     int ilast  = nsl[l.hit()->wire()->superLayerId()] - 1;
     int ilocal = l.hit()->wire()->localId();
     int jlocal = s.hit()->wire()->localId();
 
     double ddist1 = l.hit()->drift() / hcell[l.hit()->wire()->layerId()];
     double ddist2 = s.hit()->drift() / hcell[s.hit()->wire()->layerId()];
     double ddist = 0.5 * (ddist1 + ddist2);
 
     //...Case by case...
     if (ilocal > 0 && ilocal < ilast) {
     if (abs(jlocal - ilocal) > 1) {
         return -20;
     }
     else {
         if (ddist > 0.65  &&
         ((ddist1 > 0.7 && ddist2 > 0.7) ||
          (ddist1 > 0.95 || ddist2 >0.95))) {
 
         //...O.K. 2 consective hits
         if (ichg > 0) {
             return 1;
         }
         else {
             return 0;
         }
         }
         else {
         return -20;
         }
     }
     }
     else if (ilocal == 0) {
     if (jlocal > 1 && jlocal < ilast) {
         return -20;
     }
     else {
         if (ddist > 0.65 &&
         ((ddist1 > 0.7 && ddist2 > 0.7) ||
          (ddist1 > 0.95 || ddist2 > 0.95))) {
         if (jlocal == ilast) {
 
             //...O.K. 2 consective hits
             if(ichg > 0){
             return 0;
             }
             else {
             return 1;
             }
         }
         else if (jlocal == 1) {
             if (ichg > 0) {
             return 1;
             }
             else {
             return 0;
             }
         }
         }
     }
     }
     else if (ilocal == ilast) {
     if(jlocal > 0 && jlocal < ilast-1) {
         return -20;
     }
     else {
         if (ddist > 0.65 &&
         ((ddist1 > 0.7 && ddist2 > 0.7) ||
          (ddist1 > 0.95 || ddist2 > 0.95))) {
         if (jlocal == ilast-1) {
             //...O.K. 2 consective hits
             if(ichg > 0){
             return 0;
             } else {
             return 1;
             }
         } else if(jlocal == 0){
             if(ichg > 0){
             return 1;
             } else {
             return 0;
             }
         }
         } else {
         return -20;
         }
     }
     }
     
     //...fails
     return -50;
}
 
int
TBuilder0::consectiveHits(TMLink & l, TMLink & s, int ichg) const {
 
    //...Diagonal length of a cell...
    static float hcell[50] = {0.,0.,0.,0.,0.,0.,
                  0.687499,0.747198,0.806896,
                  0.,0.,0.,0.,0.,0.,
                  0.782967,0.820598,0.858229,
                  0.,0.,0.,0.,0.,
                  0.878423,0.908646,0.939845,0.970068,
                  0.,0.,0.,0.,0.,
                  0.892908,0.916188,0.940219,0.963499,
                  0.,0.,0.,0.,0.,
                  0.901912,0.920841,0.940382,0.959312,
                  0.,0.,0.,0.,0.};
    const TMDCWire & wire = * l.hit()->wire();
    const TMDCWire & next = * s.hit()->wire();
 
    //...Check same layer ?...
    if (wire.layerId() != next.layerId()) return -1;
 
    //...Are these consective ?...
    if (! wire.consective(next)) return -20;
 
    //...Drift distance...
    int ilast  = wire.layer()->nWires() - 1;
    unsigned layerId = wire.layerId();
    double ddist1 = l.hit()->drift() / hcell[layerId];
    double ddist2 = s.hit()->drift() / hcell[layerId];
    double ddist = 0.5 * (ddist1 + ddist2);
 
    //...Distance check...
    if (ddist > 0.65  &&
    ((ddist1 > 0.7 && ddist2 > 0.7) || (ddist1 > 0.95 || ddist2 > 0.95))) {
        
    //...O.K. 2 consective hits
    if (ichg > 0) return 1;
    else          return 0;
    }
    else {
    return -50;
    }
}