Merge pull request #4 from gwm17/devel

Multi-threading simulation

README.md

@@ -34,11 +34,13 @@ Detector geometry is encoded using ROOT math libraries in the `src/Detectors` fo
 
 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.
 
+To run the geometry code, one needs to provide an input file containing the following: the path of a Mask kinematics data file, the path to which data should be written, the path to a file containing a list of dead channels (optional, if not used, write None for the path), the number of threads to be used by the thread pool, and a keyword for the array type (current options are Sabre or Anasen)
+
 To run Detectors use the format
 
-`./bin/Detectors <kinematics_datafile> <new_detection_datafile> <new_detection_statsfile>`
+`./bin/Detectors <input_file>`
 
-where the detection datafile contains all of the kinematics data as well as information about which particles are detected and the statsfile is a text file containing efficiency statistics.
+An example input file is provided with 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.

detector_input.txt

@@ -0,0 +1,5 @@
+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_input.txt

@@ -1,7 +1,8 @@
 ----------Data Information----------
 OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root
+NumberOfThreads: 6
 ----------Reaction Information----------
-NumberOfSamples: 100000
+NumberOfSamples: 1000000
 begin_chain
 	begin_step
 		Type: Reaction

src/Detectors/AnasenArray.cpp

@@ -230,7 +230,7 @@ double AnasenArray::RunConsistencyCheck()
 
 }
 
-DetectorResult AnasenArray::IsRing1(Mask::Nucleus& nucleus)
+DetectorResult AnasenArray::IsRing1(const Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	double thetaIncident;
@@ -255,7 +255,7 @@ DetectorResult AnasenArray::IsRing1(Mask::Nucleus& nucleus)
 	return observation;
 }
 
-DetectorResult AnasenArray::IsRing2(Mask::Nucleus& nucleus)
+DetectorResult AnasenArray::IsRing2(const Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	double thetaIncident;
@@ -280,7 +280,7 @@ DetectorResult AnasenArray::IsRing2(Mask::Nucleus& nucleus)
 	return observation;
 }
 
-DetectorResult AnasenArray::IsQQQ(Mask::Nucleus& nucleus)
+DetectorResult AnasenArray::IsQQQ(const Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	double thetaIncident;
@@ -324,7 +324,7 @@ DetectorResult AnasenArray::IsQQQ(Mask::Nucleus& nucleus)
 	return observation;
 }
 
-DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus)
+DetectorResult AnasenArray::IsDetected(const Mask::Nucleus& nucleus)
 {
 	DetectorResult result;
 	if(nucleus.GetKE() <= s_energyThreshold)
@@ -338,235 +338,3 @@ DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus)
 		result = IsQQQ(nucleus);
 	return result;
 }
