MdcSegFinder.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: MdcSegFinder.cxx,v 1.17 2012/10/18 08:35:38 zhangy Exp $
//
// Description:
//
//
// Environment:
//      Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Authors:
//  Steve Schaffner
//  Zhang Yao(zhangyao@ihep.ac.cn) Migration for BESIII
//  2010-04-13
//  Zhang Yao(zhangyao@ihep.ac.cn) Keep best segment for 4-hit pattern
//------------------------------------------------------------------------
 
#include "MdcTrkRecon/MdcSegFinder.h"
#include "MdcData/MdcHit.h"
#include "MdcData/MdcHitMap.h"
#include "MdcData/MdcHitUse.h"
#include "MdcGeom/BesAngle.h"
#include "MdcGeom/Constants.h"
#include "MdcGeom/MdcDetector.h"
#include "MdcGeom/MdcSuperLayer.h"
#include "MdcTrkRecon/GmsList.h"
#include "MdcTrkRecon/MdcLine.h"
#include "MdcTrkRecon/MdcMap.h"
#include "MdcTrkRecon/MdcSeg.h"
#include "MdcTrkRecon/MdcSegList.h"
#include "MdcTrkRecon/MdcSegUsage.h"
#include "MdcTrkRecon/countBits.h"
#include "MdcTrkRecon/mdcWrapAng.h"
#include "MdcTrkRecon/mdcWrapWire.h"
// yzhang hist cut
#include "GaudiKernel/NTuple.h"
#include <map>
 
extern NTuple::Tuple*       g_tupleFindSeg;
extern NTuple::Item<double> g_findSegChi2;
extern NTuple::Item<double> g_findSegIntercept;
extern NTuple::Item<double> g_findSegSlope;
extern NTuple::Item<double> g_findSegChi2Refit;
extern NTuple::Item<double> g_findSegChi2Add;
extern NTuple::Item<int>    g_findSegPat;
extern NTuple::Item<int>    g_findSegNhit;
extern NTuple::Item<int>    g_findSegPat34;
extern NTuple::Item<int>    g_findSegSl;
extern NTuple::Item<double> g_findSegMc;
extern NTuple::Item<double> g_findSegAmbig;
extern int                  haveDigiTk[43][288];
extern int                  haveDigiAmbig[43][288];
// zhangy hist cut
 
using std::cout;
using std::endl;
 
