RowTripletBasedTrackFinderProcessor.cc

Go to the documentation of this file.

// #ifdef USE_ROWTRIPLETS
#include "RowTripletBasedTrackFinderProcessor.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 {

RowTripletBasedTrackFinderProcessor aRowTripletBasedTrackFinderProcessor;

RowTripletBasedTrackFinderProcessor::RowTripletBasedTrackFinderProcessor() :
                Processor("RowTripletBasedTrackFinderProcessor") {
        // the processor description
        _description = "RowTripletBasedTrackFinderProcessor 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("TripletTracks"));

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

        registerOptionalParameter("TripletPreSelectionDistCut", "Coarse cut on XY and Z residuals for triplet preselection", _distCut,
                        double(10.0));

        registerOptionalParameter("TripletDefinitionChi2Cut", "Chi2 cut for triplet definition on XY and Z residuals", _trpCut,
                        double(20.0));

        registerOptionalParameter("TripletPositionMatchingChi2Cut", "Chi2 cut for triplet position matching in XY and Z", _posCut,
                        double(20.0));

        registerOptionalParameter("TripletDirectionMatchingChi2Cut", "Chi2 cut for triplet direction matching in XY and Z", _dirCut,
                        double(20.0));

        registerOptionalParameter("SegmentMatchingChi2Cut", "Chi2/Ndf cut for segment matching (Ndf=4 for B off, 5 for B on)", _segCut,
                        double(30.0));

        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));
}

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

        const gear::TPCParameters& tpcParameters = Global::GEAR->getTPCParameters();
        bool cartesian = (tpcParameters.getCoordinateType() == PadRowLayout2D::CARTESIAN);
        const std::vector<TPCModule*> moduleVector = tpcParameters.getModules();
        m_out(DEBUG) << " Number of modules: " << moduleVector.size() << endl;
        for (std::vector<TPCModule*>::const_iterator itMod = moduleVector.begin(); itMod != moduleVector.end(); ++itMod) {
                int module = (*itMod)->getModuleID();
                m_out(DEBUG) << " module " << module << "  with " << (*itMod)->getNRows() << " rows" << endl;
                for (std::vector<TPCModule*>::const_iterator itMod2 = itMod; itMod2 != moduleVector.end(); ++itMod2) {
                        if (areNeighbourModules(*itMod, *itMod2, cartesian)) {
                                _modNeighbours[module].push_back((*itMod2)->getModuleID());
                                m_out(DEBUG) << "   neighbour " << (*itMod2)->getModuleID() << endl;
                        }
                }
        }

        // TPC center
        _Xcenter = tpcParameters.getDoubleVal("TPC_cylinder_centre_c0");
        _Ycenter = tpcParameters.getDoubleVal("TPC_cylinder_centre_c1");
}

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

