#include "ResonanceNodeGenerator.h"

deque<NodeBase*>& ResonanceNodeGenerator::Generate(int locMaxSChannelJTimes2, int locMaxDecayJTimes2, deque<ExchangeNode*>& locExchangeNodes)
{
	//if locMaxSChannelJTimes2 < 0, no s-channel (generate for other, input channels)
	dMaxSChannelJTimes2 = locMaxSChannelJTimes2;
	int locMinSChannelJTimes2 = ((locMaxSChannelJTimes2 % 2) == 0) ? 0 : 1;
	if(locMaxSChannelJTimes2 < 0)
		locMinSChannelJTimes2 = 5; //don't do loop

	dMaxDecayJTimes2 = locMaxDecayJTimes2;

	int locPsiXCharge = ParticleCharge(dBeamPID) + ParticleCharge(dTargetPID);
	double locPsiXMass = ParticleMass(dBeamPID) + ParticleMass(dTargetPID) + dDefaultBWMassOffset;
	double locPsiXWidth = dDefaultBWWidth;
	Particle_t locPsiXPID = Unknown;

	//generate all possible reaction trees for s-channel resonances
	for(int locPsiXJTimes2 = locMinSChannelJTimes2; locPsiXJTimes2 <= dMaxSChannelJTimes2; locPsiXJTimes2 += 2)
	{
		for(int locPsiXParity = -1; locPsiXParity <= 1; locPsiXParity += 2)
		{
			Particle* locPsiX = new Particle(locPsiXPID, locPsiXCharge, locPsiXJTimes2, locPsiXParity, locPsiXMass, locPsiXWidth, true);

			//create production node
				//don't use multipole basis:
				//amplitudes with different L are independent due to barrier factor (would also be different due to resonance width if not binning in sqrt(s)
			ResonanceNode* locProductionNode = new ResonanceNode(locPsiX, NULL);
			gProductionNode = locProductionNode;
			locProductionNode->dIsProductionNodeFlag = true;

			//create target & beam nodes
			ResonanceNode* locBeamNode = new ResonanceNode(dPsiBeam, locProductionNode);
			ResonanceNode* locTargetNode = new ResonanceNode(dPsiTarget, locProductionNode);
			locProductionNode->dLeftFinalNode = locBeamNode;
			locProductionNode->dRightFinalNode = locTargetNode;

			//create resonance decay node
			ResonanceNode* locParentNode = new ResonanceNode(locPsiX, locProductionNode);
			locProductionNode->dPreviousDecayNode = locParentNode;
			deque<ResonanceNode*> locOtherUndecayedNodesToProcess(1, locParentNode);

			//Pick production L & S, then decay s-channel resonance
			LoopOver_DecayAmplitudes(locProductionNode, dFinalStateParticleSet, locOtherUndecayedNodesToProcess);

			//if tree was good have saved a clone of this node, but not these nodes: delete so no memory leak
			delete locParentNode;
			delete locProductionNode;
			delete locBeamNode;
			delete locTargetNode;
		}
	}

	//GENERATE all possible reaction trees for all input exchange-nodes: including all possible J/P/S/L
	for(size_t loc_i = 0; loc_i < locExchangeNodes.size(); ++loc_i)
	{
		ExchangeNode* locExchangeNode = locExchangeNodes[loc_i];
		locExchangeNode->dBeamParticle = dPsiBeam;
		locExchangeNode->dTargetParticle = dPsiTarget;
		Generate_ExchangeTrees(locExchangeNode);
	}

	if(dDebugFlag)
		cout << dSavedTreeNodes.size()  << " decay trees generated" << endl;

	if(dDebugFlag)
	{
		for(size_t loc_i = 0; loc_i < dSavedTreeNodes.size(); ++loc_i)
		{
			cout << "DecayTree " << loc_i + 1 << ":" << endl;
			Print_Tree(dSavedTreeNodes[loc_i]);
			cout << endl;
		}
	}

	//loop through and set resonance masses
	for(size_t loc_i = 0; loc_i < dSavedTreeNodes.size(); ++loc_i)
		CalcAndSet_ResonanceMasses_Init(dSavedTreeNodes[loc_i]);

	return dSavedTreeNodes;
}

