RowBasedPadPulseRoadSearchProcessor.cc

Go to the documentation of this file.

// #ifdef USE_ROWPPRS
#include "RowBasedPadPulseRoadSearchProcessor.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>

// the common flag word definitions in MarlinTPC
#include "FlagwordDefinitions.h"

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

namespace marlintpc {

static bool comparePadRow(pb_Pulse* one, pb_Pulse* two) {
        return ((one->getRow()) < (two->getRow())); // increasing
}

static bool comparePadCol(pb_Pulse* one, pb_Pulse* two) {
        return ((one->getColumn()) < (two->getColumn())); // increasing
}

static bool compareRelCharge(pb_Pulse* one, pb_Pulse* two) {
        return ((one->getRelCharge()) > (two->getRelCharge())); // decreasing
}

static bool compareCharge(pb_Pulse* one, pb_Pulse* two) {
        return ((one->getCharge()) > (two->getCharge())); // decreasing
}

RowBasedPadPulseRoadSearchProcessor aRowBasedPadPulseRoadSearchProcessor;

RowBasedPadPulseRoadSearchProcessor::RowBasedPadPulseRoadSearchProcessor() :
                Processor("RowBasedPadPulseRoadSearchProcessor") {
        // the processor description
        _description = "RowBasedPadPulseRoadSearchProcessor is intended to simultaneously find TPC hits and tracks";

        registerInputCollection(LCIO::TRACKERPULSE, "InputTrackerPulses", "Name of the input TrackerPulses collection",
                        _inputTrackerPulsesColName, string("TPCPulses"));

        registerOutputCollection(LCIO::TRACKERHIT, "OutputTrackerHits", "Name of the output TrackerHits collection",
                        _outputTrackerHitsColName, string("TPCHits"));

        registerOutputCollection(LCIO::TRACK, "OutputTracks", "Name of the output Track Collection", _outputTrackColName,
                        string("TPCTracks"));

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

        registerOptionalParameter("DriftVelocity", "Drift velocity", _driftVelocity, double(74.0));

        registerOptionalParameter("MaxColDiffNeighbourhood", "Maximal column difference to pulse in neighbourhood", _maxColDiff, int(1));

        registerOptionalParameter("MaxTimeDiff", "Maximum drift time difference in neighboorhood/hit", _dtMax, double(80.0));

        registerOptionalParameter("MaxSeedsInRow", "Maximum number of seeding pulses in row", _maxSeedsInRow, int(10));

        registerOptionalParameter("MinRowDiff", "Minimum row difference for pulse pair seeding road", _minRowDiff, int(5));

        registerOptionalParameter("MaxDxyRoad", "Maximum interpolation deviation in XY (road search)", _maxDxyRoad, double(2.0));

        registerOptionalParameter("MaxDzRoad", "Maximum interpolation deviation in Z (road search)", _maxDzRoad, double(3.0));

        registerOptionalParameter("MinRowNum", "Minimum number of rows for segment seed", _minRows, int(6));

        registerOptionalParameter("MinRowDensity", "Minimum row 'density' for segment seed", _rowDensityCut, double(0.8));

        registerOptionalParameter("MaxVarXY", "Maximum (average) XY variance (unweighted Chi2/ndf) for segment seed", _maxVarXY,
                        double(0.25));

        registerOptionalParameter("MaxOverlapFraction", "Maximum overlap fraction of segment seeds", _maxOverlapFraction, double(0.5));

        registerOptionalParameter("MaxColDiffHit", "Maximum distance of pad to seed in hit (if > 0)", _maxColDiffHit, int(0));

        registerOptionalParameter("MinNumPulses", "Minimum number of pulses in hit", _minNumPulses, int(2));

        registerOptionalParameter("MaxColGap", "Maximum column gap in hit", _maxColGap, int(1));

        registerProcessorParameter("OffsetCol", "sigma0^2 in column direction", _offsetCol, double(0.0085));

        registerProcessorParameter("SlopeCol", "diffusion in column direction", _slopeCol, double(0.00011));

        registerProcessorParameter("OffsetDrift", "sigma0^2 in drift direction", _offsetDrift, double(0.25));

        registerProcessorParameter("SlopeDrift", "diffusion term in drift direction", _slopeDrift, double(0.));

        registerProcessorParameter("OffsetRow", "sigma0^2 in row direction", _offsetRow, double(0.));

        registerProcessorParameter("SlopeRow", "diffusion in row direction", _slopeRow, double(0.));

        registerOptionalParameter("ChargeConversionFactor", "Optional: The charge conversion factor from ADC values to primary electrons",
                        _chargeConversionFactor, double(1.));

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

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

        registerOptionalParameter("SkipRoadSearch", "Skip road search, hit finding only", _skipRoadSearch, bool(false));

        registerOptionalParameter("SkipMultiplePulses", "Skip multiple pulse candidates as road search seeds", _skipMultiplePulses,
                        bool(true));
}

void RowBasedPadPulseRoadSearchProcessor::init() {
        streamlog_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();

        // get row offsets per module layer from module offset[0]
        std::vector<double> xOffset;
        for (std::vector<TPCModule*>::const_iterator itMod = moduleVector.begin(); itMod != moduleVector.end(); ++itMod) {
                Vector2D offset = (*itMod)->getOffset();
                xOffset.push_back(offset[0]);
        }
        sort(xOffset.begin(), xOffset.end());
        // get list of (different) offsets
        double eps = 0.1; // precision of offsets
        std::vector<double> layerOffset;
        layerOffset.push_back(xOffset.front());
        for (std::vector<double>::iterator it = xOffset.begin() + 1; it != xOffset.end(); ++it)
                if (abs(layerOffset.back() - (*it)) > eps) layerOffset.push_back(*it);
        const int nLayer = layerOffset.size();
        std::vector<int> layerRows(nLayer);
        // get maximum number of rows per layer
        for (std::vector<TPCModule*>::const_iterator itMod = moduleVector.begin(); itMod != moduleVector.end(); ++itMod) {
                Vector2D offset = (*itMod)->getOffset();
                int layer = 0;
                while (abs(layerOffset[layer] - offset[0]) > eps)
                        layer++;
                layerRows[layer] = max(layerRows[layer], (*itMod)->getNRows());
        }

        streamlog_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();
                Vector2D offset = (*itMod)->getOffset();
                _moduleRowOffset[module] = 0;
                int layer = 0;
                while (abs(layerOffset[layer] - offset[0]) > eps) {
                        _moduleRowOffset[module] += layerRows[layer];
                        layer++;
                }
                streamlog_out(DEBUG) << " module " << module << "  with " << (*itMod)->getNRows() << " rows, offset "
                                << _moduleRowOffset[module] << endl;
                const std::vector<double> modExtend = (*itMod)->getLocalModuleExtent();

                Vector2D localOffset = (*itMod)->globalToLocal(offset[0], offset[1]);
                for (std::vector<TPCModule*>::const_iterator itMod2 = itMod; itMod2 != moduleVector.end(); ++itMod2) {
                        if (areNeighbourModules(*itMod, *itMod2, cartesian)) {
                                _modNeighbours[module].push_back((*itMod2)->getModuleID());
                                streamlog_out(DEBUG) << "   neighbour " << (*itMod2)->getModuleID() << endl;
                        }
                }
        }

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

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

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

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

