#include "HCA_TensorContractor.h"

//int locHelicity12 = locHelicity1 - locHelicity2;
void HCA_TensorContractor::Calc_Invariants_2BodyDecay(int locJ1, int locJ2, int locHelicity12, int locS12, int locL12, int locJ12, bool locPsi1IsPhotonFlag, bool locPsi2IsPhotonFlag, deque<HCA_ContractionResult>& locFinalAmplitudes)
{
	//build all wave functions & amplitude terms //PsiS12 * PsiL12 * PsiJ12
		//if p & epsilon or g~ are needed, will include them when contracting
	deque<AmplitudeTerm> locAmplitudeTerms;
	Build_AmplitudeTerms(locJ1, locJ2, locHelicity12, locS12, locL12, locJ12, locPsi1IsPhotonFlag, locPsi2IsPhotonFlag, locAmplitudeTerms);
	if(dDebugFlag)
	{
		cout << "# amplitude terms = " << locAmplitudeTerms.size() << endl;
		for(size_t loc_i = 0; loc_i < locAmplitudeTerms.size(); ++loc_i)
		{
			cout << "Amplitude Term " << loc_i << ":" << endl;
			Print_AmplitudeTermInfo(locAmplitudeTerms[loc_i]);
		}
	}
	if(locAmplitudeTerms.empty())
		return;

	//now, determine the number of contractions
		//each scheme will result in a new fit term
	deque<ContractionScheme> locContractionSchemes;
	Build_ContractionSchemes(locJ1, locJ2, locL12, locJ12, locContractionSchemes);
	if(dDebugFlag)
	{
		cout << "# contraction schemes = " << locContractionSchemes.size() << endl;
		for(size_t loc_i = 0; loc_i < locContractionSchemes.size(); ++loc_i)
		{
			cout << "Contraction Scheme " << loc_i << ":" << endl;
			Print_ContractionSchemeInfo(locContractionSchemes[loc_i]);
		}
	}
	if(locContractionSchemes.empty())
		return; //shouldn't be possible

	//finally, contract the tensors to get the final amplitude terms
	Contract_Tensors(locContractionSchemes, locAmplitudeTerms, locFinalAmplitudes);
	if(dDebugFlag)
		cout << "# final amplitudes = " << locFinalAmplitudes.size() << endl;

	//set remaining amplitude members
	for(size_t loc_i = 0; loc_i < locFinalAmplitudes.size(); ++loc_i)
	{
		locFinalAmplitudes[loc_i].dJ1 = locJ1;
		locFinalAmplitudes[loc_i].dJ2 = locJ2;
		locFinalAmplitudes[loc_i].dHelicity12 = locHelicity12;
		locFinalAmplitudes[loc_i].dS12 = locS12;
		locFinalAmplitudes[loc_i].dL12 = locL12;
		locFinalAmplitudes[loc_i].dJ12 = locJ12;
		locFinalAmplitudes[loc_i].dPsi1IsPhotonFlag = locPsi1IsPhotonFlag;
		locFinalAmplitudes[loc_i].dPsi2IsPhotonFlag = locPsi2IsPhotonFlag;
	}

	//print result to screen
	if(dDebugFlag)
	{
		cout << "# final amplitudes = " << locFinalAmplitudes.size() << endl;
		for(size_t loc_i = 0; loc_i < locFinalAmplitudes.size(); ++loc_i)
		{
			cout << "Final Amplitude " << loc_i << ":" << endl;
			Print_FinalAmplitude(locFinalAmplitudes[loc_i]);
		}
	}
}

void HCA_TensorContractor::Filter_Amplitudes(deque<HCA_ContractionResult>& locFinalAmplitudes)
{
	//see if any of the final amplitude terms are identical: if so, remove them
		//ONLY use this to recduce number of print statements
		//e.g. if two different contraction schemes gave the same answer for a given helicity, keep them both: they may give different answers for different helicities
	deque<HCA_ContractionResult>::iterator locIterator = locFinalAmplitudes.begin();
	while(locIterator != locFinalAmplitudes.end())
	{
		bool locDupeFoundFlag = false;
		HCA_ContractionResult& locFinalAmplitude1 = *locIterator;
		deque<HCA_ContractionResult>::iterator loc2ndIterator = locFinalAmplitudes.begin();
		for(; loc2ndIterator != locIterator; ++loc2ndIterator)
		{
			HCA_ContractionResult& locFinalAmplitude2 = *loc2ndIterator;
			if(locFinalAmplitude1 != locFinalAmplitude2)
				continue;
			locIterator = locFinalAmplitudes.erase(locIterator);
			locDupeFoundFlag = true;
			break;
		}
		if(!locDupeFoundFlag)
			++locIterator;
	}
}

void HCA_TensorContractor::Build_AmplitudeTerms(int locJ1, int locJ2, int locHelicity12, int locS12, int locL12, int locJ12, bool locPsi1IsPhotonFlag, bool locPsi2IsPhotonFlag, deque<AmplitudeTerm>& locAmplitudeTerms)
{
	if((abs(locHelicity12) > locS12) || (abs(locHelicity12) > locJ12))
		return;

	//PsiL12
	WaveFunction locPsiL12;
	Build_Psi(locL12, 0, false, false, locPsiL12); //returning false isn't possible
	double locFactorialL12 = Factorial(locL12);
	locPsiL12.dCoefficient *= sqrt(pow(2.0, locL12)*locFactorialL12*locFactorialL12/Factorial(2*locL12)); //SU 2008 Eq 4.10
	if(dDebugFlag)
	{
		cout << "Psi_L12:" << endl;
		Print_WaveFunctionInfo(locPsiL12);
	}

	//PsiJ12
	WaveFunction locPsiJ12;
	Build_Psi(locJ12, locHelicity12, true, false, locPsiJ12); //returning false isn't possible
	if(dDebugFlag)
	{
		cout << "Psi_J12:" << endl;
		Print_WaveFunctionInfo(locPsiJ12);
	}

	//PsiS12 = sum_helicities( clebsch * PsiJ1 * PsiJ2 )
		//Helicities for S12: H12 = H1 - H2: H2 = H1 - H12
	//loop over internal helicities //can range from +/- J1/J2
	for(int locS12Helicity1 = -1*locJ1; locS12Helicity1 <= locJ1; ++locS12Helicity1)
	{
		int locS12Helicity2 = locS12Helicity1 - locHelicity12; //internal helicity constraint
		if(abs(locS12Helicity2) > locJ2)
			continue; //invalid helicity

		if(dDebugFlag)
			cout << "Psi_J12 component helicities: " << locS12Helicity1 << ", " << locS12Helicity2 << endl;
		WaveFunction locPsiJ1;
		if(!Build_Psi(locJ1, locS12Helicity1, false, locPsi1IsPhotonFlag, locPsiJ1))
		{
			if(dDebugFlag)
				cout << "bad helicity for J1" << endl;
			continue; //bad helicity for photon
		}
		if(dDebugFlag)
		{
			cout << "Psi_J1:" << endl;
			Print_WaveFunctionInfo(locPsiJ1);
		}

		WaveFunction locPsiJ2;
		if(!Build_Psi(locJ2, -1*locS12Helicity2, false, locPsi2IsPhotonFlag, locPsiJ2))
		{
			if(dDebugFlag)
				cout << "bad helicity for J2" << endl;
			continue; //bad helicity for photon
		}
		if(dDebugFlag)
		{
			cout << "Psi_J2:" << endl;
			Print_WaveFunctionInfo(locPsiJ2);
		}

		//Clebsch-Gordan Coefficient
		//double clebschGordan(int ij1, int ij2, int im1, int im2, int ij, int im)
		double locPsiS12ClebschTerm = 1.0;
		if((locJ1 != 0) && (locJ2 != 0))
			locPsiS12ClebschTerm = clebschGordan(locJ1, locJ2, locS12Helicity1, -1*locS12Helicity2, locS12, locHelicity12);
		if(dDebugFlag)
			cout << "Clebsch: J1, H1, J2, -1*H2, S12, H12, Clebsch = " << locJ1 << ", " << locS12Helicity1 << ", " << locJ2 << ", " << -1*locS12Helicity2 << ", " << locS12 << ", " << locHelicity12 << ", " << locPsiS12ClebschTerm << endl << endl;

		deque<AmplitudeTerm> locNewAmplitudeTerms;
		Build_AmplitudeTerms(locPsiL12, locPsiJ12, locPsiJ1, locPsiJ2, locPsiS12ClebschTerm, locNewAmplitudeTerms);
		locAmplitudeTerms.insert(locAmplitudeTerms.end(), locNewAmplitudeTerms.begin(), locNewAmplitudeTerms.end());
	}
}