void ResonanceNodeGenerator::Generate_ExchangeTrees(ExchangeNode* locExchangeNode)
{
	gProductionNode = locExchangeNode;

	set<Particle*> locUnusedFinalStateParticles = dFinalStateParticleSet;
	deque<ResonanceNode*> locUndecayedNodesToProcess;

	//Build Left Node
	Particle_t locLeftPID = locExchangeNode->dLeftFinalPID;
	map<Particle_t, deque<Particle*> >::iterator locIterator = dFinalStateParticleMap.find(locLeftPID);
	if(locIterator != dFinalStateParticleMap.end())
	{
		//Final State Particle
		Particle* locParticle = (locIterator->second)[0];
		locUnusedFinalStateParticles.erase(locParticle);

		ResonanceNode* locParticleNode = new ResonanceNode(locParticle, locExchangeNode);
		locExchangeNode->dLeftFinalNode = locParticleNode;
	}
	else if(locLeftPID != Unknown)
	{
		//Resonance with PID given: use specific J^P
		Particle* locPsiX = new Particle(locLeftPID);
		ResonanceNode* locParticleNode = new ResonanceNode(locPsiX, locExchangeNode);
		locExchangeNode->dLeftFinalNode = locParticleNode;
		locUndecayedNodesToProcess.push_back(locParticleNode);
	}

	//Build Right Node
	Particle_t locRightPID = locExchangeNode->dRightFinalPID;
	locIterator = dFinalStateParticleMap.find(locRightPID);
	if(locIterator != dFinalStateParticleMap.end())
	{
		//Final State Particle
		size_t locIndex = (locRightPID == locLeftPID) ? 1 : 0;
		Particle* locParticle = (locIterator->second)[locIndex];
		locUnusedFinalStateParticles.erase(locParticle);

		ResonanceNode* locParticleNode = new ResonanceNode(locParticle, locExchangeNode);
		locExchangeNode->dRightFinalNode = locParticleNode;
	}
	else if(locRightPID != Unknown)
	{
		//Resonance with PID given: use specific J^P
		Particle* locPsiX = new Particle(locRightPID);
		ResonanceNode* locParticleNode = new ResonanceNode(locPsiX, locExchangeNode);
		locExchangeNode->dRightFinalNode = locParticleNode;
		locUndecayedNodesToProcess.push_back(locParticleNode);
	}

	//Decay Resonances
	if((locLeftPID != Unknown) && (locRightPID != Unknown))
	{
		ResonanceNode* locResonanceNode = locUndecayedNodesToProcess.back();
		locUndecayedNodesToProcess.pop_back();
		Generate_DecayTrees(locResonanceNode, locUnusedFinalStateParticles, locUndecayedNodesToProcess);
	}
	else if((locLeftPID == Unknown) && (locRightPID == Unknown)) //both unknown: you poor fool
	{
		//pick charge
		for(int locPsiXCharge = dMinParticleCharge; locPsiXCharge <= dMaxParticleCharge; ++locPsiXCharge)
		{
			int locPsiYCharge = ParticleCharge(dBeamPID) + ParticleCharge(dTargetPID) - locPsiXCharge;
			if((locPsiYCharge < dMinParticleCharge) || (locPsiYCharge > dMaxParticleCharge))
				continue; //bad charge

			//left particle (X)
			double locPsiXMass = ParticleMass(dBeamPID) + ParticleMass(dTargetPID) + dDefaultBWMassOffset;
			for(int locPsiXJTimes2 = 0; locPsiXJTimes2 <= dMaxDecayJTimes2; ++locPsiXJTimes2)
			{
				for(int locPsiXParity = -1; locPsiXParity <= 1; locPsiXParity += 2)
				{
					Particle* locPsiX = new Particle(Unknown, locPsiXCharge, locPsiXJTimes2, locPsiXParity, locPsiXMass, dDefaultBWWidth, true);
					ResonanceNode* locParticleNode = new ResonanceNode(locPsiX, locExchangeNode);
					locExchangeNode->dLeftFinalNode = locParticleNode;
					locUndecayedNodesToProcess.push_back(locParticleNode);

					//right particle (Y)
					for(int locPsiYJTimes2 = 0; locPsiYJTimes2 <= dMaxDecayJTimes2; ++locPsiYJTimes2)
					{
						for(int locPsiYParity = -1; locPsiYParity <= 1; locPsiYParity += 2)
						{
							Particle* locPsiY = new Particle(Unknown, locPsiYCharge, locPsiYJTimes2, locPsiYParity, locPsiXMass, dDefaultBWWidth, true);
							ResonanceNode* locParticleNode = new ResonanceNode(locPsiY, locExchangeNode);
							locExchangeNode->dRightFinalNode = locParticleNode;
							Generate_DecayTrees(locParticleNode, locUnusedFinalStateParticles, locUndecayedNodesToProcess);
						}
					}
				}
			}
		}
	}
	else //one is unknown, the other isn't
	{
		Particle_t locKnownPID = (locLeftPID == Unknown) ? locRightPID : locLeftPID;

		int locPsiXCharge = ParticleCharge(dBeamPID) + ParticleCharge(dTargetPID) - ParticleCharge(locKnownPID);
		double locPsiXMass = ParticleMass(dBeamPID) + ParticleMass(dTargetPID) - ParticleMass(locKnownPID) + dDefaultBWMassOffset;
		for(int locPsiXJTimes2 = 0; locPsiXJTimes2 <= dMaxDecayJTimes2; ++locPsiXJTimes2)
		{
			for(int locPsiXParity = -1; locPsiXParity <= 1; locPsiXParity += 2)
			{
				Particle* locPsiX = new Particle(Unknown, locPsiXCharge, locPsiXJTimes2, locPsiXParity, locPsiXMass, dDefaultBWWidth, true);
				ResonanceNode* locParticleNode = new ResonanceNode(locPsiX, locExchangeNode);
				if(locRightPID == Unknown)
					locExchangeNode->dRightFinalNode = locParticleNode;
				else
					locExchangeNode->dLeftFinalNode = locParticleNode;
				Generate_DecayTrees(locParticleNode, locUnusedFinalStateParticles, locUndecayedNodesToProcess);
			}
		}
	}
}