//-------------------------------------------------------------------------
MdcSegFinder::MdcSegFinder( int useAllAmbig ) : patternList( useAllAmbig ) {
  //-------------------------------------------------------------------------
}
//-------------------------------------------------------------------------
int MdcSegFinder::createSegs( const MdcDetector* gm, MdcSegList& segs,
                              const MdcMap<const MdcHit*, MdcSegUsage*>& usedHits,
                              const MdcHitMap* map, double bunchTime ) {
  //-------------------------------------------------------------------------
  // Forms all possible Segs of hits to act as candidate tracks
  // Construct group of 8 hits; translate hit wires into bit-mapped word
  // and look up valid patterns (if any) for this group.
 
  _addHits           = segs.segPar()->addHits;
  int          nSegs = 0;
  int          newSegs;
  unsigned int groupWord;
  MdcHit*      groupHits[8];
  int          lPrevHit = 0; // flags that 2nd layer was hit in last group
  // which layer/wire to look at for each hit in group
  static const int layerWire[8][2] = { { 0, -1 }, { 0, 0 },  { 1, 0 }, { 2, -1 },
                                       { 2, 0 },  { 3, -1 }, { 3, 0 }, { 3, 1 } };
 
  // Loop through the superlayers
  const MdcSuperLayer* slayer = gm->firstSlay();
  while ( slayer != 0 )
  {
    const MdcLayer* layArray[4];
    int             wireIndex[4];
 
    //    layArray[0] = slayer->firstLayer();
    //    for (int i = 1; i < 4; i++) layArray[i] = gm->nextLayer(layArray[i-1]);
    int nslay = slayer->nLayers();
    for ( int i = 0; i < nslay; i++ ) layArray[i] = slayer->layer( i );
    if ( nslay != 4 ) layArray[3] = 0;
 
    //    std::cout<<"-------------super layer----- "<<slayer->Id() << " nlayer
    //    "<<nslay<<std::endl;
    // for (int i = 0; i < nslay; i++) std::cerr<<layArray[i]->Id()<<"
    // "<<layArray[i]->nWires()<<"\n"; slayer = slayer->next();continue;
 
    // Loop over cells (index refers to cells in 2nd layer of superlayer)
    for ( int cellIndex = 0; cellIndex < layArray[1]->nWires(); cellIndex++ )
    {
      // yzhang change FIXME
      double phi = layArray[1]->getWire( cellIndex )->phiE();
      for ( int i = 0; i < 4; i++ )
      {
        wireIndex[i] = cellIndex;
        if ( slayer->slayNum() > 4 ) continue; // nWires over No.4 slayer is uniform
        if ( ( slayer->slayNum() > 9 ) && ( i == 3 ) ) break; // stop when get last layer
        if ( i == 1 ) continue;
        if ( i == 3 ) phi = layArray[2]->getWire( wireIndex[2] )->phiE();
        BesAngle tmp( phi - layArray[i]->phiEPOffset() );
        int      wlow = (int)floor( layArray[i]->nWires() * tmp.rad() / twopi );
        int      wbig = (int)ceil( layArray[i]->nWires() * tmp.rad() / twopi );
        if ( tmp == 0 )
        {
          wlow = -1;
          wbig = 1;
        }
        if ( i != 3 ) { wireIndex[i] = wbig; }
        else { wireIndex[i] = wlow; }
        // zhangy change
        wireIndex[i] = mdcWrapWire( wireIndex[i], layArray[i]->nWires() );
      }
      // Form group of 8 wires
      groupWord        = 0u;
      unsigned bitWord = 1u;
      int      nGroup  = 0;
      for ( int ii = 0; ii < 8; ii++ )
      {
        groupHits[ii] = 0;
        //-----------------------------
        if ( layArray[layerWire[ii][0]] == 0 ) continue;
        //-----------------------------
        int laynum = layArray[layerWire[ii][0]]->layNum();
        int wire   = wireIndex[layerWire[ii][0]] + layerWire[ii][1];
        wire = mdcWrapWire( wire, layArray[layerWire[ii][0]]->nWires() ); // handle wrap-around
        MdcHit* thisHit = map->hitWire( laynum, wire );
        if ( thisHit != 0 )
        {
          if ( !usedHits.get( thisHit )->deadHit() )
          {
            groupHits[ii] = thisHit;
            groupWord |= bitWord;
            nGroup++;
          }
          else
          {
            //      cout << "Skipping hit " << endl;
          }
        }
        // Quit if no hits in 1st 2 layers or 1 hit in 3 layers
        if ( ( ii == 2 && nGroup < 1 ) || ( ii == 4 && nGroup < 2 ) ) break;
        bitWord <<= 1;
      }
      if ( nGroup < 3 ) continue;
 
      // int lastWire = mdcWrapWire(cellIndex - 1, layArray[1]->nWires());//2011-08-10
      int lastWire = mdcWrapWire( cellIndex - 1, layArray[0]->nWires() );
      if ( map->hitWire( layArray[1]->layNum(), lastWire ) != 0 ) lPrevHit = 1;
      //      if (layArray[1]->hitWire(lastWire) != 0) lPrevHit = 1;
      else lPrevHit = 0;
 
      // Try all allowed 4-hit patterns for this group (if any)
      int nPatt = patternList.npatt4[groupWord];
 
      if ( ( layArray[1]->layNum() == 5 ) && ( cellIndex == 46 ) )
      {
        for ( int i = 0; i < 4; i++ )
        {
          // std::cout<<__FILE__<<" "<<__LINE__<<"====("<<layArray[i]->layNum()<<","<<
          // wireIndex[i]<<")" << std::endl;
        }
        // std::cout<<__FILE__<<" "<<__LINE__<< " groupWord:"<<groupWord<<" nPatt:"<<nPatt<<
        // std::endl;
      }
 
      if ( nPatt > 0 )
      {
        newSegs = tryPatterns( groupHits, groupWord, 4, lPrevHit, nPatt,
                               patternList.allowedPatt4[groupWord], slayer, segs, usedHits,
                               map, bunchTime );
        nSegs += newSegs;
      }
 
      // If unsuccessful, try 3-hit patterns ???? may want to try 2nd pass here
      if ( !segs.segPar()->find3[slayer->index()] ) continue;
      int nUnused = 0;
      for ( int i = 0; i < 8; i++ )
      {
        if ( groupHits[i] != 0 )
        {
          if ( usedHits.get( groupHits[i] )->usedSeg() == 0 ) nUnused++;
        }
      }
      if ( nUnused > 0 )
      {
        nPatt = patternList.npatt3[groupWord];
        if ( nPatt > 0 )
        {
 
          newSegs = tryPatterns( groupHits, groupWord, 3, lPrevHit, nPatt,
                                 patternList.allowedPatt3[groupWord], slayer, segs, usedHits,
                                 map, bunchTime );
          nSegs += newSegs;
        }
      }
    } // cellIndex loop
 
    slayer = slayer->next();
  } // superlayer loop
 
  if ( nSegs > 0 )
  {
    segs.massageSegs();
    segs.resetSeed( gm );
  }
 
  nSegs = segs.count();
 
  if ( 5 == segs.segPar()->lPrint ) { cout << "Number of Segs found: " << nSegs << "\n"; }
 
  return nSegs;
}
 
