/*
 * ProfilerPlaneSliceFitter.cpp
 *
 *  Created on: Oct 17, 2015
 *      Author: hovanes
 */

#include "ProfilerPlaneSliceFitter.hh"

using namespace ROOT::Minuit2;
using namespace std;

double ProfilerPlaneSliceFitter::ppsfRange = 15.0 ; // range for fitting in mm
double ProfilerPlaneSliceFitter::ppsfMaxWidthGuess = 40; 	// maximum width for initial width parameter

// Main constructor
ProfilerPlaneSliceFitter::ProfilerPlaneSliceFitter( const std::string name,
		const std::vector<double>& xAxis, const std::vector<double>& yAxis ) :
		MutexedClass(), ProfilerPlaneSlice( name, xAxis, yAxis ), ProfilerPlaneFitter( name ),
				ppsfDistFunc( 0 ) {
	ppsfDistFunc = new CauchyDistribution( &ppsAxis, &ppsValues );
	gErrorIgnoreLevel = 1001;
}

// Copy constructor
ProfilerPlaneSliceFitter::ProfilerPlaneSliceFitter( const ProfilerPlaneSliceFitter& f ) :
		MutexedClass( f ), ProfilerPlaneSlice( f ), ProfilerPlaneFitter( f ),
				ppsfDistFunc( f.ppsfDistFunc ) {
	return;
}

// Destructor
ProfilerPlaneSliceFitter::~ProfilerPlaneSliceFitter() {
	// Delete the distribution function created in the constructor
	if ( ppsfDistFunc != 0 ) {
		delete ppsfDistFunc;
		ppsfDistFunc = 0;
	}
}

// Assignment operator
ProfilerPlaneSliceFitter& ProfilerPlaneSliceFitter::operator=( const ProfilerPlaneSliceFitter& f ) {
	MutexedClass::operator=( f );
	ProfilerPlaneSlice::operator=( f );
	ProfilerPlaneFitter::operator =( f );
	MtxLock objLock( mcMutex );
	ppsfDistFunc = f.ppsfDistFunc;
	return *this;
}

// Check if the fit was good and return the parameter state reference. If it was not, redo the fit.
// This method was created because UserState return const reference whic cannot be reassigned.
const MnUserParameterState& ProfilerPlaneSliceFitter::checkFit( FunctionMinimum& fitMin,
		MnUserParameters&  initialValues ) {
//	if( ! fitMin.IsValid() ) return false;
	// Get the parameter object from which we can pull out the values.
	const MnUserParameterState& minStateInit = fitMin.UserState();
	double chi2df = minStateInit.Fval();
	if( chi2df > 3 ) {
		// If Chi2 is not so good redo the fit using the new parameters
		// estimates as the starting point. Also Change the error definition
		// to chi2.
		ppsfDistFunc->setErrorDef( chi2df );
		initialValues.SetValue("Amplitude", minStateInit.Value( "Amplitude" ) );
		initialValues.SetValue("Mean", minStateInit.Value( "Mean" ) );
		initialValues.SetValue("Width", minStateInit.Value( "Width" ) );
		initialValues.SetValue("Background", minStateInit.Value( "Background" ) );

		fitMin = this->Minimize( *ppsfDistFunc, initialValues );
		return fitMin.UserState();
	} else {
		ppsfDistFunc->setErrorDef( 1.0 );
		return minStateInit;
	}
}

