#include "Particle.h"
#include <sstream>
#include <string>

/////////////////////////////////////////////////////////
//
// Particle
//
// This class enables you to set
// a mass and an impulse of a particle,
// change values and get them back.
// fA and fZ variables are still ambiguous
// and need to be discussed.
// Each particle also stores information about its energy states
//
////////////////////////////////////////////////////////

ExcitationState::ExcitationState() {
	//Default constructor
	fMean = 0.;
	fWidth = 0.1;
	fShape = "gauss";
	fStrength = 1;
}
;

ExcitationState::ExcitationState(Double_t mean, Double_t width, TString shape,
		Int_t strength) :
		fMean(mean), fWidth(width), fShape(shape), fStrength(strength) {
	//Constructor initializing all necessary variables for exc. state
}

TString ExcitationState::CreateConfigString() {
	ConfigDictionary CD;
	CD.SetDouble("mean", fMean);
	CD.SetDouble("width", fWidth);
	CD.SetString("shape", fShape.Data());
	CD.SetInt("strength", fStrength);
	return TString(CD.ToString().c_str());
}
;

int ExcitationState::ReadConfigString(TString cs) {
	ConfigDictionary CD(cs.Data());
	try {
		fMean = CD.GetDouble("mean");
		fWidth = CD.GetDouble("width");
		fShape = CD.GetString("shape").c_str();
		fStrength = CD.GetInt("strength");
		return SUCCESS;
	} catch (std::string & e) {
		Error("ExcitationState::ReadConfigString",
				"Couldn't find parameter: %s", e.c_str());
		return NOTFOUND;
	}
	if (fWidth <= 0.) {
		Error("ExcitationState::ReadConfigString", "Width parameter <= 0!");
		return FAILURE;
	}
	if (fMean <= 0.) {
		Error("ExcitationState::ReadConfigString", "Mean parameter <= 0!");
		return FAILURE;
	}
	if (fStrength <= 0) {
		Error("ExcitationState::ReadConfigString", "Strength parameter <= 0!");
		return FAILURE;
	}

	return SUCCESS;
}
;

/////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////

ClassImp(Particle);

Int_t Particle::numberOfParticles = 0;

Particle::Particle() {
// default constructor
// set Z = 0, A = 0
	Reset();
	fZ = 0;
	fA = 0;
	numberOfParticles++;
	std::stringstream ss;
	ss << "unnamed particle no. ";
	ss << numberOfParticles;
	SetName(ss.str().c_str());

//	cout << "kjadbhfkjasdf" << endl;

}
;

Particle::Particle(const char *name, Double_t mass, Int_t A, Int_t Z,
		Bool_t obs) {
	//name: name of the object
	//mass: rest-mass of the object, in MeV
	//A: mass number, number of nucleons
	//Z: charge, equals number in periodic table
	//set total energy equal to the mass of the particle

	fName = name;
	fMass = mass;
	fGroundStateMass = mass;
	Info("Particle::Particle", "Ground state mass was set to %f",
			fGroundStateMass);
	fA = A;
	fZ = Z;
	fObservable = obs;
	fImpulse.SetE(mass);

}

Particle::Particle(TString cs) {
//config constructor
	ReadConfigString(cs);
	CreateStateMassFunctions();
}

Particle::~Particle() {
//destructor
	for (UInt_t i = 0; i < fState.size(); i++) {
		delete fState[i];
//		Info("DetManager::~DetManager", "Detector %s deleting", det->GetName());
//		delete det;
	}
	fState.clear();

}

//public functions

//_____________________________________________________________________________
void Particle::AddExcitationState(ExcitationState exstate) {
//add new excitation states to the particle
	fExcitationStates.push_back(exstate);
}

//_____________________________________________________________________________
void Particle::ClearExcitationStates() {
//clear excitation states
	fExcitationStates.clear();
}

//_____________________________________________________________________________
Int_t Particle::GetNumberOfExStates() {
//return number of excitation states
	return fExcitationStates.size();
}

