/*
 * CauchyDistribution.cpp
 *
 *  Created on: Oct 13, 2015
 *      Author: Hovanes Egiyan
 */

#include "CauchyDistribution.hh"

using namespace ROOT::Minuit2;
using namespace std;

const unsigned CauchyDistribution::cdNumberOfParameters = 6;

CauchyDistribution::CauchyDistribution( const std::vector<double>* xAxis,
		const std::vector<double>* yAxis ) :
		FCNGradientBase(), cdErrorDef( 1.0 ), cdAxisX( xAxis ), cdAxisY( yAxis ) {
	return;
}

CauchyDistribution::~CauchyDistribution() {
	return;
}

CauchyDistribution& CauchyDistribution::operator =( const CauchyDistribution& dist ) {
	FCNGradientBase::operator=( dist );
	cdErrorDef = dist.cdErrorDef;
//	cdAxisX = new vector<double>( *dist.cdAxisX );
//	cdAxisY = new vector<double>( *dist.cdAxisY );
	// Need to copy pointers since this class is not supopeed to have the control
	// over the axis values.
	cdAxisX = dist.cdAxisX;
	cdAxisY = dist.cdAxisY;

	return *this;
}

// Calculates the chi2 for the Cauchy distribution .
// The function is of the form alpha/(1+((x_i-x0)/gamma)^2) ,
// where alpha is the amplitude in the median point, gamma is the
// full width, and the x0 is the median. The standard deviation is
// assumed to be follow poison statistics, e.g. sigma^2 = y .
double CauchyDistribution::operator()( const vector<double>& x ) const {
	double alpha 	= x[0]; // amplitude
	double x0 		= x[1]; // median
	double gamma 	= x[2]/2.0; // half width, second parameter for the function is the full width
	double bkg		= x[3]; // flat background
	double xMin     = x[4];
	double xMax     = x[5];

	double chi2 = 0; // running sum starts at 0
	double x_i, y_i, sigma2, mainPart, funVal;
	for ( unsigned iPoint = 0; iPoint < (*cdAxisX).size(); iPoint++ ) {
		x_i = (*cdAxisX)[iPoint];
		if( xMin < x_i && x_i < xMax ) {
		  y_i = (*cdAxisY)[iPoint];
                  if( y_i > 1.0e-3 ) {
                    sigma2 = fabs( y_i );
                    mainPart = 1.0 / (1.0 + pow( (x_i - x0) / gamma, 2.0 ));
                    funVal = alpha * mainPart + bkg;
                    chi2 += ( funVal - y_i ) * ( funVal - y_i ) / sigma2;
                  }
		}
	}
	return ( chi2 / (*cdAxisX).size() );
}

// Calculates the gradient of the chi2 with respect to the parameters for Cauchy distribution.
// The function is of the form alpha/(1+((x_i-x0)/gamma)^2) ,
// where alpha is the amplitude in the median point, gamma is the
// full width, and the x0 is the median. The standard deviation is
// assumed to be follow poison statistics, e.g. sigma^2 = y .
vector<double> CauchyDistribution::Gradient( const std::vector<double>& x ) const {
	double alpha 	= x[0]; // amplitude
	double x0 		= x[1]; // median
	double gamma 	= x[2]/2.0; // full width, second parameter for the function is the full width
	double bkg		= x[3]; // flat background
	double xMin     = x[4];
	double xMax     = x[5];

	double dChi2_dAlpha = 0; // running sum starts at 0
	double dChi2_dGamma = 0; // running sum starts at 0
	double dChi2_dX0 = 0; // running sum starts at 0
	double dChi2_dBkg = 0; // running sum starts at 0

	double x_i, y_i, sigma2, mainPart;
	for ( unsigned iPoint = 0; iPoint < cdAxisX->size(); iPoint++ ) {
		x_i = (*cdAxisX)[iPoint];
		if( xMin < x_i && x_i < xMax ) {
		  y_i = (*cdAxisY)[iPoint];
		  sigma2 = fabs( y_i );
		  mainPart = 1.0 / (1.0 + pow( (x_i - x0) / gamma, 2.0 ));
		  // Derivative with respect to amplitude
		  dChi2_dAlpha += 2.0 * mainPart * (alpha * mainPart + bkg - y_i);
		  //  Derivative with respect to the half-width
		  dChi2_dGamma += ( 4.0 * alpha / (sigma2 * gamma*gamma*gamma) ) * (alpha * mainPart + bkg - y_i) *
				mainPart * mainPart * pow( (x_i - x0) / gamma, 2 );
		  // Derivative with respect to mean value
		  dChi2_dX0 += ( 4.0 * alpha / (sigma2 * gamma*gamma) ) * (alpha * mainPart + bkg - y_i) *
				( mainPart * mainPart ) * (x_i - x0);
		  // derivative with respect to flat background
		  dChi2_dBkg += ( 2.0 / sigma2 ) * (alpha * mainPart + bkg - y_i) * mainPart;
		}

	}
	// Create and fill the vector of derivatives
	vector<double> vectorOfDerivatives( cdNumberOfParameters );
	vectorOfDerivatives[0] = dChi2_dAlpha;
	vectorOfDerivatives[1] = dChi2_dX0;
	vectorOfDerivatives[2] = dChi2_dGamma / 2.0; // Derivative with respect two full width
	vectorOfDerivatives[3] = dChi2_dBkg;
	vectorOfDerivatives[4] = 0;
	vectorOfDerivatives[5] = 0;

	return vectorOfDerivatives;
}