/*
 * BinaryReaction.cpp
 *
 *  Created on: 24.5.2010
 *      Author: Vratislav
 */

#include "BinaryReaction.h"

BinaryReaction::BinaryReaction() {

}

BinaryReaction::~BinaryReaction() {

}


ClassImp(BinaryReaction);

//TRandom3 *ranpokus = new TRandom3(1);

TRandom3 BinaryReaction::ranPhi(1);


Int_t BinaryReaction::FillReaction(Double_t _fThetaIn, Double_t _fT1, Double_t _fM[], Double_t _fThetaCM, Double_t _fPhiIn)
{
	//_fThetaIn in rad
	//_fT1 in MeV
	//_fM[] in MeV
	//_fThetaCM in rad
	//_fPhiIn in rad


	fThetaIn = _fThetaIn;
	fT1 = _fT1;
	for (Int_t i = 0; i < 4; i++) {
		fM[i] = _fM[i];
	}
	fThetaCM = _fThetaCM;
	fPhiIn = _fPhiIn;

	if (FlatCalculate() != 1) { return -1; }
	SetPhi();
	//SetOrientation();
	SetOrientationV();

	return 1;

}

Int_t BinaryReaction::FlatCalculate()
{
	const Double_t p1 = Sqrt(Power(fT1, 2) + 2*fT1*fM[0]);		//in MeV
	const Double_t E = fT1 + fM[0] + fM[1];						//in MeV

	//beta and gamma for Lor. trans.
	const Double_t beta = p1/E;
	const Double_t gamma = 1/Sqrt(1 - Power(beta, 2));

	//calculating of T3cm and T4cm
	//T3cm = ( Tcm^2 + 2*Tcm*m4*c^2 )/(2*( Tcm + m4*c^2 + m3*c^2 )), where Tcm = T3cm + T4cm
		//Ecm = Tcm + m3*c^2 + m4*c^2
			//Ecm^2 = E^2 - (p1*c)^2, where E = T1 + m1*c^2 + m2*c^2 and (p1*c)^2 = T1^2 + 2T1*m1*c^2

	Double_t Ecm = 0.;					//energy in CMS, MeV
	Double_t Tcm = 0.;					//summary kinetic energy in CMS, MeV
	Double_t TaCM = 0., TbCM = 0.;
	Double_t pcm = 0.;					//impuls in CMS, MeV

	Ecm = Sqrt( Power(E, 2) - Power(p1, 2) );
//	Tcm = Sqrt( Power(fM[0]+fM[1], 2) + 2*fT1*fM[1] ) - fM[2] - fM[3];
//	cout << Tcm << endl;
	Tcm = Ecm - fM[2] - fM[3];
//	cout << Tcm << endl;
	static Int_t counter = 0;
//	cout << Tcm << endl;

	if (Tcm < 0) {
		counter++;
		//cout << "\nThere is no solution, Tcm is smaller than threshold energy\n";
		Warning("FlatCalculate", "%s: No solution, Tcm is smaller than threshold E", GetName());
		cout << "Tcm is " << Tcm << endl
			<< "Ecm is " << Ecm << endl
			<< "fTa is " << fTa << endl
			<< "fTb is " << fTb << endl
			<< "fM[0] is " << fM[0] << endl
			<< "m_dp  is " << 2*fM[2] << endl
			<< "fT1 is " << fT1 << endl
			<< "E   is " << E << endl;
		printf("\n\n!!!!!!!!!\t%d\t!!!!!!!!!!!\n\n", counter);
		return -1;
	}

	TaCM = ( Power(Tcm, 2) + 2*Tcm*fM[3] )/( 2*(Tcm + fM[3] + fM[2]) );
	TbCM = ( Power(Tcm, 2) + 2*Tcm*fM[2] )/( 2*(Tcm + fM[3] + fM[2]) );

	//p3cm = -p4cm ==> |p3cm| = |p4cm| = pcm, (pcm*c)^2 = T3cm^2 + 2*T3cm*m3*c^2
	pcm = Sqrt( Power(TaCM, 2) + 2*TaCM*fM[2] );

	//Lorentz transformation of impuls from CMS to Lab:
	Double_t pax, pay, pbx, pby;		//in MeV
	Double_t EaCM, EbCM;				//in MeV
	Double_t paq = 0., pbq = 0.;		//in MeV

	//evaluate of T3
	EaCM = TaCM + fM[2];
	pay = pcm*Sin(fThetaCM);
	pax = gamma*(pcm*Cos(fThetaCM) + beta*EaCM);
	paq = Power(pax, 2) + Power(pay, 2);
	fTa = Sqrt( Power(fM[2], 2) + paq ) - fM[2];

	//evaluate of T4
	EbCM = TbCM + fM[3];
	pby = pcm*Sin(fThetaCM + TMath::Pi());
	pbx = gamma*(pcm*Cos(fThetaCM + TMath::Pi()) + beta*EbCM);
	pbq = Power(pbx, 2) + Power(pby, 2);
	fTb = Sqrt( Power(fM[3], 2) + pbq ) - fM[3];

	//calculating of theta3 and theta4
		//cotg(theta3) = 0 ==> theta3 = pi/2, cotg(alfa) = 1/tg(alfa)
		//tg(theta3) = ( pcm*sin(thetacm) )/( gamma*(pcm*cos(thetacm) + beta*Ecm/c) );
	if ( ( pcm*Sin(fThetaCM) )/( gamma*(pcm*Cos(fThetaCM) + beta*Sqrt(Power(pcm, 2) + Power(fM[2], 2))) ) > 10000000 ) {		//the upper member of the cotg(theta)
		fThetaA = TMath::Pi()/2;
	}
	else {
		fThetaA = ATan( (pcm*Sin(fThetaCM))/(gamma*(pcm*Cos(fThetaCM) + beta*Sqrt(Power(pcm, 2) + Power(fM[2], 2)))) );
		if (fThetaA < 0) {
			fThetaA = fThetaA + TMath::Pi();
		}
	}


	if ( ( pcm*Sin(TMath::Pi() - fThetaCM) )/( gamma*(pcm*Cos(TMath::Pi() - fThetaCM) + beta*Sqrt(Power(pcm, 2) + Power(fM[3], 2))) ) > 10000000 ) {		//the upper member of the cotg(theta)
		fThetaB = TMath::Pi()/2;
	}
	else {
		fThetaB = ATan( (pcm*Sin(TMath::Pi() - fThetaCM))/(gamma*(pcm*Cos(TMath::Pi() - fThetaCM) + beta*Sqrt(Power(pcm, 2) + Power(fM[3], 2)))) );
		if (fThetaB < 0) {
			fThetaB = fThetaB + TMath::Pi();
		}
	}

	return 1;

}