void HCA_TensorContractor::Build_AmplitudeTerms(WaveFunction& locPsiL12, WaveFunction& locPsiJ12, WaveFunction& locPsiJ1, WaveFunction& locPsiJ2, double locPsiS12ClebschTerm, deque<AmplitudeTerm>& locNewAmplitudeTerms)
{
	//create initial term
	AmplitudeTerm locAmplitudeTerm;
	double locCoefficient = locPsiS12ClebschTerm * locPsiJ1.dCoefficient * locPsiJ2.dCoefficient * locPsiL12.dCoefficient * locPsiJ12.dCoefficient;
	locAmplitudeTerm.dCoefficient = locCoefficient;
	//order is important!
	if(!locPsiJ1.dComponentPermutations.empty())
		locAmplitudeTerm.dComplexConjugateFlag[PsiJ1] = locPsiJ1.dComplexConjugateFlag;
	if(!locPsiJ2.dComponentPermutations.empty())
		locAmplitudeTerm.dComplexConjugateFlag[PsiJ2] = locPsiJ2.dComplexConjugateFlag;
	if(!locPsiL12.dComponentPermutations.empty())
		locAmplitudeTerm.dComplexConjugateFlag[PsiL12] = locPsiL12.dComplexConjugateFlag;
	if(!locPsiJ12.dComponentPermutations.empty())
		locAmplitudeTerm.dComplexConjugateFlag[PsiJ12] = locPsiJ12.dComplexConjugateFlag;

	deque<AmplitudeTerm> locOldAmplitudeTerms(1, locAmplitudeTerm);
	if(locPsiJ12.dRank > 0)
		locNewAmplitudeTerms.clear();
	else
		locNewAmplitudeTerms = locOldAmplitudeTerms;

	//PsiJ12
	for(size_t loc_i = 0; loc_i < locPsiJ12.dComponentPermutations.size(); ++loc_i)
	{
		ComponentPermutation& locComponentPermutation = locPsiJ12.dComponentPermutations[loc_i];
		for(size_t loc_j = 0; loc_j < locOldAmplitudeTerms.size(); ++loc_j)
		{
			AmplitudeTerm locNewAmplitudeTerm = AmplitudeTerm(locOldAmplitudeTerms[loc_j]);
			locNewAmplitudeTerm.dComponentHelicities[PsiJ12] = locComponentPermutation.dComponentHelicities;
			locNewAmplitudeTerm.dCoefficient *= locComponentPermutation.dCoefficient;
			locNewAmplitudeTerms.push_back(locNewAmplitudeTerm);
		}
	}
	locOldAmplitudeTerms = locNewAmplitudeTerms;

	//PsiL12
	if(locPsiL12.dRank > 0)
		locNewAmplitudeTerms.clear();
	for(size_t loc_i = 0; loc_i < locPsiL12.dComponentPermutations.size(); ++loc_i)
	{
		ComponentPermutation& locComponentPermutation = locPsiL12.dComponentPermutations[loc_i];
		for(size_t loc_j = 0; loc_j < locOldAmplitudeTerms.size(); ++loc_j)
		{
			AmplitudeTerm locNewAmplitudeTerm = AmplitudeTerm(locOldAmplitudeTerms[loc_j]);
			locNewAmplitudeTerm.dComponentHelicities[PsiL12] = locComponentPermutation.dComponentHelicities;
			locNewAmplitudeTerm.dCoefficient *= locComponentPermutation.dCoefficient;
			locNewAmplitudeTerms.push_back(locNewAmplitudeTerm);
		}
	}
	locOldAmplitudeTerms = locNewAmplitudeTerms;

	//PsiJ1
	if(locPsiJ1.dRank > 0)
		locNewAmplitudeTerms.clear();
	for(size_t loc_i = 0; loc_i < locPsiJ1.dComponentPermutations.size(); ++loc_i)
	{
		ComponentPermutation& locComponentPermutation = locPsiJ1.dComponentPermutations[loc_i];
		for(size_t loc_j = 0; loc_j < locOldAmplitudeTerms.size(); ++loc_j)
		{
			AmplitudeTerm locNewAmplitudeTerm = AmplitudeTerm(locOldAmplitudeTerms[loc_j]);
			locNewAmplitudeTerm.dComponentHelicities[PsiJ1] = locComponentPermutation.dComponentHelicities;
			locNewAmplitudeTerm.dCoefficient *= locComponentPermutation.dCoefficient;
			locNewAmplitudeTerms.push_back(locNewAmplitudeTerm);
		}
	}
	locOldAmplitudeTerms = locNewAmplitudeTerms;

	//PsiJ2
	if(locPsiJ2.dRank > 0)
		locNewAmplitudeTerms.clear();
	for(size_t loc_i = 0; loc_i < locPsiJ2.dComponentPermutations.size(); ++loc_i)
	{
		ComponentPermutation& locComponentPermutation = locPsiJ2.dComponentPermutations[loc_i];
		for(size_t loc_j = 0; loc_j < locOldAmplitudeTerms.size(); ++loc_j)
		{
			AmplitudeTerm locNewAmplitudeTerm = AmplitudeTerm(locOldAmplitudeTerms[loc_j]);
			locNewAmplitudeTerm.dComponentHelicities[PsiJ2] = locComponentPermutation.dComponentHelicities;
			locNewAmplitudeTerm.dCoefficient *= locComponentPermutation.dCoefficient;
			locNewAmplitudeTerms.push_back(locNewAmplitudeTerm);
		}
	}
}

