MdcSegGrouperSt.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: MdcSegGrouperSt.cxx,v 1.14 2011/09/26 01:06:37 zhangy Exp $
//
// Description:
//
//
// Environment:
//      Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Authors:
//  Zhang Yao(zhangyao@ihep.ac.cn)  Migrate to BESIII
//
// Copyright (C)  1996  The Board of Trustees of
// The Leland Stanford Junior University.  All Rights Reserved.
//------------------------------------------------------------------------
#include "MdcTrkRecon/MdcSegGrouperSt.h"
#include "CLHEP/Alist/AList.h"
#include "MdcData/MdcHit.h"
#include "MdcGeom/BesAngle.h"
#include "MdcGeom/MdcDetector.h"
#include "MdcGeom/MdcSuperLayer.h"
#include "MdcTrkRecon/GmsList.h"
#include "MdcTrkRecon/MdcSeg.h"
#include "MdcTrkRecon/MdcSegInfoSterO.h"
#include "MdcTrkRecon/MdcSegList.h"
#include "MdcTrkRecon/MdcSegUsage.h"
#include "MdcTrkRecon/MdcSegWorks.h"
#include "MdcTrkRecon/MdcTrack.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrkBase/TrkFit.h"
#include "TrkBase/TrkRecoTrk.h"
#include "TrkFitter/TrkHelixMaker.h"
 
// yzhang hist cut
#include "AIDA/IHistogram1D.h"
extern AIDA::IHistogram1D* g_z0Cut;
extern AIDA::IHistogram1D* g_ctCut;
extern AIDA::IHistogram1D* g_delCt;
extern AIDA::IHistogram1D* g_delZ0;
// zhangy hist cut
 
// Constructor
//------------------------------------------------------------------------
MdcSegGrouperSt::MdcSegGrouperSt( const MdcDetector* gm, int debug )
    : MdcSegGrouper( gm, gm->nStereoSuper( -1 ) + gm->nStereoSuper( 1 ), debug ) {
  //------------------------------------------------------------------------
 
  lTestGroup  = true;
  lTestSingle = false;
 
  isValid = new bool*[nPly()];
  for ( int j = 0; j < nPly(); j++ ) { isValid[j] = 0; }
}
 
//------------------------------------------------------------------------
void MdcSegGrouperSt::resetList() {
  //------------------------------------------------------------------------
  // Clear existing lists for stereo finding
  int nsuper = _gm->nSuper();
  for ( int i = 0; i < nsuper; i++ ) { segList[i].removeAll(); }
}
 
