BeReaction.cpp

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

#include "BeReaction.h"

ClassImp(BeReaction);

Int_t ranseed = 1;
//TRandom3	BeReaction::ranTheta(1277372118);
//TRandom3	BeReaction::ranTheta(0);
//TRandom3	BeReaction::ranMass(1277372118);
TRandom3	BeReaction::ranMass(0);
TRandom3	BeReaction::ranTheta(0);
Double_t BeReaction::LipM[] = {0};
Double_t BeReaction::BeDecayM[] = {0};
Double_t BeReaction::BeStateEnergy[] = {0};
Double_t BeReaction::DipDecayM[] = {0};
TF1 BeReaction::theta("thetaCM", "TMath::Sin(x)", 0., TMath::Pi());
TF1 BeReaction::thetaUniform("thetaCM", "1", 0., TMath::Pi());
TF1 BeReaction::dipEnergy("dpexcitation", "TMath::Sqrt( x*([0] - x) )", 0., 1.);

BeReaction::BeReaction()
{
	ReadParameters();
//	PrintParameters();

}

BeReaction::~BeReaction() {

}

void BeReaction::ReadParameters(const char *parameterfile)
{
	ifstream parfile;

	parfile.open(parameterfile);												//parameter file opening
	if (!parfile.is_open()) {
		Warning("BeReaction::ReadParameters", "File %s opening error\n", parameterfile);
		return;
	}//if

	Double_t value;
	TString line;

	while (parfile.good()) {
		parfile >> value;
		line.ReadLine(parfile);
		if (line.Contains("//")) line.Resize(line.Index("//"));
		line.ToLower();
		if (line.Contains("6li")) {
			LipM[0] = value;
			continue;
		}
		if (line.Contains("proton")) {
			LipM[1] = value;
			BeDecayM[2] = 2* value;
			DipDecayM[2] = value;
			DipDecayM[3] = value;
			continue;
		}
		if (line.Contains("neutron")) {
			LipM[3] = value;
			continue;
		}
		if (line.Contains("4he") || line.Contains("alpha")) {
			BeDecayM[3] = value;
			continue;
		}
		if (line.Contains("6be_g.s._e")) {
			BeStateEnergy[0] = value;
			continue;
		}
		if (line.Contains("6be_g.s._gamma")) {
			BeStateEnergy[1] = value;
			continue;
		}
		if (line.Contains("6be_e.s._e")) {
			BeStateEnergy[2] = value;
			continue;
		}
		if (line.Contains("6be_e.s._gamma")) {
			BeStateEnergy[3] = value;
			continue;
		}
		if (line.Contains("6be_(g.s.)_prob")) {
			BeStateEnergy[4] = value;
			continue;
		}
		if (line.Contains("6be_3body_thr")) {
			BeStateEnergy[5] = value;
			continue;
		}
	}//while

	//set 6Be mass
	LipM[2] = BeDecayM[3] /*4He*/ + 2*DipDecayM[3] /*proton*/ - BeStateEnergy[5] /*threshold*/;
	BeDecayM[0] = BeDecayM[3] + 2*DipDecayM[3] - BeStateEnergy[5];

	parfile.close();
	if (parfile.is_open()) {												//parameter file closing
		Warning("BeReaction::ReadParameters", "File %s closing error\n", parameterfile);
		return;
	}

	Info("BeReaction::ReadParameters", "Mass and other parameters were set");
	return;
}

