#ifndef _HCAmplitude_Step_
#define _HCAmplitude_Step_

#include <string>
#include <complex>
#include <vector>
#include <sstream>
#include <stdlib.h>

#include <TLorentzVector.h>
#include <TLorentzRotation.h>
#include <TVector3.h>

#include "IUAmpTools/Amplitude.h"
#include "IUAmpTools/AmpParameter.h"
#include "IUAmpTools/UserAmplitude.h"

#include "AMPTOOLS_AMPS/wignerD.h"
#include "AMPTOOLS_AMPS/clebschGordan.h"
#include "AMPTOOLS_AMPS/barrierFactor.h"

#include "HCA_TensorContraction/HCA_TensorContractor.h"

using namespace std;

class Kinematics;

class HCAmplitude_Step
{
	public:
		HCAmplitude_Step() : dStepDebugFlag(false) {}

		void Calculate_TensorContraction(void);

		complex<GDouble> calcAmplitude(GDouble** locKinematics) const;
		complex<GDouble> Calc_CovariantDecayAmplitude(TLorentzVector locP4_1_CM, TLorentzVector locP4_2_CM) const;

		GDouble Calc_BarrierFactor(GDouble locX, int locL) const;
		int Factorial(int locN, int locMinN = 2) const;

		bool dStepDebugFlag;

		//particle 4-vectors are sum of the particle indices set. if num == -1, is target
		int dNumParticleIndices_Product1;
		int* dParticleIndices_Product1;

		int dNumParticleIndices_Product2;
		int* dParticleIndices_Product2;

		int dNumParticleIndices_XsParent;
		int* dParticleIndices_XsParent;

		int dNumParticleIndices_YsParent;
		int* dParticleIndices_YsParent;

		bool dProduct1IsPhotonFlag;
		bool dProduct2IsPhotonFlag;

		int dJ1Times2;
		int dHelicity1Times2;

		int dJ2Times2;
		int dHelicity2Times2;

		int dJXTimes2;
		int dHelicityXTimes2;

		int dL12;
		int dS12Times2;

		int dContractionIndex;

		int dHelicity12Times2;
		GDouble dTargetMass;
		GDouble dClebschGordanCoefficient;

		HCA_TensorContractor dHCABuilder;
		HCA_ContractionResult dFinalAmplitude;
};

inline int HCAmplitude_Step::Factorial(int locN, int locMinN) const
{
	if(locN < 0)
		abort();
	if((locN == 1) || (locN == 0) || (locN < locMinN))
		return 1;
	return locN * Factorial(locN - 1);
}

#endif