MdcSeg.cxx

Go to the documentation of this file.

// MdcSeg.cxx
 
#include "MdcTrkRecon/MdcSeg.h"
#include "MdcData/MdcHit.h"
#include "MdcData/MdcHitMap.h"
#include "MdcData/MdcHitUse.h"
#include "MdcGeom/BesAngle.h"
#include "MdcGeom/MdcLayer.h"
#include "MdcGeom/MdcSuperLayer.h"
#include "MdcTrkRecon/MdcLine.h"
#include "MdcTrkRecon/MdcMap.h"
#include "MdcTrkRecon/MdcSegInfoSterO.h"
#include "MdcTrkRecon/MdcSegParams.h"
#include "MdcTrkRecon/MdcSegUsage.h"
#include "MdcTrkRecon/mdcWrapAng.h"
#include "MdcTrkRecon/mdcWrapWire.h"
#include <math.h>
#include <stdlib.h>
 
// yzhang hist cut
#include "AIDA/IHistogram1D.h"
extern AIDA::IHistogram1D* g_nSigAdd;
// zhangy hist cut
 
extern int    haveDigiTk[43][288];
extern double haveDigiPt[43][288];
extern double haveDigiTheta[43][288];
extern double haveDigiPhi[43][288];
extern int    haveDigiAmbig[43][288];
// Initialize static pointer
MdcSegParams* MdcSeg::segParam = 0;
const double  twoPi            = Constants::twoPi;
//------------------------------------------------------------------------
MdcSeg::MdcSeg( double bt ) {
  //------------------------------------------------------------------------
  _info      = 0;
  _bunchTime = bt;
}
 
//------------------------------------------------------------------------
MdcSeg::~MdcSeg() {
  //------------------------------------------------------------------------
  if ( _info != 0 ) delete _info;
  reset(); // delete Hots
}
 
//------------------------------------------------------------------------
MdcSeg::MdcSeg( const MdcSeg& other )
    : GmsListLink()
    , _slayer( other._slayer )
    , _phi( other._phi )
    , _slope( other._slope )
    , _chisq( other._chisq )
    , _qual( other._qual )
    , _pattern( other._pattern )
    , _info( other._info )
    , _bunchTime( other._bunchTime )
//------------------------------------------------------------------------
{
  HepAListDeleteAll( _theList );
  for ( int i = 0; i < other.nHit(); i++ ) { _theList.append( other.hit( i ) ); }
  for ( int j = 0; j < 3; j++ )
  {
    _errmat[0] = other._errmat[0];
    _errmat[1] = other._errmat[1];
    _errmat[2] = other._errmat[2];
  }
  segParam = other.segParam;
}
 
//------------------------------------------------------------------------
MdcSeg& MdcSeg::operator=( const MdcSeg& other ) {
  //------------------------------------------------------------------------
  if ( &other != this )
  {
 
    HepAListDeleteAll( _theList );
    for ( int i = 0; i < other.nHit(); i++ ) { _theList.append( other.hit( i ) ); }
    _slayer    = other._slayer;
    _phi       = other._phi;
    _slope     = other._slope;
    _errmat[0] = other._errmat[0];
    _errmat[1] = other._errmat[1];
    _errmat[2] = other._errmat[2];
    _chisq     = other._chisq;
    _qual      = other._qual;
    _pattern   = other._pattern;
    _info      = other._info;
    _bunchTime = other._bunchTime;
    segParam   = other.segParam;
  }
 
  return *this;
}
 
//------------------------------------------------------------------------
void MdcSeg::setInfo( MdcSegInfo* newInfo ) {
  //------------------------------------------------------------------------
  delete _info; // if any
  _info = newInfo;
}
 
//------------------------------------------------------------------------
void MdcSeg::setValues( int nInPatt, int nhit, MdcHit* hits[], MdcLine* span,
                        const MdcSuperLayer* slay, int ambig[] ) {
  //------------------------------------------------------------------------
  _qual = 0;
  if ( nInPatt == 4 ) _qual |= segFullFlag;
  _phi       = BesAngle( span->intercept );
  _slope     = span->slope;
  _chisq     = span->chisq;
  _errmat[0] = span->errmat[0];
  _errmat[1] = span->errmat[1];
  _errmat[2] = span->errmat[2];
  reset();
  _slayer = slay;
  for ( int i = 0; i < nhit; i++ )
  {
    MdcHitUse* alink = new MdcHitUse( *( hits[i] ), superlayer()->rad0(), ambig[i] );
    append( alink );
  }
 
  return;
}
 