void BeReaction::FillProcess(Double_t _LiT, Double_t _LiThetaIn, Double_t _LiPhi)
{

	//simple MC generator using phase volume

	Double_t	operatingMarray[4] = {0};		//auxiliary array for mass setting
//	SetLipMasses(operatingMarray);				//to simulate in zero approximation
//	SetLipMasses_uniformBeMass(operatingMarray);	//to determine energy resolution
	SetLipMasses_discreteBeMass(operatingMarray);	//to determine energy resolution
	fLip.SetName("fLip");
//	fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, ThetaCMdistr(0., TMath::Pi()), _LiPhi);
//	fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, ThetaCMdistrUniform(0., TMath::Pi()), _LiPhi);
	fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, ThetaCMdistrDiscrete(), _LiPhi);
	SetBeDecayMasses(operatingMarray);
	fBeDecay.SetName("fBeDecay");
	fBeDecay.FillReaction(fLip.GetThetaA(), fLip.GetTa(), operatingMarray, ThetaCMdistr(0., TMath::Pi()), fLip.GetPhiA() );
	SetDipDecayMasses(operatingMarray);
	fDipDecay.SetName("fDipDecay");
	fcorrect = fDipDecay.FillReaction( fBeDecay.GetThetaA(), fBeDecay.GetTa(), operatingMarray, ThetaCMdistr(0., TMath::Pi()), fBeDecay.GetPhiA() );
	if (fcorrect != 1) {
		cout << operatingMarray[0] - 2*DipDecayM[2] << endl;
	}

}

void BeReaction::FillProcess(Double_t _LiT, Double_t _LiThetaIn,
		Double_t _LiPhi, Double_t *_p_alpha, Double_t *_p_p1, Double_t *_p_p2,
		Double_t _thetaCMmin, Double_t _thetaCMmax)
{

	Double_t thetaCM = ThetaCMdistr(_thetaCMmin, _thetaCMmax);
	FillProcess(_LiT, _LiThetaIn, _LiPhi, _p_alpha, _p_p1, _p_p2, thetaCM);
}

void BeReaction::FillProcess(Double_t _LiT, Double_t _LiThetaIn,
		Double_t _LiPhi, Double_t *_p_alpha, Double_t *_p_p1, Double_t *_p_p2,
		Double_t _thetaCM)
{
	//function variables: LiT, LiThetaIn, LiPhiIn, BeThetaCM, impulses array
	//_LiT:
	//_thetaCMmin: in rad
	//_thetaCMmax: in rad


	//physical generator based on theoretical calculations

		//set masses for 6Li + p --> 6Be + n: mass of 6Be taken from generator
	TLorentzVector pa(_p_alpha[0], _p_alpha[1], _p_alpha[2]);	//alpha in the a-p-p CM
	pa.SetE( E( Power(pa.P(), 2), BeDecayM[3] ) );
	TLorentzVector pp1(_p_p1[0], _p_p1[1], _p_p1[2]);			//proton in the a-p-p CM
	pp1.SetE( E( Power(pp1.P(), 2), LipM[1] ) );
	TLorentzVector pp2(_p_p2[0], _p_p2[1], _p_p2[2]);			//proton in the a-p-p CM
	pp2.SetE( E( Power(pp2.P(), 2), LipM[1] ) );
	TLorentzVector pbe;											//6Be in the a-p-p CM
	pbe = pa + pp1 + pp2;

	Double_t	operatingMarray[4] = {0};		//auxilliary array for mass setting
	operatingMarray[0] = LipM[0];
	operatingMarray[1] = LipM[1];
	operatingMarray[2] = pbe.E();
	operatingMarray[3] = LipM[3];

		//reaction 6Li + p --> 6Be + n: BeThetaCM taken from generator
	fLip.SetName("fLip");
//	fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, _thetaCM*TMath::DegToRad(), _LiPhi);
	fLip.FillReaction(_LiThetaIn, _LiT, operatingMarray, _thetaCM/*ThetaCMdistr(_thetaCMmin, _thetaCMmax)*/, _LiPhi);	//thetaCM from generator is not used

	//transformation from a-p-p CM system into lab system
	TRotation r1;
	r1.SetZAxis(-GetNeutronP(), GetBeP());
	TLorentzRotation rot(r1);

	//rotation to CM system where axis are parallel with lab system (CMaux)
	pa.Transform(rot);
	pp1.Transform(rot);
	pp2.Transform(rot);

	TVector3 beta = GetBe().BoostVector();

	//CMaux --> lab
	pa.Boost(beta);
	pp1.Boost(beta);
	pp2.Boost(beta);


	//set the variables fBeDecay.fTb(.ThetaA, .PhiA) (alpha),
	//					fDipDecay.fTa(.ThetaA, .PhiA) (one of the protons),
	//					fDipDecay.fTb(.ThetaB, .PhiB) (one of the protons)
		//set variables for the alpha particle
	fBeDecay.SetM(3, BeDecayM[3]);
	fBeDecay.SetPhiB(pa.Phi());
	fBeDecay.SetThetaB(pa.Theta());
	fBeDecay.SetTb(pa.E()-BeDecayM[3]);

		//set variables for the first proton
	fDipDecay.SetM(2, DipDecayM[2]);
	fDipDecay.SetPhiA(pp1.Phi());
	fDipDecay.SetThetaA(pp1.Theta());
	fDipDecay.SetTa(pp1.E()-DipDecayM[2]);

		//set variables for the second proton
	fDipDecay.SetM(3, DipDecayM[3]);
	fDipDecay.SetPhiB(pp2.Phi());
	fDipDecay.SetThetaB(pp2.Theta());
	fDipDecay.SetTb(pp2.E()-DipDecayM[3]);
//	Info("BeReaction::FillProcess", "%f\t%f\t%f", pp2.E(), pp2.Phi(), pp2.Theta());

	return;

}

