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    : [email protected]
//-----------------------------------------------------------------------------
// 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/TMDCWire.h"
#include "TrkReco/TMDCLayer.h"
#include "TrkReco/TMDCTsf.h"
#include "TrkReco/TMDCWireHit.h"
 
#include <math.h>
#include "CLHEP/Matrix/Vector.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*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*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;
}