#include "Intensity.h"
#include "TMath.h"

#include <iostream>
#include <stdio.h>
#include <string>
#include <cstring>
#include <sstream>


//ClassImp(Intensity);

using namespace std;

Intensity::Intensity(){
  //  std::cout<<"new object Intensity is created!"<<std::endl;
  //  _a1 = new FullAmp();
}

Intensity::Intensity(vector<Resonance> Rf1, vector<Resonance> Rf2, vector<Resonance> Rf3,  Particle pi1,  Particle pi2,  Particle pi3){
  const int sizeRf1 = Rf1.size();
  double _maxJr1=0;
  double _maxJr2=0;
  double _maxJr3=0;
  _maxJr1n=0;
  _maxJr2n=0;
  _maxJr3n=0;
  for(int ir=0;ir<sizeRf1;ir++){
    _Rf1.push_back(Rf1[ir]);
    if(_maxJr1<_Rf1[ir].J()){ 
      _maxJr1=_Rf1[ir].J();
      _maxJr1n=ir;
    }
  }
  const int sizeRf2 = Rf2.size();
  for(int ir=0;ir<sizeRf2;ir++){
    _Rf2.push_back(Rf2[ir]);
    if(_maxJr2<_Rf2[ir].J()){
      _maxJr2=_Rf2[ir].J();
      _maxJr2n=ir;
    }
  }
  const int sizeRf3 = Rf3.size();
  for(int ir=0;ir<sizeRf3;ir++){
    _Rf3.push_back(Rf3[ir]);
    if(_maxJr3<_Rf3[ir].J()){
      _maxJr3=_Rf3[ir].J();
      _maxJr3n=ir;
    }
  }
  _pi1 = pi1;
  _pi2 = pi2;
  _pi3 = pi3;
  //_s = s;

  // //TODO: think how to properly use dencity matrix!
  // if(_Rf1.size()!=0){
  //   _J = _maxJr1;
  //   _Mj = _Rf1[_maxJr1n].Mj();
  // }
  // if(_Rf2.size()!=0){
  //   _J = _maxJr2;
  //   _Mj = _Rf2[_maxJr2n].Mj();
  // }
  // if(_Rf3.size()!=0){
  //   _J = _maxJr3;
  //   _Mj = _Rf3[_maxJr3n].Mj();
  // }


  // _J = -1;//TEST
  // _Mj = _Rf1[0].Mj();//TEST

  _s1 = _pi1.S();
 _lya1 = _pi1.lya();
  _s2 = _pi2.S();
  _lya2 = _pi2.lya();
  _s3 = _pi3.S();
  _lya3 = _pi3.lya();

  _mass1 = _pi1.M();
  _mass2 = _pi2.M();
  _mass3 = _pi3.M();

  // _massR = _Rf.M();
  // _widthR = _Rf.W();
 
  int sizelya1 = _lya1.size();
  int sizelya2 = _lya2.size();
  int sizelya3 = _lya3.size();

  for(int ilya1=0;ilya1<sizelya1;ilya1++){
    for(int ilya2=0;ilya2<sizelya2;ilya2++){
      _lya.push_back(_lya1[ilya1]-_lya2[ilya2]);
    }
  }

  for(int ilya2=0;ilya2<sizelya2;ilya2++){
    for(int ilya3=0;ilya3<sizelya3;ilya3++){
      _lya.push_back(_lya2[ilya2]-_lya3[ilya3]);
    }
  }

  for(int ilya1=0;ilya1<sizelya1;ilya1++){
    for(int ilya3=0;ilya3<sizelya3;ilya3++){
      _lya.push_back(_lya1[ilya1]-_lya3[ilya3]);
    }
  }
  //  _a1 = new FullAmp();
}

void Intensity::PrintPar(){
  std::cout<<"Parameters for intensity calculation: "<<std::endl;
  //  std::cout<<" mass = "<<_massR<<", width = "<<_widthR<<std::endl;
  std::cout<<" mass1 = "<<_mass1<<" mass2 = "<<_mass2<<" mass3 = "<<_mass3<<std::endl;
 std::cout<<" sizelya = "<<_lya.size()<<std::endl;
  // std::cout<<" max Spin = "<<_J<<std::endl;
  //  std::cout<<" sizemj ="<<_Mj.size()<<" sizelya = "<<_lya.size()<<std::endl;
}