void HCA_TensorContractor::Build_ContractionSchemes(int locJ1, int locJ2, int locL12, int locJ12, deque<ContractionScheme>& locContractionSchemes)
{
	int locN = (locJ1 + locJ2 + locL12 - locJ12) % 2;
	if((locN == 0) && (locJ2 == 0))
	{
		ContractionScheme locContractionScheme;
		locContractionScheme.dNumContractions[PsiJ12][PsiL12] = (locJ12 + locL12 - locJ1)/2;
		locContractionScheme.dNumContractions[PsiJ12][PsiJ1] = (locJ12 - locL12 + locJ1)/2;
		locContractionScheme.dNumContractions[PsiL12][PsiJ1] = (locL12 + locJ1 - locJ12)/2;
		locContractionSchemes.push_back(locContractionScheme);
	}
	else if((locN == 0) && (locJ1 == 0))
	{
		ContractionScheme locContractionScheme;
		locContractionScheme.dNumContractions[PsiJ12][PsiL12] = (locJ12 + locL12 - locJ2)/2;
		locContractionScheme.dNumContractions[PsiJ12][PsiJ2] = (locJ12 - locL12 + locJ2)/2;
		locContractionScheme.dNumContractions[PsiL12][PsiJ2] = (locL12 + locJ2 - locJ12)/2;
		locContractionSchemes.push_back(locContractionScheme);
	}
	else if((locN == 1) && (locJ2 == 0))
	{
		ContractionScheme locContractionScheme;
		//order is important!!
		locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiJ1);
		locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiL12);
		locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiJ12);
		locContractionScheme.dNumContractions[PsiJ12][PsiL12] = (locJ12 + locL12 - locJ1 - 1)/2;
		locContractionScheme.dNumContractions[PsiJ12][PsiJ1] = (locJ12 - locL12 + locJ1 - 1)/2;
		locContractionScheme.dNumContractions[PsiL12][PsiJ1] = (locL12 + locJ1 - locJ12 - 1)/2;
		locContractionSchemes.push_back(locContractionScheme);
	}
	else if((locN == 1) && (locJ1 == 0))
	{
		ContractionScheme locContractionScheme;
		//order is important!!
		locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiJ2);
		locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiL12);
		locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiJ12);
		locContractionScheme.dNumContractions[PsiJ12][PsiL12] = (locJ12 + locL12 - locJ2 - 1)/2;
		locContractionScheme.dNumContractions[PsiJ12][PsiJ2] = (locJ12 - locL12 + locJ2 - 1)/2;
		locContractionScheme.dNumContractions[PsiL12][PsiJ2] = (locL12 + locJ2 - locJ12 - 1)/2;
		locContractionSchemes.push_back(locContractionScheme);
	}
	else if(locN == 0)
	{
		//neither J1 or J2 has spin 0
		for(int locN_J12_J1 = 0; locN_J12_J1 <= min(locJ12, locJ1); ++locN_J12_J1)
		{
			for(int locN_J12_J2 = 0; locN_J12_J2 <= min(locJ12 - locN_J12_J1, locJ2); ++locN_J12_J2)
			{
				int locN_J12_L12 = locJ12 - locN_J12_J1 - locN_J12_J2;
				if(locN_J12_L12 > locL12)
					continue; //invalid combo

				int locN_J1_J2 = (locJ1 + locJ2 - locL12 - locJ12)/2 + locN_J12_L12;
				if(locN_J1_J2 < 0)
					continue; //if too large for J1/J2, will fail below checks

				int locN_L12_J1 = locJ1 - locN_J12_J1 - locN_J1_J2;
				if(locN_L12_J1 < 0)
					continue; //if too large for L12, will catch below (cannot be too large for J1)

				int locN_L12_J2 = locJ2 - locN_J12_J2 - locN_J1_J2;
				if(locN_L12_J2 < 0)
					continue; //if too large for L12, will catch below (cannot be too large for J2)

				if((locN_L12_J1 + locN_L12_J2 + locN_J12_L12) > locL12)
					continue; //invalid combo

				ContractionScheme locContractionScheme;
				locContractionScheme.dNumContractions[PsiJ12][PsiJ1] = locN_J12_J1;
				locContractionScheme.dNumContractions[PsiJ12][PsiJ2] = locN_J12_J2;
				locContractionScheme.dNumContractions[PsiJ12][PsiL12] = locN_J12_L12;
				locContractionScheme.dNumContractions[PsiJ1][PsiJ2] = locN_J1_J2;
				locContractionScheme.dNumContractions[PsiL12][PsiJ1] = locN_L12_J1;
				locContractionScheme.dNumContractions[PsiL12][PsiJ2] = locN_L12_J2;
				locContractionSchemes.push_back(locContractionScheme);
			}
		}
	}
	else
	{
		//n = 1, neither J1 or J2 has spin 0
			//one of the contractions with the Levi-Civita tensor is 0, the others are 1
		for(int loc_i = 0; loc_i < 4; ++loc_i)
		{
			int locN_J12_LC = (loc_i == 0) ? 0 : 1;
			if((locJ12 == 0) && (locN_J12_LC == 1))
				continue; //don't waste time

			int locN_L12_LC = (loc_i == 1) ? 0 : 1;
			if((locL12 == 0) && (locN_L12_LC == 1))
				continue; //don't waste time

			int locN_J1_LC = (loc_i == 2) ? 0 : 1;
			int locN_J2_LC = (loc_i == 3) ? 0 : 1;

			for(int locN_J12_J1 = 0; locN_J12_J1 <= min(locJ12 - locN_J12_LC, locJ1 - locN_J1_LC); ++locN_J12_J1)
			{
				for(int locN_J12_J2 = 0; locN_J12_J2 <= min(locJ12 - locN_J12_J1 - locN_J12_LC, locJ2 - locN_J2_LC); ++locN_J12_J2)
				{
					int locN_J12_L12 = locJ12 - locN_J12_J1 - locN_J12_J2 - locN_J12_LC;
					if(locN_J12_L12 > locL12)
						continue; //invalid combo

					int locTermLCSum = locN_J12_LC + locN_L12_LC - locN_J1_LC - locN_J2_LC;
					int locN_J1_J2 = (locJ1 + locJ2 - locL12 - locJ12 + locTermLCSum)/2 + locN_J12_L12;
					if(locN_J1_J2 < 0)
						continue; //if too large for J1/J2, will fail below checks

					int locN_L12_J1 = locJ1 - locN_J12_J1 - locN_J1_J2 - locN_J1_LC;
					if(locN_L12_J1 < 0)
						continue; //if too large for L12, will catch below (cannot be too large for J1)

					int locN_L12_J2 = locJ2 - locN_J12_J2 - locN_J1_J2 - locN_J2_LC;
					if(locN_L12_J2 < 0)
						continue; //if too large for L12, will catch below (cannot be too large for J2)

					if((locN_L12_J1 + locN_L12_J2 + locN_J12_L12 + locN_L12_LC) > locL12)
						continue; //invalid combo

					ContractionScheme locContractionScheme;
					//order is important!
					if(locN_J1_LC == 1)
						locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiJ1);
					if(locN_J2_LC == 1)
						locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiJ2);
					if(locN_L12_LC == 1)
						locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiL12);
					if(locN_J12_LC == 1)
						locContractionScheme.dLeviCivitaContractionTensors.push_back(PsiJ12);
					locContractionScheme.dNumContractions[PsiJ12][PsiJ1] = locN_J12_J1;
					locContractionScheme.dNumContractions[PsiJ12][PsiJ2] = locN_J12_J2;
					locContractionScheme.dNumContractions[PsiJ12][PsiL12] = locN_J12_L12;
					locContractionScheme.dNumContractions[PsiJ1][PsiJ2] = locN_J1_J2;
					locContractionScheme.dNumContractions[PsiL12][PsiJ1] = locN_L12_J1;
					locContractionScheme.dNumContractions[PsiL12][PsiJ2] = locN_L12_J2;
					locContractionSchemes.push_back(locContractionScheme);
				}
			}
		}
	}
}

