#ifndef _HCAmplitude_Chain_
#define _HCAmplitude_Chain_

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

#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Vector/LorentzRotation.h"
#include "CLHEP/Vector/ThreeVector.h"

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

#include "particleType.h"
#include "resonanceinfo.h"

#include "HCAmplitude_Step.h"
#include "TChannel_Basic.h"

using namespace std;

class Kinematics;

class HCAmplitude_Chain : public UserAmplitude<HCAmplitude_Chain>
{
	public:
		HCAmplitude_Chain() : UserAmplitude<HCAmplitude_Chain>(){}
		HCAmplitude_Chain(const vector<string>& locArguments);

		string name() const{return "HCAmplitude_Chain";}

		complex<GDouble> calcAmplitude(GDouble** locKinematics) const;

	private:

		void ExtractInts_SingleUnderscore(string locInfoString, vector<int>& locInfoVector);
		void ExtractInts_DoubleUnderscore(string locInfoString, vector<vector<int> >& locInfoVector);
		int BuildArray(const vector<int>& locInputVector, int*& locArray) const;

		void Build_DecayAmplitudes(int locStepIndex, int locKidIndex, vector<int> locCurrentHelicitiesTimes2, vector<int> locCurrentJsTimes2, vector<HCAmplitude_Step*>& locStepAmplitudes);
		void Create_DecayAmplitude(int locStepIndex, vector<int> locCurrentHelicitiesTimes2, vector<int> locCurrentJsTimes2, vector<HCAmplitude_Step*>& locStepAmplitudes);
		bool Handle_HelicityDecursion(int& locStepIndex, deque<int>& locResonanceHelicityIndices, deque<int>& locResumeAtValues);
		void Find_AmplitudeIndices(const deque<int>& locStepHelicityIndices);

		//For creating amplitude objects
		vector<int> dExternalHelicitiesTimes2;
		vector<Particle_t> dPIDs;

		vector<int> dJXTimes2Vector;
		vector<int> dL12Vector;
		vector<int> dS12Times2Vector;
		vector<int> dContractionIndexVector;
		vector<vector<vector<int> > > dChainParticleIndices;

		deque<deque<int> > dChainAmplitudeIndices;

		//for HC-amplitude calculation
		int dNumSteps;
		int dNumAmplitudeCombos;
		int** dAmplitudeHelicityComboIndices; //2D: combo (size = dNumAmplitudeCombos), step (size = dNumSteps)
		int* dAmplitudeArrayPermutationSizes; //sizes of 2nd dimension of dDecayAmplitudes 
		HCAmplitude_Step*** dDecayAmplitudes; //2D array of pointers: step (size = dNumSteps), permutation
		complex<GDouble>** dAmplitudeValues; //results of dDecayAmplitudes computations (has same size)

		//for other amplitude component calculations
		Amplitude* dOtherAmplitudeComponent;
};

inline void HCAmplitude_Chain::ExtractInts_SingleUnderscore(string locInfoString, vector<int>& locInfoVector)
{
	locInfoVector.clear();
	while(true)
	{
		size_t locUnderscoreIndex = locInfoString.find("_");
		string locSubString = locInfoString.substr(0, locUnderscoreIndex);
		istringstream locTempStream(locSubString);
		int locInfoInt = -100;
		locTempStream >> locInfoInt;
		locInfoVector.push_back(locInfoInt);
		if(locUnderscoreIndex == string::npos)
			return;
		locInfoString = locInfoString.substr(locUnderscoreIndex + 1);
	}
}

inline void HCAmplitude_Chain::ExtractInts_DoubleUnderscore(string locInfoString, vector<vector<int> >& locTotalInfoVector)
{
	locTotalInfoVector.clear();
	while(true)
	{
		size_t locDoubleUnderscoreIndex = locInfoString.find("__");
		string locSubString = locInfoString.substr(0, locDoubleUnderscoreIndex);
		vector<int> locInfoVector;
		ExtractInts_SingleUnderscore(locSubString, locInfoVector);
		locTotalInfoVector.push_back(locInfoVector);
		if(locDoubleUnderscoreIndex == string::npos)
			return;
		locInfoString = locInfoString.substr(locDoubleUnderscoreIndex + 2);
	}
}

inline int HCAmplitude_Chain::BuildArray(const vector<int>& locInputVector, int*& locArray) const
{
	locArray = new int[locInputVector.size()];
	for(size_t loc_i = 0; loc_i < locInputVector.size(); ++loc_i)
		locArray[loc_i] = locInputVector[loc_i];
	return locInputVector.size();
}

#endif