RowBasedHitFinderProcessor.cc

Go to the documentation of this file.

#include "RowBasedHitFinderProcessor.h"

//C++
#include <iostream>
#include <algorithm>

//LCIO
#include <EVENT/LCFloatVec.h>
#include <IMPL/TrackerDataImpl.h>
#include <IMPL/TrackerPulseImpl.h>
#include <IMPL/TrackerHitImpl.h>
#include <IMPL/LCFlagImpl.h>

// GEAR
#include <gear/TPCParameters.h>
#include <gear/TPCModule.h>
#include <gearimpl/GlobalPadIndex.h>

//MARLIN
#include <marlin/Global.h>

//Streamlog
#include <streamlog/streamlog.h>

//lccd
#include <lccd_exceptions.h>

// the exception by TPCCondData in MarlinTPC
#include "TPCCondDataExeptions.h"
#include "TPCConditions.h"

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

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

namespace marlintpc
{

RowBasedHitFinderProcessor aRowBasedHitFinderProcessor;


RowBasedHitFinderProcessor::RowBasedHitFinderProcessor() : Processor("RowBasedHitFinderProcessor"), _channelCorrectionListener(NULL), _tpcConditionsListener(NULL)

{
        _description = "RowBasedHitFinderProcessor calculates TrackerHits from TrackerPulses" ;

        // register steering parameters: name, description, class-variable, default value
        registerInputCollection(LCIO::TRACKERPULSE,
                        "InputTrackerPulses" ,
                        "Name of the input TrackerPulses collection"  ,
                        _inputTrackerPulsesCollectionName ,
                        string("TPCPulses")) ;

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

        registerOptionalParameter("InputTPCConditions",
                        "Optional: Name of the input collection containing the TPC conditions data",
                        _inputTPCConditionsCollectionName,
                        TPCConditions::getDefaultColName()) ;

        registerProcessorParameter("WriteOutputToStorage" ,
                        "if true, the output hits collection is written to storage (default: true)"  ,
                        _outputHitsPersistent,
                        true);

        registerProcessorParameter("ChannelCorrectionCollectionName" ,
                        "Name of the conditions collection in which the channel correction is stored",
                        _channelCorrectionCollectionName,
                        tpcconddata::ChannelCorrection::getDefaultColName());

        registerOptionalParameter("MaxEmptyPadsOverride" ,
                        "This parameter is depreciated and will be removed in the next release."  ,
                        _deprecated ,
                        1) ;

        registerOptionalParameter("MaxDeadChannelsOverride" ,
                        "Maximum number of empty consecutive pads in hit (default: 1)"  ,
                        _maxDeadChannelsOverride ,
                        1) ;

        registerOptionalParameter("MinHitSizeOverride" ,
                        "Minimum size of hit / Minimum number of Pads (default: 1)"  ,
                        _minHitSizeOverride ,
                        1) ;

        registerOptionalParameter("MinMaxPulseSliceHeightOverride",
                        "Minimum size of the maximum slice value of the maximal pulse (in charge) in a hit.",
                        _minMaxPulseSliceHeightOverride,
                        float(12.)) ;

        registerOptionalParameter("MaxTimeSpreadOverride" ,
                        "Maximum time between pulses in a hit in ns (default: 200.)"  ,
                        _maxTimeSpreadOverride ,
                        float(200.)) ;

        registerOptionalParameter("VDriftOverride" ,
                        "Optional: Set drift velocity in case there is no conditions data in mm/mu s"  ,
                        _driftVelocityOverride ,
                        float(0.)) ;

        registerOptionalParameter("DeadPadsOverride",
                        "Optional: the gear pad indices of dead channels",
                        _deadChannelIndicesOverride,
                        _deadChannelIndicesOverride) ;

        registerProcessorParameter("ChargeConversionIsCalibrated" ,
                        "Set true, if the charge conversion factor from ADC values to primary electrons is known (default:true).",
                        _chargeIsCalibrated,
                        bool(true)) ;

        registerOptionalParameter("ChargeConversionFactorOverride" ,
                        "Optional: The charge conversion factor from ADC values to primary electrons",
                        _chargeConversionFactorOverride ,
                        float(1.)) ;
}



void RowBasedHitFinderProcessor::init()
{
        printParameters();
        // get drift velocity, either measured value from calibration
        // or nominal value given as process parameter
        if(parameterSet("VDriftOverride"))
        {
                streamlog_out(WARNING) << "the drift velocity is taken from the steering file. please change" << endl;
        }
        else // CAUTION! this is not only true for this value, but currently it's the only one in use
        {
                try
                {
                        _tpcConditionsListener = new TPCConditionsListener(_inputTPCConditionsCollectionName);
                }

                catch(lccd::LCCDException &)
                {
                        streamlog_out(ERROR) << "Collection named \"" << _inputTPCConditionsCollectionName << "\" with TPCConditions is not available." << std::endl;
                        streamlog_out(ERROR) << "Use the ConditionsProcessor to provide it, or set the VDriftOverride parameter if this is not possible." << std::endl;
                        throw;
                }
        }

        if(parameterSet("DeadPadsOverride"))
        {
                streamlog_out(WARNING) << "the indices of dead pads are taken from the steering file. please change" << endl;
        }
        else
        {
                try
                {
                        _channelCorrectionListener = new ChannelCorrectionListener(_channelCorrectionCollectionName);
                }
                catch(lccd::LCCDException &)
                {
                        streamlog_out(ERROR) << "Collection named \"" << _channelCorrectionCollectionName << "\" with ChannelCorrections is not available." << std::endl;
                        streamlog_out(ERROR) << "Use the ConditionsProcessor to provide it, or set the DeadPadsOverride parameter if this is not possible." << std::endl;
                        throw;
                }
        }

        if(parameterSet("ChargeConversionFactorOverride"))
        {
                streamlog_out(WARNING) << "the conversion factor from adc counts to charge is taken from the steering file. please change" << endl;
        }
        else if(_tpcConditionsListener == 0)
        {
                try
                {
                        _tpcConditionsListener = new TPCConditionsListener(_inputTPCConditionsCollectionName);
                }
                catch(lccd::LCCDException &)
                {
                        streamlog_out(ERROR) << "Collection named \"" << _inputTPCConditionsCollectionName << "\" with TPCConditions is not available." << std::endl;
                        streamlog_out(ERROR) << "Use the ConditionsProcessor to provide it, or set the ChargeConversionFactorOverride parameter if this is not possible." << std::endl;
                        throw;
                }
        }

        if(parameterSet("MaxEmptyPadsOverride"))
        {
                streamlog_out(WARNING) << "The parameter [MaxEmptyPadsOverride] is deprecated and has no effect! It will be removed in the next release. Please use [MaxDeadChannelsOverride]." << endl;
        }
}



void RowBasedHitFinderProcessor::processRunHeader(LCRunHeader* run)
{
        run->parameters().setValue(_processorName + "_revision", "$Rev: 1842 $");

        for(ProcParamMap::iterator i = _map.begin(); i != _map.end(); i++)
        {
                if(!i->second->isOptional() || i->second->valueSet())
                        run->parameters().setValue(_processorName + "_" + i->second->name(), i->second->value());
        }
}



void RowBasedHitFinderProcessor::processEvent(LCEvent* evt)
{
        // get the input pulse collection
        LCCollection* inputPulses;
        try
        {
                inputPulses = evt->getCollection(_inputTrackerPulsesCollectionName);
        }
        catch(DataNotAvailableException &e)
        {
                // if no pulses are available, skip to next event
                streamlog_out(DEBUG) << "No pulses available in event " << evt->getEventNumber() << endl;
                return;
        }

        // create output hit collection and set the persistency flag
        LCCollectionVec* tpcHitCollection = new LCCollectionVec(LCIO::TRACKERHIT);
        tpcHitCollection->setTransient(!_outputHitsPersistent);


        // set the flags that cellID1 should be stored
        lcio::LCFlagImpl trkFlag(0);
        trkFlag.setBit(lcio::LCIO::RTHBIT_ID1);
        tpcHitCollection->setFlag(trkFlag.getFlag());

        // the main parameter for the geometry description is the gear parameter:
        const gear::TPCParameters& tpcParameters = Global::GEAR->getTPCParameters();

        map<int, map<int, vector<TrackerPulse*>*>*> theRowSortedPulsesPerModule;
        map<int, vector<hitCandidate> >             theHitCandidatesPerModule;


        // (A): sort pulses
        _sortPulses(inputPulses, tpcParameters, theRowSortedPulsesPerModule);

        // (B)
        _findHitCandidates(theRowSortedPulsesPerModule, tpcParameters, theHitCandidatesPerModule);

        // Clean up the maps
        map<int, map<int, vector<TrackerPulse*>* >* >::iterator aModuleWithPulses;
        map<int, vector<TrackerPulse*>* >::iterator aRowWithPulses;
        // outer loop:: the modules
        for(aModuleWithPulses = theRowSortedPulsesPerModule.begin();
                        aModuleWithPulses != theRowSortedPulsesPerModule.end(); ++aModuleWithPulses)
        {
                // inner loop:: the rows
                for(aRowWithPulses  = aModuleWithPulses->second->begin();
                                aRowWithPulses != aModuleWithPulses->second->end(); ++aRowWithPulses)
                {
                        delete aRowWithPulses->second;
                }
                delete aModuleWithPulses->second;
        }

        // (C)
        _createHitCollection(theHitCandidatesPerModule, tpcParameters, tpcHitCollection);

        evt->addCollection(tpcHitCollection, _outputTrackerHitsCollectionName);

        streamlog_out(DEBUG4) << "Found " << tpcHitCollection->getNumberOfElements() << " hits in the event." << endl;

        streamlog_out(DEBUG3) << "There are " << inputPulses->getNumberOfElements()
                                  << " pulses to begin with, and a total of "
                                  << tpcHitCollection->getNumberOfElements() << " hits found." << endl;
}



void RowBasedHitFinderProcessor::check(LCEvent* evt)
{
}



void RowBasedHitFinderProcessor::end()
{
        if(_channelCorrectionListener != 0)  // release memory, if needed
                delete _channelCorrectionListener;

}



void RowBasedHitFinderProcessor::_sortPulses(LCCollection* inputPulses, const gear::TPCParameters& tpcParameters,
                map<int, map<int, vector<TrackerPulse*>* >* >& rowSortedPulsesPerModule)
{
        pair<map<int, map<int, vector<TrackerPulse*>*>*>::iterator, bool> moduleInserter;
        pair<map<int, vector<TrackerPulse*>*>::iterator, bool> rowInserter;

        for(int aPulse = 0; aPulse < inputPulses->getNumberOfElements(); ++aPulse)
        {
                TrackerPulse* inputPulse = dynamic_cast<TrackerPulse*>(inputPulses->getElementAt(aPulse));
                if(inputPulse&&(!_padIsDead(inputPulse->getCellID0(),inputPulse->getCellID1())))  // the cast worked
                                {
                        const int moduleID         = inputPulse->getCellID1();
                        const int padID            = inputPulse->getCellID0();

                        const TPCModule& theModule = tpcParameters.getModule(moduleID);
                        const int rowID            = theModule.getRowNumber(padID);

                        // first create/use the module insertion iterator
                        map<int, vector<TrackerPulse*>*>* rmap = new map<int, vector<TrackerPulse*>*>;
                        moduleInserter =
                                        rowSortedPulsesPerModule.insert(pair<int,  map<int, vector<TrackerPulse*>* >* >
                        (moduleID, rmap));
                        if(!moduleInserter.second)
                                delete rmap;
                        // use the newly created/already existing map pointer to insert the row
                        vector<TrackerPulse*>* pmap = new vector<TrackerPulse*>;
                        rowInserter = ((moduleInserter.first)->second)->insert
                                        (pair<int, vector<TrackerPulse*>* >(rowID, pmap));
                        if(!rowInserter.second)
                                delete pmap;

                        // and finally add the Pulse itself
                        ((rowInserter.first)->second)->push_back(inputPulse);
                                }
        } // end input pulse loop
}



void RowBasedHitFinderProcessor::_findHitCandidates
(map<int, map<int, vector<TrackerPulse*>* >* >& rowPulsesPerModule, const gear::TPCParameters& tpcParameters,
                map<int, vector<hitCandidate> >& theHitCandidateMap)
{
        // two iterator types are needed: for the modules and for the rows
        map<int, map<int, vector<TrackerPulse*>* >* >::iterator aModuleWithPulses;
        map<int, vector<TrackerPulse*>* >::iterator aRowWithPulses;

        // outer loop:: the modules
        for(aModuleWithPulses = rowPulsesPerModule.begin(); aModuleWithPulses != rowPulsesPerModule.end(); ++aModuleWithPulses)
        {
                int moduleID = 0;
                vector<hitCandidate> theHitCandidates;

                // inner loop:: the rows
                for(aRowWithPulses  = aModuleWithPulses->second->begin(); aRowWithPulses != aModuleWithPulses->second->end(); ++aRowWithPulses)
                {
                        streamlog_out(DEBUG) << "Starting search of neighbouring pulses in in row: " << aRowWithPulses->first << endl;

                        // third loop: now row wise, find hit candidates and remove the pulses consecutively
                        while(aRowWithPulses->second->size() != 0)
                        {
                                // first find and remove the pulse with the highest charge from the list
                                TrackerPulse* theSeedPulse = _findHighestPulseInRow(*(aRowWithPulses->second));
                                moduleID = theSeedPulse->getCellID1();

                                streamlog_out(DEBUG) << "Found a highest pulse with charge = " << theSeedPulse->getCharge()
                                                         << " on pad " << theSeedPulse->getCellID0() << " and with time " << theSeedPulse->getTime() << endl;
                                // time saver, conclude the search,
                                // if the maximal pulse height is below the minimum requirement
                                if(!_maxPulseHeightCheck(theSeedPulse))
                                {
                                        streamlog_out(DEBUG) << "Pulse is rejected, since it is too small. Continue in next row. " << endl;
                                        break;
                                }

                                // take the found pulse as prototype hit
                                hitCandidate theProtoHit(theSeedPulse);

                                // look in the adjacent pads for further pulses, store them in a vector
                                // this function modifies the initial list by removing pulses!
                                // check if it's needed before
                                if(aRowWithPulses->second->size() != 0)
                                        _addNeighbourPulses(theProtoHit, tpcParameters, *(aRowWithPulses->second));

                                theHitCandidates.push_back(theProtoHit);
                        }
                }
                theHitCandidateMap.insert(pair<int, vector<hitCandidate> >(moduleID, theHitCandidates));

                streamlog_out(DEBUG) << "On module " << aModuleWithPulses->first
                                << " there are " << theHitCandidates.size() << " hits found." << endl;

        } // end outer loop  of modules
}



TrackerPulse* RowBasedHitFinderProcessor::_findHighestPulseInRow(vector<TrackerPulse*>& pulsesInARow)
{
        TrackerPulse* highestPulse = pulsesInARow.at(0);

        vector<TrackerPulse*>::iterator eraser = pulsesInARow.begin();
        for(vector<TrackerPulse*>::iterator iter = pulsesInARow.begin(); iter < pulsesInARow.end(); ++iter)
        {
                if(fabs((*iter)->getCharge()) > fabs(highestPulse->getCharge()))
                {
                        highestPulse = *iter;
                        eraser = iter;
                }
        }
        pulsesInARow.erase(eraser);
        return highestPulse;
}



bool RowBasedHitFinderProcessor::_maxPulseHeightCheck(TrackerPulse* theSeedPulse)
{
        EVENT::TrackerData*  trackerData = theSeedPulse->getTrackerData();
        vector<float> theADCSpectrum = trackerData->getChargeValues();

        return *max_element(theADCSpectrum.begin(), theADCSpectrum.end()) >= _getMinMaxPulseSliceHeight();
}



void RowBasedHitFinderProcessor::_addNeighbourPulses
(hitCandidate& theCandidate,  const gear::TPCParameters& tpcParameters, vector<TrackerPulse*>& pulsesInARow)
{
        // the initial values for searching
        const int theStartPadIndex      = (theCandidate.getMaximumPulse())->getCellID0();
        const int theModuleIndex        = (theCandidate.getMaximumPulse())->getCellID1();
        const TPCModule& theModule      = tpcParameters.getModule(theModuleIndex);
        const float theStartTime        = (theCandidate.getMaximumPulse())->getTime();
        const float theHighestCharge    = theCandidate.getMaxPulseHeight();

        // the temporary vars for searching
        unsigned int numberOfDeadChannels = 0;
        float tempCharge            = theHighestCharge;
        int tempPadIndex            = theStartPadIndex;
        float tempTime              = theStartTime;
        bool finishedLeft           = false; // arbitrary choice: start with searching "left", used to change search direction
        bool stopSearch             = false; // loop control

        // start loop until no further pulses are found
        while(!stopSearch)
        {
                try // to get the index of the next pad
                {
                        if(!finishedLeft)
                                tempPadIndex = theModule.getLeftNeighbour(tempPadIndex);
                        else
                                tempPadIndex = theModule.getRightNeighbour(tempPadIndex);
                }
                catch(...)  // FIXME: what kind of exception is thrown?
                {
                        streamlog_out(DEBUG) << "The hit candidate is at the border!" << endl;
                        hitflag::setAtModuleBorder(theCandidate.qualityFlag);

                        if(!finishedLeft)
                        {
                                finishedLeft  = true;
                                tempPadIndex  = theStartPadIndex;
                        }
                        else
                                stopSearch = true;
                }


                if(_padIsDead(tempPadIndex, theModuleIndex))  // no pulse found because of dead pad
                {
                        numberOfDeadChannels++; // increase the counter
                        //std::cout<<"found dead pad: module="<<theModuleIndex<<"; row="<< theModule.getRowNumber(tempPadIndex)<<"; pad="<<theModule.getPadNumber(tempPadIndex)<<std::endl;
                }
                else
                {
                        // if more than one pulse matches the requirements, the best match needs to be found
                        // store this in a temporary vector that is searched for the best match
                        vector<vector<TrackerPulse*>::iterator> tempPulses;
                        tempPulses.clear();

                        // check the pulses if (a) they are on the right pad and (b) in the right time window
                        for(vector<TrackerPulse*>::iterator iter = pulsesInARow.begin(); iter < pulsesInARow.end(); ++iter)
                        {
                                if(((*iter)->getCellID0() == tempPadIndex) && (fabs((*iter)->getTime() - tempTime) <  _getMaxTimeSpread()))
                                        tempPulses.push_back(iter);
                        }

                        // now make distinction on the number of matched pulses
                        if(tempPulses.size() > 0)  // if several pulses can be matched, find the best match!
                        {
                                if(tempPulses.size() > 1)  // possibly a multiple hit in z !?
                                        hitflag::setMultipleHitCandidate(theCandidate.qualityFlag);

                                // this function returns the iterator directly, if only one pulse is passed to it
                                vector<TrackerPulse*>::iterator bestPulseMatch = _findBestNeighbourPulseMatch(tempPulses, tempTime);

                                theCandidate.addPulse(*bestPulseMatch); // store the pulse into hit candidate
                                tempTime = (*bestPulseMatch)->getTime(); // set the new time reference
                                pulsesInARow.erase(bestPulseMatch); // delete the stored pulse from list to disallow double-counting
                                tempPulses.clear();

                                // do a quick check if there is a rise in charge -> indication for multiple hits
                                float thisCharge = (*bestPulseMatch)->getCharge();
                                if(fabs(thisCharge) > fabs(tempCharge))
                                        hitflag::setMultipleHitCandidate(theCandidate.qualityFlag);
                                tempCharge = thisCharge;

                                if(numberOfDeadChannels > 0)  // there was at least one dead channel in between
                                {
                                        hitflag::setDeadChannel(theCandidate.qualityFlag);
                                        numberOfDeadChannels = 0; // reset the counter, if needed
                                }
                        }
                        else // nothing found
                        {
                                tempCharge = 0;

                                /* in two cases the search will be abandoned:
                        /            a) the pad isn't dead
                        /            b) the allowed number of consecutive dead channels is exceeded */
                                // case a) the pad isn't dead -> "end"

                                if(numberOfDeadChannels > 0)
                                        hitflag::setNextToDeadChannel(theCandidate.qualityFlag);

                                if(!finishedLeft)  // but if the other direction isn't searched yet -> start from scratch
                                {
                                        finishedLeft   = true;
                                        tempTime       = theStartTime;
                                        numberOfDeadChannels = 0;
                                        tempPadIndex   = theStartPadIndex;
                                        tempCharge            = theHighestCharge;
                                }
                                else // and end
                                        stopSearch = true;

                        }
                }

                // cover case b); unfortunately this can only be done *after* evaluation if the pad is dead
                // the "end" code is duplicated; the only other way is a goto (as far as i can think now)
                if(numberOfDeadChannels > _getMaxDeadChannels())  // this is a dead end!
                {
                        hitflag::setNextToDeadChannel(theCandidate.qualityFlag);
                        if(!finishedLeft)  // but if the other direction isn't searched yet -> start from scratch
                        {
                                finishedLeft   = true;
                                tempTime       = theStartTime;
                                numberOfDeadChannels = 0;
                                tempPadIndex   = theStartPadIndex;
                                tempCharge      = theHighestCharge;
                        }
                        else // and end
                                stopSearch = true;
                }
        }// search loop
        streamlog_out(DEBUG) << "A hit candidate with a width of " << theCandidate.getWidth() << " is found." << endl;
}



vector<TrackerPulse*>::iterator RowBasedHitFinderProcessor::_findBestNeighbourPulseMatch
(vector<vector<TrackerPulse*>::iterator>& pulseMatches, float theTime)
{
        // searches pulses on neighboring pads
        // if only one is found, this is taken
        // if there are multiple candidates, the one with the least time difference is taken
        if(pulseMatches.size() == 1)
                return pulseMatches.at(0);
        else
        {
                vector<TrackerPulse*>::iterator bestMatch = *(pulseMatches.begin());
                float minTime = fabs((*bestMatch)->getTime() - theTime);

                for(vector<vector<TrackerPulse*>::iterator>::iterator searchIter = pulseMatches.begin();
                                searchIter < pulseMatches.end(); ++searchIter)
                {
                        if(fabs((**searchIter)->getTime() - theTime) < minTime)
                        {
                                minTime = fabs((**searchIter)->getTime() - theTime);
                                bestMatch = *searchIter;
                        }
                }
                streamlog_out(DEBUG) << "In the hit more than one pulse was found to match to the protohit,"
                                "the best match in time was chosen." << endl;
                return bestMatch;
        }
}



void RowBasedHitFinderProcessor::_createHitCollection
(map<int, vector<hitCandidate> >& theHitCandidateMap, const gear::TPCParameters& tpcParameters, LCCollectionVec* tpcHitCollection)
{
        map<int, vector<hitCandidate> >::iterator hitCandidatesPerModule;

        // loop over modules
        for(hitCandidatesPerModule = theHitCandidateMap.begin();
                        hitCandidatesPerModule != theHitCandidateMap.end(); ++hitCandidatesPerModule)
        {
                //loop over hitCandidates:
                for(vector<hitCandidate>::iterator candidate = hitCandidatesPerModule->second.begin();
                                candidate < hitCandidatesPerModule->second.end(); ++candidate)
                {
                        // check validity in terms of what is defined in the steering file
                        if(_hitIsValid(*candidate) == false)
                                continue;

                        // calculate the center of gravity in both measurement plane coordinates
                        float covMatrix[6] = {0., 0., 0., 0., 0., 0.};
                        Vector3D positionGlobal3D = _calculatePositionsAndErrors(*candidate, tpcParameters, covMatrix);

                        const double position[3] = { positionGlobal3D[0], positionGlobal3D[1], positionGlobal3D[2] };

                        TrackerHitImpl* theHit = new TrackerHitImpl();
                        //      theHit->setType(int type) ???
                        theHit->setPosition(position);
                        theHit->setCovMatrix(covMatrix);
                        theHit->setEDep(_convertCharge(*candidate));
                        //  streamlog_out(WARNING)<< " Position standard deviations [" << sqrt(theHit->getCovMatrix().at(0)) << "|" << sqrt(theHit->getCovMatrix().at(2)) << "|" <<  sqrt(theHit->getCovMatrix().at(5))  << "]" << endl;
                        if(_chargeIsCalibrated)
                        {
                                theHit->setEDep(_convertCharge(*candidate));
                                theHit->setEDepError(sqrt(pow(_getChargeConversionFactor(),2)*candidate->getChargeMeasError()+_convertCharge(*candidate)));
                        }
                        else
                        {
                                theHit->setEDep(candidate->getCharge());
                                theHit->setEDepError(sqrt(candidate->getChargeMeasError()));
                        }
                        //use the cellIDs of the pad, which is closest to the hit position, for the hit
                        theHit->setCellID0((tpcParameters.getNearestPad( position[0],position[1])).getModuleID());
                        theHit->setCellID1((tpcParameters.getNearestPad( position[0],position[1])).getPadIndex());

                        vector<TrackerPulse*> hitPulses = candidate->getPulses();
                        for(vector<TrackerPulse*>::const_iterator aHitPulse = hitPulses.begin();
                                        aHitPulse < hitPulses.end(); ++aHitPulse)
                        {
                                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))  // split pulse?
                                        hitflag::setSplitPulse(candidate->qualityFlag);
                                if(pulseflag::isOverflow(pulseQuality))  // overflow
                                        hitflag::setOverflowPulse(candidate->qualityFlag);
                                if(pulseflag::isAnomalousShape(pulseQuality))  // anomalous shape
                                        hitflag::setAnomalousPulseShape(candidate->qualityFlag);
                                if(pulseflag::isPlateauCutoff(pulseQuality))  // anomalous shape
                                        hitflag::setPlateauCutoffPulse(candidate->qualityFlag);

                        }
                        theHit->setQuality(candidate->qualityFlag);

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

                        streamlog_out(DEBUG) << "The hit coordinates are [" << position[0] << "|" << position[1] << "|" << position[2] << "]" << endl
                                        << " with standard deviations [" << sqrt(theHit->getCovMatrix().at(0)) << "|" << sqrt(theHit->getCovMatrix().at(2)) << "|" <<  sqrt(theHit->getCovMatrix().at(5))  << "]" << endl
                                        << " with variances [" << theHit->getCovMatrix().at(0) << "|" << theHit->getCovMatrix().at(2) << "|" <<  theHit->getCovMatrix().at(5)  << "]" << endl
                                        << ", the charge is " << theHit->getEDep() << " +- " << theHit->getEDepError() << " and has the quality byte: ["  << theHit->getQuality() << "]" << endl;
                        streamlog_out(DEBUG) << "The detailed list of hit flags: " <<  " has dead channel [" <<  hitflag::hasDeadChannel(theHit->getQuality())
                        << "], has noisy channel [" <<  hitflag::hasNoisyChannel(theHit->getQuality())
                        << "], is at the border of a module [" <<  hitflag::isAtModuleBorder(theHit->getQuality())
                        << "], has an overflow pulse [" <<  hitflag::hasOverflowPulse(theHit->getQuality())
                        << "], contains a split pulse [" <<  hitflag::hasSplitPulse(theHit->getQuality())
                        << "], has an anomalous shaped pulse [" <<  hitflag::hasAnomalousPulseShape(theHit->getQuality())
                        << "], has an plateau pulse [" <<  hitflag::hasPlateauCutoffPulse(theHit->getQuality())
                        << "], is next to a dead channel [" <<  hitflag::isNextToDeadChannel(theHit->getQuality())
                        << "], and is a multiple hit candidate [" <<  hitflag::isMultipleHitCandidate(theHit->getQuality()) << "]." << endl;


                        tpcHitCollection->addElement(theHit);
                }
        }
        return;
}



