diff --git a/CMakeLists.txt b/CMakeLists.txt
index 289801a..55d5672 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,13 @@
 cmake_minimum_required(VERSION 3.16)
 
+
+if(NOT CMAKE_BUILD_TYPE STREQUAL "Debug")
+    set(CMAKE_BUILD_TYPE "Release")
+    message("Building release")
+else()
+    message("Building debug")
+endif()
+
 project(Mask)
 
 set(MASK_BINARY_DIR ${CMAKE_CURRENT_SOURCE_DIR}/bin)
diff --git a/README.md b/README.md
index ec05426..1842173 100644
--- a/README.md
+++ b/README.md
@@ -1,55 +1,58 @@
 # MASK: Monte cArlo Simulation of Kinematics
 MASK is a Monte Carlo simulation of reaction kinematics for use detector systems at Florida State University.
-MASK is capable of simulating multi-step kinematic reaction-decay sequences, storing data in a lightweight binary format, after which the kinematic data can be fed to a detector geometry for efficiency testing. Currently geometries for ANASEN and SABRE are included in the code.
+MASK is capable of simulating multi-step kinematic reaction-decay sequences, storing data in a lightweight binary format, after which the kinematic data can be fed to a detector geometry for efficiency testing. Currently geometries for ANASEN and SABRE are included in the code. Mask depends on the CERN ROOT analysis framework. All data is stored in the ROOT format using a custom dictionary.
 
 ## Building MASK
-Dowload the repository from github. The code is to be built using Premake5, an open source and free project building software. To build on Linux and MacOSX run
-`premake5 gmake2`
-to generate the makefiles. Then execute
-`make`
-to build the program. To build on Windows run
-`premake5 vs20xx`
-where xx should be replaced with your version of Visual Studio, and then build as a normal Visual Studio project. For more documentation see the Premake wiki. By default `make` runs in Debug mode. For release mode use `make config=release`
+To clone MASK and all submodules use `git clone --recursive https://github.com/gwm17/Kinematics.git`.
+Mask is built using the CMake build tool. On most systems this can be done using the following commands inside the Mask repository:
 
-### Building RootPlot
-One of the executables, RootPlot, requires linking against external libraries from the ROOT cern analysis package. Path to the necessary header files and libraries must be set by the user on a machine-by-machine basis in the premake5.lua file. Set the ROOTIncludepath and ROOTLibpath to match your install. An easy way to check for the paths on a unix system is through the root-config tool. root-config --cflags has the include path (the path after -I) and root-config --glibs has the lib path (after -L). If you do not have ROOT installed, you will not be able to compile RootPlot, and you will not be able to use the generic `make` command, which will try to build all executables. Instead use `make Mask` or `make DetectEff` (or `make config=release Mask`, etc).
+- `mkdir build`
+- `cd build`
+- `cmake ..`
+- `make`
 
-## Running MASK
+By default Mask builds for release. To build for debug replace `cmake ..` with `cmake -DCMAKE_BUILD_TYPE=Debug ..`. Mask uses CMake to find the installed ROOT libraries and headers.
+
+## Using the kinematics simulation
 By default MASK is capable of simulating reactions of up to three steps. Here is a brief outline of each type:
 
-0. A reaction of type 0 is a pure decay. It is assumed isotropic by default; any other case will require the modification of the code.
-1. A reaction of type 1 is a pure reaction. It can incorporate all of the input file sampling parameters.
-2. A reaction of type 2 is a reaction followed by a subsequent decay of the residual nucleus. Again, all sampling is allowed.
-3. A reaction of type 3 is a reaction followed by a subsequent decay of the residual, followed by a decay of one of the products. Again, all sampling is allowed
+0. A pure decay involves a "target" decaying into an ejectile and residual. It is assumed isotropic by default; any other case will require the modification of the code.
+1. A single step reaction involves a target being hit by a projectile and emitting an ejectile and a residual. It can incorporate all of the input file sampling parameters.
+2. A two step reaction is a single step reaction where the residual subsequently decays into two daughters (called breakups). Again, all sampling is allowed.
+3. A three step reaction is a two step reaction where one of the breakup particles subsequently decays into two more breakup particles. Again, all sampling is allowed
 