//------------------------------------------------------------------------
void MdcSeg::setValues( int nInPatt, double inPhi, double inSlope, double chi2,
                        double inError[3], const MdcSuperLayer* slay ) {
  //------------------------------------------------------------------------
  // Sets segment values with no associated hits
  _qual = 0;
  if ( nInPatt == 4 ) _qual |= segFullFlag;
  _phi       = inPhi;
  _slope     = inSlope;
  _chisq     = chi2;
  _errmat[0] = inError[0];
  _errmat[1] = inError[1];
  _errmat[2] = inError[2];
  _slayer    = slay;
  reset(); // clears hit list
 
  return;
}
 
//------------------------------------------------------------------------
void MdcSeg::markHits( const MdcMap<const MdcHit*, MdcSegUsage*>& usedHits ) const {
  //------------------------------------------------------------------------
  for ( int i = 0; i < nHit(); i++ )
  {
    MdcHitUse*   alink = _theList[i];
    MdcSegUsage* x;
    if ( usedHits.get( alink->mdcHit(), x ) ) x->setUsedAmbig( alink->ambig() );
  }
}
 
//------------------------------------------------------------------------
void MdcSeg::plotSegAll() const {
  //------------------------------------------------------------------------
  // print hit
  // if(superlayer()!=NULL) std::cout<<"sl"<<superlayer()->slayNum()<<" p"<<segPattern()<<"
  // st"<<quality();
  for ( int ihit = 0; ihit < nHit(); ihit++ )
  {
    hit( ihit )->mdcHit()->print( std::cout );
    std::cout << setw( 2 ) << hit( ihit )->ambig() << " ";
  }
 
  // cout << setiosflags( ios::fixed );
  if ( info() != NULL )
  {
    std::cout << " phi " << setprecision( 2 ) << phi() << " slope " << std::setw( 2 )
              << setprecision( 2 ) << slope() << " ";
    if ( superlayer()->whichView() == 0 )
    {
      std::cout << setprecision( 2 ) << "phi0=" << info()->par( 0 );
      cout << setprecision( 5 ) << " cpa=" << info()->par( 1 );
    }
    else
    {
      std::cout << setprecision( 2 ) << "z0=" << info()->par( 0 ) << setprecision( 2 )
                << " ct=" << info()->par( 1 );
    }
    if ( fabs( info()->arc() ) > 0.0001 )
    { std::cout << setprecision( 2 ) << " arc=" << info()->arc(); }
    std::cout << setprecision( 3 ) << " chi2=" << _chisq;
  }
  else
  {
    std::cout << " phi " << setprecision( 2 ) << phi() << " slope " << std::setw( 2 )
              << setprecision( 2 ) << slope() << " chi2 " << setprecision( 3 ) << chisq();
  }
 
  std::cout << std::defaultfloat;
  // cout << setprecision( 6 );
  // cout << setiosflags( ios::scientific );
}
//------------------------------------------------------------------------
void MdcSeg::plotSeg() const {
  //------------------------------------------------------------------------
  std::cout << superlayer()->slayNum() << " pat" << segPattern() << " ";
  for ( int ihit = 0; ihit < nHit(); ihit++ )
  {
    hit( ihit )->mdcHit()->print( std::cout ); // print hit
    std::cout << hit( ihit )->ambig() << " ";
  }
  if ( info() != NULL ) { cout << " . "; }
  // std::cout << std::endl;
}
 
//------------------------------------------------------------------------
// void
// MdcSeg::plotZ(int openIt, int color) const {
//------------------------------------------------------------------------
/*
// Plot line to indicate segment
double start[2], finish[2];
 
const double length = 0.018 * display->width(windowZ);
double ct = ((MdcSegInfoSterO *) info())->ct();   // !!! this cast stinks
double z0 = ((MdcSegInfoSterO *) info())->z0();
double arc = ((MdcSegInfoSterO *) info())->arc();
double dels = length / sqrt(1. + ct*ct);
double z = z0 + ((MdcSegInfoSterO *) info())->arc() * ct;
start[0] = arc - dels;
start[1] = z - ct * dels;
finish[0] = arc + dels;
finish[1] = z + ct * dels;
*/
//}
 
