diff --git a/include/DecaySystem.h b/include/DecaySystem.h
new file mode 100644
index 0000000..d435561
--- /dev/null
+++ b/include/DecaySystem.h
@@ -0,0 +1,41 @@
+#ifndef DECAYSYSTEM_H
+#define DECAYSYSTEM_H
+
+#include "ReactionSystem.h"
+#include "AngularDistribution.h"
+
+namespace Mask {
+
+class DecaySystem: public ReactionSystem {
+public:
+	DecaySystem();
+	DecaySystem(std::vector<int>& z, std::vector<int>& a);
+	~DecaySystem();
+
+	bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
+	void RunSystem() override;
+	void SetRandomGenerator(TRandom3* gen) override;
+
+	inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
+
+	inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
+	inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
+	inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
+
+	inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
+	inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
+
+
+private:
+	void LinkTarget() override;
+	void SetSystemEquation() override;
+
+	Reaction step1;
+
+	AngularDistribution decay1dist;
+	
+};
+
+};
+
+#endif
\ No newline at end of file
diff --git a/include/Kinematics.h b/include/Kinematics.h
index 24dc017..a015cee 100644
--- a/include/Kinematics.h
+++ b/include/Kinematics.h
@@ -2,6 +2,8 @@
 #define KINEMATICS_H
 
 #include "ReactionSystem.h"
+#include "DecaySystem.h"
+#include "OneStepSystem.h"
 #include "TwoStepSystem.h"
 #include "ThreeStepSystem.h"
 #include "Plotter.h"
@@ -32,7 +34,7 @@ public:
 	inline void SetTreeFlag() { save_tree_flag = true; };
 	inline void SetPlotterFlag() { do_plotter_flag = true; };
 	inline int GetNumberOfSamples() { return m_nsamples; };
-	inline const char* GetSystemName() { return sys == nullptr ? "" : sys->GetSystemEquation(); };
+	inline const char* GetSystemName() { return sys == nullptr ? "" : sys->GetSystemEquation().c_str(); };
 	inline const char* GetOutputName() { return m_outfile_name.c_str(); };
 	inline int GetReactionType() { return m_rxn_type; };
 	void Run();
diff --git a/include/OneStepSystem.h b/include/OneStepSystem.h
new file mode 100644
index 0000000..6b58fc4
--- /dev/null
+++ b/include/OneStepSystem.h
@@ -0,0 +1,33 @@
+#ifndef ONESTEPSYSTEM_H
+#define ONESTEPSYSTEM_H
+
+#include "ReactionSystem.h"
+
+namespace Mask {
+
+class OneStepSystem: public ReactionSystem {
+public:
+	OneStepSystem();
+	OneStepSystem(std::vector<int>& z, std::vector<int>& a);
+	~OneStepSystem();
+
+	bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
+	void RunSystem() override;
+
+	inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
+	inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
+	inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
+	inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
+	inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
+
+private:
+	void LinkTarget() override;
+	void SetSystemEquation() override;
+
+	Reaction step1;
+
+};
+
+};
+
+#endif
\ No newline at end of file
diff --git a/include/ReactionSystem.h b/include/ReactionSystem.h
index 10d30dc..4aa885f 100644
--- a/include/ReactionSystem.h
+++ b/include/ReactionSystem.h
@@ -10,7 +10,7 @@
 #define REACTIONSYSTEM_H
 
 #include "Reaction.h"
-#include "AngularDistribution.h"
+#include "KinematicsExceptions.h"
 #include <vector>
 #include <TRandom3.h>
 
@@ -19,47 +19,32 @@ namespace Mask {
 class ReactionSystem {
 public:
 	ReactionSystem();
-	ReactionSystem(std::vector<int>& z, std::vector<int>& a);
 	virtual ~ReactionSystem();
 
-	virtual bool SetNuclei(std::vector<int>& z, std::vector<int>& a);
+	virtual bool SetNuclei(std::vector<int>& z, std::vector<int>& a) = 0;
+	virtual void RunSystem() = 0;
+
 	void AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness);
 
 	/*Set sampling parameters*/
-	void SetRandomGenerator(TRandom3* gen);
+	virtual void SetRandomGenerator(TRandom3* gen);
 	inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); };
-	inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
 	inline void SetTheta1Range(double min, double max) { m_theta1Range = std::make_pair(min*deg2rad, max*deg2rad); };
 	inline void SetPhi1Range(double min, double max) { m_phi1Range = std::make_pair(min*deg2rad, max*deg2rad); };
 	inline void SetExcitationDistro(double mean, double sigma) { m_exDist = std::make_pair(mean, sigma); };
