Merge pull request #8 from gwm17/devel

Devel

.gitignore

@@ -18,6 +18,8 @@ Makefile
 build/
 lib/
 bin/
+scripts/
+*.yaml
 
 
 ###Keep this file###

.gitmodules

@@ -1,3 +1,6 @@
 [submodule "src/vendor/catima"]
 	path = src/vendor/catima
 	url = https://github.com/gwm17/catima.git
+[submodule "src/vendor/yaml-cpp"]
+	path = src/vendor/yaml-cpp
+	url = https://github.com/jbeder/yaml-cpp.git

CMakeLists.txt

@@ -19,6 +19,8 @@ set(CMAKE_CXX_STANDARD 17)
 find_package(ROOT REQUIRED COMPONENTS GenVector)
 
 add_subdirectory(src/vendor/catima)
+set(YAML_CPP_BUILD_TOOLS Off CACHE BOOL "Enable parse tools")
+add_subdirectory(src/vendor/yaml-cpp)
 add_subdirectory(src/Mask)
 add_subdirectory(src/Kinematics)
 add_subdirectory(src/Detectors)

README.md

@@ -1,9 +1,12 @@
 # Mask: Monte cArlo Simulation of Kinematics
+
 Mask is a Monte Carlo simulation of reaction kinematics for use detector systems at Florida State University.
 Mask is capable of simulating multi-step kinematic reaction-decay sequences, storing data in ROOT Trees, after which the kinematic data can be fed to a detector geometry for efficiency testing. Currently geometries for ANASEN and SABRE are included in the code.
 
 ## Building Mask
-Dowload the repository from github. CMake is use to build the project; in most environments you can build Mask using the following methods:
+
+Dowload the repository from GitHub using `git clone --recursive https://github.com/gwm17/Mask.git`. CMake is use to build the project; in most environments you can build Mask using the following methods:
+
 - `mkdir build`
 - `cd build`
 - `cmake ..`
@@ -14,6 +17,7 @@ Executables will be installed to the repostiory's `bin` directory. Libraries wil
 By default Mask builds for release. To build for debug replace `cmake ..` with `cmake -DCMAKE_BUILD_TYPE=Debug ..`. Mask uses CMake to find the installed ROOT libraries and headers.
 
 ## Using the kinematics simulation
+
 By default Mask is capable of simulating reactions of up to three steps. Here is a brief outline of each type:
 
 0. A pure decay involves a "target" decaying into an ejectile and residual. It is assumed isotropic by default; any other case will require the modification of the code.
@@ -25,11 +29,12 @@ For decays, a specific angular distribution can be given as input as a text file
 
 To run Mask simply do the following from the Mask repository:
 
-`./bin/Kinematics input.txt`
+`./bin/Kinematics <your_config>.yaml`
 
-Input.txt can be replaced by any text file with the correct format.
+`<your_config.yaml>` is a YAML configuration file. An example is given in the repository named `kinematics.yaml` and can be replaced by any yaml file with the correct format.
 
 ## Using the detector geometry simulation
+
 Detector geometry is encoded using ROOT math libraries in the `src/Detectors` folder. Two different detector geometries are already present: SPS-SABRE and ANASEN. To add a new geometry, follow the guidelines outlined by each of these cases.
 
 To choose which detector scheme is run, modify the main function in `src/Detectors/main.cpp`. The included geometries also have options to do an internal geometry consistency check and print out coordinates for drawing the detector arrays, which can be useful for testing.
@@ -38,11 +43,12 @@ To run the geometry code, one needs to provide an input file containing the foll
 
 To run Detectors use the format
 
-`./bin/Detectors <input_file>`
+`./bin/Detectors <your_config>.yaml`
 
-An example input file is provided with the repository.
+`<your_config>.yaml` is a YAML configuration file. An example, `detector.yaml` is included in the repository.
 
 ## Data visualization
+
 All data is saved as ROOT trees of std::vectors of Mask::Nucleus classes. To enable this, a ROOT dictionary is generated and linked into a shared library found in the `lib` directory of the repository. This allows the user to link to the shared library for accessing and analyzing the data generated by Mask.
 
 Mask also provides a default visualization tool called RootPlot. RootPlot is run as
@@ -52,6 +58,6 @@ Mask also provides a default visualization tool called RootPlot. RootPlot is run
 where the datafile can be either the datafile from Mask or the datafile from DetEff. The outputfile is saved in the ROOT file format.
 
 ## Requirements
+
 ROOT version 6.22 or greater is required
 CMake version 3.0 or greater is required
-

detector.yaml

@@ -0,0 +1,5 @@
+InputDataFile: /media/data/gwm17/mask_tests/temp.root
+OutputDataFile: /media/data/gwm17/mask_tests/temp_det.root
+DeadChannelFile: etc/sabreDeadChannels_May2022.txt
+NumberOfThreads: 5
+ArrayType: Sabre

detector_input.txt

@@ -1,5 +0,0 @@
-InputDataFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
-OutputDataFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B_Sabre.root
-DeadChannelFile: etc/sabreDeadChannels_May2022.txt
-NumberOfThreads: 5
-ArrayType: Sabre

kinematics.yaml

@@ -0,0 +1,49 @@
+OutputFile: /media/data/gwm17/mask_tests/temp.root
+Threads: 5
+ReactionSamples: 1000000
+ReactionChain:
+  - Type: Reaction
+    Reactants:
+      - Z: 5
+        A: 10
+      - Z: 2
+        A: 3
+      - Z: 2
+        A: 4
+    BeamEnergyMean(MeV): 24.0
+    BeamEnergySigma(MeV): 0.0
+    ThetaType: Lab
+    ThetaMin(deg): 15.0
+    ThetaMax(deg): 15.0
+    PhiMin(deg): 0.0
+    PhiMax(deg): 0.0
+    ResidualExcitationMean(MeV): 16.798
+    ResidualExcitationSigma(MeV): 0.033
+  - Type: Decay
+    Reactants:
+      - Z: 5
+        A: 9
+      - Z: 2
+        A: 4
+    PhiMin(deg): 0.0
+    PhiMax(deg): 360.0
+    ResidualExcitationMean(MeV): 0.0
+    ResidualExcitationSigma(MeV): 0.522
+    AngularDistributionFile: ./etc/isotropic_dist.txt
+  - Type: Decay
+    Reactants:
+      - Z: 3
+        A: 5
+      - Z: 2
+        A: 4
+    PhiMin(deg): 0.0
+    PhiMax(deg): 360.0
+    ResidualExcitationMean(MeV): 0.0
+    ResidualExcitationSigma(MeV): 0.0
+    AngularDistributionFile: ./etc/isotropic_dist.txt
+TargetLayers:
+  - Thickness(ug/cm^2): 74
+    Elements:
+      - Z: 5
+        A: 10
+        S: 1

kinematics_input.txt

@@ -1,45 +0,0 @@
-----------Data Information----------
-OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
-NumberOfThreads: 6
-----------Reaction Information----------
-NumberOfSamples: 1000000
-begin_chain
-	begin_step
-		Type: Reaction
-		begin_nuclei
-			5 10
-			2 3
-			2 4
-		end_nuclei
-		BeamEnergyMean(MeV): 24.0
-		BeamEnergySigma(MeV): 0.0
-		ThetaType: Lab
-		ThetaMin(deg): 15.0
-		ThetaMax(deg): 15.0
-		PhiMin(deg): 0.0
-		PhMax(deg): 0.0
-		ResidualExcitationMean(MeV): 16.8
-		ResidualExcitationSigma(MeV): 0.023
-	end_step
-	begin_step
-		Type: Decay
-		begin_nuclei
-			5 9
-			2 4
-		end_nuclei
-		PhiMin(deg): 0.0
-		PhMax(deg): 360.0
-		ResidualExcitationMean(MeV): 0.0
-		ResidualExcitationSigma(MeV): 0.0
-		AngularDistributionFile: ./etc/isotropic_dist.txt
-	end_step
-end_chain
-----------Target Information----------
-NumberOfLayers: 1
-begin_layer
-	Thickness(ug/cm^2): 74
-	begin_elements (Z, A, Stoich.)
-		element	5 10 1
-	end_elements
-end_layer
---------------------------------------

src/Detectors/CMakeLists.txt

@@ -1,5 +1,9 @@
 add_executable(Detectors)