bool RowBasedHitFinderProcessor::_hitIsValid(const hitCandidate& cand)
{
        return (int)cand.getWidth() >= _getMinHitSize();
}



Vector3D RowBasedHitFinderProcessor::_calculatePositionsAndErrors(const hitCandidate& cand, const gear::TPCParameters& tpcParameters, float* covMatrix)
{
        // calculate the center of gravity for the measurement coordinate
        const TrackerPulse* maxPulse              = cand.getMaximumPulse();
        const gear::TPCModule& module             = tpcParameters.getModule(maxPulse->getCellID1());
        const gear::PadRowLayout2D& localPadPlane = module.getLocalPadLayout();

        // check whether the module geometry is described in cartesian coordinates
        // if not, then throw an exception; this isn't in the manual -- or: the convention is to use cartesian coordinates in LCIO and the global geometry
        if(module.getCoordinateType() != PadRowLayout2D::CARTESIAN)
                throw EVENT::Exception("The module does not return cartesian coordinates; this is against the general convention. Better fix it and catch the exception.");

        // now check the coordinate system type of the pad plane; this defines the natural system for all calculations
        // define two bools for easier reading
        bool cartesianPadPlane = (localPadPlane.getCoordinateType()  == PadRowLayout2D::CARTESIAN);
        bool polarPadPlane = !cartesianPadPlane;

        /*** NOTE: if the pad plane is polar, then the globalToLocal transformation also changes the coordinate type in GEAR !!! ***/

        vector<TrackerPulse*> thePulses = cand.getPulses();

        double chargeInPulse  = 0.;
        double chargeErrorInPulse  = 0.;
        double totalCharge    = 0.;

        double localPosition[2]  = { 0., 0.};
        double localPositionError[2]  = { 0., 0.};

        Vector2D localPadCenter, globalPadCenter;

        for(vector<TrackerPulse*>::const_iterator aPulse = thePulses.begin(); aPulse < thePulses.end(); ++aPulse)
        {
                globalPadCenter = module.getPadCenter((*aPulse)->getCellID0());
                localPadCenter  = module.globalToLocal(globalPadCenter[0], globalPadCenter[1]);
                chargeInPulse = (*aPulse)->getCharge();
                totalCharge += chargeInPulse;
                for(unsigned int i = 0; i < 2; ++i)
                {
                        localPosition[i] += chargeInPulse * localPadCenter[i];
                }
        }

        Vector2D globalPosition2D =  module.localToGlobal(localPosition[0] / totalCharge, localPosition[1] / totalCharge);
        double xPosition = globalPosition2D[0];
        double yPosition = globalPosition2D[1];

        //  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 = Global::GEAR->getTPCParameters().getMaxDriftLength()
                                                                           -  _getDriftVelocity() * (maxPulse->getTime()) / 1000. ; // time is in [ns], drift velocity in [mm/mus]

        // error calculation :: distinguish two cases!
        // as the variance is the square of the standard deviation, one needs to be careful with the squaring
        covMatrix[5]  =  pow(_getDriftVelocity(),2) * (maxPulse->getCovMatrix()).at(2) / 1000000.; // error propagation from maximum pulse
        double tempDebug1 =0.0;
        double tempDebug2 =0.0;
        if(polarPadPlane)
        {
                for(vector<TrackerPulse*>::const_iterator aPulse = thePulses.begin(); aPulse < thePulses.end(); ++aPulse)
                {
                        globalPadCenter = module.getPadCenter((*aPulse)->getCellID0());
                        localPadCenter  = module.globalToLocal(globalPadCenter[0], globalPadCenter[1]);
                        chargeInPulse = (*aPulse)->getCharge();
                        chargeErrorInPulse = ((*aPulse)->getCovMatrix())[0];
                        streamlog_out(DEBUG4) << "chargeErrorInPulse =" << chargeErrorInPulse <<endl;
                        localPositionError[1] += pow(chargeInPulse*module.getPadWidth((*aPulse)->getCellID0())/totalCharge,2)/12
                                        +pow((localPadCenter[1]*totalCharge - localPosition[1]) / (totalCharge*totalCharge),2)*chargeErrorInPulse;
                        tempDebug1+=pow(chargeInPulse*module.getPadWidth((*aPulse)->getCellID0())/totalCharge,2) / 12;
                        tempDebug2+=pow((localPadCenter[1]*totalCharge - localPosition[1]) / (totalCharge*totalCharge),2)*chargeErrorInPulse;

                }
                //      localPositionError[1]=tempDebug2;
                streamlog_out(DEBUG4) << "totalCharge =" << totalCharge <<endl;
                streamlog_out(DEBUG4) << "localPosition[0]/totalCharge =" << localPosition[0]/totalCharge <<endl;
                streamlog_out(DEBUG4) << "localPosition[1]/totalCharge =" << localPosition[1]/totalCharge <<endl;
                streamlog_out(DEBUG4) << "error from pad width =" << sqrt(tempDebug1) <<endl;
                streamlog_out(DEBUG4) << "error from charge=" << sqrt(tempDebug2) <<endl;
                double meanPhi = localPosition[1] / totalCharge;
                streamlog_out(DEBUG4) << "meanPhi="<< meanPhi<<endl;
                double meanRadius = localPosition[0] / totalCharge;
                double radiusErrorSquare = pow(module.getRowHeight(module.getRowNumber(maxPulse->getCellID0())), 2) / 12.;
                covMatrix[0] = pow(cos(meanPhi), 2) * radiusErrorSquare + pow(meanRadius, 2) * pow(sin(meanPhi), 2) * localPositionError[1];
                covMatrix[1] = sin(meanPhi) * cos(meanPhi) * (radiusErrorSquare - pow(meanRadius, 2) * localPositionError[1]);
                covMatrix[2] = pow(sin(meanPhi), 2) * radiusErrorSquare + pow(meanRadius, 2) * pow(cos(meanPhi), 2) * localPositionError[1];

        }
        else if(cartesianPadPlane)
        {
                for(vector<TrackerPulse*>::const_iterator aPulse = thePulses.begin(); aPulse < thePulses.end(); ++aPulse)
                {
                        globalPadCenter = module.getPadCenter((*aPulse)->getCellID0());
                        localPadCenter  = module.globalToLocal(globalPadCenter[0], globalPadCenter[1]);
                        chargeInPulse = (*aPulse)->getCharge();
                        chargeErrorInPulse = ((*aPulse)->getCovMatrix())[0];
                        localPositionError[0] += pow(chargeInPulse*module.getPadWidth((*aPulse)->getCellID0())/totalCharge,2)/12
                                        +pow((localPadCenter[0]*totalCharge-localPosition[0])/(totalCharge*totalCharge)*chargeErrorInPulse,2);

                }
                covMatrix[0] = localPositionError[0];
                covMatrix[2] = pow(module.getRowHeight(module.getRowNumber(maxPulse->getCellID0())),2) / 12.;
        }

        /*
        double padWidth = module.getPadWidth(maxPulse->getCellID0());
        if(polarPadPlane)
            {
                double radiusErrorSquare = pow(module.getRowHeight(module.getRowNumber(maxPulse->getCellID0())), 2) / 12.;
                double meanPhi = localPosition[1] / totalCharge;
                double meanRadius = localPosition[0] / totalCharge;
                double phiErrorSquare = 0.0;
                double charge = 0.0;
                const double singlePadPhiError = padWidth / sqrt(12.);
                if(cand.getWidth() == 1)
                    phiErrorSquare = pow(singlePadPhiError, 2);
                else
                    {
                        for(vector<TrackerPulse*>::const_iterator aPulse = thePulses.begin(), endPulse = thePulses.end(); aPulse != endPulse; ++aPulse)
                            {
                                charge = (*aPulse)->getCharge();
                                phiErrorSquare += pow((charge * singlePadPhiError) , 2);

                            }
                        phiErrorSquare /= pow(totalCharge, 2);
                    }
                // OK, the diagonal errors have been determined; now do the coordinate transformation, since the geometry can't do it for matrices
                covMatrix[0] = pow(cos(meanPhi), 2) * radiusErrorSquare + pow(meanRadius, 2) * pow(sin(meanPhi), 2) * phiErrorSquare;
                covMatrix[1] = sin(meanPhi) * cos(meanPhi) * (radiusErrorSquare - pow(meanRadius, 2) * phiErrorSquare);
                covMatrix[2] = pow(sin(meanPhi), 2) * radiusErrorSquare + pow(meanRadius, 2) * pow(cos(meanPhi), 2) * phiErrorSquare;

            }
        else if(cartesianPadPlane)
            {
                streamlog_out(WARNING) << "Please make sure that the global x and y are aligned with the local one; in addition: y is the row coordinate, x within row. Otherwise the errors are non-sensical." << endl;
                double xError = 0.0;
                double meanX = localPosition[0] / totalCharge;
                if(cand.getWidth() == 1){
                        xError = pow(padWidth,2) / 12.;
                }

                else
                    {
                        for(vector<TrackerPulse*>::const_iterator aPulse = thePulses.begin(), endPulse = thePulses.end(); aPulse != endPulse; ++aPulse)
                            {
                                xError += pow((localPadPlane.getPadCenter((*aPulse)->getCellID0())[0] - meanX) , 2);
                            }
                    }
                covMatrix[0] = xError;
                covMatrix[2] = pow(module.getRowHeight(module.getRowNumber(maxPulse->getCellID0())),2) / 12.;

            }
         */
        return Vector3D(xPosition, yPosition, zPosition);
}