// Fit this slice with Cauchy distribution
bool ProfilerPlaneSliceFitter::fitWithCauchy() {
	// Set the initial parameters
	MnUserParameters initialValues = this->findInitialParameters();
	this->setStrategy( MnStrategy(2) );
	// Try to fit the minimum using Minuit
	ROOT::Minuit2::FunctionMinimum fitMin = this->Minimize( *ppsfDistFunc, initialValues );

	const MnUserParameterState& minState = this->checkFit( fitMin, initialValues );

	ppsAmpl = minState.Value( "Amplitude" );
	ppsMean = minState.Value( "Mean" );
	ppsWidth = minState.Value( "Width" );
	ppsBkg = minState.Value( "Background" );
	ppsSignal = ppsAmpl * ppsWidth * 2.0 * acos( 0 );

//	cout << "Fit Chi2/N for " << ProfilerPlaneSlice::getName() << " is " << minState.Fval() << endl;

	ppsAmplErr = minState.Error( "Amplitude" );
	ppsMeanErr = minState.Error( "Mean" );
	ppsWidthErr = minState.Error( "Width" );
	ppsBkgErr = minState.Error( "Background" );
	if ( ppsAmpl > 0 && ppsWidth > 0 ) {
		ppsSignalErr = ppsSignal
				* ((ppsAmplErr / ppsAmpl) * (ppsAmplErr / ppsAmpl)
						+ (ppsWidthErr / ppsWidth) * (ppsWidthErr / ppsWidth));
	} else {
		ppsSignalErr = 0;
	}

//	cout << "Resulting amplitude for " << ProfilerPlaneSlice::getName() << " is " <<
//			ppsAmpl << " , mean is " << ppsMean << " , width is " << ppsWidth << endl;

	return fitMin.IsValid();
//	return true;

}

// This method tries to guess the initial parameters for the fit function.
// The width is twice the RMS, the amplitude is the largest value, and the
// median is the mean values of the x-distribution using y-values as weights.
ROOT::Minuit2::MnUserParameters ProfilerPlaneSliceFitter::findInitialParameters() {
	// Calculate the sum of all Y SUM f_i
	double norm = std::accumulate( ppsValues.begin(), ppsValues.end(), 0.0 );
	double mean = 0;
	double rms = 0;

//	cout << "Normalization sum is " << norm << endl;

	// Find the mean and standard deviation
	if ( norm > 1.0e-20 ) {
		// multiply elements of the x-vector by the values of the y-vector to get x_i x y_i vector
		vector<double> singleProduct( ppsAxis.size() );
		std::transform( ppsAxis.begin(), ppsAxis.end(), ppsValues.begin(), singleProduct.begin(),
				std::multiplies<double>() );
		// the sum of the new vector divided by the sum of all weights (y-s) is the mean
		mean = std::accumulate( singleProduct.begin(), singleProduct.end(), 0.0 ) / norm;

		// Find sigma^2 as Average(X^2) - Average(X)^2
		vector<double> doubleProduct( ppsAxis.size() );
		// multiply the x_i x y_i by another x_i to get y_i x x_i^2
		std::transform( singleProduct.begin(), singleProduct.end(), ppsAxis.begin(),
				doubleProduct.begin(), std::multiplies<double>() );
		// RMS (standard deviation) is the square root of difference between average x^2 and average x squared
		rms = sqrt(
				std::accumulate( doubleProduct.begin(), doubleProduct.end(), 0.0 ) / norm
						- mean * mean );
	}
	// Full width is approximately twice the sigma divided by two since the variance for cauchy distribution is
	// not defined. This is very rough.
	double width = 2.0 * rms / 4.0 ;
	width = std::min( ppsfMaxWidthGuess, width );
	// Estimate the amplitude as the largest number in the array of values.
	// max_element method returns the iterator for the largest element of the vector, needs dereferencing.
	double ampl = *std::max_element( ppsValues.begin(), ppsValues.end() );
	double bkg = 0.01 * ampl;

	double relErr = 0.3;
	MnUserParameters pars;
	pars.Add( "Amplitude", ampl, relErr * ampl );
	pars.Add( "Mean", mean, relErr * mean );
	pars.Add( "Width", width, relErr * width );
	pars.Add( "Background", bkg, relErr * bkg );

	pars.Add( "FitRangeLow", mean - ProfilerPlaneSliceFitter::ppsfRange, 0.1 );
	pars.Fix( "FitRangeLow" );
	pars.Add( "FitRangeHigh", mean + ProfilerPlaneSliceFitter::ppsfRange, 0.1 );
	pars.Fix( "FitRangeHigh" );

//	pars.SetLowerLimit( "Amplitude", 0 );
//	pars.SetLowerLimit( "Width", 0 );
//	pars.SetLowerLimit( "Background", 0 );

//	cout << "Initial amplitude is " << pars.Parameter(0).Value() << " , mean " << pars.Parameter(1).Value() <<
//			" , width is " << pars.Parameter(2).Value() << endl;

	return pars;
}