/*
 * Intensity.cpp
 *
 *  Created on: Jan 24, 2012
 *      Author: Anastasia Karavdina
 */

#include "Intensity.h"

Intensity::Intensity(list<Amplitude> lAmpl) {
	cout<<"size of Ampl fo Intencity calc: "<<lAmpl.size()<<endl;
	_lAmpl = lAmpl;
	double maxJ=0;
	int maxLya=0;
	for (list<Amplitude>::iterator it = _lAmpl.begin(); it != _lAmpl.end(); it++){
		Amplitude curA =(*it);
		Chain ch = curA.GetChain();
		double fmaxJ = ch.Jc();
		if(maxJ<fmaxJ) maxJ=fmaxJ;

		vector<double> curLya = ch.LyaFP();
		for(int ci=0;ci<curLya.size();ci++)
			_lya.push_back(curLya[ci]);
		for(int ji=0;ji<_lya.size();ji++){
			for(int jk=0;jk<_lya.size();jk++){
				if(ji==jk) continue;
				if(_lya[ji]==_lya[jk]) _lya.erase(_lya.begin()+ji,_lya.begin()+ji+1);
			}
		}
		int fMaxlya = _lya.size();
		if(maxLya<fMaxlya) maxLya=fMaxlya;
	//	cout<<" fmaxJ = "<<fmaxJ<<" fMaxlya = "<<fMaxlya<<endl;
	}
	_sizemj=2*maxJ+1;
	_sizelya =maxLya;
	cout<<"_sizemj = "<<_sizemj<<" _sizelya = "<<_sizelya<<endl;
}

Intensity::~Intensity() {
}

//check all combination before calculation in the loop!
vector<int> Intensity::Combinations(vector<int> lfpcode){
	vector <int> allComb;// vector for checking used combination
	for(int i=0;i<lfpcode.size();i++)
		allComb.push_back(i);
	//int countcomb=0;
	for (list<Amplitude>::iterator it = _lAmpl.begin(); it != _lAmpl.end(); it++){
		Amplitude *curA =&(*it);
		vector<int> busy = curA->BusyParticle();
//		cout<<"size of busy-vector is "<<busy.size()<<endl;
//		for(int ib=0;ib<busy.size();ib++)
//			cout<<busy[ib]<<", ";
//		cout<<""<<endl;
//		cout<<"combinatorics!"<<endl;
		//of final state contain particles of one species, z.B. 1=2
		//=> coherent sum of amplitude 2(1,3) and 1(2,3)
		int tmpcode1,tmpcode2;
		bool calcAmpl=true;
		pair<int,int> identpart(-1,-1);// identical particles in decay
		if(busy.size()==2){
			if(lfpcode[busy[0]]==lfpcode[busy[1]]){//identical particles in decay
				identpart.first = busy[0];
				identpart.second = busy[1];
			}
		//	cout<<"identpart.first="<<identpart.first<<" identpart.second="<<identpart.second<<endl;
		}
		for(int ipc1=0;ipc1<lfpcode.size();ipc1++){
			for(int ipc2=0;ipc2<lfpcode.size();ipc2++){
				bool participate=false;//particle participate in decay

				if(busy.size()==2){
					if((ipc2==busy[0] && ipc1==busy[1])
							||(ipc2==busy[1] && ipc1==busy[0]) ) participate=true;
				}

				if(ipc1==ipc2 || participate) continue;

				if(lfpcode[ipc1]==lfpcode[ipc2]){
					int minv=-lfpcode.size()+ipc1;
					int val1 = allComb[(allComb.size()+minv)];
					minv=-lfpcode.size()+ipc2;
					int val2 = allComb[(allComb.size()+minv)];
					int maxv=-lfpcode.size()+ipc1+1;
					minv=-lfpcode.size()+ipc1;
					int pos1ch = lfpcode.size()+minv;//position#1
					replace((allComb.end()+minv),(allComb.end()+maxv),val1, val2);
					minv=-lfpcode.size()+ipc2;
					maxv=-lfpcode.size()+ipc2+1;
					int pos2ch = lfpcode.size()+minv;//position#2
					replace((allComb.end()+minv),(allComb.end()+maxv),val2, val1);
//cout<<"pos1ch = "<<pos1ch<<" pos2ch = "<<pos2ch<<endl;
					vector<int> curComb;
					curComb.insert(curComb.begin(),(allComb.end()-lfpcode.size()),allComb.end());
					if(allComb.size()>lfpcode.size()){
						for(int i=0;i<((allComb.size()/lfpcode.size())-1);i++){
							vector<int> tmpComb;
							tmpComb.insert(tmpComb.begin(),(allComb.begin()+i*lfpcode.size()),(allComb.begin()+(i+1)*lfpcode.size()));
							if(curComb==tmpComb) break;
							//else countcomb++;
						}
					}
//					for(int k=0;k<allComb.size();k++)
//						cout<<allComb[k]<<", ";
//						cout<<""<<endl;
					allComb.insert (allComb.end(),(allComb.end()-lfpcode.size()),allComb.end());
				}
			}
		}
	}
	//cout<<"# of combination:"<<countcomb<<endl;
//	cout<<"lfpcode.size() = "<<lfpcode.size()<<endl;
	//allComb.erase((allComb.end()-lfpcode.size()),allComb.end());
	//allComb.erase((allComb.end()-lfpcode.size()),allComb.end());//TODO: it's just test!
	//allComb.erase((allComb.end()-18*lfpcode.size()),(allComb.end()-12*lfpcode.size()));//delete equal combinations
	//cout<<"This combinations may participate in calculation:"<<endl;
	//for(int k=0;k<allComb.size();k++)
	//		cout<<allComb[k]<<", ";
	//	cout<<""<<endl;

		//TODO: delete equal combinations automaticly! (important for identical particles in final state!!)
	/*//delete equal combinations
	for(int i=0;i<((allComb.size()/lfpcode.size())-1);i++){

	}*/
	vector<int> nodouble (18,0);//TODO: this is only for 3 pi0 as f.s
	//cout<<"This combinations will participate in calculation:"<<endl;
	for(int k=0;k<nodouble.size();k++){
		nodouble[k]=allComb[k];
	//			cout<<nodouble[k]<<", ";
	}
	//		cout<<""<<endl;
	//return allComb;
	return nodouble;//TODO: this is only for 3 pi0 as f.s
}

