diff --git a/include/AnasenEfficiency.h b/include/AnasenEfficiency.h
index de2f276..fbe1ec8 100644
--- a/include/AnasenEfficiency.h
+++ b/include/AnasenEfficiency.h
@@ -2,7 +2,6 @@
 #define ANASEN_EFFICIENCY_H
 
 #include <string>
-#include <THashTable.h>
 
 #include "DetectorEfficiency.h"
 #include "StripDetector.h"
@@ -12,16 +11,11 @@ class AnasenEfficiency : public DetectorEfficiency {
 public:
 	AnasenEfficiency();
 	~AnasenEfficiency();
-	void CalculateEfficiency(const std::string& filename) override;
+	void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
 	void DrawDetectorSystem(const std::string& filename) override;
 	double RunConsistencyCheck() override;
 
 private:
-	void RunDecay(const std::string& filename) override;
-	void Run2Step(const std::string& filename) override;
-	void Run3Step(const std::string& filename) override;
-	void Run1Step(const std::string& filename) override;
-
 	bool IsRing1(double theta, double phi);
 	bool IsRing2(double theta, double phi);
 	bool IsQQQ(double theta, double phi);
diff --git a/include/AngularDistribution.h b/include/AngularDistribution.h
index ef9f3a8..393699f 100644
--- a/include/AngularDistribution.h
+++ b/include/AngularDistribution.h
@@ -1,9 +1,9 @@
 #ifndef ANGULARDISTRIBUTION_H
 #define ANGULARDISTRIBUTION_H
 
-#include <TRandom3.h>
 #include <string>
 #include <vector>
+#include <random>
 
 class AngularDistribution {
 public:
@@ -11,7 +11,7 @@ public:
 	AngularDistribution(const std::string& file);
 	~AngularDistribution();
 	void ReadDistributionFile(const std::string& file);
-	void AttachRandomNumberGenerator(TRandom3* random) { generator = random; };
+	void AttachRandomNumberGenerator(std::mt19937* random) { generator = random; };
 	double GetRandomCosTheta();
 	int GetL() { return L; };
 	double GetBranchingRatio() { return branchingRatio; };
@@ -26,7 +26,9 @@ private:
 		}
 	}
 
-	TRandom3* generator; //NOT OWNED BY ANGULAR DISTRIBUTION
+	std::mt19937* generator; //NOT OWNED BY ANGULAR DISTRIBUTION
+	std::uniform_real_distribution<double> uniform_cosine_dist;
+	std::uniform_real_distribution<double> uniform_prob_dist;
 
 	double branchingRatio;
 	int L;
diff --git a/include/DecaySystem.h b/include/DecaySystem.h
index d435561..8ab0bfd 100644
--- a/include/DecaySystem.h
+++ b/include/DecaySystem.h
@@ -6,36 +6,36 @@
 
 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;
+	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(std::mt19937* 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/DetectorEfficiency.h b/include/DetectorEfficiency.h
index 7c24e0c..e0a1385 100644
--- a/include/DetectorEfficiency.h
+++ b/include/DetectorEfficiency.h
@@ -1,52 +1,22 @@
 #ifndef DETECTOREFFICIENCY_H
 #define DETECTOREFFICIENCY_H
 
-#include <THashTable.h>
-#include <TH1.h>
-#include <TH2.h>
 #include <string>
+#include <cmath>
 
 class DetectorEfficiency {
 public:
-	DetectorEfficiency() { m_rxn_type = -1; };
+	DetectorEfficiency() {};
 	virtual ~DetectorEfficiency() {};
 
-	inline void SetReactionType(int rxntype) { m_rxn_type = rxntype; };
-	virtual void CalculateEfficiency(const std::string& filename) = 0;
+	virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) = 0;
 	virtual void DrawDetectorSystem(const std::string& filename) = 0;
 	virtual double RunConsistencyCheck() = 0;
 
 protected:
-	void MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val) {
-		TH1F* h = (TH1F*) table->FindObject(name.c_str());
-		if(h) {
-			h->Fill(val);
-		} else {
-			h = new TH1F(name.c_str(), title.c_str(), bins, min, max);
-			h->Fill(val);
-			table->Add(h);
-		}
-	}
-    void MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy) {
-    	TH2F* h = (TH2F*) table->FindObject(name.c_str());
-		if(h) {
-			h->Fill(valx, valy);
-		} else {
-			h = new TH2F(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy);
-			h->Fill(valx, valy);
-			table->Add(h);
-		}
-    }
 	inline bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < epsilon ? true : false; };
 
-
-    virtual void Run2Step(const std::string& filename) = 0;
-	virtual void Run3Step(const std::string& filename) = 0;
-	virtual void RunDecay(const std::string& filename) = 0;
-	virtual void Run1Step(const std::string& filename) = 0;
-
 	static constexpr double epsilon = 1.0e-6;
-	int m_rxn_type;
 };
 
 #endif
\ No newline at end of file
diff --git a/include/Kinematics.h b/include/Kinematics.h
index a015cee..574a59d 100644
--- a/include/Kinematics.h
+++ b/include/Kinematics.h
@@ -6,66 +6,45 @@
 #include "OneStepSystem.h"
 #include "TwoStepSystem.h"
 #include "ThreeStepSystem.h"
-#include "Plotter.h"
 
-#include <TFile.h>
-#include <TTree.h>
+#include <random>
 
 namespace Mask {
-//For tree
-struct NucData {
-	double KE = -1;
-	double E = -1;
-	double Ex = -1;
-	double p = -1;
-	double theta = -1;
-	double theta_cm = -1;
-	double phi = -1;
-	int Z = -1;
-	int A = -1;
-};
 
-class Kinematics {
-public:
-	Kinematics();
-	~Kinematics();
-	bool LoadConfig(const char* filename);
-	bool SaveConfig(const char* filename);
-	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().c_str(); };
-	inline const char* GetOutputName() { return m_outfile_name.c_str(); };
-	inline int GetReactionType() { return m_rxn_type; };
-	void Run();
-
-	enum RxnType {
-		ONESTEP_RXN,
-		ONESTEP_DECAY,
-		TWOSTEP,
-		THREESTEP
+	class Kinematics {
+	public:
+		Kinematics();
+		~Kinematics();
+		bool LoadConfig(const std::string& filename);
+		bool SaveConfig(const std::string& filename);
+		inline int GetNumberOfSamples() { return m_nsamples; };
+		inline const std::string GetSystemName() { return sys == nullptr ? "" : sys->GetSystemEquation(); };
+		inline const std::string GetOutputName() { return m_outfile_name; };
+		inline const int GetReactionType() { return m_rxn_type; };
+		void Run();
+	
+		enum RxnType {
+			ONESTEP_DECAY,
+			ONESTEP_RXN,
+			TWOSTEP,
+			THREESTEP
+		};
+	
+	private:
+		void RunOneStepRxn();
+		void RunOneStepDecay();
+		void RunTwoStep();
+		void RunThreeStep();
+	
+		ReactionSystem* sys;
+	
+		std::string m_outfile_name;
+	
+		int m_rxn_type, m_nsamples;
+	
+		std::mt19937* global_generator;
 	};
 
-private:
-	void RunOneStepRxn();
-	void RunOneStepDecay();
-	void RunTwoStep();
-	void RunThreeStep();
-
-	ReactionSystem* sys;
-	Plotter plotman;
-
-	NucData ConvertNucleus(const Nucleus& nuc);
-
-	bool save_tree_flag, do_plotter_flag;
-
-	std::string m_outfile_name;
-
-	int m_rxn_type, m_nsamples;
-
-	TRandom3* global_generator;
-};
-
-};
+}
 
 #endif
diff --git a/include/MaskFile.h b/include/MaskFile.h
index b584f89..e0a8bae 100644
--- a/include/MaskFile.h
+++ b/include/MaskFile.h
@@ -9,45 +9,67 @@
 
 namespace Mask {
 
-class MaskFile {
-public:
-	enum class FileType {
-		read,
-		write,
-		append,
-		none
+	struct MaskFileHeader {
+		int rxn_type = -1;
+		int nsamples = -1;
 	};
 
-	MaskFile();
-	MaskFile(const std::string& name, MaskFile::FileType type);
-	bool Open(const std::string& name, MaskFile::FileType type);
-	inline bool IsOpen() { return file.is_open(); }
-	void Close();
+	struct MaskFileData {
+		std::vector<double> E, KE, p, theta, phi; //ordered: target, (if not decay)projectile, ejectile, residual, break1...
+		std::vector<int> Z, A;
+		std::vector<bool> detect_flag;
+		bool eof = false; //flag on end of file
+	};
+
+	class MaskFile {
+	public:
+		enum class FileType {
+			read,
+			write,
+			append,
+			none
+		};
 	
-	void WriteHeader(std::vector<Nucleus>& data);
-	void WriteData(std::vector<Nucleus>& data);
-	void ReadHeader();
-	std::vector<Nucleus> ReadData();
-
+		MaskFile();
+		MaskFile(const std::string& name, MaskFile::FileType type);
+		bool Open(const std::string& name, MaskFile::FileType type);
+		inline bool IsOpen() { return file.is_open(); }
+		void Close();
+		
+		void WriteHeader(int rxn_type, int nsamples);
+		void WriteData(std::vector<Nucleus>& data);
+		void WriteData(MaskFileData& data);
+		MaskFileHeader ReadHeader();
+		MaskFileData ReadData();
 	
+		
+	
+	private:
 
-private:
-	void FlushBuffer();
-	void FillBuffer();
+		FileType file_type;
+		std::string filename;
+		unsigned int buffer_position;
+		unsigned int buffer_end;
+		unsigned int data_size;
+		int m_rxn_type;
+		int buffersize_bytes;
+	
+		std::fstream file;
+	
+		std::vector<char> data_buffer;
+	
+		static constexpr int onestep_rxn_n = 2;
+		static constexpr int twostep_rxn_n = 4;
+		static constexpr int threestep_rxn_n = 6;
+	
+		static constexpr int buffersize = 10000; //number of events
+		static constexpr int width = 0;
+		static constexpr int precision = 3;
 
-	FileType file_type;
-	std::string filename;
-	unsigned int buffer_position;
-	int count;
-
-	std::fstream file;
-
-	std::vector<std::vector<Nucleus>> data_buffer;
-
-	static constexpr int buffersize = 1000;
-	static constexpr int width = 0;
-	static constexpr int precision = 3;
-};
+		static constexpr std::size_t double_size = sizeof(double);
+		static constexpr std::size_t int_size = sizeof(int);
+		static constexpr std::size_t bool_size = sizeof(bool);
+	};
 
 };
 
diff --git a/include/MassLookup.h b/include/MassLookup.h
index 404116f..ba7bd52 100644
--- a/include/MassLookup.h
+++ b/include/MassLookup.h
@@ -7,6 +7,7 @@ of this map (MASS) for use throughout code it is included into.
 
 Written by G.W. McCann Aug. 2020
 
+Converted to true singleton to simplify usage -- Aug. 2021 GWM
 */
 #ifndef MASS_LOOKUP_H
 #define MASS_LOOKUP_H
@@ -20,21 +21,27 @@ Written by G.W. McCann Aug. 2020
 class MassLookup {
 
   public:
-    MassLookup();
     ~MassLookup();
     double FindMass(int Z, int A);
     std::string FindSymbol(int Z, int A);
+    static MassLookup* GetInstance() {
+        if(s_instance == nullptr) {
+          s_instance = new MassLookup();
+        }
+        return s_instance;
+    }
 
   private:
+    MassLookup();
     std::unordered_map<std::string, double> massTable;
     std::unordered_map<int, std::string> elementTable;
 
+    static MassLookup* s_instance;
+
     //constants
     static constexpr double u_to_mev = 931.4940954;
     static constexpr double electron_mass = 0.000548579909;
     
 };
 
-//static instance for use throught program
-static MassLookup MASS;
 #endif
diff --git a/include/Nucleus.h b/include/Nucleus.h
index 8a42f52..70cca6f 100644
--- a/include/Nucleus.h
+++ b/include/Nucleus.h
@@ -28,6 +28,9 @@ public:
 	inline double GetGroundStateMass() const { return m_gs_mass; };
 	inline const char* GetIsotopicSymbol() const { return m_symbol.c_str(); };
 	inline double GetThetaCM() const { return m_theta_cm; };
+	inline void SetDetected() { m_detectFlag = true; };
+	inline void SetNotDetected() { m_detectFlag = false; };
+	inline bool IsDetected() { return m_detectFlag; };
 
 	inline Nucleus& operator=(const Nucleus& rhs) {
 		SetIsotope(rhs.GetZ(), rhs.GetA());
@@ -52,6 +55,8 @@ private:
 	double m_theta_cm;
 	std::string m_symbol;
 
+	bool m_detectFlag;
+
 };
 
 };
diff --git a/include/OneStepSystem.h b/include/OneStepSystem.h
index 6b58fc4..06a276b 100644
--- a/include/OneStepSystem.h
+++ b/include/OneStepSystem.h
@@ -5,29 +5,29 @@
 
 namespace Mask {
 
-class OneStepSystem: public ReactionSystem {
-public:
-	OneStepSystem();
-	OneStepSystem(std::vector<int>& z, std::vector<int>& a);
-	~OneStepSystem();
+	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;
+	
+	};
 
-	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/Plotter.h b/include/Plotter.h
deleted file mode 100644
index 5e4bb16..0000000
--- a/include/Plotter.h
+++ /dev/null
@@ -1,49 +0,0 @@
-#ifndef PLOTTER_H
-#define PLOTTER_H
-
-#include <TROOT.h>
-#include <TH1.h>
-#include <TH2.h>
-#include <TGraph.h>
-#include <THashTable.h>
-#include "Nucleus.h"
-
-namespace Mask {
-
-struct GraphData {
-	std::string name;
-	std::string title;
-	std::vector<double> xvec;
-	std::vector<double> yvec;
-	int color;
-};
-
-class Plotter {
-public:
-	Plotter();
-	~Plotter();
-	inline void ClearTable() { table->Clear(); };
-	inline THashTable* GetTable() {
-		GenerateGraphs(); 
-		return table; 
-	};
-	void FillData(const Nucleus& nuc, const std::string& modifier = "");
-
-private:
-	THashTable* table;
-
-	void GenerateGraphs();
-	void MyFill(const char* name, const char* title, int bins, float min, float max, double val);
-	void MyFill(const char* name, const char* title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy);
-	void MyFill(const char* name, const char* title, double valx, double valy, int color); //TGraph
-
-	std::vector<TGraph*> garbage_collection;
-	std::vector<GraphData> graphs;
-
-	static constexpr double rad2deg = 180.0/M_PI;
-
-};
-
-};
-
-#endif
\ No newline at end of file
diff --git a/include/QQQDetector.h b/include/QQQDetector.h
index 3ffee56..9ee6d90 100644
--- a/include/QQQDetector.h
+++ b/include/QQQDetector.h
@@ -8,8 +8,8 @@
 #ifndef QQQDETECTOR_H
 #define QQQDETECTOR_H
 
-#include <TRandom3.h>
 #include <cmath>
+#include <vector>
 
 #include "Vec3.h"
 #include "Rotation.h"
diff --git a/include/ReactionSystem.h b/include/ReactionSystem.h
index 4aa885f..7e64138 100644
--- a/include/ReactionSystem.h
+++ b/include/ReactionSystem.h
@@ -12,45 +12,65 @@
 #include "Reaction.h"
 #include "KinematicsExceptions.h"
 #include <vector>
-#include <TRandom3.h>
+#include <random>
 
 namespace Mask {
 
-class ReactionSystem {
-public:
-	ReactionSystem();
-	virtual ~ReactionSystem();
-
-	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*/
-	virtual void SetRandomGenerator(TRandom3* gen);
-	inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); };
-	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 const std::string& GetSystemEquation() const { return m_sys_equation; };
-
-protected:
-	virtual void LinkTarget() = 0;
-	virtual void SetSystemEquation() = 0;
+	class ReactionSystem {
+	public:
+		ReactionSystem();
+		virtual ~ReactionSystem();
 	
-	LayeredTarget target;
+		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*/
+		virtual void SetRandomGenerator(std::mt19937* gen);
+		inline void SetBeamDistro(double mean, double sigma) {
+			if(m_beamDist)
+				delete m_beamDist;
+			m_beamDist = new std::normal_distribution<double>(mean, sigma); 
+		}
+	
+		inline void SetTheta1Range(double min, double max) { 
+			if(m_theta1Range)
+				delete m_theta1Range;
+			m_theta1Range = new std::uniform_real_distribution<double>(std::cos(min*deg2rad), std::cos(max*deg2rad)); 
+		}
+	
+		inline void SetPhi1Range(double min, double max) {
+			if(m_phi1Range)
+				delete m_phi1Range;
+			m_phi1Range = new std::uniform_real_distribution<double>(min*deg2rad, max*deg2rad); 
+		}
+	
+		inline void SetExcitationDistro(double mean, double sigma) {
+			if(m_exDist)
+				delete m_exDist;
+			m_exDist = new std::normal_distribution<double>(mean, sigma); 
+		}
+	
+		inline const std::string& GetSystemEquation() const { return m_sys_equation; }
+	
+	protected:
+		virtual void LinkTarget() = 0;
+		virtual void SetSystemEquation() = 0;
+		
+		LayeredTarget target;
+	
+		//Sampling information
+		std::normal_distribution<double> *m_beamDist, *m_exDist;
+		std::uniform_real_distribution<double> *m_theta1Range, *m_phi1Range; 
+		std::mt19937* generator; //not owned by ReactionSystem
+	
+		bool target_set_flag, gen_set_flag;
+		int rxnLayer;
+		std::string m_sys_equation;
+		static constexpr double deg2rad = M_PI/180.0;
+	};
 
-	//Sampling information
-	std::pair<double, double> m_beamDist, m_theta1Range, m_phi1Range, m_exDist;
-	TRandom3* generator; //not owned by ReactionSystem
-
-	bool target_set_flag, gen_set_flag;
-	int rxnLayer;
-	std::string m_sys_equation;
-	static constexpr double deg2rad = M_PI/180.0;
-};
-
-};
+}
 
 #endif
