MMHitFinderProcessor.h

#ifndef MMHITFINDERPROCESSOR_H
#define MMHITFINDERPROCESSOR_H

// MARLIN
#include "marlin/Processor.h"
#include "marlin/Global.h"

// LCIO
#include "lcio.h"
#include "EVENT/LCCollection.h"
#include "IMPL/LCCollectionVec.h"
#include "EVENT/TrackerPulse.h"

// GEAR
#include "gear/TPCParameters.h"

// C++
#include <string>
#include <vector>
#include <map>
#include <set>
#include <deque>

// MarlinTPC
#include "ChannelCorrectionListener.h"
#include "TPCConditionsListener.h"

// Aida
#include <AIDA/AIDA.h>
#include <marlin/AIDAProcessor.h>
#include <AIDA/IHistogramFactory.h>
#include <AIDA/IHistogram1D.h>

// ROOT
#include <Minuit2/FCNBase.h>
#include <TMath.h>

// Minuit2 headers used to find seed track.
#include <Minuit2/MnUserParameters.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/Minuit2Minimizer.h>

#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"

namespace marlintpc
{

  /*  Processor to find Hits from Micromegas data.
   *
   *  Basic procedure: Sort pulses by module and row, search rows for groups of associated pulses, 
   *  then either fit pulses with the Pad Response Function, or use the weighted mean if calibrating.
   *
   *  Parameters:
   *
   *    InputTrackerPulses      =  Input pulses  [string]
   *    OutputTrackerHits       =  Output hits  [string]
   *    OutputTrackerHitPulses  =  Output collection of pulses associated with output hits  [string]
   *    WriteOutputToStorage    =  Set persistency of output collections  [string]
   *    MaxEmptyPads            =  Maximum number of empty pads in a hit  [integer]
   *    MinHitSizeOverride      =  Minimum number of pulses in a hit  [integer]
   *    VDriftOverride          =  Drift velocity (not really important at this stage, but will be later)  [float]
   *    ModuleList              =  A vector of modules to be used in the analysis  [IntVec]
   *    MaxTimeSpreadOverride   =  Maximum time between pulses in a hit  [float]
   *    PRFType                 =  Name of the PRF to be used  [string]
   *    PRF_Parameters          =  PRF parameters to be used (can be empty only if it is a calibration run)  [FloatVec]
   *    RowDep                  =  Set whether or not the PRF parameters have a row dependency  [boolean]
   *    DeadRowList             =  A list of rows not to be used in the analysis. Format: MODULE ROW (e.g.  0 1 0 2 ... etc.)  [IntVec]
   *    CentreColumns           =  A list of central beam columns (same format as for the rows). Only data in columns +/- 5 of the central column is used  [IntVec]
   *    Threshold               =  Minimum amplitude of the maximum pulse in a hit  [integer]
   *    CalibRun                =  IMPORTANT! Sets whether or not this is a calibration run. Uses the weighted mean instead of the PRF, and related settings  [boolean]
   *    ExcludeRows             =  Whether or not to exclude rows in DeadRowList  [boolean]
   *    SelectCentralCols       =  Whether or not to use the central 10 columns indicated in CentreColumns  [boolean]
   *    ModuleDependence        =  Set wheter or not the PRF parameters have a module dependency  [boolean]
   *    PadHits                 =  Whether or not to pad hits with empty pulses to ensure a minimum number of pulses per hit  [boolean]
   *
   *  Author: P. Hayman, Carleton University
   */

  class MMHitFinderProcessor : public marlin::Processor
  {

    public:

      virtual Processor*  newProcessor(){
        return new MMHitFinderProcessor;
      }

      MMHitFinderProcessor() ;

      virtual void init() ;
      virtual void processRunHeader(EVENT::LCRunHeader* run) ;
      virtual void processEvent(EVENT::LCEvent * evt) ;
      virtual void check(EVENT::LCEvent * evt) ;
      virtual void end() ;