//------------------------------------------------------------------------
void MdcSegGrouperSt::fillWithSegs( const MdcSegList* inSegs, const MdcTrack* track ) {
  //------------------------------------------------------------------------
  // Prepare for stereo finding
  // Load segments consistent with input track (already in phi order)
  //  Assuming from origin
  int nsuper = _gm->nSuper();
  resetList(); // Clear existing lists
  if ( 3 == _debug )
    std::cout << std::endl << "=====MdcSegGrouperSt::fillWithSegs=====" << std::endl;
 
  const TrkFit* tkFit = track->track().fitResult();
  if ( tkFit == 0 ) return;
  TrkExchangePar par  = tkFit->helix( 0.0 );
  double         kap  = 0.5 * par.omega();
  double         phi0 = par.phi0();
  double         d0   = par.d0();
  MdcSegWorks    segStuff; // holds some calculated values to pass to segInfo
 
  bool lStraight = false;
  if ( fabs( kap ) < 1.e-9 ) { lStraight = true; }
 
  double rCurl = 99999999.;
  if ( !lStraight ) { rCurl = fabs( d0 + 1. / kap ); }
 
  // Create an info object to store the info (will be owned by seg)
  MdcSegInfoSterO* info = new MdcSegInfoSterO;
 
  // Loop over superlayers
  for ( int isuper = 0; isuper < nsuper; isuper++ )
  {
    const GmsList* inList = inSegs->oneList( isuper );
 
    if ( inList->count() == 0 ) continue;
 
    MdcSeg* inSeg = (MdcSeg*)inList->first();
    // Only load stereo segments
    if ( inSeg->superlayer()->whichView() == 0 ) continue;
 
    // Calculate r & phi (approx) of axial track at slayer
    double radius = inSeg->superlayer()->rEnd();
    if ( radius >= rCurl ) break;
    double phiArg = kap * radius;
    if ( lStraight ) phiArg = d0 / radius;
    if ( phiArg < -1.0 ) phiArg = -1.0;
    else if ( phiArg > 1.0 ) phiArg = 1.0;
    segStuff.phiArg = phiArg;
    segStuff.phiAx.setRad( phi0 + asin( phiArg ) ); // Approx??!!!!
    BesAngle delPhi( fabs( inSeg->superlayer()->delPhi() ) );
    // BesAngle maxPhiA = segStuff.phiAx + delPhi + 0.05;  // allow a little slop
    // BesAngle minPhiA = segStuff.phiAx - delPhi - 0.05;
    BesAngle maxPhiA = segStuff.phiAx + delPhi + 0.1; // allow a little slop yzhang 2011-06-15
    BesAngle minPhiA = segStuff.phiAx - delPhi - 0.1; // yzhang TEMP 2011-06-15
    // yzhang FIXME
    double maxPhi  = maxPhiA;
    double minPhi  = minPhiA;
    bool   phitest = ( maxPhi > minPhi );
    // phitest = false => the valid range straddles phi = 0, so skip
    // checking against phimin and phimax
 
    // Calc some things needed in segInfo::calcStereo
    segStuff.wirLen2inv = 1. / ( inSeg->superlayer()->zEnd() * inSeg->superlayer()->zEnd() );
    // Loop over segs in superlayer
    if ( _debug == 3 )
      std::cout << std::endl
                << "Pick segment by phi from " << minPhi << " to " << maxPhi
                << " phiAx=" << segStuff.phiAx << " delPhi=" << delPhi << std::endl;
    for ( ; inSeg != 0; inSeg = (MdcSeg*)inSeg->next() )
    {
      // Test if within allowed phi range
      BesAngle phiSeg( inSeg->phi() );
 
      if ( _debug == 3 ) { inSeg->plotSeg(); }
      if ( phitest )
      {
        if ( phiSeg < minPhi )
        {
          if ( _debug == 3 )
          {
            std::cout << " CUT by phi " << phiSeg << "<min " << minPhi << std::endl;
          }         // yzhang debug
          continue; // not up to candidates
        }
        else if ( phiSeg > maxPhi )
        {
          if ( _debug == 3 )
          {
            std::cout << " CUT by phi " << phiSeg << ">max " << maxPhi << std::endl;
          }      // yzhang debug
          break; // past candidates
        }
      }
      else
      { // !phitest
        if ( phiSeg > maxPhi && phiSeg < minPhi )
        {
          if ( _debug == 3 )
          {
            std::cout << "!phitest " << phiSeg << " max " << maxPhi << " min " << minPhi
                      << std::endl;
          } // yzhang debug
          continue;
        }
      }
 
      if ( _debug == 3 ) std::cout << " **KEEP seg phi=" << phiSeg << std::endl;
 
      // std::cout<< __FILE__ << "   " << __LINE__ << " MdcSegGrouperSt::fillWithSegs call
      // calcStereo "<<std::endl;
      int isBad = info->calcStereo( inSeg, track->track(), segStuff );
 
      if ( isBad )
      {
        if ( _debug == 3 ) { std::cout << " CUT by calcStereo isBad!" << std::endl; }
        continue;
      }
      // Test for sensible values of ct and z0
      if ( g_z0Cut ) { g_z0Cut->fill( info->z0() ); }
      if ( g_ctCut ) { g_ctCut->fill( info->ct() ); }
      if ( fabs( info->ct() ) > inSeg->segPar()->ctcut )
      {
        if ( _debug == 3 )
        {
          std::cout << " CUT by ctcut! " << fabs( info->ct() )
                    << " > cut:" << inSeg->segPar()->ctcut << " " << phiSeg << std::endl;
        }
        continue;
      }
      if ( fabs( info->z0() ) > inSeg->segPar()->z0cut )
      {
        if ( _debug == 3 )
        {
          std::cout << " CUT by z0cut! " << fabs( info->z0() ) << " >cut "
                    << inSeg->segPar()->z0cut << " " << phiSeg << std::endl;
        }
        continue;
      }
      inSeg->setInfo( info );
      info = new MdcSegInfoSterO;
      segList[isuper].append( inSeg );
      // if(_debug==3)std::cout<<" APPEND this stereo seg"<<  std::endl;
 
    } // end loop over new segs
  }   //  end loop over superlayers
  delete info;
}
 
