/*
 * TwoParticleState.cpp
 *
 *  Created on: Jan 11, 2012
 *      Author: Anastasia Karavdina
 */

#include "TwoParticleState.h"
#include "TMath.h"
#include "TLorentzVector.h"
#include "TLorentzRotation.h"
//#include "TVector3.h"

/*TwoParticleState::TwoParticleState(double p, double theta, double phi,
							Particle p1, Particle p2, Particle res){
	_p = p;
	_theta = theta;
	_phi = phi;
	_s1 = p1.Spin();
	_lya1 = p1.Lya();
	_s2 = p2.Spin();
	_lya2 = p2.Lya();
	_J = res.Spin();
	_Jm = res.SpinM();
	AngMomCalc(p1, p2, res);
}*/

TwoParticleState::TwoParticleState(const TwoParticleState& st){
	_p = st._p;
	_theta = st._theta;//theta angle of 1st particle
    _phi = st._phi;//phi angle of 1st particle
	_s1 = st._s1;//spin of 1st particle

	/*_lya1.erase(_lya1.begin(),_lya1.end());
	for(int i=0;i<st.Lya1().size();i++)
		_lya1.push_back(st.Lya1()[i]);*/

	/*_lya1.erase(_lya1.begin(),_lya1.end());
	for(int i=0;i<st._lya1.size();i++)
		_lya1.push_back(st._lya1[i]);//helicity of 1st particle*/
	_lya1 = st._lya1;
	_s2 = st._s2;//spin of 2nd particle
	//_lya2 = st.Lya2();//helicity of 2nd particle
	//_lya2 = st._lya2;
	/*_lya2.erase(_lya2.begin(),_lya2.end());
		for(int i=0;i<st._lya2.size();i++)
			_lya2.push_back(st._lya2[i]);//helicity of 2nd particle*/
		_lya2 = st._lya2;
	_J = st._J;//Total moment of system
	/*_Jm = st.Jm();//Total moment projection on z-axis*/
	_Jm = st._Jm;//Total moment projection on z-axis
	//_l = st.AngMom();//angular moment
	_l = st._l;
	_mR = st._mR;//resonance mass
	_wR = st._wR;//resonance width
	/*_codeR = st.CodeR();
	_code1 = st.Code1();
	_code2 = st.Code2();*/
	_codeR = st._codeR;
	_code1 = st._code1;
	_code2 = st._code2;
	_m1 = st._m1;
	_m2 = st._m2;
	_p1 = st._p1;
	_p2 = st._p2;
	_res = st._res;

	_lfp1 = st._lfp1;
	_lfp2 = st._lfp2;
	/*_lfp1.erase(_lfp1.begin(),_lfp1.end());
		for(int i=0;i<st._lfp1.size();i++)
			_lfp1.push_back(st._lfp1[i]);
	_lfp2.erase(_lfp2.begin(),_lfp2.end());
	for(int i=0;i<st._lfp2.size();i++)
		_lfp2.push_back(st._lfp2[i]);*/
	//cout<<"size of _lfp1="<<_lfp1.size()<<endl;
//	cout<<"size of _lfp2="<<_lfp2.size()<<endl;
}

///Particle p1, p2 and resonance res
TwoParticleState::TwoParticleState(Particle p1, Particle p2, ResonanceP res){
	_p = 0;
	_theta = 0;
    _phi = 0;
	_mR = res.Mass();
	_wR = res.Width();
	_s1 = p1.Spin();
	_lya1 = p1.Lya();
	_s2 = p2.Spin();
	_lya2 = p2.Lya();
	_J = res.Spin();
	_Jm = res.SpinM();
	_codeR = res.Code();
	_code1 = p1.Code();
	_code2 = p2.Code();
	_m1 = p1.Mass();
	_m2 = p2.Mass();
	AngMomCalc(p1, p2, res);
	_p1 = p1;
	_p2 = p2;
	_res = res;
}

///Particle p1 with momentum mom1 (in lab frame) and so on.
/// recalculate momentum in the rest frame of resonance
TwoParticleState::TwoParticleState(Particle p1, TVector3 mom1,
								   Particle p2, TVector3 mom2,
								   ResonanceP res){
	_mR = res.Mass();
	_wR = res.Width();
	ResRestFrame(mom1, mom2);
	_s1 = p1.Spin();
	_lya1 = p1.Lya();
	_s2 = p2.Spin();
	_lya2 = p2.Lya();
	_J = res.Spin();
	_Jm = res.SpinM();
	AngMomCalc(p1, p2, res);
	_m1 = p1.Mass();
	_m2 = p2.Mass();
	_p1 = p1;
	_p2 = p2;
	_res = res;

}
TwoParticleState::~TwoParticleState() {
}