      //Static function to calculate the row dependency of PRF parameters, if desired. This function can be used from anywhere that needs this information.
      static double calcRowDep(const double * x, int index, int row) 
      {
        //Currently a quadratic form: par[i] = a[i]*row^2 + b[i]*row + c[i]

        return (double)(row*row)*x[3*index] + (double)(row)*x[3*index + 1] + x[3*index + 2];
      }

      //Also accessible from anywhere, this function will fill an array compatible with the MMPadResponseFunction, taking into account both row and module 
      //dependencies, if desired.
      //static const double * calcParArray(const double * x, int nPar, int nMod, int modIndex, int rowIndex)
      static const std::vector<double> calcParArray(const double * x, int nPar, int nMod, int modIndex, int rowIndex)
      {
        //double * retVal = new double[nPar];
        std::vector<double> retVal(nPar);
        if(modIndex >= 0)
        {
          if(rowIndex >= 0)
          {
            for(int i = 0; i < nPar; i++)
            {
              retVal[i] = calcRowDep(x, i*nMod + modIndex, rowIndex);
            }
          }
          else
          {
            for(int i = 0; i < nPar; i++)
            {
              retVal[i] = x[i*nMod + modIndex];
            }
          }
        }
        else if(rowIndex >= 0)
        {
          for(int i = 0; i < nPar; i++)
          {
            retVal[i] = calcRowDep(x, i, rowIndex);
          }
        }
        else
        {
          //return x;
          for(int i = 0; i < nPar; i++)
            retVal[i] = x[i];
        }

        return retVal;
      }

      // a helper class 
      class HitCandidate
      {
        public:

          //Assume the first pulse is the largest. This could be made more dynamic, but hopefully that's not necessary.
          HitCandidate(EVENT::TrackerPulse* max)
          {
            _thePulses.push_back(max);  
            _maxPulse = max;
            qualityFlag = 0;
          };

          void addPulse(const EVENT::TrackerPulse* aPulse, bool left)
          {
            if(left)
              _thePulses.push_front(aPulse);    
            else
              _thePulses.push_back(aPulse);
          };

          const EVENT::TrackerPulse* getMaximumPulse() const{
            return _maxPulse;
          };

          double getMaxCharge() const {  return _maxPulse->getCharge();  };

          int getNumberOfPulses() const {  return _thePulses.size();  };        

          const std::deque<const EVENT::TrackerPulse*>& getPulses() const{
            return _thePulses;
          };

          int getQuality() const {  return qualityFlag; };
          void setQuality(int q) {      qualityFlag = q;        };

        private:

          int qualityFlag;
          std::deque<const EVENT::TrackerPulse*> _thePulses;    //Use a double-ended queue. First, the max pulse is added. 
                                                              //Then, the pulses to the left of it are added to the front, 
                                                              //and the pulses to the right of it are added to the end.
          const EVENT::TrackerPulse* _maxPulse;
      };

    private:

      //Processor Parameters
      std::string           _inputTrackerPulsesCollectionName ;
      std::string           _outputTrackerHitsCollectionName ;
      std::string           _outputTrackerHitPulsesCollectionName ;
      std::string           _inputTPCConditionsCollectionName;

      // A vector containing the modules which should be used in the analysis
      std::vector<int>      _moduleListParam;

      bool                  _outputHitsPersistent;
      bool                  _rowDep;

      int                   _maxEmptyChannelsOverride;
      int                   _maxDeadChannelsOverride;
      int                   _minHitSizeOverride;
      float                 _minMaxPulseSliceHeightOverride;
      float                 _maxTimeSpreadOverride;
      float                 _driftVelocity;

      std::string           _rootFileName;
      std::string           _prfType;

      std::vector<float>    _parVec;


      //Optional Parameters
      std::vector<int>                  _centreColListParam;
      std::vector<int>                  _deadRowListParam;

      int                   _threshold;

      bool                                                              _excludeRows;
      bool                                                              _selectCentralCols;
      bool                                                              _calibrationRun;
      bool                                                              _moduleDep;
      bool                                                              _padHits;


