diff --git a/README.md b/README.md index cb1950b..6e30a0a 100644 --- a/README.md +++ b/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 ` +`./bin/Detectors ` -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. diff --git a/detector_input.txt b/detector_input.txt new file mode 100644 index 0000000..775494d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/input.txt b/kinematics_input.txt similarity index 95% rename from input.txt rename to kinematics_input.txt index 7922dc8..fdbe62c 100644 --- a/input.txt +++ b/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 diff --git a/src/Detectors/AnasenArray.cpp b/src/Detectors/AnasenArray.cpp index 7e31045..56598db 100644 --- a/src/Detectors/AnasenArray.cpp +++ b/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) @@ -337,236 +337,4 @@ DetectorResult AnasenArray::IsAnasen(Mask::Nucleus& nucleus) if(!result.detectFlag) result = IsQQQ(nucleus); return result; -} - -void AnasenArray::CountCoincidences(const std::vector& data, std::vector& 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."<Get("SimTree"); - std::vector* dataHandle = new std::vector(); - intree->SetBranchAddress("nuclei", &dataHandle); - - output->cd(); - TTree* outtree = new TTree("SimTree", "SimTree"); - outtree->Branch("nuclei", dataHandle); - - input->cd(); - - stats<<"Efficiency statistics for data from "<GetEntry(1); - std::vector counts; - std::vector 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: "<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; iGetEntry(i); - if(++count == percent5) - {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<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<size() == 3) - { - stats<size() == 4) - { - stats<size() == 6) - { - stats<size() == 8) - { - stats<& data, std::vector& counts); + DetectorResult IsRing1(const Mask::Nucleus& nucleus); + DetectorResult IsRing2(const Mask::Nucleus& nucleus); + DetectorResult IsQQQ(const Mask::Nucleus& nucleus); std::vector m_Ring1; std::vector m_Ring2; diff --git a/src/Detectors/CMakeLists.txt b/src/Detectors/CMakeLists.txt index 12cb098..b556f58 100644 --- a/src/Detectors/CMakeLists.txt +++ b/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 diff --git a/src/Detectors/DetectorApp.cpp b/src/Detectors/DetectorApp.cpp new file mode 100644 index 0000000..83e3f74 --- /dev/null +++ b/src/Detectors/DetectorApp.cpp @@ -0,0 +1,144 @@ +#include "DetectorApp.h" +#include +#include + +DetectorApp::DetectorApp() : + m_resources(nullptr) +{ +} + +DetectorApp::~DetectorApp() +{ + for(std::size_t i=0; i> 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; iSetDeadChannelMap(m_deadChannelFileName); + } + std::cout << "Done" << std::endl; + + std::cout << "Allocating " << m_nthreads << " threads..." << std::endl; + m_resources = std::make_unique>(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; iPushJob({[this](DetectorArray* array, uint64_t chunkSamples) + { + if(system == nullptr) + return; + + std::vector data; + DetectorResult result; + for(uint64_t i=0; iIsDetected(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: "<IsFinished() && m_fileWriter.GetQueueSize() == 0) + break; + else if(m_fileWriter.Write()) + ++count; + } + +} \ No newline at end of file diff --git a/src/Detectors/DetectorApp.h b/src/Detectors/DetectorApp.h new file mode 100644 index 0000000..e444dc2 --- /dev/null +++ b/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 +#include +#include + +class DetectorApp +{ +public: + DetectorApp(); + ~DetectorApp(); + + bool LoadConfig(const std::string& filename); + + void Run(); + +private: + + std::vector m_detectorList; //One array per thread + std::vector 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> m_resources; +}; + +#endif \ No newline at end of file diff --git a/src/Detectors/DetectorArray.cpp b/src/Detectors/DetectorArray.cpp new file mode 100644 index 0000000..bd7ddab --- /dev/null +++ b/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; +} \ No newline at end of file diff --git a/src/Detectors/DetectorArray.h b/src/Detectors/DetectorArray.h index 5e7f37c..baf8b57 100644 --- a/src/Detectors/DetectorArray.h +++ b/src/Detectors/DetectorArray.h @@ -5,6 +5,7 @@ #include #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 \ No newline at end of file diff --git a/src/Detectors/SabreArray.cpp b/src/Detectors/SabreArray.cpp index 01a41ae..49c9a01 100644 --- a/src/Detectors/SabreArray.cpp +++ b/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----------"<Get("SimTree"); - std::vector* dataHandle = new std::vector(); - intree->SetBranchAddress("nuclei", &dataHandle); - - output->cd(); - TTree* outtree = new TTree("SimTree", "SimTree"); - outtree->Branch("nuclei", dataHandle); - - input->cd(); - - stats<<"Efficiency statistics for data from "<GetEntry(1); - std::vector counts; - std::vector 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: "<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; iGetEntry(i); - if(++count == percent5) - {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<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<& data, std::vector& 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]++; - } - } -} +} \ No newline at end of file diff --git a/src/Detectors/SabreArray.h b/src/Detectors/SabreArray.h index 2ae54bc..b32ca6c 100644 --- a/src/Detectors/SabreArray.h +++ b/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& data, std::vector& counts); std::vector m_detectors; diff --git a/src/Detectors/main.cpp b/src/Detectors/main.cpp index 8cfb3b0..6c06e87 100644 --- a/src/Detectors/main.cpp +++ b/src/Detectors/main.cpp @@ -1,5 +1,4 @@ -#include "SabreArray.h" -#include "AnasenArray.h" +#include "DetectorApp.h" #include "KinematicsExceptions.h" #include #include @@ -7,7 +6,7 @@ int main(int argc, char** argv) { - if(argc != 4) + if(argc != 2) { std::cerr<<"Incorrect number of commandline arguments! Returning."<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(); + 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& dataHandle) + { + std::scoped_lock guard(m_fileMutex); + int bytes = m_tree->GetEntry(m_currentEntry); + if(bytes != 0) + { + dataHandle = *m_branchHandle; + m_currentEntry++; + return true; + } + + return false; + } +} \ No newline at end of file diff --git a/src/Mask/FileReader.h b/src/Mask/FileReader.h new file mode 100644 index 0000000..88051bb --- /dev/null +++ b/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 +#include +#include +#include + +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& 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* m_branchHandle; + + std::mutex m_fileMutex; + std::atomic m_currentEntry; + std::atomic m_size; //in entries + }; +} + +#endif \ No newline at end of file diff --git a/src/Mask/FileWriter.cpp b/src/Mask/FileWriter.cpp new file mode 100644 index 0000000..3ad7da7 --- /dev/null +++ b/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& data) + { + std::scoped_lock 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 guard(m_queueMutex); + m_dataHandle = m_queue.front(); + m_queue.pop(); + } + --m_queueSize; + + m_tree->Fill(); + return true; + } +} \ No newline at end of file diff --git a/src/Mask/FileWriter.h b/src/Mask/FileWriter.h new file mode 100644 index 0000000..7ed188f --- /dev/null +++ b/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 +#include +#include + +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& 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 m_dataHandle; + + std::mutex m_queueMutex; + std::atomic m_queueSize; + std::queue> m_queue; + }; +} + +#endif \ No newline at end of file diff --git a/src/Mask/MaskApp.cpp b/src/Mask/MaskApp.cpp index bc27127..7e76d14 100644 --- a/src/Mask/MaskApp.cpp +++ b/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----------"<>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>(m_nthreads); + std::cout<<" Complete."< 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: "<IsValid()) + { + std::cerr<<"Failure to parse reaction system... configuration not loaded."<SetLayeredTarget(target); - std::cout<<"Outputing data to file: "<SetLayeredTarget(target); + } + + std::cout<<"Reaction equation: "<GetSystemEquation()<PushJob({[this](ReactionSystem* system, uint64_t chunkSamples) + { + if(system == nullptr) + return; + + for(uint64_t i=0; iRunSystem(); + 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: "<IsFinished() && m_fileWriter.GetQueueSize() == 0) + break; + else if(m_fileWriter.Write()) + ++count; + } + + std::cout<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 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 m_systemList; //One system for each thread + std::vector m_chunkSamples; + FileWriter m_fileWriter; + std::unique_ptr> m_resources; }; diff --git a/src/Mask/RandomGenerator.cpp b/src/Mask/RandomGenerator.cpp index 1153103..df96cf6 100644 --- a/src/Mask/RandomGenerator.cpp +++ b/src/Mask/RandomGenerator.cpp @@ -1,8 +1,7 @@ #include "RandomGenerator.h" namespace Mask { - RandomGenerator* RandomGenerator::s_instance = new RandomGenerator(); - + RandomGenerator::RandomGenerator() { std::random_device rd; diff --git a/src/Mask/RandomGenerator.h b/src/Mask/RandomGenerator.h index fd8a0b2..f23ac6b 100644 --- a/src/Mask/RandomGenerator.h +++ b/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; }; diff --git a/src/Mask/ThreadPool.h b/src/Mask/ThreadPool.h new file mode 100644 index 0000000..aa09978 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include + + +namespace Mask { + + template + class ThreadPool + { + public: + + using JobFunction = std::function; + struct Job + { + JobFunction func; + std::tuple argument; + }; + + + ThreadPool(uint32_t nthreads) : + m_isStopped(false), m_numberRunning(0), m_queueSize(0), m_initShutdown(false) + { + for(uint32_t i=0; i 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 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> m_queue; + std::vector m_pool; + std::mutex m_poolMutex; + std::condition_variable m_wakeCondition; + + bool m_isStopped; + std::atomic m_numberRunning; + std::atomic m_queueSize; + std::atomic m_initShutdown; + }; +} + +#endif \ No newline at end of file