/*
 * PndHypGePeakFitFunction.cxx
 * 
 * Copyright 2013 Marcell Steinen <steinen@kph.uni-mainz.de>
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 * MA 02110-1301, USA.
 * 
 * 
 */

	//support class for peak fitting, ROOT needs this external construction (?)
#include "PndHypGePeakFitFunction.h"


PndHypGePeakFitFunction::PndHypGePeakFitFunction()
{
	
}

Double_t PndHypGePeakFitFunction::PeakFuncGaussian (Double_t *x,Double_t *par)
{
			// Peak shape is taken from gf2 http://www.phy.anl.gov/gammasphere/doc/gf2.hlp
			// composed by the sum of:
			//	(1) a Gaussian,
			//	(2) a smoothed step function to increase the background on the low-energy
			//			side of the peak. Components (2) and/or (3) can easily be set to zero if not
			//			required, and
			//	(3)	a constant				// not in gf2
			//

			//	(2)	y = constant * ERFC( (x-c)/(SQRT(2)*sigma) )
			
			// par[0]: amplitude gaus, par[1]: x0; par[2]: sigma, par[3]: amplitude step, par[4]: gradient of lin bg
			
			// with sigma in amplitude
			//Double_t Gausian = par[0]*TMath::Sqrt(2 * TMath::Pi())/par[2] * TMath::Exp(- 0.5*pow((x[0] - par[1])/par[2],2)); // par[0]: amplitude, par[1]: x[0]0; par[2]: sigma
			// or without sigma in amplitude
			
			//Double_t Gausian = par[0]/*TMath::Sqrt(2 * TMath::Pi())/par[2]*/ * TMath::Exp(- 0.5*pow((x[0] - par[1])/par[2],2)); // par[0]: amplitude, par[1]: x[0]0; par[2]: sigma
			//Double_t SmoothedStep = par[3] * TMath::Erfc( (x[0] - par[1]) /(TMath::Sqrt(2)*par[2])); // par[3]: amplitude
			//Double_t Constant = par[4]*x[0]; // par[4]: gradient of lin b
	//return Gausian + SmoothedStep +  Constant;
	return GausOnly(x,par) + SmoothedStepOnly(x,par) + LinearOnly(x,par);
}

Double_t PndHypGePeakFitFunction::PeakFunc (Double_t *x,Double_t *par)
{
			// Peak shape is taken from gf2 http://www.phy.anl.gov/gammasphere/doc/gf2.hlp
			// composed by the sum of:
			//	(1) a Gaussian,
			//	(2) a skewed Gaussian,
			//	(3) a smoothed step function to increase the background on the low-energy
			//			side of the peak. Components (2) and/or (3) can easily be set to zero if not
			//			required, and
			//	(4)	a constant				// not in gf2
			//
			//	(2)	y = constant
			//			* EXP( (x-c)/beta )
			//			* ERFC( (x-c)/(SQRT(2)*sigma) + sigma/(SQRT(2)*beta) )
			//	(3)	y = constant * ERFC( (x-c)/(SQRT(2)*sigma) )
			
			// par[0]: amplitude gaus, par[1]: x0; par[2]: sigma, par[3]: amplitude step, par[4]: gradient of lin bg
			// par[5]: amplitude skewed gaus, par[6]: beta
			
	return PeakFuncGaussian (x,par) + SkewedOnly(x,par);
}

Double_t PndHypGePeakFitFunction::PeakFuncFreeSkewedPosition (Double_t *x,Double_t *par)
{
			// Peak shape is taken from gf2 http://www.phy.anl.gov/gammasphere/doc/gf2.hlp
			// composed by the sum of:
			//	(1) a Gaussian,
			//	(2) a skewed Gaussian,
			//	(3) a smoothed step function to increase the background on the low-energy
			//			side of the peak. Components (2) and/or (3) can easily be set to zero if not
			//			required, and
			//	(4)	a constant				// not in gf2
			//
			//	(2)	y = constant
			//			* EXP( (x-c)/beta )
			//			* ERFC( (x-c)/(SQRT(2)*sigma) + sigma/(SQRT(2)*beta) )
			//	(3)	y = constant * ERFC( (x-c)/(SQRT(2)*sigma) )
			
			// par[0]: amplitude gaus, par[1]: x0; par[2]: sigma, par[3]: amplitude step, par[4]: gradient of lin bg
			// par[5]: amplitude skewed gaus, par[6]: beta,  par[7]: x0_s, par[8]: sigma_s  
			
			//Double_t SkewedGausian = par[5] 
			//				* TMath::Exp((x[0]-par[7])/par[6]) 
			//				* TMath::Erfc( ((x[0] - par[7]) /(TMath::Sqrt(2)*par[8]) ) + par[8]/TMath::Sqrt(2)/par[6]); //par[5]: amplitude skewed gaus, par[6]: beta,  par[7]: x0_s, par[8]: sigma_s 
	//return PeakFuncGaussian (x,par) + FreeSkewedOnly(x,par);//SkewedGausian;
	return PeakFuncGaussian (x,par) + NewSkewedOnly(x,par);
	
}