-target_include_directories(Detectors PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..)
+target_include_directories(Detectors PUBLIC 
+    ${CMAKE_CURRENT_SOURCE_DIR} 
+    ${CMAKE_CURRENT_SOURCE_DIR}/..
+    ${CMAKE_CURRENT_SOURCE_DIR}/../vendor/yaml-cpp/include/
+)
 
 target_sources(Detectors PUBLIC
     AnasenDeadChannelMap.cpp

src/Detectors/DetectorApp.cpp

@@ -2,6 +2,8 @@
 #include <fstream>
 #include <iostream>
 
+#include "yaml-cpp/yaml.h"
+
 DetectorApp::DetectorApp() :
     m_resources(nullptr)
 {
@@ -16,38 +18,31 @@ DetectorApp::~DetectorApp()
 bool DetectorApp::LoadConfig(const std::string& filename)
 {
     std::cout<<"----------Detector Efficiency Calculation----------"<<std::endl;
-    std::ifstream input(filename);
-    if(!input.is_open())
+    YAML::Node data;
+    try
     {
-        std::cerr<<"Unable to open input config file "<<filename<<std::endl;
+        data = YAML::LoadFile(filename);
+    }
+    catch (YAML::ParserException& e)
+    {
+        std::cerr << "Could not load config file " << filename  << " with error: " << e.what() << std::endl;
         return false;
     }
 
-    std::string junk;
-    std::string typeName;
+    m_inputFileName = data["InputDataFile"].as<std::string>();
+    m_outputFileName = data["OutputDataFile"].as<std::string>();
+    m_deadChannelFileName = data["DeadChannelFile"].as<std::string>();
+    m_nthreads = data["NumberOfThreads"].as<uint64_t>();
+    ArrayType type = StringToArrayType(data["ArrayType"].as<std::string>());
 
-    input >> junk >> m_inputFileName;
-    input >> junk >> m_outputFileName;
-    input >> junk >> m_deadChannelFileName;
-
-    input >> junk >> m_nthreads;
-    input >> junk >> typeName;
-
-    std::cout << "Creating " << m_nthreads << " detector arrays..." << std::endl;
-    ArrayType type = StringToArrayType(typeName);
     for(uint64_t i=0; i<m_nthreads; i++)
     {
         m_detectorList.push_back(CreateDetectorArray(type));
         if(m_deadChannelFileName != "None")
             m_detectorList.back()->SetDeadChannelMap(m_deadChannelFileName);
     }
-    std::cout << "Done" << std::endl;
-
-    std::cout << "Allocating " << m_nthreads << " threads..." << std::endl;
     m_resources = std::make_unique<Mask::ThreadPool<DetectorArray*, uint64_t>>(m_nthreads);
-    std::cout << "Done" << std::endl;
 
-    std::cout << "Opening input data file " << m_inputFileName << "..." << std::endl;
     m_fileReader.Open(m_inputFileName, "SimTree");
     if(!m_fileReader.IsOpen() || !m_fileReader.IsTree())
     {
@@ -55,7 +50,6 @@ bool DetectorApp::LoadConfig(const std::string& filename)
         return false;
     }
     m_nentries = m_fileReader.GetSize();
-    std::cout << "Done. Detected " << m_nentries << " events in the file." << std::endl;
 
     //Little bit of integer division mangling to make sure we read every event in file
     uint64_t quotient = m_nentries / m_nthreads;
@@ -64,17 +58,17 @@ bool DetectorApp::LoadConfig(const std::string& filename)
     for(uint64_t i=1; i<m_nthreads; i++)
         m_chunkSamples.push_back(quotient);
 
-    std::cout << "Opening output data file " << m_outputFileName << std::endl;
     m_fileWriter.Open(m_outputFileName, "SimTree");
     if(!m_fileWriter.IsOpen() || !m_fileWriter.IsTree())
     {
         std::cerr << "Unable to open output data file " << m_outputFileName << std::endl;
         return false;
     }
-    std::cout << "Done" << std::endl;
-
-    std::cout << "Ready to launch!" << std::endl;
 
+    std::cout << "Allocating " << m_nthreads << " threads..." << std::endl;
+    std::cout << "Input data file " << m_inputFileName << "..." << std::endl;
+    std::cout << "With " << m_nentries << " events in the file..." << std::endl;
+    std::cout << "Output data file " << m_outputFileName << "..." << std::endl;
     return true;
 }
 
@@ -110,6 +104,7 @@ void DetectorApp::Run()
                             nucleus.detectedKE = result.energy_deposited;
                             nucleus.detectedTheta = result.direction.Theta();
                             nucleus.detectedPhi = result.direction.Phi();
+                            nucleus.detectedPos = result.direction;
                         }
                     }
 		            m_fileWriter.PushData(data);
@@ -141,4 +136,6 @@ void DetectorApp::Run()
 				++count;
 	}
 
+    std::cout << std::endl;
+	
 }

src/Detectors/SabreArray.cpp

@@ -20,6 +20,12 @@ SabreArray::SabreArray() :
 	m_degradedDetectors[2] = false;
 	m_degradedDetectors[3] = false;
 	m_degradedDetectors[4] = true;
+    //No degraded detectors
+	// m_degradedDetectors[0] = false;
+	// m_degradedDetectors[1] = false;
+	// m_degradedDetectors[2] = false;
+	// m_degradedDetectors[3] = false;
+	// m_degradedDetectors[4] = false;
 
 	//Choose who to look at right now. Usually switch on or off degraded/non-degraded.
 	m_activeDetectors[0] = false;

src/Detectors/SabreArray.h

@@ -39,7 +39,7 @@ private:
     static constexpr double s_detectorThickness = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density)
     static constexpr double s_degraderThickness = 70.0 * 1.0e-4 * 16.69 * 1e6; //tantalum degrader (70 um thick) 
 
-    static constexpr double s_energyThreshold = 0.2; //in MeV
+    static constexpr double s_energyThreshold = 0.25; //in MeV
 
 };
 

src/Kinematics/main.cpp

@@ -24,7 +24,6 @@ int main(int argc, char** argv)
 			return 1;
 		}
 		calculator.Run();
-		//calculator.RunSingleThread();
 	}
 	catch(const std::exception& e)
 	{

src/Mask/AngularDistribution.cpp

@@ -73,13 +73,10 @@ namespace Mask {
 		m_branchingRatio = m_constants[0]*2.0;
 		m_L = l;
 
-		std::cout<<"Angular distribution from "<<file<<" branching ratio: "<<m_branchingRatio<<std::endl;
-		std::cout<<"Angular distribution from "<<file<<" L: "<<m_L<<std::endl;
-	
 		//Renormalize distribution such that total prob is 1.0.
 		//Test branching ratio to see if we "make" a decay particle,
 		//then use re-normalized distribution to pick an angle. 
-		if(m_constants[0] < 0.5)
+		if(m_constants[0] != 0.5)
 		{
 			double norm = 0.5/m_constants[0];
 			for(auto& value : m_constants)
@@ -115,4 +112,14 @@ namespace Mask {
 		return costheta;
 	}
 
+	double AngularDistribution::GetProbability(double cosTheta)
+	{
+		double prob = 0.0;
+		for (std::size_t i=0; i<m_constants.size(); i++)
+		{
+			prob += m_constants[i] * P_l(i*2, cosTheta);
+		}
+		return prob;
+	}
+
 }

src/Mask/AngularDistribution.h

@@ -17,6 +17,7 @@ namespace Mask {
 		double GetRandomCosTheta();
 		int GetL() { return m_L; }
 		double GetBranchingRatio() { return m_branchingRatio; }
+		double GetProbability(double cosTheta);
 	
 	private:
 		bool IsIsotropic() { return m_isIsotropic; }

src/Mask/CMakeLists.txt

@@ -1,6 +1,7 @@
 add_library(MaskDict SHARED)
 target_include_directories(MaskDict
     PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+    PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../vendor/yaml-cpp/include/
     SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}
 )
 
@@ -55,10 +56,14 @@ target_sources(Mask PRIVATE
     FileWriter.cpp
     FileReader.h
     FileReader.cpp
+    CoupledThreeStepSystem.h
+    CoupledThreeStepSystem.cpp
+    ConfigSerializer.h
+    ConfigSerializer.cpp
 )
 
 set(THREADS_PREFER_PTHREAD_FLAG On)
 find_package(Threads REQUIRED)
-
-target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES} Threads::Threads)
+target_compile_definitions(Mask PRIVATE YAML_CPP_STATIC_DEFINE)
+target_link_libraries(Mask catima yaml-cpp MaskDict ${ROOT_LIBRARIES} Threads::Threads)
 set_target_properties(Mask PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})

src/Mask/ConfigSerializer.cpp

