TMDCTsf.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TMDCTsf.cxx,v 1.23 2012/05/28 05:16:29 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TMDCTsf.h
// Section  : Tracking MDC
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to represent a Track Finder Segment(TSF).
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
// #ifdef TRKRECO_DEBUG_DETAIL
// #define TRKRECO_DEBUG
// #endif
#include <iostream>
// #define HEP_SHORT_NAMES
#include "CLHEP/Vector/ThreeVector.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCLayer.h"
#include "TrkReco/TMDCTsf.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
 
#include "CLHEP/Matrix/Vector.h"
#include <math.h>
 
using CLHEP::HepVector;
 
extern TMDC* GMDC;
 
TMDC* TMDCTsf::_cdc = 0;
/*
TMDCTsf::TMDCTsf(const TMDCWire * const w) {
    //...Notice...
    // This function creates TSF of 1,2,3,4(4 wire layers) or 1,2,3(3 wire
    // layers) configuration of wires in each layer. The argument w must be
    // in the first wire layer in super layer, or must be in the second.
 
    AList<int> offset;
    offset.append(new int(0));
    offset.append(new int(-1));
    offset.append(new int(-1));
    offset.append(new int(-2));
 
    AList<int> width;
    width.append(new int(1));
    width.append(new int(2));
    width.append(new int(3));
    width.append(new int(4));
 
    create(offset, width, w);
}
 
TMDCTsf::TMDCTsf(const AList<int> & offset,
                 const AList<int> & width,
                 const TMDCWire * const w) {
    create(offset, width, w);
 
}*/
 
TMDCTsf::TMDCTsf( unsigned sl ) : _sl( sl ), _angle( 0.8 ) {}
 
TMDCTsf::~TMDCTsf() {}
 
void TMDCTsf::dump( const std::string& msg, const std::string& pre ) const {
  std::cout << pre;
  std::cout << "tsf : " << _wires.length() << " : ";
  for ( unsigned i = 0; i < _wires.length(); i++ )
  {
    TMDCWire* w = _wires[i];
    std::cout << w->layerId() << "-" << w->localId() << ", ";
  }
  std::cout << std::endl;
}
 
bool TMDCTsf::create( const AList<int>& offset, const AList<int>& width,
                      const TMDCWire* const w ) {
 
  //...Check wire position...
  unsigned localId        = w->localLayerId();
  unsigned superId        = w->superLayerId();
  unsigned availableLayer = GMDC->nLocalLayer( superId ) - localId;
  if ( availableLayer < TMDCTSF_MIN_LAYERS )
  {
    std::cout << "TMDCTsf::create !!! can not create TSF for wire ";
    std::cout << superId << "-" << localId << " because of ";
    std::cout << "lack of layer(s)." << std::endl;
    return false;
  }
  availableLayer = ( availableLayer > offset.length() ) ? offset.length() : availableLayer;
 
  //...Check offset...
  bool offsetLayer = false;
  if ( w->layer()->offset() != 0. ) offsetLayer = true;
 
  //...Loop over available layers...
  for ( unsigned l = 0; l < availableLayer; l++ )
  {
    unsigned absLayerId = w->layerId() + l;
    int      wireId[2];
    wireId[0] = w->localId() + ( *offset[l] );
    wireId[1] = w->localId() + ( *offset[l] );
    if ( offsetLayer ) ++wireId[1];
 
    //...Loop over wires...
    for ( unsigned j = 0; j < *( width[l] ); j++ )
    { _wires.append( (TMDCWire*)GMDC->wire( absLayerId, wireId[l % 2] + j ) ); }
  }
 
  return true;
}
 
