/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/TrkReco/TrkReco-00-08-59-patch4-slc6tag/src/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    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// 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;
}

---