void HCA_TensorContractor::Contract_Tensors(deque<ContractionScheme>& locContractionSchemes, deque<AmplitudeTerm>& locAmplitudeTerms, deque<HCA_ContractionResult>& locFinalAmplitudes)
{
	//ok, all contraction schemes are built.  apply them on all amplitude terms. 
	for(size_t loc_i = 0; loc_i < locContractionSchemes.size(); ++loc_i)
	{
		ContractionScheme& locContractionScheme = locContractionSchemes[loc_i];

		HCA_ContractionResult locFinalAmplitude;
		locFinalAmplitude.dContractionScheme = locContractionScheme;
		locFinalAmplitude.dHasWOverWThreshold = (!locContractionScheme.dLeviCivitaContractionTensors.empty()); //will be true of all non-zero terms

		//loop over amplitude terms, building final amplitude
		for(size_t loc_j = 0; loc_j < locAmplitudeTerms.size(); ++loc_j)
		{
			if(dDebugFlag)
				cout << "Contract Tensors for Amplitude Term " << loc_j << endl;
			Contract_Tensors(locContractionScheme, locAmplitudeTerms[loc_j], locFinalAmplitude);
			if(dDebugFlag)
				cout << endl;
		}

		//contraction scheme finished: locFinalAmplitude is final for this fit term
		//loop through, eliminate terms with coefficient = 0 (e.g. cancelling terms)
		map<pair<int, int>, double>::iterator locIterator = locFinalAmplitude.dGammaTerms.begin();
		while(locIterator != locFinalAmplitude.dGammaTerms.end())
		{
			if(fabs(locIterator->second) > 1.0E-9)
				++locIterator;
			else
				locFinalAmplitude.dGammaTerms.erase(locIterator++);
		}

		//set params of final amplitude
		locFinalAmplitude.dContractionIndex = loc_i;

		//include all contraction schemes, even if they are zero. This is needed to keep track of when to exclude them
		locFinalAmplitudes.push_back(locFinalAmplitude);
	}
}

void HCA_TensorContractor::Contract_Tensors(ContractionScheme& locContractionScheme, AmplitudeTerm& locAmplitudeTerm, HCA_ContractionResult& locFinalAmplitude)
{
	//values:
	int locNumGammaTerms_PsiJ1 = 0;
	int locNumGammaTerms_PsiJ2 = 0;
	double locAmplitudeCoefficient = locAmplitudeTerm.dCoefficient;

	//contract!!!
	map<TensorType, int> locCurrentIndices;
	locCurrentIndices[PsiJ12] = 0;
	locCurrentIndices[PsiL12] = 0;
	locCurrentIndices[PsiJ1] = 0;
	locCurrentIndices[PsiJ2] = 0;

	//handle levi-civita
	if(!locContractionScheme.dLeviCivitaContractionTensors.empty())
	{
		TensorType locTensorType1 = locContractionScheme.dLeviCivitaContractionTensors[0];
		TensorType locTensorType2 = locContractionScheme.dLeviCivitaContractionTensors[1];
		TensorType locTensorType3 = locContractionScheme.dLeviCivitaContractionTensors[2];

		//grab first index of each tensor
		int locHelicity1 = locAmplitudeTerm.dComponentHelicities[locTensorType1][0];
		int locHelicity2 = locAmplitudeTerm.dComponentHelicities[locTensorType2][0];
		int locHelicity3 = locAmplitudeTerm.dComponentHelicities[locTensorType3][0];
		bool locIsConjugateFlag1 = locAmplitudeTerm.dComplexConjugateFlag[locTensorType1];
		bool locIsConjugateFlag2 = locAmplitudeTerm.dComplexConjugateFlag[locTensorType2];
		bool locIsConjugateFlag3 = locAmplitudeTerm.dComplexConjugateFlag[locTensorType3];

		//return +1, -1, 0 for gamma*w/w_th, -1*gamma*w/w_th, and 0
			//gamma = gamma of helicity = 0 wave function
		int locLeviCivitaResult = Contract_LeviCivita(locHelicity1, locHelicity2, locHelicity3, locIsConjugateFlag1, locIsConjugateFlag2, locIsConjugateFlag3);
		if(dDebugFlag)
		{
			cout << "Levi-Civita Helicities, Conjugate Flags = " << locHelicity1 << ", " << locHelicity2 << ", " << locHelicity3 << ", " << locIsConjugateFlag1 << ", " << locIsConjugateFlag2 << ", " << locIsConjugateFlag3 << endl;
			cout << "Levi-Civita Result = " << locLeviCivitaResult << endl;
		}
		if(locLeviCivitaResult == 0)
			return; //term is zero

		locAmplitudeCoefficient *= double(locLeviCivitaResult);
		if(((locHelicity1 == 0) && (locTensorType1 == PsiJ1)) || ((locHelicity2 == 0) && (locTensorType2 == PsiJ1)) || ((locHelicity3 == 0) && (locTensorType3 == PsiJ1)))
			++locNumGammaTerms_PsiJ1;
		else if(((locHelicity1 == 0) && (locTensorType1 == PsiJ2)) || ((locHelicity2 == 0) && (locTensorType2 == PsiJ2)) || ((locHelicity3 == 0) && (locTensorType3 == PsiJ2)))
			++locNumGammaTerms_PsiJ2;

		++locCurrentIndices[locTensorType1];
		++locCurrentIndices[locTensorType2];
		++locCurrentIndices[locTensorType3];
	}

	//handle metric contractions
	map<TensorType, map<TensorType, int> >::iterator locIterator = locContractionScheme.dNumContractions.begin();
	for(; locIterator != locContractionScheme.dNumContractions.end(); ++locIterator)
	{
		TensorType locTensorType1 = locIterator->first;
		map<TensorType, int>& loc2ndContractionTypeMap = locIterator->second;
		map<TensorType, int>::iterator loc2ndIterator = loc2ndContractionTypeMap.begin();
		for(; loc2ndIterator != loc2ndContractionTypeMap.end(); ++loc2ndIterator)
		{
			TensorType locTensorType2 = loc2ndIterator->first;
			int locNumContractions = loc2ndIterator->second;
			if(dDebugFlag)
				cout << "Tensor1, Tensor2, #contractions = " << locTensorType1 << ", " << locTensorType2 << ", " << locNumContractions << endl;
			//loop over num contractions over these tensors
			for(int loc_k = 0; loc_k < locNumContractions; ++loc_k)
			{
				int locHelicity1 = locAmplitudeTerm.dComponentHelicities[locTensorType1][locCurrentIndices[locTensorType1]];
				int locHelicity2 = locAmplitudeTerm.dComponentHelicities[locTensorType2][locCurrentIndices[locTensorType2]];
				bool locComplexConjugateFlag1 = locAmplitudeTerm.dComplexConjugateFlag[locTensorType1];
				bool locComplexConjugateFlag2 = locAmplitudeTerm.dComplexConjugateFlag[locTensorType2];
				if(dDebugFlag)
					cout << "Contract Metric: Tensor1, Tensor2, H1, H2, ConjFlag1, ConjFlag2: " << locTensorType1 << ", " << locTensorType2 << ", " << locHelicity1 << ", " << locHelicity2 << ", " << locComplexConjugateFlag1 << ", " << locComplexConjugateFlag2 << endl;

				if((locHelicity1 == 0) && (locHelicity2 == 0))
				{
					//coefficient is 1
					if((locTensorType1 == PsiJ1) || (locTensorType2 == PsiJ1))
						++locNumGammaTerms_PsiJ1;
					if((locTensorType1 == PsiJ2) || (locTensorType2 == PsiJ2))
						++locNumGammaTerms_PsiJ2;
					if(dDebugFlag)
						cout << "Contract Metric: #gamma_1/2 terms = " << locNumGammaTerms_PsiJ1 << ", " << locNumGammaTerms_PsiJ2 << endl;
				}
				else
				{
					bool locOnlyOneIsComplexConjugateFlag = (locComplexConjugateFlag1 != locComplexConjugateFlag2);
					double locCoefficient = 0.0;
					bool locTermIsZeroFlag = !Contract_Components_Metric(locHelicity1, locHelicity2, locOnlyOneIsComplexConjugateFlag, locCoefficient);
					if(dDebugFlag)
						cout << "Contract Metric: OnlyOneConjFlag, IsZeroFlag, Coefficient: " << locOnlyOneIsComplexConjugateFlag << ", " << locTermIsZeroFlag << ", " << locCoefficient << endl;
					if(locTermIsZeroFlag)
						return; //term is zero
					locAmplitudeCoefficient *= locCoefficient;
				}
				++locCurrentIndices[locTensorType1];
				++locCurrentIndices[locTensorType2];
			}
		}
	}

	//add terms together (that have the same number of gammas)
	pair<int, int> locNumGammaTermsPair(locNumGammaTerms_PsiJ1, locNumGammaTerms_PsiJ2);
	if(locFinalAmplitude.dGammaTerms.find(locNumGammaTermsPair) == locFinalAmplitude.dGammaTerms.end())
		locFinalAmplitude.dGammaTerms[locNumGammaTermsPair] = locAmplitudeCoefficient;
	else
		locFinalAmplitude.dGammaTerms[locNumGammaTermsPair] += locAmplitudeCoefficient;
}