void RowTripletBasedTrackFinderProcessor::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

        modHitMapType 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;
        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[mod][row].push_back(someHit);
        }

        // map of modules with vectors of segments
        modSegMapType segData;
        const int step = 2; // row step size for triplets
        rowHitMapType::iterator itRow1, itRow2, itRow3;
        for (modHitMapType::iterator itMod = hitData.begin(); itMod != hitData.end(); ++itMod) {
                const int module = itMod->first;
                //std::cout << " Module " << module << std::endl;
                rowHitMapType modData = itMod->second; // hits in module
                const int numRows = modData.rbegin()->first - modData.begin()->first + 1;
                trpListType triplets;
                m_out(DEBUG) << "  module " << module << " with " << endl;
                m_out(DEBUG) << "    " << modData.size() << " rows " << endl;
                //cout << " row range " << modData.begin()->first << " " << modData.rbegin()->first << endl;
                // loop through row triplet
                for (int row = modData.begin()->first + step; row <= modData.rbegin()->first - step; ++row) {
                        // all 3 rows must have hits to build triplets
                        if ((itRow1 = modData.find(row - step)) == modData.end()) continue;
                        if ((itRow2 = modData.find(row)) == modData.end()) continue;
                        if ((itRow3 = modData.find(row + step)) == modData.end()) continue;
                        // 3-fold loop
                        for (hitListType::iterator itHit1 = itRow1->second.begin(); itHit1 != itRow1->second.end(); ++itHit1) {
                                for (hitListType::iterator itHit3 = itRow3->second.begin(); itHit3 != itRow3->second.end(); ++itHit3) {
                                        rb_Doublet someDoublet(*itHit1, *itHit3); // doublet from outer hits
                                        for (hitListType::iterator itHit2 = itRow2->second.begin(); itHit2 != itRow2->second.end(); ++itHit2) {
                                                if (someDoublet.match(*itHit2, _distCut, _trpCut)) triplets.push_back(new rb_Triplet(someDoublet, *itHit2));
                                        }
                                }
                        }
                }
                m_out(DEBUG) << "    " << triplets.size() << " triplets " << endl;

                segListType segments;
                // equivalence classes for triplet matching
                simpleEquiClasses trpEquiClass;
                // loop over triplet pairs, check matching
                for (trpListType::iterator itTrp1 = triplets.begin(); itTrp1 != triplets.end(); ++itTrp1) {
                        trpEquiClass.addIndex(itTrp1 - triplets.begin());
                        for (trpListType::iterator itTrp2 = itTrp1; itTrp2 != triplets.end(); ++itTrp2) {
                                if ((*itTrp2)->getRow() - (*itTrp1)->getRow() > 2) break; // max. row distance <= 2
                                // matching triplets ?
                                if ((*itTrp1)->match(*itTrp2, _posCut, _dirCut))
                                        trpEquiClass.addMatch(make_pair(itTrp1 - triplets.begin(), itTrp2 - triplets.begin()));
                        }
                }
                // get classes with indices
                std::map<int, std::vector<int> > trpClasses = trpEquiClass.getClasses();
                // order classes
                std::multimap<int, int> rowInfo;
                int row, nRows, firstRow, lastRow;
                for (std::map<int, std::vector<int> >::iterator itC = trpClasses.begin(); itC != trpClasses.end(); ++itC) {
                        nRows = firstRow = lastRow = 0;
                        for (std::vector<int>::iterator itInd = itC->second.begin(); itInd != itC->second.end(); ++itInd) {
                                row = triplets[*itInd]->getRow();
                                if (nRows == 0) {
                                        firstRow = lastRow = row;
                                        nRows++;
                                } else if (row > lastRow) {
                                        lastRow = row;
                                        nRows++;
                                }
                        }
                        if (nRows > 1) rowInfo.insert(make_pair(nRows * numRows + lastRow - firstRow, itC->first)); // use only classes with at least 2 triplets
                }
                // allocate hits to classes (segments) ordered by rowInfo ( ~ number of rows)
                for (std::multimap<int, int>::reverse_iterator rit = rowInfo.rbegin(); rit != rowInfo.rend(); ++rit) {
                        trpListType matchingTriplets;
                        for (std::vector<int>::iterator itInd = trpClasses[rit->second].begin(); itInd != trpClasses[rit->second].end(); ++itInd)
                                matchingTriplets.push_back(triplets[*itInd]);
                        rb_Segment* simpleSegment = new rb_Segment(_Bzc, matchingTriplets);
                        if (simpleSegment->getNdf() < 0)
                                delete simpleSegment;
                        else
                                segments.push_back(simpleSegment); // keep valid segments
                }
                m_out(DEBUG) << "    " << segments.size() << " segments " << endl;
                // get rid of triplets
                for (trpListType::iterator itTrp = triplets.begin(); itTrp != triplets.end(); ++itTrp) {
                        delete *itTrp;
                }
                triplets.clear();

                simpleEquiClasses segEquiClass;
                // loop over segment pairs, check matching
                for (segListType::iterator itSeg1 = segments.begin(); itSeg1 != segments.end(); ++itSeg1) {
                        segEquiClass.addIndex(itSeg1 - segments.begin());
                        for (segListType::iterator itSeg2 = itSeg1; itSeg2 != segments.end(); ++itSeg2) {
                                // matching segments ?
                                if ((*itSeg1)->match(*itSeg2, _segCut, step))
                                        segEquiClass.addMatch(make_pair(itSeg1 - segments.begin(), itSeg2 - segments.begin()));
                        }
                }
                // get classes with indices
                std::map<int, std::vector<int> > segClasses = segEquiClass.getClasses();
                for (std::map<int, std::vector<int> >::iterator itC = segClasses.begin(); itC != segClasses.end(); ++itC) {
                        //cout << " seg-link " << itC->first << ": " << itC->second.size() << endl;
                        segListType matchingSegments;
                        for (std::vector<int>::iterator itInd = itC->second.begin(); itInd != itC->second.end(); ++itInd)
                                matchingSegments.push_back(segments[*itInd]);
                        segData[module].push_back(new rb_Segment(matchingSegments)); // save combined segments
                }
                // get rid of (simple) segments
                for (segListType::iterator itSeg = segments.begin(); itSeg != segments.end(); ++itSeg) {
                        delete *itSeg;
                }
                segments.clear();

                for (segListType::iterator itSeg = segData[module].begin(); itSeg != segData[module].end(); ++itSeg) {
                        map<int, int> rowMap;
                        (*itSeg)->fillRowMap(rowMap);
                        int centerRow = ((*itSeg)->getFirstRow() + (*itSeg)->getLastRow()) / 2;
                        // look forward
                        for (int row = centerRow + 1; row <= modData.rbegin()->first; ++row) {
                                if (row > (*itSeg)->getLastRow() + _maxGap) break; // too far away
                                if (rowMap.count(row) > 0) continue; // row already contained in segment
                                hitListType hits;
                                (*itSeg)->match(modData[row], hits, _unusedCut); // match list of hits
                                //if (!hits.empty()) cout << " match-fwd " << row << " " <<(*itSeg)->getFirstRow() << " " <<(*itSeg)->getLastRow() << ": " << hits.size() << endl;
                                if (hits.size() == 1) (*itSeg)->addHit(hits[0]); // single additional hit?
                        }
                        // look backward
                        for (int row = centerRow; row >= modData.begin()->first; --row) {
                                if (row < (*itSeg)->getFirstRow() - _maxGap) break; // too far away
                                if (rowMap.count(row) > 0) continue; // row already contained in segment
                                hitListType hits;
                                (*itSeg)->match(modData[row], hits, _unusedCut); // match list of hits
                                //if (!hits.empty()) cout << " match-bwd " << row << " " <<(*itSeg)->getFirstRow() << " " <<(*itSeg)->getLastRow() << ": " << hits.size() << endl;
                                if (hits.size() == 1) (*itSeg)->addHit(hits[0]); // single additional hit?
                        }
                }
        }

        std::vector<std::pair<int, int> > segIndex;
        simpleEquiClasses segEquiClass;
        // create segment index
        for (modSegMapType::iterator itMod = segData.begin(); itMod != segData.end(); ++itMod) {
                segListType segments = itMod->second;
                for (segListType::iterator itSeg = segments.begin(); itSeg != segments.end(); ++itSeg) {
                        (*itSeg)->setId(segIndex.size());
                        segIndex.push_back(make_pair(itMod->first, itSeg - segments.begin()));
                        segEquiClass.addIndex((*itSeg)->getId());
                }
        }
        // check segment matching with neighbour modules
        for (modSegMapType::iterator itMod = segData.begin(); itMod != segData.end(); ++itMod) {
                std::vector<int> neighbours = _modNeighbours[itMod->first];
                segListType segments1 = itMod->second;
                for (std::vector<int>::iterator itN = neighbours.begin(); itN != neighbours.end(); ++itN) {
                        segListType segments2 = segData[*itN];
                        for (segListType::iterator itSeg1 = segments1.begin(); itSeg1 != segments1.end(); ++itSeg1) {
                                for (segListType::iterator itSeg2 = segments2.begin(); itSeg2 != segments2.end(); ++itSeg2) {
                                        // matching segments ?
                                        if ((*itSeg1)->match(*itSeg2, _segCut, step)) segEquiClass.addMatch(make_pair((*itSeg1)->getId(), (*itSeg2)->getId()));
                                }
                        }
                }
        }
        // get classes with indices
        std::map<int, std::vector<int> > segClasses = segEquiClass.getClasses();
        segListType trackCandidates;
        for (std::map<int, std::vector<int> >::iterator itC = segClasses.begin(); itC != segClasses.end(); ++itC) {
                //cout << " seg-link-m " << itC->first << ": " << itC->second.size() << endl;
                segListType matchingSegments;
                for (std::vector<int>::iterator itInd = itC->second.begin(); itInd != itC->second.end(); ++itInd) {
                        std::pair<int, int> index = segIndex[*itInd];
                        matchingSegments.push_back(segData[index.first][index.second]);
                }
                trackCandidates.push_back(new rb_Segment(matchingSegments)); // save combined segments
        }
        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) << "RowTripletBasedTrackFinder ... done!" << endl;

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

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

        // clean up
        for (modHitMapType::iterator itMod = hitData.begin(); itMod != hitData.end(); ++itMod) {
                for (rowHitMapType::iterator itRow = itMod->second.begin(); itRow != itMod->second.end(); ++itRow) {
                        for (hitListType::iterator itH = itRow->second.begin(); itH != itRow->second.end(); ++itH) {
                                delete *itH;
                        }
                }
        }
        hitData.clear();
        for (modSegMapType::iterator itMod = segData.begin(); itMod != segData.end(); ++itMod) {
                for (segListType::iterator itSeg = itMod->second.begin(); itSeg != itMod->second.end(); ++itSeg) {
                        delete *itSeg;
                }
        }
        segData.clear();
        for (segListType::iterator itSeg = trackCandidates.begin(); itSeg != trackCandidates.end(); ++itSeg) {
                delete *itSeg;
        }
        trackCandidates.clear();

}

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

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

}