void BeReaction::PrintParameters()
{
	cout <<endl;
	for (Int_t i = 0; i < 4; i++) {
		cout << "LipM[" << i << "] is " << LipM[i] << endl;
	}

	cout <<endl;
	for (Int_t i = 0; i < 4; i++) {
		cout << "BeDecayM[" << i << "] is " << BeDecayM[i] << endl;
	}

	cout <<endl;
	for (Int_t i = 0; i < 6; i++) {
		cout << "BeStateEnergy[" << i << "] is " << BeStateEnergy[i] << endl;
	}

	cout <<endl;
	for (Int_t i = 0; i < 4; i++) {
		cout << "DipDecayM[" << i << "] is " << DipDecayM[i] << endl;
	}

	return;
}

Double_t BeReaction::ThetaCMdistr(Double_t tmin, Double_t tmax) const {

	return theta.GetRandom(tmin, tmax);

}

Double_t BeReaction::ThetaCMdistrUniform(Double_t tmin, Double_t tmax) const {

	return thetaUniform.GetRandom(tmin, tmax);

}

Double_t BeReaction::ThetaCMdistrDiscrete() const {

	const Int_t nopoints = 6;
	const Double_t angles[nopoints] = {20*TMath::DegToRad(),
										45*TMath::DegToRad(),
										70*TMath::DegToRad(),
										95*TMath::DegToRad(),
										120*TMath::DegToRad(),
										150*TMath::DegToRad()};

	Int_t angle = ranMass.Integer(nopoints);
	return angles[angle];
}

void BeReaction::SetLipMasses(Double_t _lpm[]) {

	for (Int_t i = 0; i < 4; i++) _lpm[i] = 0.;

	//constant masses setting
	_lpm[0] = LipM[0];		//mass of 6Li
	_lpm[1] = LipM[1];		//mass of proton
	_lpm[3] = LipM[3];		//mass of neutron

	//variable 6Be mass setting
	if ( BeStateEnergy[4] >= ranMass.Uniform(0., 1.) ) {
		while (_lpm[2] <= BeDecayM[3] + 2*DipDecayM[3]) {
			_lpm[2] = LipM[2] + ranMass.Gaus(BeStateEnergy[0], BeStateEnergy[1]);
//			printf("%f\t>\t%f\n",  BeStateEnergy[0], BeStateEnergy[1]);
//			printf("%f\t>\t%f\n",  _lpm[2], fBeDecay.GetMb()+2*fDipDecay.GetMa());
			//		_lpm[2] = BeDecayM[3] + 2*DipDecayM[4] + ranMass.Gaus(BeStateEnergy[0], BeStateEnergy[1]);
		}//while
	}//if
	else {
		while (_lpm[2] <= BeDecayM[3] + 2*DipDecayM[3]) {
			_lpm[2] = LipM[2] + ranMass.Gaus(BeStateEnergy[2], BeStateEnergy[3]);
//			printf("%f\t>\t%f\n",  _lpm[2], fBeDecay.GetMb()+2*fDipDecay.GetMa());
	//		_lpm[2] = BeDecayM[3] + 2*DipDecayM[4] + ranMass.Gaus(BeStateEnergy[2], BeStateEnergy[3]);
		}//while
	}//if

//	printf("%f\t%f\t%f\t%f\n", LipM[0], LipM[1], LipM[2], LipM[3]);
	return;
}