void ResonanceNodeGenerator::Generate_DecayTrees(ResonanceNode* locCurrentNode, set<Particle*> locRemainingFinalStateParticles, deque<ResonanceNode*> locOtherUndecayedNodesToProcess)
{
	Particle* locPsiX = locCurrentNode->dNodeParticle;

	//recursive function to generate decay trees
	//each decay node must have 0 (final state particle) or 2 (decaying particle) decay products
	//each decaying particle can have anywhere between 0 to 2 resonances as decay products (except when limited by final-state particle availability)
		//loop through each possible value of num resonances: loop through all possible tree nodes: will generate all possible trees

	//min num final state particles associated with OTHER (not THIS) undecayed resonances:
	int locMinPossibleFinalStateParticles_OtherDecays = 2*locOtherUndecayedNodesToProcess.size();

	//max num final state particles that could be (eventual) decay products from the decay chain starting with THIS: 
	int locMaxPossibleFinalStateParticles_ThisDecayChain = locRemainingFinalStateParticles.size() - locMinPossibleFinalStateParticles_OtherDecays;

	if(locMaxPossibleFinalStateParticles_ThisDecayChain < 2) //shouldn't be possible
		abort(); //not enough decay products for this decaying resonance: something seriously wrong happened

	//find range of # possible product resonances
	int locMaxNumProductResonances = 2;
	int locMinNumProductResonances = 0;
	if(locMaxPossibleFinalStateParticles_ThisDecayChain == 2)
		locMaxNumProductResonances = 0; // cannot have more resonances
	else if(locMaxPossibleFinalStateParticles_ThisDecayChain == 3)
		locMaxNumProductResonances = 1; //this decay cannot go to 2 resonances: would have 4 products
	if(locOtherUndecayedNodesToProcess.empty() && (locMaxPossibleFinalStateParticles_ThisDecayChain >= 3))
		locMinNumProductResonances = 1; //this is the last resonance, and must have 3+ products: need at least one resonance
	//never need locMinNumProductResonances = 2

	//loop over # resonance particles
	for(int locNumResonanceProducts = locMinNumProductResonances; locNumResonanceProducts <= locMaxNumProductResonances; ++locNumResonanceProducts)
	{	
		//loop over the products that are final state particles
		int locNumFinalStateProducts = 2 - locNumResonanceProducts; //2-body decay for mother fit
		if(locNumFinalStateProducts > 0)
		{
			//loop over all possible particles for the first product
			set<Particle*>::iterator locIterator = locRemainingFinalStateParticles.begin();
			for(; locIterator != locRemainingFinalStateParticles.end(); ++locIterator)
			{
				if(dFixParticlePermutationFlag && (locIterator != locRemainingFinalStateParticles.begin()))
					break;
				Particle* locParticle1 = *locIterator;
				int locParticle2Charge = locPsiX->dCharge - locParticle1->dCharge;
				if((locParticle2Charge < dMinParticleCharge) || (locParticle2Charge > dMaxParticleCharge))
					continue; //particle 2 will have a bad charge (e.g. --)

				ResonanceNode* locParticle1Node = new ResonanceNode(locParticle1, locCurrentNode);
				locCurrentNode->dLeftFinalNode = locParticle1Node;

				if(locNumFinalStateProducts == 2)
				{
					//loop over all possible particles for the second product
					set<Particle*>::iterator locIterator2 = locIterator;
					++locIterator2;
					set<Particle*>::iterator locStartIterator = locIterator2;
					for(; locIterator2 != locRemainingFinalStateParticles.end(); ++locIterator2)
					{
						if(dFixParticlePermutationFlag && (locIterator2 != locStartIterator))
							break;
						Particle* locParticle2 = *locIterator2;
						//check if charge is conserved
						if(locParticle2->dCharge != locParticle2Charge)
							continue;
						ResonanceNode* locParticle2Node = new ResonanceNode(locParticle2, locCurrentNode);
						locCurrentNode->dRightFinalNode = locParticle2Node;
						//particles selected; now select L12, S12 and save:
						LoopOver_DecayAmplitudes(locCurrentNode, locRemainingFinalStateParticles, locOtherUndecayedNodesToProcess);
						delete locParticle2Node; //if tree was good have saved a clone of this node, but not this node: delete so no memory leak
					}
				}
				else //product #2 is a resonance //loop over its properties, select L12 & S12, & decay it
					LoopOver_DecayAmplitudes_PickPsi2(locCurrentNode, locRemainingFinalStateParticles, locOtherUndecayedNodesToProcess);

				delete locParticle1Node; //if tree was good have saved a clone of this node, but not this node: delete so no memory leak
			}
		}
		else //both particles are resonances //loop over psi1 properties first, the psi2, then select L12 & S12, then decay one of them and save the other to decay later
			LoopOver_DecayAmplitudes_PickPsi1(locCurrentNode, locRemainingFinalStateParticles, locOtherUndecayedNodesToProcess);
	}
}

