#include <cassert>
#include <iostream>
#include <string>
#include <sstream>
#include <cstdlib>
#include <cmath>

#include "TLorentzVector.h"
#include "TLorentzRotation.h"

#include "IUAmpTools/Kinematics.h"
#include "AMPTOOLS_AMPS/Amp_R.h"
#include "AMPTOOLS_AMPS/clebschGordan.h"
#include "AMPTOOLS_AMPS/wignerD.h"

#include "UTILITIES/BeamProperties.h"

struct anglesGJHzlm_t {
  double cosThetaGJ;
  double cosTheta_H;
  double phiGJ;
  double phi_H;
  double Phi;
};

anglesGJHzlm_t getAnglesGJHzlm_R(TLorentzVector resonance, TLorentzVector daughter1,
				 TLorentzVector daughter2,TLorentzVector beam, 
				 double polAngle, double polFraction){

  TLorentzVector target(0,0,0,0.938);
  TLorentzVector cm = beam+target;

  TVector3 cmRestBoost = cm.BoostVector();
  TLorentzVector beam_cm = beam;
  TLorentzVector resonance_cm = resonance;
  TLorentzVector daughter1_cm = daughter1;
  TLorentzVector daughter2_cm = daughter2;
  beam_cm.Boost(-1.0*cmRestBoost);
  resonance_cm.Boost(-1.0*cmRestBoost);
  daughter1_cm.Boost(-1.0*cmRestBoost);
  daughter2_cm.Boost(-1.0*cmRestBoost);

  TVector3 resRestBoost = resonance_cm.BoostVector();
  TLorentzVector beam_res = beam_cm;
  TLorentzVector resonance_res = resonance_cm;
  TLorentzVector daughter1_res = daughter1_cm;
  TLorentzVector daughter2_res = daughter2_cm;
  beam_res.Boost(-1.0*resRestBoost);
  resonance_res.Boost(-1.0*resRestBoost);
  daughter1_res.Boost(-1.0*resRestBoost);
  daughter2_res.Boost(-1.0*resRestBoost);

  //endp1                                     
  TVector3 resRestBoost2 = resonance.BoostVector();
  TLorentzVector resonance_res2 = resonance;
  resonance_res2.Boost(-1.0*resRestBoost2);

  TVector3 vRestBoost = daughter1_res.BoostVector();
  TLorentzVector daughter2_v = daughter2_res;
  daughter2_v.Boost(-1.0*vRestBoost);

  TVector3 beam_cm_vec = beam_cm.Vect().Unit();
  TVector3 daughter1_res_vec = daughter1_res.Vect().Unit();

  TVector3 daughter2_vec = daughter2_v.Vect().Unit();

  GDouble Pgamma;

  TVector3 z = resonance_cm.Vect().Unit();
  TVector3 y = (beam_cm_vec.Cross(z)).Unit();
  TVector3 x = (y.Cross(z)).Unit();
  //endp2                        
  TVector3 Angles(daughter1_res_vec.Dot(x),
                  daughter1_res_vec.Dot(y),
                  daughter1_res_vec.Dot(z));

  
  //double Phi = (resonance_cm.Phi() - eps.Phi());                                                                            
  double cosTheta = Angles.CosTheta();
  double theta = Angles.Theta();
  //double phi = Angles.Phi();                                                                                                
  //double psi = daughter1_res_vec.Phi();                                                                                     

  TVector3 z_H = daughter1_res.Vect().Unit();
  TVector3 y_H = (z.Cross(z_H)).Unit();
  TVector3 x_H = (y_H.Cross(z_H)).Unit();

  TVector3 Angles_H(daughter2_vec.Dot(x_H),
                    daughter2_vec.Dot(y_H),
                    daughter2_vec.Dot(z_H));

  double cosTheta_H = Angles_H.CosTheta();
  double theta_H = Angles_H.Theta();
  double phi_H = Angles_H.Phi();

  //////////////////////////////////////////////////////////////                                                              
  //Rest-frame of the meson-resonance: Gottfried-Jackson frame//           
  //////////////////////////////////////////////////////////////                                                              
  TVector3 zGJ = beam_res.Vect().Unit();
  TVector3 yGJ = y;
  TVector3 xGJ = (yGJ.Cross(zGJ)).Unit();

  TVector3 AnglesGJ(daughter1_res_vec.Dot(xGJ),
                    daughter1_res_vec.Dot(yGJ),
                    daughter1_res_vec.Dot(zGJ));

  double cosThetaGJ = AnglesGJ.CosTheta();
  double thetaGJ = AnglesGJ.Theta();
  double phiGJ = AnglesGJ.Phi();


  //Phi
  TVector3 eps(cos(polAngle*TMath::DegToRad()), sin(polAngle*TMath::DegToRad()), 0.0); // beam polarization vector
  double Phi = atan2(yGJ.Dot(eps), beam.Vect().Unit().Dot(eps.Cross(yGJ)));

  anglesGJHzlm_t anglesGJH;
  anglesGJH.cosThetaGJ = cosThetaGJ;
  anglesGJH.cosTheta_H = cosTheta_H;
  anglesGJH.phiGJ = phiGJ;
  anglesGJH.phi_H = phi_H;
  anglesGJH.Phi = Phi;
  //endp3                          
  return anglesGJH;
};




