#include "HCA_ConfigGenerator.h"

void HCA_ConfigGenerator::Generate(deque<NodeBase*>& locTransitionTrees)
{
	dAmplitudeCounter = 0;

	//loop over helicities & do tensor contractions
	deque<HCA_ContractionResult> locAllFinalAmplitudes;
	Generate_HelicityCouplingAmplitudes(locTransitionTrees, locAllFinalAmplitudes);

	if(dDebugFlag)
	{
		cout << "FINAL HCA's calc'ed, " << locAllFinalAmplitudes.size() << " final amplitudes." << endl;
		for(size_t loc_i = 0; loc_i < locAllFinalAmplitudes.size(); ++loc_i)
		{
			cout << "Final Amplitude " << loc_i + 1 << ":" << endl;
			dHelicityCouplingAmplitudeBuilder.Print_FinalAmplitude(locAllFinalAmplitudes[loc_i]);
			cout << endl;
		}
	}

	//Get names of exchange classes
	set<string> locExchangeClassNames;
	for(size_t loc_i = 0; loc_i < locTransitionTrees.size(); ++loc_i)
	{
		ExchangeNode* locExchangeNode = dynamic_cast<ExchangeNode*>(locTransitionTrees[loc_i]);
		if(locExchangeNode == NULL)
			continue;
		locExchangeClassNames.insert(locExchangeNode->dExchangeClassName);
	}

	//Config file
	Generate_ConfigFile(locTransitionTrees);
}

void HCA_ConfigGenerator::Print_Tree(NodeBase* locTreeBase, int locNumTabs)
{
	ResonanceNode* locResonanceNode = dynamic_cast<ResonanceNode*>(locTreeBase);
	if(locResonanceNode != NULL)
		Print_ResonanceNode(locResonanceNode, locNumTabs);

	ExchangeNode* locExchangeNode = dynamic_cast<ExchangeNode*>(locTreeBase);
	Particle* locPsi1 = locExchangeNode->dLeftFinalNode->dNodeParticle;
	Particle* locPsi2 = locExchangeNode->dRightFinalNode->dNodeParticle;

	ostringstream locCurrentLine;
	locCurrentLine << "# ";
	for(int loc_i = 0; loc_i < locNumTabs; ++loc_i)
		locCurrentLine << "\t";

	Particle_t locBeamPID = locExchangeNode->dBeamParticle->dPID;
	Particle_t locTargetPID = locExchangeNode->dTargetParticle->dPID;
	locCurrentLine << ParticleType(locBeamPID) << ", " << ParticleType(locTargetPID) << " --> ";

	if(locExchangeNode->dLeftFinalPID != Unknown)
		locCurrentLine << ParticleType(locExchangeNode->dLeftFinalPID) << ", ";
	else
		locCurrentLine << double(locPsi1->dJTimes2)/2.0 << "^" << locPsi1->dParity << " (q = " << locPsi1->dCharge << "), ";

	if(locExchangeNode->dRightFinalPID != Unknown)
		locCurrentLine << ParticleType(locExchangeNode->dRightFinalPID);
	else
		locCurrentLine << double(locPsi2->dJTimes2)/2.0 << "^" << locPsi2->dParity << " (q = " << locPsi2->dCharge << ")";

	dConfigFileStream << locCurrentLine.str() << endl;
	if(locPsi1->dIsResonanceFlag)
		Print_ResonanceNode(locExchangeNode->dLeftFinalNode, locNumTabs + 1);
	if(locPsi2->dIsResonanceFlag)
		Print_ResonanceNode(locExchangeNode->dRightFinalNode, locNumTabs + 1);
}

