/*
 * Amplitude.cpp
 *
 *  Created on: Jan 14, 2012
 *      Author: Anastasia Karavdina
 */

#include "Amplitude.h"
#include "qft++.h"
#include "TwoParticleState.h"
#include <complex>
#include <list>
#include <vector>
#include "TMath.h"
#include "Utils.h"
#include "TLorentzVector.h"
#include "TVector3.h"

Amplitude::Amplitude(Chain ch) {
	_ch = ch;
}

Amplitude::Amplitude(const Amplitude& a1){
	_ch = a1._ch;
}
Amplitude::~Amplitude() {
}

/// fill list of particles which built particle#i for TwoParticleState
/// lfpcode - list of f.p. codes
void Amplitude::Build(vector<int> lfpcode){
		list<TwoParticleState>  lstates = _ch.GetStates();
		const int stsize =  lstates.size();
		vector<bool> fbusy(stsize,false);
		int cpRpr = -1;//code of resonance from previous decay
		vector<int> posR;
		//int cpRprpr = -1;//code of resonance from previous-previous decay
		vector<int> posRpr;
	//	vector<int> usedpos;

	list<TwoParticleState> nlstates;

	for (list<TwoParticleState>::iterator it = lstates.begin(); it != lstates.end(); it++){
		//for (list<TwoParticleState>::iterator it = lstates.end(); it != lstates.begin(); --it){
	//	cout<<"usedpos.size()="<<usedpos.size()<<endl;
			posR.erase (posR.begin(),posR.end());
			TwoParticleState curst=(*it);
			int cp1= curst.Code1();
			int cp2= curst.Code2();
			int cpR= curst.CodeR();
		//	cout<<"cp1="<<cp1<<" cp2="<<cp2<<" cpR="<<cpR<<" cpRpr="<<cpRpr<<endl;
			pair<int,int> fpstab;  fpstab.first=-1;  fpstab.second=-1;
			for(int ip=0;ip<lfpcode.size();ip++){
			//	cout<<"lfpcode["<<ip<<"]="<<lfpcode[ip]<<" fbusy[ip]="<<fbusy[ip]<<endl;
			//	cout<<cp1-lfpcode[ip]<<" "<<cp2-lfpcode[ip]<<endl;
				if(cp1==lfpcode[ip] && fbusy[ip]==false)
					fpstab.first=ip;
				if(cp2==lfpcode[ip] && fbusy[ip]==false)
					fpstab.second=ip;
			//	cout<<"fpstab.first = "<<fpstab.first<<" fpstab.second = "<<fpstab.second<<endl;
			}
			if(fpstab.first>=0 && fpstab.second>=0){
				//cout<<"cp1="<<cp1<<" cp2="<<cp2<<" cpR="<<cpR<<endl;
				vector<int> pos(1);
				pos[0] = fpstab.first;
				fbusy[fpstab.first]=true;
				(*it).setListFP1(pos);
				posR.push_back(pos[0]);
				if(cp1==cp2){
					for(int ip=0;ip<lfpcode.size();ip++){
						if(cp2==lfpcode[ip] && fbusy[ip]==false && ip!=fpstab.first)
							fpstab.second=ip;
					}
				}
				pos[0] = fpstab.second;
				fbusy[fpstab.second]=true;
				(*it).setListFP2(pos);
				posR.push_back(pos[0]);
				//cout<<"pos1 = "<<pos.size()<<" pos2 = "<<pos.size()<<endl;
				vector<int> res1 = (*it).ListFP1();
				vector<int> res2 = (*it).ListFP2();
				//cout<<"res1 = "<<res1.size()<<" res2 = "<<res2.size()<<endl;
			}
			else{
				if(fpstab.first>=0 && cp2==cpRpr){
				//	cout<<"cp2="<<cp2<<" cp1="<<cp1<<" cpRpr="<<cpRpr<<endl;
					vector<int> pos(1);
					pos[0] = fpstab.first;
					fbusy[fpstab.first]=true;
					(*it).setListFP1(pos);
					(*it).setListFP2(posRpr);
				//	cout<<"pos1 = "<<pos.size()<<" pos2 = "<<posRpr.size()<<endl;

				}
				if(fpstab.second>=0 && cp1==cpRpr){
				//	cout<<"cp1="<<cp1<<" cp2="<<cp2<<" cpRpr="<<cpRpr<<endl;
					vector<int> pos(1);
					pos[0] = fpstab.second;
					fbusy[fpstab.second]=true;
					(*it).setListFP2(pos);
					(*it).setListFP1(posRpr);
					//cout<<"pos2 = "<<pos.size()<<" pos1 = "<<posRpr.size()<<endl;
				}
				//TODO: R->R1+R2, where R1 and R2 are resonances

			}
			posRpr = posR;
			cpRpr = cpR;
			nlstates.push_back((*it));
	}
	_ch.SetStates(nlstates);
}

