#include "FullAmp.h"
#include "qft++.h"
#include <iostream>
#include "TMath.h"
//ClassImp(FullAmp);

//using namespace TMath;
using namespace std;

FullAmp::FullAmp(){
  // std::cout<<"new object FULL amplitude is created!"<<std::endl;
}

///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> FullAmp::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 = nfls(j, l, s, s1, lya1, s2, lya2, mass, mass1, mass2, mass0, width);
  complex<double> nfcoef;
  const double lsize = l.size();
  //  std::cout<<"lsize = "<<lsize<<std::endl;
  for(int il=0;il<lsize;il++){
    nfcoef += nfls(j, l[il], s, s1, lya1, s2, lya2, mass, mass1, mass2, mass0, width);
    // std::cout<<"j "<<j<<", l[il] "<<l[il]<<", s "<<s<<", s1 "<<s1<<", lya1 "<<lya1
    // 	     <<", s2 "<<s2<<", lya2 "<<lya2<<", mass "<<mass<<", mass1 "<<mass1<<", mass2 "
    // 	     <<mass2<<", mass0 "<<mass0<<", width "<<width<<std::endl;
  }
  //  complex<double> amp = nfcoef*Wigner_D(phi,theta,-phi,j,m,lya);
  int sign = 1;
  if(m<0) sign =-1;
  if(m==0) sign = 0;
  complex<double> fphi(TMath::Cos(phi),sign*TMath::Sin(phi));
  //  std::cout<<"phi="<<phi<<" TMath::Cos(phi) = "<<TMath::Cos(phi)<<" TMath::Sin(phi) = "<<TMath::Sin(phi)<<std::endl;
  //  std::cout<<"fphi="<<fphi<<std::endl;
  complex<double> amp = nfcoef*Wigner_d(j,m,lya, theta)*fphi;
  //complex<double> amp = nfcoef*Wigner_d(j,m,lya, theta); //no dependence on phi!

  // std::cout<<" nfcoef = "<< nfcoef<<std::endl;
  // std::cout<<"Wigner_d("<<j<<","<<m<<","<<lya<<","<<theta<<") = "<<Wigner_d(j,m,lya, theta)<<std::endl;
  // std::cout<<" amp = "<<amp<<std::endl;
  return amp;
}

///noresonance decay
complex<double> FullAmp::jtos1s2(double phi, double theta, 
				 double j, double m,
				 double s1,double lya1,
				 double s2,double lya2, vector<double> l, double s){
  double lya = lya1-lya2;
  complex<double> nfcoef;
  const double lsize = l.size();
  for(int il=0;il<lsize;il++)
    nfcoef += nfls(j, l[il], s, s1, lya1, s2, lya2);
  //  complex<double> amp = nfcoef*Wigner_D(phi,theta,-phi,j,m,lya);
  int sign = 1;
  if(m<0) sign =-1;
  if(m==0) sign = 0;
  complex<double> fphi(TMath::Cos(phi),sign*TMath::Sin(phi));
  complex<double> amp = nfcoef*Wigner_d(j,m,lya, theta)*fphi; // without dependence on phi!
  // std::cout<<"Wigner_d("<<j<<","<<m<<","<<lya<<","<<theta<<") = "<<Wigner_d(j,m,lya, theta)<<std::endl;
  return amp;
}

/// one term of helicity amplitude for l,s
complex<double> FullAmp::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;
  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);
  return nfls_res;
}

/// one term of helicity amplitude for l,s
/// for nonresonans case
complex<double> FullAmp::nfls(double j, double l, double s, double s1, double lya1, double s2, double lya2){
  double lya = lya1-lya2;
  complex<double> nfls_res = sqrt(2*l+1)*Clebsch(l,0,s,lya,j,lya)*Clebsch(s1,lya1,s2,-lya2,s,lya);
  return nfls_res;
}

// /// Sum for l,s !
// complex<double> FullAmp::NFj(double j, std::vector<double> l,std:: vector<double> s, double s1, double lya1, double s2, double lya2){
//   complex<double> NFjres =0;
//   for(int i=0; i<l.size();i++){
//   NFjres +=nfls(j, l[i], s[i], s1, lya1, s2, lya2);
//   }
//   // cout<<"NFjres = "<<NFjres<<endl;
//   return NFjres;
// }

/// Calculation of momentum for decay mass m doubleo particles with masses m1, m2
/// through resonance with mass m0
double FullAmp::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 FullAmp::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> FullAmp::pwa(double l, double mass, double mass1, double mass2,  double mass0, double width){

   complex<double> im(0,1);
   //   return im;
   double q = momentum(mass, mass1, mass2, mass0);

   double q0 = momentum(mass0, mass1, mass2, mass0);
   // std::cout<<"q="<<q<<" q0="<<q0<<std::endl;
   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));
   //  std::cout<<"a_ls = "<<a_ls<<std::endl;
   return a_ls;
  // // //  complex<double> restest(1,0);
  // // //  return  restest; 
  // complex<double> im(0,1);
  // complex<double> a_ls = width/(mass0*mass0-mass*mass-im*mass0*width);
  // //  complex<double> a_ls(1,0);
  //  //  std::cout<<"a_ls = "<<a_ls<<std::endl;
  //  return a_ls;
  // //  complex<double> restest(1,0);
  // //  return  restest; 

}