#include "DCrossSectionData.h"

ClassImp(DCrossSectionData)

DCrossSectionData::DCrossSectionData(string locReactionName, string locBeamROOTString, string locAngleParticleROOTString, double locBeamMass, double locTargetMass) : 
dReactionName(locReactionName), dBeamROOTString(locBeamROOTString), dAngleParticleROOTString(locAngleParticleROOTString), 
dBeamMass(locBeamMass), dTargetMass(locTargetMass) {}

DCrossSectionData::DCrossSectionData(DCrossSectionData* locCrossSectionData_ConvertedFrom, string locGraphName, TObjArray* locCrossSectionArray, const vector<double>& locFirstDimValues)
{
	*this = *locCrossSectionData_ConvertedFrom; //not ideal, but whatever
	dGraphName = locGraphName;
	dCrossSectionArray = locCrossSectionArray;
	dFirstDimValues = locFirstDimValues;

//	Set_GraphTitles();
}

void DCrossSectionData::Set_Data(string locGraphName, TObjArray* locCrossSectionArray, const vector<double>& locFirstDimValues, bool locIsSqrtSFlag, bool locIsCosThetaFlag, bool locGraphPerEnergyBinFlag, bool locIsDSigmaDCosThetaFlag)
{
	dGraphName = locGraphName;

	dCrossSectionArray = locCrossSectionArray;
	dFirstDimValues = locFirstDimValues;
	dIsSqrtSFlag = locIsSqrtSFlag;
	dIsCosThetaFlag = locIsCosThetaFlag;
	dGraphPerEnergyBinFlag = locGraphPerEnergyBinFlag;
	dIsDSigmaDCosThetaFlag = locIsDSigmaDCosThetaFlag;

//	Set_GraphTitles();
}

DCrossSectionData* DCrossSectionData::Convert_Data(string locNewGraphName, bool locNewIsSqrtSFlag, bool locNewIsCosThetaFlag, bool locNewIsDSigmaDCosThetaFlag)
{
	//locNewGraphPerEnergyBinFlag IS NOT YET SUPPORTED
	DCrossSectionData* locCrossSectionData = new DCrossSectionData(dReactionName, dBeamROOTString, dAngleParticleROOTString, dBeamMass, dTargetMass);

	vector<double> locNewFirstDimValues;

	//1st Dimension
	if(dGraphPerEnergyBinFlag)
	{
		if(locNewIsSqrtSFlag != dIsSqrtSFlag)
		{
			for(size_t loc_i = 0; loc_i < dFirstDimValues.size(); ++loc_i)
			{
				if(locNewIsSqrtSFlag)
					locNewFirstDimValues.push_back(Convert_BeamEToRootS(dFirstDimValues[loc_i]));
				else
					locNewFirstDimValues.push_back(Convert_RootSToBeamE(dFirstDimValues[loc_i]));
			}
		}
		else
			locNewFirstDimValues = dFirstDimValues;
	}
	else
	{
		if(locNewIsCosThetaFlag != dIsCosThetaFlag)
		{
			for(size_t loc_i = 0; loc_i < dFirstDimValues.size(); ++loc_i)
			{
				if(locNewIsCosThetaFlag)
					locNewFirstDimValues.push_back(cos(dFirstDimValues[loc_i]*TMath::Pi()/180.0));
				else
					locNewFirstDimValues.push_back(acos(dFirstDimValues[loc_i])*180.0/TMath::Pi());
			}
		}
		else
			locNewFirstDimValues = dFirstDimValues;
	}

	TObjArray* locNewCrossSectionArray = new TObjArray();
	for(Int_t loc_i = 0; loc_i < dCrossSectionArray->GetEntriesFast(); ++loc_i)
	{
		TGraphErrors* locGraphErrors = (TGraphErrors*)dCrossSectionArray->At(loc_i);
		int locNumPoints = locGraphErrors->GetN();

		//2nd Dimension
		double *locSecondDimArray = locGraphErrors->GetX();
		double *locNewSecondDimArray = NULL;
		if(dGraphPerEnergyBinFlag)
		{
			if(locNewIsCosThetaFlag != dIsCosThetaFlag)
			{
				locNewSecondDimArray = new double[locNumPoints];
				for(int loc_j = 0; loc_j < locNumPoints; ++loc_j)
				{
					if(locNewIsCosThetaFlag)
						locNewSecondDimArray[loc_j] = cos(locSecondDimArray[loc_j]*TMath::Pi()/180.0);
					else
						locNewSecondDimArray[loc_j] = acos(locSecondDimArray[loc_j])*180.0/TMath::Pi();
				}
			}
			else
				locNewSecondDimArray = locSecondDimArray;
		}
		else
		{
			if(locNewIsSqrtSFlag != dIsSqrtSFlag)
			{
				locNewSecondDimArray = new double[locNumPoints];
				for(int loc_j = 0; loc_j < locNumPoints; ++loc_j)
				{
					if(locNewIsSqrtSFlag)
						locNewSecondDimArray[loc_j] = Convert_BeamEToRootS(locSecondDimArray[loc_j]);
					else
						locNewSecondDimArray[loc_j] = Convert_RootSToBeamE(locSecondDimArray[loc_j]);
				}
			}
			else
				locNewSecondDimArray = locSecondDimArray;
		}

		//cross value & error
		double* locCrossValueArray = locGraphErrors->GetY();
		double* locCrossErrorArray = locGraphErrors->GetEY();
		double *locNewCrossValueArray, *locNewCrossErrorArray;
		if(locNewIsDSigmaDCosThetaFlag != dIsDSigmaDCosThetaFlag)
		{
			locNewCrossValueArray = new double[locNumPoints];
			locNewCrossErrorArray = new double[locNumPoints];
			for(int loc_j = 0; loc_j < locNumPoints; loc_j++)
			{
				if(locNewIsDSigmaDCosThetaFlag)
				{
					locNewCrossValueArray[loc_j] = locCrossValueArray[loc_j]*2.0*TMath::Pi();
					locNewCrossErrorArray[loc_j] = locCrossErrorArray[loc_j]*2.0*TMath::Pi();
				}
				else
				{
					locNewCrossValueArray[loc_j] = locCrossValueArray[loc_j]/2.0*TMath::Pi();
					locNewCrossErrorArray[loc_j] = locCrossErrorArray[loc_j]/2.0*TMath::Pi();
				}
			}
		}
		else
		{
			locNewCrossValueArray = locCrossValueArray;
			locNewCrossErrorArray = locCrossErrorArray;
		}

		//reverse if needed
		if(locNewSecondDimArray[0] > locNewSecondDimArray[1])
		{
			locNewSecondDimArray = Reverse_Order(locNewSecondDimArray, locNumPoints);
			locNewCrossValueArray = Reverse_Order(locNewCrossValueArray, locNumPoints);
			locNewCrossErrorArray = Reverse_Order(locNewCrossErrorArray, locNumPoints);
		}

		TGraphErrors* locNewGraphErrors = new TGraphErrors(locNumPoints, locNewSecondDimArray, locNewCrossValueArray, NULL, locNewCrossErrorArray);
		locNewCrossSectionArray->AddLast(locNewGraphErrors);
	}

	locCrossSectionData->Set_Data(locNewGraphName, locNewCrossSectionArray, locNewFirstDimValues, locNewIsSqrtSFlag, locNewIsCosThetaFlag, dGraphPerEnergyBinFlag, locNewIsDSigmaDCosThetaFlag);
	return locCrossSectionData;
}