//_____________________________________________________________________________
ExcitationState Particle::GetExcitationState(Int_t index) {
//return "index"th excitation state
	return fExcitationStates[index];
}

void Particle::Reset() {
	// Set mass and impulse of particle as zero value

	fMass = 0.;
	fGroundStateMass = 0.;
	fImpulse.SetPxPyPzE(0., 0., 0., 0.);
	return;
}

int Particle::CopyValues(Particle * other) {
	//Copies kinetic energy and direction angles from other particle
	//Dont check if other == null to increase speed, be careful!
	fImpulse.SetPxPyPzE(other->GetPx(), other->GetPy(), other->GetPz(),
			other->GetE());
	return SUCCESS;
}

void Particle::SetMPxPyPz(Double_t m, Double_t px, Double_t py, Double_t pz) {
	// m: set mass in the MeV
	// px: x cordinate of impulse
	// py: y cordinate of impulse
	// pz: z cordinate of impulse

	fImpulse.SetPx(px);
	fImpulse.SetPy(py);
	fImpulse.SetPz(pz);
	fImpulse.SetE(CalcT(m, px, py, pz) + m);
	fMass = m;
}

void Particle::SetMTNxNyNz(Double_t m, Double_t T, Double_t nx, Double_t ny,
		Double_t nz) {
// m: set mass in MeV
// T: set kinetic energy in MeV
// set direction vector of impulse with (nx,ny,nz) components
//todo: function restored, check if it works properly

	fMass = m;
	TVector3 vect(nx, ny, nz);
	vect.SetMag(CalcP(m, T));
	fImpulse.SetVectMag(vect, m);
	fImpulse.SetE(Sqrt(vect.Mag2() + m * m));

	return;
}
;

void Particle::SetMTDir(Double_t m, Double_t T, TVector3 dir) {
	// dir: set direction vector of impulse
	// m: set mass in MeV
	// T: set kinetic energy in MeV
	// imp: need to use TVector3 dir(dir(0),dir(1),dir(2))
	// where dir(0), dir(1), dir(2) are coordinates of direction vector of impulse

	Double_t p = CalcP(m, T); // impulse magnitude
	fMass = m;
	fImpulse.SetE(T + m);
	fImpulse.SetPx(dir.Unit().x() * p);
	fImpulse.SetPy(dir.Unit().y() * p);
	fImpulse.SetPz(dir.Unit().z() * p);
}

void Particle::SetP(TVector3 P) {
	// set impulse of the particle
	fImpulse.SetXYZM(P(0), P(1), P(2), fMass);
}

void Particle::SetPx(Double_t px) {
// set cartesian x coordinate of impulse in MeV
	fImpulse.SetPx(px);
	//fImpulse.SetE(0.);
	fImpulse.SetE(Sqrt(fMass * fMass + fImpulse.Mag() * fImpulse.Mag()));
	return;
}
;

void Particle::SetPy(Double_t py) {
// set cartesian y coordinate of impulse in MeV
	fImpulse.SetPy(py);
	//fImpulse.SetE(0.);
	fImpulse.SetE(Sqrt(fMass * fMass + fImpulse.Mag() * fImpulse.Mag()));
}
;
void Particle::SetPz(Double_t pz) {
// set cartesian z coordinate of impulse in MeV
	fImpulse.SetPz(pz);
	fImpulse.SetE(0.);
	fImpulse.SetE(Sqrt(fMass * fMass + fImpulse.Mag() * fImpulse.Mag()));
}
;
void Particle::SetE(Double_t E) {
	// set total energy in MeV
	fImpulse.SetE(E);
}
;

void Particle::SetObservable(Bool_t obs) {
//whether object is observable or not
	fObservable = obs;
}
;

void Particle::SetMass(Double_t mass) {
	// Set particle mass.
	// Impulse and Kinetic energy are remained
	Double_t T = GetT();
	fMass = mass;
	fImpulse.SetE(T + mass);
}