bool RowTripletBasedTrackFinderProcessor::areNeighbourModules(gear::TPCModule* mod1, gear::TPCModule* mod2, bool cartesian) {
        // are mod1 and mod2 neighbours?
        if (mod1 == mod2) return false;
        const std::vector<double> modExtend = mod1->getModuleExtent();
        const double xm1 = 0.5 * (modExtend[1] + modExtend[0]), xs1 = 0.5 * (modExtend[1] - modExtend[0]);
        const double ym1 = 0.5 * (modExtend[3] + modExtend[2]), ys1 = 0.5 * (modExtend[3] - modExtend[2]);
        const std::vector<double> modExtend2 = mod2->getModuleExtent();
        const double xm2 = 0.5 * (modExtend2[1] + modExtend2[0]), xs2 = 0.5 * (modExtend2[1] - modExtend2[0]);
        const double ym2 = 0.5 * (modExtend2[3] + modExtend2[2]), ys2 = 0.5 * (modExtend2[3] - modExtend2[2]);
        const double dx = xm2 - xm1, dy = ym2 - ym1;
        return (pow(dx * dx + dy * dy, 2) / (pow((xs1 + xs2) * dx, 2) + pow((ys1 + ys2) * dy, 2)) < 2.);
}


rb_Hit::rb_Hit(const int iHit, const int mod, const int row, const TrackerHit& aHit) :
                _hitNum(iHit), _hitModule(mod), _hitRow(row), _hitX(aHit.getPosition()[0]), _hitY(aHit.getPosition()[1]), _hitZ(
                                aHit.getPosition()[2]), _hitPhiMeas(_calcPhiMeas(aHit)), _hitVarXY(_getVarXY(aHit.getCovMatrix())), _hitVarR(
                                _getVarR(aHit.getCovMatrix())), _hitVarZ(_getVarZ(aHit.getCovMatrix())), _used(false), _cosb(1.) {
}

void rb_Hit::print() const {
        std::cout << " rb_Hit " << _hitNum << ": " << _hitModule << ", " << _hitRow << ", " << _hitX << ", " << _hitY << ", " << _hitZ
                        << ", " << _hitPhiMeas << "; " << _hitVarXY << ", " << _hitVarZ << std::endl;
}


double rb_Hit::_calcPhiMeas(const TrackerHit& aHit) {
        const gear::TPCModule &module = Global::GEAR->getTPCParameters().getModule(_hitModule);
        double phi; // azimuthal direction
        if (module.getLocalPadLayout().getCoordinateType() == gear::PadRowLayout2D::POLAR) {
                Vector2D correction = module.getOffset();
                phi = atan2(aHit.getPosition()[0] - correction[0], correction[1] - aHit.getPosition()[1]);
                //std::cout << " Polar " << phi << std::endl;
        } else {
                phi = module.getAngle(); // + 0.5 * M_PI; ??
                //std::cout << " Cartesian " << phi << std::endl;
        }
        return phi;
}


