RowBasedFastHoughTransformationProcessor.cc

Go to the documentation of this file.

// #ifdef USE_ROWFHT
#include "RowBasedFastHoughTransformationProcessor.h"

#include "SimpleHelixTrackModel.h"

#include <iostream>
#include <string>
#include <ctime>
#include <cmath>

//GEAR
#include <gear/TPCParameters.h>
#include <gear/BField.h>
#include <gearimpl/Vector3D.h>

// include what you need from lcio
#include <EVENT/LCEvent.h>
#include <IMPL/LCCollectionVec.h>
#include <UTIL/BitField64.h>
#include <UTIL/ILDConf.h>
#include <IMPL/TrackerHitImpl.h>
#include <IMPL/TrackImpl.h>
#include <IMPL/LCFlagImpl.h>
#include <Exceptions.h>

// this is the inserted part

using namespace lcio;
using namespace marlin;
using namespace std;
using namespace gear;

namespace marlintpc {

static bool compare(pair<unsigned int, float> one, pair<unsigned int, float> two) {
        return ((one.second) > (two.second));
}

RowBasedFastHoughTransformationProcessor aRowBasedFastHoughTransformationProcessor;

RowBasedFastHoughTransformationProcessor::RowBasedFastHoughTransformationProcessor() :
                Processor("RowBasedFastHoughTransformationProcessor") {
        // the processor description
        _description = "RowBasedFastHoughTransformationProcessor is intended to ...";

        registerInputCollection(LCIO::TRACKERHIT, "InputHits", "Name of the input Hit Collection", _inputColName, string("TPCHits"));
        registerOutputCollection(LCIO::TRACK, "OutputTracks", "Name of the output Track Collection", _outputColName, string("FHTTracks"));

        registerOptionalParameter("BFieldScaleFactor", "Scales magnetic field (map), use 1.0 (default) for field ON or 0.0 for field OFF",
                        _bfieldScaleFactor, double(1.0));

        registerProcessorParameter("CenterXYMeasurement", "Center of XY (anode plane) measurement", _centerXY, double(60.));

        registerProcessorParameter("RangeXYMeasurement", "Range of XY (anode plane) measurement", _rangeXY, double(200.));

        registerProcessorParameter("CenterXZMeasurement", "Center of XZ (drift time) measurement", _centerXZ, double(300.));

        registerProcessorParameter("RangeXZMeasurement", "Range of XZ (drift time) measurement", _rangeXZ, double(600.));

        registerOptionalParameter("UseXYMeasurement", "Use XY (anode plane) measurement", _useXY, bool(true));

        registerOptionalParameter("UseXZMeasurement", "Use XZ (drift time) measurement", _useXZ, bool(true));

        registerOptionalParameter("MinimumRows", "Minimal number of (correlated) rows ", _minRow, int(8));

        registerOptionalParameter("MaximumRowDifference", "Row correation parameter ", _maxRowDiff, int(1));

        registerOptionalParameter("FractionOfRows", "Relative minimum cube content (rows >= fraction * (correlated) rows) ", _fracRow,
                        float(0.8));

        registerOptionalParameter("MaximumCubes", "Maximum number of (hyper) cubes to check ", _maxCube, int(250));

        registerOptionalParameter("MinimumLevel", "Minimum number of (hyper) cubes splittings ", _minLevel, int(5));

        registerOptionalParameter("MaximumLevel", "Maximum number of (hyper) cubes splittings ", _maxLevel, int(8));

        registerOptionalParameter("EfficienyCut", "Minimum hit density (hits/tracklength) ", _effCut, float(0.8));

        registerOptionalParameter("PurityCut", "Maximum hit density (hits/tracklength) ", _purCut, float(1.1));

        registerOptionalParameter("MaxHitsPerRow", "Maximum number of hits per row ", _maxHitsPerRow, int(1));

        registerOptionalParameter("UnusedHitMatchingChi2Cut", "Chi2 cut for matching unused hits in XY and Z", _unusedCut, double(20.0));

        registerOptionalParameter("UnusedHitMaxGapCut", "Maximal (row) gap for matching unused hits in XY and Z", _maxGap, int(4));

        registerOptionalParameter("EncodedModuleID", "Module ID encoded in CellID0", _encodedModuleID, bool(true));

        registerOptionalParameter("ReferencePointAtPca", "Use PCA as reference point, else 1. hit", _refAtPCA, bool(false));

        registerOptionalParameter("TimePixFlag", "Is timepix data", _timePixFlag, bool(false));

}

void RowBasedFastHoughTransformationProcessor::init() {
        m_out(DEBUG) << " init() called  " << endl;

}

void RowBasedFastHoughTransformationProcessor::processRunHeader(LCRunHeader* run) {
        // do something for each run
}

void RowBasedFastHoughTransformationProcessor::processEvent(LCEvent* evt) {
        std::clock_t start1 = std::clock();

        m_out(DEBUG) << "Event Number: " << evt->getEventNumber() << endl;

        //get input collection
        LCCollection* inputCol = 0;

        LCCollectionVec* outputCol = new LCCollectionVec(LCIO::TRACK);
        LCFlagImpl trkFlag(0);
        trkFlag.setBit(LCIO::TRBIT_HITS);
        outputCol->setFlag(trkFlag.getFlag());

        // do a try / catch in case there is an empty event
        try {
                inputCol = evt->getCollection(_inputColName);
        } catch (DataNotAvailableException &e) {
                // no input data, just return and process the next event
                return;
        }

        Vector3D theTPCCenter(_Xcenter, _Ycenter, 0.);
        Vector3D bfield = marlin::Global::GEAR->getBField().at(theTPCCenter);

        const bool bOff(bfield.z() * _bfieldScaleFactor == 0.);
        if (bOff) {
                m_out(DEBUG) << "BField(GEAR) is OFF" << std::endl;
        } else {
                m_out(DEBUG) << "BField(GEAR) is ON" << std::endl;
        }
        // constant magnetic field (in z direction)
        _Bzc = bfield.z() * _bfieldScaleFactor * 0.0002998; // Bz*c

        // setup (number of parameters per measurement direction)
        unsigned int setup[] = { 3, 2 };
        if (bOff) setup[0] = 2;
        if (!_useXY) setup[0] = 0;
        if (!_useXZ) setup[1] = 0;
        if (setup[0] + setup[1] == 0) {
                m_out(ERROR) << " No active parameters ";
                return;
        }
        unsigned int dimension = (_useXY and _useXZ) ? 2 : 1; // dimension of measurements
        _scaleBits = _maxLevel + (_maxLevel % 2) + 8; // even number

        rowHitMapType hitData;

        UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
        int nInput = inputCol->getNumberOfElements();
        m_out(DEBUG) << "Event " << evt->getEventNumber() << " in run " << evt->getRunNumber() << ", number of input hits: " << nInput
                        << endl;

        if (_timePixFlag) {
                for (int iInput = 0; iInput < nInput; iInput++) {
                        TrackerHit* Hit = dynamic_cast<TrackerHit*>(inputCol->getElementAt(iInput));
                        encoder.setValue(Hit->getCellID0());
                        int mod = Hit->getCellID1();
                        int rowLoc = Hit->getCellID0() % 256 + mod * 256;
                        int row = (rowLoc > 1023) ? rowLoc - 1024 : 1023 - rowLoc; // simple row count in x direction
                        rb_Hit* someHit = new rb_Hit(iInput, mod, row, *Hit);
                        hitData[row].push_back(someHit);
                }
        } else {
                for (int iInput = 0; iInput < nInput; iInput++) {
                        TrackerHit* Hit = dynamic_cast<TrackerHit*>(inputCol->getElementAt(iInput));
                        encoder.setValue(Hit->getCellID0());
                        int mod =
                                        _encodedModuleID ?
                                                        encoder[ILDCellID0::module] :
                                                        Global::GEAR->getTPCParameters().getNearestModule(Hit->getPosition()[0], Hit->getPosition()[1]).getModuleID();
                        int row = encoder[ILDCellID0::layer];
                        rb_Hit* someHit = new rb_Hit(iInput, mod, row, *Hit);
                        hitData[row].push_back(someHit);
                }
        }

        // list of track candidates (segments)
        segListType trackCandidates;
        double scale = double(1 << _scaleBits);

        bool more = true;
        while (more) {
                int numRows = 0;
                int numRows2 = 0;
                int lastRow = -1;
                double xMin = 0., xMax = 0.;
                for (rowHitMapType::iterator itRow = hitData.begin(); itRow != hitData.end(); ++itRow) {
                        const int row = itRow->first;
                        //std::cout << " row " << row << std::endl;
                        for (hitListType::iterator itHit = itRow->second.begin(); itHit != itRow->second.end(); ++itHit) {
                                if (!(*itHit)->getUsed()) {
                                        double x = (*itHit)->getX();
                                        if (numRows == 0) {
                                                xMax = xMin = x;
                                                numRows++;
                                                numRows2++;
                                        } else {
                                                xMax = max(xMax, x);
                                                xMin = min(xMin, x);
                                                if (row > lastRow) {
                                                        numRows++;
                                                        if (row <= lastRow + _maxRowDiff) numRows2++;
                                                }
                                        }
                                        lastRow = row;
                                }
                        }
                }
                m_out(DEBUG) << " ... rows " << numRows << " " << numRows2 << endl;
                if (numRows2 < _minRow) break;
                const double xCenter = 0.5 * (xMin + xMax);
                const double xRange = 0.5 * (xMax - xMin);
                intListType rootDistances;
                rootDistances.reserve(nInput * dimension);
                hyperPlaneListType rootPlanes;
                rootPlanes.reserve(nInput);
                //

                for (rowHitMapType::iterator itRow = hitData.begin(); itRow != hitData.end(); ++itRow) {
                        //int row = itRow->first;
                        //std::cout << " row " << row << std::endl;
                        for (hitListType::iterator itHit = itRow->second.begin(); itHit != itRow->second.end(); ++itHit) {
                                if (!(*itHit)->getUsed()) {
                                        const double u = ((*itHit)->getX() - xCenter) / xRange;
                                        const double v[] = { ((*itHit)->getY() - _centerXY) / _rangeXY, ((*itHit)->getZ() - _centerXZ) / _rangeXZ };
                                        // calculate unit normals
                                        directionsType dirs;
                                        bool inside = true;
                                        int dist[2];
                                        unsigned int ndim = 0;
                                        for (unsigned int idim = 0; idim < 2; idim++) {
                                                if (setup[idim] > 0) {
                                                        doubleListType vec(setup[idim]);
                                                        vec[0] = 1.;
                                                        vec[1] = u;
                                                        if (setup[idim] > 2) vec[2] = 1.5 * u * u - 0.5;
                                                        double sum2 = 0.;
                                                        for (doubleListType::iterator i = vec.begin(); i != vec.end(); ++i)
                                                                sum2 += (*i) * (*i);
                                                        const double norm = sqrt(sum2);
                                                        for (doubleListType::iterator i = vec.begin(); i != vec.end(); ++i)
                                                                *i /= norm;
                                                        dirs.push_back(vec);
                                                        double r = -v[idim] / norm;
                                                        inside = inside and (abs(r * vec[0]) < 0.5);
                                                        dist[ndim++] = int(scale * r);
                                                }
                                        }
                                        if (inside) {
                                                rb_HyperPlane* hp = new rb_HyperPlane(_scaleBits, *itHit, dirs);
                                                for (unsigned int idim = 0; idim < ndim; idim++)
                                                        rootDistances.push_back(dist[idim]);
                                                rootPlanes.push_back(hp);
                                        }
                                }
                        }
                }
                m_out(DEBUG) << "rootPlanes " << rootPlanes.size() << endl;
                if (rootPlanes.size() < static_cast<unsigned int>(_minRow)) break;

                rb_HyperCube root(setup, rootPlanes, rootDistances);
                int minCubeRows = int(float(numRows2) * _fracRow);
                m_out(DEBUG) << "minCubeRows " << minCubeRows << endl;
                hitListType candidate = root.divide(_maxCube, minCubeRows, _minLevel, _maxLevel, _effCut, _purCut);
                //root.print();
                m_out(DEBUG) << "Candidate " << candidate.size() << endl;
                more = false;
                if (candidate.size() > 0) {
                        unsigned int maxHitsPerRow = _maxHitsPerRow;
                        rb_Segment* aSegment = new rb_Segment(_Bzc, candidate, maxHitsPerRow);
                        if (aSegment->getNdf() < 0)
                                delete aSegment;
                        else {
                                if (_unusedCut > 0.) {
                                        // link unused hits
                                        map<int, int> rowMap;
                                        aSegment->fillRowMap(rowMap);
                                        int centerRow = (aSegment->getFirstRow() + aSegment->getLastRow()) / 2;
                                        // look forward
                                        for (int row = centerRow + 1; row <= hitData.rbegin()->first; ++row) {
                                                if (row > aSegment->getLastRow() + _maxGap) break; // too far away
                                                if (rowMap.count(row) > 0) continue; // row already contained in segment
                                                hitListType hits;
                                                aSegment->match(hitData[row], hits, _unusedCut); // match list of hits
                                                if (hits.size() <= maxHitsPerRow) for (hitListType::iterator itHit = hits.begin(); itHit != hits.end(); ++itHit)
                                                        aSegment->addHit(*itHit); // additional hit(s)?
                                        }
                                        // look backward
                                        for (int row = centerRow; row >= hitData.begin()->first; --row) {
                                                if (row < aSegment->getFirstRow() - _maxGap) break; // too far away
                                                if (rowMap.count(row) > 0) continue; // row already contained in segment
                                                hitListType hits;
                                                aSegment->match(hitData[row], hits, _unusedCut); // match list of hits
                                                if (hits.size() <= maxHitsPerRow) for (hitListType::iterator itHit = hits.begin(); itHit != hits.end(); ++itHit)
                                                        aSegment->addHit(*itHit); // additional hit(s)?
                                        }
                                }
                                trackCandidates.push_back(aSegment); // keep valid segments
                                m_out(DEBUG) << "Segment " << trackCandidates.size() << " " << aSegment->getHitList().size() << endl;
                                more = true;
                        }
                }
                // cleanup
                for (hyperPlaneListType::iterator itHp = rootPlanes.begin(); itHp != rootPlanes.end(); ++itHp)
                        delete *itHp;
                rootPlanes.clear();
                rootDistances.clear();
        }

        m_out(DEBUG) << "  found " << trackCandidates.size() << " track candidates " << endl;

        for (segListType::iterator itSeg = trackCandidates.begin(); itSeg != trackCandidates.end(); ++itSeg) {
                // new track
                TrackImpl* theTrack = new TrackImpl();
                // get hits
                hitListType hits = (*itSeg)->getHitList();
                // position of 1. hit
                double posFirst[3];
                hits.front()->getPos(posFirst);
                // reference point, 1. hit or PCA (point of closest approach (to origin))
                double refPos[] = { 0., 0., 0. };
                if (!_refAtPCA) {
                        refPos[0] = posFirst[0];
                        refPos[1] = posFirst[1];
                        refPos[2] = posFirst[2];
                }
                float refPoint[3];
                for (int i = 0; i < 3; ++i)
                        refPoint[i] = refPos[i];
                // track parameters and covariance matrix
                TVectorD parameters(5);
                TMatrixDSym covariance(5);
                (*itSeg)->getLCIOStateAtRefPoint(refPos, parameters, covariance); // ( d0, phi0, Omega, z0, tanLambda)

                theTrack->setOmega(parameters[2]);
                theTrack->setTanLambda(parameters[4]);
                theTrack->setPhi(parameters[1]);
                theTrack->setD0(parameters[0]);
                theTrack->setZ0(parameters[3]);
                theTrack->setChi2((*itSeg)->getChi2());
                theTrack->setNdf((*itSeg)->getNdf());
                theTrack->setReferencePoint(refPoint);
                /* prepare LCIO output; unknown parts:
                 *
                 * theTrack->setTypeBit (int index, bool val=true) // ?
                 * theTrack->setdEdx (float dEdx)
                 * theTrack->setdEdxError (float dEdxError) */
                float rInner = sqrt(pow(posFirst[0], 2) + pow(posFirst[1], 2));
                theTrack->setRadiusOfInnermostHit(rInner);
                // compressed covariance matrix
                float compressedCovariance[15];
                int ind = 0;
                for (int i = 0; i < 5; ++i)
                        for (int j = 0; j <= i; ++j)
                                compressedCovariance[ind++] = covariance[i][j];
                theTrack->setCovMatrix(compressedCovariance);
                // and finally add the hits to the track
                for (hitListType::iterator itH = hits.begin(); itH != hits.end(); ++itH) {
                        TrackerHit* Hit = dynamic_cast<TrackerHit*>(inputCol->getElementAt((*itH)->getHitNum()));
                        theTrack->addHit(Hit);
                }

                m_out(DEBUG) << " +++ the found track has a chi^2 value of " << theTrack->getChi2() << " for ndf = " << theTrack->getNdf()
                                << " and the parameters [omega, tanLambda, phi0, d0, z0] = [" << theTrack->getOmega() << "," << theTrack->getTanLambda()
                                << "," << theTrack->getPhi() << "," << theTrack->getD0() << "," << theTrack->getZ0() << "] at the reference position ["
                                << refPoint[0] << "," << refPoint[1] << "," << refPoint[2] << "]." << endl;

                outputCol->addElement(theTrack);
        }

        m_out(DEBUG) << "RowBasedFastHoughTransformation ... done!" << endl;

        if (!outputCol->empty()) {
                evt->addCollection(outputCol, _outputColName);
        }

        m_out(DEBUG) << ((std::clock() - start1) / (double) CLOCKS_PER_SEC) << " sec" << endl;

        // cleanup
        for (rowHitMapType::iterator itRow = hitData.begin(); itRow != hitData.end(); ++itRow) {
                for (hitListType::iterator itH = itRow->second.begin(); itH != itRow->second.end(); ++itH) {
                        delete *itH;
                }
        }
        hitData.clear();

        // get rid of segments
        for (segListType::iterator itSeg = trackCandidates.begin(); itSeg != trackCandidates.end(); ++itSeg) {
                delete *itSeg;
        }
        trackCandidates.clear();

}

void RowBasedFastHoughTransformationProcessor::check(LCEvent* evt) {
        // check what is needed here
}

void RowBasedFastHoughTransformationProcessor::end() {
        // do something at the end of the processing, e.g. cleaning up

}


rb_HyperPlane::rb_HyperPlane(int scaleBits, rb_Hit* hit, directionsType& dirList, double distCut) :
                _scaleBits(scaleBits), _hit(hit), _cut(), _steps() {
        double scale = double(1 << _scaleBits);
        for (directionsType::iterator itDir = dirList.begin(); itDir != dirList.end(); ++itDir) {
                _cut.push_back(int(scale * distCut / (*itDir)[0])); // first is max component
                const unsigned int nSteps = 1 << itDir->size();
                intListType aSteps(nSteps);
                for (unsigned int i = 0; i < nSteps; i++) {
                        unsigned int ii = i;
                        double step = 0.;
                        for (doubleListType::iterator itNum = itDir->begin(); itNum != itDir->end(); ++itNum) {
                                step += (ii % 2 == 1) ? 0.5 * (*itNum) : -0.5 * (*itNum);
                                ii >>= 1;
                        }
                        aSteps[i] = int(scale * step);
                }
                _steps.push_back(aSteps);
        }
}

void rb_HyperPlane::print() const {
        cout << " HyperPlane " << _hit->getHitNum() << " " << _hit->getRow() << " " << _cut.size() << " " << _cut.front() << " "
                        << _cut.back() << endl;
}

int rb_HyperPlane::getRow() const {
        return _hit->getRow();
}

int rb_HyperPlane::getScaleBits() const {
        return _scaleBits;
}

rb_Hit* rb_HyperPlane::getHit() const {
        return _hit;
}

intListType rb_HyperPlane::getCut() const {
        return _cut;
}

int rb_HyperPlane::getCut(int idim) const {
        return _cut[idim];
}

std::vector<intListType> rb_HyperPlane::getSteps() const {
        return _steps;
}

int rb_HyperPlane::getStep(int idim, int b) const {
        return _steps[idim][b];
}


rb_HyperCube::rb_HyperCube(unsigned int* setup, hyperPlaneListType& planes, intListType& distances, unsigned int level) :
                _planes(planes), _distances(distances), _level(level) {
        _setup[0] = setup[0];
        _setup[1] = setup[1];
        _numChild = 1 << (setup[0] + setup[1]);
        _dimension = (setup[0] > 0 and setup[1] > 0) ? 2 : 1;
        if (level == 0) nCube = 0;
        nCube++;
}

void rb_HyperCube::print() const {
        cout << " HyperCube " << nCube << " " << _level << " " << _numChild << endl;
}


hitListType rb_HyperCube::divide(int maxCube, int minRows, int minLevel, int maxLevel, float effCut, float purCut) {
        vector<int> numHits(_numChild, 0);
        vector<int> numRows(_numChild, 0);
        vector<int> lastRow(_numChild, -1);
        vector<int> sum2(_numChild, 0);

        hitListType hits(0);

        int scaleBits = _planes.front()->getScaleBits();
        int halfScale = scaleBits / 2;
        int nChild0 = 1 << _setup[0];
        int nChild1 = 1 << _setup[1];
        unsigned int iOff = 0;
        for (hyperPlaneListType::iterator itHp = _planes.begin(); itHp != _planes.end(); ++itHp) {
                rb_HyperPlane* hp = *itHp;
                int row = hp->getRow();
                // factorized subspaces?
                if (_dimension == 2) {
                        // two factorized subspaces
                        vector<pair<int, int> > distNew1;
                        distNew1.reserve(nChild1);
                        for (int b1 = 0; b1 < nChild1; b1++) {
                                int distChild1 = _distances[iOff + 1] * 2 + hp->getStep(1, b1);
                                if (abs(distChild1) < hp->getCut(1)) distNew1.push_back(make_pair(b1, distChild1));
                        }
                        for (int b0 = 0; b0 < nChild0; b0++) {
                                int distChild0 = _distances[iOff] * 2 + hp->getStep(0, b0);
                                if (abs(distChild0) < hp->getCut(0)) {
                                        distChild0 >>= halfScale;
                                        for (vector<pair<int, int> >::iterator itP = distNew1.begin(); itP != distNew1.end(); ++itP) {
                                                int distChild1 = itP->second >> halfScale;
                                                int b = b0 + (itP->first << _setup[0]);
                                                numHits[b]++;
                                                if (row > lastRow[b]) {
                                                        lastRow[b] = row;
                                                        numRows[b]++;
                                                }
                                                sum2[b] += distChild0 * distChild0 + distChild1 * distChild1;
                                        }
                                }
                        }
                } else {
                        // single subspace
                        for (int b = 0; b < _numChild; b++) {
                                int distChild = _distances[iOff] * 2 + hp->getStep(0, b);
                                if (abs(distChild) < hp->getCut(0)) {
                                        distChild >>= halfScale;
                                        numHits[b]++;
                                        if (row > lastRow[b]) {
                                                lastRow[b] = row;
                                                numRows[b]++;
                                        }
                                        sum2[b] += distChild * distChild;
                                }
                        }
                }
                iOff += _dimension;
        }
        vector<pair<unsigned int, int> > childList;
        for (int b = 0; b < _numChild; b++) {
                //cout << " child " << _level << " " << b << " " << numRows[b] << " " << numHits[b] << endl;
                if (numRows[b] >= minRows) {
                        int weight = (numRows[b] << scaleBits) - (sum2[b] / numHits[b]);
                        childList.push_back(make_pair(b, weight));
                }
        }
        sort(childList.begin(), childList.end(), compare);
        for (vector<pair<unsigned int, int> >::iterator itP = childList.begin(); itP != childList.end(); ++itP) {
                unsigned int b = itP->first;
                // calculate child planes and distances
                hyperPlaneListType childPlanes;
                childPlanes.reserve(numHits[b]);
                intListType childDistances;
                childDistances.reserve(numHits[b] * _dimension);
                int b0 = b % (1 << _setup[0]);
                int b1 = b >> _setup[0];
                iOff = 0;
                for (hyperPlaneListType::iterator itHp = _planes.begin(); itHp != _planes.end(); ++itHp) {
                        rb_HyperPlane* hp = *itHp;
                        if (_dimension == 2) {
                                int distNew0 = _distances[iOff] * 2 + hp->getStep(0, b0);
                                int distNew1 = _distances[iOff + 1] * 2 + hp->getStep(1, b1);
                                if (abs(distNew0) < hp->getCut(0) and abs(distNew1) < hp->getCut(1)) {
                                        childPlanes.push_back(hp);
                                        childDistances.push_back(distNew0);
                                        childDistances.push_back(distNew1);
                                }
                        } else {
                                int distNew0 = _distances[iOff] * 2 + hp->getStep(0, b0);
                                if (abs(distNew0) < hp->getCut(0)) {
                                        childPlanes.push_back(hp);
                                        childDistances.push_back(distNew0);
                                }
                        }
                        iOff += _dimension;
                }
                // calculate hit density (number of hits / length (in rows))
                int rowLength = childPlanes.back()->getRow() - childPlanes.front()->getRow() + 1;
                float density = float(numHits[b]) / float(rowLength);
                //cout << " sorted child cubes " << itP->first << " " << itP->second << " " << rowLength << " " << density << endl;
                // minimum level reached ?
                if (_level >= minLevel) {
                        if ((effCut < density) && (density < purCut)) {
                                for (hyperPlaneListType::iterator itHp = childPlanes.begin(); itHp != childPlanes.end(); ++itHp)
                                        hits.push_back((*itHp)->getHit());
                                return hits;
                        }
                }
                if (_level < maxLevel) {
                        rb_HyperCube newCube(_setup, childPlanes, childDistances, _level + 1);
                        if (nCube > maxCube) return hits;
                        hitListType candidate = newCube.divide(maxCube, minRows, minLevel, maxLevel, effCut, purCut);
                        if (candidate.size() > 0) return candidate;
                }
        }
        return hits;
}

} // close namespace brackets

// #endif //#ifdef USE_ROWFHT