-	inline void SetDecay1Distribution(const std::string& file) { decay1dist.ReadDistributionFile(file); };
-	inline void SetDecay2Distribution (const std::string& file) { decay2dist.ReadDistributionFile(file); };
 
-	virtual void RunSystem();
-
-	inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
-	inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
-	inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
-	inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
-	inline const char* GetSystemEquation() const { return m_sys_equation.c_str(); };
-	inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
-	inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); };
-	inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
-	inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); };
+	inline const std::string& GetSystemEquation() const { return m_sys_equation; };
 
 protected:
-	virtual void LinkTarget();
-	virtual void SetSystemEquation();
+	virtual void LinkTarget() = 0;
+	virtual void SetSystemEquation() = 0;
 	
-	Reaction step1;
 	LayeredTarget target;
 
 	//Sampling information
 	std::pair<double, double> m_beamDist, m_theta1Range, m_phi1Range, m_exDist;
 	TRandom3* generator; //not owned by ReactionSystem
 
-	AngularDistribution decay1dist, decay2dist;
-
 	bool target_set_flag, gen_set_flag;
 	int rxnLayer;
 	std::string m_sys_equation;
diff --git a/include/ThreeStepSystem.h b/include/ThreeStepSystem.h
index e26f1a6..4d95356 100644
--- a/include/ThreeStepSystem.h
+++ b/include/ThreeStepSystem.h
@@ -2,6 +2,7 @@
 #define THREESTEPSYSTEM_H
 
 #include "ReactionSystem.h"
+#include "AngularDistribution.h"
 
 namespace Mask {
 
@@ -10,19 +11,35 @@ public:
 	ThreeStepSystem();
 	ThreeStepSystem(std::vector<int>& z, std::vector<int>& a);
 	~ThreeStepSystem();
-	bool SetNuclei(std::vector<int>& z, std::vector<int>& a);
-	void RunSystem();
+	bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
+	void RunSystem() override;
+	void SetRandomGenerator(TRandom3* gen) override;
 
+	inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
+	inline void SetDecay2Distribution(const std::string& filename) { decay2dist.ReadDistributionFile(filename); };
+
+	inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
+	inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
+	inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
+	inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
+	inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
 	inline const Nucleus& GetBreakup1() const { return step2.GetEjectile(); };
 	inline const Nucleus& GetBreakup2() const { return step2.GetResidual(); };
 	inline const Nucleus& GetBreakup3() const { return step3.GetEjectile(); };
 	inline const Nucleus& GetBreakup4() const { return step3.GetResidual(); };
 
-protected:
-	void LinkTarget();
-	void SetSystemEquation();
+	inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
+	inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); };
+	inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
+	inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); };
 
-	Reaction step2, step3;
+protected:
+	void LinkTarget() override;
+	void SetSystemEquation() override;
+
+	Reaction step1, step2, step3;
+
+	AngularDistribution decay1dist, decay2dist;
 
 };
 
diff --git a/include/TwoStepSystem.h b/include/TwoStepSystem.h
index 16ffc1b..44c927a 100644
--- a/include/TwoStepSystem.h
+++ b/include/TwoStepSystem.h
@@ -2,6 +2,7 @@
 #define TWOSTEPSYSTEM_H
 
 #include "ReactionSystem.h"
+#include "AngularDistribution.h"
 
 namespace Mask {
 
@@ -10,17 +11,31 @@ public:
 	TwoStepSystem();
 	TwoStepSystem(std::vector<int>& z, std::vector<int>& a);
 	~TwoStepSystem();
-	bool SetNuclei(std::vector<int>& z, std::vector<int>& a);
-	void RunSystem();
+	bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
+	void RunSystem() override;
 
+	void SetRandomGenerator(TRandom3* gen) override;
+
+	inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
+
+	inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); };
+	inline const Nucleus& GetTarget() const { return step1.GetTarget(); };
+	inline const Nucleus& GetProjectile() const { return step1.GetProjectile(); };
+	inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); };
+	inline const Nucleus& GetResidual() const { return step1.GetResidual(); };
 	inline const Nucleus& GetBreakup1() const { return step2.GetEjectile(); };
 	inline const Nucleus& GetBreakup2() const { return step2.GetResidual(); };
 
+	inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
+	inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
+
 private:
-	void LinkTarget();
-	void SetSystemEquation();
+	void LinkTarget() override;
+	void SetSystemEquation() override;
 	