double rb_Hit::_getVarXY(FloatVec covMatrix) {
        double var = cos(_hitPhiMeas) * cos(_hitPhiMeas) * covMatrix[0] + 2. * sin(_hitPhiMeas) * cos(_hitPhiMeas) * covMatrix[1]
                        + sin(_hitPhiMeas) * sin(_hitPhiMeas) * covMatrix[2];
        return (var > 0.) ? var : 1.;
}


double rb_Hit::_getVarR(FloatVec covMatrix) {
        double var = sin(_hitPhiMeas) * sin(_hitPhiMeas) * covMatrix[0] - 2. * sin(_hitPhiMeas) * cos(_hitPhiMeas) * covMatrix[1]
                        + cos(_hitPhiMeas) * cos(_hitPhiMeas) * covMatrix[2];
        return var;
}


double rb_Hit::_getVarZ(FloatVec covMatrix) {
        return (covMatrix[5] > 0.) ? covMatrix[5] : 1.;
}

int rb_Hit::getHitNum() const {
        return _hitNum;
}

int rb_Hit::getMod() const {
        return _hitModule;
}

int rb_Hit::getRow() const {
        return _hitRow;
}

double rb_Hit::getX() const {
        return _hitX;
}

double rb_Hit::getY() const {
        return _hitY;
}

double rb_Hit::getZ() const {
        return _hitZ;
}


void rb_Hit::getPos(double* position) const {
        position[0] = _hitX;
        position[1] = _hitY;
        position[2] = _hitZ;
}


double rb_Hit::getVarXY(const double der2) const {
        return _hitVarXY + der2 * _hitVarR;
}


double rb_Hit::getVarZ(const double der2) const {
        return _hitVarZ + der2 * _hitVarR;
}

double rb_Hit::getPhiMeas() const {
        return _hitPhiMeas;
}

bool rb_Hit::getUsed() const {
        return _used;
}


void rb_Hit::setUsed(const double cosb = 1.) {
        _used = true;
        _cosb = cosb;
}

double rb_Hit::getCosBeta() const {
        return _cosb;
}


double rb_Hit::getDistXY(const double x, const double y, const double ex, const double ey) const {
        return (_hitX - x) * ex + (_hitY - y) * ey;
}


double rb_Hit::getDistZ(const double z) const {
        return _hitZ - z;
}


rb_Doublet::rb_Doublet(rb_Hit* hit1, rb_Hit* hit2) {
        double pos1[3], pos2[3], dx, dy, dz, ds, dr;
        hit1->getPos(pos1);
        hit2->getPos(pos2);
        // hits
        _hit[0] = hit1;
        _hit[1] = hit2;
        // measurement direction
        _phiMeas = 0.5 * (hit1->getPhiMeas() + hit2->getPhiMeas());
        _cosPhi = cos(_phiMeas);
        _sinPhi = sin(_phiMeas);
        // position
        _xav = 0.5 * (pos1[0] + pos2[0]);
        _yav = 0.5 * (pos1[1] + pos2[1]);
        _zav = 0.5 * (pos1[2] + pos2[2]);
        // slope
        dx = pos2[0] - pos1[0];
        dy = pos2[1] - pos1[1];
        dz = pos2[2] - pos1[2];
        ds = sqrt(dx * dx + dy * dy); // arc length
        dr = dx * _sinPhi - dy * _cosPhi; // radial length
        // phi
        _phi = atan2(dy, dx);
        // tanl
        _tanl = dz / ds;
        // length
        _length = ds;
        // cos(beta) (angle between measurement and normal to flight direction in XY)
        _cosb = dr / ds;
        // derivatives
        const double cosb2 = _cosb * _cosb;
        _der2XY = (1. - cosb2) / cosb2; // tan(beta)^2
        _der2Z = _tanl * _tanl / cosb2; // tan(lambda)^2/cos(beta)^2
        // variances
        _varXY = hit1->getVarXY(_der2XY) + hit2->getVarXY(_der2XY);
        _varZ = hit1->getVarZ(_der2Z) + hit2->getVarZ(_der2Z);
}


bool rb_Doublet::match(rb_Hit* hit, const double distCut, const double chi2Cut) const {
        // triplet residuals
        const double dz = hit->getDistZ(_zav);
        // coarse (pre-)cut
        if (abs(dz) > distCut) return false;
        const double dxy = hit->getDistXY(_xav, _yav, _cosPhi, _sinPhi);
        // coarse (pre-)cut
        if (abs(dxy) > distCut) return false;
        const double chi2xy = dxy * dxy / (hit->getVarXY(_der2XY) + 0.25 * _varXY);
        const double chi2z = dz * dz / (hit->getVarZ(_der2Z) + 0.25 * _varZ);
        // cout << "chi2 " << hit->getRow() << ": " << chi2xy << " " << chi2z << endl;
        return (chi2xy < chi2Cut and chi2z < chi2Cut);
}


void rb_Doublet::getParameters(double& x, double& y, double& z, double& phiMeas, double& phi, double& dzds, double& ds,
                double& cosb, double& varXY, double& varZ) const {
        x = _xav;
        y = _yav;
        z = _zav;
        phiMeas = _phiMeas;
        phi = _phi;
        dzds = _tanl;
        ds = _length;
        cosb = _cosb;
        varXY = _varXY;
        varZ = _varZ;
}


rb_Hit* rb_Doublet::getHit(const int ihit) const {
        return _hit[ihit];
}