Double_t PndHypGePeakFitFunction::PeakFuncFreeSkewedPositionWithoutFirstGausian (Double_t *x,Double_t *par)
{
	return FreeSkewedOnly(x,par) + SmoothedStepOnly(x,par) + LinearOnly(x,par);
}

Double_t PndHypGePeakFitFunction::PeakFuncDoubleGausian (Double_t *x,Double_t *par)
{
			// Peak shape is taken from gf2 http://www.phy.anl.gov/gammasphere/doc/gf2.hlp
			// composed by the sum of:
			//	(1) a Gaussian,
			//	(2) a second Gaussian,
			//	(3) a smoothed step function to increase the background on the low-energy
			//			side of the peak. Components (2) and/or (3) can easily be set to zero if not
			//			required, and
			//	(4)	a constant				// not in gf2
			//
			//	(3)	y = constant * ERFC( (x-c)/(SQRT(2)*sigma) )
			
			// par[0]: amplitude gaus, par[1]: x0; par[2]: sigma, par[3]: amplitude step, par[4]: gradient of lin bg
			// par[5]: amplitude second gaus, par[6]: x0_2,  par[7]: sigma_2  
			
	return PeakFuncGaussian (x,par) + SecondGausOnly(x,par);
}

// components
Double_t PndHypGePeakFitFunction::GausOnly (Double_t *x,Double_t *par)
{
			
			Double_t Gausian = par[0] * TMath::Exp(- 0.5*pow((x[0] - par[1])/par[2],2)); // par[0]: amplitude, par[1]: x[0]0; par[2]: sigma
			
	return Gausian ;
}

Double_t PndHypGePeakFitFunction::SmoothedStepOnly (Double_t *x,Double_t *par)
{
			Double_t SmoothedStep = par[3] * TMath::Erfc( (x[0] - par[1]) /(TMath::Sqrt(2)*par[2])); // par[3]: amplitude
	return  SmoothedStep;
}

Double_t PndHypGePeakFitFunction::LinearOnly (Double_t *x,Double_t *par)
{
			//Double_t Constant = par[4]*x[0]+par[5]; // par[4]: gradient of lin bg, par[5]: const bg
			Double_t Constant = par[5]; // par[4]: gradient of lin bg, par[5]: const bg
	return  Constant;
}
Double_t PndHypGePeakFitFunction::SkewedOnly (Double_t *x,Double_t *par)
{
	Double_t SkewedGausian = par[6] 
					* TMath::Exp((x[0]-par[1])/par[7]) 
					* TMath::Erfc( ((x[0] - par[1]) /(TMath::Sqrt(2)*par[2]) ) + par[2]/TMath::Sqrt(2)/par[7]); // par[6]: amplitude skewed gaus, par[7]: beta
	return SkewedGausian;
}
Double_t PndHypGePeakFitFunction::FreeSkewedOnly (Double_t *x,Double_t *par)
{
			Double_t SkewedGausian = par[6] 
							* TMath::Exp((x[0]-par[8])/par[7]) 
							* TMath::Erfc( ((x[0] - par[8]) /(TMath::Sqrt(2)*par[9]) ) + par[9]/TMath::Sqrt(2)/par[7]); //par[6]: amplitude skewed gaus, par[7]: beta,  par[8]: x0_s, par[9]: sigma_s 
	return  SkewedGausian;
}

Double_t PndHypGePeakFitFunction::SecondGausOnly (Double_t *x,Double_t *par)
{
			Double_t SecondGausian = par[6] * TMath::Exp(- 0.5*pow((x[0] - par[7])/par[8],2)); // par[6]: amplitude, par[7]: x0_2; par[8]: sigma_2
	return  SecondGausian;
}

Double_t PndHypGePeakFitFunction::NewSkewedOnly(Double_t *x, Double_t *par)
{
	Double_t NewSkewedGausian = par[6]
										* TMath::Exp(-1/2*pow((x[0]-par[8])/par[9],2))
										* (1 + TMath::Erf(1/sqrt(2)*par[7]*(x[0]-par[8])/par[9]));
	return NewSkewedGausian;
}

Double_t PndHypGePeakFitFunction::NewFunction(Double_t *x, Double_t *par)		// added 27.08.14 seems to fit the june 2014 dataset best, skewed gaus + gaus + constant
{
	Double_t value = par[3]*par[0]/2*TMath::Exp(-par[0]/2*(2*par[1]+par[0]*par[2]*par[2]-2*x[0])) * TMath::Erfc(-(par[1]+par[0]*par[2]*par[2]-x[0])/TMath::Sqrt(2)/par[2])+ par[4]*TMath::Gaus(x[0],par[5],par[6]) +par[7];
	return value;
}