bool HCA_TensorContractor::Build_Psi(int locJ, int locHelicity, bool locComplexConjugateFlag, bool locIsPhoton, WaveFunction& locWaveFunction)
{
	if(locIsPhoton && (locHelicity == 0))
		return false;

	if(locHelicity < 0)
	{
		//psi(J, -m) = (-1)^m * psi^*(J, m) //SU 2008 Eq 4.6
		bool locGoodFlag = Build_Psi(locJ, -1*locHelicity, !locComplexConjugateFlag, locIsPhoton, locWaveFunction);
		if(!locGoodFlag)
			return false;
		locWaveFunction.dHelicity = locHelicity;
		if(((-1*locHelicity) % 2) == 1)
			locWaveFunction.dCoefficient *= -1.0;
		return true;
	}

	locWaveFunction.dSpin = locJ;
	locWaveFunction.dRank = locJ;
	locWaveFunction.dHelicity = locHelicity;
	locWaveFunction.dComplexConjugateFlag = locComplexConjugateFlag;

	//don't waste time if J = 0
	if(locJ == 0)
	{
		locWaveFunction.dComplexConjugateFlag = false;
		locWaveFunction.dCoefficient = 1.0;
		locWaveFunction.dComponentPermutations.clear();
		return true;
	}

	//don't waste time if J = 1
	if(locJ == 1)
	{
		locWaveFunction.dCoefficient = 1.0;
		ComponentPermutation locComponentPermutation;

		locComponentPermutation.dCoefficient = 1.0;
		locComponentPermutation.dComponentHelicities = deque<int>(1, locHelicity);

		locWaveFunction.dComponentPermutations.push_back(locComponentPermutation);
		return true;
	}

	double locAJTerm = Calc_AJTerm(locJ, locHelicity);
	locWaveFunction.dCoefficient = locAJTerm;
	if(locIsPhoton)
		locWaveFunction.dCoefficient *= sqrt(4.0*3.14159265359);

	// Loop over the possible number of Helicity = 0 Components for Psi: locNumHZeroComponents
	int locJMinusHelicity = locJ - locHelicity;
	int locMinNumHZeroComponents = ((locJMinusHelicity % 2) == 0) ? 0 : 1;
	for(int locNumHZeroComponents = locMinNumHZeroComponents; locNumHZeroComponents <= locJMinusHelicity; locNumHZeroComponents += 2)
	{
		double locHZeroTermCoefficient = pow(2.0, double(locNumHZeroComponents)/2.0); //use kappa to be more efficient (SU 2008 Eq 4.1)

		int locNumHPlusComponents = (locJ + locHelicity - locNumHZeroComponents)/2;
		int locNumHMinusComponents = locNumHPlusComponents - locHelicity;

		map<int, int> locNumHelicityComponents;
		locNumHelicityComponents[1] = locNumHPlusComponents;
		locNumHelicityComponents[0] = locNumHZeroComponents;
		locNumHelicityComponents[-1] = locNumHMinusComponents;

		//build all permutations
		deque<ComponentPermutation> locPermutations;
		Build_Permutations(locJ, locNumHelicityComponents, locPermutations, locHZeroTermCoefficient);
		locWaveFunction.dComponentPermutations.insert(locWaveFunction.dComponentPermutations.end(), locPermutations.begin(), locPermutations.end());
	}

	return true;
}