double RowBasedHitFinderProcessor::_convertCharge(const hitCandidate& cand)
{
        return _getChargeConversionFactor() * (cand.getCharge());
}



RowBasedHitFinderProcessor::hitCandidate::hitCandidate(TrackerPulse* aPulse)
: qualityFlag(0), _maxPulse(NULL), _maxPulseCharge(0.0), _numberOfDeadChannels(0)
{
        _thePulses.clear();
        addPulse(aPulse);
}



void RowBasedHitFinderProcessor::hitCandidate::addPulse(TrackerPulse* aPulse)
{
        _thePulses.push_back(aPulse);
        if(fabs(aPulse->getCharge()) > fabs(_maxPulseCharge))
        {
                _maxPulseCharge = aPulse->getCharge();
                _maxPulse       = aPulse;
        }
}


void RowBasedHitFinderProcessor::hitCandidate::addDeadChannel()
{
        _numberOfDeadChannels++;
}



double RowBasedHitFinderProcessor::hitCandidate::getCharge() const
{
        double tmpCharge = 0.0;
        for(vector<TrackerPulse*>::const_iterator aPulse = _thePulses.begin(); aPulse < _thePulses.end(); ++aPulse)
                tmpCharge += (*aPulse)->getCharge();
        return tmpCharge;
}