        //get input collection
        LCCollection* inputCol = 0;

        LCCollectionVec* outputHitCol = new LCCollectionVec(LCIO::TRACKERHIT);
        //tpcHitCollection->setTransient(!_outputHitsPersistent);
        // set the flags that cellID1 should be stored
        LCFlagImpl hitFlag(0);
        hitFlag.setBit(LCIO::RTHBIT_ID1);
        outputHitCol->setFlag(hitFlag.getFlag());

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

        // do a try / catch in case there is an empty event
        try {
                inputCol = evt->getCollection(_inputTrackerPulsesColName);
        } 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 gear::TPCParameters& tpcParameters = Global::GEAR->getTPCParameters();

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

        UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);

        modPulseMapType pulseData;

        int nInput = inputCol->getNumberOfElements();
        streamlog_out(DEBUG) << "Event " << evt->getEventNumber() << " in run " << evt->getRunNumber() << ", number of input pulses: "
                        << nInput << endl;

        double maxDriftLength = tpcParameters.getMaxDriftLength();
        Vector2D globalPadCenter;
        for (int iInput = 0; iInput < nInput; iInput++) {
                TrackerPulse* aPulse = dynamic_cast<TrackerPulse*>(inputCol->getElementAt(iInput));
                const int moduleID = aPulse->getCellID1();
                const TPCModule& theModule = tpcParameters.getModule(moduleID);

                //  coordinate system: z=0 at cathode, for large prototype the anode is in negative z-direction
                //  z-position is therefore: zAnode (max drift length, positive) - driftTime * driftVelocity
                double zPosition = maxDriftLength - _driftVelocity * (aPulse->getTime()) / 1000.; // time is in [ns], drift velocity in [mm/mus]

                pb_Pulse* somePulse = new pb_Pulse(iInput, *aPulse, theModule, zPosition);
                //somePulse->print();
                pulseData[moduleID][somePulse->getRow()].push_back(somePulse);
        }

        hitListType hitData;
        int numHits = 0;
        // map of modules with vectors of segments
        modSegMapType segData;
        for (modPulseMapType::iterator itMod = pulseData.begin(); itMod != pulseData.end(); ++itMod) {
                const int moduleID = itMod->first;
                streamlog_out(DEBUG) << " === Module " << moduleID << std::endl;
                const TPCModule& theModule = tpcParameters.getModule(moduleID);
                rowPulseMapType modData = itMod->second; // rows with pulses in module
                pulseListType leadingPulses, trailingPulses, seedingPulses; // pulses in module
                for (rowPulseMapType::iterator itRow = modData.begin(); itRow != modData.end(); ++itRow) {
                        //std::cout << "  Row " << itRow->first << std::endl;
                        // sort row by column
                        sort(itRow->second.begin(), itRow->second.end(), comparePadCol);
                        for (pulseListType::iterator itP1 = itRow->second.begin(); itP1 != itRow->second.end() - 1; ++itP1) {
                                int col1 = (*itP1)->getColumn();
                                int colLast = col1; // last (neighbour) column
                                for (pulseListType::iterator itP2 = itP1 + 1; itP2 != itRow->second.end(); ++itP2) {
                                        int col2 = (*itP2)->getColumn();
                                        // out of neighbourhood?
                                        if (col2 > col1 + _maxColDiff or col2 > colLast + 1) break;
                                        // matching in time?
                                        if (abs((*itP1)->getTime() - (*itP2)->getTime()) < _dtMax) {
                                                colLast = col2;
                                                (*itP1)->addNeighbour(*itP2);
                                                (*itP2)->addNeighbour(*itP1);
                                        }
                                }
                        }
                        // split into leading and trailing pulses, set index
                        const int numSeeding = seedingPulses.size();
                        for (pulseListType::iterator itP = itRow->second.begin(); itP != itRow->second.end(); ++itP) {
                                (*itP)->setIndex(itP - itRow->second.begin());
                                if ((*itP)->getRelCharge() > 1.0) {
                                        leadingPulses.push_back(*itP);
                                        if (pulseflag::isMultiplePulseCandidate((*itP)->getQuality()) and _skipMultiplePulses) continue; // skip for road seeding
                                        if ((*itP)->isMaxPulse()) seedingPulses.push_back(*itP); // take largest pulse from neighbourhood
                                } else
                                        trailingPulses.push_back(*itP);
                        }
                        // too many seeding pulses in row?
                        if (seedingPulses.size() > static_cast<unsigned int>(numSeeding + _maxSeedsInRow)) seedingPulses.resize(numSeeding); // skip row
                }
                if (_skipRoadSearch) seedingPulses.clear(); // only hit finding
                streamlog_out(DEBUG) << "  pulses (leading, seeding, trailing) " << leadingPulses.size() << " " << seedingPulses.size() << " "
                                << trailingPulses.size() << " " << std::endl;
                // sort row by relative charge
                sort(seedingPulses.begin(), seedingPulses.end(), compareRelCharge);
                seedListType seedList;
                pulseListType insidePulses; // (leading or trailing) pulses inside road
                insidePulses.reserve(max(leadingPulses.size(), trailingPulses.size()));
                // number of segment seeds
                int numSegSeeds = 0;
                // loop over pairs (of (unused) seeding pulses)
                double pos1[3], pos2[3];
                for (pulseListType::iterator itP1 = seedingPulses.begin(); itP1 != seedingPulses.end(); ++itP1) {
                        if ((*itP1)->inSeed()) continue;
                        const int row1 = (*itP1)->getRow();
                        (*itP1)->getPos(pos1);
                        //std::cout << " -seed pulse 1 " << row1 << " " << (*itP1)->getColumn() << " " << (*itP1)->getRelCharge() << std::endl;
                        for (pulseListType::iterator itP2 = itP1 + 1; itP2 != seedingPulses.end(); ++itP2) {
                                if ((*itP2)->inSeed()) continue;
                                const int row2 = (*itP2)->getRow();
                                if (abs(row1 - row2) >= _minRowDiff) {
                                        // (straight line) road candidate
                                        //std::cout << "  seed pulse 2 " << row2 << " " << (*itP2)->getColumn() << " " << (*itP2)->getRelCharge() << std::endl;
                                        (*itP2)->getPos(pos2);
                                        const double dx(pos2[0] - pos1[0]), dy(pos2[1] - pos1[1]), dz(pos2[2] - pos1[2]), dist = sqrt(dx * dx + dy * dy);
                                        const double ex(dx / dist), ey(dy / dist), dzds(dz / dist);
                                        //std::cout << "      par " << dist << " " << ex << " " << ey << " " << dzds << std::endl;
                                        insidePulses.clear();
                                        int numNearPulses = 0;
                                        for (pulseListType::iterator itP = leadingPulses.begin(); itP != leadingPulses.end(); ++itP) {
                                                (*itP)->getPos(pos2);
                                                const double dx(pos2[0] - pos1[0]), dy(pos2[1] - pos1[1]), dxy(dy * ex - dx * ey);
                                                const double ds(dx * ex + dy * ey), dz(pos2[2] - pos1[2] - ds * dzds);
                                                //std::cout << "      dist " << (*itP)->getRow() << " " << ds << " " << dxy << " " << dz << std::endl;
                                                if (abs(dxy) < _maxDxyRoad and abs(dz) < _maxDzRoad)
                                                        insidePulses.push_back(*itP);
                                                else if (abs(dxy) < 2. * _maxDxyRoad and abs(dz) < 2. * _maxDzRoad) numNearPulses++;
                                        }
                                        //std::cout << "      inside " << insidePulses.size() << " near " << numNearPulses << std::endl;
                                        // too few pulses found?
                                        if (insidePulses.size() < static_cast<unsigned int>(_minRows)) continue;
                                        // create seed from pulses
                                        pb_Seed* someSeed = new pb_Seed(numSegSeeds, _Bzc, insidePulses, _minRows, _rowDensityCut, _maxVarXY);
                                        if (!someSeed->getValid()) {
                                                delete someSeed;
                                                continue;
                                        }
                                        // iterate seed?
                                        if (numNearPulses > 0) {
                                                someSeed->roadSearch(leadingPulses, insidePulses, _maxDxyRoad, _maxDzRoad);
                                                delete someSeed;
                                                pb_Seed* iteratedSeed = new pb_Seed(numSegSeeds, _Bzc, insidePulses, _minRows, _rowDensityCut, _maxVarXY);
                                                if (!iteratedSeed->getValid()) {
                                                        delete iteratedSeed;
                                                        continue;
                                                }
                                                // use iterated seed
                                                someSeed = iteratedSeed;
                                        }
                                        // valid seed
                                        someSeed->markPulses();
                                        //std::cout << " --> accepted " << numSegSeeds << std::endl;
                                        numSegSeeds++;
                                        seedList.push_back(someSeed);
                                        break;
                                }
                        }
                }

                while (true) {
                        seedListType::iterator itWorst = seedList.end();
                        double worstFraction = _maxOverlapFraction;
                        for (seedListType::iterator itS = seedList.begin(); itS != seedList.end(); ++itS) {
                                if (!(*itS)->getValid()) continue;
                                int numUnique = (*itS)->getNumUniqueSeedPulses();
                                int numAll = (*itS)->getNumSeedPulses();
                                if (numUnique < 4) numUnique = -1; // require minimum number of unique pulses
                                double fraction = double(numUnique) / double(numAll);
                                if (fraction < worstFraction) {
                                        worstFraction = fraction;
                                        itWorst = itS;
                                }
                        }
                        if (itWorst == seedList.end()) break;
                        // invalidate worst      
                        //std::cout << " remove overlap seed " << itWorst - seedList.begin() << " " << worstFraction << std::endl;
                        (*itWorst)->freePulses();
                        (*itWorst)->setInvalid();
                }

                // remove ambigous pulses from seeds
                for (seedListType::iterator itS = seedList.begin(); itS != seedList.end(); ++itS) {
                        if (!(*itS)->getValid()) continue;
                        int numUnique = (*itS)->getNumUniqueSeedPulses();
                        int numAll = (*itS)->getNumSeedPulses();
                        if (numUnique < numAll) {
                                const int seedId = (*itS)->getSeedId();
                                //std::cout << " remove ambigous hits " << seedId << std::endl;
                                pulseListType allPulses = (*itS)->getPulseList();
                                pulseListType uniquePulses;
                                for (pulseListType::iterator itP = allPulses.begin(); itP != allPulses.end(); ++itP) {
                                        if ((*itP)->getNumSeeds() < 2) uniquePulses.push_back(*itP);
                                }
                                pb_Seed* uniqueSeed = new pb_Seed(seedId, _Bzc, uniquePulses);
                                //uniqueSeed->print();
                                // replace seed
                                delete *itS;
                                *itS = uniqueSeed;
                        }
                }
                // (re)attach ambigous pulses to nearest seed
                double pos[3], refPoint[3], par[5], xyPos, zPos, sArc, phi, lambda;
                int numAmbigousPulses = 0;
                int numReattachedPulses = 0;
                for (pulseListType::iterator itP = leadingPulses.begin(); itP != leadingPulses.end(); ++itP) {
                        if ((*itP)->getNumSeeds() < 2) continue;
                        // resolve by distance
                        numAmbigousPulses += 1;
                        (*itP)->getCOG(pos);
                        seedIdListType seedIds = (*itP)->getSeedIdList();
                        int bestId = -1;
                        double bestDist = _maxDxyRoad;
                        for (seedIdListType::iterator itId = seedIds.begin(); itId != seedIds.end(); ++itId) {
                                pb_Seed* aSeed = seedList[*itId];
                                // construct helix
                                aSeed->getRefPoint(refPoint);
                                aSeed->getPar(par);
                                simpleHelix hlx(par, refPoint, _Bzc);
                                if (hlx.getExpectedPlanePos(pos, (*itP)->getPhiMeas(), xyPos, zPos, sArc, phi, lambda)) {
                                        //std::cout << " expected pos " << xyPos << " " << zPos << std::endl;
                                        if (abs(xyPos) < bestDist) {
                                                bestId = *itId;
                                                bestDist = abs(xyPos);
                                        }
                                }
                        }
                        // clear seeds
                        (*itP)->clearSeeds();
                        // reattach
                        //std::cout << " reattach " << bestId << " " << bestDist << std::endl;
                        if (bestId >= 0) {
                                numReattachedPulses += 1;
                                (*itP)->addSeed(bestId);
                                seedList[bestId]->addPulse(*itP);
                        }
                }
                // update seeds with additional pulses
                for (seedListType::iterator itS = seedList.begin(); itS != seedList.end(); ++itS) {
                        if (!(*itS)->getValid()) continue;
                        if ((*itS)->hasAddedPulses()) {
                                pb_Seed* uniqueSeed = new pb_Seed((*itS)->getSeedId(), _Bzc, (*itS)->getPulseList());
                                //uniqueSeed->print();
                                // replace seed
                                delete *itS;
                                *itS = uniqueSeed;
                        }
                }

                for (seedListType::iterator itS = seedList.begin(); itS != seedList.end(); ++itS) {
                        if (!(*itS)->getValid()) continue;
                        int seedId = (*itS)->getSeedId();
                        (*itS)->roadSearch(trailingPulses, insidePulses, _maxDxyRoad, _maxDzRoad * 2.);
                        // add associations to pulses
                        for (pulseListType::iterator itP = insidePulses.begin(); itP != insidePulses.end(); ++itP)
                                (*itP)->addSeed(seedId);
                }

                for (pulseListType::iterator itP = trailingPulses.begin(); itP != trailingPulses.end(); ++itP) {
                        if ((*itP)->getNumSeeds() < 2) continue;
                        // resolve by distance
                        numAmbigousPulses += 1;
                        (*itP)->getCOG(pos);
                        seedIdListType seedIds = (*itP)->getSeedIdList();
                        int bestId = -1;
                        double bestDist = _maxDxyRoad;
                        for (seedIdListType::iterator itId = seedIds.begin(); itId != seedIds.end(); ++itId) {
                                pb_Seed* aSeed = seedList[*itId];
                                // construct helix
                                aSeed->getRefPoint(refPoint);
                                aSeed->getPar(par);
                                simpleHelix hlx(par, refPoint, _Bzc);
                                if (hlx.getExpectedPlanePos(pos, (*itP)->getPhiMeas(), xyPos, zPos, sArc, phi, lambda)) {
                                        //std::cout << " expected pos " << xyPos << " " << zPos << std::endl;
                                        if (abs(xyPos) < bestDist) {
                                                bestId = *itId;
                                                bestDist = abs(xyPos);
                                        }
                                }
                        }
                        // clear seeds
                        (*itP)->clearSeeds();
                        // reattach
                        //std::cout << " reattach " << bestId << " " << bestDist << std::endl;
                        if (bestId >= 0) {
                                numReattachedPulses += 1;
                                (*itP)->addSeed(bestId);
                        }
                }

                // similar to the RowBasedHitFinderProcessor
                // loop ober rows
                for (rowPulseMapType::iterator itRow = modData.begin(); itRow != modData.end(); ++itRow) {
                        const int row = itRow->first;
                        const pulseListType rowPulses = itRow->second;
                        const int numPulses = rowPulses.size();
                        const int firstCol = theModule.getPadNumber(theModule.getPadsInRow(row).front());
                        const int lastCol = theModule.getPadNumber(theModule.getPadsInRow(row).back());
                        const int maxColDiffHit = _maxColDiffHit > 0 ? _maxColDiffHit : lastCol - firstCol; // cut on column difference ?
                        //std::cout << " Make hits in row " << row << " " << numPulses << ", " << firstCol << ".." << lastCol << std::endl;
                        // use (only) leading pulses to seed hits
                        pulseListType seedingPulses;
                        for (pulseListType::const_iterator itP = rowPulses.begin(); itP != rowPulses.end(); ++itP)
                                if ((*itP)->getRelCharge() > 1.0) seedingPulses.push_back(*itP);
                        pulseListType hitPulses;
                        hitPulses.reserve(numPulses);
                        // sort by charge
                        sort(seedingPulses.begin(), seedingPulses.end(), compareCharge);
                        for (pulseListType::iterator itP = seedingPulses.begin(); itP != seedingPulses.end(); ++itP) {
                                //(*itP)->print();
                                if ((*itP)->inHit()) continue;
                                // start hit with this pulse
                                const double tSeed = (*itP)->getTime();
                                const int indSeed = (*itP)->getIndex();
                                const int colSeed = (*itP)->getColumn();
                                hitPulses.clear();
                                hitPulses.push_back(*itP);
                                int qualityFlag = 0;
                                // segment seed ID of hit seed
                                int segSeedId = (*itP)->getSeedId();
                                //std::cout << " hitSeed " << colSeed << " " << indSeed << " " << segSeedId << ", " << (*itP)->getPadId() << std::endl;
                                // check previous/next pulses/columns
                                for (int iDir = -1; iDir <= 1; iDir += 2) {
                                        int ind = indSeed + iDir;
                                        int colLast = colSeed;
                                        double dtLast = 0.;
                                        while (ind >= 0 and ind < numPulses) {
                                                pb_Pulse* somePulse = rowPulses[ind];
                                                ind += iDir;
                                                if (somePulse->inHit()) continue; // skip used pulses
                                                const int col = somePulse->getColumn();
                                                const int colDiff = abs(col - colLast);
                                                const int segId = somePulse->getSeedId();
                                                const int maxGap = (segId >= 0) and (segId == segSeedId) ? _maxColGap : 0; // max. allowed column gap
                                                if (colDiff > maxGap + 1) break; // empty (or dead?) pad (after last accepted pad)
                                                if (abs(col - colSeed) > maxColDiffHit) break; // distance to seed too large
                                                double dTime = abs(somePulse->getTime() - tSeed);
                                                if (dTime > _dtMax) continue; // skip pulses not matching in time
                                                if ((segId >= 0) and (segId != segSeedId)) {
                                                        hitflag::setMultipleHitCandidate(qualityFlag);
                                                        continue; // skip pulses from different (segment) seed
                                                }
                                                // same column?
                                                if (colDiff == 0) {
                                                        hitflag::setMultipleHitCandidate(qualityFlag);
                                                        // take smaller time difference
                                                        if (dTime >= dtLast) continue; // skip new pulse
                                                        hitPulses.pop_back(); // remove last accepted pulse
                                                } else {
                                                        // new column
                                                        if (colDiff > 1) hitflag::setDeadChannel(qualityFlag); // missing pad detected
                                                        // check charge for rerise, may indicate multiple hit
                                                        if (segId < 0)
                                                                if (somePulse->getCharge() > hitPulses.back()->getCharge()) hitflag::setMultipleHitCandidate(qualityFlag);
                                                }
                                                // accept pulse
                                                hitPulses.push_back(somePulse);
                                                colLast = col;
                                                dtLast = dTime;
                                        }
                                        if (iDir < 0) {
                                                if (colLast == firstCol) hitflag::setAtModuleBorder(qualityFlag); // at module border
                                                // correct order
                                                std::reverse(hitPulses.begin(), hitPulses.end());
                                        } else if (colLast == lastCol) hitflag::setAtModuleBorder(qualityFlag); // at module border
                                }
                                //std::cout << " hit: " << colSeed << " " << hitPulses.front()->getColumn() << " " << hitPulses.back()->getColumn() << " " << hitPulses.size() << std::endl;
                                // to few pulses? 
                                if (hitPulses.size() < static_cast<unsigned int>(_minNumPulses)) continue;
                                // hit accepted
                                (*itP)->setHitNum(numHits);
                                // calculate hit position
                                Vector2D localPadCenter;
                                double sumCharge = 0.;
                                double varCharge = 0.;
                                double sumChargePhiMeas = 0.;
                                double localPosition[2] = { 0., 0. };
                                for (pulseListType::iterator itHitP = hitPulses.begin(); itHitP != hitPulses.end(); ++itHitP) {
                                        (*itHitP)->addToHit();
                                        const double charge = (*itHitP)->getCharge();
                                        localPadCenter = theModule.globalToLocal((*itHitP)->getX(), (*itHitP)->getY());
                                        sumCharge += charge;
                                        sumChargePhiMeas += charge * (*itHitP)->getPhiMeas();
                                        localPosition[0] += charge * localPadCenter[0];
                                        localPosition[1] += charge * localPadCenter[1];
                                }
                                Vector2D globalPosition2D = theModule.localToGlobal(localPosition[0] / sumCharge, localPosition[1] / sumCharge);
                                const double position[3] = { globalPosition2D[0], globalPosition2D[1], (*itP)->getZ() };
                                // calculate covariance matrix for hit position
                                float covMatrix[6] = { 0., 0., 0., 0., 0., 0. };
                                const double cosPhiMeas = cos(sumChargePhiMeas / sumCharge);
                                const double sinPhiMeas = sin(sumChargePhiMeas / sumCharge);
                                // variance in row direction
                                // Determine the covariance matrix according to the given resolution
                                const double zDrift = std::max(maxDriftLength - position[2], 0.);

                                // Sigma squared along Z
                                const double zSigmaSquare = _offsetDrift + _slopeDrift * zDrift;
                                // Sigma squared in RPhi plane (azimuthal (Column) direction)
                                const double phiSigmaSquare = _offsetCol + _slopeCol * zDrift;
                                // Sigma squared in RPhi plane (radial (row) direction)
                                const double rSigmaSquare = _offsetRow + _slopeRow * zDrift;

                                covMatrix[0] = pow(cosPhiMeas, 2) * phiSigmaSquare + pow(sinPhiMeas, 2) * rSigmaSquare;
                                covMatrix[1] = sinPhiMeas * cosPhiMeas * (phiSigmaSquare - rSigmaSquare);
                                covMatrix[2] = pow(sinPhiMeas, 2) * phiSigmaSquare + pow(cosPhiMeas, 2) * rSigmaSquare;
                                covMatrix[5] = float(zSigmaSquare);

                                TrackerHitImpl* theHit = new TrackerHitImpl();
                                //      theHit->setType(int type) ???
                                theHit->setPosition(position);
                                theHit->setCovMatrix(covMatrix);
                                //use the cellIDs of the seeding pad for the hit
                                theHit->setCellID1((*itP)->getPadId());
                                // encode CellID0
                                const int layer = row + _moduleRowOffset[moduleID];
                                encoder.reset();
                                encoder[ILDCellID0::subdet] = ILDDetID::TPC;
                                encoder[ILDCellID0::layer] = layer;
                                encoder[ILDCellID0::module] = moduleID;
                                theHit->setCellID0(encoder.lowWord());

                                // add pulses
                                for (pulseListType::iterator itHitP = hitPulses.begin(); itHitP != hitPulses.end(); ++itHitP) {
                                        TrackerPulse* aHitPulse = dynamic_cast<TrackerPulse*>(inputCol->getElementAt((*itHitP)->getPulseNum()));
                                        theHit->rawHits().push_back(aHitPulse);

                                        // determine the quality bits from the pulses to let hit inherit them
                                        int pulseQuality = aHitPulse->getQuality();
                                        if (pulseflag::isSplit(pulseQuality)) hitflag::setSplitPulse(qualityFlag); // split pulse?
                                        if (pulseflag::isOverflow(pulseQuality)) hitflag::setOverflowPulse(qualityFlag);  // overflow
                                        if (pulseflag::isAnomalousShape(pulseQuality)) hitflag::setAnomalousPulseShape(qualityFlag); // anomalous shape
                                        if (pulseflag::isPlateauCutoff(pulseQuality)) hitflag::setPlateauCutoffPulse(qualityFlag); // anomalous shape
                                        // variance of charge measurement
                                        varCharge += aHitPulse->getCovMatrix()[0];
                                }
                                theHit->setQuality(qualityFlag);
                                // charge
                                theHit->setEDep(_chargeConversionFactor * sumCharge);
                                theHit->setEDepError(_chargeConversionFactor * sqrt(varCharge));

                                // the time is inherited from the pulse with the largest charge
                                theHit->setTime((*itP)->getTime());

                                // create row based hit
                                rb_Hit* someHit = new rb_Hit(numHits, moduleID, layer, *theHit);
                                hitData.push_back(someHit);

                                outputHitCol->addElement(theHit);
                                numHits++;
                        }
                }
                for (seedListType::iterator itS = seedList.begin(); itS != seedList.end(); ++itS) {
                        //(*itS)->print();
                        if (!(*itS)->getValid()) continue;
                        hitListType seedHits;
                        pulseListType seedPulses = (*itS)->getPulseList();
                        for (pulseListType::iterator itP = seedPulses.begin(); itP != seedPulses.end(); ++itP) {
                                const int iHit = (*itP)->getHitNum();
                                if (iHit >= 0) seedHits.push_back(hitData[iHit]);
                        }
                        rb_Segment* aSegment = new rb_Segment(_Bzc, seedHits, 1);
                        streamlog_out(DEBUG) << " segment " << seedHits.size() << " " << aSegment->getChi2() << " " << aSegment->getNdf()
                                        << std::endl;
                        if (aSegment->getNdf() < 0)
                                delete aSegment;
                        else
                                segData[moduleID].push_back(aSegment);
                }

                // cleanup seeds
                for (seedListType::iterator itS = seedList.begin(); itS != seedList.end(); ++itS)
                        delete *itS;
        }

        std::vector<std::pair<int, int> > segIndex;
        simpleEquiClasses segEquiClass;
        // create segment index
        for (modSegMapType::iterator itMod = segData.begin(); itMod != segData.end(); ++itMod) {
                //std::cout << " Segments in module " << itMod->first << std::endl;
                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());
                        //std::cout << "  " << (*itSeg)->getId() << " " << (*itSeg)->getFirstRow() << " " << (*itSeg)->getLastRow() << " " << (*itSeg)->getHitList().size() << std::endl;
                }
        }
        // check segment matching with neighbour modules
        segMatchListType segMatches;
        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, 0)) {
                                                rb_SegmentMatch* someMatch = new rb_SegmentMatch(*itSeg1, *itSeg2);
                                                //someMatch->print();
                                                segMatches.push_back(someMatch);
                                        }
                                }
                        }
                }
        }
        // check for overlapping matches
        segMatchMapType forwardMatches, backwardMatches;
        for (segMatchListType::iterator itM = segMatches.begin(); itM != segMatches.end(); ++itM) {
                forwardMatches[(*itM)->getSegId1()].push_back((*itM));
                backwardMatches[(*itM)->getSegId2()].push_back((*itM));
        }
        // check forward matches
        for (segMatchMapType::iterator itSid = forwardMatches.begin(); itSid != forwardMatches.end(); ++itSid) {
                for (segMatchListType::iterator itM1 = itSid->second.begin(); itM1 != itSid->second.end() - 1; ++itM1) {
                        for (segMatchListType::iterator itM2 = itM1 + 1; itM2 != itSid->second.end(); ++itM2) {
                                if ((*itM1)->overlap(*itM2)) {
                                        (*itM1)->setInvalid();
                                        (*itM2)->setInvalid();
                                        streamlog_out(DEBUG) << " overlap of forward matches " << (*itM1)->getSegId1() << " " << (*itM1)->getSegId2() << ", "
                                                        << (*itM2)->getSegId1() << " " << (*itM2)->getSegId2() << std::endl;
                                }
                        }
                }
        }
        // check backward matches
        for (segMatchMapType::iterator itSid = backwardMatches.begin(); itSid != backwardMatches.end(); ++itSid) {
                for (segMatchListType::iterator itM1 = itSid->second.begin(); itM1 != itSid->second.end() - 1; ++itM1) {
                        for (segMatchListType::iterator itM2 = itM1 + 1; itM2 != itSid->second.end(); ++itM2) {
                                if ((*itM1)->overlap(*itM2)) {
                                        (*itM1)->setInvalid();
                                        (*itM2)->setInvalid();
                                        streamlog_out(DEBUG) << " overlap of backward matches " << (*itM1)->getSegId1() << " " << (*itM1)->getSegId2() << ", "
                                                        << (*itM2)->getSegId1() << " " << (*itM2)->getSegId2() << std::endl;
                                }
                        }
                }
        }
        // add valid matches
        for (segMatchListType::iterator itM = segMatches.begin(); itM != segMatches.end(); ++itM)
                if ((*itM)->getValid()) segEquiClass.addMatch(make_pair((*itM)->getSegId1(), (*itM)->getSegId2()));
        // 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) {
                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
        }
        streamlog_out(DEBUG) << "  found " << numHits << " hits, " << 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*>(outputHitCol->getElementAt((*itH)->getHitNum()));
                        theTrack->addHit(Hit);
                }

                streamlog_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;

                outputTrackCol->addElement(theTrack);
        }

        streamlog_out(DEBUG) << "RowBasedPadPulseRoadSearch ... done!" << endl;

        outputHitCol->parameters().setValue("CellIDEncoding", ILDCellID0::encoder_string);
        if (!outputHitCol->empty()) evt->addCollection(outputHitCol, _outputTrackerHitsColName);
        if (!outputTrackCol->empty()) evt->addCollection(outputTrackCol, _outputTrackColName);

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

        // clean up
        for (modPulseMapType::iterator itMod = pulseData.begin(); itMod != pulseData.end(); ++itMod) {
                for (rowPulseMapType::iterator itRow = itMod->second.begin(); itRow != itMod->second.end(); ++itRow) {
                        for (pulseListType::iterator itP = itRow->second.begin(); itP != itRow->second.end(); ++itP) {
                                delete *itP;
                        }
                }
        }
        pulseData.clear();
        for (hitListType::iterator itH = hitData.begin(); itH != hitData.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 (segMatchListType::iterator itM = segMatches.begin(); itM != segMatches.end(); ++itM)
                delete *itM;
        segMatches.clear();
}

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

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

}


