diff --git a/.gitmodules b/.gitmodules index 0f578de..cde4811 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,6 @@ [submodule "src/vendor/catima"] path = src/vendor/catima url = https://github.com/gwm17/catima.git +[submodule "src/vendor/yaml-cpp"] + path = src/vendor/yaml-cpp + url = https://github.com/jbeder/yaml-cpp.git diff --git a/CMakeLists.txt b/CMakeLists.txt index ebf7578..e02f867 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,6 +19,8 @@ set(CMAKE_CXX_STANDARD 17) find_package(ROOT REQUIRED COMPONENTS GenVector) add_subdirectory(src/vendor/catima) +set(YAML_CPP_BUILD_TOOLS Off CACHE BOOL "Enable parse tools") +add_subdirectory(src/vendor/yaml-cpp) add_subdirectory(src/Mask) add_subdirectory(src/Kinematics) add_subdirectory(src/Detectors) diff --git a/detector.yaml b/detector.yaml new file mode 100644 index 0000000..48d2a78 --- /dev/null +++ b/detector.yaml @@ -0,0 +1,5 @@ +InputDataFile: /media/data/gwm17/mask_tests/temp.root +OutputDataFile: /media/data/gwm17/mask_tests/temp_det.root +DeadChannelFile: etc/sabreDeadChannels_May2022.txt +NumberOfThreads: 5 +ArrayType: Sabre diff --git a/detector_input.txt b/detector_input.txt deleted file mode 100644 index 775494d..0000000 --- a/detector_input.txt +++ /dev/null @@ -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 \ No newline at end of file diff --git a/kinematics.yaml b/kinematics.yaml new file mode 100644 index 0000000..d743212 --- /dev/null +++ b/kinematics.yaml @@ -0,0 +1,49 @@ +OutputFile: /media/data/gwm17/mask_tests/temp.root +Threads: 5 +ReactionSamples: 1000000 +ReactionChain: + - Type: Reaction + Reactants: + - Z: 5 + A: 10 + - Z: 2 + A: 3 + - Z: 2 + A: 4 + BeamEnergyMean(MeV): 24.0 + BeamEnergySigma(MeV): 0.0 + ThetaType: Lab + ThetaMin(deg): 15.0 + ThetaMax(deg): 15.0 + PhiMin(deg): 0.0 + PhiMax(deg): 0.0 + ResidualExcitationMean(MeV): 16.798 + ResidualExcitationSigma(MeV): 0.033 + - Type: Decay + Reactants: + - Z: 5 + A: 9 + - Z: 2 + A: 4 + PhiMin(deg): 0.0 + PhiMax(deg): 360.0 + ResidualExcitationMean(MeV): 0.0 + ResidualExcitationSigma(MeV): 0.522 + AngularDistributionFile: ./etc/isotropic_dist.txt + - Type: Decay + Reactants: + - Z: 3 + A: 5 + - Z: 2 + A: 4 + PhiMin(deg): 0.0 + PhiMax(deg): 360.0 + ResidualExcitationMean(MeV): 0.0 + ResidualExcitationSigma(MeV): 0.0 + AngularDistributionFile: ./etc/isotropic_dist.txt +TargetLayers: + - Thickness(ug/cm^2): 74 + Elements: + - Z: 5 + A: 10 + S: 1 \ No newline at end of file diff --git a/kinematics_input.txt b/kinematics_input.txt deleted file mode 100644 index fdbe62c..0000000 --- a/kinematics_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: 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/CMakeLists.txt b/src/Detectors/CMakeLists.txt index b556f58..fe66023 100644 --- a/src/Detectors/CMakeLists.txt +++ b/src/Detectors/CMakeLists.txt @@ -1,5 +1,9 @@ add_executable(Detectors) -target_include_directories(Detectors PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_include_directories(Detectors PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/.. + ${CMAKE_CURRENT_SOURCE_DIR}/../vendor/yaml-cpp/include/ +) target_sources(Detectors PUBLIC AnasenDeadChannelMap.cpp diff --git a/src/Detectors/DetectorApp.cpp b/src/Detectors/DetectorApp.cpp index 766ecab..b2d3be0 100644 --- a/src/Detectors/DetectorApp.cpp +++ b/src/Detectors/DetectorApp.cpp @@ -2,6 +2,8 @@ #include #include +#include "yaml-cpp/yaml.h" + DetectorApp::DetectorApp() : m_resources(nullptr) { @@ -13,6 +15,65 @@ DetectorApp::~DetectorApp() delete m_detectorList[i]; } +bool DetectorApp::ParseConfig(const std::string& filename) +{ + std::cout<<"----------Detector Efficiency Calculation----------"<(); + m_outputFileName = data["OutputDataFile"].as(); + m_deadChannelFileName = data["DeadChannelFile"].as(); + m_nthreads = data["NumberOfThreads"].as(); + ArrayType type = StringToArrayType(data["ArrayType"].as()); + + std::cout << "Creating " << m_nthreads << " detector arrays..." << std::endl; + for(uint64_t i=0; iSetDeadChannelMap(m_deadChannelFileName); + } + + std::cout << "Allocating " << m_nthreads << " threads..." << std::endl; + m_resources = std::make_unique>(m_nthreads); + + std::cout << "Opening input data file " << m_inputFileName << "..." << std::endl; + m_fileReader.Open(m_inputFileName, "SimTree"); + if(!m_fileReader.IsOpen() || !m_fileReader.IsTree()) + { + std::cerr << "Unable to open input data file " << m_inputFileName << std::endl; + return false; + } + m_nentries = m_fileReader.GetSize(); + std::cout << "Detected " << m_nentries << " events in the file..." << std::endl; + + //Little bit of integer division mangling to make sure we read every event in file + uint64_t quotient = m_nentries / m_nthreads; + uint64_t remainder = m_nentries % m_nthreads; + m_chunkSamples.push_back(quotient + remainder); + for(uint64_t i=1; i()); + YAML::Node nuclei = yamlStream["Reactants"]; + for(auto nucleus : nuclei) + { + params.Z.push_back(nucleus["Z"].as()); + params.A.push_back(nucleus["A"].as()); + } + params.phiMin = yamlStream["PhiMin(deg)"].as(); + params.phiMax = yamlStream["PhiMax(deg)"].as(); + params.meanResidualEx = yamlStream["ResidualExcitationMean(MeV)"].as(); + params.sigmaResidualEx = yamlStream["ResidualExcitationSigma(MeV)"].as(); + params.angularDistFile = yamlStream["AngularDistributionFile"].as(); + return params; + } + + static void SerializeReaction(YAML::Emitter& yamlStream, const StepParameters& params) + { + yamlStream << YAML::BeginMap; + yamlStream << YAML::Key << "Type" << RxnTypeToString(params.rxnType); + yamlStream << YAML::Key << "Reactants" << YAML::Value << YAML::BeginSeq; + for (std::size_t i=0; i()); + YAML::Node nuclei = yamlStream["Reactants"]; + for(auto nucleus : nuclei) + { + params.Z.push_back(nucleus["Z"].as()); + params.A.push_back(nucleus["A"].as()); + } + params.meanBeamEnergy = yamlStream["BeamEnergyMean(MeV)"].as(); + params.sigmaBeamEnergy = yamlStream["BeamEnergySigma(MeV)"].as(); + params.thetaType = StringToRxnThetaType(yamlStream["ThetaType"].as()); + params.thetaMin = yamlStream["ThetaMin(deg)"].as(); + params.thetaMax = yamlStream["ThetaMax(deg)"].as(); + params.phiMin = yamlStream["PhiMin(deg)"].as(); + params.phiMax = yamlStream["PhiMax(deg)"].as(); + params.meanResidualEx = yamlStream["ResidualExcitationMean(MeV)"].as(); + params.sigmaResidualEx = yamlStream["ResidualExcitationSigma(MeV)"].as(); + return params; + } + + static void SerializeTarget(YAML::Emitter& yamlStream, const LayeredTarget& target) + { + + yamlStream << YAML::BeginSeq; + for (int i=0; i z, a, s; + double thickness; + + for (auto layer : yamlStream) + { + thickness = layer["Thickness(ug/cm^2)"].as(); + YAML::Node elements = layer["Elements"]; + z.clear(); + a.clear(); + s.clear(); + for(auto element : elements) + { + z.push_back(element["Z"].as()); + a.push_back(element["A"].as()); + s.push_back(element["S"].as()); + } + target.AddLayer(z, a, s, thickness); + } + + return target; + } + + bool ConfigSerializer::SerializeConfig(const std::string& configfile, const AppParameters& params) + { + std::ofstream output(configfile); + if (!output.is_open()) + { + std::cerr << "Could not open output file " << configfile << std::endl; + return false; + } + + YAML::Emitter yamlStream; + yamlStream << YAML::BeginMap; + yamlStream << YAML::Key << "OutputFile" << YAML::Value << params.outputFileName; + yamlStream << YAML::Key << "Threads" << YAML::Value << params.nThreads; + yamlStream << YAML::Key << "ReactionSamples" << YAML::Value << params.nSamples; + yamlStream << YAML::Key << "ReactionChain" << YAML::Value << YAML::BeginSeq; + for (auto& step : params.chainParams) + { + switch(step.rxnType) + { + case RxnType::Reaction: SerializeReaction(yamlStream, step); break; + case RxnType::Decay: SerializeDecay(yamlStream, step); break; + case RxnType::None: + { + std::cerr << "Error serializing config: none type reaction found!" << std::endl; + return false; + } + } + } + yamlStream << YAML::EndSeq; + yamlStream << YAML::Key << "TargetLayers" << YAML::Value; + SerializeTarget(yamlStream, params.target); + + output << yamlStream.c_str(); + output.close(); + return true; + } + + bool ConfigSerializer::DeserializeConfig(const std::string& configfile, AppParameters& params) + { + YAML::Node data; + try + { + data = YAML::LoadFile(configfile); + } + catch(YAML::ParserException& e) + { + std::cerr << "Error reading config " << configfile << ": " << e.what() <(); + params.nThreads = data["Threads"].as(); + params.nSamples = data["ReactionSamples"].as(); + + auto steps = data["ReactionChain"]; + for (auto step : steps) + { + RxnType type = StringToRxnType(step["Type"].as()); + switch(type) + { + case RxnType::Decay: + { + params.chainParams.push_back(DeserializeDecay(step)); + break; + } + case RxnType::Reaction: + { + params.chainParams.push_back(DeserializeReaction(step)); + break; + } + case RxnType::None: + { + std::cerr << "Error deserializing config: None-type reaction found" << std::endl; + return false; + } + } + } + auto layers = data["TargetLayers"]; + params.target = DeserializeTarget(layers); + return true; + } +} \ No newline at end of file diff --git a/src/Mask/ConfigSerializer.h b/src/Mask/ConfigSerializer.h new file mode 100644 index 0000000..7b0ddd2 --- /dev/null +++ b/src/Mask/ConfigSerializer.h @@ -0,0 +1,16 @@ +#ifndef CONFIG_SERIALIZER_H +#define CONFIG_SERIALIZER_H + +#include "MaskApp.h" + +namespace Mask { + + class ConfigSerializer + { + public: + static bool SerializeConfig(const std::string& configfile, const AppParameters& params); + static bool DeserializeConfig(const std::string& configfile, AppParameters& params); + }; +} + +#endif \ No newline at end of file diff --git a/src/Mask/LayeredTarget.h b/src/Mask/LayeredTarget.h index cf6dcf7..242cbcd 100644 --- a/src/Mask/LayeredTarget.h +++ b/src/Mask/LayeredTarget.h @@ -31,9 +31,9 @@ namespace Mask { double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth); double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth); std::size_t FindLayerContaining(int Z, int A); - std::size_t GetNumberOfLayers() { return m_layers.size(); } + std::size_t GetNumberOfLayers() const { return m_layers.size(); } void SetName(std::string& n) { m_name = n; } - const Target& GetLayerInfo(int index) { return m_layers[index]; } + const Target& GetLayerInfo(int index) const { return m_layers[index]; } const std::string& GetName() { return m_name; } private: diff --git a/src/Mask/MaskApp.cpp b/src/Mask/MaskApp.cpp index 7e76d14..5806ca1 100644 --- a/src/Mask/MaskApp.cpp +++ b/src/Mask/MaskApp.cpp @@ -1,185 +1,80 @@ #include "MaskApp.h" +#include "ConfigSerializer.h" #include #include #include "TFile.h" #include "TTree.h" +#include "yaml-cpp/yaml.h" + namespace Mask { MaskApp::MaskApp() : - m_system(nullptr), m_resources(nullptr) + m_resources(nullptr) { std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<>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) + //Initialize the system + //Reaction chain + for(uint32_t i=0; i> junk >> junk; - currentParams.rxnType = StringToRxnType(junk); - if(currentParams.rxnType == RxnType::Reaction) - { - input >> junk; - for(int i=0; i<3; i++) - { - input >> z >> a; - currentParams.Z.push_back(z); - currentParams.A.push_back(a); - } - input >> junk; - input >> junk >> currentParams.meanBeamEnergy; - input >> junk >> currentParams.sigmaBeamEnergy; - input >> junk >> junk; - currentParams.thetaType = StringToRxnThetaType(junk); - input >> junk >> currentParams.thetaMin; - input >> junk >> currentParams.thetaMax; - input >> junk >> currentParams.phiMin; - input >> junk >> currentParams.phiMax; - input >> junk >> currentParams.meanResidualEx; - input >> junk >> currentParams.sigmaResidualEx; - params.push_back(currentParams); - } - else if(currentParams.rxnType == RxnType::Decay) - { - input >> junk; - for(int i=0; i<2; i++) - { - input >> z >> a; - currentParams.Z.push_back(z); - currentParams.A.push_back(a); - } - input >> junk; - input >> junk >> currentParams.phiMin; - input >> junk >> currentParams.phiMax; - input >> junk >> currentParams.meanResidualEx; - input >> junk >> currentParams.sigmaResidualEx; - input >> junk >> currentParams.angularDistFile; - params.push_back(currentParams); - } - else - { - std::cerr << "Invalid reaction information at MaskApp::LoadConfig!" << std::endl; - return false; - } - } - } - - m_system = CreateSystem(params); - if(m_system == nullptr || !m_system->IsValid()) - { - std::cerr<<"Failure to parse reaction system... configuration not loaded."<IsValid()) { std::cerr<<"Failure to parse reaction system... configuration not loaded."< avec, zvec, svec; - int s; - input>>junk>>nlayers; - for(int i=0; i>junk>>junk>>thickness; - avec.clear(); zvec.clear(); svec.clear(); - while(input>>junk) - { - if(junk == "begin_elements") - { - input>>junk>>junk>>junk; - continue; - } - else if (junk == "end_elements") - break; - input>>z>>a>>s; - zvec.push_back(z); avec.push_back(a); svec.push_back(s); - } - target.AddLayer(zvec, avec, svec, thickness); - input>>junk; - } - m_system->SetLayeredTarget(target); - + //Link the target for(auto system : m_systemList) { - system->SetLayeredTarget(target); + system->SetLayeredTarget(m_params.target); } + //Setup threading + m_resources = std::make_unique>(m_params.nThreads); + //Little bit of integer division mangling to make sure we do the total number of samples + uint64_t quotient = m_params.nSamples / m_params.nThreads; + uint64_t remainder = m_params.nSamples % m_params.nThreads; + m_chunkSamples.push_back(quotient + remainder); + for(uint64_t i=1; iGetSystemEquation()<GetSystemEquation() << std::endl; + std::cout << "Number of samples: " << m_params.nSamples << std::endl; + std::cout << "Number of threads: " << m_params.nThreads << std::endl; + std::cout << "Outputing data to file: " << m_params.outputFileName << std::endl; return true; } - - //Not implemented... yet! + bool MaskApp::SaveConfig(const std::string& filename) - { - return true; + { + std::cout << "Writing configuration to " << filename << std::endl; + return ConfigSerializer::SerializeConfig(filename, m_params); } - + void MaskApp::Run() { std::cout<<"Running simulation..."<Branch("nuclei", m_system->GetNuclei()); - - //For progress tracking - uint64_t percent5 = 0.05*m_nsamples; - uint64_t count = 0; - uint64_t npercent = 0; - - for(uint64_t i=0; iRunSystem(); - tree->Fill(); - } - - tree->Write(tree->GetName(), TObject::kOverwrite); - output->Close(); - - std::cout< chainParams; + LayeredTarget target; + }; + class MaskApp { public: @@ -23,17 +32,9 @@ namespace Mask { bool SaveConfig(const std::string& filename); void Run(); - void RunSingleThread(); private: - void RunChunk(ReactionSystem* system); - - ReactionSystem* m_system; - - std::string m_outputName; - RxnType m_rxnType; - uint64_t m_nsamples; - uint32_t m_nthreads; + AppParameters m_params; std::vector m_systemList; //One system for each thread std::vector m_chunkSamples; diff --git a/src/Mask/RxnType.h b/src/Mask/RxnType.h index cb46d59..b868758 100644 --- a/src/Mask/RxnType.h +++ b/src/Mask/RxnType.h @@ -37,6 +37,17 @@ namespace Mask { return RxnType::None; } + static std::string RxnTypeToString(RxnType type) + { + switch(type) + { + case RxnType::Decay: return "Decay"; + case RxnType::Reaction: return "Reaction"; + case RxnType::None: return "None"; + default: return "None"; + } + } + static RxnThetaType StringToRxnThetaType(const std::string& type) { if(type == "CenterOfMass") @@ -52,6 +63,17 @@ namespace Mask { return RxnThetaType::None; } } + + static std::string RxnThetaTypeToString(RxnThetaType type) + { + switch(type) + { + case RxnThetaType::CenterOfMass: return "CenterOfMass"; + case RxnThetaType::Lab: return "Lab"; + case RxnThetaType::None: return "None"; + default: return "None"; + } + } } #endif \ No newline at end of file diff --git a/src/Mask/Target.h b/src/Mask/Target.h index be8b4a4..6f3a8b1 100644 --- a/src/Mask/Target.h +++ b/src/Mask/Target.h @@ -31,12 +31,12 @@ namespace Mask { double GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle); double GetEnergyLossFractionalDepth(int zp, int ap, double startEnergy, double angle, double percent_depth); double GetReverseEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double angle, double percent_depth); - inline const double& GetThickness() { return m_thickness; } - inline int GetNumberOfElements() { return m_Z.size(); } - inline int GetElementZ(int index) { return m_Z[index]; } - inline int GetElementA(int index) { return m_A[index]; } - inline int GetElementStoich(int index) { return m_stoich[index]; } - + const double& GetThickness() const { return m_thickness; } + int GetNumberOfElements() const { return m_Z.size(); } + int GetElementZ(int index) const { return m_Z[index]; } + int GetElementA(int index) const { return m_A[index]; } + int GetElementStoich(int index) const { return m_stoich[index]; } + private: void Init(const std::vector& z, const std::vector& a, const std::vector& stoich); diff --git a/src/vendor/yaml-cpp b/src/vendor/yaml-cpp new file mode 160000 index 0000000..0e6e28d --- /dev/null +++ b/src/vendor/yaml-cpp @@ -0,0 +1 @@ +Subproject commit 0e6e28d1a38224fc8172fae0109ea7f673c096db