Merge pull request #3 from gwm17/devel

Overhaul input and reaction system construction. Much more streamline…

input.txt

@@ -1,12 +1,38 @@
 ----------Data Information----------
 OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
 ----------Reaction Information----------
-begin_nuclei(Z,A)
+NumberOfSamples: 100000
-	5	10
+begin_chain
-	2	3
+	begin_step
-	2	4
+		Type: Reaction
-	2	4
+		begin_nuclei
-end_nuclei(Z,A)
+			5 10
+			2 3
+			2 4
+		end_nuclei
+		BeamEnergyMean(MeV): 24.0
+		BeamEnergySigma(MeV): 0.0
+		ThetaType: Lab
+		ThetaMin(deg): 15.0
+		ThetaMax(deg): 15.0
+		PhiMin(deg): 0.0
+		PhMax(deg): 0.0
+		ResidualExcitationMean(MeV): 16.8
+		ResidualExcitationSigma(MeV): 0.023
+	end_step
+	begin_step
+		Type: Decay
+		begin_nuclei
+			5 9
+			2 4
+		end_nuclei
+		PhiMin(deg): 0.0
+		PhMax(deg): 360.0
+		ResidualExcitationMean(MeV): 0.0
+		ResidualExcitationSigma(MeV): 0.0
+		AngularDistributionFile: ./etc/isotropic_dist.txt
+	end_step
+end_chain
 ----------Target Information----------
 NumberOfLayers: 1
 begin_layer
@@ -15,13 +41,4 @@ begin_layer
 		element	5 10 1
 	end_elements
 end_layer
-----------Sampling Information----------
-NumberOfSamples: 1000000
-BeamMeanEnergy(MeV): 24.0 BeamEnergySigma(MeV): 0.0
-EjectileThetaType(0=Lab,1=CM): 0
-EjectileThetaMin(deg): 15.0 EjectileThetaMax(deg): 15.0
-EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 0.0
-ResidualExMean(MeV): 16.8 ResidualExSigma(MeV): 0.023
-Decay1_DistributionFile: ./etc/isotropic_dist.txt
-Decay2_DistributionFile: ./etc/isotropic_dist.txt
 --------------------------------------

src/Mask/AngularDistribution.cpp

@@ -13,7 +13,7 @@ namespace Mask {
 	}
 	
 	AngularDistribution::AngularDistribution(const std::string& file) :
-		m_branchingRatio(1.0), m_L(0), m_isIsotropic(true)
+		m_uniformCosineDist(-1.0, 1.0), m_branchingRatio(1.0), m_L(0), m_isIsotropic(true)
 	{
 		ReadDistributionFile(file);
 	}

src/Mask/DecaySystem.cpp

@@ -5,49 +5,52 @@
 
 namespace Mask {
 
-	DecaySystem::DecaySystem() :
+	DecaySystem::DecaySystem(const std::vector<StepParameters>& params) :
 		ReactionSystem()
 	{
 		m_nuclei.resize(3);
-	}
+		Init(params);
-	
-	DecaySystem::DecaySystem(const std::vector<int>& z, const std::vector<int>& a) :
-		ReactionSystem()
-	{
-		m_nuclei.resize(3);
-		SetNuclei(z, a);
 	}
 	
 	DecaySystem::~DecaySystem() {}
 	
-	bool DecaySystem::SetNuclei(const std::vector<int>& z, const std::vector<int>& a)
+	void DecaySystem::Init(const std::vector<StepParameters>& params)
 	{
-		if(z.size() != a.size() || z.size() != 2)
+		if(params.size() != 1 || params[0].rxnType != RxnType::Decay ||
-			return false;
+		   params[0].Z.size() != 2 || params[0].A.size() != 2)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at DecaySystem::Init(), does not match Decay signature!" << std::endl;
+			return;
+		}
 
-		int zr = z[0] - z[1];
+		const StepParameters& step1Params = params[0];
-		int ar = a[0] - a[1];
 
-		m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
+		int zr = step1Params.Z[0] - step1Params.Z[1];
-		m_nuclei[1] = CreateNucleus(z[1], a[1]); //breakup1
+		int ar = step1Params.A[0] - step1Params.A[1];
+
+		m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
+		m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //breakup1
 		m_nuclei[2] = CreateNucleus(zr, ar); //breakup2
 
 		m_step1.BindNuclei(&(m_nuclei[0]), nullptr, &(m_nuclei[1]), &(m_nuclei[2]));
 		SetSystemEquation();
-		return true;
+
+		
+		AddPhiRange(step1Params.phiMin, step1Params.phiMax);
+		AddDecayAngularDistribution(step1Params.angularDistFile);
+		AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
+
+		return;
 	}
 	