void TMDCTsf::solveLeftRight( void ) {
  //...For debug...
  // this->dump();
 
  //     TMDCWire * innerWire = _wires[0];
  //     TMDCWireHit * innerHit = innerWire->hit();
 
  //     unsigned outerMostLayer = _wires.last()->layerId();
  //     unsigned i = _wires.length();
  //     while (_wires[--i]->layerId() == outerMostLayer) {
  //    TMDCWireHit * outerHit = _wires[i]->hit();
  //    if (! outerHit) continue;
 
  //    //...Find four lines by two hit points...
  //    //  AList<TMDCLine> lines = innerHit->lines(outerHit);
  //    // AList<TMDCLine> lines; // = innerHit->lines(outerHit);
 
  //    //...Line loop...
  //    for (unsigned ii = 0; ii < 4; ii++) {
  //        TMDCLine * l = lines[ii];
 
  //        //...For debug...
  //        // l->dump();
 
  //        //...Find hits on this line...
  //        unsigned j = i;
  //        while (TMDCWire * w = _wires[--j]) {
 
  //        //...Check condition...
  //        if (w == innerWire) break;
  //        if (w->layerId() == outerMostLayer) continue;
  //        TMDCWireHit * hit = w->hit();
  //        if (! hit) continue;
 
  //        //...Append this hit...
  //        l->appendHit(hit);
 
  //        //...Obtain distance of closest approach...
  //        double distance = l->distance(* l->closest().last());
 
  //        //...Cal. tolerable distance...
  //        unsigned leftRight = * (l->leftRight().last());
  //        double tolerable = 5. * hit->dDistance(leftRight);
 
  //        //...Is this OK?...
  //        if (distance > tolerable) {
  //            l->removeLastHit();
  //            continue;
  //        }
 
  //        //...For debug...
  //        // hit->dump("", "    ");
  //        // std::cout << "distance = " << distance << std::endl;
  //        }
 
  //        //...Check # of hits on this line
  //        if (l->hits().length() > 2) {
  //        unsigned i = 0;
  //        while (TMDCWireHit * h = (l->hits())[i]) {
  //            unsigned leftRight = * (l->leftRight())[i];
  //            unsigned newState = h->state() | (1 << leftRight);
  //            h->state(newState);
  //            ++i;
 
  //            //...For debug...
  //            // h->dump("", "    ");
  //        }
  //        }
 
  //        //...Delete line here...
  //        delete l;
  //    }
  //     }
}
 
//----------------the following is added by liuqg--------------------//
AList<TSegment> TMDCTsf::segments( const AList<TMLink>& links ) {
  _all.removeAll();
  _all.append( links );
 
  AList<TMLink> seeds1;
  AList<TMLink> seeds2;
  for ( int i = 0; i < _all.length(); ++i )
  {
    if ( _all[i]->wire()->localLayerId() == 0 ) seeds1.append( _all[i] );
    else if ( _all[i]->wire()->localLayerId() == 1 ) seeds2.append( _all[i] );
    else continue;
  }
  AList<TSegment> list = createTsf( 0 ); // get 1234 tsf.
  unsigned        n    = list.length();
  AList<TSegment> splitted;
  AList<TMLink>   seedNeighbors;
  for ( unsigned i = 0; i < n; ++i )
  {
    seedNeighbors.removeAll();
    TSegment* c = list[i];
 
    // append seed's neighbors
    for ( int j = 0; j < c->links().length(); ++j )
    {
      if ( c->links()[j]->wire()->localLayerId() > 0 ) continue;
      unsigned seedId = c->links()[j]->wire()->localId();
      for ( int k = 0; k < seeds1.length(); ++k )
      {
        if ( seeds1[k]->wire()->localIdForPlus() + 1 == seedId ||
             seeds1[k]->wire()->localIdForMinus() - 1 == seedId )
          seedNeighbors.append( seeds1[k] );
      }
      break;
    }
 
    AList<TSegment> newClusters = c->splitTsf( seedNeighbors );
    //  cout<<"newClusters.length = "<<newClusters.length()<<endl;
    splitted.append( c ); // remove the segment created by TMDCTsf here!!!
    if ( newClusters.length() == 0 ) continue;
    list.append( newClusters );
  }
  list.remove( splitted );
  HepAListDeleteAll( splitted );
  // cout<<"sl = "<<_sl<<"    list1234 = "<<list.length()<<endl;
 
  AList<TSegment> list2 = createTsf( 1 ); // get 234 tsf.
  n                     = list2.length();
  AList<TSegment> splitted2;
 
  for ( unsigned i = 0; i < n; ++i )
  {
    seedNeighbors.removeAll();
    TSegment* c = list2[i];
 
    // append seed's neighbors
    for ( int j = 0; j < c->links().length(); ++j )
    {
      if ( c->links()[j]->wire()->localLayerId() != 1 ) continue;
      unsigned seedId = c->links()[j]->wire()->localId();
      for ( int k = 0; k < seeds2.length(); ++k )
      {
        if ( seeds2[k]->wire()->localIdForPlus() + 1 == seedId ||
             seeds2[k]->wire()->localIdForMinus() - 1 == seedId )
          seedNeighbors.append( seeds2[k] );
      }
      break;
    }
 
    AList<TSegment> newClusters2 = c->splitTsf( seedNeighbors );
    //  cout<<"newClusters2.length = "<<newClusters2.length()<<endl;
    splitted2.append( c );
    if ( newClusters2.length() == 0 ) continue;
    list2.append( newClusters2 );
  }
  list2.remove( splitted2 );
  HepAListDeleteAll( splitted2 );
  // cout<<"sl = "<<_sl<<"    list234 = "<<list2.length()<<endl;
 
  list.append( list2 );
 
  // check segment, remove the multi segs.
  AList<TSegment> badList;
  for ( int i = 0; i < list.length(); ++i )
  {
    for ( int j = i + 1; j < list.length(); ++j )
    {
      AList<TMLink> links = list[j]->links();
      AList<TMLink> multiLinks;
      for ( int k = 0; k < links.length(); ++k )
      {
        TMLink& l = *links[k];
        if ( list[i]->links().hasMember( l ) ) multiLinks.append( l );
      }
      if ( multiLinks.length() < 2 ) continue;
      unsigned minLength = links.length();
      if ( links.length() > list[i]->links().length() ) minLength = list[i]->links().length();
      if ( minLength - multiLinks.length() > 1 ) continue;
      int nHits[4]    = { 0 }; // in each layer.
      int multiLayers = 0;
      for ( int k = 0; k < multiLinks.length(); ++k )
        ++nHits[multiLinks[k]->hit()->wire()->localLayerId()];
      for ( int k = 0; k < 4; ++k )
        if ( nHits[k] > 0 ) ++multiLayers;
      if ( multiLayers >= 2 && ( links.length() > list[i]->links().length() ) )
        badList.append( list[i] );
      else if ( multiLayers >= 2 ) badList.append( list[j] );
      else continue;
    }
  }
  list.remove( badList );
 
  // reset links in badList.
  for ( int i = 0; i < badList.length(); ++i )
  {
    AList<TMLink> bads = badList[i]->links();
    for ( int j = 0; j < bads.length(); ++j )
    {
      unsigned n = bads[j]->tsfTag();
      if ( n == 0 ) continue;
      else
      {
        --n;
        bads[j]->tsfTag( n );
      }
    }
  }
  HepAListDeleteAll( badList );
 
  return list;
}
 
