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    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// 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/TTrack.h"
#include "TrkReco/TTrackHEP.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 = nullptr; // add initialize 25-05-15
    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 = -1;// add initialize 25-05-15
 
  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;
}