\ No newline at end of file
diff --git a/include/RootPlotter.h b/include/RootPlotter.h
new file mode 100644
index 0000000..2ea5559
--- /dev/null
+++ b/include/RootPlotter.h
@@ -0,0 +1,48 @@
+#ifndef ROOTPLOTTER_H
+#define ROOTPLOTTER_H
+
+#include <vector>
+#include <string>
+
+#include "Nucleus.h"
+
+#include <THashTable.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TGraph.h>
+
+struct GraphData {
+	std::string name;
+	std::string title;
+	std::vector<double> xvec;
+	std::vector<double> yvec;
+	int color;
+};
+
+class RootPlotter {
+public:
+	RootPlotter();
+	~RootPlotter();
+	inline void ClearTable() { table->Clear(); };
+	inline THashTable* GetTable() {
+		GenerateGraphs(); 
+		return table; 
+	};
+	void FillData(const Mask::Nucleus& nuc, const std::string& modifier = "");
+
+private:
+	THashTable* table;
+
+	void GenerateGraphs();
+	void MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val);
+	void MyFill(const std::string& name, const std::string& title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy);
+	void MyFill(const std::string& name, const std::string& title, double valx, double valy, int color); //TGraph
+
+	std::vector<TGraph*> garbage_collection;
+	std::vector<GraphData> graphs;
+
+	static constexpr double rad2deg = 180.0/M_PI;
+
+};
+
+#endif
\ No newline at end of file
diff --git a/include/SabreEfficiency.h b/include/SabreEfficiency.h
index 6bee112..b795acf 100644
--- a/include/SabreEfficiency.h
+++ b/include/SabreEfficiency.h
@@ -5,24 +5,19 @@
 #include "SabreDetector.h"
 #include "Target.h"
 #include "DeadChannelMap.h"
-#include "Kinematics.h"
-#include <THashTable.h>
+#include "Nucleus.h"
 
 class SabreEfficiency : public DetectorEfficiency {
 public:
 	SabreEfficiency();
 	~SabreEfficiency();
     void SetDeadChannelMap(std::string& filename) { dmap.LoadMapfile(filename); };
-	void CalculateEfficiency(const std::string& file) override;
+	void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
     void DrawDetectorSystem(const std::string& filename) override;
     double RunConsistencyCheck() override;
 
 private:
-	std::pair<bool,double> IsSabre(Mask::NucData* nucleus);
-    void Run2Step(const std::string& filename) override;
-	void Run3Step(const std::string& filename) override;
-    void RunDecay(const std::string& filename) override;
-	void Run1Step(const std::string& filename) override;
+	std::pair<bool,double> IsSabre(Mask::Nucleus& nucleus);
 
 	std::vector<SabreDetector> detectors;
     
diff --git a/include/Stopwatch.h b/include/Stopwatch.h
new file mode 100644
index 0000000..af217c1
--- /dev/null
+++ b/include/Stopwatch.h
@@ -0,0 +1,30 @@
+/*
+	Stopwatch.h
+	Simple class designed to provide timing info on parts of the process.
+	Only for use in development.
+
+	Written by G.W. McCann Oct. 2020
+*/
+#ifndef STOPWATCH_H
+#define STOPWATCH_H
+
+#include <chrono>
+
+class Stopwatch {
+
+public:
+	Stopwatch();
+	~Stopwatch();
+	void Start();
+	void Stop();
+	double GetElapsedSeconds();
+	double GetElapsedMilliseconds();
+
+private:
+	using Time = std::chrono::high_resolution_clock::time_point;
+	using Clock = std::chrono::high_resolution_clock;
+
+	Time start_time, stop_time;
+};
+
+#endif
\ No newline at end of file
diff --git a/include/StripDetector.h b/include/StripDetector.h
index 2aaef13..7479cf0 100644
--- a/include/StripDetector.h
+++ b/include/StripDetector.h
@@ -11,9 +11,9 @@
 
 //Back strips from lowest z to highest z
 
-#include <TRandom3.h>
 #include <cmath>
 #include <vector>
+#include <random>
 
 #include "Vec3.h"
 #include "Rotation.h"
@@ -27,7 +27,7 @@ public:
 	inline Mask::Vec3 GetBackStripCoordinates(int stripch, int corner) { return back_strip_coords[stripch][corner]; };
 	inline Mask::Vec3 GetRotatedFrontStripCoordinates(int stripch, int corner) { return rotated_front_strip_coords[stripch][corner]; };
 	inline Mask::Vec3 GetRotatedBackStripCoordinates(int stripch, int corner) { return rotated_back_strip_coords[stripch][corner]; };
-	inline void SetRandomNumberGenerator(TRandom3* random) { m_random = random; };
+	inline void SetRandomNumberGenerator(std::mt19937* random) { m_random = random; };
 	Mask::Vec3 GetHitCoordinates(int front_stripch, double front_strip_ratio);
 	std::pair<int,double> GetChannelRatio(double theta, double phi);
 
@@ -51,7 +51,9 @@ private:
 
 	Mask::ZRotation zRot;
 
-	TRandom3* m_random; //Not owned by StripDetector!
+	std::uniform_real_distribution<double> m_uniform_fraction;
+
+	std::mt19937* m_random; //Not owned by StripDetector!
 
 	inline bool ValidChannel(int f) { return ((f >= 0 && f < num_strips) ? true : false); };
 	inline bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); };
diff --git a/include/ThreeStepSystem.h b/include/ThreeStepSystem.h
index 4d95356..76447c8 100644
--- a/include/ThreeStepSystem.h
+++ b/include/ThreeStepSystem.h
@@ -6,43 +6,45 @@
 
 namespace Mask {
 
-class ThreeStepSystem : public ReactionSystem {
-public:
-	ThreeStepSystem();
-	ThreeStepSystem(std::vector<int>& z, std::vector<int>& a);
-	~ThreeStepSystem();
-	bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
-	void RunSystem() override;
-	void SetRandomGenerator(TRandom3* gen) override;
+	class ThreeStepSystem : public ReactionSystem {
+	public:
+		ThreeStepSystem();
+		ThreeStepSystem(std::vector<int>& z, std::vector<int>& a);
+		~ThreeStepSystem();
+		bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
+		void RunSystem() override;
+		void SetRandomGenerator(std::mt19937* 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(); };
+	
+		inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
+		inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); };
+		inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
+		inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); };
+	
+	protected:
+		void LinkTarget() override;
+		void SetSystemEquation() override;
 
-	inline void SetDecay1Distribution(const std::string& filename) { decay1dist.ReadDistributionFile(filename); };
-	inline void SetDecay2Distribution(const std::string& filename) { decay2dist.ReadDistributionFile(filename); };
+		std::uniform_real_distribution<double> m_phi2Range;
+	
+		Reaction step1, step2, step3;
+	
+		AngularDistribution decay1dist, decay2dist;
+	
+	};
 
-	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(); };
-
-	inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); };
-	inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); };
-	inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); };
-	inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); };
-
-protected:
-	void LinkTarget() override;
-	void SetSystemEquation() override;
-
-	Reaction step1, step2, step3;
-
-	AngularDistribution decay1dist, decay2dist;
-
-};
-
-};
+}
 
 #endif
\ No newline at end of file
diff --git a/include/TwoStepSystem.h b/include/TwoStepSystem.h
index 44c927a..c985069 100644
--- a/include/TwoStepSystem.h
+++ b/include/TwoStepSystem.h
@@ -6,39 +6,41 @@
 
 namespace Mask {
 
-class TwoStepSystem : public ReactionSystem {
-public:
-	TwoStepSystem();
-	TwoStepSystem(std::vector<int>& z, std::vector<int>& a);
-	~TwoStepSystem();
-	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() override;
-	void SetSystemEquation() override;
+	class TwoStepSystem : public ReactionSystem {
+	public:
+		TwoStepSystem();
+		TwoStepSystem(std::vector<int>& z, std::vector<int>& a);
+		~TwoStepSystem();
+		bool SetNuclei(std::vector<int>& z, std::vector<int>& a) override;
+		void RunSystem() override;
 	
-	Reaction step1, step2;
+		void SetRandomGenerator(std::mt19937* 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() override;
+		void SetSystemEquation() override;
 
-	AngularDistribution decay1dist;
+		std::uniform_real_distribution<double> m_phi2Range;
+		
+		Reaction step1, step2;
+	
+		AngularDistribution decay1dist;
+	
+	};
 
-};
-
-};
+}
 
 #endif
\ No newline at end of file
diff --git a/input.txt b/input.txt
index 69ced3c..bc9234c 100644
--- a/input.txt
+++ b/input.txt
@@ -1,5 +1,5 @@
 ----------Data Information----------
-OutputFile: /data1/gwm17/mask_tests/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.root
+OutputFile: /data1/gwm17/mask_tests/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.mask
 SaveTree: yes
 SavePlots: yes
 ----------Reaction Information----------
@@ -20,8 +20,8 @@ Z	A	Stoich
 0
 ~
 ----------Sampling Information----------
-NumberOfSamples: 100000
-BeamMeanEnergy(MeV): 0.6 BeamEnergySigma(MeV): 0.001
+NumberOfSamples: 10000
+BeamMeanEnergy(MeV): 0.6 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
diff --git a/premake5.lua b/premake5.lua
index 4cd98ee..1eb4b76 100644
--- a/premake5.lua
+++ b/premake5.lua
@@ -12,14 +12,35 @@ project "Mask"
 	cppdialect "C++11"
 
 	files {
-		"src/**.cpp",
-		"include/**.h",
-		"src/**.cxx"
+		"src/*.cpp",
+		"include/*.h"
 	}
 
-	prebuildcommands {
-		"rootcint -f src/kinematics_dict.cxx include/Kinematics.h include/LinkDef_Kinematics.h",
-		"{MOVE} src/kinematics_dict_rdict.pcm bin/"
+	includedirs {
+		"include"
+	}
+
+	filter "configurations:Debug"
+		symbols "On"
+
+	filter "configurations:Release"
+		optimize "On"
+
+project	"RootPlot"
+	kind "ConsoleApp"
+	language "C++"
+	targetdir "bin"
+	objdir "objs"
+	cppdialect "c++11"
+
+	files {
+		"src/Plotters/ROOT/RootPlotter.cpp",
+		"src/MaskFile.cpp",
+		"src/Nucleus.cpp",
+		"src/Vec4.cpp",
+		"src/Vec3.cpp",
+		"src/MassLookup.cpp",
+		"include/*.h"
 	}
 
 	filter "system:windows"
@@ -54,6 +75,36 @@ project "Mask"
 			"`root-config --glibs`"
 		}
 
+	filter "configurations:Debug"
+		symbols "On"
+
+	filter "configurations:Release"
+		optimize "On"
+
+project	"DetectEff"
+	kind "ConsoleApp"
+	language "C++"
+	targetdir "bin"
+	objdir "objs"
+	cppdialect "c++11"
+
+	files {
+		"src/Detectors/*.cpp",
+		"src/MaskFile.cpp",
+		"src/Nucleus.cpp",
+		"src/Vec4.cpp",
+		"src/Vec3.cpp",
+		"src/MassLookup.cpp",
+		"src/Rotation.cpp",
+		"src/Target.cpp",
+		"src/EnergyLoss.cpp",
+		"include/*.h"
+	}
+
+	includedirs {
+		"include"
+	}
+
 	filter "configurations:Debug"
 		symbols "On"
 
diff --git a/src/AngularDistribution.cpp b/src/AngularDistribution.cpp
index 453b862..1065c71 100644
--- a/src/AngularDistribution.cpp
+++ b/src/AngularDistribution.cpp
@@ -5,7 +5,7 @@
 #include "LegendrePoly.h"
 
 AngularDistribution::AngularDistribution() :
-	generator(nullptr), branchingRatio(1.0), L(0), isoFlag(true)
+	generator(nullptr), uniform_cosine_dist(-1.0, 1.0), uniform_prob_dist(0.0, 1.0), branchingRatio(1.0), L(0), isoFlag(true)
 {
 }
 
@@ -84,18 +84,18 @@ double AngularDistribution::GetRandomCosTheta() {
 		return 0.0;
 	}
 
-	if(isoFlag) return generator->Uniform(-1.0, 1.0);
+	if(isoFlag) return uniform_cosine_dist(*generator);
 
 	double test, probability;
 	double costheta;
 
-	test = generator->Uniform(0.0, 1.0);
+	test = uniform_prob_dist(*generator);
 	if(test > branchingRatio) return -10;
 
 	do {
 		probability = 0.0;
-		costheta = generator->Uniform(-1.0, 1.0);
-		test = generator->Uniform(0.0, 1.0);
+		costheta = uniform_cosine_dist(*generator);
+		test = uniform_prob_dist(*generator);
 		for(unsigned int i=0; i<constants.size(); i++)
 			probability += constants[i]*P_l(i*2, costheta);
 	} while(test > probability);
diff --git a/src/DecaySystem.cpp b/src/DecaySystem.cpp
index b18165a..eaa0ba0 100644
--- a/src/DecaySystem.cpp
+++ b/src/DecaySystem.cpp
@@ -2,74 +2,74 @@
 
 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;
+	DecaySystem::DecaySystem() :
+		ReactionSystem()
+	{
 	}
-
-	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();
+	DecaySystem::DecaySystem(std::vector<int>& z, std::vector<int>& a) :
+		ReactionSystem()
+	{
+		SetNuclei(z, a);
 	}
-
-	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;
+	
+	DecaySystem::~DecaySystem() {}
+	
+	void DecaySystem::SetRandomGenerator(std::mt19937* 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 = (*m_phi1Range)(*generator);
+			step1.SetPolarRxnAngle(rxnTheta);
+			step1.SetAzimRxnAngle(rxnPhi);
+	
+			step1.TurnOnResidualEloss();
+			step1.Calculate();
+		} else {
+			return;
+		}
 	}
-}
 
 }
\ No newline at end of file
diff --git a/src/Detectors/AnasenEfficiency.cpp b/src/Detectors/AnasenEfficiency.cpp
index 2c43f19..117644a 100644
--- a/src/Detectors/AnasenEfficiency.cpp
+++ b/src/Detectors/AnasenEfficiency.cpp
@@ -1,11 +1,8 @@
 #include "AnasenEfficiency.h"
 #include "Kinematics.h"
-#include <TH2.h>
-#include <TH1.h>
-#include <TGraph2D.h>
-#include <TGraph.h>
-#include <TCanvas.h>
-#include <TParameter.h>
+#include "MaskFile.h"
+#include <fstream>
+#include <iomanip>
 
 AnasenEfficiency::AnasenEfficiency() :
 	DetectorEfficiency()
@@ -24,7 +21,7 @@ AnasenEfficiency::~AnasenEfficiency() {}
 
 
 void AnasenEfficiency::DrawDetectorSystem(const std::string& filename) {
-	TFile* file = TFile::Open(filename.c_str(), "RECREATE");
+	std::ofstream output(filename);
 
 	std::vector<double> x, y, z;
 	std::vector<double> cx, cy, cz;
@@ -105,24 +102,18 @@ void AnasenEfficiency::DrawDetectorSystem(const std::string& filename) {
 			}
 		}
 	}
-	TGraph2D* graph = new TGraph2D(x.size(), &(x[0]), &(y[0]), &(z[0]));
-	graph->SetName("CornerGraph");
-	graph->SetMarkerStyle(2);
-	graph->SetLineColor(1);
 
-	TGraph2D* graph2 = new TGraph2D(cx.size(), &(cx[0]), &(cy[0]), &(cz[0]));
-	graph2->SetName("CenterGraph");
-	graph2->SetMarkerStyle(2);
-	graph2->SetMarkerColor(4);
+	output<<"ANASEN Geometry File -- Coordinates for Detectors"<<std::endl;
+	output<<"Edges: x y z"<<std::endl;
+	for(unsigned int i=0; i<x.size(); i++) {
+		output<<x[i]<<" "<<y[i]<<" "<<z[i]<<std::endl;
+	}
+	output<<"Centers: x y z"<<std::endl;
+	for(unsigned int i=0; i<cx.size(); i++) {
+		output<<cx[i]<<" "<<cy[i]<<" "<<cz[i]<<std::endl;
+	}
 