bool RowBasedPadPulseRoadSearchProcessor::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.);
}


pb_Pulse::pb_Pulse(const int iPulse, const EVENT::TrackerPulse& aPulse, const gear::TPCModule& aModule, const double zPos) :
                _pulseNum(iPulse), _pulseModule(aPulse.getCellID1()), _pulsePadID(aPulse.getCellID0()), _pulseRow(
                                aModule.getRowNumber(_pulsePadID)), _pulseCol(aModule.getPadNumber(_pulsePadID)), _pulseQuality(aPulse.getQuality()), _pulseCharge(
                                aPulse.getCharge()), _pulseTime(aPulse.getTime()), _pulseX(aModule.getPadCenter(_pulsePadID)[0]), _pulseY(
                                aModule.getPadCenter(_pulsePadID)[1]), _pulseZ(zPos), _pulsePhiMeas(
                                (aModule.getLocalPadLayout().getCoordinateType() == gear::PadRowLayout2D::POLAR) ?
                                                atan2((_pulseX - aModule.getOffset()[0]), (aModule.getOffset()[1] - _pulseY)) : aModule.getAngle()), _pulseDirX(
                                cos(_pulsePhiMeas)), _pulseDirY(sin(_pulsePhiMeas)), _firstCol(_pulseCol), _lastCol(_pulseCol), _maxCharge(0.), _sumCharge(
                                _pulseCharge), _sumChargeDist(0.), _seeds(), _inHit(false), _indexCol(-1), _hitNum(-1) {
}