void Particle::SetT(Double_t T) {
//set kinetic energy in MeV
//NB! impulse should be nonzero else you'll receive message "zero vector can't be stretched"

	if (T <= 0)
		fImpulse = TLorentzVector(0, 0, 0, fMass);
	else {
		fImpulse.SetRho(CalcP(fMass, T));
		fImpulse.SetE(fMass + T);
	}
}

void Particle::SetImpulse(TLorentzVector *P) {
	//set TLorenzVector
		fImpulse.SetPx(P->Px()*1000.);
		fImpulse.SetPy(P->Py()*1000.);
		fImpulse.SetPz(P->Pz()*1000.);
		fImpulse.SetE(P->E()*1000.);
}

void Particle::SetTPhiTheta(Double_t T, Double_t phi, Double_t theta) {
//change kinetic energy and direction
//T: kinetic energy in MeV
//phi: azimuthal angle of impulse in rad
//theta: polar angle of impulse in rad
//NB! impulse should be nonzero else you'll receive message "zero vector can't be stretched"

	fImpulse.SetRho(CalcP(fMass, T));
	fImpulse.SetE(fMass + T);
	fImpulse.SetPhi(phi);
	fImpulse.SetTheta(theta);
}
;

void Particle::Print(Option_t * option){
	// Print values:
	//		particle name, whether particle is observable,
	//		mass number and charge (in units of electron charge) of the particle,
	//		rest-mass and total energy in MeV, kinetic energy in MeV,
	//		coordinates of impulse and value of impulse in MeV,
	//		phi and theta angles in rad
	if (option != NULL) {
		//just checking
	}
	cout << "Name: " << fName;
	if (fObservable == 0)
		cout << " is not observable" << endl;
	else
		cout << " is observable" << endl;
	cout << "A, Z = " << fA << ",  " << fZ << endl << "Mass = " << fMass << endl
			<< "Energy = " << fImpulse.E() << endl << "Kinetic Energy = "
			<< GetT() << endl << "Impulse (Px;Py;Pz) = " << "(" << fImpulse.Px()
			<< " " << fImpulse.Py() << " " << fImpulse.Pz() << ")" << endl
			<< "P = " << fImpulse.Rho() << endl << "Phi = " << fImpulse.Phi()
			<< endl << "Theta = " << fImpulse.Theta() << endl;
}
;

TString Particle::CreateConfigString() {
	//Creates string containing config info about particle
	//so it can be saved to file or used anywhere else (gui?)
	//It saves direction of impulse as phi and theta angles
	//Its format looks like: (for example proton) #mass is stored in MeV I presume?
	//"name"="proton" "mass"="938" "observable"="true" "a_number"="1" "z_number"="1" "phi"="0.0" "theta"="0.0"

	ConfigDictionary CD;
	CD.SetString("name", GetName());
	CD.SetDouble("mass", GetM());
	CD.SetBool("observable", fObservable);
	CD.SetInt("a_number", GetA());
	CD.SetInt("z_number", GetZ());
	CD.SetDouble("phi", GetPhi());
	CD.SetDouble("theta", GetTheta());
	CD.SetDouble("k_energy", GetT()); //kinetic energy = GetT()
	CD.SetInt("exNumber", GetNumberOfExStates());
	TString ret(CD.ToString());

	return ret;
}

void Particle::ReadConfigString(TString conf) {
	//conf - string formatted as list of "key=value" pairs
	//Reads and sets parameters from config string
	//Parameters for impulse direction are phi and theta angles

	ConfigDictionary CD(conf.Data());
	try {
		SetName(CD.GetString("name").c_str());
		SetMass(CD.GetDouble("mass"));
		fGroundStateMass = GetM();
		Info("Particle::ReadConfigString", "Ground state mass was set to %f",
				fGroundStateMass);
		SetObservable(CD.GetBool("observable"));
		SetA(CD.GetInt("a_number"));
		SetZ(CD.GetInt("z_number"));
		SetPz(1.);
		SetTPhiTheta(CD.GetDouble("k_energy"), CD.GetDouble("phi"),
				CD.GetDouble("theta"));
	} catch (...) {
		Error("Particle::ReadConfigString", "Couldn't find all parameters!");
		return;
	}
}