void ResonanceNodeGenerator::LoopOver_DecayAmplitudes_PickPsi1(ResonanceNode* locCurrentNode, set<Particle*> locRemainingFinalStateParticles, deque<ResonanceNode*> locOtherUndecayedNodesToProcess)
{
	Particle* locPsiX = locCurrentNode->dNodeParticle;

	Particle_t locPsi1PID = Unknown;
	double locPsi1Mass = -1.0;
	double locPsi1Width = dDefaultBWWidth;

	//Loop over possible charge
	for(int locPsi1Charge = dMinParticleCharge; locPsi1Charge <= dMaxParticleCharge; ++locPsi1Charge)
	{
		int locPsi2Charge = locPsiX->dCharge - locPsi1Charge;
		if((locPsi2Charge < dMinParticleCharge) || (locPsi2Charge > dMaxParticleCharge))
			continue; //psi2 would be invalid charge

		//Loop over possible P_{1} (-1 or 1)
		for(int locPsi1Parity = -1; locPsi1Parity <= 1; locPsi1Parity += 2)
		{
			//Loop over possible S_{1} (from 0 to max)
				//dunno if half-integer spin or not: loop over all possibilities
			for(int locPsi1JTimes2 = 0; locPsi1JTimes2 <= dMaxDecayJTimes2; locPsi1JTimes2 += 1) //??
			{
				Particle* locPsi1 = new Particle(locPsi1PID, locPsi1Charge, locPsi1JTimes2, locPsi1Parity, locPsi1Mass, locPsi1Width, true);
				ResonanceNode* locParticle1Node = new ResonanceNode(locPsi1, locCurrentNode);
				locCurrentNode->dLeftFinalNode = locParticle1Node;
				locOtherUndecayedNodesToProcess.push_back(locParticle1Node); //save to do the to-do list, will decay particle 2 once made
				LoopOver_DecayAmplitudes_PickPsi2(locCurrentNode, locRemainingFinalStateParticles, locOtherUndecayedNodesToProcess);
				delete locParticle1Node; //if tree was good have saved a clone of this node, but not this node: delete so no memory leak
			}
		}
	}
}