//-------------------------------------------------------------------------
void MdcSegGrouperSt::plotStereo() const {
  //-------------------------------------------------------------------------
  int nsuper = _gm->nSuper();
  int isuper;
  /*
     for (isuper = 0; isuper < nsuper; isuper++) {
     if (segList[isuper].length() == 0) continue;
 
     MdcSeg *inSeg = (MdcSeg *) segList[isuper].first();
  // Only draw stereo segments
  if (inSeg->superlayer()->whichView() == 0) continue;
  for (int j = 0 ; j < (int) segList[isuper].length(); j++) {
  segList[isuper][j]->plotSeg(0,orange);
  }
  }
  */
  for ( isuper = 0; isuper < nsuper; isuper++ )
  {
    if ( segList[isuper].length() == 0 ) continue;
 
    MdcSeg* inSeg = (MdcSeg*)segList[isuper].first();
    // Only draw stereo segments
    // if (inSeg->superlayer()->whichView() == 0) continue;
    for ( int j = 0; j < (int)segList[isuper].length(); j++ )
    { segList[isuper][j]->plotSeg(); }
  }
}
//-------------------------------------------------------------------------
int MdcSegGrouperSt::incompWithSeg( const MdcSeg* refSeg, const MdcSeg* testSeg ) {
  //-------------------------------------------------------------------------
 
  return 0;
}
 