double moZero_R(double M,double m1,double m2){
  double moZeroVal = -1;
  double alpha = M*M;
  double beta = m1*m1;
  double gamma = m2*m2;
  double lambda = pow(alpha,2) + pow(beta,2) + pow(gamma,2) -
    2*alpha*beta - 2*alpha*gamma - 2*beta*gamma;
  if (lambda < 0)lambda = 0;
  moZeroVal = sqrt(lambda)/(2*M);
  return moZeroVal;
}
double F_R(double q,int l){
  double z = pow(q/0.1973,2);
  double FVal = 1;
  if (l == 1) {
    FVal = sqrt(2*z/(z+1));
  }
  if (l == 2) {
    FVal = sqrt(13*pow(z,2)/(pow(z-3,2) + 9));
  }

  FVal = 1;
  return FVal;
}

complex <GDouble> Gamma_R(double m, double q, double m0, double q0, double Gamma0, int l){
  complex <GDouble> i(0,1);  
  complex <GDouble> GammaVal(0,0);  

  //double q = 0; //breakup momentum in restframe of decaying particle
  double z = pow(q/0.1973,2);

  GammaVal = Gamma0*m0*q*pow(F_R(q,l),2)/(m*q0*pow(F_R(q0,l),2));
  
  GammaVal = Gamma0;
  return GammaVal;

}

complex <GDouble> BW_R(double m, double q, double m0, double q0, double Gamma0, int l){
  complex <GDouble> i(0,1);  
  complex <GDouble> BWVal(0,0);  
  //BWVal = m0*Gamma0/(pow(m0,2) - pow(m,2) - i*m0*Gamma(m,q,m0,q0,Gamma0,l));
  BWVal = m0*Gamma0/(m0*m0 - m*m - i*m0*Gamma0);
  
  return BWVal;
}



Amp_R::Amp_R( const vector< string >& args ) :
UserAmplitude< Amp_R >( args )
{
  nArgs = args.size();

  cout<<"nArgs = "<<nArgs<<endl;

  assert( nArgs >= 9 );
  m_j = atoi( args[0].c_str() );//total J
  m_m = atoi( args[1].c_str() );//total Jz
  m_l = atoi( args[2].c_str() );//total L
  m_s = atoi( args[3].c_str() );//S
  m_r = atoi( args[4].c_str() );//Rxn type: 0 -> nothing special, 1-> a0(980), 2-> kpStar, 3-> kmStar
  m_e = atoi( args[5].c_str() );//Reflectivity (epsilon)

  cout<<"m_r, abs(m_r) = "<<m_r<<" , "<<abs(m_r)<<endl;
  
  polAngle = AmpParameter( args[6] );
  registerParameter( polAngle );


  
  polFraction = atof(args[7].c_str());
  
  // BeamProperties configuration file
  if (polFraction == 0){
    TString beamConfigFile = args[7].c_str();
    BeamProperties beamProp(beamConfigFile);
    polFrac_vs_E = (TH1D*)beamProp.GetPolFrac();
  }


  
  
  particle_decay = args[8];

  if ( nArgs == 11 ) {
    bwCenter = AmpParameter( args[9] );
    bwWidth = AmpParameter( args[10] );

    // need to register any free parameters so the framework knows about them                                                                      
    registerParameter( bwCenter );
    registerParameter( bwWidth );

    cout<<"bwCenter = "<<bwCenter<<endl;
    cout<<"bwWidth = "<<bwWidth<<endl;

  }




}