void ResonanceNodeGenerator::LoopOver_DecayAmplitudes_PickPsi2(ResonanceNode* locCurrentNode, set<Particle*> locRemainingFinalStateParticles, deque<ResonanceNode*> locOtherUndecayedNodesToProcess)
{
	Particle* locPsiX = locCurrentNode->dNodeParticle;
	Particle* locPsi1 = locCurrentNode->dLeftFinalNode->dNodeParticle;

	Particle_t locPsi2PID = Unknown;
	int locPsi2Charge = locPsiX->dCharge - locPsi1->dCharge;
	if((locPsi2Charge < dMinParticleCharge) || (locPsi2Charge > dMaxParticleCharge))
		return; //psi2 would be invalid charge
	double locPsi2Mass = -1.0;
	double locPsi2Width = dDefaultBWWidth;

	//Loop over possible P_{2} (-1 or 1)
	for(int locPsi2Parity = -1; locPsi2Parity <= 1; locPsi2Parity += 2)
	{
		//Loop over possible S_{2} (from 0 to max)
		bool locIsJ2IntegerSpin = ((locPsiX->dJTimes2 + locPsi1->dJTimes2) % 2 == 0); //S_{2} = J_{R} - S_{1} - L_{12}
		int locMinPsi2JTimes2 = (locIsJ2IntegerSpin ? 0 : 1);
		for(int locPsi2JTimes2 = locMinPsi2JTimes2; locPsi2JTimes2 <= dMaxDecayJTimes2; locPsi2JTimes2 += 2)
		{
			Particle* locPsi2 = new Particle(locPsi2PID, locPsi2Charge, locPsi2JTimes2, locPsi2Parity, locPsi2Mass, locPsi2Width, true);
			ResonanceNode* locParticle2Node = new ResonanceNode(locPsi2, locCurrentNode);
			locCurrentNode->dRightFinalNode = locParticle2Node;
			LoopOver_DecayAmplitudes(locCurrentNode, locRemainingFinalStateParticles, locOtherUndecayedNodesToProcess);
			delete locParticle2Node; //if tree was good have saved a clone of this node, but not this node: delete so no memory leak
		}
	}
}

