README.md

@@ -34,13 +34,11 @@ 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 <input_file>`
+`./bin/Detectors <kinematics_datafile> <new_detection_datafile> <new_detection_statsfile>`
 
-An example input file is provided with the repository.
+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.
 
 ## 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

@@ -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

input.txt

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

src/Detectors/AnasenArray.cpp

@@ -230,7 +230,7 @@ double AnasenArray::RunConsistencyCheck()
 
 }
 
-DetectorResult AnasenArray::IsRing1(const Mask::Nucleus& nucleus)
+DetectorResult AnasenArray::IsRing1(Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	double thetaIncident;
@@ -255,7 +255,7 @@ DetectorResult AnasenArray::IsRing1(const Mask::Nucleus& nucleus)
 	return observation;
 }
 
-DetectorResult AnasenArray::IsRing2(const Mask::Nucleus& nucleus)
+DetectorResult AnasenArray::IsRing2(Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	double thetaIncident;
@@ -280,7 +280,7 @@ DetectorResult AnasenArray::IsRing2(const Mask::Nucleus& nucleus)
 	return observation;
 }
 
-DetectorResult AnasenArray::IsQQQ(const Mask::Nucleus& nucleus)
+DetectorResult AnasenArray::IsQQQ(Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	double thetaIncident;
@@ -324,7 +324,7 @@ DetectorResult AnasenArray::IsQQQ(const Mask::Nucleus& nucleus)
 	return observation;
 }
 
-DetectorResult AnasenArray::IsDetected(const Mask::Nucleus& nucleus)
+DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus)
 {
 	DetectorResult result;
 	if(nucleus.GetKE() <= s_energyThreshold)
@@ -338,3 +338,235 @@ DetectorResult AnasenArray::IsDetected(const 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 "Mask/Target.h"
+#include "Target.h"
-#include "Mask/Nucleus.h"
+#include "Nucleus.h"
 #include "AnasenDeadChannelMap.h"
 
 class AnasenArray : public DetectorArray
@@ -15,15 +15,17 @@ class AnasenArray : public DetectorArray
 public:
 	AnasenArray();
 	~AnasenArray();
-	virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus);
+	virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
 	virtual void DrawDetectorSystem(const std::string& filename) override;
 	virtual double RunConsistencyCheck() override;
-	virtual void SetDeadChannelMap(const std::string& filename) override { dmap.LoadMapfile(filename); }
+	inline void SetDeadChannelMap(const std::string& filename) { dmap.LoadMapfile(filename); }
 
 private:
-	DetectorResult IsRing1(const Mask::Nucleus& nucleus);
+	DetectorResult IsRing1(Mask::Nucleus& nucleus);
-	DetectorResult IsRing2(const Mask::Nucleus& nucleus);
+	DetectorResult IsRing2(Mask::Nucleus& nucleus);
-	DetectorResult IsQQQ(const 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);
 
 	std::vector<SX3Detector> m_Ring1;
 	std::vector<SX3Detector> m_Ring2;

src/Detectors/CMakeLists.txt

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

src/Detectors/DetectorApp.cpp

@@ -1,144 +0,0 @@
-#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

@@ -1,40 +0,0 @@
-#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

@@ -1,35 +0,0 @@
-#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,7 +5,6 @@
 #include <cmath>
 
 #include "Math/Point3D.h"
-#include "Mask/Nucleus.h"
 
 struct DetectorResult
 {
@@ -15,23 +14,15 @@ struct DetectorResult
 	std::string det_name = "";
 };
 
-enum class ArrayType
-{
-	None,
-	Anasen,
-	Sabre
-};
-
 class DetectorArray
 {
 public:
 	DetectorArray() {};
 	virtual ~DetectorArray() {};
 
-	virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus) = 0;
+	virtual void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) = 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; };
@@ -39,10 +30,4 @@ 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,6 +31,161 @@ 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);
@@ -104,7 +259,7 @@ double SabreArray::RunConsistencyCheck()
 }
 
 /*Returns if detected, as well as total energy deposited in SABRE*/
-DetectorResult SabreArray::IsDetected(const Mask::Nucleus& nucleus)
+DetectorResult SabreArray::IsSabre(Mask::Nucleus& nucleus)
 {
 	DetectorResult observation;
 	if(nucleus.GetKE() <= s_energyThreshold)
@@ -143,3 +298,81 @@ DetectorResult SabreArray::IsDetected(const 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 "Mask/Target.h"
+#include "Target.h"
 #include "SabreDeadChannelMap.h"
 #include "Mask/Nucleus.h"
 
@@ -12,12 +12,14 @@ class SabreArray : public DetectorArray
 public:
 	SabreArray();
 	~SabreArray();
-    virtual void SetDeadChannelMap(const std::string& filename) override { m_deadMap.LoadMapfile(filename); };
+    void SetDeadChannelMap(const std::string& filename) { m_deadMap.LoadMapfile(filename); };
-	virtual DetectorResult IsDetected(const Mask::Nucleus& nucleus) override;
+	void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) 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,4 +1,5 @@
-#include "DetectorApp.h"
+#include "SabreArray.h"
+#include "AnasenArray.h"
 #include "KinematicsExceptions.h"
 #include <iostream>
 #include <string>
@@ -6,7 +7,7 @@
 int main(int argc, char** argv)
 {
 
-	if(argc != 2)
+	if(argc != 4)
 	{
 		std::cerr<<"Incorrect number of commandline arguments! Returning."<<std::endl;
 		return 1;
@@ -18,21 +19,26 @@ int main(int argc, char** argv)
 		return 1;
 	}
 
-	try
+	std::string inputname = argv[1];
-	{
+	std::string outputname = argv[2];
-		DetectorApp app;
+	std::string statsname = argv[3];
-		if(!app.LoadConfig(argv[1]))
-		{
-			std::cerr << "Unable to load config file " << argv[1] << ". Shutting down." << std::endl;
-			return 1;
-		}
 
-		app.Run();
+	
-	}
+	SabreArray sabre;
-	catch(std::exception& e)
+	std::string mapfile = "./etc/sabreDeadChannels_May2022.txt";
-	{
+	sabre.SetDeadChannelMap(mapfile);
-		std::cerr << e.what() << std::endl;
+	sabre.CalculateEfficiency(inputname, outputname, statsname);
-		return 1;
+	//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");
+	*/
 	return 0;
 }

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/CMakeLists.txt

@@ -50,15 +50,7 @@ target_sources(Mask PRIVATE
     ThreeStepSystem.h
     TwoStepSystem.cpp
     TwoStepSystem.h
-    ThreadPool.h
-    FileWriter.h
-    FileWriter.cpp
-    FileReader.h
-    FileReader.cpp
 )
 
-set(THREADS_PREFER_PTHREAD_FLAG On)
+target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES})
-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

@@ -1,66 +0,0 @@
-#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

@@ -1,46 +0,0 @@
-#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

@@ -1,75 +0,0 @@
-#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

@@ -1,45 +0,0 @@
-#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_resources(nullptr)
+		m_system(nullptr)
 	{
 		std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
 	}
@@ -16,8 +16,6 @@ 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) 
@@ -33,25 +31,10 @@ 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)
@@ -115,20 +98,10 @@ 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);
 
@@ -159,13 +132,9 @@ namespace Mask {
 		}
 		m_system->SetLayeredTarget(target);
 
-		for(auto system : m_systemList)
+		std::cout<<"Outputing data to file: "<<m_outputName<<std::endl;
-		{
+		std::cout<<"Reaction equation: "<<GetSystemName()<<std::endl;
-			system->SetLayeredTarget(target);
+		std::cout<<"Number of samples: "<<GetNumberOfSamples()<<std::endl;
-		}
-
-		std::cout<<"Reaction equation: "<<m_system->GetSystemEquation()<<std::endl;
-		std::cout<<"Number of samples: "<<m_nsamples<<std::endl;
 	
 		return true;
 	}
@@ -177,57 +146,6 @@ 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) 
@@ -240,11 +158,11 @@ namespace Mask {
 		tree->Branch("nuclei", m_system->GetNuclei());
 	
 		//For progress tracking
-		uint64_t percent5 = 0.05*m_nsamples;
+		uint32_t percent5 = 0.05*m_nsamples;
-		uint64_t count = 0;
+		uint32_t count = 0;
-		uint64_t npercent = 0;
+		uint32_t npercent = 0;
 	
-		for(uint64_t i=0; i<m_nsamples; i++) 
+		for(uint32_t i=0; i<m_nsamples; i++) 
 		{
 			if(++count == percent5) 
 			{

src/Mask/MaskApp.h

@@ -7,10 +7,6 @@
 #include "TwoStepSystem.h"
 #include "ThreeStepSystem.h"
 #include "RxnType.h"
-#include "ThreadPool.h"
-#include "FileWriter.h"
-
-#include <memory>
 
 namespace Mask {
 
@@ -21,24 +17,18 @@ 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;
-		uint64_t m_nsamples;
+		int 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,6 +1,7 @@
 #include "RandomGenerator.h"
 
 namespace Mask {
+	RandomGenerator* RandomGenerator::s_instance = new RandomGenerator();
 
 	RandomGenerator::RandomGenerator()
 	{

src/Mask/RandomGenerator.h

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

src/Mask/ThreadPool.h

@@ -1,118 +0,0 @@
-/*
-    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