void pb_Pulse::print() const {
        std::cout << " Pulse " << _pulseModule << " " << _pulseRow << " " << _pulseCol << " " << _pulseCharge << " " << _pulseTime
                        << std::endl;
        std::cout << "       " << _pulseX << " " << _pulseY << " " << _pulseZ << " " << _pulsePhiMeas << " " << _firstCol << " "
                        << _lastCol << " " << _sumCharge << " " << _sumChargeDist << std::endl;
}

int pb_Pulse::getPulseNum() const {
        return _pulseNum;
}

int pb_Pulse::getModule() const {
        return _pulseModule;
}

int pb_Pulse::getPadId() const {
        return _pulsePadID;
}

int pb_Pulse::getRow() const {
        return _pulseRow;
}

int pb_Pulse::getColumn() const {
        return _pulseCol;
}

int pb_Pulse::getQuality() const {
        return _pulseQuality;
}

double pb_Pulse::getCharge() const {
        return _pulseCharge;
}

double pb_Pulse::getRelCharge() const {
        return (_lastCol + 1 - _firstCol) * _pulseCharge / _sumCharge;
}

double pb_Pulse::getTime() const {
        return _pulseTime;
}

double pb_Pulse::getX() const {
        return _pulseX;
}

double pb_Pulse::getY() const {
        return _pulseY;
}