/*
Cross products:
	//helicity = +1: -1/sqrt(2) * (1, i, 0)
	//helicity = -1: 1/sqrt(2) * (1, -i, 0)
	//helicity = 0: (0, 0, gamma)
	//helicity = *+1: -1/sqrt(2) * (1, -i, 0) = - helicity -1
	//helicity = *-1: 1/sqrt(2) * (1, i, 0) = - helicity 1


//a x b = (a2b3 - a3b2, a3b1 - a1b3, a1b2 - a2b1)
//b x a = - a x b
//a x a = 0

//0: 
+1 x +1
0 x 0
-1 x -1
1* x 1*
0* x 0*
-1* x -1*
0 x 0*
0* x 0
*+1 x -1
-1 x *+1
*-1 x +1
+1 x *-1

//(0, 0, i):
+1 x -1 // -1/2 * (0, 0, -i - i)
*+1 x +1
-1 x *-1
*-1 x *+1

//(0, 0, -i):
-1 x +1
+1 x *+1
*-1 x -1
*+1 x *-1

//1/sqrt(2) * (-i*gamma, gamma, 0):
+1 x 0
+1 x 0*
0 x *-1
0* x *-1

//1/sqrt(2) * (i*gamma, -gamma, 0)
0 x +1
0* x +1
*-1 x 0
*-1 x 0*

//1/sqrt(2) * (-i*gamma, -gamma, 0)
-1 x 0
-1 x 0*
0 x *+1
0* x *+1

//1/sqrt(2) * (i*gamma, gamma, 0)
0 x -1
0* x -1
*+1 x 0
*+1 x *0
*/

/*
//Now, dot product with these, excluding the 0-cross cases:
	//helicity = +1: -1/sqrt(2) * (1, i, 0)
	//helicity = -1: 1/sqrt(2) * (1, -i, 0)
	//helicity = 0: (0, 0, gamma)
	//helicity = *+1: -1/sqrt(2) * (1, -i, 0) = - helicity -1
	//helicity = *-1: 1/sqrt(2) * (1, i, 0) = - helicity 1


//X cdot (0, 0, i) = 
0: +1, -1, *+1, *-1
i*gamma: 0, 0*

//X cdot (0, 0, -i) = 
0: +1, -1, *+1, *-1
-i*gamma: 0, 0*

//X cdot 1/sqrt(2) * (-i*gamma, gamma, 0) = 
0: 0, 0*, +1, *-1
-i*gamma: -1
i*gamma: *+1

//X cdot 1/sqrt(2) * (i*gamma, -gamma, 0) = 
0: 0, 0*, +1, *-1
i*gamma: -1
-i*gamma: *+1

//X cdot 1/sqrt(2) * (-i*gamma, -gamma, 0) = 
0: 0, 0*, -1, *+1
i*gamma: +1 //-1/2 * (-i*gamma -gamma*i) = i*gamma
-i*gamma = *-1

//X cdot 1/sqrt(2) * (i*gamma, gamma, 0) = 
0: 0, 0*, -1, *+1
-i*gamma: +1 //-1/2 * (i*gamma + i*gamma) = -i*gamma
i*gamma: *-1
*/

/*
So, the combined results for bcd:

//if (0, 0, i):
i*gamma:
0/0*, +1, -1
0/0*, *+1, +1
0/0*, -1, *-1
0/0*, *-1, *+1

//if (0, 0, -i): 
-i*gamma:
0/0*, -1, +1
0/0*, +1, *+1
0/0*, *-1, -1
0/0*, *+1, *-1

if 1/sqrt(2) * (-i*gamma, gamma, 0): 
-i*gamma: 
-1, +1, 0/0*
-1, 0/0*, *-1
i*gamma:
*+1, +1, 0/0*
*+1, 0/0*, *-1

//if 1/sqrt(2) * (i*gamma, -gamma, 0):
i*gamma:
-1, 0/0*, +1
-1, *-1 x 0/0*
-i*gamma:
*+1, 0/0*, +1
*+1, *-1 x 0/0*

//if 1/sqrt(2) * (-i*gamma, -gamma, 0):
i*gamma: 
+1, -1 x 0/0*
+1, 0/0* x *+1
-i*gamma:
*-1, -1 x 0/0*
*-1, 0/0* x *+1

//if 1/sqrt(2) * (i*gamma, gamma, 0):
-i*gamma:
+1, 0/0*, -1
+1, *+1 x 0/0*
i*gamma:
*-1, 0/0*, -1
*-1, *+1 x 0/0*
*/

/*
Organizing the combined results for bcd:

//i*gamma:
0/0*, +1, -1
0/0*, *+1, +1
0/0*, -1, *-1
0/0*, *-1, *+1
*+1, +1, 0/0*
*+1, 0/0*, *-1
-1, 0/0*, +1
-1, *-1 x 0/0*
+1, -1 x 0/0*
+1, 0/0* x *+1
*-1, 0/0*, -1
*-1, *+1 x 0/0*

//-i*gamma:
0/0*, -1, +1
0/0*, +1, *+1
0/0*, *-1, -1
0/0*, *+1, *-1
-1, +1, 0/0*
-1, 0/0*, *-1
*+1, 0/0*, +1
*+1, *-1 x 0/0*
*-1, -1 x 0/0*
*-1, 0/0* x *+1
+1, 0/0*, -1
+1, *+1 x 0/0*
*/

/*
Summarizing the combined results for bcd:
//gamma = gamma of helicity = 0 wave function
//0 can also be 0*

//i*gamma:
0, +1, -1
0, *+1, +1
0, -1, *-1
0, *-1, *+1
*+1, +1, 0
*+1, 0, *-1
-1, 0, +1
-1, *-1, 0
+1, -1, 0
+1, 0, *+1
*-1, 0, -1
*-1, *+1, 0

//-i*gamma:
0, -1, +1
0, +1, *+1
0, *-1, -1
0, *+1, *-1
-1, +1, 0
-1, 0, *-1
*+1, 0, +1
*+1, *-1, 0
*-1, -1, 0
*-1, 0, *+1
+1, 0, -1
+1, *+1, 0
*/

int HCA_TensorContractor::Contract_LeviCivita(int locHelicity1, int locHelicity2, int locHelicity3, bool locIsConjugateFlag1, bool locIsConjugateFlag2, bool locIsConjugateFlag3)
{
	//return +1, -1, 0 for gamma*w/w_th, -1*gamma*w/w_th, and 0
		//gamma = gamma of helicity = 0 wave function

	//Helicity_0* = Helicity_0
	//Helicity_+1* = -1*Helicity_-1
	//Helicity_-1* = -1*Helicity_+1
	int locSignMultiplier = 1;
	if((locIsConjugateFlag1) && (locHelicity1 != 0))
	{
		locHelicity1 *= -1;
		locSignMultiplier *= -1;
	}
	if((locIsConjugateFlag2) && (locHelicity2 != 0))
	{
		locHelicity2 *= -1;
		locSignMultiplier *= -1;
	}
	if((locIsConjugateFlag3) && (locHelicity3 != 0))
	{
		locHelicity3 *= -1;
		locSignMultiplier *= -1;
	}

	if((locHelicity1 == locHelicity2) || (locHelicity1 == locHelicity3) || (locHelicity2 == locHelicity3))
		return 0;

	//Determine cyclical order
	if(((locHelicity1 + 1) == locHelicity2) || ((locHelicity1 - 2) == locHelicity2)) 
		return locSignMultiplier; //0, +1, -1  OR  -1, 0, +1  OR  +1, -1, 0
	else
		return -1*locSignMultiplier; //0, -1, +1  OR  -1, +1, 0  OR  +1, 0, -1
}

