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

#include "ProfilerPlaneSliceSimulator.hh"

using namespace std;

// Main constructor
ProfilerPlaneSliceSimulator::ProfilerPlaneSliceSimulator( const std::string name,
		const std::vector<double>& xAxis, const std::vector<double>& yAxis ) :
		MutexedClass(), ProfilerPlaneSlice( name, xAxis, yAxis ), ppssHist( 0 ), ppssFunc( 0 ),
				ppssRandGenerator() {
	struct timespec tms;
	clock_gettime( CLOCK_REALTIME, &tms );
	ppssRandGenerator.SetSeed( tms.tv_nsec );
	return;
}

// Copy constructor
ProfilerPlaneSliceSimulator::ProfilerPlaneSliceSimulator( const ProfilerPlaneSliceSimulator& s ) :
		MutexedClass( s ), ProfilerPlaneSlice( s ), ppssHist( 0 ), ppssFunc( 0 ),
				ppssRandGenerator() {
	ppssRandGenerator.SetSeed( s.ppssRandGenerator.GetSeed() );
	// Clone the histogram and the function
	std::string newHistName = std::string( ppssHist->GetName() ) + ":clone";
	ppssHist = reinterpret_cast<TH1D*>( s.ppssHist->Clone( newHistName.c_str() ) );
	std::string newFuncName = std::string( ppssFunc->GetName() ) + ":clone";
	ppssFunc = reinterpret_cast<TF1*>( s.ppssFunc->Clone( newFuncName.c_str() ) );
	return;
}

// destructor
ProfilerPlaneSliceSimulator::~ProfilerPlaneSliceSimulator() {
	MtxLock objLock( mcMutex );
	// delete the histogram and the function if they exist
	if ( ppssHist != 0 && !ppssHist->IsZombie() ) {
		delete ppssHist;
		ppssHist = 0;
	}
	if ( ppssFunc != 0 && !ppssFunc->IsZombie() ) {
		delete ppssFunc;
		ppssFunc = 0;
	}
	return;
}

// Assignment operator
ProfilerPlaneSliceSimulator& ProfilerPlaneSliceSimulator::operator =(
		const ProfilerPlaneSliceSimulator& s ) {
	MutexedClass::operator =( s );
	ProfilerPlaneSlice::operator =( s );
	MtxLock objLock( mcMutex );
	// Delete the old histogram and function
	if ( ppssHist != 0 && !ppssHist->IsZombie() ) {
		delete ppssHist;
		ppssHist = 0;
	}
	if ( ppssFunc != 0 && !ppssFunc->IsZombie() ) {
		delete ppssFunc;
		ppssFunc = 0;
	}
	ppssRandGenerator.SetSeed( s.ppssRandGenerator.GetSeed() );
	// Clone the histogram and function
	std::string newHistName = std::string( ppssHist->GetName() ) + ":clone";
	ppssHist = dynamic_cast<TH1D*>( s.ppssHist->Clone( newHistName.c_str() ) );
	std::string newFuncName = std::string( ppssFunc->GetName() ) + ":clone";
	ppssFunc = dynamic_cast<TF1*>( s.ppssFunc->Clone( newFuncName.c_str() ) );
	return *this;
}

// Generate the values and fill the vector
void ProfilerPlaneSliceSimulator::simulateValues( const vector<double>& effcCorr ) {
	MtxLock objLock( mcMutex );

	// Loop over all bins/slices of the histogram and generate a random value in that bin.
	// The number is generated in the center of that bin using the function member of
	// this class. The number is distributed according to Poisson distribution but scaled
	// by the value of efficiency vector requested.
	for ( Int_t iBin = 1; iBin <= ppssHist->GetNbinsX(); iBin++ ) {
//		double xVal = ppssHist->GetBinCenter( iBin );
		double xMin = ppssHist->GetBinLowEdge( iBin );
		double xMax = xMin + ppssHist->GetBinWidth( iBin );
//		cout << "Simulating between " << xMin << " and " << xMax << " as bin " << iBin << " for " << getName() << endl;
		double meanValue = ppssFunc->Integral( xMin, xMax ) / ppssHist->GetBinWidth( iBin ) ;
//		double meanValue = ppssFunc->Eval( xVal );
		double effcCorrVal = 1.0;
		if( effcCorr[iBin-1] > 1.0e-20 ) effcCorrVal = effcCorr[iBin-1];
		double rndmNumber = ppssRandGenerator.PoissonD( meanValue ) / effcCorrVal;
		ppssHist->SetBinContent( iBin, rndmNumber );
	}

//	cout << "Amplitude of the histogram is " << ppssHist->GetMaximum() <<
//			" , mean is " << ppssHist->GetMean() << " , RMS is " << ppssHist->GetRMS() << endl;
//
//	cout << "Here is the content of  the simulation function for " << getName() << endl;
//	for( Int_t iBin = 1; iBin <= ppssHist->GetNbinsX(); iBin++ ) {
//		double xVal = ppssHist->GetBinCenter( iBin );
//		cout << ppssFunc->Eval( xVal ) << " " ;
//	}
//	cout << endl;
//	cout << "Here is the content of simulated histogram for " << getName() << " with address " <<
//			hex << showbase << ppssHist << dec << endl;
//	for( Int_t iBin = 1; iBin <= ppssHist->GetNbinsX(); iBin++ ) {
//		cout << ppssHist->GetBinContent( iBin ) << " " ;
//	}
//	cout << endl;
	return;
}