-	TCanvas* canvas = new TCanvas();
-	canvas->SetName("ANASEN Detector");
-	graph->Draw("A|P");
-	graph2->Draw("same|P");
-	canvas->Write();
-	graph->Write();
-	graph2->Write();
-	file->Close();
+	output.close();
 }
 
 double AnasenEfficiency::RunConsistencyCheck() {
@@ -248,305 +239,86 @@ bool AnasenEfficiency::IsQQQ(double theta, double phi) {
 	return false;
 }
 
-void AnasenEfficiency::CalculateEfficiency(const std::string& file) {
+void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) {
 	std::cout<<"----------ANASEN Efficiency Calculation----------"<<std::endl;
-	std::cout<<"Loading in output from kinematics simulation: "<<file<<std::endl;
+	std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
 	std::cout<<"Running efficiency calculation..."<<std::endl;
 
-	switch(m_rxn_type) {
-		case Mask::Kinematics::ONESTEP_DECAY:
-		{
-			RunDecay(file);
+	Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
+	Mask::MaskFile output(outputname, Mask::MaskFile::FileType::write);
+	std::ofstream stats(statsname);
+	stats<<std::setprecision(5);
+
+	Mask::MaskFileHeader header = input.ReadHeader();
+	output.WriteHeader(header.rxn_type, header.nsamples);
+	stats<<"Efficiency statistics for data from "<<inputname<<" using the ANASEN geometry"<<std::endl;
+	stats<<"Given in order of target=0, projectile=1, ejectile=2, residual=3, .... etc."<<std::endl;
+
+	Mask::MaskFileData data;
+
+	std::vector<int> counts;
+	switch(header.rxn_type) {
+		case 0:
+			counts.resize(3, 0);
 			break;
-		}
-		case Mask::Kinematics::ONESTEP_RXN:
-		{
-			Run1Step(file);
-		}
-		case Mask::Kinematics::TWOSTEP:
-		{
-			Run2Step(file);
+		case 1:
+			counts.resize(4, 0);
 			break;
-		}
-		case Mask::Kinematics::THREESTEP:
-		{
-			Run3Step(file);
+		case 2:
+			counts.resize(6, 0);
 			break;
+		case 3:
+			counts.resize(8, 0);
+			break;
+		default:
+		{
+			std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<header.rxn_type<<"). Quiting..."<<std::endl;
+			input.Close();
+			output.Close();
+			stats.close();
+			return;
 		}
 	}
+
+	int percent5 = header.nsamples*0.05;
+	int count = 0;
+	int npercent = 0;
+
+	while(true) {
+		if(++count == percent5) {//Show progress every 5%
+			npercent++;
+			count = 0;
+			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
+		}
+
+		data = input.ReadData();
+		if(data.eof)
+			break;
+
+		for(unsigned int i=0; i<data.Z.size(); i++) {
+			if(data.KE[i] >= threshold && (IsRing1(data.theta[i], data.phi[i]) || IsRing2(data.theta[i], data.phi[i]) || IsQQQ(data.theta[i], data.phi[i]))) {
+				data.detect_flag[i] = true;
+				counts[i]++;
+			} else if(data.detect_flag[i] == true) {
+				data.detect_flag[i] = false;
+			}
+		}
+
+		output.WriteData(data);
+	}
+
+	input.Close();
+	output.Close();
+
+	stats<<std::setw(10)<<"Index"<<"|"<<std::setw(10)<<"Efficiency"<<std::endl;
+	stats<<"---------------------"<<std::endl;
+	for(unsigned int i=0; i<counts.size(); i++) {
+		stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+	}
+	stats.close();
+
 	std::cout<<std::endl;
 	std::cout<<"Complete."<<std::endl;
 	std::cout<<"---------------------------------------------"<<std::endl;
-}
-
-void AnasenEfficiency::RunDecay(const std::string& filename) {
-	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
-	TTree* tree = (TTree*) input->Get("DataTree");
-
-	THashTable* table = new THashTable();
-	Mask::NucData* eject = new Mask::NucData();
-	Mask::NucData* resid = new Mask::NucData();
-
-	tree->SetBranchAddress("ejectile", &eject);
-	tree->SetBranchAddress("residual", &resid);
-
-	double nevents = tree->GetEntries();
-
-	//Progress tracking
-	int percent5 = nevents*0.05;
-	int count = 0;
-	int npercent = 0;
-
-	int detected_eject=0;
-	int detected_resid=0;
-	for(int i=0; i<tree->GetEntries(); i++) {
-		if(++count == percent5) {//Show progress every 5%
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		tree->GetEntry(i);
-
-		if(eject->KE >= threshold && (IsRing1(eject->theta, eject->phi) || IsRing2(eject->theta, eject->phi) || IsQQQ(eject->theta, eject->phi))) {
-			detected_eject++;
-			MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/deg2rad,300,0.0,30.0,eject->KE);
-		}
-
-		if(resid->KE > threshold && (IsRing1(resid->theta, resid->phi) || IsRing2(resid->theta, resid->phi) || IsQQQ(resid->theta, resid->phi))) {
-			detected_resid++;
-			MyFill(table, "detected_resid_theta_ke","detected residual theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,resid->theta/deg2rad,300,0.0,30.0,resid->KE);
-		}
-
-	}
-
-	double ejecteff = ((double) detected_eject)/nevents;
-	double resideff = ((double) detected_resid)/nevents;
-	TParameter<double> eject_eff("Light Breakup Efficiency", ejecteff);
-	TParameter<double> resid_eff("Heavy Breakup Efficiency", resideff);
-
-	input->cd();
-	eject_eff.Write();
-	resid_eff.Write();
-	table->Write();
-	input->Close();
-}
-
-void AnasenEfficiency::Run1Step(const std::string& filename) {
-	std::cout<<"SabreEfficiency::Run1Step Not Implemented!"<<std::endl;
-	return;
-}
-
-void AnasenEfficiency::Run2Step(const std::string& filename) {
-
-	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
-	TTree* tree = (TTree*) input->Get("DataTree");
-
-	Mask::NucData* eject = new Mask::NucData();
-	Mask::NucData* break1 = new Mask::NucData();
-	Mask::NucData* break2 = new Mask::NucData();
-
-	THashTable* table = new THashTable();
-
-	tree->SetBranchAddress("ejectile", &eject);
-	tree->SetBranchAddress("breakup1", &break1);
-	tree->SetBranchAddress("breakup2", &break2);
-
-	double nevents = tree->GetEntries();
-
-	//Progress tracking
-	int percent5 = nevents*0.05;
-	int count = 0;
-	int npercent = 0;
-
-	double costheta_cm =0;
-
-	bool break1_flag, break2_flag, eject_flag;
-	int b1_count=0, b2_count=0, eject_count=0;
-	int b1b2_count=0, b1e_count=0, b2e_count=0, b1b2e_count=0;
-	for(int i=0; i<tree->GetEntries(); i++) {
-		if(++count == percent5) {//Show progress every 5%
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		break1_flag = false;
-		break2_flag = false;
-		eject_flag = false;
-
-		tree->GetEntry(i);
-
-		if(eject->KE >= threshold && (IsRing1(eject->theta, eject->phi) || IsRing2(eject->theta, eject->phi) || IsQQQ(eject->theta, eject->phi))) {
-			eject_count++;
-			eject_flag = true;
-			MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/deg2rad,300,0.0,30.0,eject->KE);
-		}
-
-		if(break1->KE >= threshold && (IsRing1(break1->theta, break1->phi) || IsRing2(break1->theta, break1->phi) || IsQQQ(break1->theta, break1->phi))) {
-			b1_count++;
-			break1_flag = true;
-			costheta_cm = std::cos(break1->theta_cm);
-			MyFill(table,"detected_break1_cm_theta","detected breakup1 CM theta;cos(#theta_{CM})",20,-1.0,1.0,costheta_cm);
-			MyFill(table, "detected_break1_theta_ke","detected breakup1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/deg2rad,300,0.0,30.0,break1->KE);
-		}
-
-		if(break2->KE >= threshold && (IsRing1(break2->theta, break2->phi) || IsRing2(break2->theta, break2->phi) || IsQQQ(break2->theta, break2->phi))) {
-			b2_count++;
-			break2_flag = true;
-			MyFill(table, "detected_break2_theta_ke","detected breakup2 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break2->theta/deg2rad,300,0.0,30.0,break2->KE);
-		}
-
-		if(break1_flag && break2_flag && eject_flag) {
-			b1b2e_count++;
-			b1b2_count++;
-			b1e_count++;
-			b2e_count++;
-			MyFill(table,"b1b2_theta_theta","b1b2_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break1->theta/deg2rad,180,0.0,180.0,break2->theta/deg2rad);
-			MyFill(table,"b1b2_ke_ke","b1b2_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break1->KE,300,0.0,30.0,break2->KE);
-			MyFill(table,"b1b2_phi_phi","b1b2_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break1->phi/deg2rad,360,0.0,360.0,break2->phi/deg2rad);
-		} else if(break1_flag && eject_flag) {
-			b1e_count++;
-		} else if(break2_flag && eject_flag) {
-			b2e_count++;
-		} else if(break1_flag && break2_flag) {
-			b1b2_count++;
-			MyFill(table,"b1b2_theta_theta","b1b2_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break1->theta/deg2rad,180,0.0,180.0,break2->theta/deg2rad);
-			MyFill(table,"b1b2_ke_ke","b1b2_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break1->KE,300,0.0,30.0,break2->KE);
-			MyFill(table,"b1b2_phi_phi","b1b2_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break1->phi/deg2rad,360,0.0,360.0,break2->phi/deg2rad);
-		}
-	}
-
-	double eeff = ((double) eject_count)/nevents;
-	double b1eff = ((double) b1_count)/nevents;
-	double b2eff = ((double) b2_count)/nevents;
-	double b1b2eff = b1b2_count/nevents;
-	double b1eeff = b1e_count/nevents;
-	double b2eeff = b2e_count/nevents;
-	double b1b2eeff = b1b2e_count/nevents;
-	TParameter<double> eject_eff("Ejectile Efficiency", eeff);
-	TParameter<double> break1_eff("Breakup 1 Efficiency", b1eff);
-	TParameter<double> break2_eff("Breakup 2 Efficiency", b2eff);
-	TParameter<double> break1break2_eff("Breakup1 Coincident with Breakup2", b1b2eff);
-	TParameter<double> break1eject_eff("Breakup1 Coincident with Ejectile", b1eeff);
-	TParameter<double> break2eject_eff("Breakup2 Coincident with Ejectile", b2eeff);
-	TParameter<double> break1break2eject_eff("All Three Particles Coincident", b1b2eeff);
-
-
-	input->cd();
-	table->Write();
-	eject_eff.Write();
-	break1_eff.Write();
-	break2_eff.Write();
-	break1break2_eff.Write();
-	break1eject_eff.Write();
-	break2eject_eff.Write();
-	break1break2eject_eff.Write();
-
-	input->Close();
-}
-
-void AnasenEfficiency::Run3Step(const std::string& filename) {
-	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
-	TTree* tree = (TTree*) input->Get("DataTree");
-
-	THashTable* table = new THashTable();
-
-	Mask::NucData* break1 = new Mask::NucData();
-	Mask::NucData* break3 = new Mask::NucData();
-	Mask::NucData* break4 = new Mask::NucData();
-
-
-	tree->SetBranchAddress("breakup1", &break1);
-	tree->SetBranchAddress("breakup3", &break3);
-	tree->SetBranchAddress("breakup4", &break4);
-
-	double nevents = tree->GetEntries();
-
-	//Progress tracking
-	int percent5 = nevents*0.05;
-	int count = 0;
-	int npercent = 0;
-
-	bool break1_flag, break3_flag, break4_flag;
-	int b1_count=0, b3_count=0, b4_count=0;
-	int b1b3_count=0, b1b4_count=0, b3b4_count=0, b1b3b4_count=0;
-
-	for(int i=0; i<tree->GetEntries(); i++) {
-		if(++count == percent5) {//Show progress every 5%
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		break1_flag = false;
-		break3_flag = false;
-		break4_flag = false;
-
-		tree->GetEntry(i);
-
-		if(break1->KE >= threshold && (IsRing1(break1->theta, break1->phi) || IsRing2(break1->theta, break1->phi) || IsQQQ(break1->theta, break1->phi))) {
-			b1_count++;
-			break1_flag = true;
-			MyFill(table, "detected_break1_theta_ke","detected breakup1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/deg2rad,300,0.0,30.0,break1->KE);
-		}
-
-
-		if(break3->KE >= threshold && (IsRing1(break3->theta, break3->phi) || IsRing2(break3->theta, break3->phi) || IsQQQ(break3->theta, break3->phi))) {
-			b3_count++;
-			break3_flag = true;
-			MyFill(table, "detected_break3_theta_ke","detected breakup3 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break3->theta/deg2rad,300,0.0,30.0,break3->KE);
-		}
-
-		if(break4->KE >= threshold && (IsRing1(break4->theta, break4->phi) || IsRing2(break4->theta, break4->phi) || IsQQQ(break4->theta, break4->phi))) {
-			b4_count++;
-			break4_flag = true;
-			MyFill(table, "detected_break4_theta_ke","detected breakup4 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break4->theta/deg2rad,300,0.0,30.0,break4->KE);
-		}
-
-		if(break1_flag && break3_flag && break4_flag) {
-			b1b3b4_count++;
-			b1b3_count++;
-			b1b4_count++;
-			b3b4_count++;
-			MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/deg2rad,180,0.0,180.0,break4->theta/deg2rad);
-			MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break3->KE,300,0.0,30.0,break4->KE);
-			MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/deg2rad,360,0.0,360.0,break4->phi/deg2rad);
-		} else if(break1_flag && break3_flag) {
-			b1b3_count++;
-		} else if(break1_flag && break4_flag) {
-			b1b4_count++;
-		} else if(break3_flag && break4_flag) {
-			b3b4_count++;
-			MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/deg2rad,180,0.0,180.0,break4->theta/deg2rad);
-			MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break3->KE,300,0.0,30.0,break4->KE);
-			MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/deg2rad,360,0.0,360.0,break4->phi/deg2rad);
-		}
-	}
-
-	double b1eff = ((double) b1_count)/nevents;
-	double b3eff = ((double) b3_count)/nevents;
-	double b4eff = ((double) b4_count)/nevents;
-	double b1b3eff = b1b3_count/nevents;
-	double b1b4eff = b1b4_count/nevents;
-	double b3b4eff = b3b4_count/nevents;
-	double b1b3b4eff = b1b3b4_count/nevents;
-	TParameter<double> break1_eff("Breakup1 Efficiency", b1eff);
-	TParameter<double> break3_eff("Breakup3 Efficiency", b3eff);
-	TParameter<double> break4_eff("Breakup4 Efficiency", b4eff);
-	TParameter<double> b1b3_eff("Break1 Coincident with Break3", b1b3eff);
-	TParameter<double> b1b4_eff("Break1 Coincident with Break4", b1b4eff);
-	TParameter<double> b3b4_eff("Break3 Coincident with Break4", b3b4eff);
-	TParameter<double> b1b3b4_eff("All Breakups Coincident", b1b3b4eff);
-
-	input->cd();
-	break1_eff.Write();
-	break3_eff.Write();
-	break4_eff.Write();
-	b1b3_eff.Write();
-	b1b4_eff.Write();
-	b3b4_eff.Write();
-	b1b3b4_eff.Write();
-	table->Write();
-	input->Close();
 }
\ No newline at end of file
diff --git a/src/Detectors/DetectorEfficiency.cpp b/src/Detectors/DetectorEfficiency.cpp
new file mode 100644
index 0000000..5d84933
--- /dev/null
+++ b/src/Detectors/DetectorEfficiency.cpp
@@ -0,0 +1,33 @@
+#include "SabreEfficiency.h"
+#include "AnasenEfficiency.h"
+#include <iostream>
+#include <string>
+
+int main(int argc, char** argv) {
+
+	if(argc != 4) {
+		std::cerr<<"Incorrect number of commandline arguments! Returning."<<std::endl;
+		return 1;
+	}
+
+	std::string inputname = argv[1];
+	std::string outputname = argv[2];
+	std::string statsname = argv[3];
+
+	/*
+	SabreEfficiency sabre;
+	std::string mapfile = "./etc/DeadChannels.txt";
+	sabre.SetDeadChannelMap(mapfile);
+	sabre.CalculateEfficiency(inputname, outputname, statsname);
+	//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;
+	//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt");
+	*/
+
+	
+ 	AnasenEfficiency anasen;
+ 	anasen.CalculateEfficiency(inputname, outputname, statsname);
+ 	//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
+ 	//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_test.root");
+	
+	return 0;
+}
\ No newline at end of file
diff --git a/src/Detectors/SabreEfficiency.cpp b/src/Detectors/SabreEfficiency.cpp
index 4272d6c..6c9ff3d 100644
--- a/src/Detectors/SabreEfficiency.cpp
+++ b/src/Detectors/SabreEfficiency.cpp
@@ -1,12 +1,9 @@
 #include "SabreEfficiency.h"
+#include "MaskFile.h"
 #include <fstream>
 #include <iostream>
-#include <TFile.h>
-#include <TTree.h>
-#include <TParameter.h>
-#include <TGraph.h>
-#include <TGraph2D.h>
-#include <TCanvas.h>
+#include <iomanip>
+
 
 SabreEfficiency::SabreEfficiency() : 
 	DetectorEfficiency(), deadlayer(DEADLAYER_THIN), sabre_eloss(SABRE_THICKNESS)
@@ -28,9 +25,9 @@ SabreEfficiency::SabreEfficiency() :
 
 SabreEfficiency::~SabreEfficiency() {}
 
-void SabreEfficiency::CalculateEfficiency(const std::string& file) {
+void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) {
 	std::cout<<"----------SABRE Efficiency Calculation----------"<<std::endl;
-	std::cout<<"Loading in output from kinematics simulation: "<<file<<std::endl;
+	std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
 	std::cout<<"Running efficiency calculation..."<<std::endl;
 
 	if(!dmap.IsValid()) {
@@ -40,35 +37,92 @@ void SabreEfficiency::CalculateEfficiency(const std::string& file) {
 		std::cout<<"---------------------------------------------"<<std::endl;
 	}
 
-	switch(m_rxn_type) {
-		case Mask::Kinematics::ONESTEP_DECAY:
-		{
-			RunDecay(file);
+
+	Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
+	Mask::MaskFile output(outputname, Mask::MaskFile::FileType::read);
+	std::ofstream stats(statsname);
+	stats<<std::setprecision(5);
+
+	Mask::MaskFileHeader header = input.ReadHeader();
+	output.WriteHeader(header.rxn_type, header.nsamples);
+	stats<<"Efficiency statistics for data from "<<inputname<<" using the SPS-SABRE geometry"<<std::endl;
+	stats<<"Given in order of target=0, projectile=1, ejectile=2, residual=3, .... etc."<<std::endl;
+
+	Mask::MaskFileData data;
+	Mask::Nucleus nucleus;
+
+	std::vector<int> counts;
+	switch(header.rxn_type) {
+		case 0:
+			counts.resize(3, 0);
 			break;
-		}
-		case Mask::Kinematics::ONESTEP_RXN:
-		{
-			Run1Step(file);
+		case 1:
+			counts.resize(4, 0);
 			break;
-		}
-		case Mask::Kinematics::TWOSTEP:
-		{
-			Run2Step(file);
+		case 2:
+			counts.resize(6, 0);
 			break;
-		}
-		case Mask::Kinematics::THREESTEP:
-		{
-			Run3Step(file);
+		case 3:
+			counts.resize(8, 0);
 			break;
+		default:
+		{
+			std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<header.rxn_type<<"). Quiting..."<<std::endl;
+			input.Close();
+			output.Close();
+			stats.close();
+			return;
 		}
 	}
+
+	int percent5 = header.nsamples*0.05;
+	int count = 0;
+	int npercent = 0;
+
+	while(true) {
+		if(++count == percent5) {//Show progress every 5%
+			npercent++;
+			count = 0;
+			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
+		}
+
+		data = input.ReadData();
+		if(data.eof)
+			break;
+
+		for(unsigned int i=0; i<data.Z.size(); i++) {
+			nucleus.SetIsotope(data.Z[i], data.A[i]);
+			nucleus.SetVectorSpherical(data.theta[i], data.phi[i], data.p[i], data.E[i]);
+			auto result = IsSabre(nucleus);
+			if(result.first) {
+				data.detect_flag[i] = true;
+				counts[i]++;
+			} else if(data.detect_flag[i] == true) {
+				data.detect_flag[i] = false;
+			}
+		}
+
+		output.WriteData(data);
+	}
+
+	input.Close();
+	output.Close();
+
+	stats<<std::setw(10)<<"Index"<<"|"<<std::setw(10)<<"Efficiency"<<std::endl;
+	stats<<"---------------------"<<std::endl;
+	for(unsigned int i=0; i<counts.size(); i++) {
+		stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+	}
+	stats.close();
+
 	std::cout<<std::endl;
 	std::cout<<"Complete."<<std::endl;
 	std::cout<<"---------------------------------------------"<<std::endl;
 }
 
 void SabreEfficiency::DrawDetectorSystem(const std::string& filename) {
-	TFile* output = TFile::Open(filename.c_str(), "RECREATE");
+	std::ofstream output(filename);
 
 	std::vector<double> ringxs, ringys, ringzs;
     std::vector<double> wedgexs, wedgeys, wedgezs;
@@ -95,27 +149,16 @@ void SabreEfficiency::DrawDetectorSystem(const std::string& filename) {
  		}
  	}
 
- 	TGraph2D* gr = new TGraph2D(ringxs.size(), &(ringxs[0]), &(ringys[0]), &(ringzs[0]));
-	gr->SetName("ring_detector_edges");
-	gr->SetTitle("SABRE Detector; x(m); y(m); z(m)");
-	gr->SetMarkerStyle(1);
+ 	output<<"SABRE Geometry File -- Coordinates for Detectors"<<std::endl;
+	output<<"Edges: x y z"<<std::endl;
+	for(unsigned int i=0; i<ringxs.size(); i++) {
+		output<<ringxs[i]<<" "<<ringys[i]<<" "<<ringzs[i]<<std::endl;
+	}
+	for(unsigned int i=0; i<wedgexs.size(); i++) {
+		output<<wedgexs[i]<<" "<<wedgeys[i]<<" "<<wedgezs[i]<<std::endl;
+	}
 
-	TGraph2D* gw = new TGraph2D(wedgexs.size(), &(wedgexs[0]), &(wedgeys[0]), &(wedgezs[0]));
-	gw->SetName("wedge_detector_edges");
-	gw->SetTitle("SABRE Detector Wedges; x(m); y(m); z(m)");
-	gw->SetMarkerStyle(1);
-
-	TCanvas* canvas = new TCanvas();
-	canvas->SetName("detector_system");
-	canvas->cd();
-	gr->Draw("AP");
-	gw->Draw("same P");
-
-	canvas->Write();
-	gr->Write();
-	gw->Write();
-
-	output->Close();
+ 	output.close();
 }
 
 double SabreEfficiency::RunConsistencyCheck() {
@@ -141,263 +184,25 @@ double SabreEfficiency::RunConsistencyCheck() {
 }
 
 /*Returns if detected, as well as total energy deposited in SABRE*/
-std::pair<bool,double> SabreEfficiency::IsSabre(Mask::NucData* nucleus) {
-	if(nucleus->KE <= ENERGY_THRESHOLD) {
+std::pair<bool,double> SabreEfficiency::IsSabre(Mask::Nucleus& nucleus) {
+	if(nucleus.GetKE() <= ENERGY_THRESHOLD) {
 		return std::make_pair(false, 0.0);
 	}
 
 	Mask::Vec3 coords;
 	double thetaIncident, eloss, e_deposited;
 	for(int i=0; i<5; i++) {
-		auto chan = detectors[i].GetTrajectoryRingWedge(nucleus->theta, nucleus->phi);
+		auto chan = detectors[i].GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
 		if(chan.first != -1 && chan.second != -1) {
 			if(dmap.IsDead(i, chan.first, 0) || dmap.IsDead(i, chan.second, 1)) break; //dead channel check
-			coords = detectors[i].GetTrajectoryCoordinates(nucleus->theta, nucleus->phi);
+			coords = detectors[i].GetTrajectoryCoordinates(nucleus.GetTheta(), nucleus.GetPhi());
 			thetaIncident = std::acos(coords.Dot(detectors[i].GetNormTilted())/(coords.GetR()));
-			eloss = deadlayer.getEnergyLossTotal(nucleus->Z, nucleus->A, nucleus->KE, M_PI - thetaIncident);
-			if((nucleus->KE - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check
-			e_deposited = sabre_eloss.getEnergyLossTotal(nucleus->Z, nucleus->A, nucleus->KE - eloss, M_PI - thetaIncident);
+			eloss = deadlayer.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), M_PI - thetaIncident);
+			if((nucleus.GetKE() - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check
+			e_deposited = sabre_eloss.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE() - eloss, M_PI - thetaIncident);
 			return std::make_pair(true, e_deposited);
 		}
 	}
 
 	return std::make_pair(false,0.0);
 }
-
-void SabreEfficiency::RunDecay(const std::string& filename) {
-	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
-	TTree* tree = (TTree*) input->Get("DataTree");
-
-	THashTable* table = new THashTable();
-
-	Mask::NucData* eject = new Mask::NucData();
-	Mask::NucData* resid = new Mask::NucData();
-
-
-	tree->SetBranchAddress("ejectile", &eject);
-	tree->SetBranchAddress("residual", &resid);
-
-	double nevents = tree->GetEntries();
-
-	//Progress tracking
-	int percent5 = nevents*0.05;
-	int count = 0;
-	int npercent = 0;
-
-	std::pair<bool, double> eject_result, resid_result;
-	int detected_eject=0, detected_resid=0;
-	for(int i=0; i<tree->GetEntries(); i++) {
-		if(++count == percent5) {//Show progress every 5%
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		tree->GetEntry(i);
-
-		eject_result = IsSabre(eject);
-		resid_result = IsSabre(resid);
-
-		if(eject_result.first) {
-			detected_eject++;
-			MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/DEG2RAD,300,0.0,30.0,eject->KE);
-			MyFill(table, "detected_eject_theta_edep","detected ejectile theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,eject->theta/DEG2RAD,300,0.0,30.0,eject_result.second);
-		}
-
-		if(resid_result.first) {
-			detected_resid++;
-			MyFill(table, "detected_resid_theta_ke","detected residual theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,resid->theta/DEG2RAD,300,0.0,30.0,resid->KE);
-			MyFill(table, "detected_resid_theta_edep","detected residual theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,resid->theta/DEG2RAD,300,0.0,30.0,resid_result.second);
-		}
-	}
-
-	double ejecteff = ((double) detected_eject)/nevents;
-	double resideff = ((double) detected_resid)/nevents;
-	TParameter<double> eject_eff("Ejectile Efficiency", ejecteff);
-	TParameter<double> resid_eff("Residual Efficiency", resideff);
-
-	input->cd();
-	eject_eff.Write();
-	resid_eff.Write();
-	table->Write();
-	input->Close();
-
-}
-
-void SabreEfficiency::Run1Step(const std::string& filename) {
-	std::cout<<"SabreEfficiency::Run1Step Not Implemented!"<<std::endl;
-	return;
-}
-
-void SabreEfficiency::Run2Step(const std::string& filename) {
-
-	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
-	TTree* tree = (TTree*) input->Get("DataTree");
-
-	Mask::NucData* break1 = new Mask::NucData();
-	Mask::NucData* break2 = new Mask::NucData();
-
-	THashTable* table = new THashTable();
-
-	tree->SetBranchAddress("breakup1", &break1);
-	tree->SetBranchAddress("breakup2", &break2);
-
-	double nevents = tree->GetEntries();
-
-	//Progress tracking
-	int percent5 = nevents*0.05;
-	int count = 0;
-	int npercent = 0;
-
-	double costheta_cm =0;
-
-	std::pair<bool, double> break1_result, break2_result;
-	int detected_b1=0, detected_b2=0;
-
-	for(int i=0; i<tree->GetEntries(); i++) {
-		if(++count == percent5) {//Show progress every 5%
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		tree->GetEntry(i);
-
-		break1_result = IsSabre(break1);
-		break2_result = IsSabre(break2);
-
-		if(break1_result.first) {
-			detected_b1++;
-			costheta_cm = std::cos(break1->theta_cm);
-			MyFill(table,"detected_break1_cm_theta","cos(#theta_{CM})",20,-1.0,1.0,costheta_cm);
-			MyFill(table, "detected_break1_theta_ke","detected break1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1->KE);
-			MyFill(table, "detected_break1_theta_edep","detected break1 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1_result.second);
-		}
-
-		if(break2_result.first) {
-			detected_b2++;
-			MyFill(table, "detected_break2_theta_ke","detected break2 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break2->theta/DEG2RAD,300,0.0,30.0,break2->KE);
-			MyFill(table, "detected_break2_theta_edep","detected break2 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break2->theta/DEG2RAD,300,0.0,30.0,break2_result.second);
-		
-		}
-	}
-
-	double b1eff = ((double) detected_b1)/nevents;
-	double b2eff = ((double) detected_b2)/nevents;
-	TParameter<double> break1_eff("Breakup1 Efficiency", b1eff);
-	TParameter<double> break2_eff("Breakup2 Efficiency", b2eff);
-
-	input->cd();
-	table->Write();
-	break1_eff.Write();
-	break2_eff.Write();
-	input->Close();
-
-}
-
-void SabreEfficiency::Run3Step(const std::string& filename) {
-	TFile* input = TFile::Open(filename.c_str(), "UPDATE");
-	TTree* tree = (TTree*) input->Get("DataTree");
-
-	THashTable* table = new THashTable();
-
-	Mask::NucData* break1 = new Mask::NucData();
-	Mask::NucData* break3 = new Mask::NucData();
-	Mask::NucData* break4 = new Mask::NucData();
-
-
-	tree->SetBranchAddress("breakup1", &break1);
-	tree->SetBranchAddress("breakup3", &break3);
-	tree->SetBranchAddress("breakup4", &break4);
-
-	double nevents = tree->GetEntries();
-
-	//Progress tracking
-	int percent5 = nevents*0.05;
-	int count = 0;
-	int npercent = 0;
-
-	std::pair<bool, double> break1_result, break3_result, break4_result;
-	int detected_b1=0, detected_b3=0, detected_b4=0;
-	int b1b3_count=0, b1b4_count=0, b3b4_count=0, b1b3b4_count=0;
-
-
-	for(int i=0; i<tree->GetEntries(); i++) {
-		if(++count == percent5) {//Show progress every 5%
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		break1_result = IsSabre(break1);
-		break3_result = IsSabre(break3);
-		break4_result = IsSabre(break4);
-
-
-		tree->GetEntry(i);
-
-		if(break1_result.first) {
-			detected_b1++;
-			MyFill(table, "detected_break1_theta_ke","detected break1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1->KE);
-			MyFill(table, "detected_break1_theta_edep","detected break1 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1_result.second);
-		}
-
-		if(break3_result.first) {
-			detected_b3++;
-			MyFill(table, "detected_break3_theta_ke","detected break3 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break3->theta/DEG2RAD,300,0.0,30.0,break3->KE);
-			MyFill(table, "detected_break3_theta_edep","detected break3 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break3->theta/DEG2RAD,300,0.0,30.0,break3_result.second);
-		}
-
-		if(break4_result.first) {
-			detected_b1++;
-			MyFill(table, "detected_break4_theta_ke","detected break4 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break4->theta/DEG2RAD,300,0.0,30.0,break4->KE);
-			MyFill(table, "detected_break4_theta_edep","detected break4 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break4->theta/DEG2RAD,300,0.0,30.0,break4_result.second);
-		}
-
-		if(break1_result.first && break3_result.first && break4_result.first) {
-			b1b3b4_count++;
-			b1b3_count++;
-			b1b4_count++;
-			b3b4_count++;
-			MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/DEG2RAD,180,0,180,break4->theta/DEG2RAD);
-			MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0.0,15.0,break3->KE,150,0.0,15.0,break4->KE);
-			MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/DEG2RAD,360,0.0,360.0,break4->phi/DEG2RAD);
-		} else if(break1_result.first && break3_result.first) {
-			b1b3_count++;
-		} else if(break1_result.first && break4_result.first) {
-			b1b4_count++;
-		} else if(break3_result.first && break4_result.first) {
-			b3b4_count++;
-			MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/DEG2RAD,180,0,180,break4->theta/DEG2RAD);
-			MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0.0,15.0,break3->KE,150,0.0,15.0,break4->KE);
-			MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/DEG2RAD,360,0.0,360.0,break4->phi/DEG2RAD);
-		}
-	}
-
-	double b1eff = ((double) detected_b1)/nevents;
-	double b3eff = ((double) detected_b3)/nevents;
-	double b4eff = ((double) detected_b4)/nevents;
-	double b1b3eff = b1b3_count/nevents;
-	double b1b4eff = b1b4_count/nevents;
-	double b3b4eff = b3b4_count/nevents;
-	double b1b3b4eff = b1b3b4_count/nevents;
-	TParameter<double> break1_eff("Breakup1 Efficiency", b1eff);
-	TParameter<double> break3_eff("Breakup3 Efficiency", b3eff);
-	TParameter<double> break4_eff("Breakup4 Efficiency", b4eff);
-	TParameter<double> b1b3_eff("Break1 Coincident with Break3", b1b3eff);
-	TParameter<double> b1b4_eff("Break1 Coincident with Break4", b1b4eff);
-	TParameter<double> b3b4_eff("Break3 Coincident with Break4", b3b4eff);
-	TParameter<double> b1b3b4_eff("All Breakups Coincident", b1b3b4eff);
-
-	input->cd();
-	break1_eff.Write();
-	break3_eff.Write();
-	break4_eff.Write();
-	b1b3_eff.Write();
-	b1b4_eff.Write();
-	b3b4_eff.Write();
-	b1b3b4_eff.Write();
-	table->Write();
-
-	input->Close();
-}
diff --git a/src/Detectors/StripDetector.cpp b/src/Detectors/StripDetector.cpp
index 80ad21b..b08a346 100644
--- a/src/Detectors/StripDetector.cpp
+++ b/src/Detectors/StripDetector.cpp
@@ -13,7 +13,7 @@
 */
 
 StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) :
-	m_random(nullptr)
+	m_random(nullptr), m_uniform_fraction(0.0, 1.0)
 {
 
 	num_strips = ns;
@@ -88,7 +88,7 @@ Mask::Vec3 StripDetector::GetHitCoordinates(int front_stripch, double front_stri
 	double y;
 	//If we have a random number generator, randomize y position within pixel. Otherwise take halfway.
 	if(m_random) {
-		y = -total_width/2.0 + (front_stripch + m_random->Uniform(0.0, 1.0))*front_strip_width;
+		y = -total_width/2.0 + (front_stripch + m_uniform_fraction(*m_random))*front_strip_width;
 	} else {
 		y = -total_width/2.0 + (front_stripch+0.5)*front_strip_width;
 	}
diff --git a/src/EnergyLoss.cpp b/src/EnergyLoss.cpp
index 266d942..c332ecf 100644
--- a/src/EnergyLoss.cpp
+++ b/src/EnergyLoss.cpp
@@ -46,7 +46,7 @@ double EnergyLoss::GetEnergyLoss(int zp, int ap, double e_initial, double thickn
   if( ZP != zp) { 
     ZP = zp;
     AP = ap;
-    MP = MASS.FindMass(ZP, AP)*MEV2U;
+    MP = MassLookup::GetInstance()->FindMass(ZP, AP)*MEV2U;
   }
 
   double e_final = e_initial;
@@ -98,7 +98,7 @@ double EnergyLoss::GetReverseEnergyLoss(int zp, int ap, double e_final, double t
   if( ZP != zp) { 
     ZP = zp;
     AP = ap;
-    MP = MASS.FindMass(ZP, AP)*MEV2U;
+    MP = MassLookup::GetInstance()->FindMass(ZP, AP)*MEV2U;
   }
 
   double e_initial = e_final;
diff --git a/src/Kinematics.cpp b/src/Kinematics.cpp
index 50d25a6..55b737e 100644
--- a/src/Kinematics.cpp
+++ b/src/Kinematics.cpp
@@ -5,471 +5,355 @@
 
 namespace Mask {
 
-Kinematics::Kinematics() :
-	sys(nullptr), save_tree_flag(false), do_plotter_flag(false), global_generator(new TRandom3(0))
-{
-	std::cout<<"----------GWM Kinematics Simulation----------"<<std::endl;
-}
-
-Kinematics::~Kinematics() {
-	delete global_generator;
-	if(sys) delete sys;
-}
-
-bool Kinematics::LoadConfig(const char* filename) {
-	std::cout<<"Loading configuration in "<<filename<<"..."<<std::endl;
-
-	std::ifstream input(filename);
-	if(!input.is_open()) {
-		std::cerr<<"Unable to load configuration in "<<filename<<", check that it exists"<<std::endl;
-		return false;
+	Kinematics::Kinematics() :
+		sys(nullptr)
+	{
+		std::cout<<"----------GWM Kinematics Simulation----------"<<std::endl;
+	
+		std::random_device rd;
+		global_generator = new std::mt19937(rd());
 	}
-
-	std::string junk;
-	getline(input, junk);
-	input>>junk>>m_outfile_name;
-	input>>junk>>junk;
-	if(junk == "yes") save_tree_flag = true;
-	input>>junk>>junk;
-	if(junk == "yes") do_plotter_flag = true;
-
-	std::vector<int> avec, zvec, svec;
-	int z, a, s;
-	getline(input, junk);
-	getline(input, junk);
-	input>>junk>>m_rxn_type;
-	getline(input, junk);
-	getline(input, junk);
-	switch(m_rxn_type) {
-		case 0:
-		{
-			sys = new DecaySystem();
-			m_rxn_type = ONESTEP_DECAY;
-			for(int i=0; i<2; i++) {
-				input>>z>>a;
-				avec.push_back(a);
-				zvec.push_back(z);
-			}
-			break;
-		}
-		case 1:
-		{
-			sys = new OneStepSystem();
-			m_rxn_type = ONESTEP_RXN;
-			for(int i=0; i<3; i++) {
-				input>>z>>a;
-				avec.push_back(a);
-				zvec.push_back(z);
-			}
-			break;
-		}
-		case 2:
-		{
-			sys = new TwoStepSystem();
-			m_rxn_type = TWOSTEP;
-			for(int i=0; i<4; i++) {
-				input>>z>>a;
-				avec.push_back(a);
-				zvec.push_back(z);
-			}
-			break;
-		}
-		case 3:
-		{
-			sys = new ThreeStepSystem();
-			m_rxn_type = THREESTEP;
-			for(int i=0; i<5; i++) {
-				input>>z>>a;
-				avec.push_back(a);
-				zvec.push_back(z);
-			}
-			break;
-		}
-		default:
+	
+	Kinematics::~Kinematics() {
+		delete global_generator;
+		if(sys) delete sys;
+	}
+	
+	bool Kinematics::LoadConfig(const std::string& filename) {
+		std::cout<<"Loading configuration in "<<filename<<"..."<<std::endl;
+	
+		std::ifstream input(filename);
+		if(!input.is_open()) {
+			std::cerr<<"Unable to load configuration in "<<filename<<", check that it exists"<<std::endl;
 			return false;
-	}
-	sys->SetNuclei(zvec, avec);
-
-	int nlayers;
-	double thickness;
-	getline(input, junk);
-	getline(input, junk);
-	input>>junk>>junk;
-	input>>junk>>nlayers;
-	for(int i=0; i<nlayers; i++) {
-		input>>junk>>junk>>thickness;
+		}
+	
+		std::string junk;
+		getline(input, junk);
+		input>>junk>>m_outfile_name;
+		input>>junk>>junk;
+		input>>junk>>junk;
+	
+		std::vector<int> avec, zvec, svec;
+		int z, a, s;
 		getline(input, junk);
 		getline(input, junk);
-		avec.clear(); zvec.clear(); svec.clear();
-		while(input>>z) {
-			if(z == 0) break;
-			input>>a>>s;
-			zvec.push_back(z); avec.push_back(a); svec.push_back(s);
+		input>>junk>>m_rxn_type;
+		getline(input, junk);
+		getline(input, junk);
+		switch(m_rxn_type) {
+			case 0:
+			{
+				sys = new DecaySystem();
+				m_rxn_type = ONESTEP_DECAY;
+				for(int i=0; i<2; i++) {
+					input>>z>>a;
+					avec.push_back(a);
+					zvec.push_back(z);
+				}
+				break;
+			}
+			case 1:
+			{
+				sys = new OneStepSystem();
+				m_rxn_type = ONESTEP_RXN;
+				for(int i=0; i<3; i++) {
+					input>>z>>a;
+					avec.push_back(a);
+					zvec.push_back(z);
+				}
+				break;
+			}
+			case 2:
+			{
+				sys = new TwoStepSystem();
+				m_rxn_type = TWOSTEP;
+				for(int i=0; i<4; i++) {
+					input>>z>>a;
+					avec.push_back(a);
+					zvec.push_back(z);
+				}
+				break;
+			}
+			case 3:
+			{
+				sys = new ThreeStepSystem();
+				m_rxn_type = THREESTEP;
+				for(int i=0; i<5; i++) {
+					input>>z>>a;
+					avec.push_back(a);
+					zvec.push_back(z);
+				}
+				break;
+			}
+			default:
+				return false;
 		}
-		sys->AddTargetLayer(zvec, avec, svec, thickness);
-		input>>junk;
-	}
-	std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
-
-	double par1, par2;
-	std::string dfile1, dfile2;
-	getline(input, junk);
-	getline(input, junk);
-
-	input>>junk>>m_nsamples;
-	input>>junk>>par1>>junk>>par2;
-	sys->SetBeamDistro(par1, par2);
-	input>>junk>>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;
-	sys->SetPhi1Range(par1, par2);
-	input>>junk>>par1>>junk>>par2;
-	sys->SetExcitationDistro(par1, par2);
-	input>>junk>>dfile1;
-	input>>junk>>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<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
-
-	return true;
-}
-
-bool Kinematics::SaveConfig(const char* filename) { return true; }
-
-NucData Kinematics::ConvertNucleus(const Nucleus& nuc) {
-	NucData datum;
-	datum.E = nuc.GetE();
-	datum.KE = nuc.GetKE();
-	datum.p = nuc.GetP();
-	datum.theta = nuc.GetTheta();
-	datum.theta_cm = nuc.GetThetaCM();
-	datum.phi = nuc.GetPhi();
-	datum.Ex = nuc.GetExcitationEnergy();
-	datum.Z = nuc.GetZ();
-	datum.A = nuc.GetA();
-	return datum;
-}
-
-void Kinematics::Run() {
-	std::cout<<"Running simulation..."<<std::endl;
-	switch(m_rxn_type) {
-		case ONESTEP_DECAY:
-		{
-			RunOneStepDecay();
-			break;
-		}
-		case ONESTEP_RXN:
-		{
-			RunOneStepRxn();
-			break;
-		}
-		case TWOSTEP:
-		{
-			RunTwoStep();
-			break;
-		}
-		case THREESTEP:
-		{
-			RunThreeStep();
-			break;
-		}
-	}
-	std::cout<<std::endl;
-	std::cout<<"Complete."<<std::endl;
-	std::cout<<"---------------------------------------------"<<std::endl;
-}
-
-void Kinematics::RunOneStepRxn() {
-	OneStepSystem* this_sys = dynamic_cast<OneStepSystem*>(sys);
-	if(this_sys == nullptr) {
-		return;
-	}
-
+		sys->SetNuclei(zvec, avec);
 	
-	TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
-	TTree* tree;
-	NucData targ, proj, eject, residual;
-	if(save_tree_flag) {
-		tree = new TTree("DataTree","DataTree");
-		tree->Branch("target", &targ);
-		tree->Branch("projectile", &proj);
-		tree->Branch("ejectile", &eject);
-		tree->Branch("residual", &residual);
+		int nlayers;
+		double thickness;
+		getline(input, junk);
+		getline(input, junk);
+		input>>junk>>junk;
+		input>>junk>>nlayers;
+		for(int i=0; i<nlayers; i++) {
+			input>>junk>>junk>>thickness;
+			getline(input, junk);
+			getline(input, junk);
+			avec.clear(); zvec.clear(); svec.clear();
+			while(input>>z) {
+				if(z == 0) break;
+				input>>a>>s;
+				zvec.push_back(z); avec.push_back(a); svec.push_back(s);
+			}
+			sys->AddTargetLayer(zvec, avec, svec, thickness);
+			input>>junk;
+		}
+		std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
+	
+		double par1, par2;
+		std::string dfile1, dfile2;
+		getline(input, junk);
+		getline(input, junk);
+	
+		input>>junk>>m_nsamples;
+		input>>junk>>par1>>junk>>par2;
+		sys->SetBeamDistro(par1, par2);
+		input>>junk>>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;
+		sys->SetPhi1Range(par1, par2);
+		input>>junk>>par1>>junk>>par2;
+		sys->SetExcitationDistro(par1, par2);
+		input>>junk>>dfile1;
+		input>>junk>>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<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
+	
+		return true;
 	}
 	
-
-	//For progress tracking
-	int percent5 = 0.05*m_nsamples;
-	int count = 0;
-	int npercent = 0;
-
-	std::vector<Nucleus> data;
-	data.resize(4);
-	for(int i=0; i<m_nsamples; i++) {
-		if(++count == percent5) {//Show update every 5 percent
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		this_sys->RunSystem();
-		
-		if(save_tree_flag) {
-			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(this_sys->GetTarget(), "targ");
-			plotman.FillData(this_sys->GetProjectile(), "proj");
-			plotman.FillData(this_sys->GetEjectile(), "eject");
-			plotman.FillData(this_sys->GetResidual(), "resid");
-		}
-
-	}
-
+	bool Kinematics::SaveConfig(const std::string& filename) { return true; }
 	
-	output->cd();
-	if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite);
-	if(do_plotter_flag) {
-		plotman.GetTable()->Write();
-		plotman.ClearTable();
-	}
-	output->Close();
-	
-}
-
-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;
-	NucData targ, eject, residual;
-	if(save_tree_flag) {
-		tree = new TTree("DataTree","DataTree");
-		tree->Branch("target", &targ);
-		tree->Branch("ejectile", &eject);
-		tree->Branch("residual", &residual);
-	}
-
-	//For progress tracking
-	int percent5 = 0.05*m_nsamples;
-	int count = 0;
-	int npercent = 0;
-
-	for(int i=0; i<m_nsamples; i++) {
-		if(++count == percent5) {//Show update every 5 percent
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
+	void Kinematics::Run() {
+		std::cout<<"Running simulation..."<<std::endl;
+		switch(m_rxn_type) {
+			case ONESTEP_DECAY:
+			{
+				RunOneStepDecay();
+				break;
+			}
+			case ONESTEP_RXN:
+			{
+				RunOneStepRxn();
+				break;
+			}
+			case TWOSTEP:
+			{
+				RunTwoStep();
+				break;
+			}
+			case THREESTEP:
+			{
+				RunThreeStep();
+				break;
+			}
 		}
-
-		this_sys->RunSystem();
-		if(save_tree_flag) {
-			targ = ConvertNucleus(this_sys->GetTarget());
-			eject = ConvertNucleus(this_sys->GetEjectile());
-			residual = ConvertNucleus(this_sys->GetResidual());
-			tree->Fill();
-		}
-		if(do_plotter_flag) {
-			plotman.FillData(this_sys->GetTarget(), "targ");
-			plotman.FillData(this_sys->GetEjectile(), "eject");
-			plotman.FillData(this_sys->GetResidual(), "resid");
-		}
-
-	}
-
-	output->cd();
-	if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite);
-	if(do_plotter_flag) {
-		plotman.GetTable()->Write();
-		plotman.ClearTable();
-	}
-	output->Close();
-}
-
-void Kinematics::RunTwoStep() {
-
-	TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(sys);
-	if(this_sys == nullptr) {
-		return;
-	}
-
-	
-	TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
-	TTree* tree;
-	NucData targ, proj, eject, residual, break1, break2;
-	if(save_tree_flag) {
-		tree = new TTree("DataTree","DataTree");
-		tree->Branch("target", &targ);
-		tree->Branch("projectile", &proj);
-		tree->Branch("ejectile", &eject);
-		tree->Branch("residual", &residual);
-		tree->Branch("breakup1", &break1);
-		tree->Branch("breakup2", &break2);
+		std::cout<<std::endl;
+		std::cout<<"Complete."<<std::endl;
+		std::cout<<"---------------------------------------------"<<std::endl;
 	}
 	
-
-
-	//For progress tracking
-	int percent5 = 0.05*m_nsamples;
-	int count = 0;
-	int npercent = 0;
-
-	for(int i=0; i<m_nsamples; i++) {
-		if(++count == percent5) {//Show update every 5 percent
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		this_sys->RunSystem();
-		if(save_tree_flag) {
-			targ = ConvertNucleus(this_sys->GetTarget());
-			proj = ConvertNucleus(this_sys->GetProjectile());
-			eject = ConvertNucleus(this_sys->GetEjectile());
-			residual = ConvertNucleus(this_sys->GetResidual());
-			break1 = ConvertNucleus(this_sys->GetBreakup1());
-			break2 = ConvertNucleus(this_sys->GetBreakup2());
-			tree->Fill();
-		}
-		if(do_plotter_flag) {
-			plotman.FillData(this_sys->GetTarget(), "targ");
-			plotman.FillData(this_sys->GetProjectile(), "proj");
-			plotman.FillData(this_sys->GetEjectile(), "eject");
-			plotman.FillData(this_sys->GetResidual(), "resid");
-			plotman.FillData(this_sys->GetBreakup1(), "break1");
-			plotman.FillData(this_sys->GetBreakup2(), "break2");
+	void Kinematics::RunOneStepRxn() {
+		OneStepSystem* this_sys = dynamic_cast<OneStepSystem*>(sys);
+		if(this_sys == nullptr) {
+			return;
 		}
 		
+		MaskFile output(m_outfile_name, MaskFile::FileType::write);
+		std::vector<Nucleus> data;
+		data.resize(4);
+		output.WriteHeader(m_rxn_type, m_nsamples);
+		
+	
+		//For progress tracking
+		int percent5 = 0.05*m_nsamples;
+		int count = 0;
+		int npercent = 0;
+	
+		for(int i=0; i<m_nsamples; i++) {
+			if(++count == percent5) {//Show update every 5 percent
+				npercent++;
+				count = 0;
+				std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
+			}
+	
+			this_sys->RunSystem();
+			data[0] = this_sys->GetTarget();
+			data[1] = this_sys->GetProjectile();
+			data[2] = this_sys->GetEjectile();
+			data[3] = this_sys->GetResidual();
+			output.WriteData(data);
+			
+		}
+	
+		output.Close();
+	}
+	
+	void Kinematics::RunOneStepDecay() {
+		DecaySystem* this_sys = dynamic_cast<DecaySystem*>(sys);
+		if(this_sys == nullptr) {
+			return;
+		}
+	
+		MaskFile output(m_outfile_name, MaskFile::FileType::write);
+		std::vector<Nucleus> data;
+		data.resize(3);
+		output.WriteHeader(m_rxn_type, m_nsamples);
+	
+		//For progress tracking
+		int percent5 = 0.05*m_nsamples;
+		int count = 0;
+		int npercent = 0;
+	
+		for(int i=0; i<m_nsamples; i++) {
+			if(++count == percent5) {//Show update every 5 percent
+				npercent++;
+				count = 0;
+				std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
+			}
+	
+			this_sys->RunSystem();
+			data[0] = this_sys->GetTarget();
+			data[1] = this_sys->GetEjectile();
+			data[2] = this_sys->GetResidual();
+			output.WriteData(data);
+		}
+	
+		output.Close();
+	}
+	
+	void Kinematics::RunTwoStep() {
+	
+		TwoStepSystem* this_sys = dynamic_cast<TwoStepSystem*>(sys);
+		if(this_sys == nullptr) {
+			return;
+		}
+	
+		MaskFile output(m_outfile_name, MaskFile::FileType::write);
+		std::vector<Nucleus> data;
+		data.resize(6);
+		output.WriteHeader(m_rxn_type, m_nsamples);
+	
+		//For progress tracking
+		int percent5 = 0.05*m_nsamples;
+		int count = 0;
+		int npercent = 0;
+	
+		for(int i=0; i<m_nsamples; i++) {
+			if(++count == percent5) {//Show update every 5 percent
+				npercent++;
+				count = 0;
+				std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
+			}
+	
+			this_sys->RunSystem();
+			data[0] = this_sys->GetTarget();
+			data[1] = this_sys->GetProjectile();
+			data[2] = this_sys->GetEjectile();
+			data[3] = this_sys->GetResidual();
+			data[4] = this_sys->GetBreakup1();
+			data[5] = this_sys->GetBreakup2();
+			output.WriteData(data);
+		}
+	
+		output.Close();
+	}
+	
+	void Kinematics::RunThreeStep() {
+	
+		ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(sys);
+		if(this_sys == nullptr) {
+			return;
+		}
+	
+		MaskFile output(m_outfile_name, MaskFile::FileType::write);
+		std::vector<Nucleus> data;
+		data.resize(8);
+		output.WriteHeader(m_rxn_type, m_nsamples);
+	
+		//For progress updating
+		int percent5 = 0.05*m_nsamples;
+		int count = 0;
+		int npercent = 0;
+	
+		for(int i=0; i<m_nsamples; i++) {
+			if(++count == percent5) {//Show update every 5 percent
+				npercent++;
+				count = 0;
+				std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
+			}
+	
+			this_sys->RunSystem();
+			data[0] = this_sys->GetTarget();
+			data[1] = this_sys->GetProjectile();
+			data[2] = this_sys->GetEjectile();
+			data[3] = this_sys->GetResidual();
+			data[4] = this_sys->GetBreakup1();
+			data[5] = this_sys->GetBreakup2();
+			data[6] = this_sys->GetBreakup3();
+			data[7] = this_sys->GetBreakup4();
+			output.WriteData(data);
+		}
+	
+		output.Close();
 	}
 
-	
-	output->cd();
-	if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite);
-	if(do_plotter_flag) {
-		plotman.GetTable()->Write();
-		plotman.ClearTable();
-	}
-	output->Close();
-	
 }
-
-void Kinematics::RunThreeStep() {
-
-	ThreeStepSystem* this_sys = dynamic_cast<ThreeStepSystem*>(sys);
-	if(this_sys == nullptr) {
-		return;
-	}
-
-	TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE");
-	TTree* tree;
-	NucData targ, proj, eject, residual, break1, break2, break3, break4;
-	if(save_tree_flag) {
-		tree = new TTree("DataTree","DataTree");
-		tree->Branch("target", &targ);
-		tree->Branch("projectile", &proj);
-		tree->Branch("ejectile", &eject);
-		tree->Branch("residual", &residual);
-		tree->Branch("breakup1", &break1);
-		tree->Branch("breakup2", &break2);
-		tree->Branch("breakup3", &break3);
-		tree->Branch("breakup4", &break4);
-	}
-	
-	//For progress updating
-	int percent5 = 0.05*m_nsamples;
-	int count = 0;
-	int npercent = 0;
-
-	for(int i=0; i<m_nsamples; i++) {
-		if(++count == percent5) {//Show update every 5 percent
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		this_sys->RunSystem();
-		if(save_tree_flag) {
-			targ = ConvertNucleus(this_sys->GetTarget());
-			proj = ConvertNucleus(this_sys->GetProjectile());
-			eject = ConvertNucleus(this_sys->GetEjectile());
-			residual = ConvertNucleus(this_sys->GetResidual());
-			break1 = ConvertNucleus(this_sys->GetBreakup1());
-			break2 = ConvertNucleus(this_sys->GetBreakup2());
-			break3 = ConvertNucleus(this_sys->GetBreakup3());
-			break4 = ConvertNucleus(this_sys->GetBreakup4());
-			tree->Fill();
-		}
-		if(do_plotter_flag) {
-			plotman.FillData(this_sys->GetTarget(), "targ");
-			plotman.FillData(this_sys->GetProjectile(), "proj");
-			plotman.FillData(this_sys->GetEjectile(), "eject");
-			plotman.FillData(this_sys->GetResidual(), "resid");
-			plotman.FillData(this_sys->GetBreakup1(), "break1");
-			plotman.FillData(this_sys->GetBreakup2(), "break2");
-			plotman.FillData(this_sys->GetBreakup3(), "break3");
-			plotman.FillData(this_sys->GetBreakup4(), "break4");
-		}
-	}
-
-	output->cd();
-	if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite);
-	if(do_plotter_flag) {
-		plotman.GetTable()->Write();
-		plotman.ClearTable();
-	}
-	output->Close();
-}
-
-};
diff --git a/src/MaskFile.cpp b/src/MaskFile.cpp
index af9b279..5af094a 100644
--- a/src/MaskFile.cpp
+++ b/src/MaskFile.cpp
@@ -3,177 +3,336 @@
 #include <iomanip>
 #include <sstream>
 
+/*
+
+	FORMAT
+
+	HEADER (contains rxntype & nuclei numbers & beam kinetic energy) (64 bits, 8 bytes)
+		NSAMPLES(32bit) RXNTYPE(32bit)
+	end HEADER
+
+	There are NSAMPLES * (number of saved nuclei) data in the file. The number of nuclei saved is related to the 
+	RXNTYPE. All nuclei (target, projectile, ejectile, residual, break1, etc...) are saved. A datum is as follows:
+
+	DATUM (contains kinematic data for a nucleus) (384 bits, 48 bytes)
+		Z(32bit) A(32bit) DETECTFLAG(32bit) E(64bit) KE(64bit) P(64bit) THETA(64bit) PHI(64bit)
+	end DATUM
+
+
+*/
+
 namespace Mask {
 
-MaskFile::MaskFile() : 
-	file_type(MaskFile::FileType::none), filename(""), buffer_position(0), count(0), file()
-{
-}
-
-MaskFile::MaskFile(const std::string& name, MaskFile::FileType type) :
-	file_type(type), filename(name), buffer_position(0), count(0), file()
-{
-	Open(filename, type);
-}
-
-bool MaskFile::Open(const std::string& name, MaskFile::FileType type) {
-	if(IsOpen()) {
-		std::cerr<<"Attempted to open file that is already open!"<<std::endl;
-		return false;
+	MaskFile::MaskFile() : 
+		file_type(MaskFile::FileType::none), filename(""), buffer_position(0), buffer_end(0), data_size(0), buffersize_bytes(0), file()
+	{
 	}
-
-	if(type == FileType::read) {
-		file.open(name, std::ios::in);
-		buffer_position = 0;
-		if(IsOpen()) ReadHeader();
-	} else if (type == FileType::write) {
-		file.open(name, std::ios::out | std::ios::trunc);
-		file<<std::setprecision(precision);
-		file<<std::scientific;
-	} else if (type == FileType::append) {
-		file.open(name, std::ios::out | std::ios::app);
-		file<<std::setprecision(precision);
-		file<<std::scientific;
-	} else {
-		std::cerr<<"Invalid FileType at MaskFile::Open()"<<std::endl;
+	
+	MaskFile::MaskFile(const std::string& name, MaskFile::FileType type) :
+		file_type(type), filename(name), buffer_position(0), buffer_end(0), data_size(0), buffersize_bytes(0), file()
+	{
+		Open(filename, type);
+	}
+	
+	bool MaskFile::Open(const std::string& name, MaskFile::FileType type) {
+		if(IsOpen()) {
+			std::cerr<<"Attempted to open file that is already open!"<<std::endl;
+			return false;
+		}
+	
+		if(type == FileType::read) {
+			file.open(name, std::ios::in);
+			buffer_position = 0;
+		} else if (type == FileType::write) {
+			file.open(name, std::ios::out | std::ios::trunc);
+		} else if (type == FileType::append) {
+			file.open(name, std::ios::out | std::ios::app);
+		} else {
+			std::cerr<<"Invalid FileType at MaskFile::Open()"<<std::endl;
+			return IsOpen();
+		}
+	
+	
 		return IsOpen();
 	}
-
-	count = 0;
-
-	return IsOpen();
-}
-
-void MaskFile::Close() {
-	if(data_buffer.size() > 0) {
-		FlushBuffer();
-	}
-	file.close();
-}
-
-void MaskFile::WriteHeader(std::vector<Nucleus>& nuclei) {
-	file<<std::setw(width)<<"ZT"<<","
-		<<std::setw(width)<<"AT"<<","
-		<<std::setw(width)<<"ZP"<<","
-		<<std::setw(width)<<"AP"<<","
-		<<std::setw(width)<<"ZE"<<","
-		<<std::setw(width)<<"AE"<<","
-		<<std::setw(width)<<"ZR"<<","
-		<<std::setw(width)<<"AR"<<","
-		<<std::setw(width)<<"ZB1"<<","
-		<<std::setw(width)<<"AB1"<<","
-		<<std::setw(width)<<"ZB2"<<","
-		<<std::setw(width)<<"AB2"<<","
-		<<std::setw(width)<<"ZB3"<<","
-		<<std::setw(width)<<"AB3"<<","
-		<<std::setw(width)<<"ZB4"<<","
-		<<std::setw(width)<<"AB4"<<","
-		<<std::setw(width)<<"KEP"<<","
-		<<std::endl;
-
-	file<<std::setw(width)<<nuclei[0].GetZ()<<","
-		<<std::setw(width)<<nuclei[0].GetA()<<","
-		<<std::setw(width)<<nuclei[1].GetZ()<<","
-		<<std::setw(width)<<nuclei[1].GetA()<<","
-		<<std::setw(width)<<nuclei[2].GetZ()<<","
-		<<std::setw(width)<<nuclei[2].GetA()<<","
-		<<std::setw(width)<<nuclei[3].GetZ()<<","
-		<<std::setw(width)<<nuclei[3].GetA()<<",";
-	if(nuclei.size() == 6) {
-		file<<std::setw(width)<<nuclei[4].GetZ()<<","
-			<<std::setw(width)<<nuclei[4].GetA()<<","
-			<<std::setw(width)<<nuclei[5].GetZ()<<","
-			<<std::setw(width)<<nuclei[5].GetA()<<",";
-	} else if(nuclei.size() == 8) {
-		file<<std::setw(width)<<nuclei[4].GetZ()<<","
-			<<std::setw(width)<<nuclei[4].GetA()<<","
-			<<std::setw(width)<<nuclei[5].GetZ()<<","
-			<<std::setw(width)<<nuclei[5].GetA()<<","
-			<<std::setw(width)<<nuclei[6].GetZ()<<","
-			<<std::setw(width)<<nuclei[6].GetA()<<","
-			<<std::setw(width)<<nuclei[7].GetZ()<<","
-			<<std::setw(width)<<nuclei[7].GetA()<<",";
-	}
-	file<<std::setw(width)<<nuclei[1].GetKE()<<",";
-	file<<std::endl;
-
-	file<<std::setw(width)<<"Event"<<","
-		<<std::setw(width)<<"EE"<<","<<std::setw(width)<<"KEE"<<","<<std::setw(width)<<"PE"<<","<<std::setw(width)<<"ThetaE"<<","<<std::setw(width)<<"PhiE"<<","<<std::setw(width)<<"ThetaCME"<<","
-		<<std::setw(width)<<"ER"<<","<<std::setw(width)<<"KER"<<","<<std::setw(width)<<"PR"<<","<<std::setw(width)<<"ThetaR"<<","<<std::setw(width)<<"PhiR"<<","<<std::setw(width)<<"ThetaCMR"<<","
-		<<std::setw(width)<<"EB1"<<","<<std::setw(width)<<"KEB1"<<","<<std::setw(width)<<"PB1"<<","<<std::setw(width)<<"ThetaB1"<<","<<std::setw(width)<<"PhiB1"<<","<<std::setw(width)<<"ThetaCMB1"<<","
-		<<std::setw(width)<<"EB2"<<","<<std::setw(width)<<"KEB2"<<","<<std::setw(width)<<"PB2"<<","<<std::setw(width)<<"ThetaB2"<<","<<std::setw(width)<<"PhiB2"<<","<<std::setw(width)<<"ThetaCMB2"<<","
-		<<std::setw(width)<<"EB3"<<","<<std::setw(width)<<"KEB3"<<","<<std::setw(width)<<"PB3"<<","<<std::setw(width)<<"ThetaB3"<<","<<std::setw(width)<<"PhiB3"<<","<<std::setw(width)<<"ThetaCMB3"<<","
-		<<std::setw(width)<<"EB4"<<","<<std::setw(width)<<"KEB4"<<","<<std::setw(width)<<"PB4"<<","<<std::setw(width)<<"ThetaB4"<<","<<std::setw(width)<<"PhiB4"<<","<<std::setw(width)<<"ThetaCMB4"<<","
-		<<std::endl;
-}
-
-void MaskFile::ReadHeader() {
-	std::string junk;
-	std::getline(file, junk);
-}
-
-void MaskFile::WriteData(std::vector<Nucleus>& data) {
-	if(count == 0 && data_buffer.size() == 0) {
-		WriteHeader(data);
+	
+	void MaskFile::Close() {
+		//Final flush (if necessary)
+		if(buffer_position > 0 && buffer_position < buffersize_bytes && (file_type == FileType::write || file_type == FileType::append)) {
+			file.write(data_buffer.data(), buffer_position);
+		}
+		file.close();
 	}
 	
-	data_buffer.push_back(data);
-	if(data_buffer.size() == buffersize) {
-		FlushBuffer();
-	}
-}
-
-void MaskFile::FlushBuffer() {
-	for(auto& event : data_buffer) {
-		file<<std::setw(width)<<count<<",";
-		for(unsigned int i=2; i< event.size(); i++) {
-			Nucleus& nuc = event[i];
-			file<<std::setw(width)<<nuc.GetE()<<","
-				<<std::setw(width)<<nuc.GetKE()<<","
-				<<std::setw(width)<<nuc.GetP()<<","
-				<<std::setw(width)<<nuc.GetTheta()<<","
-				<<std::setw(width)<<nuc.GetPhi()<<","
-				<<std::setw(width)<<nuc.GetThetaCM()<<",";
+	void MaskFile::WriteHeader(int rxn_type, int nsamples) {
+		//size of a datum = data nuclei per event * (# doubles per nucleus * sizeof double + # ints per nucleus * sizeof int + sizeof bool)
+		if(rxn_type == 0) {
+			m_rxn_type = 0;
+			data_size = 3 * ( 5 * double_size + 2 * int_size + bool_size); 
 		}
-		file<<std::endl;
-		count++;
-	}
-	data_buffer.clear();
-}
-
-std::vector<Nucleus> MaskFile::ReadData() {
-	if(buffer_position == data_buffer.size()) {
-		FillBuffer();
-		buffer_position = 0;
-	}
-	std::vector<Nucleus> data = data_buffer[buffer_position];
-	buffer_position++;
+		else if (rxn_type == 1) {
+			m_rxn_type = 1;
+			data_size = 4 * ( 5 * double_size + 2 * int_size + bool_size);
+		}
+		else if (rxn_type == 2) {
+			m_rxn_type = 2;
+			data_size = 6 * ( 5 * double_size + 2 * int_size + bool_size);
+		}
+		else if (rxn_type == 3) {
+			m_rxn_type = 3;
+			data_size = 8 * ( 5 * double_size + 2 * int_size + bool_size);
+		} else {
+			std::cerr<<"Invalid number of nuclei at MaskFile::WriteHeader! Returning"<<std::endl;
+			return;
+		}
+		buffersize_bytes = buffersize * data_size; //buffersize_bytes = # of events * size of an event
 	
-	return data;
-}
-
-void MaskFile::FillBuffer() {
-	std::vector<Nucleus> data;
-	std::string line;
-	std::stringstream linebuf;
-	std::string junk;
-	int Z, A;
-	double E, p, theta, phi, tcm;
-	data_buffer.clear();
-
-	while(data_buffer.size() <= buffersize) {
-		if(!std::getline(file, line)) break;
-		linebuf.str(line);
-		while(linebuf>>Z) {
-			linebuf>>A;
-			linebuf>>E>>junk>>junk>>p>>theta>>phi>>tcm;
-			data.emplace_back(Z, A);
-			data[data.size()-1].SetVectorSpherical(theta, phi, p, E);
-			data[data.size()-1].SetThetaCM(tcm);
-		}
-		data_buffer.push_back(data);
-		data.clear();
+		data_buffer.resize(buffersize_bytes);
+	
+		file.write((char*) &nsamples, int_size);
+		file.write((char*) &m_rxn_type, int_size);
+	}
+	
+	MaskFileHeader MaskFile::ReadHeader() {
+		MaskFileHeader header;
+		std::vector<char> temp_buffer(4);
+		file.read(temp_buffer.data(), 4);
+		header.nsamples = *(int*)(&temp_buffer[0]);
+		file.read(temp_buffer.data(), 4);
+		m_rxn_type = *(int*)(&temp_buffer[0]);
+	
+		//size of a datum = data nuclei per event * (# doubles per nucleus * sizeof double + # ints per nucleus * sizeof int + sizeof bool)
+		if(m_rxn_type == 0) {
+			data_size = 3 * ( 5 * double_size + 2 * int_size + bool_size);
+		}
+		else if (m_rxn_type == 1) {
+			data_size = 4 * ( 5 * double_size + 2 * int_size + bool_size);
+		}
+		else if (m_rxn_type == 2) {
+			data_size = 6 * ( 5 * double_size + 2 * int_size + bool_size);
+		}
+		else if (m_rxn_type == 3) {
+			data_size = 8 * ( 5 * double_size + 2 * int_size + bool_size);
+		} else {
+			std::cerr<<"Unexpected reaction type at MaskFile::ReadHeader (rxn type = "<<m_rxn_type<<")! Returning"<<std::endl;
+			return header;
+		}
+		buffersize_bytes = buffersize * data_size;//buffersize_bytes = size of a datum * # of events
+	
+		header.rxn_type = m_rxn_type;
+	
+		data_buffer.resize(buffersize_bytes);
+	
+		return header;
+	}
+	
+	void MaskFile::WriteData(std::vector<Nucleus>& data) {
+	
+		char* data_pointer;
+		double datum;
+		int number;
+		bool flag;
+		std::size_t j;
+		for(unsigned int i=0; i<data.size(); i++) {
+			number = data[i].GetZ();
+			data_pointer = (char*) &number;
+			for(j=0; j<int_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+
+			number = data[i].GetA();
+			data_pointer = (char*) &number;
+			for(j=0; j<int_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+
+			flag = data[i].IsDetected();
+			data_pointer = (char*) &flag;
+			for(j=0; j<bool_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+
+			datum = data[i].GetE();
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+
+			datum = data[i].GetKE();
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+	
+			datum = data[i].GetP();
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+	
+			datum = data[i].GetTheta();
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+	
+			datum = data[i].GetPhi();
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+		}
+	
+		//Flush the buffer when it is full, and reset the position.
+		if(buffer_position == buffersize_bytes) {
+			file.write(data_buffer.data(), data_buffer.size());
+			buffer_position = 0;
+		}
+	}
+
+	void MaskFile::WriteData(MaskFileData& data) {
+	
+		char* data_pointer;
+		double datum;
+		int number;
+		bool flag;
+		std::size_t j;
+		for(unsigned int i=0; i<data.Z.size(); i++) {
+			number = data.Z[i];
+			data_pointer = (char*) &number;
+			for(j=0; j<int_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+
+			number = data.A[i];
+			data_pointer = (char*) &number;
+			for(j=0; j<int_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+
+			flag = data.detect_flag[i];
+			data_pointer = (char*) &flag;
+			for(j=0; j<bool_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+
+			datum = data.E[i];
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+
+			datum = data.KE[i];
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+	
+			datum = data.p[i];
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+	
+			datum = data.theta[i];
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+	
+			datum = data.phi[i];
+			data_pointer = (char*) &datum;
+			for(j=0; j<double_size; j++) {
+				data_buffer[buffer_position] = *(data_pointer + j);
+				buffer_position++;
+			}
+		}
+	
+		//Flush the buffer when it is full, and reset the position.
+		if(buffer_position == buffersize_bytes) {
+			file.write(data_buffer.data(), data_buffer.size());
+			buffer_position = 0;
+		}
+	}
+	
+	/*
+		Read data from the buffer and submit it to the client side as a MaskFileData struct.
+		When file reaches the end of the file (no more data to read), an empty MaskFileData with 
+		eof == true is sent out signaling that the file is finished.
+
+		Should be used like
+
+		Mask::MaskFile input(file, Mask::MaskFile::FileType::read);
+		Mask::MaskFileHeader header = input.ReadHeader();
+		Mask::MaskFileData data;
+		while(true) {
+			data = input.ReadData();
+			if(data.eof) break;
+
+			Do some stuff...
+
+		}
+		input.Close();
+
+		Good luck
+	*/
+	MaskFileData MaskFile::ReadData() {
+	
+		MaskFileData data;
+	
+		//Fill the buffer when needed, reset the positon, and find the end
+		if(buffer_position == data_buffer.size() || buffer_position == buffer_end) {
+			file.read(data_buffer.data(), buffersize_bytes);
+			buffer_position = 0;
+			buffer_end = file.gcount();
+			if(buffer_end == 0 && file.eof()) {
+				data.eof = true;
+				return data;
+			}
+		}
+	
+		unsigned int local_end = buffer_position + data_size;
+		if(local_end > buffer_end) {
+			std::cerr<<"Attempting to read past end of file at MaskFile::ReadData! Returning empty"<<std::endl;
+			data.eof = true;
+			return data;
+		}
+
+		while(buffer_position < local_end) {
+			data.Z.push_back(*(int*)(&data_buffer[buffer_position]));
+			buffer_position += int_size;
+			data.A.push_back(*(int*)(&data_buffer[buffer_position]));
+			buffer_position += int_size;
+			data.detect_flag.push_back(*(bool*)(&data_buffer[buffer_position]));
+			buffer_position += bool_size;
+			data.E.push_back(*(double*)(&data_buffer[buffer_position]));
+			buffer_position += double_size;
+			data.KE.push_back(*(double*)(&data_buffer[buffer_position]));
+			buffer_position += double_size;
+			data.p.push_back(*(double*)(&data_buffer[buffer_position]));
+			buffer_position += double_size;
+			data.theta.push_back(*(double*)(&data_buffer[buffer_position]));
+			buffer_position += double_size;
+			data.phi.push_back(*(double*)(&data_buffer[buffer_position]));
+			buffer_position += double_size;
+		}
+	
+		return data;
 	}
-}
 
 };
\ No newline at end of file
diff --git a/src/MassLookup.cpp b/src/MassLookup.cpp
index 0de274a..14e3171 100644
--- a/src/MassLookup.cpp
+++ b/src/MassLookup.cpp
@@ -16,61 +16,47 @@ Written by G.W. McCann Aug. 2020
   Read in AMDC mass file, preformated to remove excess info. Here assumes that by default
   the file is in a local directory etc/
 */
+
+MassLookup* MassLookup::s_instance = nullptr;
+
 MassLookup::MassLookup() {
-  std::ifstream massfile("./etc/mass.txt");
-  if(massfile.is_open()) {
-    std::string junk, A, element;
-    int Z;
-    double atomicMassBig, atomicMassSmall, isotopicMass;
-    getline(massfile,junk);
-    getline(massfile,junk);
-    while(massfile>>junk) {
-      massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall;
-      isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev;
-      std::string key = "("+std::to_string(Z)+","+A+")";
-      massTable[key] = isotopicMass;
-      elementTable[Z] = element;
-    }
-  } else {
-    throw MassFileException();
-  }
+	std::ifstream massfile("./etc/mass.txt");
+	if(massfile.is_open()) {
+		std::string junk, A, element;
+		int Z;
+		double atomicMassBig, atomicMassSmall, isotopicMass;
+		getline(massfile,junk);
+		getline(massfile,junk);
+		while(massfile>>junk) {
+			massfile>>Z>>A>>element>>atomicMassBig>>atomicMassSmall;
+			isotopicMass = (atomicMassBig + atomicMassSmall*1e-6 - Z*electron_mass)*u_to_mev;
+			std::string key = "("+std::to_string(Z)+","+A+")";
+			massTable[key] = isotopicMass;
+			elementTable[Z] = element;
+		}
+	} else {
+	  throw MassFileException();
+	}
 }
 
 MassLookup::~MassLookup() {}
 
 //Returns nuclear mass in MeV
 double MassLookup::FindMass(int Z, int A) {
-  std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")";
-  auto data = massTable.find(key);
-  if(data == massTable.end()) {
-    throw MassException();
-  }
-  return data->second;
-  /*try {
-    double mass = massTable.at(key);
-    return mass;
-  } catch (std::out_of_range& oor) {
-    std::cerr<<"Mass of "<<key<<" (Z,A) not found in Mass Table! Returning 1"<<std::endl;
-    return 1;
-  }*/
+  	std::string key = "("+std::to_string(Z)+","+std::to_string(A)+")";
+  	auto data = massTable.find(key);
+  	if(data == massTable.end()) {
+  	  	throw MassException();
+  	}
+  	return data->second;
 }
 
 //returns element symbol
 std::string MassLookup::FindSymbol(int Z, int A) {
-  auto data = elementTable.find(Z);
-  if(data == elementTable.end()) {
-    throw MassException();
-  }
-  std::string fullsymbol = std::to_string(A) + data->second;
-  return fullsymbol;
-
-  /*try {
-    std::string element = elementTable.at(Z);
-    std::string fullsymbol = std::to_string(A) + element;
-    return fullsymbol;
-  } catch (std::out_of_range& oor) {
-    std::cerr<<"Atomic number: "<<Z<<" not found in Element Table! Returning void."<<std::endl;
-    std::string element = "void";
-    return element;
-  }*/
+  	auto data = elementTable.find(Z);
+  	if(data == elementTable.end()) {
+  	  	throw MassException();
+  	}
+  	std::string fullsymbol = std::to_string(A) + data->second;
+  	return fullsymbol;
 }
diff --git a/src/Nucleus.cpp b/src/Nucleus.cpp
index f98d241..d5d621f 100644
--- a/src/Nucleus.cpp
+++ b/src/Nucleus.cpp
@@ -12,23 +12,23 @@
 namespace Mask {
 
 Nucleus::Nucleus () :
-	Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_theta_cm(0), m_symbol("")
+	Vec4(), m_z(0), m_a(0), m_gs_mass(0), m_theta_cm(0), m_symbol(""), m_detectFlag(false)
 {
 }
 
 Nucleus::Nucleus(int Z, int A) :
-	Vec4(), m_z(Z), m_a(A), m_theta_cm(0)
+	Vec4(), m_z(Z), m_a(A), m_theta_cm(0), m_detectFlag(false)
 {
-	m_gs_mass = MASS.FindMass(Z, A);
-	m_symbol = MASS.FindSymbol(Z, A);
+	m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A);
+	m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A);
 	SetVectorCartesian(0,0,0,m_gs_mass); //by defualt a nucleus has mass given by the g.s.
 }
 
 Nucleus::Nucleus(int Z, int A, double px, double py, double pz, double E) :
 	Vec4(px, py, pz, E), m_z(Z), m_a(A)
 {
-	m_gs_mass = MASS.FindMass(Z, A);
-	m_symbol = MASS.FindSymbol(Z, A);
+	m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A);
+	m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A);
 }
 
 Nucleus::~Nucleus() {}
@@ -38,8 +38,8 @@ bool Nucleus::SetIsotope(int Z, int A) {
 	
 	m_z = Z;
 	m_a = A;
-	m_gs_mass = MASS.FindMass(Z, A);
-	m_symbol = MASS.FindSymbol(Z, A);
+	m_gs_mass = MassLookup::GetInstance()->FindMass(Z, A);
+	m_symbol = MassLookup::GetInstance()->FindSymbol(Z, A);
 	SetVectorCartesian(0,0,0,m_gs_mass);
 	return true;
 }
diff --git a/src/OneStepSystem.cpp b/src/OneStepSystem.cpp
index d9b2b69..ad2e257 100644
--- a/src/OneStepSystem.cpp
+++ b/src/OneStepSystem.cpp
@@ -2,76 +2,76 @@
 
 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;
+	OneStepSystem::OneStepSystem() :
+		ReactionSystem()
+	{
 	}
-
-	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();
+	OneStepSystem::OneStepSystem(std::vector<int>& z, std::vector<int>& a) :
+		ReactionSystem()
+	{
+		SetNuclei(z, a);
 	}
-
-	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);
+	OneStepSystem::~OneStepSystem() {}
 	
-		step1.TurnOnResidualEloss();
-		step1.Calculate();
-	} else {
-		return;
+	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 = (*m_beamDist)(*generator);
+			double rxnTheta = std::acos((*m_theta1Range)(*generator));
+			double rxnPhi = (*m_phi1Range)(*generator);
+			double residEx = (*m_exDist)(*generator);
+		
+			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/Plotter.cpp b/src/Plotter.cpp
deleted file mode 100644
index f9a5753..0000000
--- a/src/Plotter.cpp
+++ /dev/null
@@ -1,89 +0,0 @@
-#include "Plotter.h"
-#include <iostream>
-
-namespace Mask {
-
-Plotter::Plotter() :
-	table(new THashTable())
-{
-}
-
-Plotter::~Plotter() {
-	for(unsigned int i=0; i<garbage_collection.size(); i++) {
-		delete garbage_collection[i];
-	}
-	garbage_collection.clear();
-	delete table;
-}
-
-void Plotter::FillData(const Nucleus& nuc, const std::string& modifier) {
-	std::string sym = nuc.GetIsotopicSymbol();
-	std::string ke_vs_th_name = sym + modifier + "_ke_vs_theta";
-	std::string ke_vs_th_title = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
-	std::string ke_vs_ph_name = sym + modifier + "_ke_vs_phi";
-	std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
-	std::string ex_name = sym + modifier + "_ex";
-	std::string ex_title = ex_name + ";E_{ex} (MeV);counts";
-	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.GetTheta()*rad2deg, nuc.GetKE(), 2);
-	MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, nuc.GetKE(), 4);
-	MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
-	MyFill(angdist_name.c_str(), angdist_title.c_str(),100,-1.0,1.0,std::cos(nuc.GetThetaCM()));
-	
-}
-
-void Plotter::MyFill(const char* name, const char* title, int bins, float min, float max, double val) {
-	TH1F* h = (TH1F*) table->FindObject(name);
-	if(h) {
-		h->Fill(val);
-	} else {
-		h = new TH1F(name, title, bins, min, max);
-		h->Fill(val);
-		table->Add(h);
-	}
-}
-
-void Plotter::MyFill(const char* name, const char* title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy) {
-	TH2F* h = (TH2F*) table->FindObject(name);
-	if(h) {
-		h->Fill(valx, valy);
-	} else {
-		h = new TH2F(name, title, binsx, minx, maxx, binsy, miny, maxy);
-		h->Fill(valx, valy);
-		table->Add(h);
-	}
-}
-
-void Plotter::MyFill(const char* name, const char* title, double valx, double valy, int color) {
-	for(auto& g : graphs) {
-		if(g.name == name) {
-			g.xvec.push_back(valx);
-			g.yvec.push_back(valy);
-			return;
-		}
-	}
-
-	GraphData new_g;
-	new_g.name = name;
-	new_g.title = title;
-	new_g.xvec.push_back(valx);
-	new_g.yvec.push_back(valy);
-	new_g.color = color;
-
-	graphs.push_back(new_g);
-}
-
-void Plotter::GenerateGraphs() {
-	for(auto& g : graphs) {
-		TGraph* graph = new TGraph(g.xvec.size(), &(g.xvec[0]), &(g.yvec[0]));
-		graph->SetName(g.name.c_str());
-		graph->SetTitle(g.title.c_str());
-		graph->SetMarkerColor(g.color);
-		table->Add(graph);
-		garbage_collection.push_back(graph);
-	}
-}
-
-};
\ No newline at end of file
diff --git a/src/Plotters/ROOT/RootPlotter.cpp b/src/Plotters/ROOT/RootPlotter.cpp
new file mode 100644
index 0000000..1996863
--- /dev/null
+++ b/src/Plotters/ROOT/RootPlotter.cpp
@@ -0,0 +1,137 @@
+#include "RootPlotter.h"
+#include "MaskFile.h"
+#include <TFile.h>
+
+#include <iostream>
+
+RootPlotter::RootPlotter() :
+	table(new THashTable())
+{
+}
+
+RootPlotter::~RootPlotter() {}
+
+void RootPlotter::FillData(const Mask::Nucleus& nuc, const std::string& modifier) {
+	std::string sym = nuc.GetIsotopicSymbol();
+	std::string ke_vs_th_name = sym + modifier + "_ke_vs_theta";
+	std::string ke_vs_th_title = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
+	std::string ke_vs_ph_name = sym + modifier + "_ke_vs_phi";
+	std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
+	std::string ex_name = sym + modifier + "_ex";
+	std::string ex_title = ex_name + ";E_{ex} (MeV);counts";
+	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.GetTheta()*rad2deg, nuc.GetKE(), 2);
+	MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, nuc.GetKE(), 4);
+	MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy());
+	MyFill(angdist_name.c_str(), angdist_title.c_str(),100,-1.0,1.0,std::cos(nuc.GetThetaCM()));
+	
+}
+
+void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val) {
+	TH1F* h = (TH1F*) table->FindObject(name.c_str());
+	if(h) {
+		h->Fill(val);
+	} else {
+		h = new TH1F(name.c_str(), title.c_str(), bins, min, max);
+		h->Fill(val);
+		table->Add(h);
+	}
+}
+
+void RootPlotter::MyFill(const std::string& name, const std::string& title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy) {
+	TH2F* h = (TH2F*) table->FindObject(name.c_str());
+	if(h) {
+		h->Fill(valx, valy);
+	} else {
+		h = new TH2F(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy);
+		h->Fill(valx, valy);
+		table->Add(h);
+	}
+}
+
+void RootPlotter::MyFill(const std::string& name, const std::string& title, double valx, double valy, int color) {
+	for(auto& g : graphs) {
+		if(g.name == name) {
+			g.xvec.push_back(valx);
+			g.yvec.push_back(valy);
+			return;
+		}
+	}
+
+	GraphData new_g;
+	new_g.name = name;
+	new_g.title = title;
+	new_g.xvec.push_back(valx);
+	new_g.yvec.push_back(valy);
+	new_g.color = color;
+
+	graphs.push_back(new_g);
+}
+
+void RootPlotter::GenerateGraphs() {
+	for(auto& g : graphs) {
+		TGraph* graph = new TGraph(g.xvec.size(), &(g.xvec[0]), &(g.yvec[0]));
+		graph->SetName(g.name.c_str());
+		graph->SetTitle(g.title.c_str());
+		graph->SetMarkerColor(g.color);
+		table->Add(graph);
+		garbage_collection.push_back(graph);
+	}
+}
+
+int main(int argc, char** argv) {
+	if(argc != 3) {
+		std::cout<<"Unable to run ROOT plotting tool with incorrect number of arguments! Expected 2 args, given: "<<argc<<" Exiting."<<std::endl;
+		return 1;
+	}
+
+	std::string inputname = argv[1];
+	std::string outputname = argv[2];
+
+	Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
+
+	TFile* root_out = TFile::Open(outputname.c_str(), "RECREATE");
+
+	RootPlotter plotter;
+
+	Mask::MaskFileHeader header = input.ReadHeader();
+
+	std::cout<<"File Header -- rxn type: "<<header.rxn_type<<" nsamples: "<<header.nsamples<<std::endl;
+
+	Mask::MaskFileData data;
+	Mask::Nucleus nucleus;
+
+	double flush_frac = 0.05;
+	int count = 0, flush_val = flush_frac*header.nsamples, flush_count = 0;
+	while(true) {
+		if(count == flush_val) {
+			count = 0;
+			flush_count++;
+			std::cout<<"\rPercent of file processed: "<<flush_frac*flush_count*100<<"%"<<std::flush;
+		}
+
+		data = input.ReadData();
+		if(data.eof)
+			break;
+		for(unsigned int i=0; i<data.Z.size(); i++) {
+			nucleus.SetIsotope(data.Z[i], data.A[i]);
+			nucleus.SetVectorSpherical(data.theta[i], data.phi[i], data.p[i], data.E[i]);
+			plotter.FillData(nucleus);
+			if(data.detect_flag[i] == true) {
+				plotter.FillData(nucleus, "detected");
+			}
+		}
+		count++;
+	}
+	std::cout<<std::endl;
+
+	input.Close();
+
+	root_out->cd();
+	plotter.GetTable()->Write();
+	root_out->Close();
+
+	return 0;
+}
\ No newline at end of file
diff --git a/src/ReactionSystem.cpp b/src/ReactionSystem.cpp
index 4034dff..d282832 100644
--- a/src/ReactionSystem.cpp
+++ b/src/ReactionSystem.cpp
@@ -4,20 +4,26 @@
 
 namespace Mask {
 
-ReactionSystem::ReactionSystem() :
-	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("")
-{
-}
+	ReactionSystem::ReactionSystem() :
+		m_beamDist(nullptr), m_theta1Range(nullptr), m_phi1Range(nullptr), m_exDist(nullptr), generator(nullptr), target_set_flag(false),
+		gen_set_flag(false), rxnLayer(0), m_sys_equation("")
+	{
+	}
+	
+	ReactionSystem::~ReactionSystem() {
+		delete m_beamDist;
+		delete m_theta1Range;
+		delete m_phi1Range;
+		delete m_exDist;
+	}
+	
+	void ReactionSystem::SetRandomGenerator(std::mt19937* gen) {
+		generator = gen;
+		gen_set_flag = true;
+	}
+	
+	void ReactionSystem::AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness) {
+		target.AddLayer(zt, at, stoich, thickness);
+	}
 
-ReactionSystem::~ReactionSystem() {}
-
-void ReactionSystem::SetRandomGenerator(TRandom3* gen) {
-	generator = gen;
-	gen_set_flag = true;
-}
-
-void ReactionSystem::AddTargetLayer(std::vector<int>& zt, std::vector<int>& at, std::vector<int>& stoich, double thickness) {
-	target.AddLayer(zt, at, stoich, thickness);
-}
-
-};
\ No newline at end of file
+}
\ No newline at end of file
diff --git a/src/Stopwatch.cpp b/src/Stopwatch.cpp
new file mode 100644
index 0000000..337098e
--- /dev/null
+++ b/src/Stopwatch.cpp
@@ -0,0 +1,31 @@
+/*
+	Stopwatch.cpp
+	Simple class designed to provide timing info on parts of the process.
+	Only for use in development.
+
+	Written by G.W. McCann Oct. 2020
+*/
+#include "Stopwatch.h"
+
+Stopwatch::Stopwatch() {
+	start_time = Clock::now();
+	stop_time = start_time;
+}
+
+Stopwatch::~Stopwatch() {}
+
+void Stopwatch::Start() {
+	start_time = Clock::now();
+}
+
+void Stopwatch::Stop() {
+	stop_time = Clock::now();
+}
+
+double Stopwatch::GetElapsedSeconds() {
+	return std::chrono::duration_cast<std::chrono::duration<double>>(stop_time-start_time).count();
+}
+
+double Stopwatch::GetElapsedMilliseconds() {
+	return std::chrono::duration_cast<std::chrono::duration<double>>(stop_time-start_time).count()*1000.0;
+}
\ No newline at end of file
diff --git a/src/ThreeStepSystem.cpp b/src/ThreeStepSystem.cpp
index 88c5256..83b5113 100644
--- a/src/ThreeStepSystem.cpp
+++ b/src/ThreeStepSystem.cpp
@@ -3,131 +3,129 @@
 
 namespace Mask {
 
-ThreeStepSystem::ThreeStepSystem() :
-	ReactionSystem()
-{
-}
-
-ThreeStepSystem::ThreeStepSystem(std::vector<int>& z, std::vector<int>& a) :
-	ReactionSystem()
-{
-	SetNuclei(z, a);
-}
-
-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;
+	ThreeStepSystem::ThreeStepSystem() :
+		ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
+	{
 	}
-	int zr = z[0] + z[1] - z[2];
-	int zb2 = zr - z[3];
-	int ar = a[0] + a[1] - a[2];
-	int ab2 = ar - a[3];
-
-	step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
-	step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]);
-	step3.SetNuclei(zb2, ab2, 0, 0, z[4], a[4]);
-	SetSystemEquation();
-	return true;
-}
-
-void ThreeStepSystem::LinkTarget() {
-	step1.SetLayeredTarget(&target);
-	step2.SetLayeredTarget(&target);
-	step3.SetLayeredTarget(&target);
-
-	rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
-	if(rxnLayer != -1) {
-		step1.SetRxnLayer(rxnLayer);
-		step2.SetRxnLayer(rxnLayer);
-		step3.SetRxnLayer(rxnLayer);
-		target_set_flag = true;
-	} else {
-		throw ReactionLayerException();
-	}
-}
-
-void ThreeStepSystem::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();
-	m_sys_equation += "-> ";
-	m_sys_equation += step2.GetEjectile().GetIsotopicSymbol();
-	m_sys_equation += " + ";
-	m_sys_equation += step2.GetResidual().GetIsotopicSymbol();
-	m_sys_equation += "-> ";
-	m_sys_equation += step3.GetEjectile().GetIsotopicSymbol();
-	m_sys_equation += " + ";
-	m_sys_equation += step3.GetResidual().GetIsotopicSymbol();
-}
-
-void ThreeStepSystem::RunSystem() {
-	if(!gen_set_flag) return;
 	
-	//Link up the target if it hasn't been done yet
-	if(!target_set_flag) {
-		LinkTarget();
+	ThreeStepSystem::ThreeStepSystem(std::vector<int>& z, std::vector<int>& a) :
+		ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
+	{
+		SetNuclei(z, a);
+	}
+	
+	ThreeStepSystem::~ThreeStepSystem() {}
+	
+	void ThreeStepSystem::SetRandomGenerator(std::mt19937* 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;
+		}
+		int zr = z[0] + z[1] - z[2];
+		int zb2 = zr - z[3];
+		int ar = a[0] + a[1] - a[2];
+		int ab2 = ar - a[3];
+	
+		step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
+		step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]);
+		step3.SetNuclei(zb2, ab2, 0, 0, z[4], a[4]);
+		SetSystemEquation();
+		return true;
+	}
+	
+	void ThreeStepSystem::LinkTarget() {
+		step1.SetLayeredTarget(&target);
+		step2.SetLayeredTarget(&target);
+		step3.SetLayeredTarget(&target);
+	
+		rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
+		if(rxnLayer != -1) {
+			step1.SetRxnLayer(rxnLayer);
+			step2.SetRxnLayer(rxnLayer);
+			step3.SetRxnLayer(rxnLayer);
+			target_set_flag = true;
+		} else {
+			throw ReactionLayerException();
+		}
+	}
+	
+	void ThreeStepSystem::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();
+		m_sys_equation += "-> ";
+		m_sys_equation += step2.GetEjectile().GetIsotopicSymbol();
+		m_sys_equation += " + ";
+		m_sys_equation += step2.GetResidual().GetIsotopicSymbol();
+		m_sys_equation += "-> ";
+		m_sys_equation += step3.GetEjectile().GetIsotopicSymbol();
+		m_sys_equation += " + ";
+		m_sys_equation += step3.GetResidual().GetIsotopicSymbol();
+	}
+	
+	void ThreeStepSystem::RunSystem() {
+		if(!gen_set_flag) return;
+		
+		//Link up the target if it hasn't been done yet
+		if(!target_set_flag) {
+			LinkTarget();
+		}
+	
+		//Sample parameters
+		double bke = (*m_beamDist)(*generator);
+		double rxnTheta = acos((*m_theta1Range)(*generator));
+		double rxnPhi = (*m_phi1Range)(*generator);
+		double decay1costheta = decay1dist.GetRandomCosTheta();
+		double decay1Theta = std::acos(decay1costheta);
+		double decay1Phi = m_phi2Range(*generator);
+		double decay2costheta = decay2dist.GetRandomCosTheta();
+		double decay2Theta = std::acos(decay2costheta);
+		double decay2Phi = m_phi2Range(*generator);
+		double residEx = (*m_exDist)(*generator);
+	
+		step1.SetBeamKE(bke);
+		step1.SetPolarRxnAngle(rxnTheta);
+		step1.SetAzimRxnAngle(rxnPhi);
+		step1.SetExcitation(residEx);
+	
+		step2.SetPolarRxnAngle(decay1Theta);
+		step2.SetAzimRxnAngle(decay1Phi);
+	
+		step3.SetPolarRxnAngle(decay2Theta);
+		step3.SetAzimRxnAngle(decay2Phi);
+	
+		step1.Calculate();
+	
+		step2.SetTarget(step1.GetResidual());
+		if(decay1costheta == -10) {
+			step2.ResetEjectile();
+			step2.ResetResidual();
+			step3.ResetTarget();
+			step3.ResetEjectile();
+			step3.ResetResidual();
+			return;
+		}
+		step2.Calculate();
+	
+		step3.SetTarget(step2.GetResidual());
+		if(decay2costheta == -10) {
+			step3.ResetEjectile();
+			step3.ResetResidual();
+			return;
+		}
+		step3.TurnOnResidualEloss();
+		step3.Calculate();
+	
 	}
 
-	//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 decay1costheta = decay1dist.GetRandomCosTheta();
-	double decay1Theta = std::acos(decay1costheta);
-	double decay1Phi = generator->Uniform(0, M_PI*2.0);
-	double decay2costheta = decay2dist.GetRandomCosTheta();
-	double decay2Theta = std::acos(decay2costheta);
-	double decay2Phi = generator->Uniform(0, M_PI*2.0);
-	double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
-
-	step1.SetBeamKE(bke);
-	step1.SetPolarRxnAngle(rxnTheta);
-	step1.SetAzimRxnAngle(rxnPhi);
-	step1.SetExcitation(residEx);
-
-	step2.SetPolarRxnAngle(decay1Theta);
-	step2.SetAzimRxnAngle(decay1Phi);
-
-	step3.SetPolarRxnAngle(decay2Theta);
-	step3.SetAzimRxnAngle(decay2Phi);
-
-	step1.Calculate();
-
-	step2.SetTarget(step1.GetResidual());
-	if(decay1costheta == -10) {
-		step2.ResetEjectile();
-		step2.ResetResidual();
-		step3.ResetTarget();
-		step3.ResetEjectile();
-		step3.ResetResidual();
-		return;
-	}
-	step2.Calculate();
-
-	step3.SetTarget(step2.GetResidual());
-	if(decay2costheta == -10) {
-		step3.ResetEjectile();
-		step3.ResetResidual();
-		return;
-	}
-	step3.TurnOnResidualEloss();
-	step3.Calculate();
-
-}
-
-};
\ No newline at end of file
+}
\ No newline at end of file
diff --git a/src/TwoStepSystem.cpp b/src/TwoStepSystem.cpp
index c4f69d9..575cbc6 100644
--- a/src/TwoStepSystem.cpp
+++ b/src/TwoStepSystem.cpp
@@ -3,105 +3,105 @@
 
 namespace Mask {
 
-TwoStepSystem::TwoStepSystem() :
-	ReactionSystem()
-{
-}
-
-TwoStepSystem::TwoStepSystem(std::vector<int>& z, std::vector<int>& a) :
-	ReactionSystem()
-{
-	SetNuclei(z, a);
-}
-
-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) {
-		return false;
+	TwoStepSystem::TwoStepSystem() :
+		ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
+	{
 	}
-	int zr = z[0] + z[1] - z[2];
-	int ar = a[0] + a[1] - a[2];
-
-	step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
-	step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]);
-	SetSystemEquation();
-	return true;
-}
-
-void TwoStepSystem::LinkTarget() {
-	step1.SetLayeredTarget(&target);
-	step2.SetLayeredTarget(&target);
-
-	rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
-	if(rxnLayer != -1) {
-		step1.SetRxnLayer(rxnLayer);
-		step2.SetRxnLayer(rxnLayer);
-		target_set_flag = true;
-	} else {
-		throw ReactionLayerException();
-	}
-}
-
-void TwoStepSystem::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();
-	m_sys_equation += "-> ";
-	m_sys_equation += step2.GetEjectile().GetIsotopicSymbol();
-	m_sys_equation += "+";
-	m_sys_equation += step2.GetResidual().GetIsotopicSymbol();
-}
-
-void TwoStepSystem::RunSystem() {
-	if(!gen_set_flag) return;
 	
-	//Link up the target if it hasn't been done yet
-	if(!target_set_flag) {
-		LinkTarget();
+	TwoStepSystem::TwoStepSystem(std::vector<int>& z, std::vector<int>& a) :
+		ReactionSystem(), m_phi2Range(0, 2.0*M_PI)
+	{
+		SetNuclei(z, a);
 	}
-
-	//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 decay1costheta = decay1dist.GetRandomCosTheta();
-	double decay1Theta = std::acos(decay1costheta);
-	double decay1Phi = generator->Uniform(0, M_PI*2.0);
-	double residEx = generator->Gaus(m_exDist.first, m_exDist.second);
-
-	step1.SetBeamKE(bke);
-	step1.SetPolarRxnAngle(rxnTheta);
-	step1.SetAzimRxnAngle(rxnPhi);
-	step1.SetExcitation(residEx);
-
-	step2.SetPolarRxnAngle(decay1Theta);
-	step2.SetAzimRxnAngle(decay1Phi);
 	
-	step1.Calculate();
-
-	step2.SetTarget(step1.GetResidual());
-	if(decay1costheta == -10) {
-		step2.ResetEjectile();
-		step2.ResetResidual();
-		return;
+	TwoStepSystem::~TwoStepSystem() {
+	
+	}
+	
+	void TwoStepSystem::SetRandomGenerator(std::mt19937* 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) {
+			return false;
+		}
+		int zr = z[0] + z[1] - z[2];
+		int ar = a[0] + a[1] - a[2];
+	
+		step1.SetNuclei(z[0], a[0], z[1], a[1], z[2], a[2]);
+		step2.SetNuclei(zr, ar, 0, 0, z[3], a[3]);
+		SetSystemEquation();
+		return true;
+	}
+	
+	void TwoStepSystem::LinkTarget() {
+		step1.SetLayeredTarget(&target);
+		step2.SetLayeredTarget(&target);
+	
+		rxnLayer = target.FindLayerContaining(step1.GetTarget().GetZ(), step1.GetTarget().GetA());
+		if(rxnLayer != -1) {
+			step1.SetRxnLayer(rxnLayer);
+			step2.SetRxnLayer(rxnLayer);
+			target_set_flag = true;
+		} else {
+			throw ReactionLayerException();
+		}
+	}
+	
+	void TwoStepSystem::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();
+		m_sys_equation += "-> ";
+		m_sys_equation += step2.GetEjectile().GetIsotopicSymbol();
+		m_sys_equation += "+";
+		m_sys_equation += step2.GetResidual().GetIsotopicSymbol();
+	}
+	
+	void TwoStepSystem::RunSystem() {
+		if(!gen_set_flag) return;
+		
+		//Link up the target if it hasn't been done yet
+		if(!target_set_flag) {
+			LinkTarget();
+		}
+	
+		//Sample parameters
+		double bke = (*m_beamDist)(*generator);
+		double rxnTheta = acos((*m_theta1Range)(*generator));
+		double rxnPhi = (*m_phi1Range)(*generator);
+		double decay1costheta = decay1dist.GetRandomCosTheta();
+		double decay1Theta = std::acos(decay1costheta);
+		double decay1Phi = m_phi2Range(*generator);
+		double residEx = (*m_beamDist)(*generator);
+	
+		step1.SetBeamKE(bke);
+		step1.SetPolarRxnAngle(rxnTheta);
+		step1.SetAzimRxnAngle(rxnPhi);
+		step1.SetExcitation(residEx);
+	
+		step2.SetPolarRxnAngle(decay1Theta);
+		step2.SetAzimRxnAngle(decay1Phi);
+		
+		step1.Calculate();
+	
+		step2.SetTarget(step1.GetResidual());
+		if(decay1costheta == -10) {
+			step2.ResetEjectile();
+			step2.ResetResidual();
+			return;
+		}
+		step2.TurnOnResidualEloss();
+		step2.Calculate();
+	
+	
 	}
-	step2.TurnOnResidualEloss();
-	step2.Calculate();
 
-
-}
-
-};
\ No newline at end of file
+}
\ No newline at end of file
diff --git a/src/main.cpp b/src/main.cpp
index 197e8af..3fd5446 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,21 +1,20 @@
 #include <iostream>
 #include <string>
+#include "Stopwatch.h"
+#include "MaskFile.h"
 #include "Kinematics.h"
-#include "SabreEfficiency.h"
-#include "AnasenEfficiency.h"
 #include "KinematicsExceptions.h"
 
-#include <TGraph2D.h>
-#include <TCanvas.h>
-
 int main(int argc, char** argv) {
 	if(argc<2) {
 		std::cerr<<"Incorrect number of arguments!"<<std::endl;
 		return 1;
 	}
-	
-	
+	Stopwatch sw;
+
+
 	Mask::Kinematics calculator;
+	sw.Start();
 	try {
 		if(!calculator.LoadConfig(argv[1])) {
 			return 1;
@@ -26,26 +25,26 @@ int main(int argc, char** argv) {
 		std::cerr<<"Terminating process."<<std::endl;
 		return 1;
 	}
+	sw.Stop();
+	std::cout<<"Time elapsed(seconds): "<<sw.GetElapsedSeconds()<<std::endl;
 	
-	
+	Mask::MaskFile input(calculator.GetOutputName(), Mask::MaskFile::FileType::read);
+	Mask::MaskFileData data;
+	Mask::MaskFileHeader header = input.ReadHeader();
 
-	/*
-	SabreEfficiency sabre;
-	std::string mapfile = "./etc/DeadChannels.txt";
-	sabre.SetReactionType(calculator.GetReactionType());
-	sabre.SetDeadChannelMap(mapfile);
-	sabre.CalculateEfficiency(calculator.GetOutputName());
-	//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;
-	//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.root");
-	*/
-
-
- 	AnasenEfficiency anasen;
- 	anasen.SetReactionType(calculator.GetReactionType());
- 	anasen.CalculateEfficiency(calculator.GetOutputName());
- 	//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
- 	//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_test.root");
+	std::cout<<"Header Found -- rxn type: "<<header.rxn_type<<" nsamples: "<<header.nsamples;
+	std::cout<<std::endl;
 
+	int counter=0;
+	while(!data.eof) {
+		data = input.ReadData();
+		for(unsigned int i=0; i<data.E.size(); i++)
+			counter++;
+			//std::cout<<"Data Found -- E: "<<data.E[i]<<" KE: "<<data.KE[i]<<" p: "<<data.p[i]<<" theta: "<<data.theta[i]<<" phi: "<<data.phi[i]<<std::endl;
+	}
+	std::cout<<"events found: "<<counter<<std::endl;
+	input.Close();
+	std::cout<<"File closed."<<std::endl;
 
 	return 0;