-
-void AnasenArray::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
-{
-	if (data.size() == 3 && data[1].isDetected && data[2].isDetected)
-	{
-		counts[0]++;
-	}
-	else if (data.size() == 4 && data[2].isDetected && data[3].isDetected)
-	{
-		counts[0]++;
-	}
-	else if(data.size() == 6)
-	{
-		if(data[2].isDetected && data[4].isDetected) 
-		{
-			counts[0]++;
-		}
-		if(data[2].isDetected && data[5].isDetected)
-		{
-			counts[1]++;
-		}
-		if(data[4].isDetected && data[5].isDetected)
-		{
-			counts[2]++;
-		}
-		if(data[2].isDetected && data[4].isDetected && data[5].isDetected)
-		{
-			counts[3]++;
-		}
-	}
-	else if(data.size() == 8)
-	{
-		if(data[2].isDetected && data[4].isDetected) 
-		{
-			counts[0]++;
-		}
-		if(data[2].isDetected && data[6].isDetected)
-		{
-			counts[1]++;
-		}
-		if(data[2].isDetected && data[7].isDetected)
-		{
-			counts[2]++;
-		}
-		if(data[4].isDetected && data[6].isDetected)
-		{
-			counts[3]++;
-		}
-		if(data[4].isDetected && data[7].isDetected)
-		{
-			counts[4]++;
-		}
-		if(data[6].isDetected && data[7].isDetected)
-		{
-			counts[5]++;
-		}
-		if(data[2].isDetected && data[4].isDetected && data[6].isDetected)
-		{
-			counts[6]++;
-		}
-		if(data[2].isDetected && data[4].isDetected && data[7].isDetected)
-		{
-			counts[7]++;
-		}
-		if(data[2].isDetected && data[6].isDetected && data[7].isDetected)
-		{
-			counts[8]++;
-		}
-		if(data[4].isDetected && data[6].isDetected && data[7].isDetected)
-		{
-			counts[9]++;
-		}
-		if(data[2].isDetected && data[4].isDetected && data[6].isDetected && data[7].isDetected)
-		{
-			counts[10]++;
-		}
-	}
-}
-
-void AnasenArray::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) {
-
-	if(!dmap.IsValid())
-	{
-		std::cerr<<"Invalid Dead Channel Map at AnasenArray::CalculateEfficiency()! If you want to run with all possible";
-		std::cerr<<"channels active, simply load an empty file."<<std::endl;
-		return;
-	}
-
-	std::cout<<"----------ANASEN Efficiency Calculation----------"<<std::endl;
-	std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
-	std::cout<<"Running efficiency calculation..."<<std::endl;
-
-	TFile* input = TFile::Open(inputname.c_str(), "READ");
-	TFile* output = TFile::Open(outputname.c_str(), "RECREATE");
-	std::ofstream stats(statsname);
-	stats<<std::setprecision(5);
-
-	TTree* intree = (TTree*) input->Get("SimTree");
-	std::vector<Mask::Nucleus>* dataHandle = new std::vector<Mask::Nucleus>();
-	intree->SetBranchAddress("nuclei", &dataHandle);
-
-	output->cd();
-	TTree* outtree = new TTree("SimTree", "SimTree");
-	outtree->Branch("nuclei", dataHandle);
-
-	input->cd();
-
-	stats<<"Efficiency statistics for data from "<<inputname<<" using the ANASEN geometry"<<std::endl;
-	stats<<"Given in order of target=0, projectile=1, ejectile=2, residual=3, .... etc."<<std::endl;
-
-	intree->GetEntry(1);
-	std::vector<int> counts;
-	std::vector<int> coinc_counts;
-	counts.resize(dataHandle->size());
-	switch(counts.size())
-	{
-		case 3: coinc_counts.resize(1, 0); break;
-		case 4: coinc_counts.resize(1, 0); break;
-		case 6: coinc_counts.resize(4, 0); break;
-		case 8: coinc_counts.resize(11, 0); break;
-		default:
-		{
-			std::cerr<<"Bad reaction type at AnasenArray::CalculateEfficiency (given value: "<<counts.size()<<"). Quiting..."<<std::endl;
-			input->Close();
-			output->Close();
-			stats.close();
-			return;
-		}
-	}
-
-	uint64_t nentries = intree->GetEntries();
-	uint64_t percent5 = nentries*0.05;
-	uint64_t count = 0;
-	uint64_t npercent = 0;
-
-	Mask::Nucleus nucleus;
-	for(uint64_t i=0; i<nentries; i++)
-	{
-		intree->GetEntry(i);
-		if(++count == percent5)
-		{//Show progress every 5%
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		for(std::size_t j=0; j<dataHandle->size(); j++)
-		{
-			Mask::Nucleus& nucleus = (*dataHandle)[j];
-			DetectorResult result = IsAnasen(nucleus);
-			if(result.detectFlag)
-			{
-				nucleus.isDetected = true;
-				nucleus.detectedKE = result.energy_deposited;
-				nucleus.detectedTheta = result.direction.Theta();
-				nucleus.detectedPhi = result.direction.Phi();
-				counts[j]++;
-			}
-		}
-
-		CountCoincidences(*dataHandle, coinc_counts);
-
-		outtree->Fill();
-	}
-	input->Close();
-	output->cd();
-	outtree->Write(outtree->GetName(), TObject::kOverwrite);
-	output->Close();
-
-	delete dataHandle;
-
-	stats<<std::setw(10)<<"Index"<<"|"<<std::setw(10)<<"Efficiency"<<std::endl;
-	stats<<"---------------------"<<std::endl;
-	for(unsigned int i=0; i<counts.size(); i++)
-	{
-		stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	stats<<"Coincidence Efficiency"<<std::endl;
-	stats<<"---------------------"<<std::endl;
-	if(dataHandle->size() == 3)
-	{
-		stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	else if(dataHandle->size() == 4)
-	{
-		stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	else if(dataHandle->size() == 6)
-	{
-		stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	else if(dataHandle->size() == 8)
-	{
-		stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	stats.close();
-
-	std::cout<<std::endl;
-	std::cout<<"Complete."<<std::endl;
-	std::cout<<"---------------------------------------------"<<std::endl;
-}

src/Detectors/AnasenArray.h

@@ -6,8 +6,8 @@
 #include "DetectorArray.h"
 #include "SX3Detector.h"
 #include "QQQDetector.h"
-#include "Target.h"
-#include "Nucleus.h"
+#include "Mask/Target.h"
+#include "Mask/Nucleus.h"
 #include "AnasenDeadChannelMap.h"
 
 class AnasenArray : public DetectorArray
@@ -15,17 +15,15 @@ class AnasenArray : public DetectorArray
 public:
 	AnasenArray();
 	~AnasenArray();
-	virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
+	virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus);
 	virtual void DrawDetectorSystem(const std::string& filename) override;
 	virtual double RunConsistencyCheck() override;
-	inline void SetDeadChannelMap(const std::string& filename) { dmap.LoadMapfile(filename); }
+	virtual void SetDeadChannelMap(const std::string& filename) override { dmap.LoadMapfile(filename); }
 
 private:
-	DetectorResult IsRing1(Mask::Nucleus& nucleus);
-	DetectorResult IsRing2(Mask::Nucleus& nucleus);
-	DetectorResult IsQQQ(Mask::Nucleus& nucleus);
-	DetectorResult IsAnasen(Mask::Nucleus& nucleus);
-	void CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts);
+	DetectorResult IsRing1(const Mask::Nucleus& nucleus);
+	DetectorResult IsRing2(const Mask::Nucleus& nucleus);
+	DetectorResult IsQQQ(const Mask::Nucleus& nucleus);
 
 	std::vector<SX3Detector> m_Ring1;
 	std::vector<SX3Detector> m_Ring2;

src/Detectors/CMakeLists.txt

@@ -18,6 +18,9 @@ target_sources(Detectors PUBLIC
     SabreArray.h
     SX3Detector.cpp
     SX3Detector.h
+    DetectorApp.h
+    DetectorApp.cpp
+    DetectorArray.cpp
 )
 
 target_link_libraries(Detectors

src/Detectors/DetectorApp.cpp

@@ -0,0 +1,144 @@
+#include "DetectorApp.h"
+#include <fstream>
+#include <iostream>
+
+DetectorApp::DetectorApp() :
+    m_resources(nullptr)
+{
+}
+
+DetectorApp::~DetectorApp()
+{
+    for(std::size_t i=0; i<m_detectorList.size(); i++)
+        delete m_detectorList[i];
+}
+
+bool DetectorApp::LoadConfig(const std::string& filename)
+{
+    std::cout<<"----------Detector Efficiency Calculation----------"<<std::endl;
+    std::ifstream input(filename);
+    if(!input.is_open())
+    {
+        std::cerr<<"Unable to open input config file "<<filename<<std::endl;
+        return false;
+    }
+
+    std::string junk;
+    std::string typeName;
+
+    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())
+    {
+        std::cerr << "Unable to open input data file " << m_inputFileName << std::endl;
+        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;
+    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;
+    }
+    std::cout << "Done" << std::endl;
+
+    std::cout << "Ready to launch!" << std::endl;
+
+    return true;
+}
+
+void DetectorApp::Run()
+{
+	std::cout<<"Running efficiency calculation..."<<std::endl;
+
+    if(m_detectorList.size() != m_nthreads)
+    {
+        std::cerr << "Detector list not equal to number of threads" << std::endl;
+        return;
+    }
+
+    for(uint64_t i=0; i<m_detectorList.size(); i++)
+    {
+        //Create a job for the thread pool, using a lambda and providing a tuple of the arguments
+        m_resources->PushJob({[this](DetectorArray* array, uint64_t chunkSamples)
+            {
+                if(system == nullptr)
+		    	    return;
+
+                std::vector<Mask::Nucleus> data;
+                DetectorResult result;
+	            for(uint64_t i=0; i<chunkSamples; i++)
+	            {
+                    m_fileReader.Read(data);
+                    for(auto& nucleus : data)
+                    {
+                        result = array->IsDetected(nucleus);
+                        if(result.detectFlag)
+                        {
+                            nucleus.isDetected = true;
+                            nucleus.detectedKE = result.energy_deposited;
+                            nucleus.detectedTheta = result.direction.Theta();
+                            nucleus.detectedPhi = result.direction.Phi();
+                        }
+                    }
+		            m_fileWriter.PushData(data);
+	            }
+            },
+        {m_detectorList[i], m_chunkSamples[i]} //arguments to function, in order
+        }
+        );
+    }
+
+	uint64_t size = m_fileReader.GetSize();
+	uint64_t count = 0;
+	double percent = 0.05;
+	uint64_t flushVal = size*percent;
+	uint64_t flushCount = 0;
+
+	while(true)
+	{
+        if(count == flushVal)
+	    {
+	    	count = 0;
+	    	++flushCount;
+	    	std::cout<<"\rPercent of data written to disk: "<<percent*flushCount*100<<"%"<<std::flush;
+	    }
+
+		if(m_resources->IsFinished() && m_fileWriter.GetQueueSize() == 0)
+				break;
+		else if(m_fileWriter.Write())
+				++count;
+	}
+	
+}

src/Detectors/DetectorApp.h

@@ -0,0 +1,40 @@
+#ifndef DETECTOR_APP_H
+#define DETECTOR_APP_H
+
+#include "DetectorArray.h"
+#include "Mask/FileWriter.h"
+#include "Mask/FileReader.h"
+#include "Mask/ThreadPool.h"
+
+#include <string>
+#include <vector>
+#include <memory>
+
+class DetectorApp
+{
+public:
+    DetectorApp();
+    ~DetectorApp();
+
+    bool LoadConfig(const std::string& filename);
+    
+    void Run();
+
+private:
+
+    std::vector<DetectorArray*> m_detectorList; //One array per thread
+    std::vector<uint64_t> m_chunkSamples;
+    Mask::FileWriter m_fileWriter;
+    Mask::FileReader m_fileReader;
+
+    std::string m_inputFileName;
+    std::string m_outputFileName;
+    std::string m_deadChannelFileName;
+
+    uint64_t m_nthreads;
+    uint64_t m_nentries;
+
+    std::unique_ptr<Mask::ThreadPool<DetectorArray*, uint64_t>> m_resources;
+};
+
+#endif

src/Detectors/DetectorArray.cpp

@@ -0,0 +1,35 @@
+#include "DetectorArray.h"
+#include "AnasenArray.h"
+#include "SabreArray.h"
+
+DetectorArray* CreateDetectorArray(ArrayType type)
+{
+    switch(type)
+    {
+        case ArrayType::None: return nullptr;
+        case ArrayType::Anasen: return new AnasenArray();
+        case ArrayType::Sabre: return new SabreArray();
+    }
+    return nullptr;
+}
+
+std::string ArrayTypeToString(ArrayType type)
+{
+    switch(type)
+    {
+        case ArrayType::None: return "None";
+        case ArrayType::Anasen: return "Anasen";
+        case ArrayType::Sabre: return "Sabre";
+    }
+    return "None";
+}
+
+ArrayType StringToArrayType(const std::string& value)
+{
+    if (value == "Anasen")
+        return ArrayType::Anasen;
+    else if (value == "Sabre")
+        return ArrayType::Sabre;
+    else
+        return ArrayType::None;
+}

src/Detectors/DetectorArray.h

@@ -5,6 +5,7 @@
 #include <cmath>
 
 #include "Math/Point3D.h"
+#include "Mask/Nucleus.h"
 
 struct DetectorResult
 {
@@ -14,15 +15,23 @@ struct DetectorResult
 	std::string det_name = "";
 };
 
+enum class ArrayType
+{
+	None,
+	Anasen,
+	Sabre
+};
+
 class DetectorArray
 {
 public:
 	DetectorArray() {};
 	virtual ~DetectorArray() {};
 
-	virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) = 0;
+	virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus) = 0;
 	virtual void DrawDetectorSystem(const std::string& filename) = 0;
 	virtual double RunConsistencyCheck() = 0;
+	virtual void SetDeadChannelMap(const std::string& filename) = 0;
 
 protected:
 	bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < s_epsilon ? true : false; };
@@ -30,4 +39,10 @@ protected:
 	static constexpr double s_epsilon = 1.0e-6;
 };
 
+DetectorArray* CreateDetectorArray(ArrayType type);
+
+std::string ArrayTypeToString(ArrayType type);
+
+ArrayType StringToArrayType(const std::string& value);
+
 #endif

src/Detectors/SabreArray.cpp

@@ -31,161 +31,6 @@ SabreArray::SabreArray() :
 
 SabreArray::~SabreArray() {}
 
-void SabreArray::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname)
-{
-	std::cout<<"----------SABRE Efficiency Calculation----------"<<std::endl;
-	std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
-	std::cout<<"Running efficiency calculation..."<<std::endl;
-
-	if(!m_deadMap.IsValid())
-	{
-		std::cerr<<"Unable to run SABRE Efficiency without a dead channel map."<<std::endl;
-		std::cerr<<"If you have no dead channels, simply make a file that's empty"<<std::endl;
-		std::cerr<<"Exiting."<<std::endl;
-		std::cout<<"---------------------------------------------"<<std::endl;
-	}
-
-
-	TFile* input = TFile::Open(inputname.c_str(), "READ");
-	TFile* output = TFile::Open(outputname.c_str(), "RECREATE");
-	std::ofstream stats(statsname);
-	stats<<std::setprecision(5);
-
-	TTree* intree = (TTree*) input->Get("SimTree");
-	std::vector<Mask::Nucleus>* dataHandle = new std::vector<Mask::Nucleus>();
-	intree->SetBranchAddress("nuclei", &dataHandle);
-
-	output->cd();
-	TTree* outtree = new TTree("SimTree", "SimTree");
-	outtree->Branch("nuclei", dataHandle);
-
-	input->cd();
-
-	stats<<"Efficiency statistics for data from "<<inputname<<" using the ANASEN geometry"<<std::endl;
-	stats<<"Given in order of target=0, projectile=1, ejectile=2, residual=3, .... etc."<<std::endl;
-
-	intree->GetEntry(1);
-	std::vector<int> counts;
-	std::vector<int> coinc_counts;
-	counts.resize(dataHandle->size());
-	switch(counts.size())
-	{
-		case 3: coinc_counts.resize(1, 0); break;
-		case 4: coinc_counts.resize(1, 0); break;
-		case 6: coinc_counts.resize(4, 0); break;
-		case 8: coinc_counts.resize(11, 0); break;
-		default:
-		{
-			std::cerr<<"Bad reaction type at AnasenEfficiency::CalculateEfficiency (given value: "<<counts.size()<<"). Quiting..."<<std::endl;
-			input->Close();
-			output->Close();
-			stats.close();
-			return;
-		}
-	}
-
-	uint64_t nentries = intree->GetEntries();
-	uint64_t percent5 = nentries*0.05;
-	uint64_t count = 0;
-	uint64_t npercent = 0;
-
-	Mask::Nucleus nucleus;
-	for(uint64_t i=0; i<nentries; i++)
-	{
-		intree->GetEntry(i);
-		if(++count == percent5)
-		{//Show progress every 5%
-			npercent++;
-			count = 0;
-			std::cout<<"\rPercent completed: "<<npercent*5<<"%"<<std::flush;
-		}
-
-		for(std::size_t j=0; j<dataHandle->size(); j++)
-		{
-			Mask::Nucleus& nucleus = (*dataHandle)[j];
-			DetectorResult result = IsSabre(nucleus);
-			if(result.detectFlag)
-			{
-				nucleus.isDetected = true;
-				nucleus.detectedKE = result.energy_deposited;
-				nucleus.detectedTheta = result.direction.Theta();
-				nucleus.detectedPhi = result.direction.Phi();
-				counts[j]++;
-			}
-		}
-
-		CountCoincidences(*dataHandle, coinc_counts);
-
-		outtree->Fill();
-	}
-	input->Close();
-	output->cd();
-	outtree->Write(outtree->GetName(), TObject::kOverwrite);
-	output->Close();
-
-	delete dataHandle;
-
-	stats<<std::setw(10)<<"Index"<<"|"<<std::setw(10)<<"Efficiency"<<std::endl;
-	stats<<"---------------------"<<std::endl;
-	for(unsigned int i=0; i<counts.size(); i++) {
-		stats<<std::setw(10)<<i<<"|"<<std::setw(10)<<((double)counts[i]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	stats<<"Coincidence Efficiency"<<std::endl;
-	stats<<"---------------------"<<std::endl;
-	if(counts.size() == 3)
-	{
-		stats<<std::setw(10)<<"1 + 2"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	else if(counts.size() == 4)
-	{
-		stats<<std::setw(10)<<"2 + 3"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	else if(counts.size() == 6)
-	{
-		stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 4 + 5"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	else if(counts.size() == 8)
-	{
-		stats<<std::setw(10)<<"2 + 4"<<"|"<<std::setw(10)<<((double)coinc_counts[0]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[1]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[2]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[3]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[4]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[5]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 4 + 6"<<"|"<<std::setw(10)<<((double)coinc_counts[6]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 4 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[7]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[8]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[9]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-		stats<<std::setw(10)<<"2 + 4 + 6 + 7"<<"|"<<std::setw(10)<<((double)coinc_counts[10]/nentries)<<std::endl;
-		stats<<"---------------------"<<std::endl;
-	}
-	stats.close();
-
-	std::cout<<std::endl;
-	std::cout<<"Complete."<<std::endl;
-	std::cout<<"---------------------------------------------"<<std::endl;
-}
-
 void SabreArray::DrawDetectorSystem(const std::string& filename)
 {
 	std::ofstream output(filename);
@@ -259,7 +104,7 @@ double SabreArray::RunConsistencyCheck()
 }
 
 /*Returns if detected, as well as total energy deposited in SABRE*/
-DetectorResult SabreArray::IsSabre(Mask::Nucleus& nucleus)
+DetectorResult SabreArray::IsDetected(const Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	if(nucleus.GetKE() <= s_energyThreshold)
@@ -298,81 +143,3 @@ DetectorResult SabreArray::IsSabre(Mask::Nucleus& nucleus)
 	observation.detectFlag = false;
 	return observation;
 }
-
-void SabreArray::CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts)
-{
-	if (data.size() == 3 && data[1].isDetected && data[2].isDetected)
-	{
-		counts[0]++;
-	}
-	else if (data.size() == 4 && data[2].isDetected && data[3].isDetected)
-	{
-		counts[0]++;
-	}
-	else if(data.size() == 6)
-	{
-		if(data[2].isDetected && data[4].isDetected) 
-		{
-			counts[0]++;
-		}
-		if(data[2].isDetected && data[5].isDetected)
-		{
-			counts[1]++;
-		}
-		if(data[4].isDetected && data[5].isDetected)
-		{
-			counts[2]++;
-		}
-		if(data[2].isDetected && data[4].isDetected && data[5].isDetected)
-		{
-			counts[3]++;
-		}
-	}
-	else if(data.size() == 8)
-	{
-		if(data[2].isDetected && data[4].isDetected) 
-		{
-			counts[0]++;
-		}
-		if(data[2].isDetected && data[6].isDetected)
-		{
-			counts[1]++;
-		}
-		if(data[2].isDetected && data[7].isDetected)
-		{
-			counts[2]++;
-		}
-		if(data[4].isDetected && data[6].isDetected)
-		{
-			counts[3]++;
-		}
-		if(data[4].isDetected && data[7].isDetected)
-		{
-			counts[4]++;
-		}
-		if(data[6].isDetected && data[7].isDetected)
-		{
-			counts[5]++;
-		}
-		if(data[2].isDetected && data[4].isDetected && data[6].isDetected)
-		{
-			counts[6]++;
-		}
-		if(data[2].isDetected && data[4].isDetected && data[7].isDetected)
-		{
-			counts[7]++;
-		}
-		if(data[2].isDetected && data[6].isDetected && data[7].isDetected)
-		{
-			counts[8]++;
-		}
-		if(data[4].isDetected && data[6].isDetected && data[7].isDetected)
-		{
-			counts[9]++;
-		}
-		if(data[2].isDetected && data[4].isDetected && data[6].isDetected && data[7].isDetected)
-		{
-			counts[10]++;
-		}
-	}
-}

src/Detectors/SabreArray.h

@@ -3,7 +3,7 @@
 
 #include "DetectorArray.h"
 #include "SabreDetector.h"
-#include "Target.h"
+#include "Mask/Target.h"
 #include "SabreDeadChannelMap.h"
 #include "Mask/Nucleus.h"
 
@@ -12,14 +12,12 @@ class SabreArray : public DetectorArray
 public:
 	SabreArray();
 	~SabreArray();
-    void SetDeadChannelMap(const std::string& filename) { m_deadMap.LoadMapfile(filename); };
-	void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
+    virtual void SetDeadChannelMap(const std::string& filename) override { m_deadMap.LoadMapfile(filename); };
+	virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus) override;
     void DrawDetectorSystem(const std::string& filename) override;
     double RunConsistencyCheck() override;
 
 private:
-	DetectorResult IsSabre(Mask::Nucleus& nucleus);
-    void CountCoincidences(const std::vector<Mask::Nucleus>& data, std::vector<int>& counts);
 
 	std::vector<SabreDetector> m_detectors;
     

src/Detectors/main.cpp

@@ -1,5 +1,4 @@
-#include "SabreArray.h"
-#include "AnasenArray.h"
+#include "DetectorApp.h"
 #include "KinematicsExceptions.h"
 #include <iostream>
 #include <string>
@@ -7,7 +6,7 @@
 int main(int argc, char** argv)
 {
 
-	if(argc != 4)
+	if(argc != 2)
 	{
 		std::cerr<<"Incorrect number of commandline arguments! Returning."<<std::endl;
 		return 1;
@@ -19,26 +18,21 @@ int main(int argc, char** argv)
 		return 1;
 	}
 
-	std::string inputname = argv[1];
-	std::string outputname = argv[2];
-	std::string statsname = argv[3];
+	try
+	{
+		DetectorApp app;
+		if(!app.LoadConfig(argv[1]))
+		{
+			std::cerr << "Unable to load config file " << argv[1] << ". Shutting down." << std::endl;
+			return 1;
+		}
 
-	
-	SabreArray sabre;
-	std::string mapfile = "./etc/sabreDeadChannels_May2022.txt";
-	sabre.SetDeadChannelMap(mapfile);
-	sabre.CalculateEfficiency(inputname, outputname, statsname);
-	//std::cout<<"Running consistency check(0=success): "<<sabre.RunConsistencyCheck()<<std::endl;
-	//sabre.DrawDetectorSystem("/data1/gwm17/10B3He/Feb2021/simulation/SABREGeo.txt");
-	
-
-	/*
- 	AnasenArray anasen;
- 	std::string mapfile = "./etc/AnasenDeadChannels.txt";
- 	anasen.SetDeadChannelMap(mapfile);
- 	anasen.CalculateEfficiency(inputname, outputname, statsname);
- 	//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
- 	//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_fixedZ.txt");
-	*/
+		app.Run();
+	}
+	catch(std::exception& e)
+	{
+		std::cerr << e.what() << std::endl;
+		return 1;
+	}
 	return 0;
 }

src/Kinematics/main.cpp

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

src/Mask/CMakeLists.txt

@@ -50,7 +50,15 @@ target_sources(Mask PRIVATE
     ThreeStepSystem.h
     TwoStepSystem.cpp
     TwoStepSystem.h
+    ThreadPool.h
+    FileWriter.h
+    FileWriter.cpp
+    FileReader.h
+    FileReader.cpp
 )
 
-target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES})
+set(THREADS_PREFER_PTHREAD_FLAG On)
+find_package(Threads REQUIRED)
+
+target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES} Threads::Threads)
 set_target_properties(Mask PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})

