diff --git a/include/RootPlotter.h b/include/RootPlotter.h
index 17dc408..31cb931 100644
--- a/include/RootPlotter.h
+++ b/include/RootPlotter.h
@@ -5,6 +5,8 @@
 #include <string>
 
 #include "Nucleus.h"
+#include "RxnType.h"
+#include "MaskFile.h"
 
 #include <THashTable.h>
 #include <TH1.h>
@@ -29,6 +31,8 @@ public:
 		return table; 
 	};
 	void FillData(const Mask::Nucleus& nuc, double detKE = 0.0, const std::string& modifier = "");
+	void FillCorrelations(const Mask::MaskFileData& data, Mask::RxnType type);
+	void FillCorrelationsDetected(const Mask::MaskFileData& data, Mask::RxnType type);
 
 private:
 	THashTable* table;
diff --git a/include/SabreEfficiency.h b/include/SabreEfficiency.h
index a5346a6..ea866c6 100644
--- a/include/SabreEfficiency.h
+++ b/include/SabreEfficiency.h
@@ -1,6 +1,7 @@
 #ifndef SABREEFFICIENCY_H
 #define SABREEFFICIENCY_H
 
+#include "MaskFile.h"
 #include "DetectorEfficiency.h"
 #include "SabreDetector.h"
 #include "Target.h"
@@ -18,6 +19,7 @@ public:
 
 private:
 	std::pair<bool,double> IsSabre(Mask::Nucleus& nucleus);
+    void CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, Mask::RxnType rxn_type);
 
 	std::vector<SabreDetector> detectors;
     
diff --git a/src/Detectors/DetectorEfficiency.cpp b/src/Detectors/DetectorEfficiency.cpp
index 8671b37..5250ad7 100644
--- a/src/Detectors/DetectorEfficiency.cpp
+++ b/src/Detectors/DetectorEfficiency.cpp
@@ -15,15 +15,16 @@ int main(int argc, char** argv) {
 	std::string outputname = argv[2];
 	std::string statsname = argv[3];
 
-	/*
+	
 	SabreEfficiency sabre;
 	std::string mapfile = "./etc/SabreDeadChannels.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");
-	*/
+	
 
+	/*
 	try {
  		AnasenEfficiency anasen;
  		std::string mapfile = "./etc/AnasenDeadChannels.txt";
@@ -36,5 +37,6 @@ int main(int argc, char** argv) {
 		std::cerr<<"Terminating."<<std::endl;
 		return 1;
 	}
+	*/
 	return 0;
-}
\ No newline at end of file
+}
diff --git a/src/Detectors/SabreEfficiency.cpp b/src/Detectors/SabreEfficiency.cpp
index 854b255..06f91b8 100644
--- a/src/Detectors/SabreEfficiency.cpp
+++ b/src/Detectors/SabreEfficiency.cpp
@@ -39,7 +39,7 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
 
 
 	Mask::MaskFile input(inputname, Mask::MaskFile::FileType::read);
-	Mask::MaskFile output(outputname, Mask::MaskFile::FileType::read);
+	Mask::MaskFile output(outputname, Mask::MaskFile::FileType::write);
 	std::ofstream stats(statsname);
 	stats<<std::setprecision(5);
 
@@ -52,18 +52,23 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
 	Mask::Nucleus nucleus;
 
 	std::vector<int> counts;
+	std::vector<int> coinc_counts;
 	switch(header.rxn_type) {
 		case Mask::RxnType::PureDecay:
 			counts.resize(3, 0);
+			coinc_counts.resize(1, 0);
 			break;
 		case Mask::RxnType::OneStepRxn:
 			counts.resize(4, 0);
+			coinc_counts.resize(1, 0);
 			break;
 		case Mask::RxnType::TwoStepRxn:
 			counts.resize(6, 0);
+			coinc_counts.resize(4, 0);
 			break;
 		case Mask::RxnType::ThreeStepRxn:
 			counts.resize(8, 0);
+			coinc_counts.resize(11, 0);
 			break;
 		default:
 		{
@@ -96,12 +101,15 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
 			auto result = IsSabre(nucleus);
 			if(result.first) {
 				data.detect_flag[i] = true;
+				data.KE[i] = result.second;
 				counts[i]++;
 			} else if(data.detect_flag[i] == true) {
 				data.detect_flag[i] = false;
 			}
 		}
 
+		CountCoincidences(data, coinc_counts, header.rxn_type);
+
 		output.WriteData(data);
 	}
 