//-------------------------------------------------------------------------
int MdcSegGrouperSt::incompWithGroup( MdcSeg** segGroup, const MdcSeg* testSeg, int iply ) {
  //-------------------------------------------------------------------------
  // Test that the latest seg lies more or less in a line with the others
  //  Currently requiring line to point to IP; also require rough
  //  agreement in ct.
  // iply = depth index of current seg (not yet in group)
 
  int i;
  // Find first segment already in group
  for ( i = iply - 1; i > 0; i-- )
  {
    if ( !leaveGap[i] ) break;
  }
  const MdcSeg* first = ( *combList[i] )[currentSeg[i]];
 
  // Test ct
  MdcSegInfoSterO* firstInfo = (MdcSegInfoSterO*)first->info();
  MdcSegInfoSterO* newInfo   = (MdcSegInfoSterO*)testSeg->info();
 
  if ( g_delCt ) { g_delCt->fill( firstInfo->ct() - newInfo->ct() ); } // yzhang hist cut
  if ( fabs( firstInfo->ct() - newInfo->ct() ) > testSeg->segPar()->delCtCut )
  {
    if ( _debug == 3 )
    {
      cout << "---MdcSegGrouperSt Ct CUT!" << endl; // yzhang debug
      std::cout << "first:";
      first->plotSeg();
      std::cout << std::endl;
      std::cout << "test:";
      testSeg->plotSeg();
      std::cout << std::endl;
      std::cout << "first.ct " << firstInfo->ct() << " test.ct " << newInfo->ct()
                << " delta ct " << firstInfo->ct() - newInfo->ct() << " Cut "
                << testSeg->segPar()->delCtCut << std::endl; // yzhang debug
 
      std::cout << "--- " << std::endl;
    }
    return -1;
  }
  else {}
 
  double arcFirst = firstInfo->arc();
  double arcNew   = newInfo->arc();
  double zFirst   = firstInfo->z0() + firstInfo->ct() * arcFirst;
  double zNew     = newInfo->z0() + newInfo->ct() * arcNew;
  // project line from IP through 1st seg to new seg
  double zProj = zFirst / arcFirst * arcNew;
  if ( g_delZ0 ) { g_delZ0->fill( zProj - zNew ); } // yzhang hist cut
 
  if ( fabs( zProj - zNew ) > testSeg->segPar()->delZ0Cut )
  {
    if ( 3 == _debug )
    {
      cout << "MdcSegGrouperSt delZ0Cut  not incompWithGroup CUT!" << endl; // yzhang debug
      testSeg->plotSeg();
      std::cout << " zProj " << zProj << " zNew " << zNew << " CUT! "
                << testSeg->segPar()->delZ0Cut << std::endl;
    }
    return -1;
  }
 
  /*
     double delZ = newInfo->z0() + newInfo->ct() * arcNew - zFirst;
     double z0 = zFirst - arcFirst * delZ / (arcNew - arcFirst);
     if (fabs(z0) > testSeg->segPar()->delZ0Cut)return -1;
 
     for (i = iply - 1; i > 0; i--) {
     if (segGroup[i] != 0) break;
     }
     const MdcSeg *last = segGroup[i];
     MdcSegInfoSterO* lastInfo = (MdcSegInfoSterO *) last->info();
 
  // Test that slope from last segment to new one is roughly = slope from
  //   first seg to last
  double arcLast = lastInfo->arc();
  double arcFirst = firstInfo->arc();
  double arcNew = newInfo->arc();
 
  double delZold = lastInfo->z0() - firstInfo->z0() +
  lastInfo->ct() * arcLast - firstInfo->ct() * arcFirst;
  double delArcold = arcLast - arcFirst;
  double delZnew = newInfo->z0() - lastInfo->z0() +
  newInfo->ct() * arcNew - lastInfo->ct() * arcLast;
  double delArcnew = arcNew - arcLast;
  double p1 = delZold * delArcnew;
  double p2 = delZnew * delArcold;
 
  cout << " test: " << p1 << "  " << p2 << endl;
 
  if (fabs(p2 - p1) > 0.02) return -1;
 
*/
  return 0;
}
 
//**************************************************************************
void MdcSegGrouperSt::resetComb( const MdcSeg* seed ) {
  //**************************************************************************
 
  // Delete existing list of valid/invalid segs
  if ( isValid != 0 )
  {
    int i;
    for ( i = 0; i < nDeep; i++ )
    {
      delete[] isValid[i];
      isValid[i] = 0;
    }
  }
 
  // Grab the seglist for each slayer
  int islay  = 0;
  int iply   = 0;
  nPlyFilled = 0;
  nNull      = 0;
 
  // Set up all sorts of stuff for fast grouping of segs in nextGroup()
  for ( const MdcSuperLayer* thisSlay = _gm->firstSlay(); thisSlay != 0;
        thisSlay                      = thisSlay->next() )
  {
    // bool noGoodYet = true;
    islay++;
    if ( thisSlay->whichView() == 0 ) continue;
    firstGood[iply] = 0;
    firstBad[iply]  = segList[islay - 1].length();
    if ( firstGood[iply] > firstBad[iply] ) firstGood[iply] = firstBad[iply];
    if ( firstGood[iply] == firstBad[iply] )
    {
      // If there are no valid segs for this ply, skip it
      continue;
    }
    // Associate correct seglist with this ply
    combList[iply] = &segList[islay - 1];
    leaveGap[iply] = false;
    iply++;
  }
 
  nPlyFilled = iply;
  resetSegCounters();
  maxNull = nPlyFilled - 2;
}
 
//---------------------------------------------------------------------
MdcTrack* MdcSegGrouperSt::storePar( MdcTrack* trk, double parms[2], double chisq, TrkContext&,
                                     double ) {
  //---------------------------------------------------------------------
  assert( trk != 0 );
  TrkHelixMaker maker;
  if ( 3 == _debug )
    std::cout << "-----storePar z0 " << parms[0] << " ct " << parms[1] << std::endl;
  maker.addZValues( trk->track(), parms[0], parms[1], chisq );
  return trk;
}