#ifndef _HCA_TensorContractor_
#define _HCA_TensorContractor_

#include "AMPTOOLS_AMPS/clebschGordan.h"
#include <iostream>
#include <deque>
#include <map>
#include <math.h>
#include <stdlib.h>
#include <fstream>

#include "HCA_Components.h"
#include "HCA_ContractionResult.h"

using namespace std;

class HCA_TensorContractor
{
	//builds the tensor contraction terms
	public:
		HCA_TensorContractor() : dDebugFlag(false) {}

		int Factorial(int locN);
		double Calc_AJTerm(int locJ, int locM);

		void Calc_Invariants_2BodyDecay(int locJ1, int locJ2, int locHelicity12, int locS12, int locL12, int locJ12, bool locPsi1IsPhotonFlag, bool locPsi2IsPhotonFlag, deque<HCA_ContractionResult>& locFinalAmplitudes);
		void Filter_Amplitudes(deque<HCA_ContractionResult>& locFinalAmplitudes);

		void Print_WaveFunctionInfo(WaveFunction& locPsi);
		void Print_AmplitudeTermInfo(AmplitudeTerm& locAmplitudeTerm);
		void Print_ContractionSchemeInfo(ContractionScheme& locContractionScheme);
		void Print_FinalAmplitude(HCA_ContractionResult& locFinalAmplitude);

		bool dDebugFlag;

	private:
		void Build_AmplitudeTerms(int locJ1, int locJ2, int locHelicity12, int locS12, int locL12, int locJ12, bool locPsi1IsPhotonFlag, bool locPsi2IsPhotonFlag, deque<AmplitudeTerm>& locAmplitudeTerms);
		void Build_ContractionSchemes(int locJ1, int locJ2, int locL12, int locJ12, deque<ContractionScheme>& locContractionSchemes);
		void Contract_Tensors(deque<ContractionScheme>& locContractionSchemes, deque<AmplitudeTerm>& locAmplitudeTerms, deque<HCA_ContractionResult>& locFinalAmplitudes);
		void Contract_Tensors(ContractionScheme& locContractionScheme, AmplitudeTerm& locAmplitudeTerm, HCA_ContractionResult& locFinalAmplitude);

		bool Build_Psi(int locJ, int locHelicity, bool locComplexConjugateFlag, bool locIsPhoton, WaveFunction& locWaveFunction);
		void Build_AmplitudeTerms(WaveFunction& locPsiL12, WaveFunction& locPsiJ12, WaveFunction& locPsiJ1, WaveFunction& locPsiJ2, double locPsiS12ClebschTerm, deque<AmplitudeTerm>& locNewAmplitudeTerms);
		int Contract_LeviCivita(int locHelicity1, int locHelicity2, int locHelicity3, bool locIsConjugateFlag1, bool locIsConjugateFlag2, bool locIsConjugateFlag3);
		bool Contract_Components_Metric(int locHelicity1, int locHelicity2, bool locOnlyOneIsComplexConjugateFlag, double& locCoefficient);
		void Build_Permutations(int locRank, map<int, int>& locNumHelicityComponents, deque<ComponentPermutation>& locPermutations, double locCoefficient);
		bool Handle_PermutationDecursion(int& locIndex, deque<int>& locPermutation, deque<int>& locResumeAtValues, map<int, int>& locNumUsedHelicityComponents);
};

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

//term in front of poldeque definition //assumes locM >= 0!!
inline double HCA_TensorContractor::Calc_AJTerm(int locJ, int locM)
{
	//use kappa to be more efficient (SU 2008 Eq 4.1)
	return sqrt(double(Factorial(locJ + locM)*Factorial(locJ - locM))/double(Factorial(2*locJ)));
}

#endif //_HCA_TensorContractor_