-	Reaction step2;
+	Reaction step1, step2;
+
+	AngularDistribution decay1dist;
 
 };
 
diff --git a/input.txt b/input.txt
index 768161a..69ced3c 100644
--- a/input.txt
+++ b/input.txt
@@ -1,10 +1,10 @@
 ----------Data Information----------
-OutputFile: /data1/gwm17/TRIUMF_7Bed/simulation/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.root
+OutputFile: /data1/gwm17/mask_tests/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.root
 SaveTree: yes
 SavePlots: yes
 ----------Reaction Information----------
 ReactionType: 2
-Z	A (order is target, projectile, ejectile, break1, break3)
+Z	A (order is target, projectile, ejectile, break1, break3(if pure decay is target, ejectile))
 1	2
 4	7
 1	1
diff --git a/src/DecaySystem.cpp b/src/DecaySystem.cpp
new file mode 100644
index 0000000..b18165a
--- /dev/null
+++ b/src/DecaySystem.cpp
@@ -0,0 +1,75 @@
+#include "DecaySystem.h"
+
+namespace Mask {
+
+DecaySystem::DecaySystem() :
+	ReactionSystem()
+{
+}
+
+DecaySystem::DecaySystem(std::vector<int>& z, std::vector<int>& a) :
+	ReactionSystem()
+{
+	SetNuclei(z, a);
+}
+
+DecaySystem::~DecaySystem() {}
+
+void DecaySystem::SetRandomGenerator(TRandom3* gen) {
+	generator = gen;
+	decay1dist.AttachRandomNumberGenerator(gen);
+	gen_set_flag = true;
+}
+
+bool DecaySystem::SetNuclei(std::vector<int>& z, std::vector<int>& a) {
+	if(z.size() != a.size() || z.size() != 2) {
+		return false;
+	}
+
+	step1.SetNuclei(z[0], a[0], 0, 0, z[1], a[1]);
+	SetSystemEquation();
+	return true;
+}
+
+void DecaySystem::LinkTarget() {
+	step1.SetLayeredTarget(&target);
+
+	rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
+	if(rxnLayer != -1) {
+		step1.SetRxnLayer(rxnLayer);
+		target_set_flag = true;
+	} else {
+		throw ReactionLayerException();
+	}
+}
+
+void DecaySystem::SetSystemEquation() {
+	m_sys_equation = step1.GetTarget().GetIsotopicSymbol();
+	m_sys_equation += "-> ";
+	m_sys_equation += step1.GetEjectile().GetIsotopicSymbol();
+	m_sys_equation += "+";
+	m_sys_equation += step1.GetResidual().GetIsotopicSymbol();
+}
+
+void DecaySystem::RunSystem() {
+	if(!gen_set_flag) return;
+	
+	//Link up the target if it hasn't been done yet
+	if(!target_set_flag) {
+		LinkTarget();
+	}
+
+	if(step1.IsDecay()) {
+		double rxnTheta = std::acos(decay1dist.GetRandomCosTheta());
+		double rxnPhi = generator->Uniform(0, 2.0*M_PI);
+		step1.SetPolarRxnAngle(rxnTheta);
+		step1.SetAzimRxnAngle(rxnPhi);
+
+		step1.TurnOnResidualEloss();
+		step1.Calculate();
+	} else {
+		return;
+	}
+}
+
+}
\ No newline at end of file
diff --git a/src/Kinematics.cpp b/src/Kinematics.cpp
index b51de2a..b6d4997 100644
--- a/src/Kinematics.cpp
+++ b/src/Kinematics.cpp
@@ -43,9 +43,9 @@ bool Kinematics::LoadConfig(const char* filename) {
 	switch(m_rxn_type) {
 		case 0:
 		{
-			sys = new ReactionSystem();
+			sys = new DecaySystem();
 			m_rxn_type = ONESTEP_DECAY;
-			for(int i=0; i<3; i++) {
+			for(int i=0; i<2; i++) {
 				input>>z>>a;
 				avec.push_back(a);
 				zvec.push_back(z);
@@ -54,7 +54,7 @@ bool Kinematics::LoadConfig(const char* filename) {
 		}
 		case 1:
 		{
-			sys = new ReactionSystem();
+			sys = new OneStepSystem();
 			m_rxn_type = ONESTEP_RXN;
 			for(int i=0; i<3; i++) {
 				input>>z>>a;
@@ -109,6 +109,7 @@ bool Kinematics::LoadConfig(const char* filename) {
 		sys->AddTargetLayer(zvec, avec, svec, thickness);
 		input>>junk;
 	}
+	std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
 
 	double par1, par2;
 	std::string dfile1, dfile2;
@@ -119,7 +120,23 @@ bool Kinematics::LoadConfig(const char* filename) {
 	input>>junk>>par1>>junk>>par2;
 	sys->SetBeamDistro(par1, par2);
 	input>>junk>>par1;
-	sys->SetReactionThetaType(par1);
+	switch(m_rxn_type) {
+		case ONESTEP_RXN :
+		{
+			dynamic_cast<OneStepSystem*>(sys)->SetReactionThetaType(par1);
+			break;
+		}
+		case TWOSTEP:
+		{
+			dynamic_cast<TwoStepSystem*>(sys)->SetReactionThetaType(par1);
+			break;
+		}
+		case THREESTEP:
+		{
+			dynamic_cast<ThreeStepSystem*>(sys)->SetReactionThetaType(par1);
+			break;
+		}
+	}
 	input>>junk>>par1>>junk>>par2;
 	sys->SetTheta1Range(par1, par2);
 	input>>junk>>par1>>junk>>par2;
@@ -128,13 +145,35 @@ bool Kinematics::LoadConfig(const char* filename) {
 	sys->SetExcitationDistro(par1, par2);
 	input>>junk>>dfile1;
 	input>>junk>>dfile2;
-	sys->SetDecay1Distribution(dfile1);
-	sys->SetDecay2Distribution(dfile2);
+	switch(m_rxn_type) {
+		case ONESTEP_DECAY:
+		{
+			DecaySystem* this_sys = dynamic_cast<DecaySystem*>(sys);
+			this_sys->SetDecay1Distribution(dfile1);
+			std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<std::endl;
+			std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<std::endl;
+			break;
+		}
+		case TWOSTEP:
+		{
+			TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(sys);
+			this_sys->SetDecay1Distribution(dfile1);
+			std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<std::endl;
+			std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<std::endl;
+			break;
+		}
+		case THREESTEP:
+		{
+			ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(sys);
+			this_sys->SetDecay1Distribution(dfile1);
+			this_sys->SetDecay2Distribution(dfile2);
+			std::cout<<"Decay1 angular momentum: "<<this_sys->GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<<this_sys->GetDecay2AngularMomentum()<<std::endl;
+			std::cout<<"Decay1 total branching ratio: "<<this_sys->GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<<this_sys->GetDecay2BranchingRatio()<<std::endl;
+			break;
+		}
+	}
 	sys->SetRandomGenerator(global_generator);
 
-	std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
-	std::cout<<"Decay1 angular momentum: "<<sys->GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<<sys->GetDecay2AngularMomentum()<<std::endl;
-	std::cout<<"Decay1 total branching ratio: "<<sys->GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<<sys->GetDecay2BranchingRatio()<<std::endl;
 	std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
 
 	return true;
@@ -186,6 +225,10 @@ void Kinematics::Run() {
 }
 
 void Kinematics::RunOneStepRxn() {
+	OneStepSystem* this_sys = dynamic_cast<OneStepSystem*>(sys);
+	if(this_sys == nullptr) {
+		return;
+	}
 
 	
 	TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
@@ -214,20 +257,20 @@ void Kinematics::RunOneStepRxn() {
 			std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
 		}
 
-		sys->RunSystem();
+		this_sys->RunSystem();
 		
 		if(save_tree_flag) {
-			targ = ConvertNucleus(sys->GetTarget());
-			proj = ConvertNucleus(sys->GetProjectile());
-			eject = ConvertNucleus(sys->GetEjectile());
-			residual = ConvertNucleus(sys->GetResidual());
+			targ = ConvertNucleus(this_sys->GetTarget());
+			proj = ConvertNucleus(this_sys->GetProjectile());
+			eject = ConvertNucleus(this_sys->GetEjectile());
+			residual = ConvertNucleus(this_sys->GetResidual());
 			tree->Fill();
 		}
 		if(do_plotter_flag) {
-			plotman.FillData(sys->GetTarget());
-			plotman.FillData(sys->GetProjectile());
-			plotman.FillData(sys->GetEjectile());
-			plotman.FillData(sys->GetResidual());
+			plotman.FillData(this_sys->GetTarget());
+			plotman.FillData(this_sys->GetProjectile());
+			plotman.FillData(this_sys->GetEjectile());
+			plotman.FillData(this_sys->GetResidual());
 		}
 
 	}
@@ -244,6 +287,10 @@ void Kinematics::RunOneStepRxn() {
 }
 
 void Kinematics::RunOneStepDecay() {
+	DecaySystem* this_sys = dynamic_cast<DecaySystem*>(sys);
+	if(this_sys == nullptr) {
+		return;
+	}
 
 	TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
 	TTree* tree;
@@ -267,17 +314,17 @@ void Kinematics::RunOneStepDecay() {
 			std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
 		}
 
-		sys->RunSystem();
+		this_sys->RunSystem();
 		if(save_tree_flag) {
-			targ = ConvertNucleus(sys->GetTarget());
-			eject = ConvertNucleus(sys->GetEjectile());
-			residual = ConvertNucleus(sys->GetResidual());
+			targ = ConvertNucleus(this_sys->GetTarget());
+			eject = ConvertNucleus(this_sys->GetEjectile());
+			residual = ConvertNucleus(this_sys->GetResidual());
 			tree->Fill();
 		}
 		if(do_plotter_flag) {
-			plotman.FillData(sys->GetTarget());
-			plotman.FillData(sys->GetEjectile());
-			plotman.FillData(sys->GetResidual());
+			plotman.FillData(this_sys->GetTarget());
+			plotman.FillData(this_sys->GetEjectile());
+			plotman.FillData(this_sys->GetResidual());
 		}
 
 	}
diff --git a/src/OneStepSystem.cpp b/src/OneStepSystem.cpp
new file mode 100644
index 0000000..d9b2b69
--- /dev/null
+++ b/src/OneStepSystem.cpp
@@ -0,0 +1,77 @@
+#include "OneStepSystem.h"
+
+namespace Mask {
+
+OneStepSystem::OneStepSystem() :
+	ReactionSystem()
+{
+}
+
+OneStepSystem::OneStepSystem(std::vector<int>& z, std::vector<int>& a) :
+	ReactionSystem()
+{
+	SetNuclei(z, a);
+}
+
+OneStepSystem::~OneStepSystem() {}
+
+bool OneStepSystem::SetNuclei(std::vector<int>& z, std::vector<int>& a) {
+	if(z.size() != a.size() || z.size() != 3) {
+		return false;
+	}
+
+	step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
+	SetSystemEquation();
+	return true;
+}
+
+void OneStepSystem::LinkTarget() {
+	step1.SetLayeredTarget(&target);
+
+	rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
+	if(rxnLayer != -1) {
+		step1.SetRxnLayer(rxnLayer);
+		target_set_flag = true;
+	} else {
+		throw ReactionLayerException();
+	}
+}
+
+void OneStepSystem::SetSystemEquation() {
+	m_sys_equation = step1.GetTarget().GetIsotopicSymbol();
+	m_sys_equation += "(";
+	m_sys_equation += step1.GetProjectile().GetIsotopicSymbol();
+	m_sys_equation += ", ";
+	m_sys_equation += step1.GetEjectile().GetIsotopicSymbol();
+	m_sys_equation += ")";
+	m_sys_equation += step1.GetResidual().GetIsotopicSymbol();
+}
+
+void OneStepSystem::RunSystem() {
+	if(!gen_set_flag) return;
+	
+	//Link up the target if it hasn't been done yet
+	if(!target_set_flag) {
+		LinkTarget();
+	}
+
+	if(!step1.IsDecay()) {
+		//Sample parameters
+		double bke = generator->Gaus(m_beamDist.first, m_beamDist.second);
+		double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second)));
+		double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second);
+		double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
+	
+		step1.SetBeamKE(bke);
+		step1.SetPolarRxnAngle(rxnTheta);
+		step1.SetAzimRxnAngle(rxnPhi);
+		step1.SetExcitation(residEx);
+	
+		step1.TurnOnResidualEloss();
+		step1.Calculate();
+	} else {
+		return;
+	}
+}
+
+}
\ No newline at end of file
diff --git a/src/ReactionSystem.cpp b/src/ReactionSystem.cpp
index 3cf53d1..4034dff 100644
--- a/src/ReactionSystem.cpp
+++ b/src/ReactionSystem.cpp
@@ -9,30 +9,10 @@ ReactionSystem::ReactionSystem() :
 {
 }
 
-ReactionSystem::ReactionSystem(std::vector<int>& z, std::vector<int>& a) :
-	m_beamDist(0,0), m_theta1Range(0,0), m_phi1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("")
-{
-	SetNuclei(z, a);
-}
-
-ReactionSystem::~ReactionSystem() {
-
-}
-
-bool ReactionSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
-	if(z.size() != a.size() || z.size() < 3) {
-		return false;
-	}
-
-	step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
-	SetSystemEquation();
-	return true;
-}
+ReactionSystem::~ReactionSystem() {}
 
 void ReactionSystem::SetRandomGenerator(TRandom3* gen) {
 	generator = gen;
-	decay1dist.AttachRandomNumberGenerator(gen);
-	decay2dist.AttachRandomNumberGenerator(gen);
 	gen_set_flag = true;
 }
 
@@ -40,66 +20,4 @@ void ReactionSystem::AddTargetLayer(std::vector<int>& zt, std::vector<int>& at,
 	target.AddLayer(zt, at, stoich, thickness);
 }
 
-void ReactionSystem::LinkTarget() {
-	step1.SetLayeredTarget(&target);
-
-	rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
-	if(rxnLayer != -1) {
-		step1.SetRxnLayer(rxnLayer);
-		target_set_flag = true;
-	} else {
-		throw ReactionLayerException();
-	}
-}
-
-void ReactionSystem::SetSystemEquation()  {
-	m_sys_equation = step1.GetTarget().GetIsotopicSymbol();
-	m_sys_equation += "(";
-	m_sys_equation += step1.GetProjectile().GetIsotopicSymbol();
-	m_sys_equation += ", ";
-	m_sys_equation += step1.GetEjectile().GetIsotopicSymbol();
-	m_sys_equation += ")";
-	m_sys_equation += step1.GetResidual().GetIsotopicSymbol();
-}
-
-void ReactionSystem::RunSystem() {
-	if(!gen_set_flag) return;
-	
-	//Link up the target if it hasn't been done yet
-	if(!target_set_flag) {
-		LinkTarget();
-	}
-
-	if(step1.IsDecay()) {
-		double decaycostheta = decay1dist.GetRandomCosTheta();
-		if(decaycostheta == -10.0) {
-			step1.ResetEjectile();
-			step1.ResetResidual();
-			return;
-		}
-		double rxnTheta = std::acos(decay1dist.GetRandomCosTheta());
-		double rxnPhi = generator->Uniform(0, 2.0*M_PI);
-		step1.SetPolarRxnAngle(rxnTheta);
-		step1.SetAzimRxnAngle(rxnPhi);
-
-		step1.TurnOnResidualEloss();
-		step1.Calculate();
-	} else {
-		//Sample parameters
-		double bke = generator->Gaus(m_beamDist.first, m_beamDist.second);
-		double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second)));
-		double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second);
-		double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
-	
-		step1.SetBeamKE(bke);
-		step1.SetPolarRxnAngle(rxnTheta);
-		step1.SetAzimRxnAngle(rxnPhi);
-		step1.SetExcitation(residEx);
-	
-		step1.TurnOnResidualEloss();
-		step1.Calculate();
-	}
-
-}
-
 };