//dot products:
//contract across the metric //returned value is always real
bool HCA_TensorContractor::Contract_Components_Metric(int locHelicity1, int locHelicity2, bool locOnlyOneIsComplexConjugateFlag, double& locCoefficient)
{
	//helicity = +1: -1/sqrt(2) * (0, 1, i, 0)
	//helicity = -1: 1/sqrt(2) * (0, 1, -i, 0)
	//helicity = 0: (gamma*beta, 0, 0, gamma)

	//metric: delta_ij, except metric_00 = 0
	//metric * helicity_deque = helicity_deque with time-component = 0
	//so helicity_deque1^a * metric_ab * helicity_deque2^b = helicity_deque1 * helicity_deque2 (dot product, but ignoring time components)
		//ignoring the time component only effects (helicity1=0 * helicity2=0): otherwise term was 0 anyway

/*
	//all cases
	if(!loc2IsComplexConjugateFlag)
	{
		if((locHelicity1 == 1) && (locHelicity2 == 1))
			return 0.0; //[- 1/sqrt(2) * (0, 1, i, 0)]^2 = 1/2 * (1 - 1) = 0
		else if((locHelicity1 == 1) && (locHelicity2 == -1))
			return -1.0; //[- 1/sqrt(2) * (0, 1, i, 0)] * [1/sqrt(2) * (0, 1, -i, 0)] = -1/2 * (1 + 1) = -1
		else if((locHelicity1 == 1) && (locHelicity2 == 0))
			return 0.0; //[- 1/sqrt(2) * (0, 1, i, 0)] * [(gamma*beta, 0, 0, gamma)] = 0

		else if((locHelicity == 0) && (locHelicity2 == 1))
			return 0.0; //[(gamma*beta, 0, 0, gamma)] * [- 1/sqrt(2) * (0, 1, i, 0)] = 0
		else if((locHelicity == 0) && (locHelicity2 == -1))
			return 0.0; //[(gamma*beta, 0, 0, gamma)] * [1/sqrt(2) * (0, 1, -i, 0)] = 0
		else if((locHelicity == 0) && (locHelicity2 == 0))
			return locGamma1*locGamma2; //[(gamma*beta, 0, 0, gamma)] * [(gamma*beta, 0, 0, gamma)] = gamma1*gamma2

		else if((locHelicity1 == -1) && (locHelicity2 == 1))
			return -1.0; //[1/sqrt(2) * (0, 1, -i, 0)] * [-1/sqrt(2) * (0, 1, i, 0)] = -1/2 * (1 + 1) = -1
		else if((locHelicity1 == -1) && (locHelicity2 == -1))
			return 0.0; //[1/sqrt(2) * (0, 1, -i, 0)]^2 = 1/2 * (1 - 1) = 0
		else if((locHelicity1 == -1) && (locHelicity2 == 0))
			return 0.0; //[1/sqrt(2) * (0, 1, -i, 0)] * [(gamma*beta, 0, 0, gamma)] = 0
	}
	else
	{
		if((locHelicity1 == 1) && (locHelicity2 == 1))
			return 1.0; //[- 1/sqrt(2) * (0, 1, i, 0)] * [- 1/sqrt(2) * (0, 1, -i, 0)] = 1/2 * (1 + 1) = 1
		else if((locHelicity1 == 1) && (locHelicity2 == -1))
			return 0.0; //[- 1/sqrt(2) * (0, 1, i, 0)] * [1/sqrt(2) * (0, 1, i, 0)] = -1/2 * (1 - 1) = 0
		else if((locHelicity1 == 1) && (locHelicity2 == 0))
			return 0.0; //[- 1/sqrt(2) * (0, 1, i, 0)] * [(gamma*beta, 0, 0, gamma)] = 0

		else if((locHelicity == 0) && (locHelicity2 == 1))
			return 0.0; //[(gamma*beta, 0, 0, gamma)] * [- 1/sqrt(2) * (0, 1, -i, 0)] = 0
		else if((locHelicity == 0) && (locHelicity2 == -1))
			return 0.0; //[(gamma*beta, 0, 0, gamma)] * [1/sqrt(2) * (0, 1, i, 0)] = 0
		else if((locHelicity == 0) && (locHelicity2 == 0))
			return locGamma1*locGamma2; //[(gamma*beta, 0, 0, gamma)] * [(gamma*beta, 0, 0, gamma)] = gamma1*gamma2

		else if((locHelicity1 == -1) && (locHelicity2 == 1))
			return 0.0; //[1/sqrt(2) * (0, 1, -i, 0)] * [-1/sqrt(2) * (0, 1, -i, 0)] = -1/2 * (1 - 1) = 0
		else if((locHelicity1 == -1) && (locHelicity2 == -1))
			return 1.0; //[1/sqrt(2) * (0, 1, -i, 0)] * [1/sqrt(2) * (0, 1, i, 0)] = 1/2 * (1 + 1) = 1
		else if((locHelicity1 == -1) && (locHelicity2 == 0))
			return 0.0; //[1/sqrt(2) * (0, 1, -i, 0)] * [(gamma*beta, 0, 0, gamma)] = 0
	}

	//helicity = +1: -1/sqrt(2) * (0, 1, i, 0)
	//helicity = -1: 1/sqrt(2) * (0, 1, -i, 0)

//both complex:
	//helicity = +1: -1/sqrt(2) * (0, 1, -i, 0)
	//helicity = -1: 1/sqrt(2) * (0, 1, i, 0)
	//+1 * +1 = -1 * -1 = 0
	//+1 * -1 = -1.0
*/

//	if((locHelicity == 0) && (locHelicity2 == 0))
//		return locGamma1*locGamma2; //[(gamma*beta, 0, 0, gamma)] * [(gamma*beta, 0, 0, gamma)] = gamma1*gamma2 (across metric with time-component = 0)

	if(!locOnlyOneIsComplexConjugateFlag)
	{
		//both are complex or neither is complex
		if(locHelicity1 != -1*locHelicity2)
			return false;
		locCoefficient = -1.0;
		return true;
	}

	//one & only one of them is complex
	if(locHelicity1 != locHelicity2)
		return false;
	locCoefficient = 1.0;
	return true;
}