void BinaryReaction::SetPhi()
{
	//PhiA and PhiB assignment
	//fPhiA = ranpokus->Uniform(0., 2.*TMath::Pi());

	/*fPhiA = ranPhi.Uniform(0., 2.*TMath::Pi());
	if ( fPhiA < TMath::Pi() ) { fPhiB = fPhiA + TMath::Pi(); }
	else { fPhiB = fPhiA - TMath::Pi();	}*/

	fPhiA = ranPhi.Uniform(0., 2.*TMath::Pi());
	fPhiB = fPhiA - TMath::Pi();
}

void BinaryReaction::SetOrientation()
{
	if (fThetaA + fThetaIn <= TMath::Pi()) {
		fThetaA = fThetaA + fThetaIn;
	}
	else {
		fThetaA = 2*TMath::Pi() - fThetaA + fThetaIn;
		if ( fPhiA < TMath::Pi() ) { fPhiA = fPhiA + TMath::Pi(); }
		else { fPhiA = fPhiA - TMath::Pi(); }
	}
	if (fPhiA + fPhiIn > 2*TMath::Pi()) { fPhiA = fPhiA + fPhiIn - 2*TMath::Pi(); }
	else { fPhiA = fPhiA + fPhiIn; }

	if (fThetaB - fThetaIn >= 0) {
		fThetaB = fThetaB - fThetaIn;
	}
	else {
		fThetaB = -1*(fThetaB - fThetaIn);
		if ( fPhiB < TMath::Pi() ) { fPhiB = fPhiB + TMath::Pi(); }
		else { fPhiB = fPhiB - TMath::Pi(); }
	}
	if (fPhiB + fPhiIn > 2*TMath::Pi()) { fPhiB = fPhiB + fPhiIn - 2*TMath::Pi(); }
	else { fPhiB = fPhiB + fPhiIn; }

}