void HCA_ConfigGenerator::Print_ResonanceNode(ResonanceNode* locCurrentNode, int locNumTabs)
{
	ostringstream locCurrentLine;
	locCurrentLine << "# ";
	for(int loc_i = 0; loc_i < locNumTabs; ++loc_i)
		locCurrentLine << "\t";

	if(locCurrentNode->dIsProductionNodeFlag)
	{
		Particle_t locBeamPID = locCurrentNode->dLeftFinalNode->dNodeParticle->dPID;
		Particle_t locTargetPID = locCurrentNode->dRightFinalNode->dNodeParticle->dPID;
		int locJTimes2 = locCurrentNode->dNodeParticle->dJTimes2;
		int locParity = locCurrentNode->dNodeParticle->dParity;
		int locCharge = locCurrentNode->dNodeParticle->dCharge;
		locCurrentLine << ParticleType(locBeamPID) << ", " << ParticleType(locTargetPID) << " --> ";
		locCurrentLine << double(locJTimes2)/2.0 << "^" << locParity << " (q = " << locCharge << ")";
		locCurrentLine << " (S12x2 = " << locCurrentNode->dS12Times2 << ", L12 = " << locCurrentNode->dL12 << ")" << endl;

		dConfigFileStream << locCurrentLine.str();
		Print_ResonanceNode((ResonanceNode*)locCurrentNode->dPreviousDecayNode, locNumTabs + 1);
		return;
	}

	int locJTimes2 = locCurrentNode->dNodeParticle->dJTimes2;
	int locParity = locCurrentNode->dNodeParticle->dParity;
	int locCharge = locCurrentNode->dNodeParticle->dCharge;
	locCurrentLine << double(locJTimes2)/2.0 << "^" << locParity << " (q = " << locCharge << ") --> ";

	Particle* locPsi1 = locCurrentNode->dLeftFinalNode->dNodeParticle;
	if(locPsi1->dIsResonanceFlag)
		locCurrentLine << double(locPsi1->dJTimes2)/2.0 << "^" << locPsi1->dParity << " (q = " << locPsi1->dCharge << "), ";
	else
		locCurrentLine << ParticleType(locPsi1->dPID) << ", ";

	Particle* locPsi2 = locCurrentNode->dRightFinalNode->dNodeParticle;
	if(locPsi2->dIsResonanceFlag)
		locCurrentLine << double(locPsi2->dJTimes2)/2.0 << "^" << locPsi2->dParity << " (q = " << locPsi2->dCharge << ")";
	else
		locCurrentLine << ParticleType(locPsi2->dPID);

	locCurrentLine << " (S12x2 = " << locCurrentNode->dS12Times2 << ", L12 = " << locCurrentNode->dL12 << ")" << endl;

	dConfigFileStream << locCurrentLine.str();
	if(locPsi1->dIsResonanceFlag)
		Print_ResonanceNode(locCurrentNode->dLeftFinalNode, locNumTabs + 1);
	if(locPsi2->dIsResonanceFlag)
		Print_ResonanceNode(locCurrentNode->dRightFinalNode, locNumTabs + 1);
}

void HCA_ConfigGenerator::Generate_HelicityCouplingAmplitudes(deque<NodeBase*>& locTransitionTrees, deque<HCA_ContractionResult>& locAllFinalAmplitudes)
{
	for(size_t loc_i = 0; loc_i < locTransitionTrees.size(); ++loc_i)
	{
		//Initialize (Get starting node(s))
		deque<ResonanceNode*> locOtherUndecayedNodesToProcess;
		ResonanceNode* locCurrentNode = Get_StartingNodes(locTransitionTrees[loc_i], locOtherUndecayedNodesToProcess);
		if(locCurrentNode == NULL)
			continue;

		while(true)
		{
			Particle* locPsiX = locCurrentNode->dNodeParticle;
			Particle* locPsi1 = locCurrentNode->dLeftFinalNode->dNodeParticle;
			Particle* locPsi2 = locCurrentNode->dRightFinalNode->dNodeParticle;

			//loop over possible helicities
			int locMaxHelicity12Times2 = locPsi1->dJTimes2 + locPsi2->dJTimes2;
			if(locMaxHelicity12Times2 > locCurrentNode->dS12Times2)
				locMaxHelicity12Times2 = locCurrentNode->dS12Times2;
			if(locMaxHelicity12Times2 > locPsiX->dJTimes2)
				locMaxHelicity12Times2 = locPsiX->dJTimes2;
			for(int locHelicity12Times2 = -1*locMaxHelicity12Times2; locHelicity12Times2 <= locMaxHelicity12Times2; locHelicity12Times2 += 2)
			{
				//calculate helicity coupling amplitudes
				deque<HCA_ContractionResult> locNewFinalAmplitudes;
				dHelicityCouplingAmplitudeBuilder.Calc_Invariants_2BodyDecay(locPsi1->dJTimes2/2, locPsi2->dJTimes2/2, locHelicity12Times2/2, locCurrentNode->dS12Times2/2, locCurrentNode->dL12, locPsiX->dJTimes2/2, locPsi1->dPID == Gamma, locPsi2->dPID == Gamma, locNewFinalAmplitudes);
				locAllFinalAmplitudes.insert(locAllFinalAmplitudes.end(), locNewFinalAmplitudes.begin(), locNewFinalAmplitudes.end());

				locCurrentNode->dNumContractionSchemes = locNewFinalAmplitudes.size();
				locCurrentNode->dTensorAmplitudes[locHelicity12Times2] = locNewFinalAmplitudes;
			}

			if((!locPsi1->dIsResonanceFlag) && (!locPsi2->dIsResonanceFlag))
			{
				if(locOtherUndecayedNodesToProcess.empty())
					break; //finished this decay tree, go to the next one
				locCurrentNode = locOtherUndecayedNodesToProcess.back();
				locOtherUndecayedNodesToProcess.pop_back();
				continue;
			}
			else if(locPsi2->dIsResonanceFlag && (!locPsi1->dIsResonanceFlag))
			{
				locCurrentNode = locCurrentNode->dRightFinalNode;
				continue;
			}
			else if(locPsi1->dIsResonanceFlag && (!locPsi2->dIsResonanceFlag))
			{
				locCurrentNode = locCurrentNode->dLeftFinalNode;
				continue;
			}
			else //both are resonances
			{
				locCurrentNode = locCurrentNode->dLeftFinalNode;
				locOtherUndecayedNodesToProcess.push_back(locCurrentNode->dRightFinalNode);
				continue;
			}
		}
	}

	//Filter out duplicate amplitudes (no need to print to screen twice)
		//e.g. two decays that have identical quantum numbers
	dHelicityCouplingAmplitudeBuilder.Filter_Amplitudes(locAllFinalAmplitudes);
}