\ No newline at end of file
diff --git a/src/ThreeStepSystem.cpp b/src/ThreeStepSystem.cpp
index 8252109..88c5256 100644
--- a/src/ThreeStepSystem.cpp
+++ b/src/ThreeStepSystem.cpp
@@ -18,6 +18,13 @@ ThreeStepSystem::~ThreeStepSystem() {
 
 }
 
+void ThreeStepSystem::SetRandomGenerator(TRandom3* gen) {
+	generator = gen;
+	decay1dist.AttachRandomNumberGenerator(gen);
+	decay2dist.AttachRandomNumberGenerator(gen);
+	gen_set_flag = true;
+}
+
 bool ThreeStepSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
 	if(z.size() != a.size() || z.size() < 5) {
 		return false;
diff --git a/src/TwoStepSystem.cpp b/src/TwoStepSystem.cpp
index e7c17d2..c4f69d9 100644
--- a/src/TwoStepSystem.cpp
+++ b/src/TwoStepSystem.cpp
@@ -18,8 +18,14 @@ TwoStepSystem::~TwoStepSystem() {
 
 }
 
+void TwoStepSystem::SetRandomGenerator(TRandom3* gen) {
+	generator = gen;
+	decay1dist.AttachRandomNumberGenerator(gen);
+	gen_set_flag = true;
+}
+
 bool TwoStepSystem::SetNuclei(std::vector<int>&z, std::vector<int>& a) {
-	if(z.size() != a.size() || z.size() < 4) {
+	if(z.size() != a.size() || z.size() != 4) {
 		return false;
 	}
 	int zr = z[0] + z[1] - z[2];