src/Mask/FileReader.cpp

@@ -0,0 +1,66 @@
+#include "FileReader.h"
+
+namespace Mask {
+
+    FileReader::FileReader() :
+        m_file(nullptr), m_tree(nullptr), m_branchHandle(nullptr), m_currentEntry(0), m_size(0)
+    {
+    }
+
+    FileReader::FileReader(const std::string& filename, const std::string& treename) :
+        m_file(nullptr), m_tree(nullptr), m_branchHandle(nullptr), m_currentEntry(0), m_size(0)
+    {
+        Open(filename, treename);
+    }
+
+    FileReader::~FileReader()
+    {
+        Close();
+    }
+
+    void FileReader::Open(const std::string& filename, const std::string& treename)
+    {
+        if(m_file != nullptr || m_tree != nullptr)
+            Close();
+
+        m_file = TFile::Open(filename.c_str(), "READ");
+        if(m_file != nullptr && m_file->IsOpen())
+        {
+            m_tree = (TTree*) m_file->Get(treename.c_str());
+            if(m_tree == nullptr)
+            {
+                m_file->Close();
+                delete m_file;
+                m_file = nullptr;
+            }
+            m_branchHandle = new std::vector<Nucleus>();
+            m_tree->SetBranchAddress("nuclei", &m_branchHandle);
+            m_size = m_tree->GetEntries();
+            m_currentEntry = 0; //Reset file position
+        }
+    }
+
+    void FileReader::Close()
+    {
+        if(m_file != nullptr && m_file->IsOpen())
+        {
+            m_file->Close();
+            delete m_file;
+            m_file = nullptr;
+        }
+    }
+
+    bool FileReader::Read(std::vector<Nucleus>& dataHandle)
+    {
+        std::scoped_lock<std::mutex> guard(m_fileMutex);
+        int bytes = m_tree->GetEntry(m_currentEntry);
+        if(bytes != 0)
+        {
+            dataHandle = *m_branchHandle;
+            m_currentEntry++;
+            return true;
+        }
+
+        return false;
+    }
+}