-For decays, a specific angular distribution can be given as input as a text file with values of coefficiencts of a Legendre polynomial series. Examples can be found in the `./etc` directory, including an isotropic case. It is assumed that the decays in the center-of-mass frame are isotropic in phi (i.e. m=0). Decay1 corresponds to the first decay, if there are multiple steps, Decay2 to the second. If there are no decays, these parameters are not used (or if only one decay, Decay2_AngularMomentum is not used). The input file requires that the user include target information, which will be used to calculate energy loss for all of the reactants and reaction products. The target can contain layers, and each layer can be composed of a compound of elements with a given stoichiometry. If the user wishes to not include energy loss in the kinematics, simply give all target layers a thickness of 0. Note that more layers and more thickness = more time spent calculating energy loss. These energy loss methods are only applicable for solid targets, and should not be applied to gas or liquid targets. Energy loss calculations have a stated uncertainty of approximately five percent.
+For decays, a specific angular distribution can be given as input as a text file with values of coefficiencts of a Legendre polynomial series. Examples can be found in the `./etc` directory, including an isotropic case. It is assumed that the decays in the center-of-mass frame are isotropic in phi (i.e. m=0). Decay1 corresponds to the first decay, if there are multiple steps, Decay2 to the second. If there are no decays, these parameters are not used (or if only one decay, Decay2_AngularMomentum is not used).
 