void TwoParticleState::ResRestFrame(TVector3 mom1, TVector3 mom2){
	double E1 = TMath::Sqrt(_m1*_m1+mom1.Mag2());
	TLorentzVector p1l(mom1,E1);

	double E2 = TMath::Sqrt(_m2*_m2+mom2.Mag2());
	TLorentzVector p2l(mom2,E2);

	TVector3 respv = (mom1+mom2);
	double Er = E1+E2;
	TLorentzVector resl(respv,Er);

	TLorentzRotation lt;//Lorentz transformation to rest system of resonance
	lt.RotateZ(-(respv.Phi()));
	lt.RotateY((TMath::Pi()-respv.Theta()));//rotate to the rest frame
	resl.Transform(lt);
	lt.Boost(-resl.BoostVector());
	resl.Boost(-resl.BoostVector());
	p1l.Transform(lt);
	p2l.Transform(lt);
	TVector3 p1lnew = p1l.Vect();
	TVector3 p2lnew = p2l.Vect();
	TVector3 p12lnew= p1lnew+p2lnew;
	 if(fabs(p1lnew.Mag()-p2lnew.Mag())>1e-6){
		 cout<<"Attention! Something went wrong during boost to the rest frame of resonance!"<<endl;
		 cout<<"p1new="<<p1lnew.Mag()<<" p2lnew="<<p2lnew.Mag()<<endl;
		 cout<<"p12lnew="<<p12lnew.Mag()<<endl;
	 }
	/* else{
		 cout<<"_p = "<<p1lnew.Mag()<<" _theta = "<<p1lnew.Theta()<<" _phi = "<<p1lnew.Phi()<<endl;
	 }*/
     _p = p1lnew.Mag();
     _theta = p1lnew.Theta();
     _phi =  p1lnew.Phi();
     if(_phi<0) _phi+=TMath::Pi();
    // cout<<"_p = "<<_p<<" _theta = "<<_theta<<" _phi = "<<_phi<<endl;
}


void TwoParticleState::AngMomCalc(Particle p1, Particle p2, Particle res){
	_l.erase(_l.begin(),_l.end());
	  double J = res.Spin();
	  double s1 = p1.Spin();
	  double s2 = p2.Spin();
	  double maxl = 2*J;
	  if(J<(s1+s2)) maxl = 2*(s1+s2);
	  int Pm = res.P();
	  int Pd1 = p1.P();
	  int Pd2 = p2.P();
	//  cout<<"J="<<J<<" s1="<<s1<<" s2="<<s2<<endl;
	 // cout<<"Pm="<<Pm<<" Pd1="<<Pd1<<" Pd2="<<Pd2<<endl;
	  for(int il=0;il<=int(maxl);il++){
	    double l = il;
	    int Pres = Pd1*Pd2*TMath::Power(-1,il);
	    if(Pres!=Pm) continue;
	//    cout<<"fabs(s1+s2-l) = "<<fabs(s1+s2-l)<<" fabs(s1+s2+l) = "<<fabs(s1+s2+l)<<endl;
	    if(fabs(s1+s2-l)<=J && fabs(s1+s2+l)>=J){
	    	//cout<<"add l = "<<l<<endl;
	    	_l.push_back(l);
	    }
	  }
//	  for(int i=0;i<_l.size();i++)
	//	  cout<<"_l["<<i<<"]="<<_l[i]<<endl;
}
vector<double> TwoParticleState::AngMom(){
	AngMomCalc(_p1, _p2, _res);
	return _l;}
double TwoParticleState::Mom(){return _p;}
double TwoParticleState::Theta(){return _theta;}
double TwoParticleState::Phi(){return _phi;}
double TwoParticleState::S1(){return _s1;}
double TwoParticleState::S2(){return _s2;}
double TwoParticleState::J(){return _J;}
double TwoParticleState::mR(){return _mR;}
double TwoParticleState::wR(){return _wR;}
vector<double> TwoParticleState::Lya1() {return _lya1;}
vector<double> TwoParticleState::Lya2() {return _lya2;}
vector<double> TwoParticleState::Jm() {return _Jm;}
int TwoParticleState::CodeR(){return _codeR;}
int TwoParticleState::Code1(){return _code1;}
int TwoParticleState::Code2(){return _code2;}
double TwoParticleState::mass1(){return _m1;}
double TwoParticleState::mass2(){return _m2;}
pair<Particle,Particle> TwoParticleState::GetFP(){//get pair of final particles
	pair<Particle,Particle> fppair;
	fppair.first = _p1;
	fppair.second = _p2;
	return fppair;
}
ResonanceP TwoParticleState::GetIP(){//get initial particles
	return _res;
}

 void  TwoParticleState::setListFP1(vector<int> lfp){
	  _lfp1.erase(_lfp1.begin(),_lfp1.end());
	  for(int i=0;i<lfp.size();i++)
	  _lfp1.push_back(lfp[i]);
  }

  void  TwoParticleState::setListFP2(vector<int> lfp){
	  _lfp2.erase(_lfp2.begin(),_lfp2.end());
	  for(int i=0;i<lfp.size();i++)
		  _lfp2.push_back(lfp[i]);
  }