#include "BreitWigner.h"

BreitWigner::BreitWigner( const vector< string >& args ) : UserAmplitude< BreitWigner >( args )
{
	assert( args.size() == 5 );
	
	m_mass0 = AmpParameter( args[0] );
	m_width0 = AmpParameter( args[1] );
	m_orbitL = atoi( args[2].c_str() );

	// Extract particle indices for product 1
	string locParticleIndexString = args[3];
	vector<int> locParticleIndices;
	while(true)
	{
		size_t locUnderscoreIndex = locParticleIndexString.find_first_of("_");
		string locSubString = locParticleIndexString.substr(0, locUnderscoreIndex);
		istringstream locTempStream(locSubString);
		int locIndex = -100;
		locTempStream >> locIndex;
		locParticleIndices.push_back(locIndex);
		if(locUnderscoreIndex == string::npos)
			break;
		locParticleIndexString = locParticleIndexString.substr(locUnderscoreIndex + 1);
	}
	dNumParticleIndices_Product1 = locParticleIndices.size();
	dParticleIndices_Product1 = new int[locParticleIndices.size()];
	for(size_t loc_i = 0; loc_i < locParticleIndices.size(); ++loc_i)
		dParticleIndices_Product1[loc_i] = locParticleIndices[loc_i];

	// Extract particle indices for product 2
	locParticleIndexString = args[4];
	locParticleIndices.clear();
	while(true)
	{
		size_t locUnderscoreIndex = locParticleIndexString.find_first_of("_");
		string locSubString = locParticleIndexString.substr(0, locUnderscoreIndex);
		istringstream locTempStream(locSubString);
		int locIndex = -100;
		locTempStream >> locIndex;
		locParticleIndices.push_back(locIndex);
		if(locUnderscoreIndex == string::npos)
			break;
		locParticleIndexString = locParticleIndexString.substr(locUnderscoreIndex + 1);
	}
	dNumParticleIndices_Product2 = locParticleIndices.size();
	dParticleIndices_Product2 = new int[locParticleIndices.size()];
	for(size_t loc_i = 0; loc_i < locParticleIndices.size(); ++loc_i)
		dParticleIndices_Product2[loc_i] = locParticleIndices[loc_i];

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

	// make sure the input variables look reasonable
	assert(m_orbitL >= 0);
}

complex< GDouble > BreitWigner::calcAmplitude( GDouble** pKin ) const
{
  TLorentzVector P1, P2, Ptot, Ptemp;
  
	// Get/Build Product 1
	for(int loc_i = 0; loc_i < dNumParticleIndices_Product1; ++loc_i)
	{
		int locIndex = dParticleIndices_Product1[loc_i];
		P1 += TLorentzVector(pKin[locIndex][1], pKin[locIndex][2], pKin[locIndex][3], pKin[locIndex][0]);
	}

	// Get/Build Product 2
	for(int loc_i = 0; loc_i < dNumParticleIndices_Product2; ++loc_i)
	{
		int locIndex = dParticleIndices_Product2[loc_i];
		P2 += TLorentzVector(pKin[locIndex][1], pKin[locIndex][2], pKin[locIndex][3], pKin[locIndex][0]);
	}

  Ptot = P1 + P2;
  GDouble mass  = Ptot.M();
  GDouble mass1 = P1.M();
  GDouble mass2 = P2.M();
  
  // assert positive breakup momenta     
  GDouble q0 = fabs( breakupMomentum(m_mass0, mass1, mass2) );
  GDouble q  = fabs( breakupMomentum(mass,    mass1, mass2) );
  
  GDouble F0 = barrierFactor(q0, m_orbitL);
  GDouble F  = barrierFactor(q,  m_orbitL);
  
  GDouble width = m_width0*(m_mass0/mass)*(q/q0)*((F*F)/(F0*F0));
  //GDouble width = m_width0;
  
  // this first factor just gets normalization right for BW's that have
  // no additional s-dependence from orbital L
  complex<GDouble> bwtop( sqrt( m_mass0 * m_width0 / 3.1416 ), 0.0 );
  
  complex<GDouble> bwbottom( ( m_mass0*m_mass0 - mass*mass ) ,
                           -1.0 * ( m_mass0 * width ) );
  
  return( F * bwtop / bwbottom );
}

void BreitWigner::updatePar( const AmpParameter& par ){
 
  // could do expensive calculations here on parameter updates
  
}