      //Class Variables
      int                   _evt;
      int                         _nPar;

      double *                                            _parArray;

      std::set<int>         _moduleList;
      std::map<int, std::set<int> >             _deadRowList;
      std::map<int, int >         _centreColList;
      std::map<int, int >         _moduleMap;

      ROOT::Minuit2::Minuit2Minimizer* min;

      // Some AIDA histograms used for error checking.
      std::map<int, AIDA::IHistogram1D *>  numberOfPulsesPerHit;
      std::map<int, AIDA::IHistogram1D *>  averagePulsesPerHitPerRow;
      std::map<int, AIDA::IHistogram1D *>  numberOfHitsPerEvent;
      std::map<int, AIDA::IHistogram1D *>  chiSquareDistrubtionOfHits;
      std::map<int, AIDA::IHistogram1D *>  hitFitResults;
      std::map<int, AIDA::IHistogram1D *>  hitDistribution;
      std::map<int, AIDA::IHistogram1D *>  hitErrorDistribution;
      std::map<int, AIDA::IHistogram1D *>  hitShapes;

      //Vectors used for calculating the average number of pulse in a hit per row.
      std::map<int, std::vector<int> >    numberOfHitsInRow;
      std::map<int, std::vector<int> >    numberOfPulsesInRow;


      //Typedefs 

      //Can't declare this before HitCandidate, for some reason...
      typedef std::map<int, std::vector<MMHitFinderProcessor::HitCandidate> > SortedHitMap;
      typedef std::map<int, std::map<int, std::vector<TrackerPulse*> > > SortedPulseMap;


      /************************************************************************************************/
      // helper functions:

      // Sorts pulses into collection corresponding to modules and rows. 
      // i.e a pulse on row 5 of module 1 is put in the same container as the other pulses
      // from that row and module.
      void _sortPulses                                                                  (       EVENT::LCCollection*, const gear::TPCParameters&, SortedPulseMap&               );

      // Collects the pulses in a row into a potential hit. This is the first stage of hit determination.
      // As long as the pulses are beside each other spatially and close in time they are a potential hit.
      // Hence the name hit CANDIDATE.
      void _findHitCandidates                                                           (       SortedPulseMap&, const gear::TPCParameters&, SortedHitMap&                              );

      // Does all the neseccary checks and calculations to turn a HitCandidate into a full fledged hit!
      // It then stores this TPCHit in the output hit collection which is saved to the lcio file.
      void _createHitCollection                                                 (       SortedHitMap&, const gear::TPCParameters&, IMPL::LCCollectionVec*               );

      // All these functions, as the names suggest, calculate the position and position error of the hit in 3D.
      // The "fixed" coordinate is determined by the row position.
      // The "free" coordinate is determined by fitting the PRF to the pulse amplitudes in the row.
      // The "Z" coordinate is determined from the time measurement of the main pulse in the hit.
      // Each has a return value of a pair: first = value, second = error.
      std::pair<double, double> _calculateFixedCoordinate       (       const MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters&   );
      std::pair<double, double> _calculateFreeCoordinate        (       const MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters&   );
      std::pair<double, double> _calculateFreeCoordinatePRF     (       const MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters&   );
      std::pair<double, double> _calculateFreeCoordinateWeightedMean    (       const MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters&   );
      std::pair<double, double> _calculateZCoordinate           (       const MMHitFinderProcessor::HitCandidate&                                                               );

      //Useful helper method for finding the largest pulse in a row
      EVENT::TrackerPulse* _findHighestPulseInRow                       (       std::vector<EVENT::TrackerPulse*>&                                                                              );
      //Useful helper method for finding pulses that are near a max pulse (i.e., forming a HitCandidate)
      void _addNeighbourPulses                                                  (       MMHitFinderProcessor::HitCandidate&, const gear::TPCParameters&, 
          std::vector<EVENT::TrackerPulse*>&                                                                            );

  };//End MMHitFinderProcessor class

}//End namespace marlintpc

#endif