1
0
Fork 0
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.

21 changed files with 288 additions and 421 deletions

1
.gitignore vendored
View File

@ -19,7 +19,6 @@ build/
lib/
bin/
scripts/
*.yaml
###Keep this file###

3
.gitmodules vendored
View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -24,6 +24,7 @@ int main(int argc, char** argv)
return 1;
}
calculator.Run();
//calculator.RunSingleThread();
}
catch(const std::exception& e)
{

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -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 +0,0 @@
Subproject commit 0e6e28d1a38224fc8172fae0109ea7f673c096db