@@ -0,0 +1,223 @@
+#include "ConfigSerializer.h"
+#include "yaml-cpp/yaml.h"
+#include "RxnType.h"
+
+namespace Mask {
+
+    static void SerializeDecay(YAML::Emitter& yamlStream, const StepParameters& params)
+    {
+        yamlStream << YAML::BeginMap;
+		yamlStream << YAML::Key << "Type" << RxnTypeToString(params.rxnType); 
+		yamlStream << YAML::Key << "Reactants" << YAML::Value << YAML::BeginSeq;
+		for (std::size_t i=0; i<params.Z.size(); i++)
+		{
+			yamlStream << YAML::BeginMap;
+			yamlStream << YAML::Key << "Z" << YAML::Value << params.Z[i];
+			yamlStream << YAML::Key << "A" << YAML::Value << params.A[i];
+			yamlStream << YAML::EndMap;
+		}
+		yamlStream << YAML::EndSeq;
+        yamlStream << YAML::Key << "PhiMin(deg)" << YAML::Value << params.phiMin;
+        yamlStream << YAML::Key << "PhiMax(deg)" << YAML::Value << params.phiMax;
+        yamlStream << YAML::Key << "ResidualExcitationMean(MeV)" << YAML::Value << params.meanResidualEx;
+        yamlStream << YAML::Key << "ResidualExcitationSigma(MeV)" << YAML::Value << params.sigmaResidualEx;
+		yamlStream << YAML::Key << "AngularDistributionFile" << YAML::Value << params.angularDistFile;
+		yamlStream << YAML::EndMap;
+    }
+
+    static StepParameters DeserializeDecay(YAML::Node& yamlStream)
+    {
+        StepParameters params;
+        params.rxnType = StringToRxnType(yamlStream["Type"].as<std::string>());
+        YAML::Node nuclei = yamlStream["Reactants"];
+        for(auto nucleus : nuclei)
+        {
+            params.Z.push_back(nucleus["Z"].as<int>());
+            params.A.push_back(nucleus["A"].as<int>());
+        }
+        params.phiMin = yamlStream["PhiMin(deg)"].as<double>();
+        params.phiMax = yamlStream["PhiMax(deg)"].as<double>();
+        params.meanResidualEx = yamlStream["ResidualExcitationMean(MeV)"].as<double>();
+        params.sigmaResidualEx = yamlStream["ResidualExcitationSigma(MeV)"].as<double>();
+        params.angularDistFile = yamlStream["AngularDistributionFile"].as<std::string>();
+        return params;
+    }
+
+    static void SerializeReaction(YAML::Emitter& yamlStream, const StepParameters& params)
+    {
+        yamlStream << YAML::BeginMap;
+		yamlStream << YAML::Key << "Type" << RxnTypeToString(params.rxnType); 
+		yamlStream << YAML::Key << "Reactants" << YAML::Value << YAML::BeginSeq;
+		for (std::size_t i=0; i<params.Z.size(); i++)
+		{
+			yamlStream << YAML::BeginMap;
+			yamlStream << YAML::Key << "Z" << YAML::Value << params.Z[i];
+			yamlStream << YAML::Key << "A" << YAML::Value << params.A[i];
+			yamlStream << YAML::EndMap;
+		}
+		yamlStream << YAML::EndSeq;
+        yamlStream << YAML::Key << "BeamEnergyMean(MeV)" << YAML::Value << params.meanBeamEnergy;
+        yamlStream << YAML::Key << "BeamEnergySigma(MeV)" << YAML::Value << params.sigmaBeamEnergy;
+        yamlStream << YAML::Key << "ThetaType" << YAML::Value << RxnThetaTypeToString(params.thetaType);
+        yamlStream << YAML::Key << "ThetaMin(deg)" << YAML::Value << params.thetaMin;
+        yamlStream << YAML::Key << "ThetaMax(deg)" << YAML::Value << params.thetaMax;
+        yamlStream << YAML::Key << "PhiMin(deg)" << YAML::Value << params.phiMin;
+        yamlStream << YAML::Key << "PhiMax(deg)" << YAML::Value << params.phiMax;
+        yamlStream << YAML::Key << "ResidualExcitationMean(MeV)" << YAML::Value << params.meanResidualEx;
+        yamlStream << YAML::Key << "ResidualExcitationSigma(MeV)" << YAML::Value << params.sigmaResidualEx;
+		yamlStream << YAML::EndMap;
+    }
+
+    static StepParameters DeserializeReaction(YAML::Node& yamlStream)
+    {
+        StepParameters params;
+        params.rxnType = StringToRxnType(yamlStream["Type"].as<std::string>());
+        YAML::Node nuclei = yamlStream["Reactants"];
+        for(auto nucleus : nuclei)
+        {
+            params.Z.push_back(nucleus["Z"].as<int>());
+            params.A.push_back(nucleus["A"].as<int>());
+        }
+        params.meanBeamEnergy = yamlStream["BeamEnergyMean(MeV)"].as<double>();
+        params.sigmaBeamEnergy = yamlStream["BeamEnergySigma(MeV)"].as<double>();
+        params.thetaType = StringToRxnThetaType(yamlStream["ThetaType"].as<std::string>());
+        params.thetaMin = yamlStream["ThetaMin(deg)"].as<double>();
+        params.thetaMax = yamlStream["ThetaMax(deg)"].as<double>();
+        params.phiMin = yamlStream["PhiMin(deg)"].as<double>();
+        params.phiMax = yamlStream["PhiMax(deg)"].as<double>();
+        params.meanResidualEx = yamlStream["ResidualExcitationMean(MeV)"].as<double>();
+        params.sigmaResidualEx = yamlStream["ResidualExcitationSigma(MeV)"].as<double>();
+        return params;
+    }
+
+    static void SerializeTarget(YAML::Emitter& yamlStream, const LayeredTarget& target)
+    {
+
+        yamlStream << YAML::BeginSeq;
+        for (int i=0; i<target.GetNumberOfLayers(); i++)
+        {
+            const Target& layer = target.GetLayerInfo(i);
+            yamlStream << YAML::BeginMap;
+            yamlStream << YAML::Key << "Thickness(ug/cm^2)" << YAML::Value << layer.GetThickness();
+            yamlStream << YAML::Key << "Elements" << YAML::Value << YAML::BeginSeq;
+            for (int j=0; j<layer.GetNumberOfElements(); j++)
+            {
+                yamlStream << YAML::BeginMap;
+                yamlStream << YAML::Key << "Z" << YAML::Value << layer.GetElementZ(j);
+                yamlStream << YAML::Key << "A" << YAML::Value << layer.GetElementA(j);
+                yamlStream << YAML::Key << "S" << YAML::Value << layer.GetElementStoich(j);
+                yamlStream << YAML::EndMap;
+            }
+            yamlStream << YAML::EndSeq;
+            yamlStream << YAML::EndMap;
+        }
+        yamlStream << YAML::EndSeq;
+    }
+
+    static LayeredTarget DeserializeTarget(YAML::Node& yamlStream)
+    {
+        LayeredTarget target;
+        std::vector<int> z, a, s;
+        double thickness;
+
+        for (auto layer : yamlStream)
+        {
+            thickness = layer["Thickness(ug/cm^2)"].as<double>();
+            YAML::Node elements = layer["Elements"];
+            z.clear();
+            a.clear();
+            s.clear();
+            for(auto element : elements)
+            {
+                z.push_back(element["Z"].as<int>());
+                a.push_back(element["A"].as<int>());
+                s.push_back(element["S"].as<int>());
+            }
+            target.AddLayer(z, a, s, thickness);
+        }
+
+        return target;
+    }
+
+    bool ConfigSerializer::SerializeConfig(const std::string& configfile, const AppParameters& params)
+    {
+		std::ofstream output(configfile);
+		if (!output.is_open())
+		{
+			std::cerr << "Could not open output file " << configfile << std::endl;
+			return false;
+		}
+
+		YAML::Emitter yamlStream;
+		yamlStream << YAML::BeginMap;
+		yamlStream << YAML::Key << "OutputFile" << YAML::Value << params.outputFileName;
+		yamlStream << YAML::Key << "Threads" << YAML::Value << params.nThreads;
+		yamlStream << YAML::Key << "ReactionSamples" << YAML::Value << params.nSamples;
+		yamlStream << YAML::Key << "ReactionChain" << YAML::Value << YAML::BeginSeq;
+		for (auto& step : params.chainParams)
+		{
+			switch(step.rxnType)
+            {
+                case RxnType::Reaction: SerializeReaction(yamlStream, step); break;
+                case RxnType::Decay: SerializeDecay(yamlStream, step); break;
+                case RxnType::None:
+                {
+                    std::cerr << "Error serializing config: none type reaction found!" << std::endl;
+                    return false;
+                }
+            }
+		}
+        yamlStream << YAML::EndSeq;
+        yamlStream << YAML::Key << "TargetLayers" << YAML::Value;
+        SerializeTarget(yamlStream, params.target);
+
+        output << yamlStream.c_str();
+        output.close();
+        return true;
+    }
+
+    bool ConfigSerializer::DeserializeConfig(const std::string& configfile, AppParameters& params)
+    {
+        YAML::Node data;
+		try 
+		{
+			data = YAML::LoadFile(configfile);
+		}
+		catch(YAML::ParserException& e)
+		{
+			std::cerr << "Error reading config " << configfile << ": " << e.what() <<std::endl;
+			return false;
+		}
+
+        params.outputFileName = data["OutputFile"].as<std::string>();
+        params.nThreads = data["Threads"].as<uint32_t>();
+        params.nSamples = data["ReactionSamples"].as<uint64_t>();
+
+        auto steps = data["ReactionChain"];
+        for (auto step : steps)
+        {
+            RxnType type = StringToRxnType(step["Type"].as<std::string>());
+            switch(type)
+            {
+                case RxnType::Decay:
+                {
+                    params.chainParams.push_back(DeserializeDecay(step));
+                    break;
+                }
+                case RxnType::Reaction:
+                {
+                    params.chainParams.push_back(DeserializeReaction(step));
+                    break;
+                }
+                case RxnType::None:
+                {
+                    std::cerr << "Error deserializing config: None-type reaction found" << std::endl;
+                    return false;
+                }
+            }
+        }
+        auto layers = data["TargetLayers"];
+        params.target = DeserializeTarget(layers);
+        return true;
+    }
+}

src/Mask/ConfigSerializer.h

@@ -0,0 +1,16 @@
+#ifndef CONFIG_SERIALIZER_H
+#define CONFIG_SERIALIZER_H
+
+#include "MaskApp.h"
+
+namespace Mask {
+
+    class ConfigSerializer
+    {
+    public:
+        static bool SerializeConfig(const std::string& configfile, const AppParameters& params);
+        static bool DeserializeConfig(const std::string& configfile, AppParameters& params);
+    };
+}
+
+#endif

src/Mask/CoupledThreeStepSystem.cpp