void ResonanceNodeGenerator::LoopOver_DecayAmplitudes(ResonanceNode* locCurrentNode, set<Particle*> locRemainingFinalStateParticles, deque<ResonanceNode*> locOtherUndecayedNodesToProcess)
{
	Particle* locPsiX = locCurrentNode->dNodeParticle;
	Particle* locPsi1 = locCurrentNode->dLeftFinalNode->dNodeParticle;
	Particle* locPsi2 = locCurrentNode->dRightFinalNode->dNodeParticle;

	//Determine if legal J's input
	if(((locPsiX->dJTimes2 + locPsi1->dJTimes2 + locPsi2->dJTimes2) % 2) == 1)
		return; //not possible: J_{R} = S_{1} + S_{2} + L_{12}, so 2*L_{12} (and thus the sum above) must be even

	//Determine whether L_{12} is even or odd
	//P_{R} = P_{1}*P_{2}*(-1)^(L_{12})
	bool locL12IsEvenFlag = ((locPsi1->dParity*locPsi2->dParity) == locPsiX->dParity);

	//Loop over possible S_{12} (from S_{1} & S_{2})
	int locMinS12Times2 = abs(locPsi1->dJTimes2 - locPsi2->dJTimes2);
	int locMaxS12Times2 = locPsi1->dJTimes2 + locPsi2->dJTimes2;
	set<int> locChosenL12s;
	for(int locS12Times2 = locMinS12Times2; locS12Times2 <= locMaxS12Times2; locS12Times2 += 2)
	{
		//Loop over possible L_{12} (from J_{R} & S_{12})
		//Min L_{12}
		int locMinL12 = abs(locPsiX->dJTimes2 - locS12Times2)/2;
		if(locL12IsEvenFlag != ((locMinL12 % 2) == 0))
			++locMinL12;
		//Max L_{12}
		int locMaxL12 = (locPsiX->dJTimes2 + locS12Times2)/2;
		if((locMaxL12 + locMinL12) % 2 == 1)
			--locMaxL12;

		for(int locL12 = locMinL12; locL12 <= locMaxL12; locL12 += 2) //increments of 2: narrowed to even/odd
		{
			if(dLimitToOneSFlag && (locChosenL12s.find(locL12) != locChosenL12s.end()))
				break; //This L12 already chosen for a previous S12
			locChosenL12s.insert(locL12);

			//J/L/S chosen. 
			locCurrentNode->dL12 = locL12;
			locCurrentNode->dS12Times2 = locS12Times2;

			//Handle decaying resonances
			if((!locPsi1->dIsResonanceFlag) && (!locPsi2->dIsResonanceFlag))
			{
				if(locOtherUndecayedNodesToProcess.empty()) //this tree has finished construction: save a clone of it (so don't overwrite it
				{
					ResonanceNode* locResonanceNode = dynamic_cast<ResonanceNode*>(gProductionNode);
					ExchangeNode* locExchangeNode = dynamic_cast<ExchangeNode*>(gProductionNode);
					if(locResonanceNode != NULL)
						dSavedTreeNodes.push_back(Clone_DecayTree(locResonanceNode));
					else if(locExchangeNode != NULL)
					{
						//Register any fit params that were present
						for(size_t loc_i = 0; loc_i < locExchangeNode->dFitParams.size(); ++loc_i)
							dFitParams.insert(locExchangeNode->dFitParams[loc_i]);
						dSavedTreeNodes.push_back(Clone_DecayTree(locExchangeNode));
					}
					else
					{
						cout << "ABORTING: NODE TYPE NOT RECOGNIZED." << endl;
						abort();
					}
				}
				else
				{
					set<Particle*> locNewRemainingFinalStateParticles = locRemainingFinalStateParticles;
					if(!locCurrentNode->dIsProductionNodeFlag)
					{
						locNewRemainingFinalStateParticles.erase(locPsi1);
						locNewRemainingFinalStateParticles.erase(locPsi2);
					}

					//this tree branch has finished construction, but not the whole tree: hop over to the next node
					deque<ResonanceNode*> locNewNodesToProcess = locOtherUndecayedNodesToProcess;
					ResonanceNode* locNextNodeToProcess = locNewNodesToProcess.back();
					locNewNodesToProcess.pop_back();
					Generate_DecayTrees(locNextNodeToProcess, locNewRemainingFinalStateParticles, locNewNodesToProcess);
				}
				continue;
			}

			//decay particle 2 (if particle 1 is resonance, will be decayed later):
			set<Particle*> locNewRemainingFinalStateParticles = locRemainingFinalStateParticles;
			if(!locPsi1->dIsResonanceFlag)
				locNewRemainingFinalStateParticles.erase(locPsi1);
			Generate_DecayTrees(locCurrentNode->dRightFinalNode, locNewRemainingFinalStateParticles, locOtherUndecayedNodesToProcess);
		}
	}
}

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

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

	ostringstream 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 << "), ";

	cout << locCurrentLine.str();
	if(!locPrintChainFlag)
		return;

	if(locPsi1->dIsResonanceFlag)
		Print_ResonanceNode(locExchangeNode->dLeftFinalNode, locPrintChainFlag, locNumTabs + 1);
	if(locPsi2->dIsResonanceFlag)
		Print_ResonanceNode(locExchangeNode->dRightFinalNode, locPrintChainFlag, locNumTabs + 1);
}

