Update README

.gitignore

@@ -19,6 +19,7 @@ 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;
 }
 

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,9 +73,6 @@ 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. 

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}
 )
 
@@ -57,10 +58,12 @@ target_sources(Mask PRIVATE
     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/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/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/vendor/yaml-cpp

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