void HCA_ConfigGenerator::Generate_ConfigFile(deque<NodeBase*>& locTransitionTrees)
{
	dConfigFileStream.open(dConfigFileName.c_str());

	dConfigFileStream << "#####################################" << endl;
	dConfigFileStream << "####	THIS IS A CONFIG FILE	 ####" << endl;
	dConfigFileStream << "#####################################" << endl;
	dConfigFileStream << "##" << endl;
	dConfigFileStream << "##  Blank lines or lines beginning with a \"#\" are ignored." << endl;
	dConfigFileStream << "##" << endl;
	dConfigFileStream << "##  Double colons (::) are treated like a space." << endl;
	dConfigFileStream << "##     This is sometimes useful for grouping (for example," << endl;
	dConfigFileStream << "##     grouping strings like \"reaction::sum::amplitudeName\")" << endl;
	dConfigFileStream << "##" << endl;
	dConfigFileStream << "##  All non-comment lines must begin with one of the following keywords." << endl;
	dConfigFileStream << "##" << endl;
	dConfigFileStream << "##  (note:  <word> means necessary" << endl;
	dConfigFileStream << "##	    (word) means optional)" << endl;
	dConfigFileStream << "##" << endl;
	dConfigFileStream << "##  include	  <file>" << endl;
	dConfigFileStream << "##  define	  <word> (defn1) (defn2) (defn3) ..." << endl;
	dConfigFileStream << "##  fit 	  <fitname>" << endl;
	dConfigFileStream << "##  keyword	  <keyword> <min arguments> <max arguments>" << endl;
	dConfigFileStream << "##  reaction	  <reaction> <particle1> <particle2> (particle3) ..." << endl;
	dConfigFileStream << "##  data	  <reaction> <class> (arg1) (arg2) (arg3) ..." << endl;
	dConfigFileStream << "##  genmc	  <reaction> <class> (arg1) (arg2) (arg3) ..." << endl;
	dConfigFileStream << "##  accmc	  <reaction> <class> (arg1) (arg2) (arg3) ..." << endl;
	dConfigFileStream << "##  normintfile   <reaction> <file>" << endl;
	dConfigFileStream << "##  sum 	  <reaction> <sum> (sum2) (sum3) ..." << endl;
	dConfigFileStream << "##  amplitude	  <reaction> <sum> <amp> <class> (arg1) (arg2) ([par]) ... " << endl;
	dConfigFileStream << "##  initialize    <reaction> <sum> <amp> <\"events\"/\"polar\"/\"cartesian\">" << endl;
	dConfigFileStream << "##		    <value1> <value2> (\"fixed\"/\"real\")" << endl;
	dConfigFileStream << "##  scale	  <reaction> <sum> <amp> <value or [parameter]>" << endl;
	dConfigFileStream << "##  constrain	  <reaction1> <sum1> <amp1> <reaction2> <sum2> <amp2> ..." << endl;
	dConfigFileStream << "##  permute	  <reaction> <sum> <amp> <index1> <index2> ..." << endl;
	dConfigFileStream << "##  parameter	  <par> <value> (\"fixed\"/\"bounded\"/\"gaussian\")" << endl;
	dConfigFileStream << "##		    (lower/central) (upper/error)" << endl;
	dConfigFileStream << "##    DEPRECATED:" << endl;
	dConfigFileStream << "##  datafile	  <reaction> <file> (file2) (file3) ..." << endl;
	dConfigFileStream << "##  genmcfile	  <reaction> <file> (file2) (file3) ..." << endl;
	dConfigFileStream << "##  accmcfile	  <reaction> <file> (file2) (file3) ..." << endl;
	dConfigFileStream << "##" << endl;
	dConfigFileStream << "#####################################" << endl;
	dConfigFileStream << endl;

	//build reaction name
	for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
	{
		if(loc_i != 0)
			dReactionName += "_";
		dReactionName += ParticleType(dFinalStatePIDs[loc_i]);
	}
	dConfigFileStream << "reaction " << dReactionName << " " << ParticleType(dBeamPID) << "_0";
	for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
		dConfigFileStream << " " << ParticleType(dFinalStatePIDs[loc_i]) << "_" << loc_i + 1; //names are unique: won't permute! all permutations handled (already) manually
	dConfigFileStream << endl;

	//build fit name
	string locFitName = dFitName;
	dConfigFileStream << "fit " << locFitName << endl;
	dConfigFileStream << endl;

	//file sources
	dConfigFileStream << "genmc " << dReactionName << " ROOTDataReader " << dGenerateDataFileName << endl;
	dConfigFileStream << "accmc " << dReactionName << " ROOTDataReader " << dAcceptedDataFileName << endl;
	dConfigFileStream << "data " << dReactionName << " ROOTDataReader " << dExperimentDataFileName << endl;

	//normalization integral file name //add "input" afterwards to read-in the results rather than re-compute
	dConfigFileStream << "normintfile " << dReactionName << " " << dNormalizationIntegralFileName << endl;
	dConfigFileStream << endl;

	//build all possible permutations of external helicities
	deque<deque<int> > locHelicityTimes2Permutations;
	Build_HelicityPermutations(locHelicityTimes2Permutations);

	//build PID string
	deque<Particle_t> locPIDs;
	for(int loc_i = 0; loc_i < int(dFinalStateParticles.size() + 2); ++loc_i)
	{
		int locParticleIndex = loc_i - 1;
		if(loc_i == 0) //beam
			locParticleIndex = 0;
		else if(loc_i == 1) //target
			locParticleIndex = -1;
		Particle* locParticle = dReverseParticleIndexMap[locParticleIndex];
		locPIDs.push_back(locParticle->dPID);
	}
	ostringstream locPIDStream;
	locPIDStream << locPIDs[0];
	for(size_t loc_i = 1; loc_i < locPIDs.size(); ++loc_i)
		locPIDStream << "_" << locPIDs[loc_i];
	dPIDString = locPIDStream.str();

	//loop over helicity permutations & print sums
	deque<string> locSumNames;
	for(size_t loc_i = 0; loc_i < locHelicityTimes2Permutations.size(); ++loc_i)
	{
		deque<int>& locPermutation = locHelicityTimes2Permutations[loc_i];
		ostringstream locSumStream;
		locSumStream << "Hx2";
		for(size_t loc_j = 0; loc_j < locPermutation.size(); ++loc_j)
			locSumStream << "_" << locPermutation[loc_j];
		dConfigFileStream << "sum " << dReactionName << " " << locSumStream.str() << endl;
		locSumNames.push_back(locSumStream.str());
	}
	dConfigFileStream << endl;

	dConfigFileStream << "# HCAmplitude_Chain Arguments: external_helicities, pids, # steps, steps (JXx2 L12 S12x2 ContractionIndex ParticleX__Particle1__Particle2__...etc) other_name, other_args" << endl;
	dConfigFileStream << "#\t\tExternal_helicities & pids: by p4 array index, EXCEPT target inserted at index = 1 (with indices starting at 0)" << endl;
	dConfigFileStream << "#\t\tParticles: _-separated p4 array index locations" << endl;
	dConfigFileStream << "#\t\tOther: Name is class name (e.g. t-channel exchange), args are the arguments to pass into it. leave these blank to ignore them" << endl;
	dConfigFileStream << "# BreitWigner Arguments: MassX, WidthX, L12, P4Indices_1, P4Indices_2" << endl;
	dConfigFileStream << "# Uniform (Phase Space) Arguments: None" << endl;
	dConfigFileStream << endl;

	//include phase space amplitude
	if(dIncludePhaseSpaceFlag)
	{
		dConfigFileStream << "# Phase Space amplitudes:" << endl;
		string locFirstAmpFullName = dReactionName + string("::") + locSumNames[0] + string("::Amp0");
		dUniqueFullAmplitudeNames.push_back(locFirstAmpFullName);
		for(size_t loc_i = 0; loc_i < locSumNames.size(); ++loc_i)
		{
			string locAmpFullName = (loc_i == 0) ? locFirstAmpFullName : dReactionName + string("::") + locSumNames[loc_i] + string("::Amp0");
			dConfigFileStream << "amplitude " << locAmpFullName << " Uniform " << endl;
			if(loc_i != 0)
				dConfigFileStream << "constrain " << locAmpFullName << " " << locFirstAmpFullName << endl;
			dConfigFileStream << endl;
		}
		dAmplitudeCounter = 1;
		dConfigFileStream << endl;
	}

	//for a given J/P/L/S combo (all): constrain all params equal, regardless of int/ext-ernal helicity state
	for(size_t loc_i = 0; loc_i < locTransitionTrees.size(); ++loc_i)
	{
		//reset
		gAmplitudeNames.clear(); //new tree: clear: used for constraining fit params of different helicity combos for the same decay tree

		//Initialize (Get starting node(s))
		deque<ResonanceNode*> locOtherUndecayedNodesToProcess;
		ResonanceNode* locCurrentNode = Get_StartingNodes(locTransitionTrees[loc_i], locOtherUndecayedNodesToProcess);
		if(locCurrentNode == NULL)
			continue;

		//init
		gProductionNode = locTransitionTrees[loc_i];
		for(size_t loc_j = 0; loc_j < locSumNames.size(); ++loc_j)
		{
			//loop over external helicity states & contraction schemes: //loop through decay tree
			vector<string> locBWArgumentsVector, locAmplitudeArgumentsVector;
			vector<int> locContractionIndices;
			Loop_PrintingAmplitudes(locCurrentNode, locOtherUndecayedNodesToProcess, locSumNames[loc_j], locHelicityTimes2Permutations[loc_j], locAmplitudeArgumentsVector, locBWArgumentsVector, locContractionIndices);
		}
	}

	dConfigFileStream << "include " << dSeedFileName << endl;
	dConfigFileStream.close();

	ofstream locSeedFileStream;
	locSeedFileStream.open(dSeedFileName.c_str());

	//DEFINE FIT PARAMETERS
	if(!dFitParams.empty())
	{
		locSeedFileStream << "# FIT PARAMETERS" << endl;
		set<FitParamInfo*>::iterator locIterator = dFitParams.begin();
		for(; locIterator != dFitParams.end(); ++locIterator)
		{
			locSeedFileStream << "parameter " << (*locIterator)->dName << " " << (*locIterator)->dValue;
			if((*locIterator)->dType != "")
			{
				locSeedFileStream << " " << (*locIterator)->dType;
				if(((*locIterator)->dLowerBound < (*locIterator)->dUpperBound) || ((*locIterator)->dUpperBound >= 0.0))
					locSeedFileStream << " " << (*locIterator)->dLowerBound << " " << (*locIterator)->dUpperBound;
			}
			locSeedFileStream << endl;
		}
		locSeedFileStream << endl;
	}


	//INITIALIZE AMPLITUDES
	locSeedFileStream << "# INITIALIZE AMPLITUDES" << endl;
	locSeedFileStream << "# \"real\" is necessary to fix the arbitrary, overall phase (else fit is underconstrained)" << endl;
	for(size_t loc_i = 0; loc_i < dUniqueFullAmplitudeNames.size(); ++loc_i)
	{
		locSeedFileStream << "initialize " << dUniqueFullAmplitudeNames[loc_i] << " polar 1.0 0.0";
		if(loc_i == 0)
			locSeedFileStream << " real";
		locSeedFileStream << endl;
	}
	locSeedFileStream << endl;

	locSeedFileStream.close();
}