void ResonanceNodeGenerator::Print_ResonanceNode(ResonanceNode* locCurrentNode, bool locPrintChainFlag, int locNumTabs)
{
	ostringstream 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;
		cout << locCurrentLine.str();
		if(locPrintChainFlag)
			Print_ResonanceNode((ResonanceNode*)locCurrentNode->dPreviousDecayNode, locPrintChainFlag, 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;

	cout << locCurrentLine.str();
	if(!locPrintChainFlag)
		return;

	if(locPsi1->dIsResonanceFlag)
		Print_ResonanceNode(locCurrentNode->dLeftFinalNode, locPrintChainFlag, locNumTabs + 1);
	if(locPsi2->dIsResonanceFlag)
		Print_ResonanceNode(locCurrentNode->dRightFinalNode, locPrintChainFlag, locNumTabs + 1);
}

double ResonanceNodeGenerator::CalcAndSet_ResonanceMasses_Init(NodeBase* locInputTransitionTree)
{
	//Initialize (Get starting node(s))
	ResonanceNode* locResonanceNode = dynamic_cast<ResonanceNode*>(locInputTransitionTree);
	if(locResonanceNode != NULL)
		return CalcAndSet_ResonanceMasses(locResonanceNode);

	//is exchange node
	ExchangeNode* locExchangeNode = dynamic_cast<ExchangeNode*>(locInputTransitionTree);
	Particle* locPsi1 = locExchangeNode->dLeftFinalNode->dNodeParticle;
	if(locPsi1->dIsResonanceFlag)
		CalcAndSet_ResonanceMasses(locExchangeNode->dLeftFinalNode);

	Particle* locPsi2 = locExchangeNode->dRightFinalNode->dNodeParticle;
	if(locPsi2->dIsResonanceFlag)
		CalcAndSet_ResonanceMasses(locExchangeNode->dRightFinalNode);

	return 0.0;
}

double ResonanceNodeGenerator::CalcAndSet_ResonanceMasses(ResonanceNode* locCurrentNode)
{
	if(locCurrentNode->dIsProductionNodeFlag)
		return CalcAndSet_ResonanceMasses((ResonanceNode*)locCurrentNode->dPreviousDecayNode);

	//assign threshold mass + 0.05, unless decay recognized
	//also set width param name
	Particle* locPsiX = locCurrentNode->dNodeParticle;
	Particle* locPsi1 = locCurrentNode->dLeftFinalNode->dNodeParticle;
	Particle* locPsi2 = locCurrentNode->dRightFinalNode->dNodeParticle;

	double locResonanceMass = dDefaultBWMassOffset;
	if(locPsi1->dIsResonanceFlag)
		locResonanceMass += CalcAndSet_ResonanceMasses(locCurrentNode->dLeftFinalNode);
	else
		locResonanceMass += locPsi1->dMass;

	if(locPsi2->dIsResonanceFlag)
		locResonanceMass += CalcAndSet_ResonanceMasses(locCurrentNode->dRightFinalNode);
	else
		locResonanceMass += locPsi2->dMass;

	locPsiX->dMass = locResonanceMass;
	locPsiX->dMassParamName = "";
	locPsiX->dWidthParamName = "";

	//Loop through input resonance decays, matching to this decay
	locCurrentNode->dResonanceRecognizedFlag = false;
	for(size_t loc_i = 0; loc_i < dResonanceDecays.size(); ++loc_i)
	{
		if(!dResonanceDecays[loc_i]->IsMatch(locPsiX, locPsi1, locPsi2))
			continue;

		locPsiX->dMass = dResonanceDecays[loc_i]->dPsiXMass;
		locPsiX->dWidth = dResonanceDecays[loc_i]->dPsiXWidth;

		if(dResonanceDecays[loc_i]->dMassFitParamInfo != NULL)
		{
			locPsiX->dMassParamName = dResonanceDecays[loc_i]->dMassFitParamInfo->dName;
			dFitParams.insert(dResonanceDecays[loc_i]->dMassFitParamInfo);
		}
		if(dResonanceDecays[loc_i]->dWidthFitParamInfo != NULL)
		{
			locPsiX->dWidthParamName = dResonanceDecays[loc_i]->dWidthFitParamInfo->dName;
			dFitParams.insert(dResonanceDecays[loc_i]->dWidthFitParamInfo);
		}

		locCurrentNode->dResonanceRecognizedFlag = true;
		break;
	}

	return locPsiX->dMass;
}