/*!
* @file STools.h
* @author Nick Ryder, on behalf of the SoLid collaboration.
* @date 18 Nov 2016
*/

#ifndef __STOOLS_H__
#define __STOOLS_H__

#include <algorithm>
#include <functional>
#include <vector>

#include <TH1D.h>

//! Time sort any vector objects that have a htime().
template<class obj> void timeSortReconObjects(std::vector<obj>* unsortedObjects) {
    std::sort(unsortedObjects->begin(), 
          unsortedObjects->end(), 
          [](obj a, obj b)->bool{return a->htime() < b->htime();}
    );
}

//! pedestall finder, takes a samples histogram, returns the pedestall value
float PedestalFinder(TH1D *hin);

// functions related with gain finding
//==============================================================================
typedef std::vector<float> FloatVector;
typedef std::function<int(TH1D*, FloatVector&, FloatVector&)> PeakFindFunction;
/*!
 * \brief Finds the gain using the residuals method:
 * Peaks are found in the input histogram using f,
 * A scan is then done over a range of gain candidates:
 * all distances between the peaks are calculated and it is checked
 * how close all these distances are to an integere multiple
 * of the gain candidate. The candidate that is closest to integer
 * is considered the gain.
 * \param h  histogram of which the gain will be determined
 * \param f  peakfinding function to be applied on the histogram
 * \return the calulated gain, returns -1 if no gain was found,
 * \return returns -2 if not enough statistics were there
 */
float residualsMethodGainFinder(TH1D *h, PeakFindFunction f);

/*!
 * \brief finds the range around a peak in which a gaussian fit can be done,
 * starts from the peak position and moves out as long as the first
 * derivative is positive (left) / negative (right)
 * \param[in]  hin  histogram in which a fit range has to be found
 * \param[in]  pp   peak position in the histogram
 * \param[out] rlb  fit range low bound
 * \param[out] rhb  fit range high bound
 */
void findFitRange(TH1D *hin, float pp, double &rlb, double &rhb);

/*!
 * \brief calculates the 5 point first derivative if possible,
 * lower (3 or 2 points) if needed, returns 0 if out of range
 * \param h histogram in which the derivative is calculated
 * \param b bin in which the derivative has to ba calculated
 * \return the first derivative (or 0 if out of range)
 */
float FirstDerivative(TH1D *h, int b);

/*!
 * \brief calculates the residuals for a point
 * \param v vector of peak positions
 * \param t test point for which the residuals are calulated
 * \return residuals of the point t in a peak configuration listed in v
 */
float CalculateResiduals(std::vector<float> v, float t);

/*!
 * \brief peak finder for gain calculation,
 * uses a simple local maximum (5 bins on either side),
 * returns peaks sorted with descending amplitude
 * \param[in]  hin  histogram for which the peaks have to be found
 * \param[out] posX vector with the X positions of the peaks
 * \param[out] posY vector with the peak values
 * \return the number of peaks found
 */
int findGainPeaks_LocMax(TH1D *hin, std::vector<float> &posX,
                            std::vector<float> &posY);

//! Simple gain finding.
/*! This gain finder is a simple and robust method for finding gains, although
 * the precision is yet to be understood. The method involved 'folding' PA peaks
 * from peak amplitude (or integral) spectra back on eachother, and does this
 * for various 'trial gains'. The method postulates that at the correct gain,
 * the width of this folded distribution (taken as the RMS to avoid fitting) is
 * minimum when the trial gain equals the actaul gain.
 *
 * This has been tested for a small number of files from both Bristol and Gent
 * Note: this method is not to replace more sophisticated gain finding methods -
 * rather to be used as a cross check. It is expected that this method gives
 * less precise results, but is more robust against problems from peak finding,
 * gaussian fitting and tuning parameters.
 *
 * It seems to work for around 2000 random trigger waveforms. The only parameters
 * that need tuning is the search range, and it is already known that outside this
 * range strange things may happen.
 *
 * \param h   histogram for which the gain has to be found
 * \param amp boolean to state if the histogram is an amplitude or an integral
 *            histogram, influences the scanned range
 * \return the calculated gain
 */
double histoFoldingGainFinder(TH1D * h, bool amp);

#endif