@@ -0,0 +1,210 @@
+#include "CoupledThreeStepSystem.h"
+
+#include "Math/Boost.h"
+#include "Math/Vector3D.h"
+#include "Math/RotationY.h"
+#include "Math/RotationZ.h"
+
+namespace Mask {
+
+    CoupledThreeStepSystem::CoupledThreeStepSystem(const std::vector<StepParameters>& params) :
+		ReactionSystem()
+	{
+		m_nuclei.resize(8);
+		Init(params);
+	}
+	
+	CoupledThreeStepSystem::~CoupledThreeStepSystem() {}
+	
+	void CoupledThreeStepSystem::Init(const std::vector<StepParameters>& params)
+	{
+		if(params.size() != 3 || params[0].rxnType != RxnType::Reaction || params[1].rxnType != RxnType::Decay ||
+		  params[2].rxnType != RxnType::Decay || params[0].Z.size() != 3 || params[0].A.size() != 3 || params[1].Z.size() != 2 ||
+		  params[1].A.size() != 2 || params[2].Z.size() != 2 || params[2].A.size() != 2)
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ThreeStepSystem::Init(), does not match ThreeStep signature!" << std::endl;
+			return;
+		}
+
+		const StepParameters& step1Params = params[0];
+		const StepParameters& step2Params = params[1];
+		const StepParameters& step3Params = params[2];
+
+		//Setup nuclei
+		int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2];
+		int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2];
+		if(zr != step2Params.Z[0] || ar != step2Params.A[0])
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step one and step two are not sequential! Step one recoil (Z,A): ("
+					  << zr << "," << ar << ") Step two target (Z,A): (" << step2Params.Z[0] << "," << step2Params.A[0] << ")" <<std::endl;
+			return;  
+		}
+		int zb2 = step2Params.Z[0] - step2Params.Z[1];
+		int ab2 = step2Params.A[0] - step2Params.A[1];
+		if(zb2 != step3Params.Z[0] || ab2 != step3Params.A[0])
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step two and step three are not sequential! Step two heavy (Z,A): ("
+					  << zb2 << "," << ab2 << ") Step three target (Z,A): (" << step3Params.Z[0] << "," << step3Params.A[0] << ")" <<std::endl;
+			return;
+		}
+		int zb4 = step3Params.Z[0] - step3Params.Z[1];
+		int ab4 = step3Params.A[0] - step3Params.A[1];
+
+		m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
+		m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //projectile
+		m_nuclei[2] = CreateNucleus(step1Params.Z[2], step1Params.A[2]); //ejectile
+		m_nuclei[3] = CreateNucleus(zr, ar); //residual
+		m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
+		m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
+		m_nuclei[6] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
+		m_nuclei[7] = CreateNucleus(zb4, ab4); //breakup4
+
+		m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
+		m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
+		m_step3.BindNuclei(&(m_nuclei[5]), nullptr, &(m_nuclei[6]), &(m_nuclei[7]));
+	
+		SetSystemEquation();
+
+		//Step one sampling parameters
+		AddBeamDistribution(step1Params.meanBeamEnergy, step1Params.sigmaBeamEnergy);
+		m_step1.SetEjectileThetaType(step1Params.thetaType);
+		AddThetaRange(step1Params.thetaMin, step1Params.thetaMax);
+		AddPhiRange(step1Params.phiMin, step1Params.phiMax);
+		AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
+
+		//Step two sampling parameters
+		AddPhiRange(step2Params.phiMin, step2Params.phiMax);
+		AddDecayAngularDistribution(step2Params.angularDistFile);
+		AddExcitationDistribution(step2Params.meanResidualEx, step2Params.sigmaResidualEx);
+
+		//Step three sampling parameters
+		AddPhiRange(step3Params.phiMin, step3Params.phiMax);
+		AddDecayAngularDistribution(step3Params.angularDistFile);
+		AddExcitationDistribution(step3Params.meanResidualEx, step3Params.sigmaResidualEx);
+	}
+
+	void CoupledThreeStepSystem::SetLayeredTarget(const LayeredTarget& target)
+	{
+		m_target = target;
+		m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
+		if(m_rxnLayer != m_target.GetNumberOfLayers())
+		{
+			m_step1.SetLayeredTarget(&m_target);
+			m_step2.SetLayeredTarget(&m_target);
+			m_step3.SetLayeredTarget(&m_target);
+			m_step1.SetRxnLayer(m_rxnLayer);
+			m_step2.SetRxnLayer(m_rxnLayer);
+			m_step3.SetRxnLayer(m_rxnLayer);
+			m_isTargetSet = true;
+		}
+		else
+			throw ReactionLayerException();
+	}
+	
+	void CoupledThreeStepSystem::SetSystemEquation()
+	{
+		std::stringstream stream;
+		stream << m_nuclei[0].isotopicSymbol << "("
+			   << m_nuclei[1].isotopicSymbol << ", "
+			   << m_nuclei[2].isotopicSymbol << ")"
+			   << m_nuclei[3].isotopicSymbol << "->"
+			   << m_nuclei[4].isotopicSymbol << "+"
+			   << m_nuclei[5].isotopicSymbol << "->"
+			   << m_nuclei[6].isotopicSymbol << "+"
+			   << m_nuclei[7].isotopicSymbol;
+		m_sysEquation = stream.str();
+	}
+
+	CoupledThreeStepParameters CoupledThreeStepSystem::SampleParameters()
+	{
+		CoupledThreeStepParameters params;
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		params.beamEnergy = m_beamDistributions[0](gen);
+		params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
+		params.rxnPhi = m_phiRanges[0](gen);
+		params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
+        params.cosRelativeAngle = m_decayAngularDistributions[1].GetRandomCosTheta();
+		params.decay1Theta = std::acos(params.cosdecay1Theta);
+		params.decay1Phi = m_phiRanges[1](gen);
+		params.residEx = m_exDistributions[0](gen);
+		params.decay1Ex = m_exDistributions[1](gen);
+		params.decay2Ex = m_exDistributions[2](gen);
+		params.rxnDepth = (m_rxnDepthDist(gen));
+		return params;
+	}
+
+    //Called after running step2
+    void CoupledThreeStepSystem::SampleCoupling(CoupledThreeStepParameters& params)
+    {
+        ROOT::Math::Boost boost(m_nuclei[5].vec4.BoostToCM());
+		ROOT::Math::PxPyPzEVector a1Vec = (boost * m_nuclei[4].vec4);
+		ROOT::Math::XYZVector a1PVec(a1Vec.Px(), a1Vec.Py(), a1Vec.Pz());
+        ROOT::Math::PxPyPzEVector a2Vec;
+		ROOT::Math::XYZVector a2PVec;
+		ROOT::Math::RotationZ zRot;
+		ROOT::Math::RotationY yRot;
+		a1PVec *= 1.0/a1Vec.P();
+        zRot.SetAngle(a1PVec.Phi());
+		yRot.SetAngle(a1PVec.Theta());
+		double relAngle = std::acos(params.cosRelativeAngle);
+		ROOT::Math::XYZVector a1PVecAligned = yRot.Inverse() * (zRot.Inverse() * a1PVec);
+		a2PVec.SetXYZ(
+			std::sin(relAngle),
+			0.0,
+			std::cos(relAngle)
+		);
+		a2PVec = zRot * (yRot * a2PVec);
+		params.decay2Theta = a2PVec.Theta();
+		params.decay2Phi = a2PVec.Phi();
+    }
+	
+	void CoupledThreeStepSystem::RunSystem()
+	{
+		CoupledThreeStepParameters params;
+		do
+		{
+			params = SampleParameters();
+		}
+		while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) 
+				&& m_step2.CheckDecayThreshold(params.residEx, params.decay1Ex)
+				&& m_step3.CheckDecayThreshold(params.decay1Ex, params.decay2Ex)));
+		
+		m_step1.SetReactionDepth(params.rxnDepth);
+		m_step1.SetBeamKE(params.beamEnergy);
+		m_step1.SetPolarRxnAngle(params.rxnTheta);
+		m_step1.SetAzimRxnAngle(params.rxnPhi);
+		m_step1.SetExcitation(params.residEx);
+	
+		m_step2.SetReactionDepth(params.rxnDepth);
+		m_step2.SetPolarRxnAngle(params.decay1Theta);
+		m_step2.SetAzimRxnAngle(params.decay1Phi);
+		m_step2.SetExcitation(params.decay1Ex);
+	
+		m_step3.SetReactionDepth(params.rxnDepth);
+		
+		m_step3.SetExcitation(params.decay2Ex);
+	
+		m_step1.Calculate();
+	
+		if(params.cosdecay1Theta == -10)
+		{
+			m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
+			m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
+			m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
+			m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
+			return;
+		}
+		m_step2.Calculate();
+	
+		m_step3.SetResidualEnergyLoss(true);
+
+		SampleCoupling(params); //must be called here as frame needs to be set
+
+		m_step3.SetPolarRxnAngle(params.decay2Theta);
+		m_step3.SetAzimRxnAngle(params.decay2Phi);
+		m_step3.Calculate();
+	}
+}

src/Mask/CoupledThreeStepSystem.h

@@ -0,0 +1,46 @@
+#ifndef COUPLEDTHREESTEPSYSTEM_H
+#define COUPLEDTHREESTEPSYSTEM_H
+
+#include "ReactionSystem.h"
+#include "AngularDistribution.h"
+
+namespace Mask {
+
+	struct CoupledThreeStepParameters
+	{
+		double beamEnergy = 0.;
+		double rxnTheta = 0.;
+		double rxnPhi = 0.;
+		double cosdecay1Theta = 0.;
+		double cosRelativeAngle = 0.;
+		double decay1Theta = 0.;
+		double decay1Phi = 0.;
+		double decay2Theta = 0.;
+		double decay2Phi = 0.;
+		double residEx = 0.;
+		double decay1Ex = 0.;
+		double decay2Ex = 0.;
+		double rxnDepth = 0.;
+	};
+
+	class CoupledThreeStepSystem : public ReactionSystem
+	{
+	public:
+		CoupledThreeStepSystem(const std::vector<StepParameters>& params);
+		~CoupledThreeStepSystem();
+
+		virtual void SetLayeredTarget(const LayeredTarget& target) override;
+		virtual void RunSystem() override;
+		
+	protected:
+		void Init(const std::vector<StepParameters>& params);
+		void SetSystemEquation() override;
+		CoupledThreeStepParameters SampleParameters();
+        void SampleCoupling(CoupledThreeStepParameters& params);
+	
+		Reaction m_step1, m_step2, m_step3;
+	};
+
+}
+
+#endif

