THistogram.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: THistogram.cxx,v 1.8 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : THistogram.h
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : [email protected]
//-----------------------------------------------------------------------------
// Description : A class for a histogram used in tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include <stdio.h>
#include <stdlib.h>
 
#include "TrkReco/THistogram.h"
#include "TrkReco/TCircle.h"
 
THistogram::THistogram(unsigned nBins) : _nBins(nBins) {
    _binSize = 2. * M_PI / (float) _nBins;
    _bins = (unsigned *) malloc(_nBins * sizeof(unsigned));
    _masks = (bool *) malloc(_nBins * sizeof(bool));
    _links = (AList<TMLink> **) malloc(_nBins * sizeof(AList<TMLink> *));
    for (unsigned i = 0; i < _nBins; i++) {
    _bins[i] = 0;
    _masks[i] = false;
    _links[i] = new AList<TMLink>;
    }
}
 
THistogram::~THistogram() {
    free(_bins);
    free(_masks);
    for (unsigned i = 0; i < _nBins; i++)
    delete _links[i];
    free(_links);
}
 
void
THistogram::dump(const std::string & msg, const std::string & pre) const {
    std::cout << pre;
    std::cout << "THistogram dump:#bins=" << _nBins << std::endl;
    unsigned nLoops = _nBins / 15 + 1;
    unsigned n0 = 0;
    unsigned n1 = 0;
    for (unsigned i = 0; i < nLoops; i++) {
    for (unsigned j = 0; j < 15; j++) {
        if (n0 == _nBins) break;
        printf("%4d", n0);
        ++n0;
    }
    std::cout << std::endl;
    for (unsigned j = 0; j < 15; j++) {
        if (n1 == _nBins) break;
        if (! _masks[n1]) printf("%4d", _bins[n1]);
        else              printf("-%3d", _bins[n1]);
        ++n1;
    }   
    std::cout << std::endl;
    }
 
    if (msg.find("detail") != std::string::npos) {
    for (unsigned i = 0; i < _nBins; i++) {
        std::cout << "bin " << i << " : ";
        for (unsigned j = 0; j < _links[i]->length(); j++) {
        std::cout << (* _links[i])[j]->wire()->name() << ",";
        }
        std::cout << std::endl;
    }
    }
 
    return;
}
 
void
THistogram::fillX(const AList<TMLink> & links) {
    _all = (AList<TMLink> &) links;
    unsigned nLinks = links.length();
    double offset = _binSize / 4.;
    for (unsigned i = 0; i < nLinks; i++) {
    TMLink * l = links[i];
    const HepPoint3D & p = l->position();
    unsigned pos = (unsigned) floor((p.x() + offset) / _binSize);
 
    //...Why is this needed?...
    pos %= _nBins;
 
    ++_bins[pos];
    _links[pos]->append(l);
    }
}
 
void
THistogram::fillY(const AList<TMLink> & links) {
    _all = (AList<TMLink> &) links;
    unsigned nLinks = links.length();
    for (unsigned i = 0; i < nLinks; i++) {
    TMLink * l = links[i];
    const HepPoint3D & p = l->position();
    unsigned pos = (unsigned) floor(p.y() / _binSize);
 
    //...Why is this needed?...
    pos %= _nBins;
 
    ++_bins[pos];
    _links[pos]->append(l);
    }
}
 
void
THistogram::fillPhi(const AList<TMLink> & links) {
    _all = (AList<TMLink> &) links;
    unsigned nLinks = links.length();
    double offset = _binSize / 4.;
    for (unsigned i = 0; i < nLinks; i++) {
    TMLink * l = links[i];
    const HepPoint3D & p = l->position();
    float phi = atan2(p.y(), p.x()) + M_PI;
//  std::cout<<"atan "<<atan2(-1., -1.)<<std::endl;
    unsigned pos = (unsigned) floor((phi + offset) / _binSize);
 
//  std::cout <<"layer "<<l->wire()->layerId()<<" cell "<<l->wire()->localId()
//        <<" x "<<p.x()<<" y "<<p.y()<<" pos "<<pos<<" xypos hit "<<l->hit()->xyPosition()<<" xypos wire "<<l->wire()->xyPosition()<<" drift "<<l->hit()->drift(0)<< "+-" <<l->hit()->dDrift(0)<<std::endl;
//  std::cout<<"binsize "<<_binSize<<" offset "<<offset<<" phi "<<phi<<std::endl;
    //...Why is this needed?...
    pos %= _nBins;
 
    ++_bins[pos];
    _links[pos]->append(l);
    }
}
 