TVector3 Particle::GetBoost() {
	// return vector beta (impulse components divided by the time component)
	return fImpulse.BoostVector();
}
;

void Particle::BoostTransform(TVector3 beta) {
//perform a boost transformation
//to correct working,we input (-beta), see TLorentzVector documentation
//this bug will be corrected soon
	fImpulse.Boost(-beta);
	return;
}
;

void Particle::CreateStatesWeigthsHist() {

	fStatesWeigths.SetBins(GetNumberOfExStates(), 0, GetNumberOfExStates());
//	fStatesWeigths.SetBins(8, 1, 8+1);

	for (Int_t i = 0; i < GetNumberOfExStates(); i++) {
//		cout << "\t\tjhagvdjhavsdasvd" << endl;
//		cout << i << endl;
//		cout << GetExcitationState(i).GetStrength() << endl;
//		cout << fStatesWeigths.GetNbinsX() << endl;
		fStatesWeigths.AddBinContent(i + 1,
				GetExcitationState(i).GetStrength());
	}
}

void Particle::DrawWeigths() {
	fStatesWeigths.Draw();
	return;
}

void Particle::DrawMassDistribution(Int_t i, Option_t *option) {
	//"i" is number of excitation state
	//"option" is draw option in TF1
	if (fState.empty()) {
		Error("Particle::DrawMassDistribution",
				"There are no defined states for particle %s", GetName());
		return;
	}

	if (i > (Int_t) fState.size()) {
		Error("Particle::DrawMassDistribution", "Maximum possible index is %d",
				(Int_t) (fState.size()));
		return;
	}

	fState[i]->Draw(option);
	return;
}

void Particle::CreateStateMassFunctions() {
	//method creates TF1 function for each state as defined e.g.
	//in configuration file
	//
	//Two distributions are realized for the present in the code:
	//Gauss and Lorentz distributions.
	//In future releases other distributions will be added.
	//
	//In order to optimize the calculation process there exist
	//limits of distribution:
	//(-3*sigma ; 4*sigma) for Gauss;
	//(-3*FWHM ; 4*FWHM) for Lorentz.
	//This limits provide 98% of integral area.

	CreateStatesWeigthsHist();

	TF1 *massSpectrum;
	massSpectrum = NULL;
	Double_t xmin = 0.; //lower and
	Double_t xmax = 0.; //upper range for each state

	for (Int_t i = 0; i < GetNumberOfExStates(); i++) {
		// cout << GetNumberOfExStates()<< " WTF!!!!! " << GetName() << endl; 
		//gauss distribution of the excited state
		if (GetExcitationState(i).GetShape() == "gauss") {

			const Double_t mean = GetExcitationState(i).GetMean();
			const Double_t sigma = GetExcitationState(i).GetWidth() / 2.35;

			Info(
					"Particle::CreateStateMassFunctions",
					"Creating gauss distribution with mean %f and FWHM of %f MeV",
					mean, sigma);

//			const Double_t mean = GetExcitationState(i).GetMean();
//			const Double_t sigma = GetExcitationState(i).GetWidth()/2.35;
			xmin = mean - 4 * sigma; // important cuz for high mean and little sigma it doesnt work cuz of fNpx (see TRandom)
			if (mean - 3 * sigma < 0. || sigma < 0.) {
				Warning("Particle::CreateStateMassFunctions",
						"There is something weird in this state");
				xmin = mean - 3 * sigma;
			} //if

			xmax = mean + 4 * sigma;
			Info("Particle::CreateStateMassFunctions",
					"Range for this state was set from %f to %f", xmin, xmax);
			//	(-3*sigma, +4*sigma)
			massSpectrum = new TF1("massGauss", "TMath::Gaus(x, [0], [1])",
					xmin, xmax);
			massSpectrum->SetParameters(mean, sigma); //1st parameter is mean value
													  //2nd parameter is sigma
		} //if gauss

		if (GetExcitationState(i).GetShape() == "lorentzian") {

			const Double_t mean = GetExcitationState(i).GetMean();
			const Double_t fwhm = GetExcitationState(i).GetWidth();

			Info(
					"Particle::CreateStateMassFunctions",
					"Creating lorentz distribution with mean %f and FWHM of %f MeV",
					mean, fwhm);

			//			const Double_t mean = GetExcitationState(i).GetMean();
			//			const Double_t sigma = GetExcitationState(i).GetWidth()/2.35;

			if (-3 * fwhm < 0.) {
				Warning("Particle::CreateStateMassFunctions",
						"There is something weird in this state");
				xmin = mean - 3 * fwhm;
			} //if

			xmax = mean + 4 * fwhm;
			Info("Particle::CreateStateMassFunctions",
					"Range for this state was set from %f to %f", xmin, xmax);
			//	(-3*sigma, +4*sigma)
			massSpectrum = new TF1("massLorentz",
					"TMath::BreitWigner(x, [0], [1])", xmin, xmax);
			massSpectrum->SetParameters(mean, fwhm); //1st parameter is mean value

		} //if lorentz

		if (massSpectrum) {
			fState.push_back(massSpectrum);
			Info("Particle::CreateStateMassFunctions",
					"Function %s was created", fState[i]->GetName());
		} else {
			Warning("Particle::CreateStateMassFunctions",
					"Function for e.s. number %d of type %s was not created", i,
					GetExcitationState(i).GetShape().Data());
		}
	}

	return;
}