src/Mask/FileReader.h

@@ -0,0 +1,46 @@
+#ifndef FILE_READER_H
+#define FILE_READER_H
+
+#include "Nucleus.h"
+
+#include "TFile.h"
+#include "TTree.h"
+
+#include <vector>
+#include <string>
+#include <mutex>
+#include <atomic>
+
+namespace Mask {
+
+    class FileReader
+    {
+    public:
+        FileReader();
+        FileReader(const std::string& filename, const std::string& treename);
+        ~FileReader();
+
+        void Open(const std::string& filename, const std::string& treename); //Not thread safe
+        void Close(); //Not thread safe
+
+        /*
+            Read: fills entry to given dataHandle. Returns true if data was successfully filled, otherwise returns false
+        */
+        bool Read(std::vector<Nucleus>& dataHandle); //Thread safe
+        uint64_t GetSize() { return m_size; }//In entries (implicitly thread safe)
+        bool IsOpen() { return m_file == nullptr ? false : m_file->IsOpen(); } //Should be safe?
+        bool IsTree() { return m_tree != nullptr; } //Should be safe?
+
+    private:
+        TFile* m_file;
+        TTree* m_tree;
+
+        std::vector<Nucleus>* m_branchHandle;
+
+        std::mutex m_fileMutex;
+        std::atomic<uint64_t> m_currentEntry;
+        std::atomic<uint64_t> m_size; //in entries
+    };
+}
+
+#endif