double pb_Pulse::getZ() const {
        return _pulseZ;
}


void pb_Pulse::getPos(double* position) const {
        position[0] = _pulseX;
        position[1] = _pulseY;
        position[2] = _pulseZ;
}


void pb_Pulse::getCOG(double* position) const {
        position[0] = _pulseX + _pulseDirX * _sumChargeDist / _sumCharge;
        position[1] = _pulseY + _pulseDirY * _sumChargeDist / _sumCharge;
        position[2] = _pulseZ;
}

double pb_Pulse::getPhiMeas() const {
        return _pulsePhiMeas;
}


void pb_Pulse::getMeasDir(double* direction) const {
        direction[0] = _pulseDirX;
        direction[1] = _pulseDirY;
}


void pb_Pulse::addNeighbour(pb_Pulse* aPulse) {
        _firstCol = min(_firstCol, aPulse->getColumn());
        _lastCol = max(_lastCol, aPulse->getColumn());
        _maxCharge = max(_maxCharge, aPulse->getCharge());
        _sumCharge += aPulse->getCharge();
        _sumChargeDist += aPulse->getCharge() * (_pulseDirX * (aPulse->getX() - _pulseX) + _pulseDirY * (aPulse->getY() - _pulseY));
}

bool pb_Pulse::isMaxPulse() const {
        return (_pulseCharge >= _maxCharge);
}

