mirror of
https://github.com/gwm17/Mask.git
synced 2024-11-22 18:28:51 -05:00
Compare commits
8 Commits
a4eaeaa634
...
9628e73328
Author | SHA1 | Date | |
---|---|---|---|
9628e73328 | |||
Gordon McCann | 142a2dffd1 | ||
Gordon McCann | 47d84e8dca | ||
Gordon McCann | 49b54af286 | ||
Gordon McCann | d68f5de383 | ||
Gordon McCann | cf22fb6384 | ||
Gordon McCann | 6815d3d2e7 | ||
Gordon McCann | 29a136a14b |
2
.gitignore
vendored
2
.gitignore
vendored
|
@ -18,6 +18,8 @@ Makefile
|
||||||
build/
|
build/
|
||||||
lib/
|
lib/
|
||||||
bin/
|
bin/
|
||||||
|
scripts/
|
||||||
|
*.yaml
|
||||||
|
|
||||||
|
|
||||||
###Keep this file###
|
###Keep this file###
|
||||||
|
|
3
.gitmodules
vendored
3
.gitmodules
vendored
|
@ -1,3 +1,6 @@
|
||||||
[submodule "src/vendor/catima"]
|
[submodule "src/vendor/catima"]
|
||||||
path = src/vendor/catima
|
path = src/vendor/catima
|
||||||
url = https://github.com/gwm17/catima.git
|
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,6 +19,8 @@ set(CMAKE_CXX_STANDARD 17)
|
||||||
find_package(ROOT REQUIRED COMPONENTS GenVector)
|
find_package(ROOT REQUIRED COMPONENTS GenVector)
|
||||||
|
|
||||||
add_subdirectory(src/vendor/catima)
|
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/Mask)
|
||||||
add_subdirectory(src/Kinematics)
|
add_subdirectory(src/Kinematics)
|
||||||
add_subdirectory(src/Detectors)
|
add_subdirectory(src/Detectors)
|
||||||
|
|
18
README.md
18
README.md
|
@ -1,9 +1,12 @@
|
||||||
# Mask: Monte cArlo Simulation of Kinematics
|
# 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 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.
|
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
|
## Building Mask
|
||||||
Dowload the repository from github. CMake is use to build the project; in most environments you can build Mask using the following methods:
|
|
||||||
|
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:
|
||||||
|
|
||||||
- `mkdir build`
|
- `mkdir build`
|
||||||
- `cd build`
|
- `cd build`
|
||||||
- `cmake ..`
|
- `cmake ..`
|
||||||
|
@ -14,6 +17,7 @@ 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.
|
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
|
## 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:
|
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.
|
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.
|
||||||
|
@ -25,11 +29,12 @@ 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:
|
To run Mask simply do the following from the Mask repository:
|
||||||
|
|
||||||
`./bin/Kinematics input.txt`
|
`./bin/Kinematics <your_config>.yaml`
|
||||||
|
|
||||||
Input.txt can be replaced by any text file with the correct format.
|
`<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.
|
||||||
|
|
||||||
## Using the detector geometry simulation
|
## 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.
|
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.
|
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.
|
||||||
|
@ -38,11 +43,12 @@ To run the geometry code, one needs to provide an input file containing the foll
|
||||||
|
|
||||||
To run Detectors use the format
|
To run Detectors use the format
|
||||||
|
|
||||||
`./bin/Detectors <input_file>`
|
`./bin/Detectors <your_config>.yaml`
|
||||||
|
|
||||||
An example input file is provided with the repository.
|
`<your_config>.yaml` is a YAML configuration file. An example, `detector.yaml` is included in the repository.
|
||||||
|
|
||||||
## Data visualization
|
## 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.
|
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
|
Mask also provides a default visualization tool called RootPlot. RootPlot is run as
|
||||||
|
@ -52,6 +58,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.
|
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
|
## Requirements
|
||||||
|
|
||||||
ROOT version 6.22 or greater is required
|
ROOT version 6.22 or greater is required
|
||||||
CMake version 3.0 or greater is required
|
CMake version 3.0 or greater is required
|
||||||
|
|
||||||
|
|
5
detector.yaml
Normal file
5
detector.yaml
Normal file
|
@ -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
|
|
@ -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
|
|
49
kinematics.yaml
Normal file
49
kinematics.yaml
Normal file
|
@ -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
|
|
@ -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
|
|
||||||
--------------------------------------
|
|
|
@ -1,5 +1,9 @@
|
||||||
add_executable(Detectors)
|
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
|
target_sources(Detectors PUBLIC
|
||||||
AnasenDeadChannelMap.cpp
|
AnasenDeadChannelMap.cpp
|
||||||
|
|
|
@ -2,6 +2,8 @@
|
||||||
#include <fstream>
|
#include <fstream>
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
|
|
||||||
|
#include "yaml-cpp/yaml.h"
|
||||||
|
|
||||||
DetectorApp::DetectorApp() :
|
DetectorApp::DetectorApp() :
|
||||||
m_resources(nullptr)
|
m_resources(nullptr)
|
||||||
{
|
{
|
||||||
|
@ -16,38 +18,31 @@ DetectorApp::~DetectorApp()
|
||||||
bool DetectorApp::LoadConfig(const std::string& filename)
|
bool DetectorApp::LoadConfig(const std::string& filename)
|
||||||
{
|
{
|
||||||
std::cout<<"----------Detector Efficiency Calculation----------"<<std::endl;
|
std::cout<<"----------Detector Efficiency Calculation----------"<<std::endl;
|
||||||
std::ifstream input(filename);
|
YAML::Node data;
|
||||||
if(!input.is_open())
|
try
|
||||||
{
|
{
|
||||||
std::cerr<<"Unable to open input config file "<<filename<<std::endl;
|
data = YAML::LoadFile(filename);
|
||||||
|
}
|
||||||
|
catch (YAML::ParserException& e)
|
||||||
|
{
|
||||||
|
std::cerr << "Could not load config file " << filename << " with error: " << e.what() << std::endl;
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
|
|
||||||
std::string junk;
|
m_inputFileName = data["InputDataFile"].as<std::string>();
|
||||||
std::string typeName;
|
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>());
|
||||||
|
|
||||||
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++)
|
for(uint64_t i=0; i<m_nthreads; i++)
|
||||||
{
|
{
|
||||||
m_detectorList.push_back(CreateDetectorArray(type));
|
m_detectorList.push_back(CreateDetectorArray(type));
|
||||||
if(m_deadChannelFileName != "None")
|
if(m_deadChannelFileName != "None")
|
||||||
m_detectorList.back()->SetDeadChannelMap(m_deadChannelFileName);
|
m_detectorList.back()->SetDeadChannelMap(m_deadChannelFileName);
|
||||||
}
|
}
|
||||||
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);
|
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");
|
m_fileReader.Open(m_inputFileName, "SimTree");
|
||||||
if(!m_fileReader.IsOpen() || !m_fileReader.IsTree())
|
if(!m_fileReader.IsOpen() || !m_fileReader.IsTree())
|
||||||
{
|
{
|
||||||
|
@ -55,7 +50,6 @@ bool DetectorApp::LoadConfig(const std::string& filename)
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
m_nentries = m_fileReader.GetSize();
|
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
|
//Little bit of integer division mangling to make sure we read every event in file
|
||||||
uint64_t quotient = m_nentries / m_nthreads;
|
uint64_t quotient = m_nentries / m_nthreads;
|
||||||
|
@ -64,17 +58,17 @@ bool DetectorApp::LoadConfig(const std::string& filename)
|
||||||
for(uint64_t i=1; i<m_nthreads; i++)
|
for(uint64_t i=1; i<m_nthreads; i++)
|
||||||
m_chunkSamples.push_back(quotient);
|
m_chunkSamples.push_back(quotient);
|
||||||
|
|
||||||
std::cout << "Opening output data file " << m_outputFileName << std::endl;
|
|
||||||
m_fileWriter.Open(m_outputFileName, "SimTree");
|
m_fileWriter.Open(m_outputFileName, "SimTree");
|
||||||
if(!m_fileWriter.IsOpen() || !m_fileWriter.IsTree())
|
if(!m_fileWriter.IsOpen() || !m_fileWriter.IsTree())
|
||||||
{
|
{
|
||||||
std::cerr << "Unable to open output data file " << m_outputFileName << std::endl;
|
std::cerr << "Unable to open output data file " << m_outputFileName << std::endl;
|
||||||
return false;
|
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;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -110,6 +104,7 @@ void DetectorApp::Run()
|
||||||
nucleus.detectedKE = result.energy_deposited;
|
nucleus.detectedKE = result.energy_deposited;
|
||||||
nucleus.detectedTheta = result.direction.Theta();
|
nucleus.detectedTheta = result.direction.Theta();
|
||||||
nucleus.detectedPhi = result.direction.Phi();
|
nucleus.detectedPhi = result.direction.Phi();
|
||||||
|
nucleus.detectedPos = result.direction;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
m_fileWriter.PushData(data);
|
m_fileWriter.PushData(data);
|
||||||
|
@ -141,4 +136,6 @@ void DetectorApp::Run()
|
||||||
++count;
|
++count;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
std::cout << std::endl;
|
||||||
|
|
||||||
}
|
}
|
|
@ -20,6 +20,12 @@ SabreArray::SabreArray() :
|
||||||
m_degradedDetectors[2] = false;
|
m_degradedDetectors[2] = false;
|
||||||
m_degradedDetectors[3] = false;
|
m_degradedDetectors[3] = false;
|
||||||
m_degradedDetectors[4] = true;
|
m_degradedDetectors[4] = true;
|
||||||
|
//No degraded detectors
|
||||||
|
// m_degradedDetectors[0] = false;
|
||||||
|
// m_degradedDetectors[1] = false;
|
||||||
|
// m_degradedDetectors[2] = false;
|
||||||
|
// m_degradedDetectors[3] = false;
|
||||||
|
// m_degradedDetectors[4] = false;
|
||||||
|
|
||||||
//Choose who to look at right now. Usually switch on or off degraded/non-degraded.
|
//Choose who to look at right now. Usually switch on or off degraded/non-degraded.
|
||||||
m_activeDetectors[0] = false;
|
m_activeDetectors[0] = false;
|
||||||
|
|
|
@ -39,7 +39,7 @@ private:
|
||||||
static constexpr double s_detectorThickness = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density)
|
static constexpr double s_detectorThickness = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density)
|
||||||
static constexpr double s_degraderThickness = 70.0 * 1.0e-4 * 16.69 * 1e6; //tantalum degrader (70 um thick)
|
static constexpr double s_degraderThickness = 70.0 * 1.0e-4 * 16.69 * 1e6; //tantalum degrader (70 um thick)
|
||||||
|
|
||||||
static constexpr double s_energyThreshold = 0.2; //in MeV
|
static constexpr double s_energyThreshold = 0.25; //in MeV
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
|
@ -24,7 +24,6 @@ int main(int argc, char** argv)
|
||||||
return 1;
|
return 1;
|
||||||
}
|
}
|
||||||
calculator.Run();
|
calculator.Run();
|
||||||
//calculator.RunSingleThread();
|
|
||||||
}
|
}
|
||||||
catch(const std::exception& e)
|
catch(const std::exception& e)
|
||||||
{
|
{
|
||||||
|
|
|
@ -73,13 +73,10 @@ namespace Mask {
|
||||||
m_branchingRatio = m_constants[0]*2.0;
|
m_branchingRatio = m_constants[0]*2.0;
|
||||||
m_L = l;
|
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.
|
//Renormalize distribution such that total prob is 1.0.
|
||||||
//Test branching ratio to see if we "make" a decay particle,
|
//Test branching ratio to see if we "make" a decay particle,
|
||||||
//then use re-normalized distribution to pick an angle.
|
//then use re-normalized distribution to pick an angle.
|
||||||
if(m_constants[0] < 0.5)
|
if(m_constants[0] != 0.5)
|
||||||
{
|
{
|
||||||
double norm = 0.5/m_constants[0];
|
double norm = 0.5/m_constants[0];
|
||||||
for(auto& value : m_constants)
|
for(auto& value : m_constants)
|
||||||
|
@ -115,4 +112,14 @@ namespace Mask {
|
||||||
return costheta;
|
return costheta;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
double AngularDistribution::GetProbability(double cosTheta)
|
||||||
|
{
|
||||||
|
double prob = 0.0;
|
||||||
|
for (std::size_t i=0; i<m_constants.size(); i++)
|
||||||
|
{
|
||||||
|
prob += m_constants[i] * P_l(i*2, cosTheta);
|
||||||
|
}
|
||||||
|
return prob;
|
||||||
|
}
|
||||||
|
|
||||||
}
|
}
|
|
@ -17,6 +17,7 @@ namespace Mask {
|
||||||
double GetRandomCosTheta();
|
double GetRandomCosTheta();
|
||||||
int GetL() { return m_L; }
|
int GetL() { return m_L; }
|
||||||
double GetBranchingRatio() { return m_branchingRatio; }
|
double GetBranchingRatio() { return m_branchingRatio; }
|
||||||
|
double GetProbability(double cosTheta);
|
||||||
|
|
||||||
private:
|
private:
|
||||||
bool IsIsotropic() { return m_isIsotropic; }
|
bool IsIsotropic() { return m_isIsotropic; }
|
||||||
|
|
|
@ -1,6 +1,7 @@
|
||||||
add_library(MaskDict SHARED)
|
add_library(MaskDict SHARED)
|
||||||
target_include_directories(MaskDict
|
target_include_directories(MaskDict
|
||||||
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
|
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
|
||||||
|
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../vendor/yaml-cpp/include/
|
||||||
SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}
|
SYSTEM PUBLIC ${ROOT_INCLUDE_DIRS}
|
||||||
)
|
)
|
||||||
|
|
||||||
|
@ -55,10 +56,14 @@ target_sources(Mask PRIVATE
|
||||||
FileWriter.cpp
|
FileWriter.cpp
|
||||||
FileReader.h
|
FileReader.h
|
||||||
FileReader.cpp
|
FileReader.cpp
|
||||||
|
CoupledThreeStepSystem.h
|
||||||
|
CoupledThreeStepSystem.cpp
|
||||||
|
ConfigSerializer.h
|
||||||
|
ConfigSerializer.cpp
|
||||||
)
|
)
|
||||||
|
|
||||||
set(THREADS_PREFER_PTHREAD_FLAG On)
|
set(THREADS_PREFER_PTHREAD_FLAG On)
|
||||||
find_package(Threads REQUIRED)
|
find_package(Threads REQUIRED)
|
||||||
|
target_compile_definitions(Mask PRIVATE YAML_CPP_STATIC_DEFINE)
|
||||||
target_link_libraries(Mask catima MaskDict ${ROOT_LIBRARIES} Threads::Threads)
|
target_link_libraries(Mask catima yaml-cpp MaskDict ${ROOT_LIBRARIES} Threads::Threads)
|
||||||
set_target_properties(Mask PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})
|
set_target_properties(Mask PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})
|
223
src/Mask/ConfigSerializer.cpp
Normal file
223
src/Mask/ConfigSerializer.cpp
Normal file
|
@ -0,0 +1,223 @@
|
||||||
|
#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;
|
||||||
|
}
|
||||||
|
}
|
16
src/Mask/ConfigSerializer.h
Normal file
16
src/Mask/ConfigSerializer.h
Normal file
|
@ -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
|
210
src/Mask/CoupledThreeStepSystem.cpp
Normal file
210
src/Mask/CoupledThreeStepSystem.cpp
Normal file
|
@ -0,0 +1,210 @@
|
||||||
|
#include "CoupledThreeStepSystem.h"
|
||||||
|
|
||||||
|
#include "Math/Boost.h"
|
||||||
|
#include "Math/Vector3D.h"
|
||||||
|
#include "Math/RotationY.h"
|
||||||
|
#include "Math/RotationZ.h"
|
||||||
|
|
||||||
|
namespace Mask {
|
||||||
|
|
||||||
|
CoupledThreeStepSystem::CoupledThreeStepSystem(const std::vector<StepParameters>& params) :
|
||||||
|
ReactionSystem()
|
||||||
|
{
|
||||||
|
m_nuclei.resize(8);
|
||||||
|
Init(params);
|
||||||
|
}
|
||||||
|
|
||||||
|
CoupledThreeStepSystem::~CoupledThreeStepSystem() {}
|
||||||
|
|
||||||
|
void CoupledThreeStepSystem::Init(const std::vector<StepParameters>& params)
|
||||||
|
{
|
||||||
|
if(params.size() != 3 || params[0].rxnType != RxnType::Reaction || params[1].rxnType != RxnType::Decay ||
|
||||||
|
params[2].rxnType != RxnType::Decay || params[0].Z.size() != 3 || params[0].A.size() != 3 || params[1].Z.size() != 2 ||
|
||||||
|
params[1].A.size() != 2 || params[2].Z.size() != 2 || params[2].A.size() != 2)
|
||||||
|
{
|
||||||
|
m_isValid = false;
|
||||||
|
std::cerr << "Invalid parameters at ThreeStepSystem::Init(), does not match ThreeStep signature!" << std::endl;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
const StepParameters& step1Params = params[0];
|
||||||
|
const StepParameters& step2Params = params[1];
|
||||||
|
const StepParameters& step3Params = params[2];
|
||||||
|
|
||||||
|
//Setup nuclei
|
||||||
|
int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2];
|
||||||
|
int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2];
|
||||||
|
if(zr != step2Params.Z[0] || ar != step2Params.A[0])
|
||||||
|
{
|
||||||
|
m_isValid = false;
|
||||||
|
std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step one and step two are not sequential! Step one recoil (Z,A): ("
|
||||||
|
<< zr << "," << ar << ") Step two target (Z,A): (" << step2Params.Z[0] << "," << step2Params.A[0] << ")" <<std::endl;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
int zb2 = step2Params.Z[0] - step2Params.Z[1];
|
||||||
|
int ab2 = step2Params.A[0] - step2Params.A[1];
|
||||||
|
if(zb2 != step3Params.Z[0] || ab2 != step3Params.A[0])
|
||||||
|
{
|
||||||
|
m_isValid = false;
|
||||||
|
std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step two and step three are not sequential! Step two heavy (Z,A): ("
|
||||||
|
<< zb2 << "," << ab2 << ") Step three target (Z,A): (" << step3Params.Z[0] << "," << step3Params.A[0] << ")" <<std::endl;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
int zb4 = step3Params.Z[0] - step3Params.Z[1];
|
||||||
|
int ab4 = step3Params.A[0] - step3Params.A[1];
|
||||||
|
|
||||||
|
m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target
|
||||||
|
m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //projectile
|
||||||
|
m_nuclei[2] = CreateNucleus(step1Params.Z[2], step1Params.A[2]); //ejectile
|
||||||
|
m_nuclei[3] = CreateNucleus(zr, ar); //residual
|
||||||
|
m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
|
||||||
|
m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
|
||||||
|
m_nuclei[6] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
|
||||||
|
m_nuclei[7] = CreateNucleus(zb4, ab4); //breakup4
|
||||||
|
|
||||||
|
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
|
||||||
|
m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
|
||||||
|
m_step3.BindNuclei(&(m_nuclei[5]), nullptr, &(m_nuclei[6]), &(m_nuclei[7]));
|
||||||
|
|
||||||
|
SetSystemEquation();
|
||||||
|
|
||||||
|
//Step one sampling parameters
|
||||||
|
AddBeamDistribution(step1Params.meanBeamEnergy, step1Params.sigmaBeamEnergy);
|
||||||
|
m_step1.SetEjectileThetaType(step1Params.thetaType);
|
||||||
|
AddThetaRange(step1Params.thetaMin, step1Params.thetaMax);
|
||||||
|
AddPhiRange(step1Params.phiMin, step1Params.phiMax);
|
||||||
|
AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx);
|
||||||
|
|
||||||
|
//Step two sampling parameters
|
||||||
|
AddPhiRange(step2Params.phiMin, step2Params.phiMax);
|
||||||
|
AddDecayAngularDistribution(step2Params.angularDistFile);
|
||||||
|
AddExcitationDistribution(step2Params.meanResidualEx, step2Params.sigmaResidualEx);
|
||||||
|
|
||||||
|
//Step three sampling parameters
|
||||||
|
AddPhiRange(step3Params.phiMin, step3Params.phiMax);
|
||||||
|
AddDecayAngularDistribution(step3Params.angularDistFile);
|
||||||
|
AddExcitationDistribution(step3Params.meanResidualEx, step3Params.sigmaResidualEx);
|
||||||
|
}
|
||||||
|
|
||||||
|
void CoupledThreeStepSystem::SetLayeredTarget(const LayeredTarget& target)
|
||||||
|
{
|
||||||
|
m_target = target;
|
||||||
|
m_rxnLayer = m_target.FindLayerContaining(m_nuclei[0].Z, m_nuclei[0].A);
|
||||||
|
if(m_rxnLayer != m_target.GetNumberOfLayers())
|
||||||
|
{
|
||||||
|
m_step1.SetLayeredTarget(&m_target);
|
||||||
|
m_step2.SetLayeredTarget(&m_target);
|
||||||
|
m_step3.SetLayeredTarget(&m_target);
|
||||||
|
m_step1.SetRxnLayer(m_rxnLayer);
|
||||||
|
m_step2.SetRxnLayer(m_rxnLayer);
|
||||||
|
m_step3.SetRxnLayer(m_rxnLayer);
|
||||||
|
m_isTargetSet = true;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
throw ReactionLayerException();
|
||||||
|
}
|
||||||
|
|
||||||
|
void CoupledThreeStepSystem::SetSystemEquation()
|
||||||
|
{
|
||||||
|
std::stringstream stream;
|
||||||
|
stream << m_nuclei[0].isotopicSymbol << "("
|
||||||
|
<< m_nuclei[1].isotopicSymbol << ", "
|
||||||
|
<< m_nuclei[2].isotopicSymbol << ")"
|
||||||
|
<< m_nuclei[3].isotopicSymbol << "->"
|
||||||
|
<< m_nuclei[4].isotopicSymbol << "+"
|
||||||
|
<< m_nuclei[5].isotopicSymbol << "->"
|
||||||
|
<< m_nuclei[6].isotopicSymbol << "+"
|
||||||
|
<< m_nuclei[7].isotopicSymbol;
|
||||||
|
m_sysEquation = stream.str();
|
||||||
|
}
|
||||||
|
|
||||||
|
CoupledThreeStepParameters CoupledThreeStepSystem::SampleParameters()
|
||||||
|
{
|
||||||
|
CoupledThreeStepParameters params;
|
||||||
|
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
|
||||||
|
params.beamEnergy = m_beamDistributions[0](gen);
|
||||||
|
params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
|
||||||
|
params.rxnPhi = m_phiRanges[0](gen);
|
||||||
|
params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
|
||||||
|
params.cosRelativeAngle = m_decayAngularDistributions[1].GetRandomCosTheta();
|
||||||
|
params.decay1Theta = std::acos(params.cosdecay1Theta);
|
||||||
|
params.decay1Phi = m_phiRanges[1](gen);
|
||||||
|
params.residEx = m_exDistributions[0](gen);
|
||||||
|
params.decay1Ex = m_exDistributions[1](gen);
|
||||||
|
params.decay2Ex = m_exDistributions[2](gen);
|
||||||
|
params.rxnDepth = (m_rxnDepthDist(gen));
|
||||||
|
return params;
|
||||||
|
}
|
||||||
|
|
||||||
|
//Called after running step2
|
||||||
|
void CoupledThreeStepSystem::SampleCoupling(CoupledThreeStepParameters& params)
|
||||||
|
{
|
||||||
|
ROOT::Math::Boost boost(m_nuclei[5].vec4.BoostToCM());
|
||||||
|
ROOT::Math::PxPyPzEVector a1Vec = (boost * m_nuclei[4].vec4);
|
||||||
|
ROOT::Math::XYZVector a1PVec(a1Vec.Px(), a1Vec.Py(), a1Vec.Pz());
|
||||||
|
ROOT::Math::PxPyPzEVector a2Vec;
|
||||||
|
ROOT::Math::XYZVector a2PVec;
|
||||||
|
ROOT::Math::RotationZ zRot;
|
||||||
|
ROOT::Math::RotationY yRot;
|
||||||
|
a1PVec *= 1.0/a1Vec.P();
|
||||||
|
zRot.SetAngle(a1PVec.Phi());
|
||||||
|
yRot.SetAngle(a1PVec.Theta());
|
||||||
|
double relAngle = std::acos(params.cosRelativeAngle);
|
||||||
|
ROOT::Math::XYZVector a1PVecAligned = yRot.Inverse() * (zRot.Inverse() * a1PVec);
|
||||||
|
a2PVec.SetXYZ(
|
||||||
|
std::sin(relAngle),
|
||||||
|
0.0,
|
||||||
|
std::cos(relAngle)
|
||||||
|
);
|
||||||
|
a2PVec = zRot * (yRot * a2PVec);
|
||||||
|
params.decay2Theta = a2PVec.Theta();
|
||||||
|
params.decay2Phi = a2PVec.Phi();
|
||||||
|
}
|
||||||
|
|
||||||
|
void CoupledThreeStepSystem::RunSystem()
|
||||||
|
{
|
||||||
|
CoupledThreeStepParameters params;
|
||||||
|
do
|
||||||
|
{
|
||||||
|
params = SampleParameters();
|
||||||
|
}
|
||||||
|
while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx)
|
||||||
|
&& m_step2.CheckDecayThreshold(params.residEx, params.decay1Ex)
|
||||||
|
&& m_step3.CheckDecayThreshold(params.decay1Ex, params.decay2Ex)));
|
||||||
|
|
||||||
|
m_step1.SetReactionDepth(params.rxnDepth);
|
||||||
|
m_step1.SetBeamKE(params.beamEnergy);
|
||||||
|
m_step1.SetPolarRxnAngle(params.rxnTheta);
|
||||||
|
m_step1.SetAzimRxnAngle(params.rxnPhi);
|
||||||
|
m_step1.SetExcitation(params.residEx);
|
||||||
|
|
||||||
|
m_step2.SetReactionDepth(params.rxnDepth);
|
||||||
|
m_step2.SetPolarRxnAngle(params.decay1Theta);
|
||||||
|
m_step2.SetAzimRxnAngle(params.decay1Phi);
|
||||||
|
m_step2.SetExcitation(params.decay1Ex);
|
||||||
|
|
||||||
|
m_step3.SetReactionDepth(params.rxnDepth);
|
||||||
|
|
||||||
|
m_step3.SetExcitation(params.decay2Ex);
|
||||||
|
|
||||||
|
m_step1.Calculate();
|
||||||
|
|
||||||
|
if(params.cosdecay1Theta == -10)
|
||||||
|
{
|
||||||
|
m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
|
||||||
|
m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
|
||||||
|
m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
|
||||||
|
m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
m_step2.Calculate();
|
||||||
|
|
||||||
|
m_step3.SetResidualEnergyLoss(true);
|
||||||
|
|
||||||
|
SampleCoupling(params); //must be called here as frame needs to be set
|
||||||
|
|
||||||
|
m_step3.SetPolarRxnAngle(params.decay2Theta);
|
||||||
|
m_step3.SetAzimRxnAngle(params.decay2Phi);
|
||||||
|
m_step3.Calculate();
|
||||||
|
}
|
||||||
|
}
|
46
src/Mask/CoupledThreeStepSystem.h
Normal file
46
src/Mask/CoupledThreeStepSystem.h
Normal file
|
@ -0,0 +1,46 @@
|
||||||
|
#ifndef COUPLEDTHREESTEPSYSTEM_H
|
||||||
|
#define COUPLEDTHREESTEPSYSTEM_H
|
||||||
|
|
||||||
|
#include "ReactionSystem.h"
|
||||||
|
#include "AngularDistribution.h"
|
||||||
|
|
||||||
|
namespace Mask {
|
||||||
|
|
||||||
|
struct CoupledThreeStepParameters
|
||||||
|
{
|
||||||
|
double beamEnergy = 0.;
|
||||||
|
double rxnTheta = 0.;
|
||||||
|
double rxnPhi = 0.;
|
||||||
|
double cosdecay1Theta = 0.;
|
||||||
|
double cosRelativeAngle = 0.;
|
||||||
|
double decay1Theta = 0.;
|
||||||
|
double decay1Phi = 0.;
|
||||||
|
double decay2Theta = 0.;
|
||||||
|
double decay2Phi = 0.;
|
||||||
|
double residEx = 0.;
|
||||||
|
double decay1Ex = 0.;
|
||||||
|
double decay2Ex = 0.;
|
||||||
|
double rxnDepth = 0.;
|
||||||
|
};
|
||||||
|
|
||||||
|
class CoupledThreeStepSystem : public ReactionSystem
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
CoupledThreeStepSystem(const std::vector<StepParameters>& params);
|
||||||
|
~CoupledThreeStepSystem();
|
||||||
|
|
||||||
|
virtual void SetLayeredTarget(const LayeredTarget& target) override;
|
||||||
|
virtual void RunSystem() override;
|
||||||
|
|
||||||
|
protected:
|
||||||
|
void Init(const std::vector<StepParameters>& params);
|
||||||
|
void SetSystemEquation() override;
|
||||||
|
CoupledThreeStepParameters SampleParameters();
|
||||||
|
void SampleCoupling(CoupledThreeStepParameters& params);
|
||||||
|
|
||||||
|
Reaction m_step1, m_step2, m_step3;
|
||||||
|
};
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
|
@ -31,9 +31,9 @@ namespace Mask {
|
||||||
double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth);
|
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);
|
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 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; }
|
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; }
|
const std::string& GetName() { return m_name; }
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
|
|
@ -1,185 +1,80 @@
|
||||||
#include "MaskApp.h"
|
#include "MaskApp.h"
|
||||||
|
#include "ConfigSerializer.h"
|
||||||
#include <fstream>
|
#include <fstream>
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
|
|
||||||
#include "TFile.h"
|
#include "TFile.h"
|
||||||
#include "TTree.h"
|
#include "TTree.h"
|
||||||
|
|
||||||
|
#include "yaml-cpp/yaml.h"
|
||||||
|
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
|
|
||||||
MaskApp::MaskApp() :
|
MaskApp::MaskApp() :
|
||||||
m_system(nullptr), m_resources(nullptr)
|
m_resources(nullptr)
|
||||||
{
|
{
|
||||||
std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
|
std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<<std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
MaskApp::~MaskApp()
|
MaskApp::~MaskApp()
|
||||||
{
|
{
|
||||||
delete m_system;
|
|
||||||
for(std::size_t i=0; i<m_systemList.size(); i++)
|
for(std::size_t i=0; i<m_systemList.size(); i++)
|
||||||
delete m_systemList[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;
|
std::cout << "Loading configuration in " << filename << std::endl;
|
||||||
|
|
||||||
std::ifstream input(filename);
|
if (!ConfigSerializer::DeserializeConfig(filename, m_params))
|
||||||
if(!input.is_open())
|
|
||||||
{
|
{
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
|
|
||||||
std::string junk;
|
//Initialize the system
|
||||||
std::getline(input, junk);
|
//Reaction chain
|
||||||
input>>junk>>m_outputName;
|
for(uint32_t i=0; i<m_params.nThreads; i++)
|
||||||
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<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)
|
|
||||||
{
|
{
|
||||||
if(junk == "begin_chain")
|
m_systemList.push_back(CreateSystem(m_params.chainParams));
|
||||||
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())
|
if(m_systemList.back() == nullptr || !m_systemList.back()->IsValid())
|
||||||
{
|
{
|
||||||
std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
|
std::cerr<<"Failure to parse reaction system... configuration not loaded."<<std::endl;
|
||||||
return false;
|
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)
|
for(auto system : m_systemList)
|
||||||
{
|
{
|
||||||
system->SetLayeredTarget(target);
|
system->SetLayeredTarget(m_params.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_system->GetSystemEquation()<<std::endl;
|
|
||||||
std::cout<<"Number of samples: "<<m_nsamples<<std::endl;
|
|
||||||
|
|
||||||
|
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;
|
||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
||||||
//Not implemented... yet!
|
|
||||||
bool MaskApp::SaveConfig(const std::string& filename)
|
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()
|
void MaskApp::Run()
|
||||||
{
|
{
|
||||||
std::cout<<"Running simulation..."<<std::endl;
|
std::cout<<"Running simulation..."<<std::endl;
|
||||||
if(m_systemList.size() != m_nthreads)
|
if(m_systemList.size() != m_params.nThreads)
|
||||||
{
|
{
|
||||||
std::cerr << "System list not equal to number of threads" << std::endl;
|
std::cerr << "System list not equal to number of threads" << std::endl;
|
||||||
return;
|
return;
|
||||||
|
@ -205,7 +100,7 @@ namespace Mask {
|
||||||
|
|
||||||
uint64_t count = 0;
|
uint64_t count = 0;
|
||||||
double percent = 0.05;
|
double percent = 0.05;
|
||||||
uint64_t flushVal = m_nsamples*percent;
|
uint64_t flushVal = m_params.nSamples*percent;
|
||||||
uint64_t flushCount = 0;
|
uint64_t flushCount = 0;
|
||||||
while(true)
|
while(true)
|
||||||
{
|
{
|
||||||
|
@ -227,42 +122,4 @@ namespace Mask {
|
||||||
std::cout<<"---------------------------------------------"<<std::endl;
|
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,6 +14,15 @@
|
||||||
|
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
|
|
||||||
|
struct AppParameters
|
||||||
|
{
|
||||||
|
std::string outputFileName = "";
|
||||||
|
uint64_t nSamples = 0;
|
||||||
|
uint32_t nThreads = 1;
|
||||||
|
std::vector<StepParameters> chainParams;
|
||||||
|
LayeredTarget target;
|
||||||
|
};
|
||||||
|
|
||||||
class MaskApp
|
class MaskApp
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
|
@ -23,17 +32,9 @@ namespace Mask {
|
||||||
bool SaveConfig(const std::string& filename);
|
bool SaveConfig(const std::string& filename);
|
||||||
|
|
||||||
void Run();
|
void Run();
|
||||||
void RunSingleThread();
|
|
||||||
|
|
||||||
private:
|
private:
|
||||||
void RunChunk(ReactionSystem* system);
|
AppParameters m_params;
|
||||||
|
|
||||||
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<ReactionSystem*> m_systemList; //One system for each thread
|
||||||
std::vector<uint64_t> m_chunkSamples;
|
std::vector<uint64_t> m_chunkSamples;
|
||||||
|
|
|
@ -12,6 +12,7 @@
|
||||||
#include <string>
|
#include <string>
|
||||||
#include <vector>
|
#include <vector>
|
||||||
#include "Math/Vector4D.h"
|
#include "Math/Vector4D.h"
|
||||||
|
#include "Math/Point3D.h"
|
||||||
#include "MassLookup.h"
|
#include "MassLookup.h"
|
||||||
|
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
|
@ -48,6 +49,8 @@ namespace Mask {
|
||||||
double detectedKE = 0.0;
|
double detectedKE = 0.0;
|
||||||
double detectedTheta = 0.0;
|
double detectedTheta = 0.0;
|
||||||
double detectedPhi = 0.0;
|
double detectedPhi = 0.0;
|
||||||
|
|
||||||
|
ROOT::Math::XYZPoint detectedPos = ROOT::Math::XYZPoint(0., 0., 0.);
|
||||||
};
|
};
|
||||||
|
|
||||||
Nucleus CreateNucleus(uint32_t z, uint32_t a);
|
Nucleus CreateNucleus(uint32_t z, uint32_t a);
|
||||||
|
|
|
@ -229,7 +229,6 @@ namespace Mask {
|
||||||
m_residual->vec4 = m_target->vec4 - m_ejectile->vec4;
|
m_residual->vec4 = m_target->vec4 - m_ejectile->vec4;
|
||||||
|
|
||||||
//energy loss for the *light* break up nucleus
|
//energy loss for the *light* break up nucleus
|
||||||
double keorig = m_ejectile->GetKE();
|
|
||||||
double ejectKE = m_ejectile->GetKE() -
|
double ejectKE = m_ejectile->GetKE() -
|
||||||
m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, m_ejectile->GetKE(), m_rxnLayer, m_ejectile->vec4.Theta(), m_rxnDepth);
|
m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, m_ejectile->GetKE(), m_rxnLayer, m_ejectile->vec4.Theta(), m_rxnDepth);
|
||||||
if(ejectKE > 0.0)
|
if(ejectKE > 0.0)
|
||||||
|
@ -257,6 +256,27 @@ namespace Mask {
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
bool Reaction::CheckReactionThreshold(double beamEnergy, double residEx)
|
||||||
|
{
|
||||||
|
double Q = m_target->groundStateMass + m_projectile->groundStateMass -
|
||||||
|
(m_ejectile->groundStateMass + m_residual->groundStateMass + residEx);
|
||||||
|
double Ethresh = -Q*(m_ejectile->groundStateMass+m_residual->groundStateMass) /
|
||||||
|
(m_ejectile->groundStateMass + m_residual->groundStateMass - m_projectile->groundStateMass);
|
||||||
|
if (beamEnergy < Ethresh)
|
||||||
|
return false;
|
||||||
|
else
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
|
bool Reaction::CheckDecayThreshold(double targetEx, double residEx)
|
||||||
|
{
|
||||||
|
double Q = m_target->groundStateMass + targetEx -
|
||||||
|
(m_ejectile->groundStateMass + m_residual->groundStateMass + residEx);
|
||||||
|
if ( Q < 0.0)
|
||||||
|
return false;
|
||||||
|
else
|
||||||
|
return true;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -42,6 +42,9 @@ namespace Mask {
|
||||||
void SetRxnLayer(std::size_t layer) { m_rxnLayer = layer; };
|
void SetRxnLayer(std::size_t layer) { m_rxnLayer = layer; };
|
||||||
void SetResidualEnergyLoss(bool isEloss) { m_isResidEloss = isEloss; };
|
void SetResidualEnergyLoss(bool isEloss) { m_isResidEloss = isEloss; };
|
||||||
|
|
||||||
|
bool CheckReactionThreshold(double beamEnergy, double residualEx);
|
||||||
|
bool CheckDecayThreshold(double targetEx, double residualEx);
|
||||||
|
|
||||||
bool IsDecay() const { return m_isDecay; };
|
bool IsDecay() const { return m_isDecay; };
|
||||||
|
|
||||||
std::size_t GetRxnLayer() const { return m_rxnLayer; };
|
std::size_t GetRxnLayer() const { return m_rxnLayer; };
|
||||||
|
|
|
@ -37,6 +37,17 @@ namespace Mask {
|
||||||
return RxnType::None;
|
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)
|
static RxnThetaType StringToRxnThetaType(const std::string& type)
|
||||||
{
|
{
|
||||||
if(type == "CenterOfMass")
|
if(type == "CenterOfMass")
|
||||||
|
@ -52,6 +63,17 @@ namespace Mask {
|
||||||
return RxnThetaType::None;
|
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
|
#endif
|
|
@ -31,11 +31,11 @@ namespace Mask {
|
||||||
double GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle);
|
double GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle);
|
||||||
double GetEnergyLossFractionalDepth(int zp, int ap, double startEnergy, double angle, double percent_depth);
|
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);
|
double GetReverseEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double angle, double percent_depth);
|
||||||
inline const double& GetThickness() { return m_thickness; }
|
const double& GetThickness() const { return m_thickness; }
|
||||||
inline int GetNumberOfElements() { return m_Z.size(); }
|
int GetNumberOfElements() const { return m_Z.size(); }
|
||||||
inline int GetElementZ(int index) { return m_Z[index]; }
|
int GetElementZ(int index) const { return m_Z[index]; }
|
||||||
inline int GetElementA(int index) { return m_A[index]; }
|
int GetElementA(int index) const { return m_A[index]; }
|
||||||
inline int GetElementStoich(int index) { return m_stoich[index]; }
|
int GetElementStoich(int index) const { return m_stoich[index]; }
|
||||||
|
|
||||||
private:
|
private:
|
||||||
void Init(const std::vector<int>& z, const std::vector<int>& a, const std::vector<int>& stoich);
|
void Init(const std::vector<int>& z, const std::vector<int>& a, const std::vector<int>& stoich);
|
||||||
|
|
|
@ -2,6 +2,11 @@
|
||||||
#include "RandomGenerator.h"
|
#include "RandomGenerator.h"
|
||||||
#include "KinematicsExceptions.h"
|
#include "KinematicsExceptions.h"
|
||||||
|
|
||||||
|
#include "Math/Boost.h"
|
||||||
|
#include "Math/Vector3D.h"
|
||||||
|
#include "Math/RotationZ.h"
|
||||||
|
#include "Math/RotationY.h"
|
||||||
|
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
|
|
||||||
ThreeStepSystem::ThreeStepSystem(const std::vector<StepParameters>& params) :
|
ThreeStepSystem::ThreeStepSystem(const std::vector<StepParameters>& params) :
|
||||||
|
@ -56,8 +61,8 @@ namespace Mask {
|
||||||
m_nuclei[3] = CreateNucleus(zr, ar); //residual
|
m_nuclei[3] = CreateNucleus(zr, ar); //residual
|
||||||
m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
|
m_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
|
||||||
m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
|
m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
|
||||||
m_nuclei[5] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
|
m_nuclei[6] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
|
||||||
m_nuclei[5] = CreateNucleus(zb4, ab4); //breakup4
|
m_nuclei[7] = CreateNucleus(zb4, ab4); //breakup4
|
||||||
|
|
||||||
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
|
m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
|
||||||
m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
|
m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
|
||||||
|
@ -115,43 +120,57 @@ namespace Mask {
|
||||||
m_sysEquation = stream.str();
|
m_sysEquation = stream.str();
|
||||||
}
|
}
|
||||||
|
|
||||||
|
ThreeStepParameters ThreeStepSystem::SampleParameters()
|
||||||
|
{
|
||||||
|
ThreeStepParameters params;
|
||||||
|
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
|
||||||
|
params.beamEnergy = m_beamDistributions[0](gen);
|
||||||
|
params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
|
||||||
|
params.rxnPhi = m_phiRanges[0](gen);
|
||||||
|
params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
|
||||||
|
params.decay1Theta = std::acos(params.cosdecay1Theta);
|
||||||
|
params.decay1Phi = m_phiRanges[1](gen);
|
||||||
|
params.cosdecay2Theta = m_decayAngularDistributions[1].GetRandomCosTheta();
|
||||||
|
params.decay2Theta = std::acos(params.cosdecay2Theta);
|
||||||
|
params.decay2Phi = m_phiRanges[1](gen);
|
||||||
|
params.residEx = m_exDistributions[0](gen);
|
||||||
|
params.decay1Ex = m_exDistributions[1](gen);
|
||||||
|
params.decay2Ex = m_exDistributions[2](gen);
|
||||||
|
params.rxnDepth = (m_rxnDepthDist(gen));
|
||||||
|
return params;
|
||||||
|
}
|
||||||
|
|
||||||
void ThreeStepSystem::RunSystem()
|
void ThreeStepSystem::RunSystem()
|
||||||
{
|
{
|
||||||
//Sample parameters
|
ThreeStepParameters params;
|
||||||
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
|
do
|
||||||
double bke = m_beamDistributions[0](gen);
|
{
|
||||||
double rxnTheta = std::acos((m_thetaRanges[0])(gen));
|
params = SampleParameters();
|
||||||
double rxnPhi = m_phiRanges[0](gen);
|
}
|
||||||
double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
|
while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx)
|
||||||
double decay1Theta = std::acos(decay1costheta);
|
&& m_step2.CheckDecayThreshold(params.residEx, params.decay1Ex)
|
||||||
double decay1Phi = m_phiRanges[1](gen);
|
&& m_step3.CheckDecayThreshold(params.decay1Ex, params.decay2Ex)));
|
||||||
double decay2costheta = m_decayAngularDistributions[1].GetRandomCosTheta();
|
|
||||||
double decay2Theta = std::acos(decay2costheta);
|
|
||||||
double decay2Phi = m_phiRanges[2](gen);
|
|
||||||
double residEx = m_exDistributions[0](gen);
|
|
||||||
double decay1Ex = m_exDistributions[1](gen);
|
|
||||||
double decay2Ex = m_exDistributions[2](gen);
|
|
||||||
double rxnDepth = (m_rxnDepthDist(gen));
|
|
||||||
|
|
||||||
m_step1.SetReactionDepth(rxnDepth);
|
m_step1.SetReactionDepth(params.rxnDepth);
|
||||||
m_step1.SetBeamKE(bke);
|
m_step1.SetBeamKE(params.beamEnergy);
|
||||||
m_step1.SetPolarRxnAngle(rxnTheta);
|
m_step1.SetPolarRxnAngle(params.rxnTheta);
|
||||||
m_step1.SetAzimRxnAngle(rxnPhi);
|
m_step1.SetAzimRxnAngle(params.rxnPhi);
|
||||||
m_step1.SetExcitation(residEx);
|
m_step1.SetExcitation(params.residEx);
|
||||||
|
|
||||||
m_step2.SetReactionDepth(rxnDepth);
|
m_step2.SetReactionDepth(params.rxnDepth);
|
||||||
m_step2.SetPolarRxnAngle(decay1Theta);
|
m_step2.SetPolarRxnAngle(params.decay1Theta);
|
||||||
m_step2.SetAzimRxnAngle(decay1Phi);
|
m_step2.SetAzimRxnAngle(params.decay1Phi);
|
||||||
m_step2.SetExcitation(decay1Ex);
|
m_step2.SetExcitation(params.decay1Ex);
|
||||||
|
|
||||||
m_step3.SetReactionDepth(rxnDepth);
|
m_step3.SetReactionDepth(params.rxnDepth);
|
||||||
m_step3.SetPolarRxnAngle(decay2Theta);
|
m_step3.SetPolarRxnAngle(params.decay2Theta);
|
||||||
m_step3.SetAzimRxnAngle(decay2Phi);
|
m_step3.SetAzimRxnAngle(params.decay2Phi);
|
||||||
m_step3.SetExcitation(decay2Ex);
|
m_step3.SetResidualEnergyLoss(true);
|
||||||
|
m_step3.SetExcitation(params.decay2Ex);
|
||||||
|
|
||||||
m_step1.Calculate();
|
m_step1.Calculate();
|
||||||
|
|
||||||
if(decay1costheta == -10)
|
if(params.cosdecay1Theta == -10)
|
||||||
{
|
{
|
||||||
m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
|
m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
|
||||||
m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
|
m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
|
||||||
|
@ -161,15 +180,13 @@ namespace Mask {
|
||||||
}
|
}
|
||||||
m_step2.Calculate();
|
m_step2.Calculate();
|
||||||
|
|
||||||
if(decay2costheta == -10)
|
if(params.cosdecay2Theta == -10)
|
||||||
{
|
{
|
||||||
m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
|
m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
|
||||||
m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
|
m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
m_step3.SetResidualEnergyLoss(true);
|
|
||||||
m_step3.Calculate();
|
m_step3.Calculate();
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
|
@ -6,6 +6,23 @@
|
||||||
|
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
|
|
||||||
|
struct ThreeStepParameters
|
||||||
|
{
|
||||||
|
double beamEnergy = 0.;
|
||||||
|
double rxnTheta = 0.;
|
||||||
|
double rxnPhi = 0.;
|
||||||
|
double cosdecay1Theta = 0.;
|
||||||
|
double cosdecay2Theta = 0.;
|
||||||
|
double decay1Theta = 0.;
|
||||||
|
double decay1Phi = 0.;
|
||||||
|
double decay2Theta = 0.;
|
||||||
|
double decay2Phi = 0.;
|
||||||
|
double residEx = 0.;
|
||||||
|
double decay1Ex = 0.;
|
||||||
|
double decay2Ex = 0.;
|
||||||
|
double rxnDepth = 0.;
|
||||||
|
};
|
||||||
|
|
||||||
class ThreeStepSystem : public ReactionSystem
|
class ThreeStepSystem : public ReactionSystem
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
|
@ -18,6 +35,7 @@ namespace Mask {
|
||||||
protected:
|
protected:
|
||||||
void Init(const std::vector<StepParameters>& params);
|
void Init(const std::vector<StepParameters>& params);
|
||||||
void SetSystemEquation() override;
|
void SetSystemEquation() override;
|
||||||
|
ThreeStepParameters SampleParameters();
|
||||||
|
|
||||||
Reaction m_step1, m_step2, m_step3;
|
Reaction m_step1, m_step2, m_step3;
|
||||||
|
|
||||||
|
|
|
@ -93,34 +93,58 @@ namespace Mask {
|
||||||
m_sysEquation = stream.str();
|
m_sysEquation = stream.str();
|
||||||
}
|
}
|
||||||
|
|
||||||
|
TwoStepParameters TwoStepSystem::SampleParameters()
|
||||||
|
{
|
||||||
|
TwoStepParameters params;
|
||||||
|
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
|
||||||
|
params.beamEnergy = (m_beamDistributions[0])(gen);
|
||||||
|
params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
|
||||||
|
params.rxnPhi = (m_phiRanges[0])(gen);
|
||||||
|
params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
|
||||||
|
params.decay1Theta = std::acos(params.cosdecay1Theta);
|
||||||
|
params.decay1Phi = m_phiRanges[1](gen);
|
||||||
|
params.residEx = (m_exDistributions[0])(gen);
|
||||||
|
params.decay2Ex = m_exDistributions[1](gen);
|
||||||
|
params.rxnDepth = (m_rxnDepthDist(gen));
|
||||||
|
return params;
|
||||||
|
}
|
||||||
|
|
||||||
void TwoStepSystem::RunSystem()
|
void TwoStepSystem::RunSystem()
|
||||||
{
|
{
|
||||||
//Sample parameters
|
//Sample parameters
|
||||||
std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
|
// std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
|
||||||
double bke = (m_beamDistributions[0])(gen);
|
// double bke = (m_beamDistributions[0])(gen);
|
||||||
double rxnTheta = std::acos((m_thetaRanges[0])(gen));
|
// double rxnTheta = std::acos((m_thetaRanges[0])(gen));
|
||||||
double rxnPhi = (m_phiRanges[0])(gen);
|
// double rxnPhi = (m_phiRanges[0])(gen);
|
||||||
double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
|
// double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
|
||||||
double decay1Theta = std::acos(decay1costheta);
|
// double decay1Theta = std::acos(decay1costheta);
|
||||||
double decay1Phi = m_phiRanges[1](gen);
|
// double decay1Phi = m_phiRanges[1](gen);
|
||||||
double residEx = (m_exDistributions[0])(gen);
|
// double residEx = (m_exDistributions[0])(gen);
|
||||||
double decay2Ex = m_exDistributions[1](gen);
|
// double decay2Ex = m_exDistributions[1](gen);
|
||||||
double rxnDepth = (m_rxnDepthDist(gen));
|
// double rxnDepth = (m_rxnDepthDist(gen));
|
||||||
|
|
||||||
m_step1.SetReactionDepth(rxnDepth);
|
TwoStepParameters params;
|
||||||
m_step1.SetBeamKE(bke);
|
do
|
||||||
m_step1.SetPolarRxnAngle(rxnTheta);
|
{
|
||||||
m_step1.SetAzimRxnAngle(rxnPhi);
|
params = SampleParameters();
|
||||||
m_step1.SetExcitation(residEx);
|
}
|
||||||
|
while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx)
|
||||||
|
&& m_step2.CheckDecayThreshold(params.residEx, params.decay2Ex)));
|
||||||
|
|
||||||
m_step2.SetReactionDepth(rxnDepth);
|
m_step1.SetReactionDepth(params.rxnDepth);
|
||||||
m_step2.SetPolarRxnAngle(decay1Theta);
|
m_step1.SetBeamKE(params.beamEnergy);
|
||||||
m_step2.SetAzimRxnAngle(decay1Phi);
|
m_step1.SetPolarRxnAngle(params.rxnTheta);
|
||||||
m_step2.SetExcitation(decay2Ex);
|
m_step1.SetAzimRxnAngle(params.rxnPhi);
|
||||||
|
m_step1.SetExcitation(params.residEx);
|
||||||
|
|
||||||
|
m_step2.SetReactionDepth(params.rxnDepth);
|
||||||
|
m_step2.SetPolarRxnAngle(params.decay1Theta);
|
||||||
|
m_step2.SetAzimRxnAngle(params.decay1Phi);
|
||||||
|
m_step2.SetExcitation(params.decay2Ex);
|
||||||
|
|
||||||
m_step1.Calculate();
|
m_step1.Calculate();
|
||||||
|
|
||||||
if(decay1costheta == -10)
|
if(params.cosdecay1Theta == -10)
|
||||||
{
|
{
|
||||||
m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
|
m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
|
||||||
m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
|
m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
|
||||||
|
|
|
@ -6,6 +6,19 @@
|
||||||
|
|
||||||
namespace Mask {
|
namespace Mask {
|
||||||
|
|
||||||
|
struct TwoStepParameters
|
||||||
|
{
|
||||||
|
double beamEnergy = 0.;
|
||||||
|
double rxnTheta = 0.;
|
||||||
|
double rxnPhi = 0.;
|
||||||
|
double cosdecay1Theta = 0.;
|
||||||
|
double decay1Theta = 0.;
|
||||||
|
double decay1Phi = 0.;
|
||||||
|
double residEx = 0.;
|
||||||
|
double decay2Ex = 0.;
|
||||||
|
double rxnDepth = 0.;
|
||||||
|
};
|
||||||
|
|
||||||
class TwoStepSystem : public ReactionSystem
|
class TwoStepSystem : public ReactionSystem
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
|
@ -18,6 +31,7 @@ namespace Mask {
|
||||||
private:
|
private:
|
||||||
void Init(const std::vector<StepParameters>& params);
|
void Init(const std::vector<StepParameters>& params);
|
||||||
void SetSystemEquation() override;
|
void SetSystemEquation() override;
|
||||||
|
TwoStepParameters SampleParameters();
|
||||||
|
|
||||||
Reaction m_step1, m_step2;
|
Reaction m_step1, m_step2;
|
||||||
};
|
};
|
||||||
|
|
|
@ -1,7 +1,8 @@
|
||||||
#include "RootPlotter.h"
|
#include "RootPlotter.h"
|
||||||
#include <TFile.h>
|
#include <TFile.h>
|
||||||
#include <TTree.h>
|
#include <TTree.h>
|
||||||
#include <TVector3.h>
|
#include <Math/Vector3D.h>
|
||||||
|
#include "Math/Boost.h"
|
||||||
|
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
|
|
||||||
|
@ -10,7 +11,8 @@ static double FullPhi(double phi)
|
||||||
return phi < 0.0 ? 2.0*M_PI + phi : phi;
|
return phi < 0.0 ? 2.0*M_PI + phi : phi;
|
||||||
}
|
}
|
||||||
|
|
||||||
RootPlotter::RootPlotter()
|
RootPlotter::RootPlotter() :
|
||||||
|
m_target({5},{9},{1}, 74.0)
|
||||||
{
|
{
|
||||||
TH1::AddDirectory(kFALSE);
|
TH1::AddDirectory(kFALSE);
|
||||||
}
|
}
|
||||||
|
@ -42,10 +44,16 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam
|
||||||
std::cout<<"\rPercent of data processed: "<<flushCount*flushFrac*100<<"%"<<std::flush;
|
std::cout<<"\rPercent of data processed: "<<flushCount*flushFrac*100<<"%"<<std::flush;
|
||||||
}
|
}
|
||||||
tree->GetEntry(i);
|
tree->GetEntry(i);
|
||||||
for(Mask::Nucleus& nuc : *(dataHandle))
|
// for(Mask::Nucleus& nuc : *(dataHandle))
|
||||||
|
// {
|
||||||
|
// FillData(nuc);
|
||||||
|
// }
|
||||||
|
for(int i=0; i<dataHandle->size(); i++)
|
||||||
{
|
{
|
||||||
FillData(nuc);
|
FillData(dataHandle->at(i), i);
|
||||||
}
|
}
|
||||||
|
//Don't leave this in!
|
||||||
|
//Correlations(*dataHandle);
|
||||||
}
|
}
|
||||||
std::cout<<std::endl;
|
std::cout<<std::endl;
|
||||||
input->Close();
|
input->Close();
|
||||||
|
@ -57,40 +65,142 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam
|
||||||
output->Close();
|
output->Close();
|
||||||
}
|
}
|
||||||
|
|
||||||
void RootPlotter::FillData(const Mask::Nucleus& nuc)
|
void RootPlotter::FillData(const Mask::Nucleus& nuc, int i)
|
||||||
{
|
{
|
||||||
std::string modifier = "";
|
std::string mod = "_detected";
|
||||||
if(nuc.isDetected)
|
std::string num = "_" + std::to_string(i);
|
||||||
modifier = "detected";
|
|
||||||
|
|
||||||
std::string sym = nuc.isotopicSymbol;
|
std::string sym = nuc.isotopicSymbol;
|
||||||
std::string ke_vs_th_name = sym + modifier + "_ke_vs_theta";
|
std::string ke_vs_th_name = sym + num + "_ke_vs_theta";
|
||||||
std::string ke_vs_th_title = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
|
std::string ke_vs_th_title = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
|
||||||
std::string ke_vs_ph_name = sym + modifier + "_ke_vs_phi";
|
std::string ke_vs_th_name_det = sym + num + mod + "_ke_vs_theta";
|
||||||
|
std::string ke_vs_th_title_det = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
|
||||||
|
std::string ke_vs_ph_name = sym + num + "_ke_vs_phi";
|
||||||
std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
|
std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
|
||||||
std::string th_vs_ph_name = sym + modifier + "_theta_vs_phi";
|
std::string ke_vs_ph_name_det = sym + num + mod + "_ke_vs_phi";
|
||||||
|
std::string ke_vs_ph_title_det = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
|
||||||
|
std::string th_vs_ph_name = sym + num + "_theta_vs_phi";
|
||||||
std::string th_vs_ph_title = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
|
std::string th_vs_ph_title = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
|
||||||
std::string ex_name = sym + modifier + "_ex";
|
std::string th_vs_ph_name_det = sym + num + mod + "_theta_vs_phi";
|
||||||
|
std::string th_vs_ph_title_det = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
|
||||||
|
std::string ex_name = sym + num + "_ex";
|
||||||
std::string ex_title = ex_name + ";E_{ex} (MeV);counts";
|
std::string ex_title = ex_name + ";E_{ex} (MeV);counts";
|
||||||
std::string angdist_name = sym + modifier +"_angDist";
|
std::string ex_name_det = sym + num + mod + "_ex";
|
||||||
|
std::string ex_title_det = ex_name + ";E_{ex} (MeV);counts";
|
||||||
|
std::string angdist_name = sym + num + "_angDist";
|
||||||
std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts";
|
std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts";
|
||||||
|
std::string angdist_name_det = sym + num + mod +"_angDist";
|
||||||
|
std::string angdist_title_det = angdist_name+";cos#right(#theta_{CM}#left);counts";
|
||||||
|
std::string hist_ke_th_name = sym + num + "_hist_ke_vs_theta";
|
||||||
|
std::string hist_ke_th_title = hist_ke_th_name + ";#theta_{lab};Kinetic Energy (MeV)";
|
||||||
|
std::string hist_ke_th_name_det = sym + num + mod + "_hist_ke_vs_theta";
|
||||||
|
std::string hist_ke_th_title_det = hist_ke_th_name + ";#theta_{lab};Kinetic Energy (MeV)";
|
||||||
|
|
||||||
MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
|
MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
|
||||||
MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.GetKE(), 4);
|
MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.GetKE(), 4);
|
||||||
MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
|
MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
|
||||||
MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
|
MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
|
||||||
MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
|
MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
|
||||||
if(nuc.isDetected)
|
MyFill(hist_ke_th_name, hist_ke_th_title, 180, 0.0, 180.0, 400, 0.0, 20.0, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE());
|
||||||
|
if(nuc.isDetected && nuc.detectedKE > 0.25)
|
||||||
{
|
{
|
||||||
MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
|
MyFill(ke_vs_th_name_det, ke_vs_th_title_det, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2);
|
||||||
MyFill(ke_vs_ph_name, ke_vs_ph_title, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4);
|
MyFill(ke_vs_ph_name_det, ke_vs_ph_title_det, FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4);
|
||||||
MyFill(th_vs_ph_name, th_vs_ph_title, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
|
MyFill(th_vs_ph_name_det, th_vs_ph_title_det, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2);
|
||||||
MyFill(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy());
|
MyFill(ex_name_det, ex_title_det, 260, -1.0, 25, nuc.GetExcitationEnergy());
|
||||||
MyFill(angdist_name, angdist_title, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
|
MyFill(angdist_name_det, angdist_title_det, 20, -1.0, 1.0, std::cos(nuc.thetaCM));
|
||||||
|
MyFill(hist_ke_th_name_det, hist_ke_th_title_det, 180, 0.0, 180.0, 400, 0.0, 20.0, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE);
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
//Correlation analysis for coupled step reactions, will need to be customized
|
||||||
|
//for a given reaction. By default it is off.
|
||||||
|
void RootPlotter::Correlations(const std::vector<Mask::Nucleus>& event)
|
||||||
|
{
|
||||||
|
static constexpr double lowerBoundEx = -2.7;
|
||||||
|
static constexpr double upperBoundEx = 2.7;
|
||||||
|
static constexpr double targetDepth = 0.5;
|
||||||
|
const Mask::Nucleus& secondary = event[6];
|
||||||
|
const Mask::Nucleus& primary = event[4];
|
||||||
|
const Mask::Nucleus& parent = event[3];
|
||||||
|
const Mask::Nucleus& li = event[5];
|
||||||
|
if(secondary.vec4.P() == 0.0 || primary.vec4.P() == 0.0)
|
||||||
|
{
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
ROOT::Math::Boost boostParent(parent.vec4.BoostToCM());
|
||||||
|
ROOT::Math::Boost boostIntermediate(li.vec4.BoostToCM());
|
||||||
|
ROOT::Math::PxPyPzEVector a2Vec = (secondary.vec4);
|
||||||
|
ROOT::Math::PxPyPzEVector a1Vec = (primary.vec4);
|
||||||
|
|
||||||
|
double a2KE = secondary.GetKE() + m_target.GetReverseEnergyLossFractionalDepth(secondary.Z, secondary.A, secondary.GetKE(), a2Vec.Theta(), targetDepth);
|
||||||
|
double a2P = std::sqrt(a2KE * (a2KE + 2.0 * secondary.groundStateMass));
|
||||||
|
a2Vec.SetPxPyPzE(
|
||||||
|
a2P * std::sin(a2Vec.Theta()) * std::cos(a2Vec.Phi()),
|
||||||
|
a2P * std::sin(a2Vec.Theta()) * std::sin(a2Vec.Phi()),
|
||||||
|
a2P * std::cos(a2Vec.Theta()),
|
||||||
|
a2KE + secondary.groundStateMass
|
||||||
|
);
|
||||||
|
|
||||||
|
double a1KE = primary.GetKE() + m_target.GetReverseEnergyLossFractionalDepth(primary.Z, primary.A, primary.GetKE(), a1Vec.Theta(), targetDepth);
|
||||||
|
double a1P = std::sqrt(a1KE * (a1KE + 2.0 * primary.groundStateMass));
|
||||||
|
a1Vec.SetPxPyPzE(
|
||||||
|
a1P * std::sin(a1Vec.Theta()) * std::cos(a1Vec.Phi()),
|
||||||
|
a1P * std::sin(a1Vec.Theta()) * std::sin(a1Vec.Phi()),
|
||||||
|
a1P * std::cos(a1Vec.Theta()),
|
||||||
|
a1KE + primary.groundStateMass
|
||||||
|
);
|
||||||
|
|
||||||
|
ROOT::Math::XYZVector a1PVec;
|
||||||
|
ROOT::Math::XYZVector a2PVec;
|
||||||
|
|
||||||
|
ROOT::Math::PxPyPzEVector boostedIntA1Vec = boostIntermediate * a1Vec;
|
||||||
|
ROOT::Math::PxPyPzEVector boostedIntA2Vec = boostIntermediate * a2Vec;
|
||||||
|
ROOT::Math::PxPyPzEVector boostedParA2Vec = boostParent * a2Vec;
|
||||||
|
a1PVec.SetXYZ(boostedIntA1Vec.Px(), boostedIntA1Vec.Py(), boostedIntA1Vec.Pz());
|
||||||
|
a2PVec.SetXYZ(boostedIntA2Vec.Px(), boostedIntA2Vec.Py(), boostedIntA2Vec.Pz());
|
||||||
|
double relCosThetaCM = a1PVec.Dot(a2PVec) / (a1PVec.R() * a2PVec.R());
|
||||||
|
|
||||||
|
ROOT::Math::PxPyPzEVector badLi = parent.vec4 - a2Vec;
|
||||||
|
double secondaryExcitation = badLi.M() - li.groundStateMass;
|
||||||
|
|
||||||
|
MyFill("seoncdaryEx","secondaryEx",1000,-5.0,5.0,secondaryExcitation);
|
||||||
|
MyFill("primaryEx","primaryEx",1000,-5.0,5.0,li.GetExcitationEnergy());
|
||||||
|
MyFill("relativeDist_full","relativeDist_full",20,-1.0,1.0, relCosThetaCM);
|
||||||
|
MyFill("primaryDist_full","primaryDist_full",20,-1.0,1.0,std::cos(primary.thetaCM));
|
||||||
|
MyFill("seondaryDist_full","seconaryDist_full",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
|
||||||
|
MyFill("secondaryDist_fullOwn","seondaryDist_fullOwn",20,-1.0,1.0,std::cos(secondary.thetaCM));
|
||||||
|
|
||||||
|
bool isPrimary = false;
|
||||||
|
bool isSecondary = false;
|
||||||
|
if(primary.isDetected)
|
||||||
|
{
|
||||||
|
MyFill("primaryExMatched","primaryExMatched",1000,-5.0,5.0,li.GetExcitationEnergy());
|
||||||
|
MyFill("primaryDist","primaryDist",20,-1.0,1.0,std::cos(primary.thetaCM));
|
||||||
|
isPrimary = true;
|
||||||
|
}
|
||||||
|
if( secondary.isDetected && secondaryExcitation > lowerBoundEx && secondaryExcitation < upperBoundEx )
|
||||||
|
{
|
||||||
|
MyFill("secondaryDist","secondaryDist",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
|
||||||
|
MyFill("secondaryExMatched","secondaryExMatched",1000,-5.0,5.0,secondaryExcitation);
|
||||||
|
isSecondary = true;
|
||||||
|
}
|
||||||
|
|
||||||
|
if(isSecondary && isPrimary)
|
||||||
|
{
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
else if (isSecondary)
|
||||||
|
{
|
||||||
|
MyFill("completeDist","completeDist",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
|
||||||
|
}
|
||||||
|
else if (isPrimary)
|
||||||
|
{
|
||||||
|
MyFill("completeDist","completeDist",20,-1.0,1.0,std::cos(primary.thetaCM));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val)
|
void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val)
|
||||||
{
|
{
|
||||||
auto iter = m_map.find(name);
|
auto iter = m_map.find(name);
|
||||||
|
|
|
@ -6,6 +6,7 @@
|
||||||
#include <unordered_map>
|
#include <unordered_map>
|
||||||
|
|
||||||
#include "Nucleus.h"
|
#include "Nucleus.h"
|
||||||
|
#include "Target.h"
|
||||||
|
|
||||||
#include <TH1.h>
|
#include <TH1.h>
|
||||||
#include <TH2.h>
|
#include <TH2.h>
|
||||||
|
@ -20,7 +21,8 @@ public:
|
||||||
void Run(const std::string& inputname, const std::string& outputname);
|
void Run(const std::string& inputname, const std::string& outputname);
|
||||||
|
|
||||||
private:
|
private:
|
||||||
void FillData(const Mask::Nucleus& nuc);
|
void FillData(const Mask::Nucleus& nuc, int i);
|
||||||
|
void Correlations(const std::vector<Mask::Nucleus>& event);
|
||||||
void MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val);
|
void MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val);
|
||||||
void MyFill(const std::string& name, const std::string& title, int binsx, float minx, float maxx, int binsy, float miny, float maxy,
|
void MyFill(const std::string& name, const std::string& title, int binsx, float minx, float maxx, int binsy, float miny, float maxy,
|
||||||
double valx, double valy);
|
double valx, double valy);
|
||||||
|
@ -29,6 +31,7 @@ private:
|
||||||
std::unordered_map<std::string, std::shared_ptr<TObject>> m_map;
|
std::unordered_map<std::string, std::shared_ptr<TObject>> m_map;
|
||||||
|
|
||||||
static constexpr double s_rad2deg = 180.0/M_PI;
|
static constexpr double s_rad2deg = 180.0/M_PI;
|
||||||
|
Mask::Target m_target;
|
||||||
};
|
};
|
||||||
|
|
||||||
#endif
|
#endif
|
1
src/vendor/yaml-cpp
vendored
Submodule
1
src/vendor/yaml-cpp
vendored
Submodule
|
@ -0,0 +1 @@
|
||||||
|
Subproject commit 0e6e28d1a38224fc8172fae0109ea7f673c096db
|
Loading…
Reference in New Issue
Block a user