TBuilderCosmic.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TBuilderCosmic.cxx,v 1.6 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TBuilderCosmic.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to build a cosmic track.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#ifdef TRKRECO_DEBUG_DETAIL
#  ifndef TRKRECO_DEBUG
#    define TRKRECO_DEBUG
#  endif
#endif
#include "TrkReco/TBuilderCosmic.h"
#include "TrkReco/T3DLine.h"
#include "TrkReco/TLine0.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrack.h"
 
TBuilderCosmic::TBuilderCosmic( const std::string& name, float salvageLevel )
    : TBuilder0( name, salvageLevel ), _fitter( "TBuilderCosmic Fitter" ) {}
 
TBuilderCosmic::~TBuilderCosmic() {}
 
TTrack* TBuilderCosmic::buildStereo( TTrack& track, const AList<TMLink>& list ) const {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << name() << "(stereo) ... dump of stereo candidate hits" << std::endl;
  AList<TMLink> tmp = list;
  tmp.sort( SortByWireId );
  std::cout << "    ";
  for ( unsigned i = 0; i < tmp.length(); i++ )
  {
    TMLink* l = tmp[i];
    std::cout << l->wire()->layerId() << "-";
    std::cout << l->wire()->localId() << ",";
  }
  std::cout << std::endl;
#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];
 
    //... Require Fitting vaildation
    if ( !( l->hit()->state() & WireHitFittingValid ) ) continue;
 
    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;
  std::cout << "    ";
  tmp = forLine;
  tmp.sort( SortByWireId );
  for ( unsigned i = 0; i < tmp.length(); i++ )
  {
    TMLink* l = tmp[i];
    std::cout << l->wire()->layerId() << "-";
    std::cout << l->wire()->localId() << ",";
  }
  std::cout << std::endl;
#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.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 );
 
  //...Make a seed track again
  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.);  //Liuqg, meaningless while without magnetic field
  //    err = newLine.fit();
  //    newLine.refine(bad, 20.);
  //    err = newLine.fit();
  //    newLine.refine(bad, 10.);
  //    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[3] = newLine.b();
  a[4] = track.charge() * newLine.a();
  track._helix->a( a );
 
#ifdef LINE_COSMIC
  T3DLine* Ltrack = new T3DLine( track );
 
  //...Refine...
  err = _fitter.fit( *Ltrack );
  Ltrack->refine( bad, _trackSelector.maxSigma() * 30. );
  err = _fitter.fit( *Ltrack );
  Ltrack->refine( bad, _trackSelector.maxSigma() * 3. );
  err = _fitter.fit( *Ltrack );
  //    Ltrack->refine(bad, _trackSelector.maxSigma() * 0.21);  //liuqg, for prelimilary test
  track.refine( bad, _trackSelector.maxSigma() );
  err = _fitter.fit( *Ltrack );
 
  //...Test it...
  /*    if (! _trackSelector.select(*Ltrack)) {
          HepAListDeleteAll(forLine);
          HepAListDeleteAll(forNewLine);
          delete Ltrack;
          return NULL;
      }
  */
  //...Termination...
  HepAListDeleteAll( forLine );
  HepAListDeleteAll( forNewLine );
  track.removeLinks();
  track.append( Ltrack->links() );
  Vector a1( 5 );
  a1 = Ltrack->helix().a();
  track._helix->a( a1 );
  delete Ltrack;
  return &track;
#endif
 
  //...Refine...
  err = _fitter.fit( track );
  track.refine( bad, _trackSelector.maxSigma() * 30. );
  err = _fitter.fit( track );
  track.refine( bad, _trackSelector.maxSigma() * 3. );
  err = _fitter.fit( track );
  track.refine( bad, _trackSelector.maxSigma() ); // liuqg, for prelimilary test
  //    track.refine(bad, _trackSelector.maxSigma() * 0.21);
  err = _fitter.fit( track );
  //...Test it...
  if ( !_trackSelector.select( track ) )
  {
    HepAListDeleteAll( forLine );
    HepAListDeleteAll( forNewLine );
    return NULL;
  }
 
  //...Termination...
  HepAListDeleteAll( forLine );
  HepAListDeleteAll( forNewLine );
  return &track;
}