-	std::vector<Nucleus>* DecaySystem::GetNuclei()
+	void DecaySystem::SetLayeredTarget(const LayeredTarget& target)
 	{
-		return &m_nuclei;
+		m_target = target;
-	}
-	
-	void DecaySystem::LinkTarget() {
-		m_step1.SetLayeredTarget(&m_target);
-	
 		m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
 		if(m_rxnLayer != m_target.GetNumberOfLayers())
 		{
+			m_step1.SetLayeredTarget(&m_target);
 			m_step1.SetRxnLayer(m_rxnLayer);
 			m_isTargetSet = true;
 		}
@@ -66,14 +69,14 @@ namespace Mask {
 	
 	void DecaySystem::RunSystem()
 	{
-		//Link up the target if it hasn't been done yet
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
-		if(!m_isTargetSet)
+		double rxnTheta = m_decayAngularDistributions[0].GetRandomCosTheta();
-			LinkTarget();
+		double rxnPhi = m_phiRanges[0](gen);
+		double ex = m_exDistributions[0](gen);
 
-		double rxnTheta = std::acos(m_step1Distribution.GetRandomCosTheta());
-		double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
 		m_step1.SetPolarRxnAngle(rxnTheta);
 		m_step1.SetAzimRxnAngle(rxnPhi);
+		m_step1.SetExcitation(ex);
 		m_step1.SetResidualEnergyLoss(true);
 		m_step1.Calculate();
 	}

src/Mask/DecaySystem.h

@@ -6,29 +6,20 @@
 
 namespace Mask {
 
-	class DecaySystem: public ReactionSystem {
+	class DecaySystem: public ReactionSystem
+	{
 	public:
-		DecaySystem();
+		DecaySystem(const std::vector<StepParameters>& params);
-		DecaySystem(const std::vector<int>& z, const std::vector<int>& a);
 		~DecaySystem();
 	
-		bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
+		virtual void SetLayeredTarget(const LayeredTarget& target) override;
-		void RunSystem() override;
+		virtual void RunSystem() override;
-		std::vector<Nucleus>*GetNuclei() override;
-	
-		virtual void SetDecay1Distribution(const std::string& filename) override { m_step1Distribution.ReadDistributionFile(filename); }
-	
-		int GetDecay1AngularMomentum() { return m_step1Distribution.GetL(); }
-		double GetDecay1BranchingRatio() { return m_step1Distribution.GetBranchingRatio(); }
 	
 	private:
-		void LinkTarget() override;
+		void Init(const std::vector<StepParameters>& params);
 		void SetSystemEquation() override;
 	
 		Reaction m_step1;
-	
-		AngularDistribution m_step1Distribution;
-		
 	};
 
 }

src/Mask/MaskApp.cpp

@@ -10,7 +10,7 @@ namespace Mask {
 	MaskApp::MaskApp() :
 		m_system(nullptr)
 	{
-		std::cout<<"----------GWM Kinematics Simulation----------"<<std::endl;
+		std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
 	}
 	
 	MaskApp::~MaskApp() 
@@ -29,40 +29,87 @@ namespace Mask {
 		}
 	
 		std::string junk;
-		getline(input, junk);
+		std::getline(input, junk);
 		input>>junk>>m_outputName;
+		std::getline(input, junk);
+		std::getline(input, junk);
+		input>>junk>>m_nsamples;
 
-		std::vector<int> avec, zvec, svec;
+		std::vector<StepParameters> params;
-		int z, a, s;
+		int z, a;
-		getline(input, junk);
-		getline(input, junk);
-
 		while(input>>junk)
 		{
-			if(junk == "begin_nuclei(Z,A)")
+			if(junk == "begin_chain")
 				continue;
-			else if (junk == "end_nuclei(Z,A)")
+			else if (junk == "end_chain")
 				break;
-			else
+			else if(junk == "begin_step")
 			{
-				z = std::stoi(junk);
+				StepParameters currentParams;
-				input>>a;
+				input >> junk >> junk;
-				zvec.push_back(z);
+				currentParams.rxnType = StringToRxnType(junk);
-				avec.push_back(a);
+				if(currentParams.rxnType == RxnType::Reaction)
+				{
+					input >> junk;
+					for(int i=0; i<3; i++)
+					{
+						input >> z >> a;
+						currentParams.Z.push_back(z);
+						currentParams.A.push_back(a);
+					}
+					input >> junk;
+					input >> junk >> currentParams.meanBeamEnergy;
+					input >> junk >> currentParams.sigmaBeamEnergy;
+					input >> junk >> junk;
+					currentParams.thetaType = StringToRxnThetaType(junk);
+					input >> junk >> currentParams.thetaMin;
+					input >> junk >> currentParams.thetaMax;
+					input >> junk >> currentParams.phiMin;
+					input >> junk >> currentParams.phiMax;
+					input >> junk >> currentParams.meanResidualEx;
+					input >> junk >> currentParams.sigmaResidualEx;
+					params.push_back(currentParams);
+				}
+				else if(currentParams.rxnType == RxnType::Decay)
+				{
+					input >> junk;
+					for(int i=0; i<2; i++)
+					{
+						input >> z >> a;
+						currentParams.Z.push_back(z);
+						currentParams.A.push_back(a);
+					}
+					input >> junk;
+					input >> junk >> currentParams.phiMin;
+					input >> junk >> currentParams.phiMax;
+					input >> junk >> currentParams.meanResidualEx;
+					input >> junk >> currentParams.sigmaResidualEx;
+					input >> junk >> currentParams.angularDistFile;
+					params.push_back(currentParams);
+				}
+				else
+				{
+					std::cerr << "Invalid reaction information at MaskApp::LoadConfig!" << std::endl;
+					return false;
+				}
 			}
 		}
 
-		m_system = CreateSystem(zvec, avec);
+		m_system = CreateSystem(params);
-		if(m_system == nullptr)
+		if(m_system == nullptr || !m_system->IsValid())
 		{
+			std::cerr<<"Param size: "<<params.size()<<std::endl;
 			std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
 			return false;
 		}
+		std::getline(input, junk);
+		std::getline(input, junk);
 
+		LayeredTarget target;
 		int nlayers;
 		double thickness;
-		getline(input, junk);
+		std::vector<int> avec, zvec, svec;
-		getline(input, junk);
+		int s;
 		input>>junk>>nlayers;
 		for(int i=0; i<nlayers; i++) 
 		{
@@ -80,42 +127,26 @@ namespace Mask {
 				input>>z>>a>>s;
 				zvec.push_back(z); avec.push_back(a); svec.push_back(s);
 			}
-			m_system->AddTargetLayer(zvec, avec, svec, thickness);
+			target.AddLayer(zvec, avec, svec, thickness);
 			input>>junk;
 		}
+		m_system->SetLayeredTarget(target);
+
+		std::cout<<"Outputing data to file: "<<m_outputName<<std::endl;
 		std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
-	
-		double par1, par2;
-		std::string dfile1, dfile2;
-		std::string thetaTypeString;
-		getline(input, junk);
-		getline(input, junk);
-	
-		input>>junk>>m_nsamples;
-		input>>junk>>par1>>junk>>par2;
-		m_system->SetBeamDistro(par1, par2);
-		input>>junk>>thetaTypeString;
-		m_system->SetReactionThetaType(StringToRxnThetaType(thetaTypeString));
-		
-		input>>junk>>par1>>junk>>par2;
-		m_system->SetTheta1Range(par1, par2);
-		input>>junk>>par1>>junk>>par2;
-		m_system->SetPhi1Range(par1, par2);
-		input>>junk>>par1>>junk>>par2;
-		m_system->SetExcitationDistro(par1, par2);
-		input>>junk>>dfile1;
-		input>>junk>>dfile2;
-		m_system->SetDecay1Distribution(dfile1);
-		m_system->SetDecay2Distribution(dfile2);
-	
 		std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
 	
 		return true;
 	}
 	
-	bool MaskApp::SaveConfig(const std::string& filename) { return true; }
+	//Not implemented... yet!
+	bool MaskApp::SaveConfig(const std::string& filename)
+	{ 
+		return true;
+	}
 	
-	void MaskApp::Run() {
+	void MaskApp::Run()
+	{
 		std::cout<<"Running simulation..."<<std::endl;
 		if(m_system == nullptr) 
 		{

src/Mask/OneStepSystem.cpp

@@ -5,51 +5,56 @@
 
 namespace Mask {
 
-	OneStepSystem::OneStepSystem() :
+	OneStepSystem::OneStepSystem(const std::vector<StepParameters>& params) :
 		ReactionSystem()
 	{
 		m_nuclei.resize(4);
-	}
+		Init(params);
-	
-	OneStepSystem::OneStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
-		ReactionSystem()
-	{
-		m_nuclei.resize(4);
-		SetNuclei(z, a);
 	}
 	
 	OneStepSystem::~OneStepSystem() {}
 	
-	bool OneStepSystem::SetNuclei(const std::vector<int>& z, const std::vector<int>& a)
+	void OneStepSystem::Init(const std::vector<StepParameters>& params)
 	{
-		if(z.size() != a.size() || z.size() != 3)
+		if(params.size() != 1 || params[0].rxnType != RxnType::Reaction ||
-			return false;
+		   params[0].Z.size() != 3 || params[0].A.size() != 3)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at OneStepSystem::Init(), does not match OneStep signature!" << std::endl;
+			return;
+		}
 
-		int zr = z[0] + z[1] - z[2];
+		const StepParameters& step1Params = params[0];
-		int ar = a[0] + a[1] - a[2];
 
-		m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
+		//Set nuclei
-		m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
+
-		m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
+		int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2];
+		int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2];
+
+		m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
+		m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //projectile
+		m_nuclei[2] = CreateNucleus(step1Params.Z[2], step1Params.A[2]); //ejectile
 		m_nuclei[3] = CreateNucleus(zr, ar); //residual
 
 		m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
 		SetSystemEquation();
-		return true;
+
+		//Set sampling parameters
+
+		AddBeamDistribution(step1Params.meanBeamEnergy, step1Params.sigmaBeamEnergy);
+		m_step1.SetEjectileThetaType(step1Params.thetaType);
+		AddThetaRange(step1Params.thetaMin, step1Params.thetaMax);
+		AddPhiRange(step1Params.phiMin, step1Params.phiMax);
+		AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
 	}
 	
-	std::vector<Nucleus>* OneStepSystem::GetNuclei()
+	void OneStepSystem::SetLayeredTarget(const LayeredTarget& target)
 	{
-		return &m_nuclei;
+		m_target = target;
-	}
-	
-	void OneStepSystem::LinkTarget()
-	{
-		m_step1.SetLayeredTarget(&m_target);
-	
 		m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
 		if(m_rxnLayer != m_target.GetNumberOfLayers())
 		{
+			m_step1.SetLayeredTarget(&m_target);
 			m_step1.SetRxnLayer(m_rxnLayer);
 			m_isTargetSet = true;
 		}
@@ -67,18 +72,14 @@ namespace Mask {
 		m_sysEquation = stream.str();
 	}
 	
-	void OneStepSystem::RunSystem() {
+	void OneStepSystem::RunSystem()
-		//Link up the target if it hasn't been done yet
+	{
-		if(!m_isTargetSet)
-		{
-			LinkTarget();
-		}
-	
 		//Sample parameters
-		double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
-		double rxnTheta = std::acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator()));
+		double bke = (m_beamDistributions[0])(gen);
-		double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
+		double rxnTheta = std::acos((m_thetaRanges[0])(gen));
-		double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
+		double rxnPhi = (m_phiRanges[0])(gen);
+		double residEx = (m_exDistributions[0])(gen);
 		
 		m_step1.SetBeamKE(bke);
 		m_step1.SetPolarRxnAngle(rxnTheta);

src/Mask/OneStepSystem.h

@@ -5,24 +5,20 @@
 
 namespace Mask {
 
-	class OneStepSystem: public ReactionSystem {
+	class OneStepSystem: public ReactionSystem
+	{
 	public:
-		OneStepSystem();
+		OneStepSystem(const std::vector<StepParameters>& params);
-		OneStepSystem(const std::vector<int>& z, const std::vector<int>& a);
 		~OneStepSystem();
 	
-		bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
+		virtual void SetLayeredTarget(const LayeredTarget& target) override;
 		void RunSystem() override;
-		std::vector<Nucleus>* GetNuclei() override;
-	
-		virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
 	
 	private:
-		void LinkTarget() override;
+		void Init(const std::vector<StepParameters>& params);
 		void SetSystemEquation() override;
 	
 		Reaction m_step1;
-	
 	};
 
 }

src/Mask/Reaction.cpp

@@ -187,14 +187,15 @@ namespace Mask {
 	//Calculate in CM, where decay is isotropic
 	void Reaction::CalculateDecay()
 	{
-		double Q = m_target->vec4.M() - m_ejectile->groundStateMass - m_residual->groundStateMass;
+		double residualMass = m_residual->groundStateMass + m_ex;
+		double Q = m_target->vec4.M() - m_ejectile->groundStateMass - residualMass;
 		if(Q < 0)
 			throw QValueException();
 	
 		ROOT::Math::Boost boost(m_target->vec4.BoostToCM());
 		m_target->vec4 = boost*m_target->vec4;
 		double ejectE_cm = (m_ejectile->groundStateMass*m_ejectile->groundStateMass - 
-						   m_residual->groundStateMass*m_residual->groundStateMass + m_target->vec4.E()*m_target->vec4.E()) /
+						   residualMass*residualMass + m_target->vec4.E()*m_target->vec4.E()) /
 					       (2.0*m_target->vec4.E());
 		double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - m_ejectile->groundStateMass*m_ejectile->groundStateMass);
 	

src/Mask/ReactionSystem.cpp

@@ -10,41 +10,91 @@
 namespace Mask {
 
 	ReactionSystem::ReactionSystem() :
-		m_beamDist(nullptr), m_theta1Range(nullptr), m_phi1Range(nullptr), m_exDist(nullptr), m_isTargetSet(false),
+		m_isTargetSet(false), m_isValid(true), m_rxnLayer(0), m_sysEquation("")
-		m_rxnLayer(0), m_sysEquation("")
 	{
 	}
 	
 	ReactionSystem::~ReactionSystem()
 	{
-		delete m_beamDist;
-		delete m_theta1Range;
-		delete m_phi1Range;
-		delete m_exDist;
 	}
 
-	void ReactionSystem::AddTargetLayer(const std::vector<int>& zt, const std::vector<int>& at, const std::vector<int>& stoich, double thickness)
+	ReactionSystem* CreateSystem(const std::vector<StepParameters>& params)
 	{
-		m_target.AddLayer(zt, at, stoich, thickness);
+		switch(params.size())
-	}
-
-	ReactionSystem* CreateSystem(const std::vector<int>& z, const std::vector<int>& a)
-	{
-		if(z.size() != a.size())
 		{
-			std::cerr<<"Size of Z list does not equal size of A list!"<<std::endl;
+			case 1:
-			return nullptr;
+			{
-		}
+				if(params[0].rxnType == RxnType::Decay)
-
+					return new DecaySystem(params);
-		switch(a.size())
+				else if (params[0].rxnType == RxnType::Reaction)
-		{
+					return new OneStepSystem(params);
-			case RxnSize::DecaySize: return new DecaySystem(z, a);
+			}
-			case RxnSize::OneStepSize: return new OneStepSystem(z, a);
+			case 2: return new TwoStepSystem(params);
-			case RxnSize::TwoStepSize: return new TwoStepSystem(z, a);
+			case 3: return new ThreeStepSystem(params);
-			case RxnSize::ThreeStepSize: return new ThreeStepSystem(z, a);
 		}
 
 		return nullptr;
 	}
 
+	/*Set sampling parameters*/
+	//Should only ever really be one of these, but they cannot be set a priori
+	void ReactionSystem::AddBeamDistribution(double mean, double sigma)
+	{
+		if(mean == -1.0 || sigma == -1.0)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ReactionSystem::AddBeamDistribution(), distribution is invalid -> mean: " << mean
+					  << " sigma: " << sigma << std::endl;
+			return;
+		}
+		m_beamDistributions.emplace_back(mean, sigma);
+	}
+
+	//Again only really one of these, but same issue as beam. This is the sampling of the primary reaction/decay (step 1)
+	void ReactionSystem::AddThetaRange(double min, double max)
+	{
+		if(min == -1.0 || max == -1.0 || min > max)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ReactionSystem::AddThetaRange(), range is invalid -> min: " << min
+					  << " max: " << max << std::endl;
+			return;
+		}
+		m_thetaRanges.emplace_back(std::cos(min*s_deg2rad), std::cos(max*s_deg2rad));
+	}
+
+	void ReactionSystem::AddPhiRange(double min, double max)
+	{
+		if(min == -1.0 || max == -1.0 || min > max)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ReactionSystem::AddPhiRange(), range is invalid -> min: " << min
+					  << " max: " << max << std::endl;
+			return;
+		}
+		m_phiRanges.emplace_back(min*s_deg2rad, max*s_deg2rad);
+	}
+
+	//Each reaction step can generate an excited nucleus (for a reaction can make an excited residual, decay can make an excited
+	//breakup2 or "heavy")
+	void ReactionSystem::AddExcitationDistribution(double mean, double sigma)
+	{
+		if(mean == -1.0 || sigma == -1.0)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ReactionSystem::AddExcitationDistribution(), distribution is invalid -> mean: " << mean
+					  << " sigma: " << sigma << std::endl;
+			return;
+		}
+		m_exDistributions.emplace_back(mean, sigma);
+	}
+
+	//Apply angular distributions to decay products. Here apply them to the breakup1 or "light" fragment
+	void ReactionSystem::AddDecayAngularDistribution(const std::string& filename)
+	{
+		if(filename.empty())
+			m_decayAngularDistributions.emplace_back();
+		else
+			m_decayAngularDistributions.emplace_back(filename);
+	}
 }

src/Mask/ReactionSystem.h

@@ -12,75 +12,69 @@
 #include "Reaction.h"
 #include "KinematicsExceptions.h"
 #include "RxnType.h"
+#include "AngularDistribution.h"
 #include <vector>
 #include <random>
 
 namespace Mask {
 
+	struct StepParameters
+	{
+		RxnType rxnType = RxnType::None;
+		std::vector<int> Z;
+		std::vector<int> A;
+		double meanBeamEnergy = -1.0;
+		double sigmaBeamEnergy = -1.0;
+		RxnThetaType thetaType = RxnThetaType::None;
+		double thetaMin = -1.0;
+		double thetaMax = -1.0;
+		double phiMin = -1.0;
+		double phiMax = -1.0;
+		double meanResidualEx = -1.0;
+		double sigmaResidualEx = -1.0;
+		std::string angularDistFile;
+	};
+
 	class ReactionSystem
 	{
 	public:
 		ReactionSystem();
 		virtual ~ReactionSystem();
 
-		virtual bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) = 0;
+		virtual void SetLayeredTarget(const LayeredTarget& target) = 0;
 		virtual void RunSystem() = 0;
-		virtual std::vector<Nucleus>* GetNuclei() = 0;
-		virtual void SetReactionThetaType(RxnThetaType type) {}
-		virtual void SetDecay1Distribution(const std::string& filename) {}
-		virtual void SetDecay2Distribution(const std::string& filename) {}
-	
-		void AddTargetLayer(const std::vector<int>& zt, const std::vector<int>& at, const std::vector<int>& stoich, double thickness);
-	
-		/*Set sampling parameters*/
-		void SetBeamDistro(double mean, double sigma)
-		{
-			if(m_beamDist)
-				delete m_beamDist;
-			m_beamDist = new std::normal_distribution<double>(mean, sigma); 
-		}
-	
-		void SetTheta1Range(double min, double max)
-		{ 
-			if(m_theta1Range)
-				delete m_theta1Range;
-			m_theta1Range = new std::uniform_real_distribution<double>(std::cos(min*s_deg2rad), std::cos(max*s_deg2rad)); 
-		}
-	
-		void SetPhi1Range(double min, double max)
-		{
-			if(m_phi1Range)
-				delete m_phi1Range;
-			m_phi1Range = new std::uniform_real_distribution<double>(min*s_deg2rad, max*s_deg2rad); 
-		}
-	
-		void SetExcitationDistro(double mean, double sigma)
-		{
-			if(m_exDist)
-				delete m_exDist;
-			m_exDist = new std::normal_distribution<double>(mean, sigma); 
-		}
 
+		std::vector<Nucleus>* GetNuclei() { return &m_nuclei; }
 		const std::string& GetSystemEquation() const { return m_sysEquation; }
+		bool IsValid() const { return m_isValid; }
 
 	protected:
-		virtual void LinkTarget() = 0;
 		virtual void SetSystemEquation() = 0;
 
+		void AddBeamDistribution(double mean, double sigma);
+		void AddThetaRange(double min, double max);
+		void AddPhiRange(double min, double max);
+		void AddExcitationDistribution(double mean, double sigma);
+		void AddDecayAngularDistribution(const std::string& filename);
+		
 		LayeredTarget m_target;
 	
 		//Sampling information
-		std::normal_distribution<double> *m_beamDist, *m_exDist;
+		std::vector<std::normal_distribution<double>> m_beamDistributions, m_exDistributions;
-		std::uniform_real_distribution<double> *m_theta1Range, *m_phi1Range; 
+		std::vector<std::uniform_real_distribution<double>> m_thetaRanges, m_phiRanges;
+		std::vector<AngularDistribution> m_decayAngularDistributions;
 
 		bool m_isTargetSet;
+		bool m_isValid;
+
 		std::size_t m_rxnLayer;
 		std::string m_sysEquation;
 		std::vector<Nucleus> m_nuclei;
+
 		static constexpr double s_deg2rad = M_PI/180.0;
 	};
 
-	ReactionSystem* CreateSystem(const std::vector<int>& z, const std::vector<int>& a);
+	ReactionSystem* CreateSystem(const std::vector<StepParameters>& params);
 }
 
 #endif

src/Mask/RxnType.h

@@ -7,11 +7,9 @@ namespace Mask {
 
 	enum class RxnType
 	{
-		PureDecay=0,
+		Decay,
-		OneStepRxn=1,
+		Reaction,
-		TwoStepRxn=2,
+		None
-		ThreeStepRxn=3,
-		None=4
 	};
 
 	enum class RxnThetaType
@@ -29,44 +27,16 @@ namespace Mask {
 		ThreeStepSize = 5
 	};
 
-	static RxnType GetRxnTypeFromString(const std::string& type_str)
+	static RxnType StringToRxnType(const std::string& type)
 	{
-		if (type_str == "PureDecay")
+		if (type == "Decay")
-			return RxnType::PureDecay;
+			return RxnType::Decay;
-		else if (type_str == "OneStepRxn")
+		else if (type == "Reaction")
-			return RxnType::OneStepRxn;
+			return RxnType::Reaction;
-		else if (type_str == "TwoStepRxn")
-			return RxnType::TwoStepRxn;
-		else if (type_str == "ThreeStepRxn")
-			return RxnType::ThreeStepRxn;
 		else
 			return RxnType::None;
 	}
 
-	static std::string GetStringFromRxnType(RxnType type)
-	{
-		switch(type)
-		{
-			case RxnType::PureDecay: return "PureDecay";
-			case RxnType::OneStepRxn: return "OneStepRxn";
-			case RxnType::TwoStepRxn: return "TwoStepRxn";
-			case RxnType::ThreeStepRxn: return "ThreeStepRxn";
-			case RxnType::None : return "None";
-		}
-
-		return "None";
-	}
-
-	static RxnType GetRxnTypeFromInt(uint32_t value)
-	{
-		return static_cast<RxnType>(value);
-	}
-
-	static uint32_t GetIntFromRxnType(RxnType type)
-	{
-		return static_cast<uint32_t>(type);
-	}
-
 	static RxnThetaType StringToRxnThetaType(const std::string& type)
 	{
 		if(type == "CenterOfMass")

src/Mask/ThreeStepSystem.cpp

@@ -4,40 +4,59 @@
 
 namespace Mask {
 	
-	ThreeStepSystem::ThreeStepSystem() :
+	ThreeStepSystem::ThreeStepSystem(const std::vector<StepParameters>& params) :
-		ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
+		ReactionSystem()
 	{
 		m_nuclei.resize(8);
-	}
+		Init(params);
-	
-	ThreeStepSystem::ThreeStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
-		ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
-	{
-		m_nuclei.resize(8);
-		SetNuclei(z, a);
 	}
 	
 	ThreeStepSystem::~ThreeStepSystem() {}
 	
-	bool ThreeStepSystem::SetNuclei(const std::vector<int>& z, const std::vector<int>& a)
+	void ThreeStepSystem::Init(const std::vector<StepParameters>& params)
 	{
-		if(z.size() != a.size() || z.size() < 5)
+		if(params.size() != 3 || params[0].rxnType != RxnType::Reaction || params[1].rxnType != RxnType::Decay ||
-			return false;
+		  params[2].rxnType != RxnType::Decay || params[0].Z.size() != 3 || params[0].A.size() != 3 || params[1].Z.size() != 2 ||
+		  params[1].A.size() != 2 || params[2].Z.size() != 2 || params[2].A.size() != 2)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ThreeStepSystem::Init(), does not match ThreeStep signature!" << std::endl;
+			return;
+		}
 
-		int zr = z[0] + z[1] - z[2];
+		const StepParameters& step1Params = params[0];
-		int zb2 = zr - z[3];
+		const StepParameters& step2Params = params[1];
-		int zb4 = zb2 - z[4];
+		const StepParameters& step3Params = params[2];
-		int ar = a[0] + a[1] - a[2];
-		int ab2 = ar - a[3];
-		int ab4 = ab2 - a[4];
 
-		m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
+		//Setup nuclei
-		m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
+		int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2];
-		m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
+		int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2];
+		if(zr != step2Params.Z[0] || ar != step2Params.A[0])
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step one and step two are not sequential! Step one recoil (Z,A): ("
+					  << zr << "," << ar << ") Step two target (Z,A): (" << step2Params.Z[0] << "," << step2Params.A[0] << ")" <<std::endl;
+			return;  
+		}
+		int zb2 = step2Params.Z[0] - step2Params.Z[1];
+		int ab2 = step2Params.A[0] - step2Params.A[1];
+		if(zb2 != step3Params.Z[0] || ab2 != step3Params.A[0])
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step two and step three are not sequential! Step two heavy (Z,A): ("
+					  << zb2 << "," << ab2 << ") Step three target (Z,A): (" << step3Params.Z[0] << "," << step3Params.A[0] << ")" <<std::endl;
+			return;
+		}
+		int zb4 = step3Params.Z[0] - step3Params.Z[1];
+		int ab4 = step3Params.A[0] - step3Params.A[1];
+
+		m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
+		m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //projectile
+		m_nuclei[2] = CreateNucleus(step1Params.Z[2], step1Params.A[2]); //ejectile
 		m_nuclei[3] = CreateNucleus(zr, ar); //residual
-		m_nuclei[4] = CreateNucleus(z[3], a[3]); //breakup1
+		m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
 		m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
-		m_nuclei[5] = CreateNucleus(z[4], a[4]); //breakup3
+		m_nuclei[5] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
 		m_nuclei[5] = CreateNucleus(zb4, ab4); //breakup4
 
 		m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
@@ -45,28 +64,40 @@ namespace Mask {
 		m_step3.BindNuclei(&(m_nuclei[5]), nullptr, &(m_nuclei[6]), &(m_nuclei[7]));
 	
 		SetSystemEquation();
-		return true;
+
+		//Step one sampling parameters
+		AddBeamDistribution(step1Params.meanBeamEnergy, step1Params.sigmaBeamEnergy);
+		m_step1.SetEjectileThetaType(step1Params.thetaType);
+		AddThetaRange(step1Params.thetaMin, step1Params.thetaMax);
+		AddPhiRange(step1Params.phiMin, step1Params.phiMax);
+		AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
+
+		//Step two sampling parameters
+		AddPhiRange(step2Params.phiMin, step2Params.phiMax);
+		AddDecayAngularDistribution(step2Params.angularDistFile);
+		AddExcitationDistribution(step2Params.meanResidualEx, step2Params.sigmaResidualEx);
+
+		//Step three sampling parameters
+		AddPhiRange(step3Params.phiMin, step3Params.phiMax);
+		AddDecayAngularDistribution(step3Params.angularDistFile);
+		AddExcitationDistribution(step3Params.meanResidualEx, step3Params.sigmaResidualEx);
 	}
 
-	std::vector<Nucleus>* ThreeStepSystem::GetNuclei()
+	void ThreeStepSystem::SetLayeredTarget(const LayeredTarget& target)
 	{
-		return &m_nuclei;
+		m_target = target;
-	}
-
-	void ThreeStepSystem::LinkTarget()
-	{
-		m_step1.SetLayeredTarget(&m_target);
-		m_step2.SetLayeredTarget(&m_target);
-		m_step3.SetLayeredTarget(&m_target);
-	
 		m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
 		if(m_rxnLayer != m_target.GetNumberOfLayers())
 		{
+			m_step1.SetLayeredTarget(&m_target);
+			m_step2.SetLayeredTarget(&m_target);
+			m_step3.SetLayeredTarget(&m_target);
 			m_step1.SetRxnLayer(m_rxnLayer);
 			m_step2.SetRxnLayer(m_rxnLayer);
 			m_step3.SetRxnLayer(m_rxnLayer);
 			m_isTargetSet = true;
-		} else
+		}
+		else
 			throw ReactionLayerException();
 	}
 	
@@ -84,22 +115,22 @@ namespace Mask {
 		m_sysEquation = stream.str();
 	}
 	
-	void ThreeStepSystem::RunSystem() {
+	void ThreeStepSystem::RunSystem()
-		//Link up the target if it hasn't been done yet
+	{
-		if(!m_isTargetSet)
-			LinkTarget();
-	
 		//Sample parameters
-		double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
-		double rxnTheta = acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator()));
+		double bke = m_beamDistributions[0](gen);
-		double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
+		double rxnTheta = std::acos((m_thetaRanges[0])(gen));
-		double decay1costheta = m_step2Distribution.GetRandomCosTheta();
+		double rxnPhi = m_phiRanges[0](gen);
+		double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
 		double decay1Theta = std::acos(decay1costheta);
-		double decay1Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
+		double decay1Phi = m_phiRanges[1](gen);
-		double decay2costheta = m_step3Distribution.GetRandomCosTheta();
+		double decay2costheta = m_decayAngularDistributions[1].GetRandomCosTheta();
 		double decay2Theta = std::acos(decay2costheta);
-		double decay2Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
+		double decay2Phi = m_phiRanges[2](gen);
-		double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
+		double residEx = m_exDistributions[0](gen);
+		double decay1Ex = m_exDistributions[1](gen);
+		double decay2Ex = m_exDistributions[2](gen);
 	
 		m_step1.SetBeamKE(bke);
 		m_step1.SetPolarRxnAngle(rxnTheta);
@@ -108,9 +139,11 @@ namespace Mask {
 	
 		m_step2.SetPolarRxnAngle(decay1Theta);
 		m_step2.SetAzimRxnAngle(decay1Phi);
+		m_step2.SetExcitation(decay1Ex);
 	
 		m_step3.SetPolarRxnAngle(decay2Theta);
 		m_step3.SetAzimRxnAngle(decay2Phi);
+		m_step3.SetExcitation(decay2Ex);
 	
 		m_step1.Calculate();
 	

src/Mask/ThreeStepSystem.h

@@ -6,35 +6,21 @@
 
 namespace Mask {
 
-	class ThreeStepSystem : public ReactionSystem {
+	class ThreeStepSystem : public ReactionSystem
+	{
 	public:
-		ThreeStepSystem();
+		ThreeStepSystem(const std::vector<StepParameters>& params);
-		ThreeStepSystem(const std::vector<int>& z, const std::vector<int>& a);
 		~ThreeStepSystem();
-		bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
-		void RunSystem() override;
-		std::vector<Nucleus>* GetNuclei() override;
 
-		virtual void SetDecay1Distribution(const std::string& filename) override { m_step2Distribution.ReadDistributionFile(filename); };
+		virtual void SetLayeredTarget(const LayeredTarget& target) override;
-		virtual void SetDecay2Distribution(const std::string& filename) override { m_step3Distribution.ReadDistributionFile(filename); };
+		virtual void RunSystem() override;
-	
-		virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
-	
-		int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); };
-		int GetDecay2AngularMomentum(){ return m_step3Distribution.GetL(); };
-		double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); };
-		double GetDecay2BranchingRatio(){ return m_step3Distribution.GetBranchingRatio(); };
 		
 	protected:
-		void LinkTarget() override;
+		void Init(const std::vector<StepParameters>& params);
 		void SetSystemEquation() override;
 	
-		std::uniform_real_distribution<double> m_phi2Range;
-	
 		Reaction m_step1, m_step2, m_step3;
 	
-		AngularDistribution m_step2Distribution, m_step3Distribution;
-	
 	};
 
 }