void HCA_ConfigGenerator::Build_HelicityPermutations(deque<deque<int> >& locHelicityTimes2Permutations)
{
	deque<int> locResumeAtValues;
	deque<int> locJTimes2s;
	deque<Particle_t> locPIDs;
	for(int loc_i = 0; loc_i < int(dFinalStateParticles.size() + 2); ++loc_i)
	{
		int locParticleIndex = loc_i - 1;
		if(loc_i == 0) //beam
			locParticleIndex = 0;
		else if(loc_i == 1) //target
			locParticleIndex = -1;
		Particle* locParticle = dReverseParticleIndexMap[locParticleIndex];

		locPIDs.push_back(locParticle->dPID);
		locJTimes2s.push_back(locParticle->dJTimes2);
		locResumeAtValues.push_back(-1*locJTimes2s.back());
	}

	int locParticleIndex = 0;
	deque<int> locHelicityTimes2s;
	locHelicityTimes2Permutations.clear();
	while(true)
	{
		if(locParticleIndex == int(locPIDs.size()))
		{
			//permutation finished, save results
			locHelicityTimes2Permutations.push_back(locHelicityTimes2s);
			//go through previous indices, resetting them, until reach an index where we can resume
			if(!Handle_HelicityDecursion(locParticleIndex, locHelicityTimes2s, locResumeAtValues, locJTimes2s))
				break;
			continue; //continue on previous index
		}

		int locJTimes2 = locJTimes2s[locParticleIndex];
		int locHelicityTimes2 = locResumeAtValues[locParticleIndex];
		if((locHelicityTimes2 == 0) && (locPIDs[locParticleIndex] == Gamma))
			locHelicityTimes2 += 2; //cannot be 0 for photon
		if(locHelicityTimes2 > locJTimes2)
		{
			if(!Handle_HelicityDecursion(locParticleIndex, locHelicityTimes2s, locResumeAtValues, locJTimes2s))
				break;
			continue;
		}
		locHelicityTimes2s.push_back(locHelicityTimes2);
		locResumeAtValues[locParticleIndex] = locHelicityTimes2 + 2;
		++locParticleIndex;
	}
}