bool pb_Pulse::inSeed() const {
        return (_seeds.size() > 0);
}

int pb_Pulse::getNumSeeds() const {
        return _seeds.size();
}


void pb_Pulse::addSeed(int s) {
        _seeds.push_back(s);
}


void pb_Pulse::removeSeed(int s) {
        _seeds.remove(s);
}

void pb_Pulse::clearSeeds() {
        _seeds.clear();
}

const seedIdListType& pb_Pulse::getSeedIdList() const {
        return _seeds;
}

int pb_Pulse::getSeedId() const {
        return _seeds.size() > 0 ? _seeds.front() : -1;
}

bool pb_Pulse::inHit() const {
        return _inHit;
}

void pb_Pulse::addToHit() {
        _inHit = true;
}

int pb_Pulse::getIndex() const {
        return _indexCol;
}


void pb_Pulse::setIndex(int i) {
        _indexCol = i;
}

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


void pb_Pulse::setHitNum(int i) {
        _hitNum = i;
}


pb_Seed::pb_Seed(const int seedId, const double Bzc, const pulseListType& pulses, const int minRows, const double rowDensityCut,
                const double maxVarXY) :
                _seedId(seedId), _bzc(Bzc), _checked(true), _numPulses(pulses.size()), _pulseList(pulses), _valid(false), _npar(
                                (Bzc != 0.) ? 5 : 4), _parameters(_npar) {
        // count number of rows
        sort(_pulseList.begin(), _pulseList.end(), comparePadRow);
        int numRows = 1;
        int firstRow = _pulseList.front()->getRow();
        int lastRow = firstRow;
        for (pulseListType::iterator itP = _pulseList.begin(); itP != _pulseList.end(); ++itP) {
                const int row = (*itP)->getRow();
                if (row > lastRow) {
                        numRows++;
                        lastRow = row;
                }
        }
        //std::cout << " pb_seed_checked " << _pulseList.size() << " " << numRows << " " << firstRow << " " << lastRow << std::endl;
        if (numRows < minRows or numRows < (lastRow + 1 - firstRow) * rowDensityCut) return;
        _fitSeed();
        // valid seed?
        _valid = (_varXY < maxVarXY);
}