src/Mask/FileWriter.cpp

@@ -0,0 +1,75 @@
+#include "FileWriter.h"
+
+namespace Mask {
+
+    FileWriter::FileWriter() :
+        m_file(nullptr), m_tree(nullptr), m_queueSize(0)
+    {
+    }
+
+    FileWriter::FileWriter(const std::string& filename, const std::string& treename) :
+        m_file(nullptr), m_tree(nullptr)
+    {
+        m_file = TFile::Open(filename.c_str(), "RECREATE");
+        if(m_file != nullptr && m_file->IsOpen())
+        {
+            m_tree = new TTree(treename.c_str(), treename.c_str());
+            m_tree->Branch("nuclei", &m_dataHandle);
+        }
+    }
+
+    FileWriter::~FileWriter()
+    {
+        Close();
+    }
+
+    void FileWriter::Open(const std::string& filename, const std::string& treename)
+    {
+        if(m_file != nullptr || m_tree != nullptr)
+            Close();
+
+        m_file = TFile::Open(filename.c_str(), "RECREATE");
+        if(m_file != nullptr && m_file->IsOpen())
+        {
+            m_tree = new TTree(treename.c_str(), treename.c_str());
+            m_tree->Branch("nuclei", &m_dataHandle);
+        }
+    }
+
+    void FileWriter::Close()
+    {
+        if(m_file->IsOpen())
+        {
+            if(m_tree != nullptr)
+                m_tree->Write(m_tree->GetName(), TObject::kOverwrite);
+
+            m_file->Close();
+            delete m_file;
+            m_file = nullptr;
+        }
+    }
+
+    void FileWriter::PushData(const std::vector<Nucleus>& data)
+    {
+        std::scoped_lock<std::mutex> guard(m_queueMutex);
+        m_queue.push(data);
+        ++m_queueSize;
+    }
+
+    bool FileWriter::Write()
+    {
+        if(m_queueSize == 0)
+            return false;
+
+        //Aquire lock for as short a time as possible
+        {
+            std::scoped_lock<std::mutex> guard(m_queueMutex);
+            m_dataHandle = m_queue.front();
+            m_queue.pop();
+        }
+        --m_queueSize;
+
+        m_tree->Fill();
+        return true;
+    }
+}

