From 0456f230f50adf08ee768ecab8593b3c0c6e193a Mon Sep 17 00:00:00 2001 From: gwm17 Date: Fri, 2 Sep 2022 17:14:30 -0400 Subject: [PATCH 1/4] Overhaul for multi-threading. Needs some testing to validate, but first signs are good. --- input.txt | 1 + src/Kinematics/main.cpp | 1 + src/Mask/CMakeLists.txt | 8 ++- src/Mask/FileWriter.cpp | 75 +++++++++++++++++++++++ src/Mask/FileWriter.h | 45 ++++++++++++++ src/Mask/MaskApp.cpp | 94 +++++++++++++++++++++++++--- src/Mask/MaskApp.h | 19 ++++-- src/Mask/RandomGenerator.cpp | 3 +- src/Mask/RandomGenerator.h | 7 ++- src/Mask/ThreadPool.h | 115 +++++++++++++++++++++++++++++++++++ 10 files changed, 348 insertions(+), 20 deletions(-) create mode 100644 src/Mask/FileWriter.cpp create mode 100644 src/Mask/FileWriter.h create mode 100644 src/Mask/ThreadPool.h diff --git a/input.txt b/input.txt index 7922dc8..7579ed3 100644 --- a/input.txt +++ b/input.txt @@ -1,5 +1,6 @@ ----------Data Information---------- OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root +NumberOfThreads: 6 ----------Reaction Information---------- NumberOfSamples: 100000 begin_chain diff --git a/src/Kinematics/main.cpp b/src/Kinematics/main.cpp index 2cec2fd..03798b2 100644 --- a/src/Kinematics/main.cpp +++ b/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) { diff --git a/src/Mask/CMakeLists.txt b/src/Mask/CMakeLists.txt index ec14d3f..d08879f 100644 --- a/src/Mask/CMakeLists.txt +++ b/src/Mask/CMakeLists.txt @@ -50,7 +50,13 @@ target_sources(Mask PRIVATE ThreeStepSystem.h TwoStepSystem.cpp TwoStepSystem.h + ThreadPool.h + FileWriter.h + FileWriter.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}) \ 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..c424965 --- /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->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 thread safe! + + void Open(const std::string& filename, const std::string& treename); + 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..32fdf31 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 +108,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({std::bind(&MaskApp::RunChunk, std::ref(*this), std::placeholders::_1), system}); + + 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; iRunSystem(); + m_fileWriter.PushData(*(system->GetNuclei())); + } + } + } diff --git a/src/Mask/MaskApp.h b/src/Mask/MaskApp.h index 62cbead..8b86cf5 100644 --- a/src/Mask/MaskApp.h +++ b/src/Mask/MaskApp.h @@ -7,6 +7,10 @@ #include "TwoStepSystem.h" #include "ThreeStepSystem.h" #include "RxnType.h" +#include "ThreadPool.h" +#include "FileWriter.h" + +#include namespace Mask { @@ -17,18 +21,23 @@ 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 + 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..babfaa4 --- /dev/null +++ b/src/Mask/ThreadPool.h @@ -0,0 +1,115 @@ +/* + 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 { + + using JobFunction = std::function; + struct Job + { + JobFunction func; + ReactionSystem* argument; + }; + + class ThreadPool + { + public: + 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--; + + 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 From 2f7b0969dc467411b18045809305ddd9e4a9bd62 Mon Sep 17 00:00:00 2001 From: gwm17 Date: Wed, 7 Sep 2022 11:29:57 -0400 Subject: [PATCH 2/4] Implemented multithreading in detectors, overhaul application api for detectors, implement detector input file. Need robust testing to verify performance gain in detectors. Verified that code compiles and executes. --- detector_input.txt | 5 + kinematics_input.txt | 45 ++++++ src/Detectors/AnasenArray.cpp | 240 +------------------------------- src/Detectors/AnasenArray.h | 16 +-- src/Detectors/CMakeLists.txt | 3 + src/Detectors/DetectorApp.cpp | 144 +++++++++++++++++++ src/Detectors/DetectorApp.h | 40 ++++++ src/Detectors/DetectorArray.cpp | 35 +++++ src/Detectors/DetectorArray.h | 17 ++- src/Detectors/SabreArray.cpp | 237 +------------------------------ src/Detectors/SabreArray.h | 8 +- src/Detectors/main.cpp | 40 +++--- src/Mask/CMakeLists.txt | 2 + src/Mask/FileReader.cpp | 66 +++++++++ src/Mask/FileReader.h | 46 ++++++ src/Mask/FileWriter.h | 6 +- src/Mask/MaskApp.cpp | 42 +++--- src/Mask/MaskApp.h | 3 +- src/Mask/ThreadPool.h | 19 +-- 19 files changed, 476 insertions(+), 538 deletions(-) create mode 100644 detector_input.txt create mode 100644 kinematics_input.txt create mode 100644 src/Detectors/DetectorApp.cpp create mode 100644 src/Detectors/DetectorApp.h create mode 100644 src/Detectors/DetectorArray.cpp create mode 100644 src/Mask/FileReader.cpp create mode 100644 src/Mask/FileReader.h 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/kinematics_input.txt b/kinematics_input.txt new file mode 100644 index 0000000..fdbe62c --- /dev/null +++ b/kinematics_input.txt @@ -0,0 +1,45 @@ +----------Data Information---------- +OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root +NumberOfThreads: 6 +----------Reaction Information---------- +NumberOfSamples: 1000000 +begin_chain + begin_step + Type: Reaction + begin_nuclei + 5 10 + 2 3 + 2 4 + end_nuclei + BeamEnergyMean(MeV): 24.0 + BeamEnergySigma(MeV): 0.0 + ThetaType: Lab + ThetaMin(deg): 15.0 + ThetaMax(deg): 15.0 + PhiMin(deg): 0.0 + PhMax(deg): 0.0 + ResidualExcitationMean(MeV): 16.8 + ResidualExcitationSigma(MeV): 0.023 + end_step + begin_step + Type: Decay + begin_nuclei + 5 9 + 2 4 + end_nuclei + PhiMin(deg): 0.0 + PhMax(deg): 360.0 + ResidualExcitationMean(MeV): 0.0 + ResidualExcitationSigma(MeV): 0.0 + AngularDistributionFile: ./etc/isotropic_dist.txt + end_step +end_chain +----------Target Information---------- +NumberOfLayers: 1 +begin_layer + Thickness(ug/cm^2): 74 + begin_elements (Z, A, Stoich.) + element 5 10 1 + end_elements +end_layer +-------------------------------------- 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.h b/src/Mask/FileWriter.h index c424965..7ed188f 100644 --- a/src/Mask/FileWriter.h +++ b/src/Mask/FileWriter.h @@ -19,15 +19,15 @@ namespace Mask { FileWriter(const std::string& filename, const std::string& treename); ~FileWriter(); - bool IsOpen() const { return m_file->IsOpen(); } + 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 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); + void Open(const std::string& filename, const std::string& treename); //Not thread safe! void Close(); //Not thread safe! private: diff --git a/src/Mask/MaskApp.cpp b/src/Mask/MaskApp.cpp index 32fdf31..7e76d14 100644 --- a/src/Mask/MaskApp.cpp +++ b/src/Mask/MaskApp.cpp @@ -39,9 +39,16 @@ namespace Mask { input>>junk>>m_nsamples; std::cout<<"Allocating resources... Asking for " << m_nthreads << " threads..."; - m_resources = std::make_unique(m_nthreads); + m_resources = std::make_unique>(m_nthreads); std::cout<<" Complete."<PushJob({std::bind(&MaskApp::RunChunk, std::ref(*this), std::placeholders::_1), system}); + for(std::size_t i=0; iPushJob({[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; @@ -243,18 +265,4 @@ namespace Mask { std::cout<<"---------------------------------------------"<RunSystem(); - m_fileWriter.PushData(*(system->GetNuclei())); - } - } - } diff --git a/src/Mask/MaskApp.h b/src/Mask/MaskApp.h index 8b86cf5..d5457fd 100644 --- a/src/Mask/MaskApp.h +++ b/src/Mask/MaskApp.h @@ -36,8 +36,9 @@ namespace Mask { 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; + std::unique_ptr> m_resources; }; diff --git a/src/Mask/ThreadPool.h b/src/Mask/ThreadPool.h index babfaa4..aa09978 100644 --- a/src/Mask/ThreadPool.h +++ b/src/Mask/ThreadPool.h @@ -20,16 +20,19 @@ namespace Mask { - using JobFunction = std::function; - struct Job - { - JobFunction func; - ReactionSystem* argument; - }; - + 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) { @@ -94,7 +97,7 @@ namespace Mask { m_numberRunning++; m_queueSize--; - job.func(job.argument); + std::apply(job.func, job.argument); m_numberRunning--; } From 3db9e510f518f684956690de5acd82b048a947c9 Mon Sep 17 00:00:00 2001 From: Gordon McCann <43148247+gwm17@users.noreply.github.com> Date: Wed, 7 Sep 2022 11:31:50 -0400 Subject: [PATCH 3/4] Remove obsolete input file --- input.txt | 45 --------------------------------------------- 1 file changed, 45 deletions(-) delete mode 100644 input.txt diff --git a/input.txt b/input.txt deleted file mode 100644 index 7579ed3..0000000 --- a/input.txt +++ /dev/null @@ -1,45 +0,0 @@ -----------Data Information---------- -OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root -NumberOfThreads: 6 -----------Reaction Information---------- -NumberOfSamples: 100000 -begin_chain - begin_step - Type: Reaction - begin_nuclei - 5 10 - 2 3 - 2 4 - end_nuclei - BeamEnergyMean(MeV): 24.0 - BeamEnergySigma(MeV): 0.0 - ThetaType: Lab - ThetaMin(deg): 15.0 - ThetaMax(deg): 15.0 - PhiMin(deg): 0.0 - PhMax(deg): 0.0 - ResidualExcitationMean(MeV): 16.8 - ResidualExcitationSigma(MeV): 0.023 - end_step - begin_step - Type: Decay - begin_nuclei - 5 9 - 2 4 - end_nuclei - PhiMin(deg): 0.0 - PhMax(deg): 360.0 - ResidualExcitationMean(MeV): 0.0 - ResidualExcitationSigma(MeV): 0.0 - AngularDistributionFile: ./etc/isotropic_dist.txt - end_step -end_chain -----------Target Information---------- -NumberOfLayers: 1 -begin_layer - Thickness(ug/cm^2): 74 - begin_elements (Z, A, Stoich.) - element 5 10 1 - end_elements -end_layer --------------------------------------- From 8b81dda70e80747338eccf793219909ee5cb94a8 Mon Sep 17 00:00:00 2001 From: gwm17 Date: Wed, 7 Sep 2022 11:35:50 -0400 Subject: [PATCH 4/4] Update REAMDE. --- README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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.