double RowBasedHitFinderProcessor::hitCandidate::getChargeMeasError() const
{
        double tmpChargeErrorSqr = 0.0;
        for(vector<TrackerPulse*>::const_iterator aPulse = _thePulses.begin(); aPulse < _thePulses.end(); ++aPulse)
                tmpChargeErrorSqr += ((*aPulse)->getCovMatrix())[0];
        return tmpChargeErrorSqr;
}


double RowBasedHitFinderProcessor::hitCandidate::getMaxPulseHeight()
{
        return _maxPulseCharge;
}



unsigned int RowBasedHitFinderProcessor::hitCandidate::getWidth() const
{
        return _thePulses.size();
}



bool RowBasedHitFinderProcessor::_pulseComparator(TrackerPulse* p0, TrackerPulse* p1)
{
        return fabs(p0->getCharge()) < fabs(p1->getCharge());
}



unsigned int RowBasedHitFinderProcessor::_getMaxDeadChannels()
{
        if(parameterSet("MaxDeadChannelsOverride"))
        {
                return _maxDeadChannelsOverride;
        }
        else
        {
                throw tpcconddata::ConditionsObjectException(string("[RowBasedHitFinderProcessor]::_getDeadChannels()  Exception: number can't be found!"));
        }
}



int RowBasedHitFinderProcessor::_getMinHitSize()
{
        if(parameterSet("MinHitSizeOverride"))
        {
                return _minHitSizeOverride;
        }
        else
        {
                throw tpcconddata::ConditionsObjectException(string("[RowBasedHitFinderProcessor]::_getMinHitSize()  Exception: Number can't be found!"));
        }
}