rb_Triplet::rb_Triplet(rb_Doublet& doublet, rb_Hit* hit) {
        //hits
        _hit[0] = doublet.getHit(0);
        _hit[1] = hit;
        _hit[2] = doublet.getHit(1);
        // parameters
        double sumVarXY, sumVarZ;
        doublet.getParameters(_posX, _posY, _posZ, _phiMeas, _phi, _tanl, _length, _cosb, sumVarXY, sumVarZ);
        // covarince matrices (diagonals)
        _varXY[0] = 0.25 * sumVarXY;
        _varXY[1] = sumVarXY / (_length * _length); // var(phi) = sumVarXY / (dr)^2 = sumVarXY / (ds * cosb)^2 = varXY[1] / cosb^2
        _varZ[0] = 0.25 * sumVarZ;
        _varZ[1] = sumVarZ / (_length * _length); // var(tanl) = sumVarZ / (ds)^2 = varZ[1]
}

int rb_Triplet::getRow() const {
        return _hit[1]->getRow();
}

double rb_Triplet::getX() const {
        return _posX;
}

double rb_Triplet::getY() const {
        return _posY;
}

double rb_Triplet::getZ() const {
        return _posZ;
}


void rb_Triplet::getPos(double* position) const {
        position[0] = _posX;
        position[1] = _posY;
        position[2] = _posZ;
}


double rb_Triplet::getVarXYPos(const double ds) const {
        return _varXY[0] + _varXY[1] * ds * ds;
}

double rb_Triplet::getVarXYDir() const {
        return _varXY[1] * (_cosb * _cosb);
}


double rb_Triplet::getVarZPos(const double ds) const {
        return _varZ[0] + _varZ[1] * ds * ds;
}

double rb_Triplet::getVarZDir() const {
        return _varZ[1];
}

double rb_Triplet::getPhiMeas() const {
        return _phiMeas;
}

double rb_Triplet::getPhi() const {
        return _phi;
}

double rb_Triplet::getTanl() const {
        return _tanl;
}

double rb_Triplet::getCosBeta() const {
        return _cosb;
}


bool rb_Triplet::match(rb_Triplet* trp, const double posCut, const double dirCut) const {
        const int rowDiff = trp->getRow() - this->getRow();
        if (rowDiff == 0) return false; // min. row distance > 0
        // compare triplet directions
        const double dphi = trp->getPhi() - _phi;
        const double dtanl = trp->getTanl() - _tanl;
        const double varXYDir = this->getVarXYDir() + trp->getVarXYDir();
        const double varZDir = this->getVarZDir() + trp->getVarZDir();
        //cout << "chi2Dir " << this->getRow() << " " << trp->getRow() << " " << dphi << " " << varXYDir << " " << dtanl << " " << varZDir << endl;
        if (dphi * dphi / varXYDir > dirCut or dtanl * dtanl / varZDir > dirCut) return false;
        // compare triplet positions at mid point
        double midPos1[3], midPos2[3], ds1, ds2;
        const double xav = 0.5 * (_posX + trp->getX());
        const double yav = 0.5 * (_posY + trp->getY());
        this->getPosCloseTo(xav, yav, midPos1, ds1);
        trp->getPosCloseTo(xav, yav, midPos2, ds2);
        const double phiMeas = 0.5 * (this->getPhiMeas() + trp->getPhiMeas());
        const double dxy = (midPos2[0] - midPos1[0]) * cos(phiMeas) + (midPos2[1] - midPos1[1]) * sin(phiMeas);
        const double dz = midPos2[2] - midPos1[2];
        const double varXYPos = this->getVarXYPos(ds1) + trp->getVarXYPos(ds2);
        const double varZPos = this->getVarZPos(ds1) + trp->getVarZPos(ds2);
        const double chi2xy = dxy * dxy / varXYPos;
        const double chi2z = dz * dz / varZPos;
        //cout << "chi2Pos " << this->getRow() << " " << trp->getRow() << " " << dxy << " " << varXYPos << " " << dz << " " << varZPos << endl;
        return (chi2xy < posCut and chi2z < posCut);
}


void rb_Triplet::getPosCloseTo(const double xref, const double yref, double* position, double& ds) const {
        const double ex = cos(_phi);
        const double ey = sin(_phi);
        ds = (xref - _posX) * ex + (yref - _posY) * ey;
        position[0] = _posX + ds * ex;
        position[1] = _posY + ds * ey;
        position[2] = _posZ + _tanl * ds;
}


rb_Hit* rb_Triplet::getHit(const int ihit) const {
        return _hit[ihit];
}


rb_Segment::rb_Segment(double Bzc, trpListType triplets) :
                _segId(0), _hitList(), _bzc(Bzc), _npar((Bzc != 0.) ? 5 : 4), _parameters(_npar), _covariance(_npar) {
        rowHitMapType rowList; // map of rows with lists of hits
        for (trpListType::iterator itTrp = triplets.begin(); itTrp != triplets.end(); ++itTrp) {
                for (int i = 0; i < 3; ++i) {
                        rb_Hit* iHit = (*itTrp)->getHit(i);
                        // Take all (previously) unused hits from the triplets
                        if (not iHit->getUsed()) {
                                rowList[iHit->getRow()].push_back(iHit);
                                iHit->setUsed((*itTrp)->getCosBeta());
                        }
                }
        }
        // use only rows with single hit
        for (rowHitMapType::iterator itRow = rowList.begin(); itRow != rowList.end(); ++itRow) {
                if (itRow->second.size() == 1) _hitList.push_back(itRow->second[0]);
        }
        fitSegment(Bzc);
}