src/Mask/FileWriter.h

@@ -0,0 +1,45 @@
+#ifndef FILE_WRITER_H
+#define FILE_WRITER_H
+
+#include "Nucleus.h"
+
+#include "TFile.h"
+#include "TTree.h"
+
+#include <queue>
+#include <mutex>
+#include <atomic>
+
+namespace Mask {
+
+    class FileWriter
+    {
+    public:
+        FileWriter();
+        FileWriter(const std::string& filename, const std::string& treename);
+        ~FileWriter();
+
+        bool IsOpen() const { return m_file == nullptr ? false : m_file->IsOpen(); }
+        bool IsTree() const { return m_tree == nullptr ? false : true; }
+
+        std::size_t GetQueueSize() const { return m_queueSize; } //Implicitly thread-safe
+
+        void PushData(const std::vector<Nucleus>& data); //Thread-safe
+        bool Write(); //Not completely thread-safe, should only be used by main application loop (acess to queue is safe, but all else not)
+
+        void Open(const std::string& filename, const std::string& treename); //Not thread safe!
+        void Close(); //Not thread safe!
+
+    private:
+        TFile* m_file;
+        TTree* m_tree;
+
+        std::vector<Nucleus> m_dataHandle;
+
+        std::mutex m_queueMutex;
+        std::atomic<std::size_t> m_queueSize;
+        std::queue<std::vector<Nucleus>> m_queue;
+    };
+}
+
+#endif