bool HCA_ConfigGenerator::Handle_HelicityDecursion(int& locParticleIndex, deque<int>& locHelicityTimes2s, deque<int>& locResumeAtValues, deque<int>& locJTimes2s)
{
	//go through previous indices, resetting them, until reach an index where we can resume
	do
	{
		if(locHelicityTimes2s.empty())
			return false; //we're done

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

		//re-do (resume-on) the previous rank
		locHelicityTimes2s.pop_back();
		--locParticleIndex;
		if(locParticleIndex < 0)
			return false; //we're done
	}
	while(locResumeAtValues[locParticleIndex] > locJTimes2s[locParticleIndex]);

	return true;
}

bool HCA_ConfigGenerator::Loop_PrintingAmplitudes(ResonanceNode* locCurrentNode, deque<ResonanceNode*> locOtherUndecayedNodesToProcess, string locSumName, deque<int>& locHelicityTimes2s, vector<string> locAmplitudeArgumentsVector, vector<string> locBWArgumentsVector, vector<int> locContractionIndices)
{
	ResonanceNode* locPreviousNode = dynamic_cast<ResonanceNode*>(locCurrentNode->dPreviousDecayNode);

	Particle* locPsiX = locCurrentNode->dNodeParticle;
	Particle* locPsi1 = locCurrentNode->dLeftFinalNode->dNodeParticle;
	Particle* locPsi2 = locCurrentNode->dRightFinalNode->dNodeParticle;

	//Determine the min & max H12x2
	int locMinH12x2 = 0, locMaxH12x2 = 0;
	Calculate_PossibleHelicityRange(locCurrentNode, locHelicityTimes2s, locMinH12x2, locMaxH12x2);

	//need indices for 4 vectors
	ResonanceNode* locParentIndicesNode = locCurrentNode->dIsProductionNodeFlag ? locPreviousNode : locCurrentNode; //use products for parent
	string locIndices_Parent = Get_ProductIndices(locParentIndicesNode);
	string locIndices_Product1 = Get_ProductIndices(locCurrentNode->dLeftFinalNode);
	string locIndices_Product2 = Get_ProductIndices(locCurrentNode->dRightFinalNode);

	//breit-wigner: see if needed
	bool locIncludeBWFlag = ((!dNoBreitWignersFlag) && (!locCurrentNode->dIsProductionNodeFlag)); //don't save if no-BW's or on production node
	if(locIncludeBWFlag) //don't save if previous node was s-channel production node AND are binning in sqrt-s
		locIncludeBWFlag = (locPreviousNode == NULL) ? true : ((!locPreviousNode->dIsProductionNodeFlag) || (!dBinInSqrtSFlag));
	if(locIncludeBWFlag) //Also, only save it if it's recognized (if that requirement is in place)
		locIncludeBWFlag = locCurrentNode->dResonanceRecognizedFlag;

	//breit-wigner save
	vector<string> locNewBWArgumentsVector = locBWArgumentsVector;
	if(locIncludeBWFlag)
	{
		//Breit-Wigner Arguments: MassX, WidthX, L12, P4Indices_1, P4Indices_2
		ostringstream locBWArguments;
		if(locPsiX->dMassParamName != "")
			locBWArguments << "[" << locPsiX->dMassParamName << "] ";
		else
			locBWArguments << locPsiX->dMass << " ";

		if(locPsiX->dWidthParamName != "")
			locBWArguments << "[" << locPsiX->dWidthParamName << "] ";
		else
			locBWArguments << locPsiX->dWidth << " ";

		locBWArguments << locCurrentNode->dL12 << " " << locIndices_Product1 << " " << locIndices_Product2;
		locNewBWArgumentsVector.push_back(locBWArguments.str());
	}

	//loop over contraction schemes
	int locNumIndependentAmplitudeTerms = locCurrentNode->dNumContractionSchemes;
	int locNumContractionSchemesUsed = 0;
	for(int locContractionIndex = 0; locContractionIndex < locNumIndependentAmplitudeTerms; ++locContractionIndex)
	{
		//Now, check to see if the amplitude is zero for all possible helicities with this contraction scheme
		//For a given H12x2, amplitude is zero if either:
			//abs(H12x2) > J12x2    OR    abs(H12x2) > S12x2    OR    Tensor result from this contraction scheme is zero
			//contraction RESULTS are helicity dependent, but contraction SCHEMES are not. 1 fit parameter per contraction SCHEME. 
		bool locAllTermsAreZeroFlag = true; //will negate if != 0 found
		for(int locH12x2 = locMinH12x2; locH12x2 <= locMaxH12x2; locH12x2 += 2)
		{
			//Contraction Scheme & helicity chosen
			if(Check_IsNodeAmplitudeZero(locCurrentNode, locH12x2, locContractionIndex))
				continue; //Is Zero

			//There is at least one non-zero term for this resonance: amplitude is good so far!
			locAllTermsAreZeroFlag = false;
			break;
		}
		if(locAllTermsAreZeroFlag)
			continue; //all terms for all possible helicities are zero for this contraction scheme: does not contribute: no fit parameter for this combination

		if(dLimitToOneContractionSchemeFlag && (locNumContractionSchemesUsed != 0))
			break;
		++locNumContractionSchemesUsed;

		vector<int> locNewContractionIndices = locContractionIndices;
		locNewContractionIndices.push_back(locContractionIndex);

		//build amp arguments string
		//(JXx2 L12 S12x2 ContractionIndex ParticleX__Particle1__Particle2__...etc)
		ostringstream locAmplitudeArguments;
		locAmplitudeArguments << locPsiX->dJTimes2 << " " << locCurrentNode->dL12 << " " << locCurrentNode->dS12Times2 << " " << locContractionIndex << " ";
		locAmplitudeArguments << locIndices_Parent << "__" << locIndices_Product1 << "__" << locIndices_Product2;

		vector<string> locNewAmplitudeArgumentsVector = locAmplitudeArgumentsVector;
		locNewAmplitudeArgumentsVector.push_back(locAmplitudeArguments.str());

		if((!locPsi1->dIsResonanceFlag) && (!locPsi2->dIsResonanceFlag))
		{
			if(locOtherUndecayedNodesToProcess.empty())
			{
				//finished navigating this decay tree for this helicity combination

				//get the amplitude index:
					//see if it is unique, or if need to constrain against another amplitude (in the same tree with same contractions but different helicities)
				map<vector<int>, pair<string, int> >::iterator locAmpIndexIterator = gAmplitudeNames.find(locNewContractionIndices);
				bool locConstrainToPreviousFlag = (locAmpIndexIterator != gAmplitudeNames.end());
				int locThisAmplitudeIndex = 0;
				if(locConstrainToPreviousFlag) //constrain
					locThisAmplitudeIndex = locAmpIndexIterator->second.second;
				else
				{
					//unique
					locThisAmplitudeIndex = dAmplitudeCounter;
					gAmplitudeNames[locNewContractionIndices] = pair<string, int>(locSumName, locThisAmplitudeIndex);
					++dAmplitudeCounter;
				}

				//build full amplitude name
				ostringstream locFullAmplitudeName;
				locFullAmplitudeName << dReactionName << "::" << locSumName << "::Amp" << locThisAmplitudeIndex;
				if(!locConstrainToPreviousFlag) //unique
				{
					dUniqueFullAmplitudeNames.push_back(locFullAmplitudeName.str());
					Print_Tree(gProductionNode);
				}

				//args: external_helicities pids, # steps, steps (JXx2 L12 S12x2 ContractionIndex ParticleX__Particle1__Particle2__...etc), other_name, other_args
				dConfigFileStream << "amplitude " << locFullAmplitudeName.str() << " HCAmplitude_Chain ";
				dConfigFileStream << locSumName << " " << dPIDString << " " << locNewAmplitudeArgumentsVector.size();
				for(size_t loc_i = 0; loc_i < locNewAmplitudeArgumentsVector.size(); ++loc_i)
					dConfigFileStream << " " << locNewAmplitudeArgumentsVector[loc_i];

				//exchange node (e.g. t-channel)
				ExchangeNode* locExchangeNode = dynamic_cast<ExchangeNode*>(gProductionNode);
				if(locExchangeNode != NULL)
				{
					if(locExchangeNode->dFitParams.empty()) //doesn't contain a fit parameter
						dConfigFileStream << " " << locExchangeNode->dExchangeClassName << " " << locSumName << " " << locExchangeNode->dClassConfigArguments << endl;
					else // contains a fit parameter: make separate
						dConfigFileStream << endl << "amplitude " << locFullAmplitudeName.str() << " " << locExchangeNode->dExchangeClassName << " " << locSumName << " " << locExchangeNode->dClassConfigArguments << endl;
				}
				else
					dConfigFileStream << endl;

				//breit-wigners
				for(size_t loc_i = 0; loc_i < locNewBWArgumentsVector.size(); ++loc_i)
				{
					dConfigFileStream << "amplitude " << locFullAmplitudeName.str();
					dConfigFileStream << " BreitWigner " << locNewBWArgumentsVector[loc_i] << endl;
				}

				//constrain if necessary
				if(locConstrainToPreviousFlag)
				{
					dConfigFileStream << "constrain " << locFullAmplitudeName.str();
					dConfigFileStream << " " << dReactionName << "::" << locAmpIndexIterator->second.first << "::Amp" << locAmpIndexIterator->second.second << endl;
				}
				dConfigFileStream << endl;
				continue;
			}
			deque<ResonanceNode*> locNewOtherUndecayedNodesToProcess = locOtherUndecayedNodesToProcess;
			ResonanceNode* locNewNode = locNewOtherUndecayedNodesToProcess.back();
			locNewOtherUndecayedNodesToProcess.pop_back();
			if(!Loop_PrintingAmplitudes(locNewNode, locNewOtherUndecayedNodesToProcess, locSumName, locHelicityTimes2s, locNewAmplitudeArgumentsVector, locNewBWArgumentsVector, locNewContractionIndices))
				return false; //this amplitude is zero for this helicity sum
		}
		else if(locPsi2->dIsResonanceFlag && (!locPsi1->dIsResonanceFlag))
		{
			if(!Loop_PrintingAmplitudes(locCurrentNode->dRightFinalNode, locOtherUndecayedNodesToProcess, locSumName, locHelicityTimes2s, locNewAmplitudeArgumentsVector, locNewBWArgumentsVector, locNewContractionIndices))
				return false; //this amplitude is zero for this helicity sum
		}
		else if(locPsi1->dIsResonanceFlag && (!locPsi2->dIsResonanceFlag))
		{
			if(!Loop_PrintingAmplitudes(locCurrentNode->dLeftFinalNode, locOtherUndecayedNodesToProcess, locSumName, locHelicityTimes2s, locNewAmplitudeArgumentsVector, locNewBWArgumentsVector, locNewContractionIndices))
				return false; //this amplitude is zero for this helicity sum
		}
		else //both are resonances
		{
			deque<ResonanceNode*> locNewOtherUndecayedNodesToProcess = locOtherUndecayedNodesToProcess;
			locNewOtherUndecayedNodesToProcess.push_back(locCurrentNode->dRightFinalNode);
			if(!Loop_PrintingAmplitudes(locCurrentNode->dLeftFinalNode, locNewOtherUndecayedNodesToProcess, locSumName, locHelicityTimes2s, locNewAmplitudeArgumentsVector, locNewBWArgumentsVector, locNewContractionIndices))
				return false; //this amplitude is zero for this helicity sum
		}
	}

	if(locNumContractionSchemesUsed == 0)
		return false; //all schemes were zero: this node cannot contribute with its current characteristics
	return true;
}