void
THistogram::remove(const AList<TMLink> & links) {
    for (unsigned i = 0; i < _nBins; i++) {
    _links[i]->remove(links);
    _bins[i] = _links[i]->length();
    }
    _all.remove(links);
}
 
AList<TMLink>
THistogram::contents(unsigned center, unsigned width) const {
    AList<TMLink> links;
    for (int i = - (int) width;
     i <= (int) width;
     i++) {
    links.append(* bin((int) center + i));
    }
    return links;
}
 
AList<TMLink>
THistogram::contents(int start, int end) const {
    AList<TMLink> links;
    for (int i = start; i <= end; i++)
    links.append(* bin(i));
    return links;
}
 
AList<TSegment0>
THistogram::clusters0(void) const {
    AList<TSegment0> list;
 
    //...Serach for empty bin...
    unsigned begin = 0;
    while (_bins[begin] > 0) begin++;
    if (begin == _nBins) return list;
 
    //...Start searching...
    unsigned loop = 0;
    while (loop < _nBins) {
    ++loop;
    unsigned id = (begin + loop) % _nBins;
    if (_bins[id]) {
        unsigned size = 0;
        TSegment0 * c = new TSegment0();
        while (_bins[id]) {
        if (_bins[id]) ++size;
        c->append(* _links[id]);
        ++loop;
        id = (begin + loop) % _nBins;
        if (loop == _nBins) break;
        }
        list.append(c);
    }
    }
    return list;
}
 
AList<TSegment0>
THistogram::segments0(void) const {
 
    //...Obtain raw clusters...
    AList<TSegment0> list = clusters0();
    unsigned n = list.length();
    if (n == 0) return list;
 
    //...Examine each cluster...
    AList<TSegment0> splitted;
    for (unsigned i = 0; i < n; i++) {
    TSegment0 * c = list[i];
 
    AList<TSegment0> newClusters = c->split();
    if (newClusters.length() == 0) {
        c->solveDualHits();
        continue;
    }
 
    list.append(newClusters);
    splitted.append(c);
#ifdef TRKRECO_DEBUG_DETAIL
    c->dump("hits", "    ");
    std::cout << "    ... splitted as" << std::endl;
    for (unsigned j = 0; j < newClusters.length(); j++) {
        std::cout << "    " << j << " : ";
        newClusters[j]->dump("hits");
    }
#endif
    }
    list.remove(splitted);
    HepAListDeleteAll(splitted);
 
    return list;
 
}
 
AList<TSegment>
THistogram::clusters(void) const {
    AList<TSegment> list;
//    std::cout<<"enter clusters"<<std::endl;
    //...Serach for empty bin...
    unsigned begin = 0;
    while (_bins[begin] > 0)begin++;
    if (begin == _nBins) return list;
 
    //...Start searching...
    unsigned loop = 0;
    while (loop < _nBins) {
    ++loop;
    unsigned id = (begin + loop) % _nBins;
    if (_bins[id]) {
        unsigned size = 0;
        TSegment * c = new TSegment();
        while (_bins[id]) {
        if (_bins[id]) ++size;
        c->append(* _links[id]);
        ++loop;
        id = (begin + loop) % _nBins;
        if (loop == _nBins) break;
        }
        list.append(c);
    }
    }
    return list;
}
 
AList<TSegment>
THistogram::segments(void) const {
  //yuany
  /*
  for (unsigned i = 0; i < _nBins; i++) {
      std::cout<<"i "<<i<<" bins[i] "<<_bins[i]<<std::endl;
    }
  */
    //...Obtain raw clusters...
    AList<TSegment> list = clusters();
    unsigned n = list.length();
    if (n == 0) return list;
 
    //...Examine each cluster...
    AList<TSegment> splitted;
    for (unsigned i = 0; i < n; i++) {
    TSegment * c = list[i];
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    base segment : ";
    c->dump("hits");
#endif
 
    AList<TSegment> newClusters = c->split();
    if (newClusters.length() == 0) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "    ... Solving dual hits" << std::endl;
#endif
        c->solveDualHits();
        continue;
    }
 
    list.append(newClusters);
    splitted.append(c);
#ifdef TRKRECO_DEBUG_DETAIL
    c->dump("hits", "    ");
    std::cout << "    ... splitted as" << std::endl;
    for (unsigned j = 0; j < newClusters.length(); j++) {
        std::cout << "    " << j << " : ";
        newClusters[j]->dump("hits");
    }
#endif
    }
    list.remove(splitted);
    HepAListDeleteAll(splitted);
    
    //yuany
//    n = list.length();
//    for (unsigned i = 0; i < n; i++) {
//  TSegment * c = list[i];
//  std::cout << "    base segment : ";
//  c->dump("hits");
//    }
 
    return list;
}