src/Mask/MaskApp.cpp

@@ -8,7 +8,7 @@
 namespace Mask {
 
 	MaskApp::MaskApp() :
-		m_system(nullptr)
+		m_system(nullptr), m_resources(nullptr)
 	{
 		std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
 	}
@@ -16,6 +16,8 @@ namespace Mask {
 	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) 
@@ -31,10 +33,25 @@ namespace Mask {
 		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)
@@ -98,10 +115,20 @@ namespace Mask {
 		m_system = CreateSystem(params);
 		if(m_system == nullptr || !m_system->IsValid())
 		{
-			std::cerr<<"Param size: "<<params.size()<<std::endl;
 			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;
+			}
+		}
+
 		std::getline(input, junk);
 		std::getline(input, junk);
 
@@ -132,9 +159,13 @@ namespace Mask {
 		}
 		m_system->SetLayeredTarget(target);
 
-		std::cout<<"Outputing data to file: "<<m_outputName<<std::endl;
-		std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
-		std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
+		for(auto system : m_systemList)
+		{
+			system->SetLayeredTarget(target);
+		}
+
+		std::cout<<"Reaction equation: "<<m_system->GetSystemEquation()<<std::endl;
+		std::cout<<"Number of samples: "<<m_nsamples<<std::endl;
 	
 		return true;
 	}
@@ -146,6 +177,57 @@ namespace Mask {
 	}
 
 	void MaskApp::Run()
