#ifndef DCrossSectionData_h
#define DCrossSectionData_h

#include <vector>
#include <string>
#include <math.h>
#include <sstream>

#include "TObject.h"
#include "TMath.h"
#include "TObjArray.h"
#include "TGraphErrors.h"

using namespace std;

class DCrossSectionData : public TObject
{
	public:
		DCrossSectionData(string locReactionName = "", string locBeamROOTString = "", string locAngleParticleROOTString = "", double locBeamMass = 0.0, double locTargetMass = 0.0);
		DCrossSectionData(DCrossSectionData* locCrossSectionData_ConvertedFrom, string locGraphName, TObjArray* locCrossSectionArray, const vector<double>& locFirstDimValues);

		void Set_Data(string locGraphName, TObjArray* locCrossSectionArray, const vector<double>& locFirstDimValues, bool locIsSqrtSFlag, bool locIsCosThetaFlag, bool locGraphPerEnergyBinFlag, bool locIsDSigmaDCosThetaFlag);
		void Set_GraphTitles(void);
		DCrossSectionData* Convert_Data(string locNewGraphName, bool locNewIsSqrtSFlag, bool locNewIsCosThetaFlag, bool locNewIsDSigmaDCosThetaFlag);

		double Convert_BeamEToRootS(double locBeamE) const;
		double Convert_RootSToBeamE(double locRootS) const;

		string dReactionName;
		string dBeamROOTString;
		string dAngleParticleROOTString;
		double dBeamMass;
		double dTargetMass;

		string dGraphName;

//		bool dIsTotalCrossSectionFlag; //array is size 1 //not yet supported

		bool dIsSqrtSFlag; //else is E_beam
		bool dIsCosThetaFlag; //else is theta
		bool dGraphPerEnergyBinFlag; //true if each graph is an enegy bin. else is angle bins //irrelevant if dIsTotalCrossSectionFlag = true
		bool dIsDSigmaDCosThetaFlag; //false if is dOmega

		TObjArray* dCrossSectionArray; //Array of TGraphErrors*
		vector<double> dFirstDimValues; //e.g. if each graph is an e-bin, is central energy //empty if dIsTotalCrossSectionFlag = true

	private:
		double* Reverse_Order(double* locInputArray, size_t locSize);

	ClassDef(DCrossSectionData, 1)
};

inline double DCrossSectionData::Convert_BeamEToRootS(double locBeamE) const
{
	//invariant mass:
	//E: locBeamE + dTargetMass
	//P: sqrt(locBeamE*locBeamE - locBeamMass*locBeamMass)
	//M2 = (locBeamE + dTargetMass)*(locBeamE + dTargetMass) - locBeamE*locBeamE + locBeamMass*locBeamMass
	//M2 = locBeamE*locBeamE + 2*locBeamE*dTargetMass + dTargetMass*dTargetMass - locBeamE*locBeamE + locBeamMass*locBeamMass
	//M2 = 2*locBeamE*dTargetMass + dTargetMass*dTargetMass + locBeamMass*locBeamMass
	return sqrt(2.0*dTargetMass*locBeamE + dTargetMass*dTargetMass + dBeamMass*dBeamMass);
}

inline double DCrossSectionData::Convert_RootSToBeamE(double locRootS) const
{
	return (locRootS*locRootS - dTargetMass*dTargetMass - dBeamMass*dBeamMass)/(2.0*dTargetMass);
}

#endif //DCrossSectionData_h