void HCA_TensorContractor::Build_Permutations(int locRank, map<int, int>& locNumHelicityComponents, deque<ComponentPermutation>& locPermutations, double locCoefficient)
{
	int locIndex = 0;
	deque<int> locResumeAtValues(locRank, -1); //each ranges from -1 to 1: if is 2, then have finished it
	//each one has 3 possibilities

	map<int, int> locNumUsedHelicityComponents;
	locNumUsedHelicityComponents[-1] = 0;
	locNumUsedHelicityComponents[0] = 0;
	locNumUsedHelicityComponents[1] = 0;

	deque<int> locPermutation;
	while(true) //index-level
	{
		//see if permutation is complete
		if(locIndex == locRank)
		{
			//save it
			ComponentPermutation locComponentPermutation;
			locComponentPermutation.dCoefficient = locCoefficient;
			locComponentPermutation.dComponentHelicities = locPermutation;
			locPermutations.push_back(locComponentPermutation);

			//go through previous indices, resetting them, until reach an index where we can resume
			if(!Handle_PermutationDecursion(locIndex, locPermutation, locResumeAtValues, locNumUsedHelicityComponents))
				return;
			continue; //continue on previous index
		}

		//get helicity value
		int locHelicity = locResumeAtValues[locIndex];
		bool locValueIsGoodFlag = true;
		while(locNumUsedHelicityComponents[locHelicity] >= locNumHelicityComponents[locHelicity])
		{
			++locResumeAtValues[locIndex];
			locHelicity = locResumeAtValues[locIndex];
			if(locHelicity == 2)
			{
				locValueIsGoodFlag = false;
				break; //for the helicities currently in locPermutation, all possible permutations have been found
			}
		}

		if(!locValueIsGoodFlag)
		{
			//go through previous indices, resetting them, until reach an index where we can resume
			if(!Handle_PermutationDecursion(locIndex, locPermutation, locResumeAtValues, locNumUsedHelicityComponents))
				return;
			continue; //continue on previous index
		}

		//helicity value is good, save it
		locPermutation.push_back(locHelicity);
		++locNumUsedHelicityComponents[locHelicity];
		++locResumeAtValues[locIndex]; //next time we come back to this index, start at a larger helicity

		//advance to next rank
		++locIndex;
	}
}

bool HCA_TensorContractor::Handle_PermutationDecursion(int& locIndex, deque<int>& locPermutation, deque<int>& locResumeAtValues, map<int, int>& locNumUsedHelicityComponents)
{
	//go through previous indices, resetting them, until reach an index where we can resume
	do
	{
		if(locPermutation.empty())
			return false; //we're done

		//reset this index
		if(locIndex < int(locResumeAtValues.size()))
			locResumeAtValues[locIndex] = -1;

		//re-do (resume-on) the previous rank
		int locLastHelicity = locPermutation.back();
		locPermutation.pop_back();
		--locNumUsedHelicityComponents[locLastHelicity];
		--locIndex;
	}
	while((locResumeAtValues[locIndex] == 2) || (locIndex == int(locResumeAtValues.size() - 1)));
	//last while condition: no other possible values for the last index (it was forced): resuming it will fail: back track in advance

	return true;
}

void HCA_TensorContractor::Print_WaveFunctionInfo(WaveFunction& locPsi)
{
	cout << "Rank, Spin, Helicity = " << locPsi.dRank << ", " << locPsi.dSpin << ", " << locPsi.dHelicity << endl;
	cout << "Coefficient, conjugate flag = " << locPsi.dCoefficient << ", " << locPsi.dComplexConjugateFlag << endl;
	for(size_t loc_i = 0; loc_i < locPsi.dComponentPermutations.size(); ++loc_i)
	{
		ComponentPermutation& locPermutation = locPsi.dComponentPermutations[loc_i];
		cout << "Permutation " << loc_i << ": Coefficient, Helicities = " << locPermutation.dCoefficient << "; ";
		for(size_t loc_j = 0; loc_j < locPermutation.dComponentHelicities.size(); ++loc_j)
		{
			cout << locPermutation.dComponentHelicities[loc_j];
			if(loc_j != (locPermutation.dComponentHelicities.size() - 1))
				cout << ", ";
		}
		cout << endl;
	}
	cout << endl;
}

void HCA_TensorContractor::Print_AmplitudeTermInfo(AmplitudeTerm& locAmplitudeTerm)
{
	cout << "Amplitude coefficient = " << locAmplitudeTerm.dCoefficient << endl;
	map<TensorType, deque<int> >::iterator locHelicityIterator = locAmplitudeTerm.dComponentHelicities.begin();
	map<TensorType, bool>::iterator locConjugateIterator = locAmplitudeTerm.dComplexConjugateFlag.begin();
	for(; locHelicityIterator != locAmplitudeTerm.dComponentHelicities.end(); ++locHelicityIterator, ++locConjugateIterator)
	{
		cout << "Tensor type (2x), conjugate flag: " << locHelicityIterator->first << ", " << locConjugateIterator->first << ", " << locConjugateIterator->second << endl;
		cout << "Component Helicities: ";
		for(size_t loc_i = 0; loc_i < locHelicityIterator->second.size(); ++loc_i)
		{
			cout << locHelicityIterator->second[loc_i];
			if(loc_i != (locHelicityIterator->second.size() - 1))
				cout << ", ";
		}
		cout << endl;
	}
	cout << endl;
}

void HCA_TensorContractor::Print_ContractionSchemeInfo(ContractionScheme& locContractionScheme)
{
	cout << "Levi-Civita Tensor terms: ";
	for(size_t loc_i = 0; loc_i < locContractionScheme.dLeviCivitaContractionTensors.size(); ++loc_i)
	{
		cout << locContractionScheme.dLeviCivitaContractionTensors[loc_i];
		if(loc_i != (locContractionScheme.dLeviCivitaContractionTensors.size() - 1))
			cout << ", ";
	}
	cout << endl;

	map<TensorType, map<TensorType, int> >::iterator locTensorIterator1 = locContractionScheme.dNumContractions.begin();
	for(; locTensorIterator1 != locContractionScheme.dNumContractions.end(); ++locTensorIterator1)
	{
		cout << "# Contractions by tensor for " << locTensorIterator1->first << ": ";
		map<TensorType, int>::iterator locTensorIterator2 = locTensorIterator1->second.begin();
		for(; locTensorIterator2 != locTensorIterator1->second.end(); ++locTensorIterator2)
		{
			if(locTensorIterator2 != locTensorIterator1->second.begin())
				cout << "; ";
			cout << "Type: " << locTensorIterator2->first << ", Num: " << locTensorIterator2->second;
		}
		cout << endl;
	}
	cout << endl;
}

void HCA_TensorContractor::Print_FinalAmplitude(HCA_ContractionResult& locFinalAmplitude)
{
	if(locFinalAmplitude.dHasWOverWThreshold)
		cout << "(w/w_th)*[";
	map<pair<int, int>, double>::iterator locIterator = locFinalAmplitude.dGammaTerms.begin();
	for(; locIterator != locFinalAmplitude.dGammaTerms.end(); ++locIterator)
	{
		if(locIterator != locFinalAmplitude.dGammaTerms.begin())
			cout << " + ";
		cout << locIterator->second << " * Gamma1^" << locIterator->first.first << " * Gamma2^" << locIterator->first.second;
	}
	if(locFinalAmplitude.dHasWOverWThreshold)
		cout << "]";
	cout << endl;
}