//------------------------------------------------------------------------
int MdcSeg::addHits( MdcLine* span, MdcHit** /* hits[]*/, const MdcHitMap* map, double corr ) {
  //------------------------------------------------------------------------
 
  /* Pick up hits adjacent to an existing segment.  Output final number of
     hits.  Note: input hits are not refound.  It picks up hits
     in cells that the segment should pass through, checking that they
     have a plausible drift distance.  In event of a wraparound (i.e. track that
     passes through 2pi), all angles are relative to phiseg (i.e. if phiseg is
     just above 0, some phi's may be negative).  */
  //****
 
  int    cell[2], ambig[2]; // wire # and ambig for cells in current layer
  double phiwire[2];
  int    cellused[4] = { 0 }; // list of wire #'s of existing hits in span
  int    lAdded      = 0;
  int    nhits       = nHit();
  int    m_debug     = 0;
 
  // Note the hits already in the fit, so we don't pick them up again.
  int firstnum = superlayer()->firstLayer()->layNum();
  for ( int i = 0; i < nHit(); i++ )
  {
    const MdcHit* dHit          = hit( i )->mdcHit();
    int           laynum        = dHit->layernumber();
    cellused[laynum - firstnum] = dHit->wirenumber();
  }
 
  // Loop through the layers, predicting cells that could be hit.
  //  for (int layIndex = 0; layIndex < 4; layIndex++) {
  for ( int layIndex = 0; layIndex < superlayer()->nLayers(); layIndex++ )
  {
    const MdcLayer* layer  = superlayer()->layer( layIndex );
    int             laynum = layer->layNum();
 
    // Calc. projected phi of hit
    double rinv     = 1. / layer->rMid();
    double ncellinv = 1. / (double)layer->nWires();
    double phiproj  = phi() + ( layer->rMid() - superlayer()->rad0() ) * slope();
    // Calc. nearest wire.
    BesAngle tmp( phiproj - layer->phiOffset() );
    cell[0]    = (int)floor( layer->nWires() * tmp.rad() / twoPi + 0.5 );
    cell[0]    = mdcWrapWire( cell[0], layer->nWires() );
    phiwire[0] = mdcWrapAng( phi(), cell[0] * twoPi * ncellinv + layer->phiOffset() );
    // Which ambiguity? +1 if phi(wire) < phi(projected).
    ambig[0] = ( phiwire[0] < phiproj ) ? 1 : -1;
 
    // Now next nearest wire.
    ambig[1]   = -ambig[0];
    cell[1]    = mdcWrapWire( cell[0] + ambig[0], layer->nWires() );
    phiwire[1] = mdcWrapAng( phi(), cell[1] * twoPi * ncellinv + layer->phiOffset() );
 
    if ( m_debug ) std::cout << " loop over the two possible wires  " << std::endl;
    //*** Loop over the two possible wires.
    for ( int iroad = 0; iroad < 2; iroad++ )
    {
      if ( cellused[laynum - firstnum] == cell[iroad] ) continue;
      if ( m_debug )
        std::cout << "possible wires  " << laynum << " " << cell[iroad] << std::endl;
      if ( map->hitWire( laynum, cell[iroad] ) != 0 )
      {
        MdcHit* ahit = map->hitWire( laynum, cell[iroad] );
        // if hit is already used, skip it!
        if ( ahit->usedHit() )
        {
          if ( m_debug ) std::cout << "hit used continue   " << std::endl;
          continue;
        }
        // drift as delphi
        // FIXME: flip ambiguity for incoming tracks
        double phidrift = ahit->driftDist( bunchTime(), ambig[iroad] ) * corr * rinv;
        double phihit   = mdcWrapAng( phi(), phidrift * ambig[iroad] + ahit->phi() );
 
        // Check the drift distance.
        double sigphi = corr * ahit->sigma( bunchTime(), 0 ) * rinv;
        // Skip hit if more than n sigma away from track.
        if ( g_nSigAdd && fabs( sigphi ) > 0.0001 )
        { g_nSigAdd->fill( fabs( phihit - phiproj ) / sigphi ); } // yzhang hist cut
        if ( fabs( phihit - phiproj ) > sigphi * segParam->nsigAddHit )
        {
          if ( m_debug )
            std::cout << fabs( phihit - phiproj ) << "> add hit sigma   " << sigphi << "*"
                      << segParam->nsigAddHit << "=" << sigphi * segParam->nsigAddHit
                      << std::endl;
          continue;
        }
 
        // Load hit for refitting
        lAdded             = 1;
        span->sigma[nhits] = sigphi;
        span->x[nhits]     = layer->rMid() - superlayer()->rad0();
        span->y[nhits]     = mdcWrapAng( span->y[0], phihit );
 
        // Add this hit to segment
        MdcHitUse* alink = new MdcHitUse( *ahit, superlayer()->rad0(), ambig[iroad] );
        append( alink );
        // std::cout<< __FILE__ << "   " << __LINE__ << " addhit  "<<std::endl;
 
        nhits++;
      } // end if hit wire
 
    } // end loop over 2 cells
 
  } // end layer loop
 
  if ( lAdded )
  {
    span->fit( nhits );
    BesAngle tmpAngle( span->intercept );
    span->intercept = tmpAngle;
    _phi            = span->intercept;
    _slope          = span->slope;
    _chisq          = span->chisq;
    _errmat[0]      = span->errmat[0];
    _errmat[1]      = span->errmat[1];
    _errmat[2]      = span->errmat[2];
  }
 
  return nhits;
}
 