src/Mask/TwoStepSystem.cpp

@@ -6,57 +6,73 @@
 
 namespace Mask {
 	
-	TwoStepSystem::TwoStepSystem() :
+	TwoStepSystem::TwoStepSystem(const std::vector<StepParameters>& params) :
-		ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
+		ReactionSystem()
 	{
 		m_nuclei.resize(6);
-	}
+		Init(params);
-	
-	TwoStepSystem::TwoStepSystem(const std::vector<int>& z, const std::vector<int>& a) :
-		ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
-	{
-		m_nuclei.resize(6);
-		SetNuclei(z, a);
 	}
 	
 	TwoStepSystem::~TwoStepSystem() {}
 	
-	bool TwoStepSystem::SetNuclei(const std::vector<int>&z, const std::vector<int>& a)
+	void TwoStepSystem::Init(const std::vector<StepParameters>& params)
 	{
-		if(z.size() != a.size() || z.size() != 4)
+		if(params.size() != 2 || params[0].rxnType != RxnType::Reaction || params[1].rxnType != RxnType::Decay ||
-			return false;
+		   params[0].Z.size() != 3 || params[0].A.size() != 3 || params[1].Z.size() != 2 || params[1].A.size() != 2)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at TwoStepSystem::Init(), does not match TwoStep signature!" << std::endl;
+			return;
+		}
 
-		int zr = z[0] + z[1] - z[2];
+		const StepParameters& step1Params = params[0];
-		int ar = a[0] + a[1] - a[2];
+		const StepParameters& step2Params = params[1];
-		int zb = zr - z[3];
-		int ab = ar - a[3];
 
-		m_nuclei[0] = CreateNucleus(z[0], a[0]); //target
+		//Setup nuclei
-		m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile
+		int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2];
-		m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile
+		int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2];
+		if(zr != step2Params.Z[0] || ar != step2Params.A[0])
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at TwoStepSystem::Init(), step one and step two are not sequential! Step one recoil (Z,A): ("
+					  << zr << "," << ar << ") Step two target (Z,A): (" << step2Params.Z[0] << "," << step2Params.A[0] << ")" <<std::endl;
+			return;  
+		}
+		int zb = step2Params.Z[0] - step2Params.Z[1];
+		int ab = step2Params.A[0] - step2Params.A[1];
+
+		m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
+		m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //projectile
+		m_nuclei[2] = CreateNucleus(step1Params.Z[2], step1Params.A[2]); //ejectile
 		m_nuclei[3] = CreateNucleus(zr, ar); //residual