double Intensity::Calc(list<TLorentzVector> llfp,vector<int> lfpcode){
	vector<TLorentzVector> lfp;
	for (list<TLorentzVector>::iterator it = llfp.begin(); it != llfp.end(); it++)
		lfp.push_back((*it));

	const int sizemj=_sizemj;
	const int sizelya=_sizelya;
	complex<double> aCalc[sizemj][sizelya];
	complex<double> aCalc_conj[sizelya][sizemj];


	vector<int> allComb=Combinations(lfpcode);
	for (list<Amplitude>::iterator it = _lAmpl.begin(); it != _lAmpl.end(); it++){
		Amplitude *curA =&(*it);
		for(int i=0;i<=((allComb.size()/lfpcode.size())-1);i++){
			vector<int> tmpComb;
			tmpComb.insert(tmpComb.begin(),(allComb.begin()+i*lfpcode.size()),(allComb.begin()+(i+1)*lfpcode.size()));
			vector<TLorentzVector> curlfp(lfp.size());
			for(int ilv=0;ilv<lfp.size();ilv++){
				int newpos = tmpComb[ilv];
			//	cout<<"old pos="<<ilv<<" newpos = "<<newpos<<" lfp.size()="<<lfp.size()<<endl;
				curlfp[ilv]=lfp[newpos];
			}
			for(int Jm=0;Jm<=sizemj;Jm++){
				for(int lya=0;lya<sizelya;lya++){
					int ejm = -(sizemj-1)*0.5+Jm;
				//	int elya = -(sizelya-1)*0.5+lya;
					int elya = _lya[lya];
					list<TLorentzVector> llfpcur;
					for(int kl=0;kl<lfp.size();kl++)
						llfpcur.push_back(curlfp[kl]);
					complex<double> anowis = curA->Calc(llfpcur, ejm, elya);
				//	cout<<"anowis = "<<anowis<<endl;
					aCalc[Jm][lya] += anowis;
				}
			}
		}
	}

	for(int Jm=0;Jm<sizemj;Jm++){
		for(int lya=0;lya<sizelya;lya++){
			//cout<<"aCalc[Jm][lya] = "<<aCalc[Jm][lya]<<endl;
			real(aCalc_conj[lya][Jm]) = real(aCalc[Jm][lya]);
			imag(aCalc_conj[lya][Jm]) = -imag(aCalc[Jm][lya]);
		}
	}

	//complex<double> aaCalc(real(aCalc),-imag(aCalc));
	//complex<double> itmp =  aCalc*aaCalc;
	complex<double> a2rhoa2c[sizemj][sizelya];
	for(int ij=0;ij<sizemj;ij++)
		for(int mk=0;mk<sizelya;mk++)
			a2rhoa2c[ij][mk] += aCalc[ij][mk]*aCalc_conj[mk][ij];

	complex<double> a2rhoa2csumlya[sizemj];
	for(int ij=0;ij<sizemj;ij++)
		for(int mk=0;mk<sizelya;mk++)
			a2rhoa2csumlya[ij] += a2rhoa2c[ij][mk];

	complex<double> intensity;
	for(int ij=0;ij<sizemj;ij++)
		intensity += a2rhoa2csumlya[ij];

	if(fabs(imag(intensity))>1e-6)
		cout<<"!Attention! imaginary part of intensity="<<imag(intensity)<<endl;
	/*cout<<"intensity="<<real(intensity)<<endl;
	cout<<"=============================="<<endl;*/
	return real(intensity);
}