@@ -114,6 +122,54 @@ void SabreEfficiency::CalculateEfficiency(const std::string& inputname, const st
 		stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/header.nsamples)<<std::endl;
 		stats<<"---------------------"<<std::endl;
 	}
+	stats<<"Coincidence Efficiency"<<std::endl;
+	stats<<"---------------------"<<std::endl;
+	if(header.rxn_type == Mask::RxnType::PureDecay)
+	{
+		stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+	}
+	else if(header.rxn_type == Mask::RxnType::OneStepRxn)
+	{
+		stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+	}
+	else if(header.rxn_type == Mask::RxnType::TwoStepRxn)
+	{
+		stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+	}
+	else if(header.rxn_type == Mask::RxnType::ThreeStepRxn)
+	{
+		stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+		stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/header.nsamples)<<std::endl;
+		stats<<"---------------------"<<std::endl;
+	}
 	stats.close();
 
 	std::cout<<std::endl;
@@ -206,3 +262,80 @@ std::pair<bool,double> SabreEfficiency::IsSabre(Mask::Nucleus& nucleus) {
 
 	return std::make_pair(false,0.0);
 }
+
+void SabreEfficiency::CountCoincidences(const Mask::MaskFileData& data, std::vector<int>& counts, Mask::RxnType rxn_type) {
+	if (rxn_type == Mask::RxnType::PureDecay && data.detect_flag[1] && data.detect_flag[2])
+	{
+		counts[0]++;
+	}
+	else if (rxn_type == Mask::RxnType::OneStepRxn && data.detect_flag[2] && data.detect_flag[3])
+	{
+		counts[0]++;
+	}
+	else if(rxn_type == Mask::RxnType::TwoStepRxn)
+	{
+		if(data.detect_flag[2] && data.detect_flag[4]) 
+		{
+			counts[0]++;
+		}
+		if(data.detect_flag[2] && data.detect_flag[5])
+		{
+			counts[1]++;
+		}
+		if(data.detect_flag[4] && data.detect_flag[5])
+		{
+			counts[2]++;
+		}
+		if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[5])
+		{
+			counts[3]++;
+		}
+	}
+	else if(rxn_type == Mask::RxnType::ThreeStepRxn)
+	{
+		if(data.detect_flag[2] && data.detect_flag[4]) 
+		{
+			counts[0]++;
+		}
+		if(data.detect_flag[2] && data.detect_flag[6])
+		{
+			counts[1]++;
+		}
+		if(data.detect_flag[2] && data.detect_flag[7])
+		{
+			counts[2]++;
+		}
+		if(data.detect_flag[4] && data.detect_flag[6])
+		{
+			counts[3]++;
+		}
+		if(data.detect_flag[4] && data.detect_flag[7])
+		{
+			counts[4]++;
+		}
+		if(data.detect_flag[6] && data.detect_flag[7])
+		{
+			counts[5]++;
+		}
+		if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6])
+		{
+			counts[6]++;
+		}
+		if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[7])
+		{
+			counts[7]++;
+		}
+		if(data.detect_flag[2] && data.detect_flag[6] && data.detect_flag[7])
+		{
+			counts[8]++;
+		}
+		if(data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
+		{
+			counts[9]++;
+		}
+		if(data.detect_flag[2] && data.detect_flag[4] && data.detect_flag[6] && data.detect_flag[7])
+		{
+			counts[10]++;
+		}
+	}
+}
\ No newline at end of file
diff --git a/src/MaskApp.cpp b/src/MaskApp.cpp
index 52e593d..a4365e8 100644
--- a/src/MaskApp.cpp
+++ b/src/MaskApp.cpp
@@ -76,7 +76,7 @@ namespace Mask {
 			case RxnType::ThreeStepRxn:
 			{
 				m_sys = new ThreeStepSystem();
-				m_rxn_type = RxnType::TwoStepRxn;
+				m_rxn_type = RxnType::ThreeStepRxn;
 				for(int i=0; i<5; i++) {
 					input>>z>>a;
 					avec.push_back(a);
diff --git a/src/Plotters/ROOT/RootPlotter.cpp b/src/Plotters/ROOT/RootPlotter.cpp
index a7259da..9df6df3 100644
--- a/src/Plotters/ROOT/RootPlotter.cpp
+++ b/src/Plotters/ROOT/RootPlotter.cpp
@@ -1,5 +1,4 @@
 #include "RootPlotter.h"
-#include "MaskFile.h"
 #include <TFile.h>
 
 #include <iostream>
@@ -39,6 +38,60 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc, double detKE, const std::st
 	
 }
 
+void RootPlotter::FillCorrelations(const Mask::MaskFileData& data, Mask::RxnType type)
+{
+	std::string theta_eject_theta_resid_name = "theta_eject_theta_resid_cor";
+	std::string theta_eject_theta_resid_title = theta_eject_theta_resid_name + ";#theta_{lab} Ejectile (deg);#theta_{lab} Residual";
+	if(type == Mask::RxnType::PureDecay)
+	{
+		MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[1]*rad2deg, data.theta[2]*rad2deg, 4);
+	}
+	else
+	{
+		MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[2]*rad2deg, data.theta[3]*rad2deg, 4);
+	}
+
+	if(type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn)
+	{
+		std::string theta_break1_theta_break2_name = "theta_break1_theta_break2_cor";
+		std::string theta_break1_theta_break2_title = theta_break1_theta_break2_name + ";#theta_{lab} Breakup1 (deg);#theta_{lab} Breakup2 (deg)";
+		MyFill(theta_break1_theta_break2_name, theta_break1_theta_break2_title, data.theta[4]*rad2deg, data.theta[5]*rad2deg, 4);
+	}
+	if(type == Mask::RxnType::ThreeStepRxn)
+	{
+		std::string theta_break3_theta_break4_name = "theta_break3_theta_break4_cor";
+		std::string theta_break3_theta_break4_title = theta_break3_theta_break4_name + ";#theta_{lab} Breakup3 (deg);#theta_{lab} Breakup4 (deg)";
+		MyFill(theta_break3_theta_break4_name, theta_break3_theta_break4_title, data.theta[6]*rad2deg, data.theta[7]*rad2deg, 4);
+	}
+}
+
+void RootPlotter::FillCorrelationsDetected(const Mask::MaskFileData& data, Mask::RxnType type)
+{
+	std::string theta_eject_theta_resid_name = "theta_eject_theta_resid_cor_detected";
+	std::string theta_eject_theta_resid_title = theta_eject_theta_resid_name + ";#theta_{lab} Ejectile (deg);#theta_{lab} Residual";
+	if(type == Mask::RxnType::PureDecay && data.detect_flag[1] && data.detect_flag[2])
+	{
+		MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[1]*rad2deg, data.theta[2]*rad2deg, 4);
+	}
+	else if(data.detect_flag[2] && data.detect_flag[3])
+	{
+		MyFill(theta_eject_theta_resid_name, theta_eject_theta_resid_title, data.theta[2]*rad2deg, data.theta[3]*rad2deg, 4);
+	}
+
+	if((type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) && data.detect_flag[4] && data.detect_flag[5])
+	{
+		std::string theta_break1_theta_break2_name = "theta_break1_theta_break2_cor_detected";
+		std::string theta_break1_theta_break2_title = theta_break1_theta_break2_name + ";#theta_{lab} Breakup1 (deg);#theta_{lab} Breakup2 (deg)";
+		MyFill(theta_break1_theta_break2_name, theta_break1_theta_break2_title, data.theta[4]*rad2deg, data.theta[5]*rad2deg, 4);
+	}
+	if(type == Mask::RxnType::ThreeStepRxn && data.detect_flag[6] && data.detect_flag[7])
+	{
+		std::string theta_break3_theta_break4_name = "theta_break3_theta_break4_cor_detected";
+		std::string theta_break3_theta_break4_title = theta_break3_theta_break4_name + ";#theta_{lab} Breakup3 (deg);#theta_{lab} Breakup4 (deg)";
+		MyFill(theta_break3_theta_break4_name, theta_break3_theta_break4_title, data.theta[6]*rad2deg, data.theta[7]*rad2deg, 4);
+	}
+}
+
 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) {
@@ -224,6 +277,8 @@ int main(int argc, char** argv) {
 				plotter.FillData(nucleus, data.KE[i], "detected");
 			}
 		}
+		plotter.FillCorrelations(data, header.rxn_type);
+		plotter.FillCorrelationsDetected(data, header.rxn_type);
 		count++;
 	}
 	std::cout<<std::endl;