complex< GDouble >
Amp_R::calcAmplitude( GDouble** pKin ) const {
  GDouble Pgamma;  
  if(polFraction > 0.) { // for fitting with constant polarization 
    Pgamma = polFraction;
  }
  else{
    int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]);
    if (bin == 0 || bin > polFrac_vs_E->GetXaxis()->GetNbins()){
      Pgamma = 0.;
    }
    else Pgamma = polFrac_vs_E->GetBinContent(bin);
    //cout<<pKin[0][0]<<" "<<bin<<" "<<Pgamma<<endl;
  }
  
  TLorentzVector   beam(pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0]);
  TLorentzVector recoil(pKin[1][1], pKin[1][2], pKin[1][3], pKin[1][0]);
  TLorentzVector     p1(pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0]); 
  TLorentzVector     p2(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]);
  TLorentzVector     p3(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]);

  TLorentzVector resonance = p1+p2+p3;
  double mResonance = resonance.Mag();
  
  //Isobar KmStar                                                                                                             
  TLorentzVector daughter1 = p2+p3;
  TLorentzVector daughter2 = p2;
  double mKmPi0 = daughter1.Mag();

  anglesGJHzlm_t anglesGJH_IsoKmStar = getAnglesGJHzlm_R(resonance,daughter1,daughter2,beam,polAngle,polFraction);
  
  double cosThetaGJ_IsoKmStar = anglesGJH_IsoKmStar.cosThetaGJ;
  double cosTheta_H_IsoKmStar = anglesGJH_IsoKmStar.cosTheta_H;
  double phi_H_IsoKmStar = anglesGJH_IsoKmStar.phi_H;
  double phiGJ_IsoKmStar = anglesGJH_IsoKmStar.phiGJ;
  
  double Phi = anglesGJH_IsoKmStar.Phi;

  //Isobar KpStar                                                                                                             
  daughter1 = p1+p3;
  daughter2 = p1;
  double mKpPi0 = daughter1.Mag();

  anglesGJHzlm_t anglesGJH_IsoKpStar = getAnglesGJHzlm_R(resonance,daughter1,daughter2,beam,polAngle,polFraction);
  
  double cosThetaGJ_IsoKpStar = anglesGJH_IsoKpStar.cosThetaGJ;
  double cosTheta_H_IsoKpStar = anglesGJH_IsoKpStar.cosTheta_H;
  double phi_H_IsoKpStar = anglesGJH_IsoKpStar.phi_H;
  double phiGJ_IsoKpStar = anglesGJH_IsoKpStar.phiGJ;
  
  //Isobar KK                                                                                                                 
  daughter1 = p1+p2;
  daughter2 = p1;
  double mKK = daughter1.Mag();

  anglesGJHzlm_t anglesGJH_IsoKK = getAnglesGJHzlm_R(resonance,daughter1,daughter2,beam,polAngle,polFraction);
  
  double cosThetaGJ = anglesGJH_IsoKK.cosThetaGJ;
  double cosTheta_H = anglesGJH_IsoKK.cosTheta_H;
  double phi_H = anglesGJH_IsoKK.phi_H;
  double phiGJ = anglesGJH_IsoKK.phiGJ;


  if (abs(m_r) == 2 || abs(m_r) == 20 || abs(m_r) == 21){
    cosThetaGJ = cosThetaGJ_IsoKpStar; 
    cosTheta_H = cosTheta_H_IsoKpStar; 
    phiGJ      = phiGJ_IsoKpStar; 
    phi_H      = phi_H_IsoKpStar; 
  }

  if (abs(m_r) == 3 || abs(m_r) == 30 || abs(m_r) == 31){
    cosThetaGJ = cosThetaGJ_IsoKmStar; 
    cosTheta_H = cosTheta_H_IsoKmStar; 
    phiGJ      = phiGJ_IsoKmStar; 
    phi_H      = phi_H_IsoKmStar; 
  }

  //////////////////////////
  complex <GDouble> i(0,1);  
  complex <GDouble> amplitude(1,0);
  complex <GDouble> amplitudeTest(0,0);

  GDouble Factor = sqrt(1 + m_s * Pgamma);
  complex <GDouble> zjm = 0;
  //complex <GDouble> zjmTest = 0;

  //complex< GDouble > rotateY = polar(1., -1.*Phi);

  int j_0 = m_j;
  int j_1 = m_l;
  int j_2 = m_s;
  int m_0 = m_m;

  complex <GDouble> preFactor(1,0);
  //double pVal = d1moResFrame;
  //double qVal = d2moIsoFrame;
  //double m1Val = daughter1.Mag();//Mass kpkm
  //double m2Val = p3.Mag();//Mass pi0
  
  //cout<<"m_r = "<<m_r<<endl;
  //if (m_r == 80){//phi
  //double Gamma0 = 0.008;
  //double m0Val = 1.019461;
  //double mVal = mKK;
  //complex <GDouble> bwVal(1,0);
  //bwVal = m0Val*Gamma0/(m0Val*m0Val - mVal*mVal - i*m0Val*Gamma0);
  //preFactor *=  bwVal;
  //}
  
  if (abs(m_r) == 1){//a0(980)
    double Gamma0 = 0.05;
    double m0Val = 0.980;
    double mVal = mKK;
    complex <GDouble> bwVal(1,0);
    bwVal = m0Val*Gamma0/(m0Val*m0Val - mVal*mVal - i*m0Val*Gamma0);
    preFactor *=  bwVal;
  }

  if (abs(m_r) == 2 || abs(m_r) == 3 || abs(m_r) == 20 || abs(m_r) == 30){//Kstar
    double Gamma0 = 0.052;
    double m0Val = 0.89167;
    double mVal = mKpPi0;

    ///  WAS : if (abs(m_r) == 3) mVal = mKmPi0;
    ///  NOW : a below
    if (abs(m_r) == 3 || abs(m_r) == 30) mVal = mKmPi0;
    complex <GDouble> bwVal(0,0);
    bwVal = m0Val*Gamma0/(m0Val*m0Val - mVal*mVal - i*m0Val*Gamma0);
    
    if (abs(m_r) == 2 || abs(m_r) == 3) {
      if (mVal < 0.83) bwVal = 0; 
      if (mVal > 0.95) bwVal = 0; 
    }
    preFactor *=  bwVal;
  }
  
  if (m_r < 0){
    double Gamma0 = bwWidth;
    double m0Val = bwCenter;
    double mVal = mResonance;
    complex <GDouble> bwVal(1,1);
    bwVal = m0Val*Gamma0/(m0Val*m0Val - mVal*mVal - i*m0Val*Gamma0);
    preFactor *=  bwVal;
  }
  //if (m_r == 1 || m_r == 11){//eta(1295)
  //double Gamma0 = 0.05425;
  //double m0Val = 1.294;
  //double mVal = mResonance;
  //double q0Val = moZero(m0Val,m1Val,m2Val);
  //complex <GDouble> fpl = F(pVal,m_l);
  //complex <GDouble> fqs = F(qVal,m_s);
  //complex <GDouble> bwVal(1,1);
  //bwVal = BW(mVal,qVal,m0Val,q0Val,Gamma0,m_l);
  //preFactor *=  sqrt(2*m_l+1)*sqrt(2*m_s+1)*fpl*fqs*bwVal;
  //}
  //if (m_r == 2 || m_r == 12){//f_1(1285)
  //double Gamma0 = 0.02567;
  //double m0Val = 1.2819;
  //double mVal = mResonance;
  //double q0Val = moZero(m0Val,m1Val,m2Val);
  //complex <GDouble> fpl = F(pVal,m_l);
  //complex <GDouble> fqs = F(qVal,m_s);
  //complex <GDouble> bwVal(1,1);
  //bwVal = BW(mVal,qVal,m0Val,q0Val,Gamma0,m_l);
  //preFactor *=  sqrt(2*m_l+1)*sqrt(2*m_s+1)*fpl*fqs*bwVal;
  //double preFactorR = real(preFactor);
  //double preFactorI = imag(preFactor);
  //if (preFactorR <= 0){
  //}else{
  //if (preFactorR >= 0){
  //  }else{//ASDF
  //cout<<"preFactorR = "<<preFactorR<<endl;
  //cout<<"bwVal = "<<bwVal<<endl;
  //cout<<"GammaVal = "<<Gamma(mVal,qVal,m0Val,q0Val,Gamma0,m_l)<<endl;
  //complex <GDouble> GammaVal = Gamma0*m0Val*qVal*pow(F(qVal,m_l),2)/(mVal*q0Val*pow(F(q0Val,m_l),2));
  //cout<<"GammaVal2 = "<<GammaVal<<endl;
  //cout<<"Gamma0, m0Val, qVal, mVal, q0Val = "<<Gamma0<<" , "<<m0Val<<" , "<<qVal<<" , "<<q0Val<<endl;
  //  }
  // }
  //}
  //if (m_r == 3){//f1(1420)
  //double Gamma0 = 0.0545;
  //double m0Val = 1.4263;
  //double mVal = mResonance;
  //double q0Val = moZero(m0Val,m1Val,m2Val);
  //complex <GDouble> fpl = F(pVal,m_l);
  //complex <GDouble> fqs = F(qVal,m_s);
  //complex <GDouble> bwVal(1,1);
  //bwVal = BW(mVal,qVal,m0Val,q0Val,Gamma0,m_l);
  //preFactor *=  sqrt(2*m_l+1)*sqrt(2*m_s+1)*fpl*fqs*bwVal;
  //}
  
  amplitude = 0;

  if (m_e == -1 || m_e == +1) {

    int neg = -1;
    double thetaM =1/sqrt(2.0);
    if (m_0 == 0) thetaM = 0.5;
    if (m_0 < 0) thetaM = 0;
    int parity = pow(neg,j_1 + j_2 + 1);
    
    for (int helicity = -j_2; helicity <= j_2; helicity++){//helicity
      if (abs(helicity) <= j_0){ 
	GDouble hel_amp = clebschGordan(j_1,j_2,0,helicity,j_0,helicity);
	amplitude += preFactor*hel_amp*thetaM*
	  ( conj(wignerD( j_0,  m_0, helicity, cosThetaGJ, phiGJ ))  
	    +        m_e*parity*pow(neg,j_0 - m_0)*conj(wignerD( j_0, -m_0, helicity, cosThetaGJ, phiGJ ))
	    )*conj(wignerD( j_2, helicity, 0, cosTheta_H, phi_H ));
      }
    }
  } else {

    for (int helicity = -j_2; helicity <= j_2; helicity++){//helicity
      if (abs(helicity) <= j_0){
	GDouble hel_amp = clebschGordan(j_1,j_2,0,helicity,j_0,helicity);
	amplitude += preFactor*hel_amp*
	  conj(wignerD( j_0, m_0, helicity, cosThetaGJ, phiGJ ))*
	  conj(wignerD( j_2, helicity, 0, cosTheta_H, phi_H ));
      }
    }
    
    complex< GDouble > rotateY = polar((GDouble)1., (GDouble)(-1.*Phi));

    if (m_e ==  10) amplitude = real(amplitude * rotateY);
    if (m_e == -10) amplitude = imag(amplitude * rotateY);

  }


  return amplitude ;

}