src/Mask/LayeredTarget.h

@@ -31,9 +31,9 @@ namespace Mask {
 		double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth);
 		double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth);
 		std::size_t FindLayerContaining(int Z, int A);
-		std::size_t GetNumberOfLayers() { return m_layers.size(); }
+		std::size_t GetNumberOfLayers() const { return m_layers.size(); }
 		void SetName(std::string& n) { m_name = n; }
-		const Target& GetLayerInfo(int index) { return m_layers[index]; }
+		const Target& GetLayerInfo(int index) const { return m_layers[index]; }
 		const std::string& GetName() { return m_name; }
 	
 	private:

src/Mask/MaskApp.cpp

@@ -1,185 +1,80 @@
 #include "MaskApp.h"
+#include "ConfigSerializer.h"
 #include <fstream>
 #include <iostream>
 
 #include "TFile.h"
 #include "TTree.h"
 
+#include "yaml-cpp/yaml.h"
+
 namespace Mask {
 
 	MaskApp::MaskApp() :
-		m_system(nullptr), m_resources(nullptr)
+		m_resources(nullptr)
 	{
 		std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
 	}
 	
 	MaskApp::~MaskApp() 
 	{
-		delete m_system;
 		for(std::size_t i=0; i<m_systemList.size(); i++)
 			delete m_systemList[i];
 	}
 
 	bool MaskApp::LoadConfig(const std::string& filename)
 	{
-		std::cout<<"Loading configuration in "<<filename<<"..."<<std::endl;
+		std::cout << "Loading configuration in " << filename << std::endl;
 		
-		std::ifstream input(filename);
-		if(!input.is_open()) 
+		if (!ConfigSerializer::DeserializeConfig(filename, m_params))
 		{
 			return false;
 		}
 
-		std::string junk;
-		std::getline(input, junk);
-		input>>junk>>m_outputName;
-		input>>junk>>m_nthreads;
-		std::getline(input, junk);
-		std::getline(input, junk);
-		input>>junk>>m_nsamples;
-
-		std::cout<<"Allocating resources... Asking for " << m_nthreads << " threads...";
-		m_resources = std::make_unique<ThreadPool<ReactionSystem*, uint64_t>>(m_nthreads);
-		std::cout<<" Complete."<<std::endl;
-
-		//Little bit of integer division mangling to make sure we do the total number of samples
-    	uint64_t quotient = m_nsamples / m_nthreads;
-    	uint64_t remainder = m_nsamples % m_nthreads;
-    	m_chunkSamples.push_back(quotient + remainder);
-    	for(uint64_t i=1; i<m_nthreads; i++)
-        	m_chunkSamples.push_back(quotient);
-
-		std::cout<<"Outputing data to file: " << m_outputName <<std::endl;
-		m_fileWriter.Open(m_outputName, "SimTree");
-
-		std::vector<StepParameters> params;
-		int z, a;
-		while(input>>junk)
+		//Initialize the system
+		//Reaction chain
+		for(uint32_t i=0; i<m_params.nThreads; i++)
 		{
-			if(junk == "begin_chain")
-				continue;
-			else if (junk == "end_chain")
-				break;
-			else if(junk == "begin_step")
-			{
-				StepParameters currentParams;
-				input >> junk >> junk;
-				currentParams.rxnType = StringToRxnType(junk);
-				if(currentParams.rxnType == RxnType::Reaction)
-				{
-					input >> junk;
-					for(int i=0; i<3; i++)
-					{
-						input >> z >> a;
-						currentParams.Z.push_back(z);
-						currentParams.A.push_back(a);
-					}
-					input >> junk;
-					input >> junk >> currentParams.meanBeamEnergy;
-					input >> junk >> currentParams.sigmaBeamEnergy;
-					input >> junk >> junk;
-					currentParams.thetaType = StringToRxnThetaType(junk);
-					input >> junk >> currentParams.thetaMin;
-					input >> junk >> currentParams.thetaMax;
-					input >> junk >> currentParams.phiMin;
-					input >> junk >> currentParams.phiMax;
-					input >> junk >> currentParams.meanResidualEx;
-					input >> junk >> currentParams.sigmaResidualEx;
-					params.push_back(currentParams);
-				}
-				else if(currentParams.rxnType == RxnType::Decay)
-				{
-					input >> junk;
-					for(int i=0; i<2; i++)
-					{
-						input >> z >> a;
-						currentParams.Z.push_back(z);
-						currentParams.A.push_back(a);
-					}
-					input >> junk;
-					input >> junk >> currentParams.phiMin;
-					input >> junk >> currentParams.phiMax;
-					input >> junk >> currentParams.meanResidualEx;
-					input >> junk >> currentParams.sigmaResidualEx;
-					input >> junk >> currentParams.angularDistFile;
-					params.push_back(currentParams);
-				}
-				else
-				{
-					std::cerr << "Invalid reaction information at MaskApp::LoadConfig!" << std::endl;
-					return false;
-				}
-			}
-		}
-
-		m_system = CreateSystem(params);
-		if(m_system == nullptr || !m_system->IsValid())
-		{
-			std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
-			return false;
-		}
-
-		for(uint32_t i=0; i<m_nthreads; i++)
-		{
-			m_systemList.push_back(CreateSystem(params));
+			m_systemList.push_back(CreateSystem(m_params.chainParams));
 			if(m_systemList.back() == nullptr || !m_systemList.back()->IsValid())
 			{
 				std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
 				return false;
 			}
 		}
-
-		std::getline(input, junk);
-		std::getline(input, junk);
-
-		LayeredTarget target;
-		int nlayers;
-		double thickness;
-		std::vector<int> avec, zvec, svec;
-		int s;
-		input>>junk>>nlayers;
-		for(int i=0; i<nlayers; i++) 
-		{
-			input>>junk>>junk>>thickness;
-			avec.clear(); zvec.clear(); svec.clear();
-			while(input>>junk) 
-			{
-				if(junk == "begin_elements")
-				{
-					input>>junk>>junk>>junk;
-					continue;
-				} 
-				else if (junk == "end_elements")
-					break;
-				input>>z>>a>>s;
-				zvec.push_back(z); avec.push_back(a); svec.push_back(s);
-			}
-			target.AddLayer(zvec, avec, svec, thickness);
-			input>>junk;
-		}
-		m_system->SetLayeredTarget(target);
-
+		//Link the target
 		for(auto system : m_systemList)
 		{
-			system->SetLayeredTarget(target);
+			system->SetLayeredTarget(m_params.target);
 		}
+		//Setup threading
+		m_resources = std::make_unique<ThreadPool<ReactionSystem*, uint64_t>>(m_params.nThreads);
+		//Little bit of integer division mangling to make sure we do the total number of samples
+    	uint64_t quotient = m_params.nSamples / m_params.nThreads;
+    	uint64_t remainder = m_params.nSamples % m_params.nThreads;
+    	m_chunkSamples.push_back(quotient + remainder);
+    	for(uint64_t i=1; i<m_params.nThreads; i++)
+        	m_chunkSamples.push_back(quotient);
+		m_fileWriter.Open(m_params.outputFileName, "SimTree");
 
-		std::cout<<"Reaction equation: "<<m_system->GetSystemEquation()<<std::endl;
-		std::cout<<"Number of samples: "<<m_nsamples<<std::endl;
 
+		std::cout << "Reaction equation: " << m_systemList[0]->GetSystemEquation() << std::endl;
+		std::cout << "Number of samples: " << m_params.nSamples << std::endl;
+		std::cout << "Number of threads: " << m_params.nThreads << std::endl;
+		std::cout << "Outputing data to file: " << m_params.outputFileName << std::endl;
 		return true;
 	}
 
-	//Not implemented... yet!
 	bool MaskApp::SaveConfig(const std::string& filename)
 	{
-		return true;
+		std::cout << "Writing configuration to " << filename << std::endl;
+		return ConfigSerializer::SerializeConfig(filename, m_params);
 	}
 	
 	void MaskApp::Run()
 	{
 		std::cout<<"Running simulation..."<<std::endl;
-		if(m_systemList.size() != m_nthreads)
+		if(m_systemList.size() != m_params.nThreads)
 		{
 			std::cerr << "System list not equal to number of threads" << std::endl;
 			return;
@@ -205,7 +100,7 @@ namespace Mask {
 
 		uint64_t count = 0;
 		double percent = 0.05;
-		uint64_t flushVal = m_nsamples*percent;
+		uint64_t flushVal = m_params.nSamples*percent;
 		uint64_t flushCount = 0;
 		while(true)
 		{
@@ -227,42 +122,4 @@ namespace Mask {
 		std::cout<<"---------------------------------------------"<<std::endl;
 	}
 	
-	void MaskApp::RunSingleThread()
-	{
-		std::cout<<"Running simulation..."<<std::endl;
-		if(m_system == nullptr) 
-		{
-			return;
-		}
-	
-		TFile* output = TFile::Open(m_outputName.c_str(), "RECREATE");
-		TTree* tree = new TTree("SimTree", "SimTree");
-		tree->Branch("nuclei", m_system->GetNuclei());
-	
-		//For progress tracking
-		uint64_t percent5 = 0.05*m_nsamples;
-		uint64_t count = 0;
-		uint64_t npercent = 0;
-	
-		for(uint64_t i=0; i<m_nsamples; i++) 
-		{
-			if(++count == percent5) 
-			{
-				npercent++;
-				count = 0;
-				std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
-			}
-	
-			m_system->RunSystem();
-			tree->Fill();
-		}
-	
-		tree->Write(tree->GetName(), TObject::kOverwrite);
-		output->Close();
-		
-		std::cout<<std::endl;
-		std::cout<<"Complete."<<std::endl;
-		std::cout<<"---------------------------------------------"<<std::endl;
-	}
-
 }

src/Mask/MaskApp.h

@@ -14,6 +14,15 @@
 
 namespace Mask {
 
+	struct AppParameters
+	{
+		std::string outputFileName = "";
+		uint64_t nSamples = 0;
+		uint32_t nThreads = 1;
+		std::vector<StepParameters> chainParams;
+		LayeredTarget target;
+	};
+
 	class MaskApp
 	{
 	public:
@@ -23,17 +32,9 @@ namespace Mask {
 		bool SaveConfig(const std::string& filename);
 
 		void Run();
-		void RunSingleThread();
 
 	private:
-		void RunChunk(ReactionSystem* system);
-
-		ReactionSystem* m_system;
-
-		std::string m_outputName;
-		RxnType m_rxnType;
-		uint64_t m_nsamples;
-		uint32_t m_nthreads;
+		AppParameters m_params;
 
 		std::vector<ReactionSystem*> m_systemList; //One system for each thread
 		std::vector<uint64_t> m_chunkSamples;

src/Mask/Nucleus.h

@@ -12,6 +12,7 @@
 #include <string>
 #include <vector>
 #include "Math/Vector4D.h"
+#include "Math/Point3D.h"
 #include "MassLookup.h"
 
 namespace Mask {
@@ -48,6 +49,8 @@ namespace Mask {
 		double detectedKE = 0.0;
 		double detectedTheta = 0.0;
 		double detectedPhi = 0.0;
+
+		ROOT::Math::XYZPoint detectedPos = ROOT::Math::XYZPoint(0., 0., 0.);
 	};
 
 	Nucleus CreateNucleus(uint32_t z, uint32_t a);

src/Mask/Reaction.cpp

@@ -229,7 +229,6 @@ namespace Mask {
 		m_residual->vec4 = m_target->vec4 - m_ejectile->vec4;
 	
 		//energy loss for the *light* break up nucleus
-		double keorig = m_ejectile->GetKE();
 		double ejectKE = m_ejectile->GetKE() - 
 					m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, m_ejectile->GetKE(), m_rxnLayer, m_ejectile->vec4.Theta(), m_rxnDepth);
 		if(ejectKE > 0.0)
@@ -257,6 +256,27 @@ namespace Mask {
 		}
 	}
 
+	bool Reaction::CheckReactionThreshold(double beamEnergy, double residEx)
+	{
+		double Q = m_target->groundStateMass + m_projectile->groundStateMass -
+				   (m_ejectile->groundStateMass + m_residual->groundStateMass + residEx);
+		double Ethresh = -Q*(m_ejectile->groundStateMass+m_residual->groundStateMass) / 
+						 (m_ejectile->groundStateMass + m_residual->groundStateMass - m_projectile->groundStateMass);
+		if (beamEnergy < Ethresh)
+			return false;
+		else
+			return true;
+	}
+
+	bool Reaction::CheckDecayThreshold(double targetEx, double residEx)
+	{
+		double Q = m_target->groundStateMass + targetEx - 
+				  (m_ejectile->groundStateMass + m_residual->groundStateMass + residEx);
+		if ( Q < 0.0)
+			return false;
+		else
+			return true;
+	}
 }
 
 

src/Mask/Reaction.h

@@ -42,6 +42,9 @@ namespace Mask {
 		void SetRxnLayer(std::size_t layer) { m_rxnLayer = layer; };
 		void SetResidualEnergyLoss(bool isEloss) { m_isResidEloss = isEloss; };
 
+		bool CheckReactionThreshold(double beamEnergy, double residualEx);
+		bool CheckDecayThreshold(double targetEx, double residualEx);
+
 		bool IsDecay() const { return m_isDecay; };
 
 		std::size_t GetRxnLayer() const { return m_rxnLayer; };

src/Mask/RxnType.h

@@ -37,6 +37,17 @@ namespace Mask {
 			return RxnType::None;
 	}
 
+	static std::string RxnTypeToString(RxnType type)
+	{
+		switch(type)
+		{
+			case RxnType::Decay: return "Decay";
+			case RxnType::Reaction: return "Reaction";
+			case RxnType::None: return "None";
+			default: return "None";
+		}
+	}
+
 	static RxnThetaType StringToRxnThetaType(const std::string& type)
 	{
 		if(type == "CenterOfMass")
@@ -52,6 +63,17 @@ namespace Mask {
 			return RxnThetaType::None;
 		}
 	}
+
+	static std::string RxnThetaTypeToString(RxnThetaType type)
+	{
+		switch(type)
+		{
+			case RxnThetaType::CenterOfMass: return "CenterOfMass";
+			case RxnThetaType::Lab: return "Lab";
+			case RxnThetaType::None: return "None";
+			default: return "None";
+		}
+	}
 }
 
 #endif

src/Mask/Target.h

@@ -31,11 +31,11 @@ namespace Mask {
 	 	double GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle);
 	 	double GetEnergyLossFractionalDepth(int zp, int ap, double startEnergy, double angle, double percent_depth);
 	 	double GetReverseEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double angle, double percent_depth);
-	 	inline const double& GetThickness() { return m_thickness; }
-	 	inline int GetNumberOfElements() { return m_Z.size(); }
-	 	inline int GetElementZ(int index) { return m_Z[index]; }
-	 	inline int GetElementA(int index) { return m_A[index]; }
-	 	inline int GetElementStoich(int index) { return m_stoich[index]; }
+	 	const double& GetThickness() const { return m_thickness; }
+	 	int GetNumberOfElements() const { return m_Z.size(); }
+	 	int GetElementZ(int index) const { return m_Z[index]; }
+	 	int GetElementA(int index) const { return m_A[index]; }
+	 	int GetElementStoich(int index) const { return m_stoich[index]; }
 
 	private:
 	 	void Init(const std::vector<int>& z, const std::vector<int>& a, const std::vector<int>& stoich);

src/Mask/ThreeStepSystem.cpp

@@ -2,6 +2,11 @@
 #include "RandomGenerator.h"
 #include "KinematicsExceptions.h"
 
+#include "Math/Boost.h"
+#include "Math/Vector3D.h"
+#include "Math/RotationZ.h"
+#include "Math/RotationY.h"
+
 namespace Mask {
 	
 	ThreeStepSystem::ThreeStepSystem(const std::vector<StepParameters>& params) :
@@ -56,8 +61,8 @@ namespace Mask {
 		m_nuclei[3] = CreateNucleus(zr, ar); //residual
 		m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
 		m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
-		m_nuclei[5] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
-		m_nuclei[5] = CreateNucleus(zb4, ab4); //breakup4
+		m_nuclei[6] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
+		m_nuclei[7] = CreateNucleus(zb4, ab4); //breakup4
 
 		m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
 		m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
@@ -115,43 +120,57 @@ namespace Mask {
 		m_sysEquation = stream.str();
 	}
 
+	ThreeStepParameters ThreeStepSystem::SampleParameters()
+	{
+		ThreeStepParameters params;
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		params.beamEnergy = m_beamDistributions[0](gen);
+		params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
+		params.rxnPhi = m_phiRanges[0](gen);
+		params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
+		params.decay1Theta = std::acos(params.cosdecay1Theta);
+		params.decay1Phi = m_phiRanges[1](gen);
+		params.cosdecay2Theta = m_decayAngularDistributions[1].GetRandomCosTheta();
+		params.decay2Theta = std::acos(params.cosdecay2Theta);
+		params.decay2Phi = m_phiRanges[1](gen);
+		params.residEx = m_exDistributions[0](gen);
+		params.decay1Ex = m_exDistributions[1](gen);
+		params.decay2Ex = m_exDistributions[2](gen);
+		params.rxnDepth = (m_rxnDepthDist(gen));
+		return params;
+	}
+	
 	void ThreeStepSystem::RunSystem()
 	{
-		//Sample parameters
-		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
-		double bke = m_beamDistributions[0](gen);
-		double rxnTheta = std::acos((m_thetaRanges[0])(gen));
-		double rxnPhi = m_phiRanges[0](gen);
-		double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
-		double decay1Theta = std::acos(decay1costheta);
-		double decay1Phi = m_phiRanges[1](gen);
-		double decay2costheta = m_decayAngularDistributions[1].GetRandomCosTheta();
-		double decay2Theta = std::acos(decay2costheta);
-		double decay2Phi = m_phiRanges[2](gen);
-		double residEx = m_exDistributions[0](gen);
-		double decay1Ex = m_exDistributions[1](gen);
-		double decay2Ex = m_exDistributions[2](gen);
-		double rxnDepth = (m_rxnDepthDist(gen));
+		ThreeStepParameters params;
+		do
+		{
+			params = SampleParameters();
+		}
+		while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) 
+				&& m_step2.CheckDecayThreshold(params.residEx, params.decay1Ex)
+				&& m_step3.CheckDecayThreshold(params.decay1Ex, params.decay2Ex)));
 		
-		m_step1.SetReactionDepth(rxnDepth);
-		m_step1.SetBeamKE(bke);
-		m_step1.SetPolarRxnAngle(rxnTheta);
-		m_step1.SetAzimRxnAngle(rxnPhi);
-		m_step1.SetExcitation(residEx);
+		m_step1.SetReactionDepth(params.rxnDepth);
+		m_step1.SetBeamKE(params.beamEnergy);
+		m_step1.SetPolarRxnAngle(params.rxnTheta);
+		m_step1.SetAzimRxnAngle(params.rxnPhi);
+		m_step1.SetExcitation(params.residEx);
 	
-		m_step2.SetReactionDepth(rxnDepth);
-		m_step2.SetPolarRxnAngle(decay1Theta);
-		m_step2.SetAzimRxnAngle(decay1Phi);
-		m_step2.SetExcitation(decay1Ex);
+		m_step2.SetReactionDepth(params.rxnDepth);
+		m_step2.SetPolarRxnAngle(params.decay1Theta);
+		m_step2.SetAzimRxnAngle(params.decay1Phi);
+		m_step2.SetExcitation(params.decay1Ex);
 	
-		m_step3.SetReactionDepth(rxnDepth);
-		m_step3.SetPolarRxnAngle(decay2Theta);
-		m_step3.SetAzimRxnAngle(decay2Phi);
-		m_step3.SetExcitation(decay2Ex);
+		m_step3.SetReactionDepth(params.rxnDepth);
+		m_step3.SetPolarRxnAngle(params.decay2Theta);
+		m_step3.SetAzimRxnAngle(params.decay2Phi);
+		m_step3.SetResidualEnergyLoss(true);
+		m_step3.SetExcitation(params.decay2Ex);
 	
 		m_step1.Calculate();
 	
-		if(decay1costheta == -10)
+		if(params.cosdecay1Theta == -10)
 		{
 			m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
 			m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
@@ -161,15 +180,13 @@ namespace Mask {
 		}
 		m_step2.Calculate();
 	
-		if(decay2costheta == -10)
+		if(params.cosdecay2Theta == -10)
 		{
 			m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
 			m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
 			return;
 		}
-		m_step3.SetResidualEnergyLoss(true);
 		m_step3.Calculate();
-	
 	}
 
 }

src/Mask/ThreeStepSystem.h

@@ -6,6 +6,23 @@
 
 namespace Mask {
 
+	struct ThreeStepParameters
+	{
+		double beamEnergy = 0.;
+		double rxnTheta = 0.;
+		double rxnPhi = 0.;
+		double cosdecay1Theta = 0.;
+		double cosdecay2Theta = 0.;
+		double decay1Theta = 0.;
+		double decay1Phi = 0.;
+		double decay2Theta = 0.;
+		double decay2Phi = 0.;
+		double residEx = 0.;
+		double decay1Ex = 0.;
+		double decay2Ex = 0.;
+		double rxnDepth = 0.;
+	};
+
 	class ThreeStepSystem : public ReactionSystem
 	{
 	public:
@@ -18,6 +35,7 @@ namespace Mask {
 	protected:
 		void Init(const std::vector<StepParameters>& params);
 		void SetSystemEquation() override;
+		ThreeStepParameters SampleParameters();
 	
 		Reaction m_step1, m_step2, m_step3;
 	

src/Mask/TwoStepSystem.cpp

@@ -93,34 +93,58 @@ namespace Mask {
 		m_sysEquation = stream.str();
 	}
 
+	TwoStepParameters TwoStepSystem::SampleParameters()
+	{
+		TwoStepParameters params;
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		params.beamEnergy = (m_beamDistributions[0])(gen);
+		params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
+		params.rxnPhi = (m_phiRanges[0])(gen);
+		params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
+		params.decay1Theta = std::acos(params.cosdecay1Theta);
+		params.decay1Phi = m_phiRanges[1](gen);
+		params.residEx = (m_exDistributions[0])(gen);
+		params.decay2Ex = m_exDistributions[1](gen);
+		params.rxnDepth = (m_rxnDepthDist(gen));
+		return params;
+	}
+	
 	void TwoStepSystem::RunSystem()
 	{
 		//Sample parameters
-		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
-		double bke = (m_beamDistributions[0])(gen);
-		double rxnTheta = std::acos((m_thetaRanges[0])(gen));
-		double rxnPhi = (m_phiRanges[0])(gen);
-		double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
-		double decay1Theta = std::acos(decay1costheta);
-		double decay1Phi = m_phiRanges[1](gen);
-		double residEx = (m_exDistributions[0])(gen);
-		double decay2Ex = m_exDistributions[1](gen);
-		double rxnDepth = (m_rxnDepthDist(gen));
+		// std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		// double bke = (m_beamDistributions[0])(gen);
+		// double rxnTheta = std::acos((m_thetaRanges[0])(gen));
+		// double rxnPhi = (m_phiRanges[0])(gen);
+		// double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
+		// double decay1Theta = std::acos(decay1costheta);
+		// double decay1Phi = m_phiRanges[1](gen);
+		// double residEx = (m_exDistributions[0])(gen);
+		// double decay2Ex = m_exDistributions[1](gen);
+		// double rxnDepth = (m_rxnDepthDist(gen));
 
-		m_step1.SetReactionDepth(rxnDepth);
-		m_step1.SetBeamKE(bke);
-		m_step1.SetPolarRxnAngle(rxnTheta);
-		m_step1.SetAzimRxnAngle(rxnPhi);
-		m_step1.SetExcitation(residEx);
+		TwoStepParameters params;
+		do
+		{
+			params = SampleParameters();
+		}
+		while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) 
+				&& m_step2.CheckDecayThreshold(params.residEx, params.decay2Ex)));
 		
-		m_step2.SetReactionDepth(rxnDepth);
-		m_step2.SetPolarRxnAngle(decay1Theta);
-		m_step2.SetAzimRxnAngle(decay1Phi);
-		m_step2.SetExcitation(decay2Ex);
+		m_step1.SetReactionDepth(params.rxnDepth);
+		m_step1.SetBeamKE(params.beamEnergy);
+		m_step1.SetPolarRxnAngle(params.rxnTheta);
+		m_step1.SetAzimRxnAngle(params.rxnPhi);
+		m_step1.SetExcitation(params.residEx);
+	
+		m_step2.SetReactionDepth(params.rxnDepth);
+		m_step2.SetPolarRxnAngle(params.decay1Theta);
+		m_step2.SetAzimRxnAngle(params.decay1Phi);
+		m_step2.SetExcitation(params.decay2Ex);
 		
 		m_step1.Calculate();
 	
-		if(decay1costheta == -10)
+		if(params.cosdecay1Theta == -10)
 		{
 			m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
 			m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);

src/Mask/TwoStepSystem.h

@@ -6,6 +6,19 @@
 
 namespace Mask {
 
+	struct TwoStepParameters
+	{
+		double beamEnergy = 0.;
+		double rxnTheta = 0.;
+		double rxnPhi = 0.;
+		double cosdecay1Theta = 0.;
+		double decay1Theta = 0.;
+		double decay1Phi = 0.;
+		double residEx = 0.;
+		double decay2Ex = 0.;
+		double rxnDepth = 0.;
+	};
+
 	class TwoStepSystem : public ReactionSystem
 	{
 	public:
@@ -18,6 +31,7 @@ namespace Mask {
 	private:
 		void Init(const std::vector<StepParameters>& params);
 		void SetSystemEquation() override;
+		TwoStepParameters SampleParameters();
 
 		Reaction m_step1, m_step2;
 	};

src/Plotters/RootPlotter.cpp

@@ -1,7 +1,8 @@
 #include "RootPlotter.h"
 #include <TFile.h>
 #include <TTree.h>
-#include <TVector3.h>
+#include <Math/Vector3D.h>
+#include "Math/Boost.h"
 
 #include <iostream>
 
@@ -10,7 +11,8 @@ static double FullPhi(double phi)
 	return phi < 0.0 ? 2.0*M_PI + phi : phi; 
 }
 
-RootPlotter::RootPlotter()
+RootPlotter::RootPlotter() :
+	m_target({5},{9},{1}, 74.0)
 {
 	TH1::AddDirectory(kFALSE);
 }
@@ -42,10 +44,16 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam
 			std::cout<<"\rPercent of data processed: "<<flushCount*flushFrac*100<<"%"<<std::flush;
 		}
 		tree->GetEntry(i);
-		for(Mask::Nucleus& nuc : *(dataHandle))
+		// for(Mask::Nucleus& nuc : *(dataHandle))
+		// {
+		// 	FillData(nuc);
+		// }
+		for(int i=0; i<dataHandle->size(); i++)
 		{
-			FillData(nuc);
+			FillData(dataHandle->at(i), i);
 		}
+		//Don't leave this in!
+		//Correlations(*dataHandle);
 	}
 	std::cout<<std::endl;
 	input->Close();
@@ -57,40 +65,142 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam
 	output->Close();
 }
 
-void RootPlotter::FillData(const Mask::Nucleus& nuc)
+void RootPlotter::FillData(const Mask::Nucleus& nuc, int i)
 {
-	std::string modifier = "";
-	if(nuc.isDetected)
-		modifier = "detected";
+	std::string mod = "_detected";
+	std::string num = "_" + std::to_string(i);
 
 	std::string sym = nuc.isotopicSymbol;
-	std::string ke_vs_th_name = sym + modifier + "_ke_vs_theta";
+	std::string ke_vs_th_name = sym + num + "_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_th_name_det = sym + num + mod + "_ke_vs_theta";
+	std::string ke_vs_th_title_det = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
+	std::string ke_vs_ph_name = sym + num + "_ke_vs_phi";
 	std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
-	std::string th_vs_ph_name = sym + modifier + "_theta_vs_phi";
+	std::string ke_vs_ph_name_det = sym + num + mod + "_ke_vs_phi";
+	std::string ke_vs_ph_title_det = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
+	std::string th_vs_ph_name = sym + num + "_theta_vs_phi";
 	std::string th_vs_ph_title = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
-	std::string ex_name = sym + modifier + "_ex";
+	std::string th_vs_ph_name_det = sym + num + mod + "_theta_vs_phi";
+	std::string th_vs_ph_title_det = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
+	std::string ex_name = sym + num + "_ex";
 	std::string ex_title = ex_name + ";E_{ex} (MeV);counts";
-	std::string angdist_name = sym + modifier +"_angDist";
+	std::string ex_name_det = sym + num + mod + "_ex";
+	std::string ex_title_det = ex_name + ";E_{ex} (MeV);counts";
+	std::string angdist_name = sym + num + "_angDist";
 	std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts";
+	std::string angdist_name_det = sym + num +  mod +"_angDist";
+	std::string angdist_title_det = angdist_name+";cos#right(#theta_{CM}#left);counts";
+	std::string hist_ke_th_name = sym + num + "_hist_ke_vs_theta";
+	std::string hist_ke_th_title = hist_ke_th_name + ";#theta_{lab};Kinetic Energy (MeV)";
+	std::string hist_ke_th_name_det = sym + num + mod + "_hist_ke_vs_theta";
+	std::string hist_ke_th_title_det = hist_ke_th_name + ";#theta_{lab};Kinetic Energy (MeV)";
 	
 	MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
 	MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.GetKE(), 4);
 	MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
 	MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
 	MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
-	if(nuc.isDetected)
+	MyFill(hist_ke_th_name, hist_ke_th_title, 180, 0.0, 180.0, 400, 0.0, 20.0, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE());
+	if(nuc.isDetected && nuc.detectedKE > 0.25)
 	{
-		MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
-		MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4);
-		MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
-		MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
-		MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
+		MyFill(ke_vs_th_name_det, ke_vs_th_title_det, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
+		MyFill(ke_vs_ph_name_det, ke_vs_ph_title_det, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4);
+		MyFill(th_vs_ph_name_det, th_vs_ph_title_det, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
+		MyFill(ex_name_det, ex_title_det, 260, -1.0, 25, nuc.GetExcitationEnergy());
+		MyFill(angdist_name_det, angdist_title_det, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
+		MyFill(hist_ke_th_name_det, hist_ke_th_title_det, 180, 0.0, 180.0, 400, 0.0, 20.0, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE);
 	}
 	
 }
 
+//Correlation analysis for coupled step reactions, will need to be customized
+//for a given reaction. By default it is off.
+void RootPlotter::Correlations(const std::vector<Mask::Nucleus>& event)
+{
+	static constexpr double lowerBoundEx = -2.7;
+    static constexpr double upperBoundEx = 2.7;
+	static constexpr double targetDepth = 0.5;
+	const Mask::Nucleus& secondary = event[6];
+	const Mask::Nucleus& primary = event[4];
+	const Mask::Nucleus& parent = event[3];
+    const Mask::Nucleus& li = event[5];
+	if(secondary.vec4.P() == 0.0 || primary.vec4.P() == 0.0)
+    {
+        return;
+    }
+	ROOT::Math::Boost boostParent(parent.vec4.BoostToCM());
+	ROOT::Math::Boost boostIntermediate(li.vec4.BoostToCM());
+	ROOT::Math::PxPyPzEVector a2Vec = (secondary.vec4);
+    ROOT::Math::PxPyPzEVector a1Vec = (primary.vec4);
+
+	double a2KE = secondary.GetKE() + m_target.GetReverseEnergyLossFractionalDepth(secondary.Z, secondary.A, secondary.GetKE(), a2Vec.Theta(), targetDepth);
+	double a2P = std::sqrt(a2KE * (a2KE + 2.0 * secondary.groundStateMass));
+	a2Vec.SetPxPyPzE(
+		a2P * std::sin(a2Vec.Theta()) * std::cos(a2Vec.Phi()),
+		a2P * std::sin(a2Vec.Theta()) * std::sin(a2Vec.Phi()),
+		a2P * std::cos(a2Vec.Theta()),
+		a2KE + secondary.groundStateMass
+	);
+
+	double a1KE = primary.GetKE() + m_target.GetReverseEnergyLossFractionalDepth(primary.Z, primary.A, primary.GetKE(), a1Vec.Theta(), targetDepth);
+	double a1P = std::sqrt(a1KE * (a1KE + 2.0 * primary.groundStateMass));
+	a1Vec.SetPxPyPzE(
+		a1P * std::sin(a1Vec.Theta()) * std::cos(a1Vec.Phi()),
+		a1P * std::sin(a1Vec.Theta()) * std::sin(a1Vec.Phi()),
+		a1P * std::cos(a1Vec.Theta()),
+		a1KE + primary.groundStateMass
+	);
+
+	ROOT::Math::XYZVector a1PVec;
+	ROOT::Math::XYZVector a2PVec;
+
+	ROOT::Math::PxPyPzEVector boostedIntA1Vec = boostIntermediate * a1Vec;
+	ROOT::Math::PxPyPzEVector boostedIntA2Vec = boostIntermediate * a2Vec;
+	ROOT::Math::PxPyPzEVector boostedParA2Vec = boostParent * a2Vec;
+    a1PVec.SetXYZ(boostedIntA1Vec.Px(), boostedIntA1Vec.Py(), boostedIntA1Vec.Pz());
+    a2PVec.SetXYZ(boostedIntA2Vec.Px(), boostedIntA2Vec.Py(), boostedIntA2Vec.Pz());
+    double relCosThetaCM = a1PVec.Dot(a2PVec) / (a1PVec.R() * a2PVec.R());
+
+    ROOT::Math::PxPyPzEVector badLi = parent.vec4 - a2Vec;
+	double secondaryExcitation = badLi.M() - li.groundStateMass;
+
+	MyFill("seoncdaryEx","secondaryEx",1000,-5.0,5.0,secondaryExcitation);
+	MyFill("primaryEx","primaryEx",1000,-5.0,5.0,li.GetExcitationEnergy());
+	MyFill("relativeDist_full","relativeDist_full",20,-1.0,1.0, relCosThetaCM);
+	MyFill("primaryDist_full","primaryDist_full",20,-1.0,1.0,std::cos(primary.thetaCM));
+	MyFill("seondaryDist_full","seconaryDist_full",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
+	MyFill("secondaryDist_fullOwn","seondaryDist_fullOwn",20,-1.0,1.0,std::cos(secondary.thetaCM));
+
+    bool isPrimary = false;
+    bool isSecondary = false;
+    if(primary.isDetected)
+    {
+		MyFill("primaryExMatched","primaryExMatched",1000,-5.0,5.0,li.GetExcitationEnergy());
+		MyFill("primaryDist","primaryDist",20,-1.0,1.0,std::cos(primary.thetaCM));
+        isPrimary = true;
+    }
+    if( secondary.isDetected && secondaryExcitation > lowerBoundEx && secondaryExcitation < upperBoundEx )
+    {
+		MyFill("secondaryDist","secondaryDist",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
+		MyFill("secondaryExMatched","secondaryExMatched",1000,-5.0,5.0,secondaryExcitation);
+        isSecondary = true;
+    }
+
+    if(isSecondary && isPrimary)
+    {
+		return;
+    }
+    else if (isSecondary)
+    {
+        MyFill("completeDist","completeDist",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
+    }
+    else if (isPrimary)
+    {
+        MyFill("completeDist","completeDist",20,-1.0,1.0,std::cos(primary.thetaCM));
+    }
+}
+
 void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val)
 {
 	auto iter = m_map.find(name);

src/Plotters/RootPlotter.h

@@ -6,6 +6,7 @@
 #include <unordered_map>
 
 #include "Nucleus.h"
+#include "Target.h"
 
 #include <TH1.h>
 #include <TH2.h>
@@ -20,7 +21,8 @@ public:
 	void Run(const std::string& inputname, const std::string& outputname);
 	
 private:
-	void FillData(const Mask::Nucleus& nuc);
+	void FillData(const Mask::Nucleus& nuc, int i);
+	void Correlations(const std::vector<Mask::Nucleus>& event);
 	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);
@@ -29,6 +31,7 @@ private:
 	std::unordered_map<std::string, std::shared_ptr<TObject>> m_map;
 
 	static constexpr double s_rad2deg = 180.0/M_PI;
+	Mask::Target m_target;
 };
 
 #endif

src/vendor/yaml-cpp

@@ -0,0 +1 @@
+Subproject commit 0e6e28d1a38224fc8172fae0109ea7f673c096db