rb_Segment::rb_Segment(double Bzc, hitListType hits, unsigned int maxHitsPerRow) :
                _segId(0), _hitList(), _bzc(Bzc), _npar((Bzc != 0.) ? 5 : 4), _parameters(_npar), _covariance(_npar) {
        rowHitMapType rowList; // map of rows with lists of hits
        for (hitListType::iterator itHit = hits.begin(); itHit != hits.end(); ++itHit) {
                rb_Hit* iHit = (*itHit);
                // Take all (previously) unused hits from the triplets
                if (not iHit->getUsed()) {
                        rowList[iHit->getRow()].push_back(iHit);
                }
        }
        // use only rows with single hit
        for (rowHitMapType::iterator itRow = rowList.begin(); itRow != rowList.end(); ++itRow) {
                if (itRow->second.size() <= maxHitsPerRow)
                        for (hitListType::iterator itHit = itRow->second.begin(); itHit != itRow->second.end(); ++itHit) {
                                rb_Hit* iHit = *itHit;
                                _hitList.push_back(iHit);
                                iHit->setUsed(1.);
                        }
        }
        fitSegment(Bzc);
}


rb_Segment::rb_Segment(segListType segments) :
                _segId(0), _hitList(), _bzc(segments.front()->getBzc()), _npar((_bzc != 0.) ? 5 : 4), _parameters(_npar), _covariance(_npar) {
        if (segments.size() > 1) {
                // combine lists of hits
                rowHitMapType rowList; // map of rows with lists of hits
                for (segListType::iterator itSeg = segments.begin(); itSeg != segments.end(); ++itSeg) {
                        hitListType hits = (*itSeg)->getHitList();
                        for (hitListType::iterator itH = hits.begin(); itH != hits.end(); ++itH)
                                rowList[(*itH)->getRow()].push_back(*itH);
                }
                // use only rows with single hit
                for (rowHitMapType::iterator itRow = rowList.begin(); itRow != rowList.end(); ++itRow) {
                        if (itRow->second.size() == 1) _hitList.push_back(itRow->second[0]);
                }
                // new reference point
                calcRefPoint();
                double oldRef[3], oldPar[5];
                // weighted average of track parameters;
                _ndf = 0;
                _chi2 = 0.;
                for (segListType::iterator itSeg = segments.begin(); itSeg != segments.end(); ++itSeg) {
                        (*itSeg)->getPar(oldPar);
                        (*itSeg)->getRefPoint(oldRef);
                        simpleHelix hlx(oldPar, oldRef, _bzc);
                        TVectorD newPar(_npar);
                        TMatrixDSym cov((*itSeg)->getCov()), newCov(_npar);
                        // track state at reference point
                        hlx.getStateAt(_refX, _refY, _refZ, cov, newPar, newCov);
                        if (itSeg == segments.begin()) {
                                // first segment
                                _parameters = newPar;
                                _covariance = newCov;
                        } else {
                                // weighted average
                                TMatrixDSym weight1(_covariance);
                                weight1.Invert(); // weight from previous segments
                                TMatrixDSym weight2(newCov);
                                weight2.Invert();  // weight from current segment
                                TMatrixDSym combCov(weight1 + weight2);
                                combCov.Invert(); // combined variance
                                TVectorD combPar(combCov * (weight1 * _parameters + weight2 * newPar)); // combined parameters
                                TVectorD resPrev(_parameters - combPar);
                                _chi2 += weight1.Similarity(resPrev); // chi2 contribution from prev. seg
                                TVectorD resCurr(newPar - combPar);
                                _chi2 += weight2.Similarity(resCurr); // chi2 contribution from curr.. seg
                                _ndf += _npar; // new degrees of freedom
                                _parameters = combPar;
                                _covariance = combCov;
                        }
                        _ndf += (*itSeg)->getNdf();
                        _chi2 += (*itSeg)->getChi2();
                }
                //cout << " combined segment " << _hitList.size() << " " << _chi2 << " " << _ndf << endl;
        } else {
                // 'copy' input segment
                rb_Segment* seg(segments.front());
                _hitList = seg->getHitList();
                double refPoint[3];
                seg->getRefPoint(refPoint);
                _refX = refPoint[0];
                _refY = refPoint[1];
                _refZ = refPoint[2];
                _parameters = seg->getPar();
                _covariance = seg->getCov();
                _ndf = seg->getNdf();
                _chi2 = seg->getChi2();
                //cout << " single segment " << _hitList.size() << " " << _chi2 << " " << _ndf << endl;
        }
}

double rb_Segment::getBzc() const {
        return _bzc;
}

int rb_Segment::getFirstRow() const {
        return _hitList.front()->getRow();
}

int rb_Segment::getLastRow() const {
        return _hitList.back()->getRow();
}


void rb_Segment::getRefPoint(double* position) const {
        position[0] = _refX;
        position[1] = _refY;
        position[2] = _refZ;
}

int rb_Segment::getNdf() const {
        return _ndf;
}

double rb_Segment::getChi2() const {
        return _chi2;
}


void rb_Segment::getPar(double* par) const {
        par[0] = 0.;
        for (int i = 0; i < _npar; ++i)
                par[5 - _npar + i] = _parameters[i];
}

const hitListType& rb_Segment::getHitList() const {
        return _hitList;
}

const TVectorD& rb_Segment::getPar() const {
        return _parameters;
}

const TMatrixDSym& rb_Segment::getCov() const {
        return _covariance;
}

int rb_Segment::getId() const {
        return _segId;
}