void DCrossSectionData::Set_GraphTitles(void)
{
	for(Int_t loc_i = 0; loc_i < dCrossSectionArray->GetEntriesFast(); ++loc_i)
	{
		TGraphErrors* locGraphErrors = (TGraphErrors*)dCrossSectionArray->At(loc_i);
		if(locGraphErrors == NULL)
			continue;

		ostringstream locGraphNameStream, locGraphTitleStream;
		locGraphNameStream << dGraphName << "_" << loc_i;
		locGraphErrors->SetName(locGraphNameStream.str().c_str());

		//Title
		locGraphTitleStream << dReactionName << ", ";
		if(dGraphPerEnergyBinFlag)
		{
			if(dIsSqrtSFlag)
				locGraphTitleStream << "#sqrt{s} = " << dFirstDimValues[loc_i] << " (GeV);";
			else
				locGraphTitleStream << "E_{" << dBeamROOTString << "} = " << dFirstDimValues[loc_i] << " (GeV);";
		}
		else
		{
			if(dIsCosThetaFlag)
				locGraphTitleStream << "cos(#theta_{" << dAngleParticleROOTString << "}) = " << dFirstDimValues[loc_i] << ";";
			else
				locGraphTitleStream << "#theta_{" << dAngleParticleROOTString << "} = " << dFirstDimValues[loc_i] << "#circ;";
		}

		//X-axis
		if(dGraphPerEnergyBinFlag)
		{
			if(dIsCosThetaFlag)
				locGraphTitleStream << "cos(#theta_{" << dAngleParticleROOTString << "});";
			else
				locGraphTitleStream << "#theta_{" << dAngleParticleROOTString << "}#circ;";
		}
		else
		{
			if(dIsSqrtSFlag)
				locGraphTitleStream << "#sqrt{s} (GeV);";
			else
				locGraphTitleStream << "E_{" << dBeamROOTString << "} (GeV);";
		}

		//Y-axis
		if(dIsDSigmaDCosThetaFlag)
			locGraphTitleStream << "d#sigma / dcos(#theta_{" << dAngleParticleROOTString << "}) (#mub);";
		else
			locGraphTitleStream << "d#sigma / d#Omega (#mub/sr);";

		locGraphErrors->SetTitle(locGraphTitleStream.str().c_str());
	}
}

double* DCrossSectionData::Reverse_Order(double* locInputArray, size_t locSize)
{
	double* locOutputArray = new double[locSize];
	for(size_t loc_i = 0; loc_i < locSize; ++loc_i)
		locOutputArray[locSize - loc_i - 1] = locInputArray[loc_i];
	return locOutputArray;
}