void HCA_ConfigGenerator::Calculate_PossibleHelicityRange(ResonanceNode* locCurrentNode, deque<int>& locHelicityTimes2s, int& locMinH12x2, int& locMaxH12x2)
{
	//locHelicityTimes2s order: beam, target, dFinalStateParticles ()
	Particle* locPsi1 = locCurrentNode->dLeftFinalNode->dNodeParticle;
	Particle* locPsi2 = locCurrentNode->dRightFinalNode->dNodeParticle;

	if(locCurrentNode->dIsProductionNodeFlag) //psi1 is beam, psi2 is target: first two of locHelicityTimes2s
	{
		locMinH12x2 = locHelicityTimes2s[0] - locHelicityTimes2s[1];
		locMaxH12x2 = locMinH12x2;
	}
	else if((!locPsi1->dIsResonanceFlag) && (!locPsi2->dIsResonanceFlag))
	{
		int locParticleIndex1 = dParticleIndexMap[locPsi1] + 1; //+1: target
		int locParticleIndex2 = dParticleIndexMap[locPsi2] + 1; //+1: target
		locMinH12x2 = locHelicityTimes2s[locParticleIndex1] - locHelicityTimes2s[locParticleIndex2];
		locMaxH12x2 = locMinH12x2;
	}
	else if(locPsi1->dIsResonanceFlag && locPsi2->dIsResonanceFlag) //both are resonances
	{
		locMaxH12x2 = locPsi1->dJTimes2 + locPsi2->dJTimes2;
		locMinH12x2 = -1*locMaxH12x2;
	}
	else if((!locPsi1->dIsResonanceFlag) && locPsi2->dIsResonanceFlag) //Psi2 is a resonance, but not Psi1
	{
		int locParticleIndex1 = dParticleIndexMap[locPsi1] + 1; //+1: target
		locMinH12x2 = locHelicityTimes2s[locParticleIndex1] - locPsi2->dJTimes2;
		locMaxH12x2 = locHelicityTimes2s[locParticleIndex1] + locPsi2->dJTimes2;
	}
	else //Psi1 is a resonance, but not Psi2
	{
		int locParticleIndex2 = dParticleIndexMap[locPsi2] + 1; //+1: target
		locMinH12x2 = -1*locPsi1->dJTimes2 - locHelicityTimes2s[locParticleIndex2];
		locMaxH12x2 = locPsi1->dJTimes2 - locHelicityTimes2s[locParticleIndex2];
	}
}

bool HCA_ConfigGenerator::Check_IsNodeAmplitudeZero(ResonanceNode* locCurrentNode, int locH12x2, int locContractionIndex)
{
	//Must check Jx2 S12x2 first: if out of range, would not have calc'd tensor results
	if((abs(locH12x2) > locCurrentNode->dNodeParticle->dJTimes2) || (abs(locH12x2) > locCurrentNode->dS12Times2))
		return true; //0: helicity out of range
	if(locCurrentNode->dTensorAmplitudes[locH12x2][locContractionIndex].dGammaTerms.empty())
		return true; //0: tensor contraction result is zero
	return false;
}