Switch to yaml for configuration files, add yaml-cpp as submodule

.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)

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)
 {
@@ -13,6 +15,65 @@ DetectorApp::~DetectorApp()
         delete m_detectorList[i];
 }
 
+bool DetectorApp::ParseConfig(const std::string& filename)
+{
+    std::cout<<"----------Detector Efficiency Calculation----------"<<std::endl;
+    YAML::Node data;
+    try
+    {
+        data = YAML::LoadFile(filename);
+    }
+    catch (YAML::ParserException& e)
+    {
+        std::cerr << "Could not load config file " << filename  << " with error: " << e.what() << std::endl;
+        return false;
+    }
+
+    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>());
+
+    std::cout << "Creating " << m_nthreads << " detector arrays..." << std::endl;
+    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 << "Allocating " << m_nthreads << " threads..." << std::endl;
+    m_resources = std::make_unique<Mask::ThreadPool<DetectorArray*, uint64_t>>(m_nthreads);
+
+    std::cout << "Opening input data file " << m_inputFileName << "..." << std::endl;
+    m_fileReader.Open(m_inputFileName, "SimTree");
+    if(!m_fileReader.IsOpen() || !m_fileReader.IsTree())
+    {
+        std::cerr << "Unable to open input data file " << m_inputFileName << std::endl;
+        return false;
+    }
+    m_nentries = m_fileReader.GetSize();
+    std::cout << "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;
+    uint64_t remainder = m_nentries % m_nthreads;
+    m_chunkSamples.push_back(quotient + remainder);
+    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;
+    }
+
+    return true;
+}
+
 bool DetectorApp::LoadConfig(const std::string& filename)
 {
     std::cout<<"----------Detector Efficiency Calculation----------"<<std::endl;

src/Detectors/DetectorApp.h

@@ -17,6 +17,7 @@ public:
     ~DetectorApp();
 
     bool LoadConfig(const std::string& filename);
+    bool ParseConfig(const std::string& filename);
     
     void Run();
 

src/Detectors/main.cpp

@@ -21,7 +21,7 @@ int main(int argc, char** argv)
 	try
 	{
 		DetectorApp app;
-		if(!app.LoadConfig(argv[1]))
+		if(!app.ParseConfig(argv[1]))
 		{
 			std::cerr << "Unable to load config file " << argv[1] << ". Shutting down." << std::endl;
 			return 1;

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

@@ -40,7 +40,7 @@ namespace Mask {
 	
 		if(!input.is_open())
 		{
-			std::cerr<<"Unable to open distribution file. All values reset to default."<<std::endl;
+			std::cerr << "Unable to open distribution file. All values reset to default." << std::endl;
 			m_L=0;
 			m_branchingRatio=1.0;
 			m_constants.clear();
@@ -59,8 +59,8 @@ namespace Mask {
 	
 		if(m_constants.size() != ((std::size_t) l+1))
 		{
-			std::cerr<<"Unexpected number of constants for given angular momentum! Expected "<<l+1<<" and given "<<m_constants.size()<<std::endl;
+			std::cerr << "Unexpected number of constants for given angular momentum! Expected " << l+1 << " and given " << m_constants.size() << std::endl;
-			std::cerr<<"Setting all values to default."<<std::endl;
+			std::cerr << "Setting all values to default." << std::endl;
 			m_L=0;
 			m_branchingRatio=1.0;
 			m_constants.clear();
@@ -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_compile_definitions(Mask PRIVATE YAML_CPP_STATIC_DEFINE)
-target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES} Threads::Threads)
+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) 
+	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 (!ConfigSerializer::DeserializeConfig(filename, m_params))
-		if(!input.is_open()) 
 		{
 			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...";
+		//Initialize the system
-		m_resources = std::make_unique<ThreadPool<ReactionSystem*, uint64_t>>(m_nthreads);
+		//Reaction chain
-		std::cout<<" Complete."<<std::endl;
+		for(uint32_t i=0; i<m_params.nThreads; i++)
-
-		//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)
 		{
-			if(junk == "begin_chain")
+			m_systemList.push_back(CreateSystem(m_params.chainParams));
-				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));
 			if(m_systemList.back() == nullptr || !m_systemList.back()->IsValid())
 			{
 				std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
 				return false;
 			}
 		}
-
+		//Link the target
-		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);
-
 		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);
+		AppParameters m_params;
-
-		ReactionSystem* m_system;
-
-		std::string m_outputName;
-		RxnType m_rxnType;
-		uint64_t m_nsamples;
-		uint32_t m_nthreads;
 
 		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,12 +31,12 @@ 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; }
+	 	const double& GetThickness() const { return m_thickness; }
-	 	inline int GetNumberOfElements() { return m_Z.size(); }
+	 	int GetNumberOfElements() const { return m_Z.size(); }
-	 	inline int GetElementZ(int index) { return m_Z[index]; }
+	 	int GetElementZ(int index) const { return m_Z[index]; }
-	 	inline int GetElementA(int index) { return m_A[index]; }
+	 	int GetElementA(int index) const { return m_A[index]; }
-	 	inline int GetElementStoich(int index) { return m_stoich[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