pb_Seed::pb_Seed(const int seedId, const double Bzc, const pulseListType& pulses) :
                _seedId(seedId), _bzc(Bzc), _checked(false), _numPulses(pulses.size()), _pulseList(pulses), _valid(true), _npar(
                                (Bzc != 0.) ? 5 : 4), _parameters(_npar) {
        // count number of rows
        sort(_pulseList.begin(), _pulseList.end(), comparePadRow);
        int numRows = 1;
        int firstRow = _pulseList.front()->getRow();
        int lastRow = firstRow;
        for (pulseListType::iterator itP = _pulseList.begin(); itP != _pulseList.end(); ++itP) {
                const int row = (*itP)->getRow();
                if (row > lastRow) {
                        numRows++;
                        lastRow = row;
                }
        }
        //std::cout << " pb_seed_unchecked " << _pulseList.size() << " " << numRows << " " << firstRow << " " << lastRow << std::endl;
        _fitSeed();
}

void pb_Seed::_fitSeed() {
        // simple unweighted (sigma = 1.) fit with COG of largest pulses per row
        pulseListType fitPulses;
        pb_Pulse* largestPulse = _pulseList.front();
        for (pulseListType::iterator itP = _pulseList.begin(); itP != _pulseList.end(); ++itP) {
                if ((*itP)->getRow() == largestPulse->getRow()) {
                        if ((*itP)->getCharge() > largestPulse->getCharge()) largestPulse = *itP;
                } else {
                        fitPulses.push_back(largestPulse);
                        largestPulse = *itP;
                }
        }
        fitPulses.push_back(largestPulse);
        // reference point
        double pos1[3], pos2[3];
        fitPulses.front()->getPos(pos1);
        fitPulses.back()->getPos(pos2);
        _refX = 0.5 * (pos1[0] + pos2[0]);
        _refY = 0.5 * (pos1[1] + pos2[1]);
        _refZ = 0.5 * (pos1[2] + pos2[2]);
        double refPoint[] = { _refX, _refY, _refZ };
        int nParXY, nPointXY, nParZS, nPointZS;
        double pos[3], sArc, chi2XY, chi2ZS;
        double parameters[] = { 0., 0., 0., 0., 0. };
        // coarse unweighted xy fit
        simpleFitXY xyFit((_bzc != 0.), _refX, _refY);
        for (pulseListType::iterator itP = fitPulses.begin(); itP != fitPulses.end(); ++itP) {
                (*itP)->getCOG(pos);
                xyFit.addPoint(pos[0], pos[1], 1.0);
        }
        nParXY = xyFit.fit(chi2XY, nPointXY);
        //std::cout << " xyFit " << chi2XY << " " << nPointXY << " " << nParXY << std::endl;
        _varXY = chi2XY / (nPointXY - nParXY);
        _parameters.SetSub(0, xyFit.getPar());
        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 zsFit;
        for (pulseListType::iterator itP = fitPulses.begin(); itP != fitPulses.end(); ++itP) {
                (*itP)->getCOG(pos);
                sArc = hlx.getArcLengthXY(pos);
                zsFit.addPoint(sArc, pos[2] - _refZ, 1.0);
        }
        nParZS = zsFit.fit(chi2ZS, nPointZS);
        //std::cout << " zsFit " << chi2ZS << " " << nPointZS << std::endl;
        _varZS = chi2ZS / (nPointZS - nParZS);
        _parameters.SetSub(nParXY, zsFit.getPar());
}