//list of participate particle in this amplitude
vector<int> Amplitude::BusyParticle(){
	vector<int> busy;
//	cout<<"size of chain for ampl calc="<<_ch.GetStates().size()<<endl;
	list<TwoParticleState>  lstates = _ch.GetStates();
			for (list<TwoParticleState>::iterator it = lstates.begin(); it != lstates.end(); it++){
				vector<int> pos1 = (*it).ListFP1();
				vector<int> pos2 = (*it).ListFP2();
			//	cout<<"busy: pos1.size()="<<pos1.size()<<" pos2.size()="<<pos2.size()<<endl;
				for(int ip1=0;ip1<pos1.size();ip1++){
					int curp = pos1[ip1];
					busy.push_back(curp);
				//	cout<<curp<<", ";
				}
				for(int ip2=0;ip2<pos2.size();ip2++){
					int curp = pos2[ip2];
					busy.push_back(curp);
				//	cout<<curp<<", ";
				}
			}
			//cout<<""<<endl;
	return busy;
}

///Amplitude calculation lfp - list of 4-momentum f.p (final particles)
complex<double> Amplitude::Calc(list<TLorentzVector> llfp, int eJm, int elya){
	vector<TLorentzVector> lfp;
	for (list<TLorentzVector>::iterator it = llfp.begin(); it != llfp.end(); it++){
		//(*it).Print();
		lfp.push_back((*it));
	}
	complex<double> Ares(1,0);
	//cout<<"Ares = "<<Ares<<endl;
		list<TwoParticleState>  lstates = _ch.GetStates();
		for (list<TwoParticleState>::iterator it = lstates.begin(); it != lstates.end(); it++){
			TwoParticleState curst=(*it);
			vector<int> pos1 = curst.ListFP1();
			vector<int> pos2 = curst.ListFP2();
			TLorentzVector sumFP;
			TVector3 mom1;
			TVector3 mom2;
			//cout<<"particle#1 build from ";
			for(int ip1=0;ip1<pos1.size();ip1++){
				int curp = pos1[ip1];
			//cout<<curp<<", ";
			//lfp[curp].Vect().Print();
				sumFP +=lfp[curp];
				mom1 += lfp[curp].Vect();
			}
			//cout<<""<<endl;
			//cout<<"particle#2 build from ";
			for(int ip2=0;ip2<pos2.size();ip2++){
				int curp = pos2[ip2];
				//cout<<curp<<", ";
				//lfp[curp].Vect().Print();
				sumFP +=lfp[curp];
				mom2 += lfp[curp].Vect();
			}
		//	cout<<""<<endl;

			curst.ResRestFrame(mom1,mom2);//TODO: don't calculate boost for each projection separately
			vector<double> lya1 = curst.Lya1();
			vector<double> lya2 = curst.Lya2();
			vector<double> Jm = curst.Jm();
			double invmass = sumFP.M();
			//cout<<"current invar, mass = "<<invmass<<endl;
			complex<double> a_add;
			for(int im=0;im<Jm.size();im++){
				for(int ly1=0;ly1<lya1.size();ly1++){
					for(int ly2=0;ly2<lya2.size();ly2++){
						if(Jm[im]==eJm && (lya1[ly1]+lya2[ly2])==elya){
						a_add += jtos1s2(curst.Phi(),curst.Theta(),
								curst.J(), Jm[im],
								curst.S1(), lya1[ly1],
								curst.S2(),lya2[ly2],
								curst.AngMom(),(curst.S1()+curst.S2()),
								invmass, curst.mass1(), curst.mass2(),
								curst.mR(),curst.wR());
						//cout<<"a_add = "<<a_add<<" curst.AngMom().size() = "<<curst.AngMom().size()<<endl;
						}
					}
				}
			}
		//	cout<<"a_add_NEW = "<<a_add<<endl;
			Ares*=a_add;
		}
	//	cout<<"curr ampl = "<<Ares<<endl;
		return Ares;
}