// Creates a new TF1 object and return the pointer
TF1* ProfilerPlaneSliceSimulator::defineDistFunc() {
//	cout << "In " << __FUNCTION__ << endl;
	MtxLock objLock( mcMutex );
	std::string funName = ppsName + "_f";
	// Create a ROOT TF1
	double xMin = 0;
	double xMax = 0;
	if ( ppsAxis.size() > 1 ) {
		double deltaX = (ppsAxis[ppsAxis.size() - 1] - ppsAxis[0]) / (ppsAxis.size() - 1);
		xMin = ppsAxis[0] - deltaX / 2.0;
		xMax = ppsAxis[ppsAxis.size() - 1] + deltaX / 2.0;
	}
//	cout << "Creating TF1 function between " << xMin << " and " << xMax << " for " << getName()
//			<< endl;
	TF1* rootFunPtr = new TF1( funName.c_str(),
			(Double_t (*)( Double_t*, Double_t* ) ) ProfilerPlaneSliceSimulator::cauchyDistribution, xMin, xMax, 5 );

	// Set parameter names
	rootFunPtr->SetParName( 0, "Amplitude" );
	rootFunPtr->SetParName( 1, "Mean" );
	rootFunPtr->SetParName( 2, "Width" );
	rootFunPtr->SetParName( 3, "Background" );
	rootFunPtr->SetParName( 4, "BkgSlope" );
	// Set parameter values
	rootFunPtr->SetParameter( 0, ppsAmpl );
	rootFunPtr->SetParameter( 1, ppsMean );
	rootFunPtr->SetParameter( 2, ppsWidth );
	rootFunPtr->SetParameter( 3, ppsBkg );
	rootFunPtr->SetParameter( 4, ppsBkgSlope );

//	cout << "Leaving " << __FUNCTION__ << endl;
	ppssFunc = rootFunPtr;
//	ppssFunc->Dump();
	return ppssFunc;
}

// Creates a new TH1D object and returns the pointer
TH1D* ProfilerPlaneSliceSimulator::defineHistogram() {
//	cout << "In " << __FUNCTION__ << endl;
	MtxLock objLock( mcMutex );
	std::string histName = ppsName + "_h";
	std::string histTitle = ppsName;
	double xMin = 0;
	double xMax = 0;
	if ( ppsAxis.size() > 1 ) {
		double deltaX = (ppsAxis[ppsAxis.size() - 1] - ppsAxis[0]) / (ppsAxis.size() - 1);
		xMin = ppsAxis[0] - deltaX / 2.0;
		xMax = ppsAxis[ppsAxis.size() - 1] + deltaX / 2.0;
	}
	TH1D* rootHistPtr = new TH1D( histName.c_str(), histTitle.c_str(), ppsAxis.size(), xMin, xMax );
//	rootHistPtr->Dump();
//	cout << "Leaving " << __FUNCTION__ << endl;
//	cout << "Histogram address is " << rootHistPtr << endl;
	ppssHist = rootHistPtr;
	return ppssHist;
}

// Cauchy (Breit-Wigner) distribution. First argument is the pointer to the
// x-value where the function needs to be evaluated, the second argument is the
// pointer to the parameters.
double ProfilerPlaneSliceSimulator::cauchyDistribution( double* x, double* par ) {
	double alpha = par[0]; // amplitude
	double x0 = par[1]; // median
	double gamma = par[2]/2.0; // half width
	double bkg = par[3]; // flat background
	double slope = par[4]; // background slope

	double x_i = x[0];

//	cout << "Calculating function at " << x_i << " with parameters " <<
//			" alpha= " << alpha << " , gamma= " << gamma << " , x0= " << x0 << endl;

	double value = alpha / (1.0 + pow( (x_i - x0) / gamma, 2.0 )) + bkg * ( 1.0 + slope * x_i ) ;
//	cout << "Value is " << value << endl;
	return value;

}