void pb_Seed::print() const {
        std::cout << "  seed " << _seedId << " " << _checked << ":" << _valid << " " << _pulseList.size() << " "
                        << this->getNumUniqueSeedPulses() << " " << _pulseList.front()->getRow() << " " << _pulseList.front()->getColumn() << " "
                        << _pulseList.back()->getRow() << " " << _pulseList.back()->getColumn() << " " << _varXY << " " << _varZS << std::endl;
}

int pb_Seed::getSeedId() const {
        return _seedId;
}

bool pb_Seed::getValid() const {
        return _valid;
}

void pb_Seed::setInvalid() {
        _valid = false;
}

int pb_Seed::getNumSeedPulses() const {
        return _numPulses;
}

int pb_Seed::getNumUniqueSeedPulses() const {
        int num = 0;
        for (pulseListType::const_iterator itP = _pulseList.begin(); itP != _pulseList.end(); ++itP) {
                if ((*itP)->getNumSeeds() < 2) num++;
        }
        return num;
}

bool pb_Seed::hasAddedPulses() const {
        return (_pulseList.size() > _numPulses);
}


void pb_Seed::addPulse(pb_Pulse* aPulse) {
        _pulseList.push_back(aPulse);
}


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


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

double pb_Seed::getVarXY() const {
        return _varXY;
}

void pb_Seed::markPulses() {
        for (pulseListType::iterator itP = _pulseList.begin(); itP != _pulseList.end(); ++itP) {
                (*itP)->addSeed(_seedId);
        }
}

void pb_Seed::freePulses() {
        for (pulseListType::iterator itP = _pulseList.begin(); itP != _pulseList.end(); ++itP) {
                (*itP)->removeSeed(_seedId);
        }
}

const pulseListType& pb_Seed::getPulseList() const {
        return _pulseList;
}


void pb_Seed::roadSearch(pulseListType& pulses, pulseListType& insidePulses, const double maxDxy, const double maxDz) const {
        insidePulses.clear();
        if (!_valid) return;
        // construct helix
        double refPoint[3], par[5];
        this->getRefPoint(refPoint);
        this->getPar(par);
        simpleHelix hlx(par, refPoint, _bzc);
        // test pulses
        double pos[3], xyPos, zPos, sArc, phi, lambda;
        for (pulseListType::iterator itP = pulses.begin(); itP != pulses.end(); ++itP) {
                (*itP)->getCOG(pos);
                if (hlx.getExpectedPlanePos(pos, (*itP)->getPhiMeas(), xyPos, zPos, sArc, phi, lambda)) {
                        //std::cout << " expected pos " << xyPos << " " << zPos << std::endl;
                        if (abs(xyPos) < maxDxy and abs(zPos) < maxDz) insidePulses.push_back(*itP);
                }
        }
        //std::cout << " road search " << pulses.size() << " " << insidePulses.size() << std::endl;
}


rb_SegmentMatch::rb_SegmentMatch(const rb_Segment* segment1, const rb_Segment* segment2) :
                _segId1(segment1->getId()), _segId2(segment2->getId()), _firstRow1(segment1->getFirstRow()), _firstRow2(
                                segment2->getFirstRow()), _lastRow1(segment1->getLastRow()), _lastRow2(segment2->getLastRow()), _valid(true) {
}

void rb_SegmentMatch::print() const {
        std::cout << "   Segment match " << _segId1 << " " << _segId2 << ", " << _firstRow1 << " " << _firstRow2 << ", " << _lastRow1
                        << " " << _lastRow2 << ", " << _valid << std::endl;
}

int rb_SegmentMatch::getSegId1() const {
        return _segId1;
}

int rb_SegmentMatch::getSegId2() const {
        return _segId2;
}

int rb_SegmentMatch::getFirstRow1() const {
        return _firstRow1;
}

int rb_SegmentMatch::getFirstRow2() const {
        return _firstRow2;
}

int rb_SegmentMatch::getLastRow1() const {
        return _lastRow1;
}

int rb_SegmentMatch::getLastRow2() const {
        return _lastRow2;
}

bool rb_SegmentMatch::getValid() const {
        return _valid;
}

bool rb_SegmentMatch::overlap(rb_SegmentMatch* aMatch) const {
        const int first1 = max(_firstRow1, aMatch->getFirstRow1());
        const int first2 = max(_firstRow2, aMatch->getFirstRow2());
        const int last1 = min(_lastRow1, aMatch->getLastRow1());
        const int last2 = min(_lastRow2, aMatch->getLastRow2());
        return first1 <= last1 and first2 <= last2;
}

void rb_SegmentMatch::setInvalid() {
        _valid = false;
}

} // close namespace brackets

// #endif //#ifdef USE_ROWPPRS