void BeReaction::SetLipMasses_uniformBeMass(Double_t _lpm[]) {

	//constant masses setting
	_lpm[0] = LipM[0];		//mass of 6Li
	_lpm[1] = LipM[1];		//mass of proton
	_lpm[3] = LipM[3];		//mass of neutron


	//uniform distrubution 6Be mass setting in range (Emin;Emax)
	const Double_t Emin = 0.;
	const Double_t Emax = 20.;
//	while (_lpm[2] <= fBeDecay.GetMb()+2*fDipDecay.GetMa()) {
		_lpm[2] = BeDecayM[3] + 2*DipDecayM[3] + ranMass.Uniform(Emin, Emax);
//		_lpm[2] = LipM[2] + BeStateEnergy[5] + ranMass.Uniform(Emin, Emax);
//	}

	return;
}

void BeReaction::SetLipMasses_discreteBeMass(Double_t _lpm[]) {

	//constant masses setting
	_lpm[0] = LipM[0];		//mass of 6Li
	_lpm[1] = LipM[1];		//mass of proton
	_lpm[3] = LipM[3];		//mass of neutron


	//discrete distrubution 6Be mass 1.4, 3, 6, 9, 13
	const Int_t nopoints = 5;
	const Double_t masses[nopoints] = {1.4, 3, 6, 9, 13};
	Int_t mass = ranMass.Integer(nopoints);

	_lpm[2] = BeDecayM[3] + 2*DipDecayM[3] + masses[mass];

	return;
}

void BeReaction::SetBeDecayMasses(Double_t _bdm[]) {

	//constant masses setting
	_bdm[1] = BeDecayM[1];		//zero mass
	_bdm[3] = BeDecayM[3];		//mass of 4He

	//variable 6Be mass setting
	_bdm[0] = fLip.GetMa();

	//variable diproton mass setting
	const Double_t freeEnergy = fLip.GetMa() - (BeDecayM[0] + BeStateEnergy[5]);
	//TF1 dipEnergy("dpexcitation", "TMath::Sqrt( x*([0] - x) )", 0., freeEnergy);
	dipEnergy.SetRange(0., freeEnergy);
	dipEnergy.SetParameter(0, freeEnergy);
	const Double_t random = dipEnergy.GetRandom();
	_bdm[2] = BeDecayM[2] + random;		//problem
//	if (BeDecayM[2] > _bdm[2]) {
//		printf("//////////////////\t%f\t%f\t%f\n", BeDecayM[2], random, freeEnergy);
//		printf("//////////////////\t%f\t%f\t%f\t%f\n", fLip.GetMa(), BeDecayM[0], fLip.GetMa() -BeDecayM[0], BeStateEnergy[5]);
//	}

}

void BeReaction::SetDipDecayMasses(Double_t _2pm[])
{

	if (fBeDecay.GetMa()/2. < 938.272) {
		printf("\n\n%f\t%f\n\n", fBeDecay.GetMa()/2., fBeDecay.GetMb());
	}
	_2pm[0] = fBeDecay.GetMa();
	DipDecayM[0] = _2pm[0];
	_2pm[1] = 0;
	_2pm[2] = DipDecayM[2];
	_2pm[3] = DipDecayM[3];
}

TVector3 BeReaction::GetAlphaP()
{
	TVector3 vAlpha(1., 1., 1.);
	Double_t pc = TMath::Sqrt(TMath::Power(GetAlphaT(), 2) + 2*GetAlphaT()*fBeDecay.GetMb());
	vAlpha.SetMagThetaPhi(pc, GetAlphaTheta(), GetAlphaPhi());

	return vAlpha;
}