double RowBasedHitFinderProcessor::_getMinMaxPulseSliceHeight()
{
        if(parameterSet("MinMaxPulseSliceHeightOverride"))
        {
                return _minMaxPulseSliceHeightOverride;
        }
        else
        {
                throw tpcconddata::ConditionsObjectException(string("[RowBasedHitFinderProcessor]::_getMinMaxPulseSliceHeight()  Exception: Number can't be found!"));
                return 0.;
        }
}



double RowBasedHitFinderProcessor::_getMaxTimeSpread()
{
        if(parameterSet("MaxTimeSpreadOverride"))
        {
                return _maxTimeSpreadOverride;
        }
        else
        {
                throw tpcconddata::ConditionsObjectException(string("[RowBasedHitFinderProcessor]::_getMaxTimeSpread()  Exception: Number can't be found!"));
                return 0.;
        }
}



double RowBasedHitFinderProcessor::_getDriftVelocity()
{
        if(parameterSet("VDriftOverride"))
        {
                return _driftVelocityOverride;
        }
        else
        {
                try
                {
                        return _tpcConditionsListener->getDriftVelocity();
                }
                catch(...)
                {
                        streamlog_out(ERROR) << "[getDriftVelocity failed] The Listener for TPCConditions could not be used. Maybe the conditions object is not available?" << endl;
                        streamlog_out(ERROR) <<  "                          Possibly you need to use the override parameterin the steering file called (\"VDriftOverride\")." << endl;
                        abort();
                }
        }
}