///calculation of coherent sum for 12, 23 and 13 decay
std::complex<double> Intensity::Calc(double phi13, double theta13, double invmass12, 
				     double phi21, double theta21, double invmass23,
				     double phi12, double theta12, double invmass13){
  int sizemjc;
  if(_Rf1.size()!=0) sizemjc = _Rf1[_maxJr1n].Mj().size();
  if(_Rf2.size()!=0)
    if(sizemjc<_Rf2[_maxJr2n].Mj().size())
      sizemjc = _Rf2[_maxJr2n].Mj().size();
  if(_Rf3.size()!=0)
    if(sizemjc<_Rf3[_maxJr3n].Mj().size())
      sizemjc = _Rf3[_maxJr3n].Mj().size();
  const int sizemj = sizemjc;
  const int sizelya = _lya.size();
  std::complex<double> a2calc[sizemj][sizelya];
  std::complex<double> a2calc_conj[sizelya][sizemj];

  const int sizerecvec1 =  _Rf1.size();
  const int sizerecvec2 =  _Rf2.size();
  const int sizerecvec3 =  _Rf3.size();
  // std::cout<<"sizevec 1,2,3: "<<sizerecvec1<<", "<<sizerecvec2<<", "<<sizerecvec3
  // 	   <<"  sizemj="<<sizemj<<" sizelya="<<sizelya<<std::endl;

  ///index schift for correct sum around Jm
  std::vector<int> addind1;
  int maxJ1 = _Rf1[_maxJr1n].J();
  for(int ir=0;ir<sizerecvec1;ir++){
    if(_Rf1[ir].Mj().size()<sizerecvec1) addind1.push_back(maxJ1);
    else  addind1.push_back(0);
  }
  std::vector<int> addind2;
  int maxJ2 = _Rf2[_maxJr2n].J();
  for(int ir=0;ir<sizerecvec2;ir++){
    if(_Rf2[ir].Mj().size()<sizerecvec2) addind2.push_back(maxJ2);
    else  addind2.push_back(0);
  }
  std::vector<int> addind3;
  int maxJ3 = _Rf3[_maxJr3n].J();
  for(int ir=0;ir<sizerecvec3;ir++){
    if(_Rf3[ir].Mj().size()<sizerecvec3) addind3.push_back(maxJ3);
    else  addind3.push_back(0);
  }
  int sizelya12 = _lya1.size()*_lya2.size();
  int sizelya23 = _lya2.size()*_lya3.size();
  int sizelya13 = _lya3.size()*_lya1.size();
  //  std::cout<<" sizelya12="<<sizelya12<<" sizelya23="<<sizelya23<<" sizelya13="<<sizelya13<<std::endl;
  for(int ij=0;ij<sizemj;ij++){
    for(int mk=0;mk<sizelya;mk++){
      a2calc[ij][mk] =0;
      for(int ir=0;ir<sizerecvec1;ir++){
	if(ij<_Rf1[ir].Mj().size() && mk<sizelya12){ 
	  int addind = addind1[ir];
	  //  cout<<"_Rf1[ir].J() = "<<_Rf1[ir].J()<<" M("<<(ij+addind)<<")="<<_Rf1[ir].Mj(ij+addind)<<endl;
	  // std::cout<<"_lya1["<<mk<<"]="<<_lya1[mk]<<std::endl;
	  // std::cout<<"_lya2["<<mk<<"]="<<_lya2[mk]<<std::endl;
	  std::complex<double> a_add = _a1.jtos1s2(phi13,theta13, _Rf1[ir].J(), _Rf1[ir].Mj(ij), 
						   _s1,_lya1[mk], 
						   _s2,_lya2[mk], 
						   _Rf1[ir].l(), (_s1+_s2),					  
						   invmass12, _mass1, _mass2, _Rf1[ir].M(), _Rf1[ir].W());
	  //	      std::cout<<"add="<<a_add<<std::endl;
	  a2calc[ij+addind][mk] += a_add;
	  //     std::cout<<"mk="<<mk<<std::endl;
	}
      }
      
      for(int ir=0;ir<sizerecvec2;ir++){
	if(ij<_Rf2[ir].Mj().size() && mk<(sizelya12+sizelya23) && mk>=sizelya12){ 
	    int addind = addind2[ir];
	    // std::cout<<"_lya2["<<mk-sizelya12<<"]="<<_lya2[mk-sizelya12]<<std::endl;
	    // std::cout<<"_lya3["<<mk-sizelya12<<"]="<<_lya3[mk-sizelya12]<<std::endl;
	    std::complex<double> a_add = _a1.jtos1s2(phi21,theta21, _Rf2[ir].J(), _Rf2[ir].Mj(ij), 
						     _s2,_lya2[mk-sizelya12], 
						     _s3,_lya3[mk-sizelya12], 
						     _Rf2[ir].l(), (_s2+_s3),
						     invmass23, _mass2, _mass3, _Rf2[ir].M(), _Rf2[ir].W());
	    //  std::cout<<"add="<<a_add<<std::endl;
	    a2calc[ij+addind][mk] += a_add;
	    // std::cout<<"mk="<<mk<<std::endl;
	  }
	  }
	
	
	
	for(int ir=0;ir<sizerecvec3;ir++){
	  if(ij<_Rf3[ir].Mj().size() && mk>=(sizelya23+sizelya12)){ 
	    int addind = addind3[ir];
	    // std::cout<<"_lya1["<<mk-sizelya12-sizelya23<<"]="<<_lya1[mk-sizelya12-sizelya23]<<std::endl;
	    // std::cout<<"_lya3["<<mk-sizelya12-sizelya23<<"]="<<_lya3[mk-sizelya12-sizelya23]<<std::endl;

	    std::complex<double> a_add = _a1.jtos1s2(phi12,theta12, _Rf3[ir].J(), _Rf3[ir].Mj(ij), 
						     _s1,_lya1[mk-sizelya12-sizelya23], 
						     _s3,_lya3[mk-sizelya12-sizelya23], 
						     _Rf3[ir].l(), (_s1+_s3), 
						     invmass13, _mass1, _mass3, _Rf3[ir].M(), _Rf3[ir].W());
	    //   std::cout<<"add="<<a_add<<std::endl;
	    a2calc[ij+addind][mk] += a_add;
	    //std::cout<<"mk="<<mk<<std::endl;
	  }
	}
      	
	
      }   
    }
    
    for(int ij=0;ij<sizemj;ij++){
      for(int mk=0;mk<sizelya;mk++){ 
	//	std::cout<<a2calc[ij][mk]<<", ";
      real(a2calc_conj[mk][ij]) = real(a2calc[ij][mk]);
      imag(a2calc_conj[mk][ij]) = -imag(a2calc[ij][mk]); 
    }
      //  std::cout<<""<<std::endl;
  }
  
  std::complex<double> a2rhoa2c[sizemj][sizelya];
  for(int ij=0;ij<sizemj;ij++){
    for(int mk=0;mk<sizelya;mk++){
      //  a2rhoa2c[ij][mk] += a2calc[ij][mk]*a2calc_conj[mk][ij]/(2*_J+1.);
      a2rhoa2c[ij][mk] += a2calc[ij][mk]*a2calc_conj[mk][ij];//TEST
    }
  }

  // for(int ij=0;ij<sizemj;ij++){
  //   for(int mk=0;mk<sizelya;mk++){
  //     std::cout<<a2rhoa2c[ij][mk]<<", ";
  //   }
  //   std::cout<<""<<std::endl;
  // }

  std::complex<double> a2rhoa2csumlya[sizemj];
  for(int ij=0;ij<sizemj;ij++){
    for(int mk=0;mk<sizelya;mk++){
      a2rhoa2csumlya[ij] += a2rhoa2c[ij][mk];
    }
  }
  
  std::complex<double> intensity2;
  for(int ij=0;ij<sizemj;ij++){
    intensity2 += a2rhoa2csumlya[ij];
  }
 
  return intensity2;
  //  return  abs(a2calc[0][0]);
}