TLorentzVector BeReaction::GetAlpha()
{
//	TLorentzVector *vAlpha = new TLorentzVector(GetAlphaP(), GetAlphaT()+fBeDecay.GetMb());
//	cout << "getalfa Zacatek" << endl;
//	cout << GetAlphaT()+fBeDecay.GetMb() << endl;
//	cout << fBeDecay.GetMb() << endl;
//	cout << GetAlphaT() << endl;
//	cout << "getalfa konec" << endl;
	return TLorentzVector(GetAlphaP(), GetAlphaT()+fBeDecay.GetMb());
}

TVector3 BeReaction::GetP1P()
{
	TVector3 vP1(1., 1., 1.);
	Double_t pc = TMath::Sqrt(TMath::Power(GetP1T(), 2) + 2*GetP1T()*fDipDecay.GetMa());
	vP1.SetMagThetaPhi(pc, GetP1Theta(), GetP1Phi());

	return vP1;
}

TLorentzVector BeReaction::GetP1()
{
//	TLorentzVector *vP1 = new TLorentzVector(*GetP1P(), GetP1T()+fDipDecay.GetMa());

	return TLorentzVector(GetP1P(), GetP1T()+fDipDecay.GetMa());
}

TVector3 BeReaction::GetP2P()
{
	TVector3 vP2(1., 1., 1.);
	Double_t pc = TMath::Sqrt(TMath::Power(GetP2T(), 2) + 2*GetP2T()*fDipDecay.GetMb());
	vP2.SetMagThetaPhi(pc, GetP2Theta(), GetP2Phi());

	return vP2;
}

TLorentzVector BeReaction::GetP2()
{
//	TLorentzVector *vP2 = new TLorentzVector(*GetP2P(), GetP2T()+fDipDecay.GetMb());
//	Info("BeReaction::GetP2", "\t%f\t%f\t%f", GetP2P().Mag(), GetP2T(), GetP2Theta());
	return TLorentzVector(GetP2P(), GetP2T()+fDipDecay.GetMb());
}

//TVector3* BeReaction::GetNeutronP()
TVector3 BeReaction::GetNeutronP()
{
//	TVector3 *vN = new TVector3(1., 1., 1.);
	TVector3 vN(1., 1., 1.);
	Double_t pc = TMath::Sqrt(TMath::Power(GetNeutronT(), 2) + 2*GetNeutronT()*fLip.GetMb());
//	vN->SetMagThetaPhi(pc, GetNeutronTheta(), GetNeutronPhi());
	vN.SetMagThetaPhi(pc, GetNeutronTheta(), GetNeutronPhi());

	return vN;
}

TLorentzVector BeReaction::GetNeutron()
{
//	TLorentzVector *vN = new TLorentzVector(*GetNeutronP(), GetNeutronT()+fLip.GetMb());
//	TLorentzVector vN = new TLorentzVector(GetNeutronP(), GetNeutronT()+fLip.GetMb());

//	return vN;
	return TLorentzVector(GetNeutronP(), GetNeutronT()+fLip.GetMb());
}

TVector3 BeReaction::GetBeP()
{
	TVector3 vBe(1., 1., 1.);
	Double_t pc = TMath::Sqrt(TMath::Power(GetBeT(), 2) + 2*GetBeT()*fLip.GetMa());
	vBe.SetMagThetaPhi(pc, GetBeTheta(), GetBePhi());

	return vBe;
}

TLorentzVector BeReaction::GetBe()
{
	return TLorentzVector(GetBeP(), GetBeT()+fLip.GetMa());
}

Double_t BeReaction::T(Double_t pc2, Double_t mc2)
{
	return TMath::Power(mc2*mc2 + pc2, 0.5) - mc2;
}

Double_t BeReaction::E(Double_t pc2, Double_t mc2)
{
	return TMath::Power(mc2*mc2 + pc2, 0.5);
}

void BeReaction::Reset()
{
	fLip.Reset();
	fBeDecay.Reset();
	fDipDecay.Reset();
	fcorrect = 0;

	return;
}