AList<TSegment> TMDCTsf::createTsf( unsigned seedlayer ) const {
  if ( _cdc == 0 ) _cdc = TMDC::getTMDC();
 
  AList<TMLink>   seeds;
  AList<TMLink>   exhits;
  AList<TSegment> list;
 
  unsigned n = _all.length();
  if ( n < 3 ) return list;
  if ( seedlayer == 1 && _sl == 10 ) return list;
 
  for ( unsigned i = 0; i < n; ++i )
  {
    TMLink* ll = _all[i];
    //  if (ll->wire()->localLayerId() == seedlayer && ll->tsfTag() == 0) seeds.append(ll);
    if ( ll->wire()->localLayerId() == seedlayer ) seeds.append( ll );
    if ( ll->wire()->localLayerId() > seedlayer ) exhits.append( ll );
  } // get seeds and exhits.
 
  for ( unsigned i = 0; i < seeds.length(); ++i )
  {
    AList<TMLink> forSegment;
    TMLink*       seed  = seeds[i];
    unsigned      lId   = seed->wire()->layerId();
    int           local = (int)seed->wire()->localId();
    // cout<<"lId: "<<lId<<" loId: "<<local<<endl;
    float phi0 = _cdc->layer( lId )->offset();
    int   nWir = (int)_cdc->layer( lId )->nWires();
    float phi  = phi0 + local * 2 * M_PI / nWir; // phi angle of seed.
    //  double angle = alpha(sl, seed);
    forSegment.append( seed );
 
    // for check
    int hitsOnLayer[3] = { 0, 0, 0 };
 
    for ( unsigned m = 0; m < exhits.length(); ++m )
    {
      TMLink*  l        = exhits[m];
      unsigned lLayer   = l->wire()->localLayerId();
      unsigned tmpId    = l->wire()->layerId();
      int      tmpLocal = l->wire()->localId();
      float    phi0tmp  = _cdc->layer( tmpId )->offset();
      int      nWirtmp  = (int)_cdc->layer( tmpId )->nWires();
      int      localOffset =
          (int)( ( phi - phi0tmp ) / ( 2 * M_PI / nWirtmp ) ); // get the corresponding localId
                                                               // of exhits.
 
      int diff = abs( tmpLocal - localOffset );
 
      if ( diff > nWirtmp / 2 ) diff = nWirtmp - diff;
      // cout<<"lId: "<<tmpId<<" loId: "<<tmpLocal<<" diff wire: "<<diff
      //<<" cal value = "<<fabs(acos((l->position() -
      // seed->position()).unit().dot(-seed->position().unit())))<<endl;
      if ( diff > lLayer + 2 ) continue; // check in natural phase..
 
      //        double DRIFT = l->cDrift();
      //        const HepPoint3D dirZ(0.0, 0.0, 1.0);
      //        HepPoint3D pos1 = l->positionD() - DRIFT *
      // dirZ.cross(l->positionD().unit());         HepPoint3D pos2 = l->positionD() + DRIFT *
      // dirZ.cross(l->positionD().unit());         if (fabs(acos((pos1 -
      // seed->position()).unit().dot(-seed->position().unit()))) < _angle
      //        || fabs(acos((pos2 -
      // seed->position()).unit().dot(-seed->position().unit()))) < _angle) {
      if ( fabs( acos( ( l->position() - seed->position() )
                           .unit()
                           .dot( -seed->position().unit() ) ) ) < _angle )
      {
        forSegment.append( l );
        ++hitsOnLayer[lLayer - 1 - seedlayer];
      }
    }
 
    // check
    int fireLayers = 0;
    for ( unsigned j = 0; j < 3; ++j )
    {
      if ( hitsOnLayer[j] > 0 ) ++fireLayers;
    }
    if ( fireLayers < 2 ) { continue; }
    list.append( new TSegment( forSegment ) );
  }
  return list;
}
 