//------------------------------------------------------------------------
void MdcSeg::reset() {
  //------------------------------------------------------------------------
  HepAListDeleteAll( _theList );
}
 
//------------------------------------------------------------------------
void MdcSeg::append( MdcHitUse* theHitUse ) {
  //------------------------------------------------------------------------
  _theList.append( theHitUse );
}
 
//------------------------------------------------------------------------
void MdcSeg::remove( MdcHitUse* theHitUse ) {
  //------------------------------------------------------------------------
  _theList.remove( theHitUse );
  delete theHitUse;
}
 
//------------------------------------------------------------------------
int MdcSeg::nHit() const {
  //------------------------------------------------------------------------
  return _theList.length();
}
 
//------------------------------------------------------------------------
double MdcSeg::testCombSeg( const MdcSeg* testSeg ) const {
  //------------------------------------------------------------------------
  int tkId = -1;
  for ( int i = 0; i < nHit(); i++ )
  {
    const MdcHit* h = hit( i )->mdcHit();
    unsigned int  l = h->layernumber();
    unsigned int  w = h->wirenumber();
    // std::cout<< __FILE__ << " ref  " << i << "
    // haveDigiTk("<<l<<","<<w<<")"<<haveDigiTk[l][w]<<std::endl;
    if ( haveDigiTk[l][w] < 0 || haveDigiTk[l][w] > 100 ) continue;
    if ( tkId < 0 ) { tkId = haveDigiTk[l][w]; }
    else if ( haveDigiTk[l][w] != tkId )
    {
      return -1; // hits in this seg not in same mc track
    }
  } // end for
 
  double nSame = 0.;
  for ( int i = 0; i < testSeg->nHit(); i++ )
  {
    const MdcHit* h = testSeg->hit( i )->mdcHit();
    unsigned int  l = h->layernumber();
    unsigned int  w = h->wirenumber();
    if ( haveDigiTk[l][w] == tkId ) { ++nSame; }
  }
 
  return nSame / testSeg->nHit();
}
 
//------------------------------------------------------------------------
double MdcSeg::testCombSegPt() const {
  //------------------------------------------------------------------------
  double truthPt = -1.;
  for ( int i = 0; i < nHit(); i++ )
  {
    const MdcHit* h = hit( i )->mdcHit();
    unsigned int  l = h->layernumber();
    unsigned int  w = h->wirenumber();
    if ( haveDigiPt[l][w] < 0 ) continue;
    // std::cout<< __FILE__ << "   " << __LINE__ << "
    // haveDigiPt("<<l<<","<<w<<")"<<haveDigiPt[l][w]<<std::endl;
    if ( truthPt < 0 ) { truthPt = haveDigiPt[l][w]; }
  } // end for
 
  return truthPt;
}
 
//------------------------------------------------------------------------
double MdcSeg::testCombSegTheta() const {
  //------------------------------------------------------------------------
  double truthTheta = -999.;
  for ( int i = 0; i < nHit(); i++ )
  {
    const MdcHit* h = hit( i )->mdcHit();
    unsigned int  l = h->layernumber();
    unsigned int  w = h->wirenumber();
    if ( haveDigiTheta[l][w] < -998. ) continue;
    // std::cout<< __FILE__ << "   " << __LINE__ << "
    // haveDigiTheta("<<l<<","<<w<<")"<<haveDigiTheta[l][w]<<std::endl;
    if ( truthTheta < -998. ) { truthTheta = haveDigiTheta[l][w]; }
  } // end for
 
  return truthTheta;
}
 
//------------------------------------------------------------------------
double MdcSeg::testCombSegPhi() const {
  //------------------------------------------------------------------------
  double truthPhi = -999.;
  for ( int i = 0; i < nHit(); i++ )
  {
    const MdcHit* h = hit( i )->mdcHit();
    unsigned int  l = h->layernumber();
    unsigned int  w = h->wirenumber();
    if ( haveDigiPhi[l][w] < -998. ) continue;
    // std::cout<< __FILE__ << "   " << __LINE__ << "
    // haveDigiPhi("<<l<<","<<w<<")"<<haveDigiPhi[l][w]<<std::endl;
    if ( truthPhi < -998. ) { truthPhi = haveDigiPhi[l][w]; }
  } // end for
 
  return truthPhi;
}
 
//------------------------------------------------------------------------
double MdcSeg::testCombSegAmbig() const {
  //------------------------------------------------------------------------
  double ambigOk = 0.;
  for ( int i = 0; i < nHit(); i++ )
  {
    const MdcHit* h = hit( i )->mdcHit();
    unsigned int  l = h->layernumber();
    unsigned int  w = h->wirenumber();
    if ( hit( i )->ambig() == haveDigiAmbig[l][w] ) ambigOk++;
  } // end for
 
  return ambigOk / nHit();
}