-To choose which detector scheme is run, modify the main function in DetectorEfficiency.cpp. The included geometries also have options to do an internal geometry consistency check and print out coordinates for drawing the detector arrays.
+The input file requires that the user include target information, which will be used to calculate energy loss for all of the reactants and reaction products. The target can contain layers, and each layer can be composed of a compound of elements with a given stoichiometry. If the user wishes to not include energy loss in the kinematics, simply give all target layers a thickness of 0. Note that more layers and more thickness = more time spent calculating energy loss. These energy loss methods are only applicable for solid targets, and should not be applied to gas or liquid targets. Energy loss calculations have a stated uncertainty of approximately five percent. The energy loss library used is called `catima` and can be found [here](https://github.com/gwm17/catima).
 
-To run MASK simply do the following from the MASK directory:
+To run MASK simply do the following from the MASK repository:
 
-`./bin/Mask input.txt`
+`./bin/MaskApp input.txt`
 
 Input.txt can be replaced by any text file with the correct format.
 
+## Using the detector geometry simulation
+Detector geometry is encoded using ROOT math libraries in the `src/Detectors` folder. Two different detector geometries are already present: SPS-SABRE and ANASEN. To add a new geometry, follow the guidelines outlined by each of these cases.
+
+To choose which detector scheme is run, modify the main function in `src/Detectors/main.cpp`. The included geometries also have options to do an internal geometry consistency check and print out coordinates for drawing the detector arrays, which can be useful for testing.
+
+
+
 To run DetEff use the format
 
 `./bin/DetEff <kinematics_datafile> <new_detection_datafile> <new_detection_statsfile>`
 
 where the detection datafile contains all of the kinematics data as well as information about which particles are detected (this is in the mask file format) and the statsfile is a text file containing efficiency statistics.
 
-RootPlot is run as
+## Data visualization
+All data is saved as ROOT trees of std::vectors of Mask::Nucleus classes. To enable this, a ROOT dictionary is generated and linked into a shared library found in the `lib` directory of the repository. This allows the user to link to the shared library for accessing and analyzing the data generated by MASK.
+
+Mask also provides a default visualization tool called RootPlot. RootPlot is run as
 
 `./bin/RootPlot <datafile> <outputfile>`
 
-where the datafile can be either the datafile from Mask or the datafile from DetEff. The outputfile is saved in the ROOT file format.
-
-PyPlotter is run as
-
-`./bin/PyPlotter <datafile> <outputfile>`
-
-where again datafile can come from either Mask or DetEff. The outputfile is saved in the Python pickle file format. Pickled files can be reopened using the PyPlotViewer script run as
-`./bin/PyPlotViewer <picklefile>`
+where the datafile can be either the datafile from Mask or the datafile from DetEff.s
 
 ## Requirements
-MASK, the kinematics simulation, and DetEff, the detection efficiency simulation, require no external dependancies. The RootPlot plotting tool requires the ROOT cern data analysis package, and the PyPlotter tool requires Python3 and the matplotlib and numpy packages. 
\ No newline at end of file
+- Requires CMake > 3.16
+- Requires ROOT > 6.16
\ No newline at end of file
diff --git a/input.txt b/input.txt
index 215185c..5ca223b 100644
--- a/input.txt
+++ b/input.txt
@@ -1,25 +1,26 @@
 ----------Data Information----------
-OutputFile: /media/data/gwm17/mask_tests/10B3Hea_test.root
+OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
 ----------Reaction Information----------
-ReactionType: OneStepRxn
-Z	A (order is target, projectile, ejectile, break1, break3(if pure decay is target, ejectile))
-5	10
-2	3
-2	4
+begin_nuclei(Z,A)
+	5	10
+	2	3
+	2	4
+	2	4
+end_nuclei(Z,A)
 ----------Target Information----------
 NumberOfLayers: 1
 begin_layer
-	Thickness(ug/cm^2): 50
+	Thickness(ug/cm^2): 74
 	begin_elements (Z, A, Stoich.)
 		element	5 10 1
 	end_elements
 end_layer
 ----------Sampling Information----------
-NumberOfSamples: 1000
+NumberOfSamples: 1000000
 BeamMeanEnergy(MeV): 24.0 BeamEnergySigma(MeV): 0.0
-EjectileThetaType(0=Lab,1=CM): 1
-EjectileThetaMin(deg): 0.0 EjectileThetaMax(deg): 180.0
-EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 360.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
diff --git a/src/Mask/AngularDistribution.cpp b/src/Mask/AngularDistribution.cpp
index 54ec7a9..0884e80 100644
--- a/src/Mask/AngularDistribution.cpp
+++ b/src/Mask/AngularDistribution.cpp
@@ -72,6 +72,9 @@ namespace Mask {
 		//Total branching ratio
 		m_branchingRatio = m_constants[0]*2.0;
 		m_L = l;
+
+		std::cout<<"Angular distribution from "<<file<<" branching ratio: "<<m_branchingRatio<<std::endl;
+		std::cout<<"Angular distribution from "<<file<<" L: "<<m_L<<std::endl;
 	
 		//Renormalize distribution such that total prob is 1.0.
 		//Test branching ratio to see if we "make" a decay particle,
diff --git a/src/Mask/DecaySystem.h b/src/Mask/DecaySystem.h
index 5a2aaab..c0a1122 100644
--- a/src/Mask/DecaySystem.h
+++ b/src/Mask/DecaySystem.h
@@ -16,7 +16,7 @@ namespace Mask {
 		void RunSystem() override;
 		std::vector<Nucleus>*GetNuclei() override;
 	
-		void SetDecay1Distribution(const std::string& filename) { m_step1Distribution.ReadDistributionFile(filename); }
+		virtual void SetDecay1Distribution(const std::string& filename) override { m_step1Distribution.ReadDistributionFile(filename); }
 	
 		int GetDecay1AngularMomentum() { return m_step1Distribution.GetL(); }
 		double GetDecay1BranchingRatio() { return m_step1Distribution.GetBranchingRatio(); }
diff --git a/src/Mask/MaskApp.cpp b/src/Mask/MaskApp.cpp
index 490b3e4..07a92e8 100644
--- a/src/Mask/MaskApp.cpp
+++ b/src/Mask/MaskApp.cpp
@@ -36,62 +36,29 @@ namespace Mask {
 		int z, a, s;
 		getline(input, junk);
 		getline(input, junk);
-		input>>junk>>junk;
-		m_rxnType = GetRxnTypeFromString(junk);
-		getline(input, junk);
-		getline(input, junk);
-		switch(m_rxnType) 
+
+		while(input>>junk)
 		{
-			case RxnType::PureDecay:
-			{
-				m_system = new DecaySystem();
-				for(int i=0; i<2; i++)
-				{
-					input>>z>>a;
-					avec.push_back(a);
-					zvec.push_back(z);
-				}
+			if(junk == "begin_nuclei(Z,A)")
+				continue;
+			else if (junk == "end_nuclei(Z,A)")
 				break;
-			}
-			case RxnType::OneStepRxn:
+			else
 			{
-				m_system = new OneStepSystem();
-				for(int i=0; i<3; i++)
-				{
-					input>>z>>a;
-					avec.push_back(a);
-					zvec.push_back(z);
-				}
-				std::cout<<"here"<<std::endl;
-				break;
+				z = std::stoi(junk);
+				input>>a;
+				zvec.push_back(z);
+				avec.push_back(a);
 			}
-			case RxnType::TwoStepRxn:
-			{
-				m_system = new TwoStepSystem();
-				for(int i=0; i<4; i++)
-				{
-					input>>z>>a;
-					avec.push_back(a);
-					zvec.push_back(z);
-				}
-				break;
-			}
-			case RxnType::ThreeStepRxn:
-			{
-				m_system = new ThreeStepSystem();
-				for(int i=0; i<5; i++)
-				{
-					input>>z>>a;
-					avec.push_back(a);
-					zvec.push_back(z);
-				}
-				break;
-			}
-			default:
-				return false;
 		}
-		m_system->SetNuclei(zvec, avec);
-	
+
+		m_system = CreateSystem(zvec, avec);
+		if(m_system == nullptr)
+		{
+			std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
+			return false;
+		}
+
 		int nlayers;
 		double thickness;
 		getline(input, junk);
@@ -120,33 +87,16 @@ namespace Mask {
 	
 		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>>par1;
-		switch(m_rxnType) 
-		{
-			case RxnType::PureDecay : break;
-			case RxnType::None : break;
-			case RxnType::OneStepRxn :
-			{
-				dynamic_cast<OneStepSystem*>(m_system)->SetReactionThetaType(par1);
-				break;
-			}
-			case RxnType::TwoStepRxn :
-			{
-				dynamic_cast<TwoStepSystem*>(m_system)->SetReactionThetaType(par1);
-				break;
-			}
-			case RxnType::ThreeStepRxn :
-			{
-				dynamic_cast<ThreeStepSystem*>(m_system)->SetReactionThetaType(par1);
-				break;
-			}
-		}
+		input>>junk>>thetaTypeString;
+		m_system->SetReactionThetaType(StringToRxnThetaType(thetaTypeString));
+		
 		input>>junk>>par1>>junk>>par2;
 		m_system->SetTheta1Range(par1, par2);
 		input>>junk>>par1>>junk>>par2;
@@ -155,36 +105,8 @@ namespace Mask {
 		m_system->SetExcitationDistro(par1, par2);
 		input>>junk>>dfile1;
 		input>>junk>>dfile2;
-		switch(m_rxnType) 
-		{
-			case RxnType::PureDecay :
-			{
-				DecaySystem* this_sys = dynamic_cast<DecaySystem*>(m_system);
-				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 RxnType::None : break;
-			case RxnType::OneStepRxn : break;
-			case RxnType::TwoStepRxn :
-			{
-				TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(m_system);
-				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 RxnType::ThreeStepRxn :
-			{
-				ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(m_system);
-				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;
-			}
-		}
+		m_system->SetDecay1Distribution(dfile1);
+		m_system->SetDecay2Distribution(dfile2);
 	
 		std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
 	
diff --git a/src/Mask/Nucleus.h b/src/Mask/Nucleus.h
index d5b4b47..7baef5c 100644
--- a/src/Mask/Nucleus.h
+++ b/src/Mask/Nucleus.h
@@ -37,8 +37,8 @@ namespace Mask {
 			return vec4.M() - groundStateMass;
 		}
 
-		int Z = 0; 
-		int A = 0;
+		uint32_t Z = 0; 
+		uint32_t A = 0;
 		double groundStateMass = 0.0;
 		std::string isotopicSymbol = "";
 		double thetaCM = 0.0;
diff --git a/src/Mask/OneStepSystem.h b/src/Mask/OneStepSystem.h
index 6e7a352..fe2248e 100644
--- a/src/Mask/OneStepSystem.h
+++ b/src/Mask/OneStepSystem.h
@@ -15,7 +15,7 @@ namespace Mask {
 		void RunSystem() override;
 		std::vector<Nucleus>* GetNuclei() override;
 	
-		void SetReactionThetaType(int type) { m_step1.SetEjectileThetaType(type); };
+		virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
 	
 	private:
 		void LinkTarget() override;
diff --git a/src/Mask/Reaction.cpp b/src/Mask/Reaction.cpp
index e77e144..fbc1d07 100644
--- a/src/Mask/Reaction.cpp
+++ b/src/Mask/Reaction.cpp
@@ -15,13 +15,13 @@ namespace Mask {
 
 	Reaction::Reaction() :
 		m_target(nullptr), m_projectile(nullptr), m_ejectile(nullptr), m_residual(nullptr), m_layeredTarget(nullptr), 
-		m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_ejectThetaType(s_lab), m_isInit(false), m_isResidEloss(false)
+		m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_ejectThetaType(RxnThetaType::None), m_isInit(false), m_isResidEloss(false)
 	{
 	}
 	
 	Reaction::Reaction(Nucleus* target, Nucleus* projectile, Nucleus* ejectile, Nucleus* residual) :
 		m_target(nullptr), m_projectile(nullptr), m_ejectile(nullptr), m_residual(nullptr),
-		m_layeredTarget(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_ejectThetaType(s_lab), m_isResidEloss(false)
+		m_layeredTarget(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_ejectThetaType(RxnThetaType::None), m_isResidEloss(false)
 	{
 		BindNuclei(target, projectile, ejectile, residual);
 	}
@@ -70,13 +70,11 @@ namespace Mask {
 		m_bke = bke - m_layeredTarget->GetProjectileEnergyLoss(m_projectile->Z, m_projectile->A, bke, m_rxnLayer, 0);
 	}
 	
-	void Reaction::SetEjectileThetaType(int type)
+	void Reaction::SetEjectileThetaType(RxnThetaType type)
 	{
 		if(m_isDecay)
 			return;
-		if(type != s_centerOfMass && type != s_lab)
-			return;
-	
+
 		m_ejectThetaType = type;
 	}
 	
@@ -180,8 +178,9 @@ namespace Mask {
 	{
 		switch(m_ejectThetaType)
 		{
-			case s_centerOfMass: CalculateReactionThetaCM(); break;
-			case s_lab: CalculateReactionThetaLab(); break;
+			case RxnThetaType::CenterOfMass: CalculateReactionThetaCM(); break;
+			case RxnThetaType::Lab: CalculateReactionThetaLab(); break;
+			case RxnThetaType::None: CalculateReactionThetaCM(); break; //default behavior
 		}
 	}
 	
diff --git a/src/Mask/Reaction.h b/src/Mask/Reaction.h
index a1c166c..7923139 100644
--- a/src/Mask/Reaction.h
+++ b/src/Mask/Reaction.h
@@ -11,6 +11,7 @@
 
 #include "Nucleus.h"
 #include "LayeredTarget.h"
+#include "RxnType.h"
 
 namespace Mask {
 
@@ -24,7 +25,7 @@ namespace Mask {
 
 		void BindNuclei(Nucleus* target, Nucleus* projectile, Nucleus* ejectile, Nucleus* residual);
 		void SetBeamKE(double bke);
-		void SetEjectileThetaType(int type);
+		void SetEjectileThetaType(RxnThetaType type);
 	
 		void SetLayeredTarget(LayeredTarget* targ) { m_layeredTarget = targ; };
 
@@ -61,12 +62,9 @@ namespace Mask {
 		double m_bke, m_theta, m_phi, m_ex;
 	
 		int m_rxnLayer;
-		int m_ejectThetaType; 
+		RxnThetaType m_ejectThetaType; 
 	
 		bool m_isDecay, m_isInit, m_isResidEloss;
-	
-		static constexpr int s_lab = 0;
-		static constexpr int s_centerOfMass = 1;
 	};
 
 }
diff --git a/src/Mask/ReactionSystem.cpp b/src/Mask/ReactionSystem.cpp
index dd18de2..f8e9cb9 100644
--- a/src/Mask/ReactionSystem.cpp
+++ b/src/Mask/ReactionSystem.cpp
@@ -1,6 +1,11 @@
 #include "ReactionSystem.h"
+#include "RxnType.h"
 #include "KinematicsExceptions.h"
 #include "LegendrePoly.h"
+#include "DecaySystem.h"
+#include "OneStepSystem.h"
+#include "TwoStepSystem.h"
+#include "ThreeStepSystem.h"
 
 namespace Mask {
 
@@ -23,4 +28,23 @@ namespace Mask {
 		m_target.AddLayer(zt, at, stoich, thickness);
 	}
 
+	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;
+			return nullptr;
+		}
+
+		switch(a.size())
+		{
+			case RxnSize::DecaySize: return new DecaySystem(z, a);
+			case RxnSize::OneStepSize: return new OneStepSystem(z, a);
+			case RxnSize::TwoStepSize: return new TwoStepSystem(z, a);
+			case RxnSize::ThreeStepSize: return new ThreeStepSystem(z, a);
+		}
+
+		return nullptr;
+	}
+
 }
\ No newline at end of file
diff --git a/src/Mask/ReactionSystem.h b/src/Mask/ReactionSystem.h
index 43a6aca..02f9fbf 100644
--- a/src/Mask/ReactionSystem.h
+++ b/src/Mask/ReactionSystem.h
@@ -11,6 +11,7 @@
 
 #include "Reaction.h"
 #include "KinematicsExceptions.h"
+#include "RxnType.h"
 #include <vector>
 #include <random>
 
@@ -25,6 +26,9 @@ namespace Mask {
 		virtual bool SetNuclei(const std::vector<int>& z, const std::vector<int>& a) = 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);
 	
@@ -56,7 +60,7 @@ namespace Mask {
 				delete m_exDist;
 			m_exDist = new std::normal_distribution<double>(mean, sigma); 
 		}
-	
+
 		const std::string& GetSystemEquation() const { return m_sysEquation; }
 	
 	protected:
@@ -76,6 +80,7 @@ namespace Mask {
 		static constexpr double s_deg2rad = M_PI/180.0;
 	};
 
+	ReactionSystem* CreateSystem(const std::vector<int>& z, const std::vector<int>& a);
 }
 
 #endif
\ No newline at end of file
diff --git a/src/Mask/RxnType.h b/src/Mask/RxnType.h
index bd2f8bc..9b4f841 100644
--- a/src/Mask/RxnType.h
+++ b/src/Mask/RxnType.h
@@ -14,6 +14,21 @@ namespace Mask {
 		None=4
 	};
 
+	enum class RxnThetaType
+	{
+		CenterOfMass,
+		Lab,
+		None
+	};
+
+	enum RxnSize
+	{
+		DecaySize = 2,
+		OneStepSize = 3,
+		TwoStepSize = 4,
+		ThreeStepSize = 5
+	};
+
 	static RxnType GetRxnTypeFromString(const std::string& type_str)
 	{
 		if (type_str == "PureDecay")
@@ -52,6 +67,21 @@ namespace Mask {
 		return static_cast<uint32_t>(type);
 	}
 
+	static RxnThetaType StringToRxnThetaType(const std::string& type)
+	{
+		if(type == "CenterOfMass")
+		{
+			return RxnThetaType::CenterOfMass;
+		}
+		else if(type == "Lab")
+		{
+			return RxnThetaType::Lab;
+		}
+		else
+		{
+			return RxnThetaType::None;
+		}
+	}
 }
 
 #endif
\ No newline at end of file
diff --git a/src/Mask/ThreeStepSystem.h b/src/Mask/ThreeStepSystem.h
index 9e800a2..a79017e 100644
--- a/src/Mask/ThreeStepSystem.h
+++ b/src/Mask/ThreeStepSystem.h
@@ -15,10 +15,10 @@ namespace Mask {
 		void RunSystem() override;
 		std::vector<Nucleus>* GetNuclei() override;
 	
-		inline void SetDecay1Distribution(const std::string& filename) { m_step2Distribution.ReadDistributionFile(filename); };
-		inline void SetDecay2Distribution(const std::string& filename) { m_step3Distribution.ReadDistributionFile(filename); };
+		virtual void SetDecay1Distribution(const std::string& filename) override { m_step2Distribution.ReadDistributionFile(filename); };
+		virtual void SetDecay2Distribution(const std::string& filename) override { m_step3Distribution.ReadDistributionFile(filename); };
 	
-		void SetReactionThetaType(int type) { m_step1.SetEjectileThetaType(type); };
+		virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
 	
 		int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); };
 		int GetDecay2AngularMomentum(){ return m_step3Distribution.GetL(); };
diff --git a/src/Mask/TwoStepSystem.h b/src/Mask/TwoStepSystem.h
index 55623b4..9b1249b 100644
--- a/src/Mask/TwoStepSystem.h
+++ b/src/Mask/TwoStepSystem.h
@@ -16,9 +16,9 @@ namespace Mask {
 		void RunSystem() override;
 		std::vector<Nucleus>* GetNuclei() override;
 	
-		void SetDecay1Distribution(const std::string& filename) { m_step2Distribution.ReadDistributionFile(filename); };
+		virtual void SetDecay1Distribution(const std::string& filename) override { m_step2Distribution.ReadDistributionFile(filename); };
 	
-		void SetReactionThetaType(int type) { m_step1.SetEjectileThetaType(type); };
+		virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); };
 
 		int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); };
 		double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); };
diff --git a/src/Plotters/RootPlotter.cpp b/src/Plotters/RootPlotter.cpp
index 6bb7310..476f934 100644
--- a/src/Plotters/RootPlotter.cpp
+++ b/src/Plotters/RootPlotter.cpp
@@ -34,13 +34,20 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam
 
 	for(uint64_t i=0; i<nentries; i++)
 	{
+		count++;
+		if(count == flushVal)
+		{
+			count = 0;
+			flushCount++;
+			std::cout<<"\rPercent of data processed: "<<flushCount*flushFrac*100<<"%"<<std::flush;
+		}
 		tree->GetEntry(i);
 		for(Mask::Nucleus& nuc : *(dataHandle))
 		{
 			FillData(nuc);
 		}
 	}
-
+	std::cout<<std::endl;
 	input->Close();
 	delete dataHandle;
 
@@ -68,15 +75,12 @@ 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_ph_name.c_str(), ke_vs_ph_title.c_str(), 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, nuc.vec4.Phi()*s_rad2deg, 2);
+	MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
+	MyFill(angdist_name.c_str(), angdist_title.c_str(),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.GetKE(), 2);
-		MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), 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, nuc.vec4.Phi()*s_rad2deg, 2);
-		MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
-		MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM));
-	}
-	else
 	{
 		MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
 		MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.vec4.Phi()*s_rad2deg, nuc.detectedKE, 4);