+	{
+		std::cout<<"Running simulation..."<<std::endl;
+		if(m_systemList.size() != m_nthreads)
+		{
+			std::cerr << "System list not equal to number of threads" << std::endl;
+			return;
+		}
+
+		//Give our thread pool some tasks
+		for(std::size_t i=0; i<m_systemList.size(); i++)
+		{
+			//bind a lambda to the job, taking in a ReactionSystem, and then provide a reaction system as the tuple arguments.
+			m_resources->PushJob({[this](ReactionSystem* system, uint64_t chunkSamples) 
+				{
+					if(system == nullptr)
+						return;
+
+					for(uint64_t i=0; i<chunkSamples; i++)
+					{
+						system->RunSystem();
+						m_fileWriter.PushData(*(system->GetNuclei()));
+					}
+				}, 
+			{m_systemList[i], m_chunkSamples[i]}});
+		}
+
+		uint64_t count = 0;
+		double percent = 0.05;
+		uint64_t flushVal = m_nsamples*percent;
+		uint64_t flushCount = 0;
+		while(true)
+		{
+			if(count == flushVal)
+			{
+				count = 0;
+				++flushCount;
+				std::cout<<"\rPercent of data written to disk: "<<percent*flushCount*100<<"%"<<std::flush;
+			}
+
+			if(m_resources->IsFinished() && m_fileWriter.GetQueueSize() == 0)
+				break;
+			else if(m_fileWriter.Write())
+				++count;
+		}
+
+		std::cout<<std::endl;
+		std::cout<<"Complete."<<std::endl;
+		std::cout<<"---------------------------------------------"<<std::endl;
+	}
+	
+	void MaskApp::RunSingleThread()
 	{
 		std::cout<<"Running simulation..."<<std::endl;
 		if(m_system == nullptr) 
@@ -158,11 +240,11 @@ namespace Mask {
 		tree->Branch("nuclei", m_system->GetNuclei());
 	
 		//For progress tracking
-		uint32_t percent5 = 0.05*m_nsamples;
-		uint32_t count = 0;
-		uint32_t npercent = 0;
+		uint64_t percent5 = 0.05*m_nsamples;
+		uint64_t count = 0;
+		uint64_t npercent = 0;
 	
-		for(uint32_t i=0; i<m_nsamples; i++) 
+		for(uint64_t i=0; i<m_nsamples; i++) 
 		{
 			if(++count == percent5) 
 			{

src/Mask/MaskApp.h

@@ -7,6 +7,10 @@
 #include "TwoStepSystem.h"
 #include "ThreeStepSystem.h"
 #include "RxnType.h"
+#include "ThreadPool.h"
+#include "FileWriter.h"
+
+#include <memory>
 
 namespace Mask {
 
@@ -17,18 +21,24 @@ namespace Mask {
 		~MaskApp();
 		bool LoadConfig(const std::string& filename);
 		bool SaveConfig(const std::string& filename);
-		int GetNumberOfSamples() const { return m_nsamples; }
-		const std::string GetSystemName() const { return m_system == nullptr ? "Invalid System" : m_system->GetSystemEquation(); }
-		const std::string GetOutputName() const { return m_outputName; }
-		const RxnType GetReactionType() const { return m_rxnType; }
 
 		void Run();
+		void RunSingleThread();
 
 	private:
+		void RunChunk(ReactionSystem* system);
+
 		ReactionSystem* m_system;
+
 		std::string m_outputName;
 		RxnType m_rxnType;
-		int m_nsamples;
+		uint64_t m_nsamples;
+		uint32_t m_nthreads;
+
+		std::vector<ReactionSystem*> m_systemList; //One system for each thread
+		std::vector<uint64_t> m_chunkSamples;
+		FileWriter m_fileWriter;
+		std::unique_ptr<ThreadPool<ReactionSystem*, uint64_t>> m_resources;
 	
 	};
 

src/Mask/RandomGenerator.cpp

@@ -1,7 +1,6 @@
 #include "RandomGenerator.h"
 
 namespace Mask {
-	RandomGenerator* RandomGenerator::s_instance = new RandomGenerator();
 	
 	RandomGenerator::RandomGenerator()
 	{

src/Mask/RandomGenerator.h

@@ -10,12 +10,15 @@ namespace Mask {
 	public:
 		~RandomGenerator();
 		std::mt19937& GetGenerator() { return rng; }
-		static RandomGenerator& GetInstance() { return *s_instance; }
+		static RandomGenerator& GetInstance()
+		{
+			static thread_local RandomGenerator s_instance;
+			return s_instance;
+		}
 
 	private:
 		RandomGenerator();
 		
-		static RandomGenerator* s_instance;
 		std::mt19937 rng;
 	};
 

src/Mask/ThreadPool.h

@@ -0,0 +1,118 @@
+/*
+    ThreadPool.h
+    Simple thread pool implementation. Jobs are q'd, and as threads are made available, jobs are dispatched to these threads.
+
+    GWM -- Aug 2022
+*/
+#ifndef THREAD_POOL_H
+#define THREAD_POOL_H
+
+#include "ReactionSystem.h"
+
+#include <vector>
+#include <atomic>
+#include <thread>
+#include <mutex>
+#include <functional>
+#include <queue>
+#include <condition_variable>
+
+
+namespace Mask {
+
+	template<typename... Types>
+	class ThreadPool
+	{
+	public:
+
+		using JobFunction = std::function<void(Types...)>;
+		struct Job
+		{
+			JobFunction func;
+			std::tuple<Types...> argument;
+		};
+
+
+		ThreadPool(uint32_t nthreads) :
+			m_isStopped(false), m_numberRunning(0), m_queueSize(0), m_initShutdown(false)
+		{
+			for(uint32_t i=0; i<nthreads; i++)
+			{
+				m_pool.push_back(std::thread(std::bind(&ThreadPool::ExecuteJob, std::ref(*this))));
+			}
+		}
+		~ThreadPool()
+		{
+			if(!m_isStopped)
+				Shutdown();
+		}
+
+		void PushJob(Job job)
+		{
+			{
+				std::unique_lock<std::mutex> guard(m_poolMutex);
+				m_queue.push(job);
+				m_queueSize++;
+			}
+			m_wakeCondition.notify_one();
+		}
+
+		bool IsFinished()
+		{
+			return m_numberRunning == 0 && m_queueSize == 0;
+		}
+
+		void Shutdown()
+		{
+			m_initShutdown = true;
+			m_wakeCondition.notify_all();
+			for(auto& thread : m_pool)
+			{
+				thread.join();
+			}
+
+			m_pool.clear();
+			m_isStopped = true;
+		}
+
+	private:
+		void ExecuteJob()
+		{
+			while(true)
+			{
+				Job job;
+				{
+					std::unique_lock<std::mutex> guard(m_poolMutex);
+					m_wakeCondition.wait(guard, [this](){ 
+						return (!m_queue.empty() || m_initShutdown); 
+					});
+					if(m_initShutdown)
+						return;
+	
+					job = m_queue.front();
+					m_queue.pop();
+				}
+
+				//Change number running first so that no crash 
+				m_numberRunning++;
+				m_queueSize--;
+
+				std::apply(job.func, job.argument);
+				m_numberRunning--;
+			}
+			
+		}
+
+		std::queue<Job, std::deque<Job>> m_queue;
+		std::vector<std::thread> m_pool;
+		std::mutex m_poolMutex;
+		std::condition_variable m_wakeCondition;
+
+		bool m_isStopped;
+		std::atomic<int> m_numberRunning;
+		std::atomic<int> m_queueSize;
+		std::atomic<bool> m_initShutdown;
+	};
+}
+
+#endif