//---------------------------------------------------------------------------
int MdcSegFinder::tryPatterns( MdcHit* groupHits[8], unsigned groupWord, int nInPatt,
                               int lPrevHit, int npatt, int* allowedPat,
                               const MdcSuperLayer* slay, MdcSegList& segs,
                               const MdcMap<const MdcHit*, MdcSegUsage*>& usedHits,
                               const MdcHitMap* map, double t0 ) {
  //---------------------------------------------------------------------------
  int nSegs = 0;
  int nbit  = 8;
 
  unsigned* patterns;
  int**     ambigPatt;
  if ( nInPatt == 3 )
  {
    patterns  = patternList.patt3;
    ambigPatt = patternList.ambigPatt3;
  }
  else
  {
    patterns  = patternList.patt4;
    ambigPatt = patternList.ambigPatt4;
  }
 
  MdcSeg* trialSeg = new MdcSeg( t0 ); // start w/ 1 seg active
 
  // Create temporary array of hit structures for fitting segment
  MdcLine span( 12 );
  int     spanAmbig[12];
  MdcHit* spanHits[12];
  MdcHit* ahit;
 
  // *** Loop over the allowed pattern for this group
  for ( int iPatt = 0; iPatt < npatt; iPatt++ )
  {
    unsigned thisPat  = patterns[allowedPat[iPatt]];
    unsigned testWord = thisPat & groupWord;
 
    if ( countbits( testWord ) < nInPatt ) continue;
    if ( lPrevHit && nInPatt == 3 && ( thisPat == 051 || thisPat == 0111 ) ) continue;
 
    // *** Load the hits into hit structure (radius only at this point)
    // radius is rel. to avg. for superlayer; resolution is converted to delphi
    unsigned bitadd = 1u;
    int      nhit   = 0;
    for ( int ibit = 0; ibit < nbit; ibit++ )
    {
      // Is this cell in the pattern?
      if ( bitadd & thisPat )
      {
        const MdcLayer* layer = groupHits[ibit]->layer();
        if ( layer == 0 ) cout << "huh?" << endl;
        span.x[nhit]   = layer->rMid() - slay->rad0();
        spanHits[nhit] = groupHits[ibit];
        nhit++;
        if ( nhit == nInPatt ) break;
      }
      bitadd <<= 1;
    }
 
    // *** Loop through all ambiguity combinations; the loop counter also
    //    serves as a bit-mapped pattern of ambiguities
    //    Discard combo if all hits have been used w/ this ambig before
 
    std::map<int, MdcSeg*> bestTrialSegs; // yzhang 2012-09-18
 
    int namb = 1 << nInPatt;
    // std::cout<< __FILE__ << "   " << __LINE__ << " namb  "<<namb<<" nInPatt
    // "<<nInPatt<<std::endl;
    for ( int iamb = 0; iamb < namb; iamb++ )
    {
 
      // Skip all but allowed ambiguity patterns
      if ( ambigPatt[allowedPat[iPatt]][iamb] != 1 ) continue;
 
      // Convert loop counter into physical ambiguity (0=> -1; 1=>+1)
      int nUsed = 0;
      int ihit;
      for ( ihit = 0; ihit < nInPatt; ihit++ )
      {
        spanAmbig[ihit] = ( ( 1u << ihit ) & iamb ) ? 1 : -1;
        nUsed += usedHits.get( spanHits[ihit] )->usedAmbig( spanAmbig[ihit] );
      }
      if ( nUsed >= nInPatt ) continue;
 
      // ***Calculate phi for each hit, correcting approx. for angle of track
      /* First calculate a correction factor for the drift distance (since
         we're treating hits as being at radius of wire).  This may slow
         things down.*/
      /* corr = 1/cos(theta), where theta = del(d)/del(r), taken over the
         span.  Use 1/cos = sqrt(1+tan**2) for calculating correction.  */
 
      double rIn   = slay->firstLayer()->rMid();
      double rOut  = slay->lastLayer()->rMid();
      double phiIn = mdcWrapAng( spanHits[nInPatt - 1]->phi(), spanHits[0]->phi() );
      double dphidr =
          ( ( spanHits[nInPatt - 1]->phi() +
              spanAmbig[nhit - 1] *
                  spanHits[nInPatt - 1]->driftDist( t0, spanAmbig[nhit - 1] ) / rOut ) -
            ( phiIn + spanAmbig[0] * spanHits[0]->driftDist( t0, spanAmbig[0] ) / rIn ) ) /
          ( rOut - rIn ); // yzhang temp FIXME
 
      double corr = sqrt( 1 + slay->rad0() * slay->rad0() * dphidr * dphidr );
 
      if ( 5 == segs.segPar()->lPrint )
      {
        std::cout << __FILE__ << " " << __LINE__ << " corr" << corr << " phi_n "
                  << spanHits[nInPatt - 1]->phi() << " dd_n "
                  << spanAmbig[nhit - 1] *
                         spanHits[nInPatt - 1]->driftDist( t0, spanAmbig[nhit - 1] )
                  << " phiIn " << phiIn << " dd_in "
                  << spanAmbig[0] * spanHits[0]->driftDist( t0, spanAmbig[0] ) << std::endl;
      }
      // yzhang add fix phi accross -x axis
      bool   crossAxis = false;
      double yOld      = 0;
      // zhangy add
      double sigmaSum = 0.;
      double driftSum = 0.;
      double driftHit[12];
      // Actual phi calc
      for ( ihit = 0; ihit < nInPatt; ihit++ )
      {
        ahit         = spanHits[ihit];
        double rinv  = 1. / ahit->layer()->rMid();
        double drift = ahit->driftDist( t0, spanAmbig[ihit] );
        span.y[ihit] =
            ahit->phi() + spanAmbig[ihit] * drift * ( corr * rinv ); // yzhang temp FIXME
 
        if ( 5 == segs.segPar()->lPrint )
        {
          std::cout << " in segment finding (" << ahit->layernumber() << ","
                    << ahit->wirenumber() << ")"
                    << " |driftDist| " << fabs( drift ) << " ambig " << spanAmbig[ihit]
                    << " corr " << corr << " rinv " << rinv << " sigma "
                    << ahit->sigma( drift, spanAmbig[ihit] ) << std::endl;
        }
        // span.sigma[ihit] = ahit->sigma(fabs(drift), spanAmbig[ihit]);
        // yzhang 2010-05-20 , FIXME
        span.sigma[ihit] = ahit->sigma( fabs( drift ), spanAmbig[ihit] ) * ( corr * rinv );
 
        // yzhang add temp FIXME
        sigmaSum += span.sigma[ihit];
        driftSum += drift;
        driftHit[ihit] = drift;
        // zhangy add temp FIXME
 
        // yzhang add fix phi accross -x axis,set flag
        if ( ( span.y[ihit] != 0 ) && ( !crossAxis ) )
        {
          if ( ( yOld / span.y[ihit] ) < 0 ) crossAxis = true;
          yOld = span.y[ihit];
        }
        // zhangy add
      }
 
      // yzhang add temp FIXME
      // for (ihit = 0; ihit < nInPatt; ihit++) {
      // span.sigma[ihit] = sigmaSum/nInPatt*
      //( fabs(driftHit[ihit]-(driftSum/nInPatt))/(driftSum/nInPatt) );
      // }
      // zhangy add temp FIXME
 
      //--yzhang add, fix phi accross -x axis , if cross +twopi
      if ( crossAxis )
      {
        for ( int ihit = 0; ihit < nInPatt; ihit++ )
        {
          // Approx,for max phi of a cell is 2*pi/40, max cross 4 cell=0.628<0.7
          if ( abs( span.y[ihit] ) < 0.7 ) break;
          if ( span.y[ihit] < 0 ) span.y[ihit] += twopi;
        }
      }
      // zhangy add
 
      // Fit the segment
      if ( 5 == segs.segPar()->lPrint )
        std::cout << " first fit(" << nInPatt << ")" << std::endl;
      span.fit( nInPatt );
      BesAngle temp    = span.intercept;
      span.intercept   = temp.posRad();
      double t_segchi2 = span.chisq / ( nInPatt - 2 );
 
      // yzhang test 2011-06-14
      if ( 5 == segs.segPar()->lPrint )
      {
        for ( int ii = 0; ii < nInPatt; ii++ )
        {
          std::cout << __FILE__ << " " << __LINE__ << " " << ii << " span.x " << setw( 10 )
                    << span.x[ii] << " y " << setw( 10 ) << span.y[ii] << " sigma "
                    << setw( 10 ) << span.sigma[ii] << std::endl;
        }
      }
      if ( g_tupleFindSeg != NULL )
      {
        g_findSegSlope     = span.slope;
        g_findSegIntercept = span.intercept;
        g_findSegChi2      = t_segchi2;
      }
      // zhangy test
 
      if ( 5 == segs.segPar()->lPrint )
      {
        std::cout << __FILE__ << " " << __LINE__ << " pattern " << thisPat << " nHit "
                  << nInPatt << " chi2/dof " << t_segchi2 << " chisq " << span.chisq
                  << " maxChisq=" << segs.segPar()->maxChisq << std::endl;
        for ( int jj = 0; jj < nInPatt; jj++ )
        {
          std::cout << "resid[" << jj << "] " << span.resid[jj] << " sigma " << span.sigma[jj]
                    << " chisq add "
                    << span.resid[jj] * span.resid[jj] / ( span.sigma[jj] * span.sigma[jj] )
                    << std::endl; // yzhang debug
        }
      }
      if ( span.chisq / ( nInPatt - 2 ) > segs.segPar()->maxChisq )
      {
        if ( 5 == segs.segPar()->lPrint )
        {
          std::cout << __FILE__ << " " << __LINE__ << "CUT! chi2/(N-2) "
                    << span.chisq / ( nInPatt - 2 ) << " > " << segs.segPar()->maxChisq
                    << std::endl;
          for ( std::map<int, MdcSeg*>::iterator iter = bestTrialSegs.begin();
                iter != bestTrialSegs.end(); iter++ )
          { cout << " bestTrialSeg nHit=" << iter->first << endl; }
        }
        continue;
      }
 
      if ( span.slope > Constants::pi || span.slope < ( -1. ) * Constants::pi )
      {
        if ( 5 == segs.segPar()->lPrint )
          std::cout << " CUT by span.|slope| " << span.slope << " > pi" << std::endl;
        continue; // yzhang add 2012-08-07
      }
 
      int nInSeg = nInPatt;
 
      // *** Pick up adjacent hits and refit
      // Recalculate correction factor (optional)
      if ( segs.segPar()->segRefit )
      {
        dphidr    = span.slope;
        corr      = sqrt( 1 + slay->rad0() * slay->rad0() * dphidr * dphidr );
        crossAxis = false;
        for ( ihit = 0; ihit < nInSeg; ihit++ )
        {
          ahit         = spanHits[ihit];
          double rinv  = 1. / ahit->layer()->rMid();
          double drift = ahit->driftDist( t0, spanAmbig[ihit] );
          span.y[ihit] =
              ahit->phi() + spanAmbig[ihit] * drift * ( corr * rinv ); // yzhang temp FIXME
          // yzhang add
          // span.y[ihit] = ahit->phi();
          // zhangy add
 
          // yzhang 2010-05-20 , FIXME
          // span.sigma[ihit] = ahit->sigma(fabs(drift), spanAmbig[ihit]);
          span.sigma[ihit] = ahit->sigma( fabs( drift ), spanAmbig[ihit] ) * ( corr * rinv );
          // yzhang add fix phi accross -x axis,set flag
          if ( ( span.y[ihit] != 0 ) && ( !crossAxis ) )
          {
            if ( ( yOld / span.y[ihit] ) < 0 ) crossAxis = true;
            yOld = span.y[ihit];
          }
          // zhangy add
        }
        // yzhang add, fix phi accross -x axis , if cross +twopi
        if ( crossAxis )
        {
          for ( int ihit = 0; ihit < nInPatt; ihit++ )
          {
            // Approx,for max phi of a cell is 2*pi/40, max cross 4 cell=0.628<0.7
            if ( abs( span.y[ihit] ) < 0.7 ) break;
            if ( span.y[ihit] < 0 ) span.y[ihit] += twopi;
          }
        }
        if ( 5 == segs.segPar()->lPrint )
          std::cout << " second fit(" << nInSeg << ")" << std::endl;
        // zhangy add
        span.fit( nInSeg );
        BesAngle temp  = span.intercept;
        span.intercept = temp.posRad();
 
        if ( span.chisq / ( nInPatt - 2 ) > segs.segPar()->maxChisq )
        {
          if ( 5 == segs.segPar()->lPrint )
          {
            std::cout << __FILE__ << " " << __LINE__ << "CUT! chi2/(N-2) "
                      << span.chisq / ( nInPatt - 2 ) << " > " << segs.segPar()->maxChisq
                      << std::endl;
            for ( std::map<int, MdcSeg*>::iterator iter = bestTrialSegs.begin();
                  iter != bestTrialSegs.end(); iter++ )
            { cout << " bestTrialSeg nHit=" << iter->first << endl; }
          }
          continue;
        }
      }
 
      if ( g_tupleFindSeg != NULL ) { g_findSegChi2Refit = span.chisq / ( nInSeg - 2 ); }
 
      trialSeg->setValues( nInPatt, nInSeg, spanHits, &span, slay, spanAmbig );
 
      if ( 5 == segs.segPar()->lPrint )
      {
        std::cout << " try:  " << std::endl;
        trialSeg->plotSegAll(); // yzhang temp 2012-09-18
        std::cout << "   " << std::endl;
      }
      if ( _addHits )
      {
        // Look for adjacent hits
        nInSeg = trialSeg->addHits( &span, spanHits, map, corr );
      }
 
      if ( 5 == segs.segPar()->lPrint )
      {
        for ( std::map<int, MdcSeg*>::iterator iter = bestTrialSegs.begin();
              iter != bestTrialSegs.end(); iter++ )
        {
          cout << " bestTrialSeg nHit=" << iter->first << " chi2 " << iter->second->chisq()
               << endl;
        }
        std::cout << "trialSeg chisq " << trialSeg->chisq()
                  << std::endl; //<<" best chi2 " <<bestTrialSeg->chisq() <<  std::endl;
        cout << "segment phi: " << trialSeg->phi() << "  slope: " << trialSeg->slope()
             << "  nhit: " << trialSeg->nHit() << " chi2 " << trialSeg->chisq() << endl;
        // cout << "layer, wire, ambig, used, drift: " << endl;
        for ( int i = 0; i < trialSeg->nHit(); i++ )
        {
          int           ambi   = trialSeg->hit( i )->ambig();
          const MdcHit* theHit = trialSeg->hit( i )->mdcHit();
          theHit->print( cout );
          cout << " ambig " << ambi;
        }
        cout << endl;
      }
      trialSeg->setPattern( iPatt );
      trialSeg->markHits( usedHits );
      if ( nInPatt == 4 ) trialSeg->setFull();
 
      if ( g_tupleFindSeg != NULL )
      {
        double t_mcRadio     = -1.;
        double t_nAmbigRight = 0.;
        int    t_tkId        = -1;
        for ( int ii = 0; ii < trialSeg->nHit(); ii++ )
        {
          unsigned int l = trialSeg->hit( ii )->mdcHit()->layernumber();
          unsigned int w = trialSeg->hit( ii )->mdcHit()->wirenumber();
          if ( haveDigiTk[l][w] < 0 || haveDigiTk[l][w] > 100 ) continue;
          if ( t_tkId == -1 ) { t_tkId = haveDigiTk[l][w]; }
          else if ( haveDigiTk[l][w] != t_tkId )
          {
            t_mcRadio = -999; // hits in this seg not in same mc track
          }
 
          // test ambig right ratio
          if ( haveDigiAmbig[l][w] == trialSeg->hit( ii )->ambig() ) t_nAmbigRight++;
        }
 
        double t_nSame = 0.;
        for ( int ii = 0; ii < trialSeg->nHit(); ii++ )
        {
          unsigned int l = trialSeg->hit( ii )->mdcHit()->layernumber();
          unsigned int w = trialSeg->hit( ii )->mdcHit()->wirenumber();
          if ( haveDigiTk[l][w] == t_tkId ) { ++t_nSame; }
          if ( t_mcRadio > -998. ) t_mcRadio = t_nSame / nInSeg;
        }
        g_findSegPat = iPatt;
        if ( 3 == nInPatt ) g_findSegPat34 = 3;
        else g_findSegPat34 = 4;
        g_findSegSl      = slay->slayNum();
        g_findSegMc      = t_mcRadio;
        g_findSegAmbig   = t_nAmbigRight / nInSeg;
        g_findSegChi2Add = span.chisq / ( nInSeg - 2 );
        g_findSegNhit    = nInSeg;
        g_tupleFindSeg->write();
      }
 
      //--keep only one segment for ambig
      if ( 3 == nInPatt )
      {
        segs.append( trialSeg ); // Add to list of Segs
        nSegs++;
      }
      else
      {
        // debug
        // cout<<" before insert nBest="<<bestTrialSegs.size()<<endl;
        // for(std::map<int,MdcSeg*>::iterator iter = bestTrialSegs.begin(); iter !=
        // bestTrialSegs.end(); iter++) {
        //   cout<<"after insert bestTrialSeg nHit="<<iter->first<<endl;
        //   iter->second->plotSeg();
        // }
        std::map<int, MdcSeg*>::iterator iter = bestTrialSegs.find( trialSeg->nHit() );
        if ( iter == bestTrialSegs.end() )
        {
          bestTrialSegs.insert(
              std::map<int, MdcSeg*>::value_type( trialSeg->nHit(), trialSeg ) );
          if ( 5 == segs.segPar()->lPrint )
          {
            std::cout << "   ----insert  " << trialSeg->nHit() << std::endl;
            trialSeg->plotSegAll();
            std::cout << " \n  " << std::endl;
          }
        }
        else
        {
          if ( trialSeg->chisq() < iter->second->chisq() )
          {
            MdcSeg* tempSeg = iter->second;
            delete tempSeg;
            iter->second = trialSeg;
            if ( 5 == segs.segPar()->lPrint )
            {
              std::cout << " ----update  " << iter->first << std::endl;
              trialSeg->plotSegAll();
              std::cout << " \n  " << std::endl;
            }
          }
          else
          {
            if ( 5 == segs.segPar()->lPrint )
            {
              std::cout << "DROP TrialSeg  " << std::endl;
              trialSeg->plotSegAll();
              std::cout << " \n  " << std::endl;
            }
            delete trialSeg;
          }
        }
        ////debug
        // cout<<" after insert nBest="<<bestTrialSegs.size()<<endl;
        // for(std::map<int,MdcSeg*>::iterator iter = bestTrialSegs.begin(); iter !=
        // bestTrialSegs.end(); iter++) {
        //   cout<<"after insert bestTrialSeg nHit="<<iter->first<<endl;
        //   iter->second->plotSeg();
        // }
        // if((bestTrialSeg->nHit() == trialSeg->nHit()) && (bestTrialSeg->chisq() >
        // trialSeg->chisq())){
        //   if(5 == segs.segPar()->lPrint) std::cout<< "KEEP trial as best. orignal
        //   bestTrialSeg  "
        //     << bestTrialSeg->chisq()<<"/ndof "<<bestTrialSeg->nHit()-2<< ">
        //     trialSeg->chisq()  "
        //       <<trialSeg->chisq()<<"/ndof "<<trialSeg->nHit()-2<<std::endl;
        //   delete bestTrialSeg;
        //   bestTrialSeg = trialSeg;
        //   }
      }
 
      //--mark best trialSeg
      // std::cout<< __LINE__<<" best " <<bestTrialSeg->chisq()<< " trial  " <<
      // trialSeg->chisq() <<std::endl;
 
      trialSeg = new MdcSeg( t0 );
    } // end ambiguity loop
 
    //--keep only one segment/Hit for ambig
    for ( std::map<int, MdcSeg*>::iterator iter = bestTrialSegs.begin();
          iter != bestTrialSegs.end(); iter++ )
    {
      segs.append( iter->second ); // Add to list of Segs
      if ( 5 == segs.segPar()->lPrint )
      {
        std::cout << " append bestTrialSeg of nHit " << iter->first << std::endl;
        iter->second->plotSeg();
        std::cout << std::endl;
      }
    }
 
    bestTrialSegs.clear(); // empty bestTrialSegs
 
    nSegs++;
 
  } // end pattern loop
 
  delete trialSeg;
  return nSegs;
}