-		m_nuclei[4] = CreateNucleus(z[3], a[3]); //breakup1
+		m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
 		m_nuclei[5] = CreateNucleus(zb, ab); //breakup2
 
 		m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
 		m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
 		SetSystemEquation();
-		return true;
+
+		//Step one sampling parameters
+		AddBeamDistribution(step1Params.meanBeamEnergy, step1Params.sigmaBeamEnergy);
+		m_step1.SetEjectileThetaType(step1Params.thetaType);
+		AddThetaRange(step1Params.thetaMin, step1Params.thetaMax);
+		AddPhiRange(step1Params.phiMin, step1Params.phiMax);
+		AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
+
+		//Step two sampling parameters
+		AddPhiRange(step2Params.phiMin, step2Params.phiMax);
+		AddDecayAngularDistribution(step2Params.angularDistFile);
+		AddExcitationDistribution(step2Params.meanResidualEx, step2Params.sigmaResidualEx);
 	}
 	
-	std::vector<Nucleus>* TwoStepSystem::GetNuclei()
+	void TwoStepSystem::SetLayeredTarget(const LayeredTarget& target)
 	{
-		return &m_nuclei;
+		m_target = target;
-	}
-	
-	void TwoStepSystem::LinkTarget()
-	{
-		m_step1.SetLayeredTarget(&m_target);
-		m_step2.SetLayeredTarget(&m_target);
-	
 		m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
 		if(m_rxnLayer != m_target.GetNumberOfLayers())
 		{
+			m_step1.SetLayeredTarget(&m_target);
+			m_step2.SetLayeredTarget(&m_target);
 			m_step1.SetRxnLayer(m_rxnLayer);
 			m_step2.SetRxnLayer(m_rxnLayer);
 			m_isTargetSet = true;
@@ -79,18 +95,16 @@ namespace Mask {
 	
 	void TwoStepSystem::RunSystem()
 	{
-		//Link up the target if it hasn't been done yet
-		if(!m_isTargetSet)
-			LinkTarget();
-	
 		//Sample parameters
-		double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator());
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
-		double rxnTheta = acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator()));
+		double bke = (m_beamDistributions[0])(gen);
-		double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator());
+		double rxnTheta = std::acos((m_thetaRanges[0])(gen));
-		double decay1costheta = m_step2Distribution.GetRandomCosTheta();
+		double rxnPhi = (m_phiRanges[0])(gen);
+		double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
 		double decay1Theta = std::acos(decay1costheta);
