mirror of
https://github.com/gwm17/Mask.git
synced 2024-11-22 18:28:51 -05:00
Compare commits
No commits in common. "142a2dffd103875c01b7466d0d369e4511409208" and "d68f5de383e326a65f88f03c34150c12d85608fb" have entirely different histories.
142a2dffd1
...
d68f5de383
1
.gitignore
vendored
1
.gitignore
vendored
|
@ -19,7 +19,6 @@ build/
|
|||
lib/
|
||||
bin/
|
||||
scripts/
|
||||
*.yaml
|
||||
|
||||
|
||||
###Keep this file###
|
||||
|
|
3
.gitmodules
vendored
3
.gitmodules
vendored
|
@ -1,6 +1,3 @@
|
|||
[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
|
||||
|
|
|
@ -19,8 +19,6 @@ 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)
|
||||
|
|
18
README.md
18
README.md
|
@ -1,12 +1,9 @@
|
|||
# Mask: Monte cArlo Simulation of Kinematics
|
||||
|
||||
Mask is a Monte Carlo simulation of reaction kinematics for use detector systems at Florida State University.
|
||||
Mask is capable of simulating multi-step kinematic reaction-decay sequences, storing data in ROOT Trees, after which the kinematic data can be fed to a detector geometry for efficiency testing. Currently geometries for ANASEN and SABRE are included in the code.
|
||||
|
||||
## Building Mask
|
||||
|
||||
Dowload the repository from GitHub using `git clone --recursive https://github.com/gwm17/Mask.git`. CMake is use to build the project; in most environments you can build Mask using the following methods:
|
||||
|
||||
Dowload the repository from github. CMake is use to build the project; in most environments you can build Mask using the following methods:
|
||||
- `mkdir build`
|
||||
- `cd build`
|
||||
- `cmake ..`
|
||||
|
@ -17,7 +14,6 @@ Executables will be installed to the repostiory's `bin` directory. Libraries wil
|
|||
By default Mask builds for release. To build for debug replace `cmake ..` with `cmake -DCMAKE_BUILD_TYPE=Debug ..`. Mask uses CMake to find the installed ROOT libraries and headers.
|
||||
|
||||
## Using the kinematics simulation
|
||||
|
||||
By default Mask is capable of simulating reactions of up to three steps. Here is a brief outline of each type:
|
||||
|
||||
0. A pure decay involves a "target" decaying into an ejectile and residual. It is assumed isotropic by default; any other case will require the modification of the code.
|
||||
|
@ -29,12 +25,11 @@ For decays, a specific angular distribution can be given as input as a text file
|
|||
|
||||
To run Mask simply do the following from the Mask repository:
|
||||
|
||||
`./bin/Kinematics <your_config>.yaml`
|
||||
`./bin/Kinematics input.txt`
|
||||
|
||||
`<your_config.yaml>` is a YAML configuration file. An example is given in the repository named `kinematics.yaml` and can be replaced by any yaml file with the correct format.
|
||||
Input.txt can be replaced by any text file with the correct format.
|
||||
|
||||
## Using the detector geometry simulation
|
||||
|
||||
Detector geometry is encoded using ROOT math libraries in the `src/Detectors` folder. Two different detector geometries are already present: SPS-SABRE and ANASEN. To add a new geometry, follow the guidelines outlined by each of these cases.
|
||||
|
||||
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.
|
||||
|
@ -43,12 +38,11 @@ To run the geometry code, one needs to provide an input file containing the foll
|
|||
|
||||
To run Detectors use the format
|
||||
|
||||
`./bin/Detectors <your_config>.yaml`
|
||||
`./bin/Detectors <input_file>`
|
||||
|
||||
`<your_config>.yaml` is a YAML configuration file. An example, `detector.yaml` is included in the repository.
|
||||
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.
|
||||
|
||||
Mask also provides a default visualization tool called RootPlot. RootPlot is run as
|
||||
|
@ -58,6 +52,6 @@ Mask also provides a default visualization tool called RootPlot. RootPlot is run
|
|||
where the datafile can be either the datafile from Mask or the datafile from DetEff. The outputfile is saved in the ROOT file format.
|
||||
|
||||
## Requirements
|
||||
|
||||
ROOT version 6.22 or greater is required
|
||||
CMake version 3.0 or greater is required
|
||||
|
||||
|
|
|
@ -1,5 +0,0 @@
|
|||
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
|
5
detector_input.txt
Normal file
5
detector_input.txt
Normal file
|
@ -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
|
|
@ -1,49 +0,0 @@
|
|||
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
|
45
kinematics_input.txt
Normal file
45
kinematics_input.txt
Normal file
|
@ -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
|
||||
--------------------------------------
|
|
@ -1,9 +1,5 @@
|
|||
add_executable(Detectors)
|
||||
target_include_directories(Detectors PUBLIC
|
||||
${CMAKE_CURRENT_SOURCE_DIR}
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/..
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/../vendor/yaml-cpp/include/
|
||||
)
|
||||
target_include_directories(Detectors PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..)
|
||||
|
||||
target_sources(Detectors PUBLIC
|
||||
AnasenDeadChannelMap.cpp
|
||||
|
|
|
@ -2,8 +2,6 @@
|
|||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
#include "yaml-cpp/yaml.h"
|
||||
|
||||
DetectorApp::DetectorApp() :
|
||||
m_resources(nullptr)
|
||||
{
|
||||
|
@ -18,31 +16,38 @@ DetectorApp::~DetectorApp()
|
|||
bool DetectorApp::LoadConfig(const std::string& filename)
|
||||
{
|
||||
std::cout<<"----------Detector Efficiency Calculation----------"<<std::endl;
|
||||
YAML::Node data;
|
||||
try
|
||||
std::ifstream input(filename);
|
||||
if(!input.is_open())
|
||||
{
|
||||
data = YAML::LoadFile(filename);
|
||||
}
|
||||
catch (YAML::ParserException& e)
|
||||
{
|
||||
std::cerr << "Could not load config file " << filename << " with error: " << e.what() << std::endl;
|
||||
std::cerr<<"Unable to open input config file "<<filename<<std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
m_inputFileName = data["InputDataFile"].as<std::string>();
|
||||
m_outputFileName = data["OutputDataFile"].as<std::string>();
|
||||
m_deadChannelFileName = data["DeadChannelFile"].as<std::string>();
|
||||
m_nthreads = data["NumberOfThreads"].as<uint64_t>();
|
||||
ArrayType type = StringToArrayType(data["ArrayType"].as<std::string>());
|
||||
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);
|
||||
}
|
||||
m_resources = std::make_unique<Mask::ThreadPool<DetectorArray*, uint64_t>>(m_nthreads);
|
||||
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())
|
||||
{
|
||||
|
@ -50,6 +55,7 @@ bool DetectorApp::LoadConfig(const std::string& filename)
|
|||
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;
|
||||
|
@ -58,17 +64,17 @@ bool DetectorApp::LoadConfig(const std::string& filename)
|
|||
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;
|
||||
|
||||
std::cout << "Allocating " << m_nthreads << " threads..." << std::endl;
|
||||
std::cout << "Input data file " << m_inputFileName << "..." << std::endl;
|
||||
std::cout << "With " << m_nentries << " events in the file..." << std::endl;
|
||||
std::cout << "Output data file " << m_outputFileName << "..." << std::endl;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
|
|
@ -24,6 +24,7 @@ int main(int argc, char** argv)
|
|||
return 1;
|
||||
}
|
||||
calculator.Run();
|
||||
//calculator.RunSingleThread();
|
||||
}
|
||||
catch(const std::exception& e)
|
||||
{
|
||||
|
|
|
@ -40,7 +40,7 @@ namespace Mask {
|
|||
|
||||
if(!input.is_open())
|
||||
{
|
||||
std::cerr << "Unable to open distribution file. All values reset to default." << std::endl;
|
||||
std::cerr<<"Unable to open distribution file. All values reset to default."<<std::endl;
|
||||
m_L=0;
|
||||
m_branchingRatio=1.0;
|
||||
m_constants.clear();
|
||||
|
@ -59,8 +59,8 @@ namespace Mask {
|
|||
|
||||
if(m_constants.size() != ((std::size_t) l+1))
|
||||
{
|
||||
std::cerr << "Unexpected number of constants for given angular momentum! Expected " << l+1 << " and given " << m_constants.size() << std::endl;
|
||||
std::cerr << "Setting all values to default." << std::endl;
|
||||
std::cerr<<"Unexpected number of constants for given angular momentum! Expected "<<l+1<<" and given "<<m_constants.size()<<std::endl;
|
||||
std::cerr<<"Setting all values to default."<<std::endl;
|
||||
m_L=0;
|
||||
m_branchingRatio=1.0;
|
||||
m_constants.clear();
|
||||
|
@ -73,6 +73,9 @@ namespace Mask {
|
|||
m_branchingRatio = m_constants[0]*2.0;
|
||||
m_L = l;
|
||||
|
||||
std::cout<<"Angular distribution from "<<file<<" branching ratio: "<<m_branchingRatio<<std::endl;
|
||||
std::cout<<"Angular distribution from "<<file<<" L: "<<m_L<<std::endl;
|
||||
|
||||
//Renormalize distribution such that total prob is 1.0.
|
||||
//Test branching ratio to see if we "make" a decay particle,
|
||||
//then use re-normalized distribution to pick an angle.
|
||||
|
|
|
@ -1,7 +1,6 @@
|
|||
add_library(MaskDict SHARED)
|
||||
target_include_directories(MaskDict
|
||||
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
|
||||
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../vendor/yaml-cpp/include/
|
||||
SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}
|
||||
)
|
||||
|
||||
|
@ -58,12 +57,10 @@ target_sources(Mask PRIVATE
|
|||
FileReader.cpp
|
||||
CoupledThreeStepSystem.h
|
||||
CoupledThreeStepSystem.cpp
|
||||
ConfigSerializer.h
|
||||
ConfigSerializer.cpp
|
||||
)
|
||||
|
||||
set(THREADS_PREFER_PTHREAD_FLAG On)
|
||||
find_package(Threads REQUIRED)
|
||||
target_compile_definitions(Mask PRIVATE YAML_CPP_STATIC_DEFINE)
|
||||
target_link_libraries(Mask catima yaml-cpp MaskDict ${ROOT_LIBRARIES} Threads::Threads)
|
||||
|
||||
target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES} Threads::Threads)
|
||||
set_target_properties(Mask PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})
|
|
@ -1,223 +0,0 @@
|
|||
#include "ConfigSerializer.h"
|
||||
#include "yaml-cpp/yaml.h"
|
||||
#include "RxnType.h"
|
||||
|
||||
namespace Mask {
|
||||
|
||||
static void SerializeDecay(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<params.Z.size(); i++)
|
||||
{
|
||||
yamlStream << YAML::BeginMap;
|
||||
yamlStream << YAML::Key << "Z" << YAML::Value << params.Z[i];
|
||||
yamlStream << YAML::Key << "A" << YAML::Value << params.A[i];
|
||||
yamlStream << YAML::EndMap;
|
||||
}
|
||||
yamlStream << YAML::EndSeq;
|
||||
yamlStream << YAML::Key << "PhiMin(deg)" << YAML::Value << params.phiMin;
|
||||
yamlStream << YAML::Key << "PhiMax(deg)" << YAML::Value << params.phiMax;
|
||||
yamlStream << YAML::Key << "ResidualExcitationMean(MeV)" << YAML::Value << params.meanResidualEx;
|
||||
yamlStream << YAML::Key << "ResidualExcitationSigma(MeV)" << YAML::Value << params.sigmaResidualEx;
|
||||
yamlStream << YAML::Key << "AngularDistributionFile" << YAML::Value << params.angularDistFile;
|
||||
yamlStream << YAML::EndMap;
|
||||
}
|
||||
|
||||
static StepParameters DeserializeDecay(YAML::Node& yamlStream)
|
||||
{
|
||||
StepParameters params;
|
||||
params.rxnType = StringToRxnType(yamlStream["Type"].as<std::string>());
|
||||
YAML::Node nuclei = yamlStream["Reactants"];
|
||||
for(auto nucleus : nuclei)
|
||||
{
|
||||
params.Z.push_back(nucleus["Z"].as<int>());
|
||||
params.A.push_back(nucleus["A"].as<int>());
|
||||
}
|
||||
params.phiMin = yamlStream["PhiMin(deg)"].as<double>();
|
||||
params.phiMax = yamlStream["PhiMax(deg)"].as<double>();
|
||||
params.meanResidualEx = yamlStream["ResidualExcitationMean(MeV)"].as<double>();
|
||||
params.sigmaResidualEx = yamlStream["ResidualExcitationSigma(MeV)"].as<double>();
|
||||
params.angularDistFile = yamlStream["AngularDistributionFile"].as<std::string>();
|
||||
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<params.Z.size(); i++)
|
||||
{
|
||||
yamlStream << YAML::BeginMap;
|
||||
yamlStream << YAML::Key << "Z" << YAML::Value << params.Z[i];
|
||||
yamlStream << YAML::Key << "A" << YAML::Value << params.A[i];
|
||||
yamlStream << YAML::EndMap;
|
||||
}
|
||||
yamlStream << YAML::EndSeq;
|
||||
yamlStream << YAML::Key << "BeamEnergyMean(MeV)" << YAML::Value << params.meanBeamEnergy;
|
||||
yamlStream << YAML::Key << "BeamEnergySigma(MeV)" << YAML::Value << params.sigmaBeamEnergy;
|
||||
yamlStream << YAML::Key << "ThetaType" << YAML::Value << RxnThetaTypeToString(params.thetaType);
|
||||
yamlStream << YAML::Key << "ThetaMin(deg)" << YAML::Value << params.thetaMin;
|
||||
yamlStream << YAML::Key << "ThetaMax(deg)" << YAML::Value << params.thetaMax;
|
||||
yamlStream << YAML::Key << "PhiMin(deg)" << YAML::Value << params.phiMin;
|
||||
yamlStream << YAML::Key << "PhiMax(deg)" << YAML::Value << params.phiMax;
|
||||
yamlStream << YAML::Key << "ResidualExcitationMean(MeV)" << YAML::Value << params.meanResidualEx;
|
||||
yamlStream << YAML::Key << "ResidualExcitationSigma(MeV)" << YAML::Value << params.sigmaResidualEx;
|
||||
yamlStream << YAML::EndMap;
|
||||
}
|
||||
|
||||
static StepParameters DeserializeReaction(YAML::Node& yamlStream)
|
||||
{
|
||||
StepParameters params;
|
||||
params.rxnType = StringToRxnType(yamlStream["Type"].as<std::string>());
|
||||
YAML::Node nuclei = yamlStream["Reactants"];
|
||||
for(auto nucleus : nuclei)
|
||||
{
|
||||
params.Z.push_back(nucleus["Z"].as<int>());
|
||||
params.A.push_back(nucleus["A"].as<int>());
|
||||
}
|
||||
params.meanBeamEnergy = yamlStream["BeamEnergyMean(MeV)"].as<double>();
|
||||
params.sigmaBeamEnergy = yamlStream["BeamEnergySigma(MeV)"].as<double>();
|
||||
params.thetaType = StringToRxnThetaType(yamlStream["ThetaType"].as<std::string>());
|
||||
params.thetaMin = yamlStream["ThetaMin(deg)"].as<double>();
|
||||
params.thetaMax = yamlStream["ThetaMax(deg)"].as<double>();
|
||||
params.phiMin = yamlStream["PhiMin(deg)"].as<double>();
|
||||
params.phiMax = yamlStream["PhiMax(deg)"].as<double>();
|
||||
params.meanResidualEx = yamlStream["ResidualExcitationMean(MeV)"].as<double>();
|
||||
params.sigmaResidualEx = yamlStream["ResidualExcitationSigma(MeV)"].as<double>();
|
||||
return params;
|
||||
}
|
||||
|
||||
static void SerializeTarget(YAML::Emitter& yamlStream, const LayeredTarget& target)
|
||||
{
|
||||
|
||||
yamlStream << YAML::BeginSeq;
|
||||
for (int i=0; i<target.GetNumberOfLayers(); i++)
|
||||
{
|
||||
const Target& layer = target.GetLayerInfo(i);
|
||||
yamlStream << YAML::BeginMap;
|
||||
yamlStream << YAML::Key << "Thickness(ug/cm^2)" << YAML::Value << layer.GetThickness();
|
||||
yamlStream << YAML::Key << "Elements" << YAML::Value << YAML::BeginSeq;
|
||||
for (int j=0; j<layer.GetNumberOfElements(); j++)
|
||||
{
|
||||
yamlStream << YAML::BeginMap;
|
||||
yamlStream << YAML::Key << "Z" << YAML::Value << layer.GetElementZ(j);
|
||||
yamlStream << YAML::Key << "A" << YAML::Value << layer.GetElementA(j);
|
||||
yamlStream << YAML::Key << "S" << YAML::Value << layer.GetElementStoich(j);
|
||||
yamlStream << YAML::EndMap;
|
||||
}
|
||||
yamlStream << YAML::EndSeq;
|
||||
yamlStream << YAML::EndMap;
|
||||
}
|
||||
yamlStream << YAML::EndSeq;
|
||||
}
|
||||
|
||||
static LayeredTarget DeserializeTarget(YAML::Node& yamlStream)
|
||||
{
|
||||
LayeredTarget target;
|
||||
std::vector<int> z, a, s;
|
||||
double thickness;
|
||||
|
||||
for (auto layer : yamlStream)
|
||||
{
|
||||
thickness = layer["Thickness(ug/cm^2)"].as<double>();
|
||||
YAML::Node elements = layer["Elements"];
|
||||
z.clear();
|
||||
a.clear();
|
||||
s.clear();
|
||||
for(auto element : elements)
|
||||
{
|
||||
z.push_back(element["Z"].as<int>());
|
||||
a.push_back(element["A"].as<int>());
|
||||
s.push_back(element["S"].as<int>());
|
||||
}
|
||||
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() <<std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
params.outputFileName = data["OutputFile"].as<std::string>();
|
||||
params.nThreads = data["Threads"].as<uint32_t>();
|
||||
params.nSamples = data["ReactionSamples"].as<uint64_t>();
|
||||
|
||||
auto steps = data["ReactionChain"];
|
||||
for (auto step : steps)
|
||||
{
|
||||
RxnType type = StringToRxnType(step["Type"].as<std::string>());
|
||||
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;
|
||||
}
|
||||
}
|
|
@ -1,16 +0,0 @@
|
|||
#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
|
|
@ -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() const { return m_layers.size(); }
|
||||
std::size_t GetNumberOfLayers() { return m_layers.size(); }
|
||||
void SetName(std::string& n) { m_name = n; }
|
||||
const Target& GetLayerInfo(int index) const { return m_layers[index]; }
|
||||
const Target& GetLayerInfo(int index) { return m_layers[index]; }
|
||||
const std::string& GetName() { return m_name; }
|
||||
|
||||
private:
|
||||
|
|
|
@ -1,80 +1,185 @@
|
|||
#include "MaskApp.h"
|
||||
#include "ConfigSerializer.h"
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
#include "TFile.h"
|
||||
#include "TTree.h"
|
||||
|
||||
#include "yaml-cpp/yaml.h"
|
||||
|
||||
namespace Mask {
|
||||
|
||||
MaskApp::MaskApp() :
|
||||
m_resources(nullptr)
|
||||
m_system(nullptr), m_resources(nullptr)
|
||||
{
|
||||
std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
|
||||
}
|
||||
|
||||
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)
|
||||
|
||||
bool MaskApp::LoadConfig(const std::string& filename)
|
||||
{
|
||||
std::cout << "Loading configuration in " << filename << std::endl;
|
||||
|
||||
if (!ConfigSerializer::DeserializeConfig(filename, m_params))
|
||||
std::cout<<"Loading configuration in "<<filename<<"..."<<std::endl;
|
||||
|
||||
std::ifstream input(filename);
|
||||
if(!input.is_open())
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
//Initialize the system
|
||||
//Reaction chain
|
||||
for(uint32_t i=0; i<m_params.nThreads; i++)
|
||||
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)
|
||||
{
|
||||
m_systemList.push_back(CreateSystem(m_params.chainParams));
|
||||
if(junk == "begin_chain")
|
||||
continue;
|
||||
else if (junk == "end_chain")
|
||||
break;
|
||||
else if(junk == "begin_step")
|
||||
{
|
||||
StepParameters currentParams;
|
||||
input >> 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."<<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;
|
||||
}
|
||||
}
|
||||
//Link the target
|
||||
|
||||
std::getline(input, junk);
|
||||
std::getline(input, junk);
|
||||
|
||||
LayeredTarget target;
|
||||
int nlayers;
|
||||
double thickness;
|
||||
std::vector<int> avec, zvec, svec;
|
||||
int s;
|
||||
input>>junk>>nlayers;
|
||||
for(int i=0; i<nlayers; i++)
|
||||
{
|
||||
input>>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);
|
||||
|
||||
for(auto system : m_systemList)
|
||||
{
|
||||
system->SetLayeredTarget(m_params.target);
|
||||
system->SetLayeredTarget(target);
|
||||
}
|
||||
//Setup threading
|
||||
m_resources = std::make_unique<ThreadPool<ReactionSystem*, uint64_t>>(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; i<m_params.nThreads; i++)
|
||||
m_chunkSamples.push_back(quotient);
|
||||
m_fileWriter.Open(m_params.outputFileName, "SimTree");
|
||||
|
||||
|
||||
std::cout << "Reaction equation: " << m_systemList[0]->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;
|
||||
std::cout<<"Reaction equation: "<<m_system->GetSystemEquation()<<std::endl;
|
||||
std::cout<<"Number of samples: "<<m_nsamples<<std::endl;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
//Not implemented... yet!
|
||||
bool MaskApp::SaveConfig(const std::string& filename)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
bool MaskApp::SaveConfig(const std::string& filename)
|
||||
{
|
||||
std::cout << "Writing configuration to " << filename << std::endl;
|
||||
return ConfigSerializer::SerializeConfig(filename, m_params);
|
||||
}
|
||||
|
||||
void MaskApp::Run()
|
||||
{
|
||||
std::cout<<"Running simulation..."<<std::endl;
|
||||
if(m_systemList.size() != m_params.nThreads)
|
||||
if(m_systemList.size() != m_nthreads)
|
||||
{
|
||||
std::cerr << "System list not equal to number of threads" << std::endl;
|
||||
return;
|
||||
|
@ -100,7 +205,7 @@ namespace Mask {
|
|||
|
||||
uint64_t count = 0;
|
||||
double percent = 0.05;
|
||||
uint64_t flushVal = m_params.nSamples*percent;
|
||||
uint64_t flushVal = m_nsamples*percent;
|
||||
uint64_t flushCount = 0;
|
||||
while(true)
|
||||
{
|
||||
|
@ -122,4 +227,42 @@ namespace Mask {
|
|||
std::cout<<"---------------------------------------------"<<std::endl;
|
||||
}
|
||||
|
||||
void MaskApp::RunSingleThread()
|
||||
{
|
||||
std::cout<<"Running simulation..."<<std::endl;
|
||||
if(m_system == nullptr)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
TFile* output = TFile::Open(m_outputName.c_str(), "RECREATE");
|
||||
TTree* tree = new TTree("SimTree", "SimTree");
|
||||
tree->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; i<m_nsamples; i++)
|
||||
{
|
||||
if(++count == percent5)
|
||||
{
|
||||
npercent++;
|
||||
count = 0;
|
||||
std::cout<<"\rPercent complete: "<<npercent*5<<"%"<<std::flush;
|
||||
}
|
||||
|
||||
m_system->RunSystem();
|
||||
tree->Fill();
|
||||
}
|
||||
|
||||
tree->Write(tree->GetName(), TObject::kOverwrite);
|
||||
output->Close();
|
||||
|
||||
std::cout<<std::endl;
|
||||
std::cout<<"Complete."<<std::endl;
|
||||
std::cout<<"---------------------------------------------"<<std::endl;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
@ -14,15 +14,6 @@
|
|||
|
||||
namespace Mask {
|
||||
|
||||
struct AppParameters
|
||||
{
|
||||
std::string outputFileName = "";
|
||||
uint64_t nSamples = 0;
|
||||
uint32_t nThreads = 1;
|
||||
std::vector<StepParameters> chainParams;
|
||||
LayeredTarget target;
|
||||
};
|
||||
|
||||
class MaskApp
|
||||
{
|
||||
public:
|
||||
|
@ -32,9 +23,17 @@ namespace Mask {
|
|||
bool SaveConfig(const std::string& filename);
|
||||
|
||||
void Run();
|
||||
void RunSingleThread();
|
||||
|
||||
private:
|
||||
AppParameters m_params;
|
||||
void RunChunk(ReactionSystem* system);
|
||||
|
||||
ReactionSystem* m_system;
|
||||
|
||||
std::string m_outputName;
|
||||
RxnType m_rxnType;
|
||||
uint64_t m_nsamples;
|
||||
uint32_t m_nthreads;
|
||||
|
||||
std::vector<ReactionSystem*> m_systemList; //One system for each thread
|
||||
std::vector<uint64_t> m_chunkSamples;
|
||||
|
|
|
@ -37,17 +37,6 @@ 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")
|
||||
|
@ -63,17 +52,6 @@ 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
|
|
@ -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);
|
||||
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]; }
|
||||
|
||||
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]; }
|
||||
|
||||
private:
|
||||
void Init(const std::vector<int>& z, const std::vector<int>& a, const std::vector<int>& stoich);
|
||||
|
||||
|
|
1
src/vendor/yaml-cpp
vendored
1
src/vendor/yaml-cpp
vendored
|
@ -1 +0,0 @@
|
|||
Subproject commit 0e6e28d1a38224fc8172fae0109ea7f673c096db
|
Loading…
Reference in New Issue
Block a user