MMScatterProcessor.h

#ifndef MMSCATTERPROCESSOR_H 
#define MMSCATTERPROCESSOR_H 

#include <marlin/Processor.h>
#include <lcio.h>
#include <string>
#include <map>
#include <vector>

#include "MMPadResponseFunction.h"

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

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

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

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

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

//ROOT
#include <Minuit2/MnUserParameters.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/Minuit2Minimizer.h>

#include "Math/Minimizer.h"
#include "Math/Functor.h"
#include <TMath.h>
#include <TTree.h>
#include <TH1.h>
#include <TF1.h>



namespace marlintpc
{
  class MMScatterProcessor : public marlin::Processor
  {
    public:

      //Mid-calibration parameters for use by the MMHitFinderProcessor
      static std::vector<double> CALIB_PARAMS;

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

      MMScatterProcessor();

      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 std::pair<double, double> calculateTrackPosition   (       const gear::TPCParameters&, const IMPL::TrackImpl*, int mod, int row                                                            );
      static std::pair<double, double> calculateTrackPosition   (       const gear::TPCParameters&, double d0, double phi0, double omega, const float * ref, int mod, int row   );
      static std::pair<double, double> calculateStraightTrackPosition ( const gear::TPCParameters&, double d0, double phi0, const float * ref, int mod, int row );

    protected:
      typedef std::map<int, std::map<int, std::vector<TrackerPulse*> > > SortedPulseMap;
      typedef std::map<std::pair<int, int>, std::vector<std::pair<double, float> > > PointsMap;

      lcio::StringVec    _inputColNames; // Name of the input collection
      bool                                                       _fitPRF; //Choose whether or not to fit the PRF to the scatter-plot
      bool                                                       _rowDep; //Choose whether or not to store the scatter-plot for each row
      bool                                                       _moduleDep; //Choose whether or not to store the scatter-plot for each row
      bool                                                       _singleTrack;
      bool                                                       _smallAngle;
      bool                                                       _profile;
      int                                                        _curEvt; //Keep track of the current event.
      int                _normType;
      int                _nPar;
      double *                                   _parArray;
      std::string                                _output_path;
      std::string                                _prfType;
      std::vector<float> _parVec;
      std::vector<int>   _moduleList;
      std::map<int, int> _moduleMap;
      PointsMap                                  _scatterMap;

      std::vector<std::vector<double> > _iterationPar; //Store the fit parameters for each run
      std::vector<double> _iterationNumbers; //Store the iteration numbers 
      int _nEvt;
      int _nIterations;
      int _calibIter;
      bool _calibRun;
      bool _calibPad;

      //A Histogram to hold the collective scatterplot.
      std::map<std::pair<int,int>, TH1*>                        _scatterPlots;

      //Histograms and a tuple to hold the result of fitting the PRF
      AIDA::ITuple*      fitData;

      /************************************************************************************************/
      // 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   (      const EVENT::TrackerHitVec&, const gear::TPCParameters&, SortedPulseMap& );

      // Fills the data points
      void _fillData             (      const SortedPulseMap&, const gear::TPCParameters&, const IMPL::TrackImpl*, PointsMap& );

      void calibrate();

      //Make a small, private class to use to fit the PRF to the scatter plots
      class PRFFcn
      {

        public:

          PRFFcn(std::string prfType, int nPar)
          {
            _prf = MMPadResponseFunction(prfType);
            _nPar = nPar;
          }

          double operator()(const double * x, const double * par) const
          {

            return par[_nPar]*_prf.prf(x[0] - par[_nPar+1], par);
          }

        private:
          MMPadResponseFunction _prf;
          int _nPar; 
      };

  };//End MMScatterProcessor class

}//End namespace marlintpc

#endif // MMSCATTERPROCESSOR_H