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

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

ClassImp( ScanPlot )

using namespace std;

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

// Main constructor
ScanPlot::ScanPlot( ScanDetector* det, ScanDetector* posX, ScanDetector* posY, string sName,
		string date, bool statErrFlag ) :
		scanGraph( 0 ), scanRawGraph( 0 ), scanName( sName ), scanPlotName( "" ), scanDate( date ),
				scanDet( det ), scanPosX( posX ), scanPosY( posY ), scanNormDet( 0 ),
				scanNormalized( false ) {
	cout << "Creating new scan plot " << det->getName() << " vs " << posY->getName() << " and "
			<< posX->getName() << endl;
	stringstream ssName;
	ssName << det->getName() << "_vs_" << posY->getName() << "_vs_" << posX->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
	cout << "X data size is " << scanPosX->getData().size() << " " << "Y data size is "
			<< scanPosY->getData().size() << " " << "Z data size is " << scanDet->getData().size()
			<< endl;
	if ( (scanDet->getData().size() != scanPosX->getData().size())
			|| (scanDet->getData().size() != scanPosY->getData().size())
			|| (scanPosX->getData().size() != scanPosY->getData().size()) ) {
		cout << "X data size is " << scanPosX->getData().size() << " " << "Y data size is "
				<< scanPosY->getData().size() << " " << "Z data size is "
				<< scanDet->getData().size() << endl;
		throw std::runtime_error(
				string( "ERROR: Data " ) + scanDet->getName()
						+ string( " and positions have different lengths :" ) );
	}

	cout << "Setting axes for plot " << scanDet->getName() << " vs " << scanPosX->getName()
			<< " vs " << scanPosX->getName() << endl;

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

	cout << "Got vectors " << endl;

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

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

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

	// Create the graph
	scanGraph = new TGraph2DErrors( scanPosX->getData().size(), xVal, yVal, zVal, xErr, yErr,
			zErr );
	scanGraph->SetTitle( scanPlotName.c_str() );
	scanGraph->SetName( scanPlotName.c_str() );
	scanGraph->SetFillColor( kYellow );
	scanGraph->GetXaxis()->SetTitle(
			(scanPosX->getName() + string( " (" ) + scanPosX->getUnit() + string( ")" )).c_str() );
	scanGraph->GetYaxis()->SetTitle(
			(scanPosY->getName() + string( " (" ) + scanPosY->getUnit() + string( ")" )).c_str() );
//	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 );
//	}

//	scanGraph->SetMinimum( 0 );
//	scanGraph->SetMaximum( 1000 );
	// At start both raw and normalized graphs are the same
	scanRawGraph = new TGraph2DErrors( *scanGraph );
	scanRawGraph->SetTitle( scanPlotName.c_str() );
	scanRawGraph->SetFillColor( kYellow );
	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.
	double* xRaw = scanRawGraph->GetX();
	double* yRaw = scanRawGraph->GetY();
	double* zRaw = scanRawGraph->GetZ();
	for ( int iPt = 0; iPt < scanRawGraph->GetN(); iPt++ ) {
		double zTemp = tempVec[iPt];
		scanGraph->SetPoint( iPt, xRaw[iPt], yRaw[iPt], zRaw[iPt] / zTemp );
		scanGraph->SetPointError( iPt, scanGraph->GetEX()[iPt], scanGraph->GetEY()[iPt],
				scanGraph->GetZ()[iPt] / zTemp );
	}

	this->setNormailzied();
	return true;
}

void ScanPlot::applyCorrections( ScanCorrection& corr ) {
	for ( int iPt = 0; iPt < scanGraph->GetN(); iPt++ ) {
		double xOld = scanRawGraph->GetX()[iPt];
		double yOld = scanRawGraph->GetY()[iPt];
		double zOld = scanRawGraph->GetZ()[iPt];

		double exOld = scanGraph->GetEX()[iPt];
		double eyOld = scanGraph->GetEY()[iPt];
		double ezOld = scanGraph->GetEZ()[iPt];

		double correctiomFactor = corr.getCorrection( iPt );

		double zNew = zOld * correctiomFactor;
		double ezNew = ezOld * correctiomFactor;
//		cout << "Correction factor for point " << iPt << " is " << correctiomFactor << endl;
//
//		cout << "Old y is " << yOld << " , yNew is " << yNew << endl;

		scanGraph->SetPoint( iPt, xOld, yOld, zNew );
		scanGraph->SetPointError( iPt, exOld, eyOld, ezNew );
	}
	return;
}