bool RowBasedHitFinderProcessor::_padIsDead(const int padID, const int moduleID)
{
        if(parameterSet("DeadPadsOverride"))  // this is just for a single module!
        {
                return (_deadChannelIndicesOverride.end() != find(_deadChannelIndicesOverride.begin(), _deadChannelIndicesOverride.end(), padID));
        }
                else
        {
                pair<int, int> tmpPair;
                tmpPair.first  = padID;    // CellID0
                tmpPair.second = moduleID; // CellID1
                bool returnvalue=false;
                try
                {
                        if (_channelCorrectionListener->isDeadFromGearID(tmpPair)||_channelCorrectionListener->isNoisyFromGearID(tmpPair))
                        //if (_channelCorrectionListener->isDeadFromGearID(tmpPair))
                        {
                        returnvalue=true;
                        }


                }
                catch (CondDataNotAvailableException &e)
                {
                }
                return returnvalue;
        }
}



bool RowBasedHitFinderProcessor::_padIsNoisy(const int moduleID, const int padID)
{

        if(parameterSet("DeadPadsOverride"))  // this is just for a single module!
                return (_deadChannelIndicesOverride.end() != find(_deadChannelIndicesOverride.begin(), _deadChannelIndicesOverride.end(), padID));
        else
        {
                pair<int, int> tmpPair;
                tmpPair.first  = padID;    // CellID0
                tmpPair.second = moduleID; // CellID1
                bool returnvalue;
                try
                {
                        returnvalue=_channelCorrectionListener->isNoisyFromGearID(tmpPair);
                        streamlog_out(MESSAGE)<<"found noisy pad"<<endl;
                }
                catch (CondDataNotAvailableException &e)
                {
                        returnvalue=false;
                }
                return returnvalue;
        }

}



double RowBasedHitFinderProcessor::_getChargeConversionFactor()
{
        if(parameterSet("ChargeConversionFactorOverride"))
        {
                return _chargeConversionFactorOverride;
        }
        else
        {
                return _tpcConditionsListener->getADCtoPrimaryElectronsFactor();
        }
}
}// namespace marlintpc