void rb_Segment::setId(int id) {
        _segId = id;
}


bool rb_Segment::match(rb_Segment* seg, const double chi2Cut, const int trpStepSize) const {
        int f1 = this->getFirstRow();
        int l1 = this->getLastRow();
        int f2 = seg->getFirstRow();
        int l2 = seg->getLastRow();
        // segments must not overlap (within triplet granularity)
        if (max(f1, f2) <= min(l1, l2) - trpStepSize) return false;
        // (row) gap must smaller than average size
        if (min(abs(l1 - f2), abs(l2 - f1)) > (l1 + l2 - f1 - f2) / 2) return false;
        // helices describing segments
        double refPoint1[3], refPoint2[3], par1[5], par2[5];
        this->getRefPoint(refPoint1);
        this->getPar(par1);
        simpleHelix hlx1(par1, refPoint1, _bzc);
        seg->getRefPoint(refPoint2);
        seg->getPar(par2);
        simpleHelix hlx2(par2, refPoint2, _bzc);
        // mid point
        const double xav = 0.5 * (refPoint1[0] + refPoint2[0]);
        const double yav = 0.5 * (refPoint1[1] + refPoint2[1]);
        const double zav = 0.5 * (refPoint1[2] + refPoint2[2]);
        TVectorD newPar1(_npar), newPar2(_npar);
        TMatrixDSym cov1(this->getCov()), cov2(seg->getCov()), newCov1(_npar), newCov2(_npar);
        // track states at midpoint
        hlx1.getStateAt(xav, yav, zav, cov1, newPar1, newCov1);
        hlx2.getStateAt(xav, yav, zav, cov2, newPar2, newCov2);
        TVectorD diffPar(newPar1 - newPar2);
        TMatrixDSym diffPrec(newCov1 + newCov2);
        diffPrec.Invert();
        // chi2
        double chi2 = diffPrec.Similarity(diffPar);
        //cout << " seg-match chi2 " << chi2 << endl;
        return (chi2 < chi2Cut * _npar);
}


bool rb_Segment::match(rb_Hit* hit, const double chi2Cut) const {
        double refPoint[3], par[5], pos[3], phiMeas, xyPos, zPos, sArc, phi, lambda;
        this->getRefPoint(refPoint);
        this->getPar(par);
        simpleHelix hlx(par, refPoint, _bzc);
        hit->getPos(pos);
        phiMeas = hit->getPhiMeas();
        // get prediction from helix
        if (!hlx.getExpectedPlanePos(pos, phiMeas, xyPos, zPos, sArc, phi, lambda)) return false;
        // cos(beta) (angle between measurement and normal to flight direction in XY)
        const double cosb2 = pow(sin(phiMeas - phi), 2);
        const double der2XY = (1. - cosb2) / cosb2;
        const double der2Z = pow(tan(lambda), 2) / cosb2;
        TVectorD newPar(_npar);
        TMatrixDSym cov(this->getCov()), newCov(_npar);
        // track state at predicted point
        pos[0] += xyPos * cos(phiMeas);
        pos[1] += xyPos * sin(phiMeas);
        pos[2] += zPos;
        hlx.getStateAt(pos[0], pos[1], pos[2], cov, newPar, newCov);
        double chi2xy = xyPos * xyPos / (newCov[_npar - 3][_npar - 3] / (cosb2) + hit->getVarXY(der2XY));
        double chi2z = zPos * zPos / (newCov[_npar - 1][_npar - 1] + hit->getVarZ(der2Z));
        //cout << " match hit, chi2 " << hit->getRow() << ": " << chi2xy << " " << chi2z << "/ " << zPos << endl;
        return (chi2xy < chi2Cut and chi2z < chi2Cut);
}


void rb_Segment::match(hitListType& hits, hitListType& matchingHits, const double chi2Cut) const {
        for (hitListType::iterator itHit = hits.begin(); itHit != hits.end(); ++itHit) {
                if ((*itHit)->getUsed()) continue; // hit already used
                if (this->match(*itHit, chi2Cut)) matchingHits.push_back(*itHit);
        }
}


void rb_Segment::getLCIOStateAtRefPoint(const double* refPos, TVectorD& newPar, TMatrixDSym& newCov) const {
        double refPoint[3], par[5];
        this->getRefPoint(refPoint);
        this->getPar(par);
        simpleHelix hlx(par, refPoint, _bzc);
        TVectorD tmpPar(_npar);
        TMatrixDSym cov(this->getCov()), tmpCov(_npar);
        // track state at reference point
        hlx.getStateAt(refPos[0], refPos[1], refPos[2], cov, tmpPar, tmpCov);
        TMatrixD hlx2LCIO = hlx.helixToLCIOJacobian().GetSub(0, 4, 5 - _npar, 4);
        newPar = hlx2LCIO * tmpPar;
        newCov = tmpCov.Similarity(hlx2LCIO);
}

void rb_Segment::calcRefPoint() {
        double pos1[3], pos2[3];
        _hitList.front()->getPos(pos1);
        _hitList.back()->getPos(pos2);
        _refX = 0.5 * (pos1[0] + pos2[0]);
        _refY = 0.5 * (pos1[1] + pos2[1]);
        _refZ = 0.5 * (pos1[2] + pos2[2]);
}