void BinaryReaction::SetOrientationV()
{
	//z' parallel to rdirection; x' is in theta_rdirection plane; y' perpendicular
	TVector3 rdirection(0,0,1);
	rdirection.SetTheta(fThetaIn);
	rdirection.SetPhi(fPhiIn);
	//rdirection.SetMag(1.);

//	TVector3	a(1);
	TVector3	a(1., 1., 1.);
	a.SetTheta(fThetaA);
	a.SetPhi(fPhiA);
	//a.SetMag(1.);

//	TVector3	b(1);
	TVector3	b(1., 1., 1.);
	b.SetTheta(fThetaB);
	b.SetPhi(fPhiB);
	//b.SetMag(1.);

	a.RotateUz(rdirection);
	b.RotateUz(rdirection);

	fThetaA = a.Theta();
	fPhiA = a.Phi();
	fThetaB = b.Theta();
	fPhiB = b.Phi();

}

void BinaryReaction::Reset()
{
	for (Int_t i = 0; i <= 3; i++) {
		fM[i] = 0.;
	}
	fThetaIn = 0.;
	fPhiIn = 0.;
	fT1 = 0.;
	fThetaCM = 0.;
	fTa = 0.;
	fTb = 0.;
	fThetaA = 0.;
	fThetaB = 0.;
	fPhiA = 0.;
	fPhiB = 0.;

	return;
}

Double_t BinaryReaction::GetXa(Double_t z)
{

	//univerzalni cast
	TVector3 v(1., 1., 1.);
	v.SetPhi(fPhiA);
	v.SetTheta(fThetaA);

//	pokud je Theta == pi/2 ==> nejaky jiny vystup
	if (v.CosTheta()*z <= 0) { return 0.; }

	v.SetMag(TMath::Abs(z/v.CosTheta()));

	return v.X();
}

Double_t BinaryReaction::GetYa(Double_t z)
{
	//univerzalni cast
	TVector3 v(1., 1., 1.);
	v.SetPhi(fPhiA);
	v.SetTheta(fThetaA);

//	pokud je Theta == pi/2 ==> nejaky jiny vystup
	if (v.CosTheta()*z <= 0) { return 0.; }

	v.SetMag(TMath::Abs(z/v.CosTheta()));

	return v.Y();
}

Double_t BinaryReaction::GetXb(Double_t z)
{

	//univerzalni cast
	TVector3 v(1., 1., 1.);
	v.SetPhi(fPhiB);
	v.SetTheta(fThetaB);

//	pokud je Theta == pi/2 ==> nejaky jiny vystup
	if (v.CosTheta()*z <= 0) { return 0.; }

	v.SetMag(TMath::Abs(z/v.CosTheta()));

	return v.X();
}

Double_t BinaryReaction::GetYb(Double_t z)
{
	//univerzalni cast
	TVector3 v(1., 1., 1.);
	v.SetPhi(fPhiB);
	v.SetTheta(fThetaB);

//	pokud je Theta == pi/2 ==> nejaky jiny vystup
	if (v.CosTheta()*z <= 0) { return 0.; }

	v.SetMag(TMath::Abs(z/v.CosTheta()));

	return v.Y();
}