///decay of particle j(m) with mass mass to s1(with helicity lya1) with mass1 and s2(with helicity lya2) with mass2
/// throught resonans with mass0 and width;
complex<double> Amplitude::jtos1s2(double phi, double theta,
				 double j, double m,
				 double s1,double lya1,
				 double s2,double lya2, vector<double> l, double s,
				 double mass,
				 double mass1, double mass2, double mass0, double width){
  double lya = lya1-lya2;
  complex<double> nfcoef;
  const double lsize = l.size();
 // cout<<"l.size() = "<<lsize<<endl;
  for(int il=0;il<lsize;il++){
    nfcoef += nfls(j, l[il], s, s1, lya1, s2, lya2, mass, mass1, mass2, mass0, width);
  }
  int sign = 1;
  if(m<0) sign =-1;
  if(m==0) sign = 0;
  complex<double> fphi(TMath::Cos(phi),sign*TMath::Sin(phi));
 // cout<<"fphi="<<fphi<<" "<<TMath::Cos(phi)<<" "<<sign*TMath::Sin(phi)<<endl;
  complex<double> amp = nfcoef*Wigner_d(j,m,lya, theta)*fphi;
 ///complex<double> amp = nfcoef*Wigner_d(j,m,lya, theta);
  //cout<<j<<", "<<m<<", "<<lya<<", "<<theta<<" "<<", "<<s<<", "<<lya1<<", "<<j<<", "<<lya2<<endl;
//  cout<<"nfcoef = "<<nfcoef<<" Wigner_d(j,m,lya, theta)="<<Wigner_d(j,m,lya, theta)<<" fphi="<<fphi<<endl;
  return amp;
}

/// one term of helicity amplitude for l,s
complex<double> Amplitude::nfls(double j, double l, double s, double s1, double lya1, double s2, double lya2,
			      double mass, double mass1, double mass2, double mass0, double width){
  double lya = lya1-lya2;
 // cout<<"Clebsch("<<l<<", 0, "<<s<<", "<<lya<<", "<<j<<", "<<lya<<") = "<<Clebsch(l,0,s,lya,j,lya)<<endl;
  //cout<<"Clebsch("<<s1<<", "<<lya1<<", "<<s2<<", "<<-lya2<<", "<<s<<", "<<lya<<") = "<<Clebsch(s1,lya1,s2,-lya2,s,lya)<<endl;
  complex<double> nfls_res = sqrt(2*l+1)*pwa(l, mass, mass1, mass2, mass0, width)*Clebsch(l,0,s,lya,j,lya)*Clebsch(s1,lya1,s2,-lya2,s,lya);
 // cout<<"pwa = "<<pwa(l, mass, mass1, mass2, mass0, width)<<" Clebsch(l,0,s,lya,j,lya) = "<<Clebsch(l,0,s,lya,j,lya)
//		  <<Clebsch(s1,lya1,s2,-lya2,s,lya)<<endl;
  return nfls_res;
}


/// Calculation of momentum for decay mass m doubleo particles with masses m1, m2
/// through resonance with mass m0
double Amplitude::momentum(double m, double m1, double m2, double m0){
  double mom = TMath::Sqrt(TMath::Power(m0,4)+TMath::Power(m1,4)+TMath::Power(m2,4)
			   -2*(TMath::Power(m0*m1,2)+TMath::Power(m1*m2,2)+TMath::Power(m2*m0,2)))/(2*m);

  return mom;
}

/// Calculation of barrier function
double Amplitude::Fl(const int n, const double p)
{
#define Pr 0.1973 // Gev/c corresponds to 1 fermi
  double ret;
  double z = (p/Pr) * (p/Pr);
  int m = n/2;
  switch (m) {
  case 0:
    ret = 1.0;
    break;
  case 1:
    ret = sqrt((2.0 * z)/(z + 1));
    break;
  case 2:
    ret = sqrt((13.0 * z * z)/(pow(z - 3.0,2.0) + 9.0 * z));
    break;
  case 3:
    ret = sqrt( (277.0 * pow(z,3.0))/(z * pow(z - 15.0,2.0) + 9.0 * pow(2.0 * z - 5.0,2.0)));
    break;
  case 4:
    ret = sqrt( (12746.0 * pow(z,4.0))/( pow(z * z - 45.0 * z + 105.0,2.0) + 25.0 * z * pow(2.0 * z - 21.0,2.0)));
    break;
  case 5:
    ret =  sqrt(z*z*z*z*z/(893025.0 +99225.0*z +6300.0*z*z +315.0*z*z*z +15.0*z*z*z*z +z*z*z*z*z));
    break;
  default:
    cerr << "Blatt-Weisskopf called for undefined L = " << n/2 << endl;
    ret = 1.0;
    break;
  }
  return ret;
}



///Partial wave amplitude a_ls, responceble for dynamics
/// 1st approximation is Breit-Wigner
complex<double> Amplitude::pwa(double l, double mass, double mass1, double mass2,  double mass0, double width){
   complex<double> im(0,1);
   double q = momentum(mass, mass1, mass2, mass0);
   double q0 = momentum(mass0, mass1, mass2, mass0);
   double fwidth =  width*mass0*q/(mass*q0);
   complex<double> a_ls = mass0*width*Fl(l,q)/(Fl(l,q0)*(mass0*mass0-mass*mass-im*mass0*fwidth));
   return a_ls;

}