double TMDCTsf::alpha( unsigned sl, TMLink* seed ) const {
  double     radius[11] = { 11.525, 16.160, 24.555, 31.060, 37.455, 44.760,
                            51.370, 57.855, 64.165, 71.570, 76.345 };
  HepPoint3D pa, pb, pc;
  double     rcell;
  if ( sl < 2 ) rcell = 0.6;
  else rcell = 0.8;
  pa.set( 0, 2 / radius[sl], 0 );
  pb.set( -2 * rcell / ( rcell * rcell + ( radius[sl] + rcell ) * ( radius[sl] + rcell ) ),
          2 * ( radius[sl] + rcell ) /
              ( rcell * rcell + ( radius[sl] + rcell ) * ( radius[sl] + rcell ) ),
          0 );
  pc.set( 2 * rcell / ( rcell * rcell + ( radius[sl] + rcell ) * ( radius[sl] + rcell ) ),
          2 * ( radius[sl] + rcell ) /
              ( rcell * rcell + ( radius[sl] + rcell ) * ( radius[sl] + rcell ) ),
          0 );
  cout << "sl = " << sl
       << "alpha2 = " << acos( ( pb - pa ).unit().dot( ( pc - pa ).unit() ) ) / 2 << endl;
  HepPoint3D pa2, pb2, pc2;
  pa2.set( 0, 2 / radius[sl], 0 );
  pb2.set( -6 * rcell /
               ( 9 * rcell * rcell + ( radius[sl] + rcell ) * ( radius[sl] + rcell ) ),
           2 * ( radius[sl] + rcell ) /
               ( 9 * rcell * rcell + ( radius[sl] + rcell ) * ( radius[sl] + rcell ) ),
           0 );
  pc2.set( 6 * rcell / ( 9 * rcell * rcell + ( radius[sl] + rcell ) * ( radius[sl] + rcell ) ),
           2 * ( radius[sl] + rcell ) /
               ( 9 * rcell * rcell + ( radius[sl] + rcell ) * ( radius[sl] + rcell ) ),
           0 );
  cout << "alpha22 = " << acos( ( pb2 - pa2 ).unit().dot( ( pc2 - pa2 ).unit() ) ) / 2 << endl;
  return 1.0;
}