/*
 * ScanPlot.cpp
 *
 *  Created on: Dec 24, 2014
 *      Author: Hovanes Egiyan
 */

#include "ScanPlot.hh"
#include "ScanCorrection.hh"

using namespace std;

// Default constructor
ScanPlot::ScanPlot()  {
	return;
}

// Main constructor
ScanPlot::ScanPlot( ScanDetector* det, ScanDetector* pos, string sName, bool statErrFlag ) :
		scanGraph(0), scanRawGraph(0), scanName(sName), scanDet(det), scanPos(pos),
		scanNormDet(0),  scanNormalized(false) {
	cout << "Creating new scan plot " << det->getName() << " vs " << pos->getName() << endl ;
	stringstream ssName;
	ssName << det->getName() << "_vs_" << pos->getName();
	scanPlotName = ssName.str();
	initialize( statErrFlag );
	return;
}


ScanPlot::~ScanPlot() {
	if( scanGraph != 0 )
			delete scanGraph;
	return;
}


void ScanPlot::initialize( bool statErrFlag ) {
	// Check if the detector and positioner have the same size
	// If not throw an exception
	if( scanDet->getData().size() != scanPos->getData().size() )
		throw std::runtime_error( string( "ERROR: Data ") + scanDet->getName() +
				string(" positions ") + scanPos->getName() + string(" have different lengths :")  );
	cout << "Setting axes for plot " << scanDet->getName() << " vs " << scanPos->getName() << endl ;

	vector<double> xVec =  scanPos->getData();
	vector<double> yVec =  scanDet->getData();
	unsigned NumberOfPoints = xVec.size();

	cout << "Got vectors " << endl;

	// Create and assign double arrays
	double xVal[NumberOfPoints];
	double yVal[NumberOfPoints];
	double xErr[NumberOfPoints];
	double yErr[NumberOfPoints];

	cout << "Created double arrays" << endl;
	for( unsigned iPt = 0; iPt < scanPos->getData().size(); iPt++ ) {
		xVal[iPt] = xVec[iPt];
		yVal[iPt] = yVec[iPt];
		xErr[iPt] = 0.0;
		yErr[iPt] = 0.0 ;
		if( statErrFlag ) yErr[iPt] = sqrt( yVec[iPt] );
	}

	cout << "Assigned double arrays " << endl;

	// Create the graph
	scanGraph = new TGraphErrors( scanPos->getData().size(), xVal, yVal, xErr, yErr );
	scanGraph->SetTitle( scanPlotName.c_str() );
	scanGraph->SetName( scanPlotName.c_str() );
	scanGraph->SetFillColor(kYellow);
	scanGraph->GetXaxis()->SetTitle( (scanPos->getName() + string(" (") + scanPos->getUnit() + string(")")).c_str() );
	scanGraph->GetYaxis()->SetTitle( (scanDet->getName() + string(" (") + scanDet->getUnit() + string(")")).c_str() );
	cout << "Set axes titles as " << scanGraph->GetXaxis()->GetTitle() << " , " << scanGraph->GetYaxis()->GetTitle() << endl;
	double xMin, xMax, yMin, yMax;
	scanGraph->ComputeRange( xMin, yMin, xMax, yMax );
	if( yMax > 0 ) {
		scanGraph->SetMaximum( 2.0 * yMax );
		scanGraph->SetMinimum( 2.0e-06 * yMax );
	}

	// At start both raw and normalized graphs are the same
	scanRawGraph = new TGraphErrors( *scanGraph );
	scanRawGraph->GetXaxis()->SetTitle( scanGraph->GetXaxis()->GetTitle() );
	scanRawGraph->GetYaxis()->SetTitle( scanGraph->GetYaxis()->GetTitle() );

	return;
}

bool ScanPlot::normalize() {

	if( scanRawGraph->GetN() == 0 ) {
		cerr << "WARNIGN: Number of points in the detector is zero. Will not normalize." << endl;
		return false;
	}

	if ( static_cast<int>( scanNormDet->getData().size() ) != scanRawGraph->GetN() ) {
		cerr << "WARNIGN: Number of points in the normalization detector is "
				<< scanNormDet->getData().size() << " instead of " << scanRawGraph->GetN()
				<< " . Cannot normalize." << endl;
		return false;
	}
	// Calculate the normalized vector from normalization detector
	vector<double> tempVec( scanNormDet->getData() );
	double minVal = *std::min_element( tempVec.begin(), tempVec.end() );
	if ( minVal < 1.0e-10 ) {
		cerr << "WARNING: normalization detector contains at least one very small element "
				<< minVal << " . Cannot not normalize " << endl;
		return false;
	}
	double vectorSum = std::accumulate( tempVec.begin(), tempVec.end(), 0.0 );
	double scaleFactor = vectorSum / tempVec.size();

	// Transform tempVec into normalized vector
	std::transform( tempVec.begin(), tempVec.end(), tempVec.begin(),
			std::bind1st( std::divides<double>(), scaleFactor ) );

	// Loop over all points and recalculate the value for each data point
	// using a newly calculated normalization factor.
	string xAxis = scanGraph->GetXaxis()->GetTitle();
	string yAxis = scanGraph->GetYaxis()->GetTitle();
	for ( int iPt = 0; iPt < scanRawGraph->GetN(); iPt++ ) {
		double yTemp = tempVec[iPt];
		double xRaw, yRaw;
		if( yTemp > 1.0e-20 ) {
			scanRawGraph->GetPoint( iPt, xRaw, yRaw );
			scanGraph->SetPoint( iPt, xRaw, yRaw / yTemp );
			scanGraph->SetPointError( iPt, scanGraph->GetErrorX( iPt ),
					scanGraph->GetErrorY( iPt ) / yTemp );
		}
	}

	scanGraph->GetXaxis()->SetTitle( xAxis.c_str() );
	scanGraph->GetYaxis()->SetTitle( yAxis.c_str() );

	this->setNormailzied();
	return true;
}

void ScanPlot::applyCorrections( ScanCorrection& corr ) {
	string xAxis = scanGraph->GetXaxis()->GetTitle();
	string yAxis = scanGraph->GetYaxis()->GetTitle();

	for( int iPt = 0; iPt < scanGraph->GetN(); iPt++ ) {
		double xOld, yOld;

		scanGraph->GetPoint( iPt, xOld, yOld );
		double exOld = scanGraph->GetErrorX( iPt );
		double eyOld = scanGraph->GetErrorY( iPt );

		double correctiomFactor = corr.getCorrection( iPt ) ;

		double yNew = yOld * correctiomFactor ;
		double eyNew = eyOld * correctiomFactor;
//		cout << "Correction factor for point " << iPt << " is " << correctiomFactor << endl;
//
//		cout << "Old y is " << yOld << " , yNew is " << yNew << endl;

		scanGraph->SetPoint( iPt, xOld, yNew );
		scanGraph->SetPointError( iPt, exOld, eyNew );
	}
	scanGraph->GetXaxis()->SetTitle( xAxis.c_str() );
	scanGraph->GetYaxis()->SetTitle( yAxis.c_str() );

	return;
}