void rb_Segment::fitSegment(double Bzc) {
        unsigned int nHit = _hitList.size();
        if (nHit < 4) {
                // to few hits
                _ndf = -1;
                _chi2 = -1.0;
                return;
        }
        //cout << " segment " << _hitList.size() << endl;
        calcRefPoint();
        double refPoint[] = { _refX, _refY, _refZ };
        //
        // fit segment
        //
        double pos[3], sArc, chi2XY, chi2ZS, cosb2, der2XY, derZ, phi, dzds;
        double parameters[] = { 0., 0., 0., 0., 0. };
        int nParXY, nPointXY, nPointZS, nParZS;
        // coarse unweighted xy prefit
        simpleFitXY xyPreFit((Bzc != 0.), _refX, _refY);
        for (hitListType::iterator itH = _hitList.begin(); itH != _hitList.end(); ++itH) {
                xyPreFit.addPoint((*itH)->getX(), (*itH)->getY(), 1.0);
        }
        nParXY = xyPreFit.fit(chi2XY, nPointXY);
        _parameters.SetSub(0, xyPreFit.getPar());
        for (int i = 0; i < nParXY; ++i)
                parameters[3 - nParXY + i] = _parameters[i];
        // simple helix
        simpleHelix hlx0(parameters, refPoint, Bzc);
        // xy fit (with weights)
        simpleFitXY xyFit((Bzc != 0.), _refX, _refY);
        for (hitListType::iterator itH = _hitList.begin(); itH != _hitList.end(); ++itH) {
                (*itH)->getPos(pos);
                sArc = hlx0.getArcLengthXY(pos);
                phi = parameters[0] * sArc + parameters[1];
                cosb2 = pow(sin(phi - (*itH)->getPhiMeas()), 2);
                der2XY = (1. - cosb2) / cosb2;
                xyFit.addPoint(pos[0], pos[1], 1.0 / ((*itH)->getVarXY(der2XY) * cosb2));
        }
        nParXY = xyFit.fit(chi2XY, nPointXY);
        _parameters.SetSub(0, xyFit.getPar());
        _covariance.SetSub(0, 0, xyFit.getCov());
        for (int i = 0; i < nParXY; ++i)
                parameters[3 - nParXY + i] = _parameters[i];
        // simple helix
        simpleHelix hlx(parameters, refPoint, Bzc);
        // coarse unweighted zs prefit
        simpleFitZS zsPreFit;
        for (hitListType::iterator itH = _hitList.begin(); itH != _hitList.end(); ++itH) {
                (*itH)->getPos(pos);
                sArc = hlx.getArcLengthXY(pos);
                zsPreFit.addPoint(sArc, pos[2] - _refZ, 1.0);
        }
        nParZS = zsPreFit.fit(chi2ZS, nPointZS);
        _parameters.SetSub(nParXY, zsPreFit.getPar());
        dzds = _parameters[nParXY];
        // zs fit (with weights)
        simpleFitZS zsFit;
        for (hitListType::iterator itH = _hitList.begin(); itH != _hitList.end(); ++itH) {
                (*itH)->getPos(pos);
                sArc = hlx.getArcLengthXY(pos);
                phi = parameters[0] * sArc + parameters[1];
                cosb2 = pow(sin(phi - (*itH)->getPhiMeas()), 2);
                derZ = dzds * dzds / cosb2;
                zsFit.addPoint(sArc, pos[2] - _refZ, 1.0 / (*itH)->getVarZ(derZ));
        }
        nParZS = zsFit.fit(chi2ZS, nPointZS);
        _parameters.SetSub(nParXY, zsFit.getPar());
        _covariance.SetSub(nParXY, nParXY, zsFit.getCov());
        // combine
        _ndf = nPointXY + nPointZS - _npar;
        _chi2 = chi2XY + chi2ZS;
        // cout << " cov " << endl; _covariance.Print();
}


void rb_Segment::fillRowMap(std::map<int, int>& rowMap) const {
        for (hitListType::const_iterator itH = _hitList.begin(); itH != _hitList.end(); ++itH)
                rowMap[(*itH)->getRow()] = 1;
}


void rb_Segment::addHit(rb_Hit* hit) {
        const int row = hit->getRow();
        if (row > _hitList.back()->getRow())
                _hitList.push_back(hit); // add to end
        else {
                for (hitListType::iterator itH = _hitList.begin(); itH != _hitList.end(); ++itH)
                        if ((*itH)->getRow() > row) {
                                _hitList.insert(itH, hit);
                                break;
                        } // insert
        }
        hit->setUsed();
}

simpleEquiClasses::simpleEquiClasses() :
                _index(), _matches() {
}


void simpleEquiClasses::addIndex(int aIndex) {
        _index[aIndex] = aIndex;
}


void simpleEquiClasses::addMatch(std::pair<int, int> aMatch) {
        _matches.push_back(aMatch);
}


std::map<int, std::vector<int> > simpleEquiClasses::getClasses() {
        for (std::vector<std::pair<int, int> >::iterator it = _matches.begin(); it != _matches.end(); ++it) {
                int i1 = it->first; // trace first element up to is ancestor
                while (_index[i1] != i1)
                        i1 = _index[i1];
                int i2 = it->second; // trace second element up to is ancestor
                while (_index[i2] != i2)
                        i2 = _index[i2];
                if (i1 != i2) _index[i1] = i2;
        }
        // final sweep-up to highest ancestor
        for (std::map<int, int>::iterator it = _index.begin(); it != _index.end(); ++it) {
                int i1 = it->first;
                while (_index[i1] != _index[_index[i1]])
                        _index[i1] = _index[_index[i1]];
        }
        std::map<int, std::vector<int> > classes;
        for (std::map<int, int>::iterator it = _index.begin(); it != _index.end(); ++it) {
                classes[it->second].push_back(it->first);
        }
        return classes;
}

} // close namespace brackets

// #endif //#ifdef USE_ROWTRIPLETS