void Particle::GenerateMass() {
	//calculate and set mass according to mass distribution functions

	if (!GetNumberOfExStates())
		return;

	//choose state:
	Int_t state = (Int_t) fStatesWeigths.GetRandom();
	//fixme Vratislav: random mass
	//for a few reaction in chain the same random numbers are obtained in each event
	//same situation as for angles
	Double_t mass = fGroundStateMass + fState[state]->GetRandom();
	SetMass(mass);

	return;
}

Bool_t Particle::IsObservable() {
//whether particle is observable (1) or not (0)
	return fObservable;
	/*
	 if(fObservable==0)
	 {
	 cout << fName << " is not observable" << endl;
	 return fObservable;
	 }
	 else
	 {
	 cout << fName << " is observable" << endl;
	 return fObservable;
	 }*/
}

//private functions
Double_t Particle::CalcT(Double_t m, Double_t px, Double_t py, Double_t pz) {
// Calculate kinetic energy
	TVector3 p(px, py, pz);
	return (Sqrt(p.Mag2() + m * m) - m);
}
;

Double_t Particle::CalcP(Double_t m, Double_t T) {
// Calculate value of impulse
	return (Sqrt(T * T + 2 * T * m));
}
;

Double_t Particle::CalcE() {
	return GetE();
}
;

Bool_t Particle::CheckEnergyConservation() {
// Check energy conservation
	Double_t ls = Power(fImpulse.Energy(), 2);
	Double_t ps = Power(fImpulse.Rho(), 2) + fMass * fMass;
	printf("%2.18f\t\n", ls);
	printf("%2.18f\t\n", ps);
	printf("%2.18f\t\n", TMath::Abs(ls / ps - 1.0));
	printf("%2.18f\t\n", 1.0e-15);
	printf("%f\t%f\n", Power(fImpulse.Rho(), 2), fMass * fMass);

	if (TMath::Abs(ls / ps - 1.0) < 1.0e-14) {
		cout << "Energy is conserved" << endl;
		return true;
	} else {
		cout << "Energy's not conserved" << endl;
		return false;
	}
}