-		double decay1Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator());
+		double decay1Phi = m_phiRanges[1](gen);
-		double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator());
+		double residEx = (m_exDistributions[0])(gen);
+		double decay2Ex = m_exDistributions[1](gen);
 	
 		m_step1.SetBeamKE(bke);
 		m_step1.SetPolarRxnAngle(rxnTheta);
@@ -99,6 +113,7 @@ namespace Mask {
 	
 		m_step2.SetPolarRxnAngle(decay1Theta);
 		m_step2.SetAzimRxnAngle(decay1Phi);
+		m_step2.SetExcitation(decay2Ex);
 		
 		m_step1.Calculate();
 	

src/Mask/TwoStepSystem.h

@@ -9,30 +9,17 @@ namespace Mask {
 	class TwoStepSystem : public ReactionSystem
 	{
 	public:
-		TwoStepSystem();
+		TwoStepSystem(const std::vector<StepParameters>& params);
-		TwoStepSystem(const std::vector<int>& z, const std::vector<int>& a);
 		~TwoStepSystem();
-		bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) override;
-		void RunSystem() override;
-		std::vector<Nucleus>* GetNuclei() override;
 
-		virtual void SetDecay1Distribution(const std::string& filename) override { m_step2Distribution.ReadDistributionFile(filename); };
+		virtual void SetLayeredTarget(const LayeredTarget& target) override;
-	
+		virtual void RunSystem() override;
-		virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
-
-		int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); };
-		double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); };
 	
 	private:
-		void LinkTarget() override;
+		void Init(const std::vector<StepParameters>& params);
 		void SetSystemEquation() override;
 
-		std::uniform_real_distribution<double> m_phi2Range;
-		
 		Reaction m_step1, m_step2;
-	
-		AngularDistribution m_step2Distribution;
-	
 	};
 
 }

src/Plotters/RootPlotter.cpp

@@ -75,18 +75,18 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc)
 	std::string angdist_name = sym + modifier +"_angDist";
 	std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts";
 	
-	MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
+	MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
-	MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.GetKE(), 4);
+	MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.GetKE(), 4);
-	MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
+	MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
-	MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
+	MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
-	MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM));
+	MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
 	if(nuc.isDetected)
 	{
-		MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
+		MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
-		MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4);
+		MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4);
-		MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
+		MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
-		MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
+		MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
-		MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM));
+		MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
 	}
 	
 }