diff --git a/.gitignore b/.gitignore index f0eaf4d..4e94ac6 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,8 @@ Makefile build/ lib/ bin/ +scripts/ +*.yaml ###Keep this file### diff --git a/.gitmodules b/.gitmodules index 0f578de..cde4811 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,6 @@ [submodule "src/vendor/catima"] path = src/vendor/catima url = https://github.com/gwm17/catima.git +[submodule "src/vendor/yaml-cpp"] + path = src/vendor/yaml-cpp + url = https://github.com/jbeder/yaml-cpp.git diff --git a/CMakeLists.txt b/CMakeLists.txt index ebf7578..e02f867 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,6 +19,8 @@ set(CMAKE_CXX_STANDARD 17) find_package(ROOT REQUIRED COMPONENTS GenVector) add_subdirectory(src/vendor/catima) +set(YAML_CPP_BUILD_TOOLS Off CACHE BOOL "Enable parse tools") +add_subdirectory(src/vendor/yaml-cpp) add_subdirectory(src/Mask) add_subdirectory(src/Kinematics) add_subdirectory(src/Detectors) diff --git a/README.md b/README.md index 6e30a0a..5c68cb9 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,12 @@ # 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. 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` - `cd build` - `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. ## 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. @@ -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: -`./bin/Kinematics input.txt` +`./bin/Kinematics .yaml` -Input.txt can be replaced by any text file with the correct format. +`` 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 + 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. @@ -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 -`./bin/Detectors ` +`./bin/Detectors .yaml` -An example input file is provided with the repository. +`.yaml` is a YAML configuration file. An example, `detector.yaml` is included in 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 @@ -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. ## Requirements + ROOT version 6.22 or greater is required CMake version 3.0 or greater is required - diff --git a/detector.yaml b/detector.yaml new file mode 100644 index 0000000..48d2a78 --- /dev/null +++ b/detector.yaml @@ -0,0 +1,5 @@ +InputDataFile: /media/data/gwm17/mask_tests/temp.root +OutputDataFile: /media/data/gwm17/mask_tests/temp_det.root +DeadChannelFile: etc/sabreDeadChannels_May2022.txt +NumberOfThreads: 5 +ArrayType: Sabre diff --git a/detector_input.txt b/detector_input.txt deleted file mode 100644 index 775494d..0000000 --- a/detector_input.txt +++ /dev/null @@ -1,5 +0,0 @@ -InputDataFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root -OutputDataFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B_Sabre.root -DeadChannelFile: etc/sabreDeadChannels_May2022.txt -NumberOfThreads: 5 -ArrayType: Sabre \ No newline at end of file diff --git a/kinematics.yaml b/kinematics.yaml new file mode 100644 index 0000000..d743212 --- /dev/null +++ b/kinematics.yaml @@ -0,0 +1,49 @@ +OutputFile: /media/data/gwm17/mask_tests/temp.root +Threads: 5 +ReactionSamples: 1000000 +ReactionChain: + - Type: Reaction + Reactants: + - Z: 5 + A: 10 + - Z: 2 + A: 3 + - Z: 2 + A: 4 + BeamEnergyMean(MeV): 24.0 + BeamEnergySigma(MeV): 0.0 + ThetaType: Lab + ThetaMin(deg): 15.0 + ThetaMax(deg): 15.0 + PhiMin(deg): 0.0 + PhiMax(deg): 0.0 + ResidualExcitationMean(MeV): 16.798 + ResidualExcitationSigma(MeV): 0.033 + - Type: Decay + Reactants: + - Z: 5 + A: 9 + - Z: 2 + A: 4 + PhiMin(deg): 0.0 + PhiMax(deg): 360.0 + ResidualExcitationMean(MeV): 0.0 + ResidualExcitationSigma(MeV): 0.522 + AngularDistributionFile: ./etc/isotropic_dist.txt + - Type: Decay + Reactants: + - Z: 3 + A: 5 + - Z: 2 + A: 4 + PhiMin(deg): 0.0 + PhiMax(deg): 360.0 + ResidualExcitationMean(MeV): 0.0 + ResidualExcitationSigma(MeV): 0.0 + AngularDistributionFile: ./etc/isotropic_dist.txt +TargetLayers: + - Thickness(ug/cm^2): 74 + Elements: + - Z: 5 + A: 10 + S: 1 \ No newline at end of file diff --git a/kinematics_input.txt b/kinematics_input.txt deleted file mode 100644 index fdbe62c..0000000 --- a/kinematics_input.txt +++ /dev/null @@ -1,45 +0,0 @@ -----------Data Information---------- -OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root -NumberOfThreads: 6 -----------Reaction Information---------- -NumberOfSamples: 1000000 -begin_chain - begin_step - Type: Reaction - begin_nuclei - 5 10 - 2 3 - 2 4 - end_nuclei - BeamEnergyMean(MeV): 24.0 - BeamEnergySigma(MeV): 0.0 - ThetaType: Lab - ThetaMin(deg): 15.0 - ThetaMax(deg): 15.0 - PhiMin(deg): 0.0 - PhMax(deg): 0.0 - ResidualExcitationMean(MeV): 16.8 - ResidualExcitationSigma(MeV): 0.023 - end_step - begin_step - Type: Decay - begin_nuclei - 5 9 - 2 4 - end_nuclei - PhiMin(deg): 0.0 - PhMax(deg): 360.0 - ResidualExcitationMean(MeV): 0.0 - ResidualExcitationSigma(MeV): 0.0 - AngularDistributionFile: ./etc/isotropic_dist.txt - end_step -end_chain -----------Target Information---------- -NumberOfLayers: 1 -begin_layer - Thickness(ug/cm^2): 74 - begin_elements (Z, A, Stoich.) - element 5 10 1 - end_elements -end_layer --------------------------------------- diff --git a/src/Detectors/CMakeLists.txt b/src/Detectors/CMakeLists.txt index b556f58..fe66023 100644 --- a/src/Detectors/CMakeLists.txt +++ b/src/Detectors/CMakeLists.txt @@ -1,5 +1,9 @@ add_executable(Detectors) -target_include_directories(Detectors PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_include_directories(Detectors PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/.. + ${CMAKE_CURRENT_SOURCE_DIR}/../vendor/yaml-cpp/include/ +) target_sources(Detectors PUBLIC AnasenDeadChannelMap.cpp diff --git a/src/Detectors/DetectorApp.cpp b/src/Detectors/DetectorApp.cpp index 83e3f74..351a8e7 100644 --- a/src/Detectors/DetectorApp.cpp +++ b/src/Detectors/DetectorApp.cpp @@ -2,6 +2,8 @@ #include #include +#include "yaml-cpp/yaml.h" + DetectorApp::DetectorApp() : m_resources(nullptr) { @@ -16,38 +18,31 @@ DetectorApp::~DetectorApp() bool DetectorApp::LoadConfig(const std::string& filename) { std::cout<<"----------Detector Efficiency Calculation----------"<(); + m_outputFileName = data["OutputDataFile"].as(); + m_deadChannelFileName = data["DeadChannelFile"].as(); + m_nthreads = data["NumberOfThreads"].as(); + ArrayType type = StringToArrayType(data["ArrayType"].as()); - 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; iSetDeadChannelMap(m_deadChannelFileName); } - std::cout << "Done" << std::endl; - - std::cout << "Allocating " << m_nthreads << " threads..." << std::endl; m_resources = std::make_unique>(m_nthreads); - std::cout << "Done" << std::endl; - std::cout << "Opening input data file " << m_inputFileName << "..." << std::endl; m_fileReader.Open(m_inputFileName, "SimTree"); if(!m_fileReader.IsOpen() || !m_fileReader.IsTree()) { @@ -55,7 +50,6 @@ 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; @@ -64,17 +58,17 @@ bool DetectorApp::LoadConfig(const std::string& filename) for(uint64_t i=1; i()); + YAML::Node nuclei = yamlStream["Reactants"]; + for(auto nucleus : nuclei) + { + params.Z.push_back(nucleus["Z"].as()); + params.A.push_back(nucleus["A"].as()); + } + params.phiMin = yamlStream["PhiMin(deg)"].as(); + params.phiMax = yamlStream["PhiMax(deg)"].as(); + params.meanResidualEx = yamlStream["ResidualExcitationMean(MeV)"].as(); + params.sigmaResidualEx = yamlStream["ResidualExcitationSigma(MeV)"].as(); + params.angularDistFile = yamlStream["AngularDistributionFile"].as(); + return params; + } + + static void SerializeReaction(YAML::Emitter& yamlStream, const StepParameters& params) + { + yamlStream << YAML::BeginMap; + yamlStream << YAML::Key << "Type" << RxnTypeToString(params.rxnType); + yamlStream << YAML::Key << "Reactants" << YAML::Value << YAML::BeginSeq; + for (std::size_t i=0; i()); + YAML::Node nuclei = yamlStream["Reactants"]; + for(auto nucleus : nuclei) + { + params.Z.push_back(nucleus["Z"].as()); + params.A.push_back(nucleus["A"].as()); + } + params.meanBeamEnergy = yamlStream["BeamEnergyMean(MeV)"].as(); + params.sigmaBeamEnergy = yamlStream["BeamEnergySigma(MeV)"].as(); + params.thetaType = StringToRxnThetaType(yamlStream["ThetaType"].as()); + params.thetaMin = yamlStream["ThetaMin(deg)"].as(); + params.thetaMax = yamlStream["ThetaMax(deg)"].as(); + params.phiMin = yamlStream["PhiMin(deg)"].as(); + params.phiMax = yamlStream["PhiMax(deg)"].as(); + params.meanResidualEx = yamlStream["ResidualExcitationMean(MeV)"].as(); + params.sigmaResidualEx = yamlStream["ResidualExcitationSigma(MeV)"].as(); + return params; + } + + static void SerializeTarget(YAML::Emitter& yamlStream, const LayeredTarget& target) + { + + yamlStream << YAML::BeginSeq; + for (int i=0; i z, a, s; + double thickness; + + for (auto layer : yamlStream) + { + thickness = layer["Thickness(ug/cm^2)"].as(); + YAML::Node elements = layer["Elements"]; + z.clear(); + a.clear(); + s.clear(); + for(auto element : elements) + { + z.push_back(element["Z"].as()); + a.push_back(element["A"].as()); + s.push_back(element["S"].as()); + } + target.AddLayer(z, a, s, thickness); + } + + return target; + } + + bool ConfigSerializer::SerializeConfig(const std::string& configfile, const AppParameters& params) + { + std::ofstream output(configfile); + if (!output.is_open()) + { + std::cerr << "Could not open output file " << configfile << std::endl; + return false; + } + + YAML::Emitter yamlStream; + yamlStream << YAML::BeginMap; + yamlStream << YAML::Key << "OutputFile" << YAML::Value << params.outputFileName; + yamlStream << YAML::Key << "Threads" << YAML::Value << params.nThreads; + yamlStream << YAML::Key << "ReactionSamples" << YAML::Value << params.nSamples; + yamlStream << YAML::Key << "ReactionChain" << YAML::Value << YAML::BeginSeq; + for (auto& step : params.chainParams) + { + switch(step.rxnType) + { + case RxnType::Reaction: SerializeReaction(yamlStream, step); break; + case RxnType::Decay: SerializeDecay(yamlStream, step); break; + case RxnType::None: + { + std::cerr << "Error serializing config: none type reaction found!" << std::endl; + return false; + } + } + } + yamlStream << YAML::EndSeq; + yamlStream << YAML::Key << "TargetLayers" << YAML::Value; + SerializeTarget(yamlStream, params.target); + + output << yamlStream.c_str(); + output.close(); + return true; + } + + bool ConfigSerializer::DeserializeConfig(const std::string& configfile, AppParameters& params) + { + YAML::Node data; + try + { + data = YAML::LoadFile(configfile); + } + catch(YAML::ParserException& e) + { + std::cerr << "Error reading config " << configfile << ": " << e.what() <(); + params.nThreads = data["Threads"].as(); + params.nSamples = data["ReactionSamples"].as(); + + auto steps = data["ReactionChain"]; + for (auto step : steps) + { + RxnType type = StringToRxnType(step["Type"].as()); + switch(type) + { + case RxnType::Decay: + { + params.chainParams.push_back(DeserializeDecay(step)); + break; + } + case RxnType::Reaction: + { + params.chainParams.push_back(DeserializeReaction(step)); + break; + } + case RxnType::None: + { + std::cerr << "Error deserializing config: None-type reaction found" << std::endl; + return false; + } + } + } + auto layers = data["TargetLayers"]; + params.target = DeserializeTarget(layers); + return true; + } +} \ No newline at end of file diff --git a/src/Mask/ConfigSerializer.h b/src/Mask/ConfigSerializer.h new file mode 100644 index 0000000..7b0ddd2 --- /dev/null +++ b/src/Mask/ConfigSerializer.h @@ -0,0 +1,16 @@ +#ifndef CONFIG_SERIALIZER_H +#define CONFIG_SERIALIZER_H + +#include "MaskApp.h" + +namespace Mask { + + class ConfigSerializer + { + public: + static bool SerializeConfig(const std::string& configfile, const AppParameters& params); + static bool DeserializeConfig(const std::string& configfile, AppParameters& params); + }; +} + +#endif \ No newline at end of file diff --git a/src/Mask/CoupledThreeStepSystem.cpp b/src/Mask/CoupledThreeStepSystem.cpp new file mode 100644 index 0000000..0c3e210 --- /dev/null +++ b/src/Mask/CoupledThreeStepSystem.cpp @@ -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& params) : + ReactionSystem() + { + m_nuclei.resize(8); + Init(params); + } + + CoupledThreeStepSystem::~CoupledThreeStepSystem() {} + + void CoupledThreeStepSystem::Init(const std::vector& 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] << ")" <" + << 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(); + } +} \ No newline at end of file diff --git a/src/Mask/CoupledThreeStepSystem.h b/src/Mask/CoupledThreeStepSystem.h new file mode 100644 index 0000000..fed700e --- /dev/null +++ b/src/Mask/CoupledThreeStepSystem.h @@ -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& params); + ~CoupledThreeStepSystem(); + + virtual void SetLayeredTarget(const LayeredTarget& target) override; + virtual void RunSystem() override; + + protected: + void Init(const std::vector& params); + void SetSystemEquation() override; + CoupledThreeStepParameters SampleParameters(); + void SampleCoupling(CoupledThreeStepParameters& params); + + Reaction m_step1, m_step2, m_step3; + }; + +} + +#endif \ No newline at end of file diff --git a/src/Mask/LayeredTarget.h b/src/Mask/LayeredTarget.h index cf6dcf7..242cbcd 100644 --- a/src/Mask/LayeredTarget.h +++ b/src/Mask/LayeredTarget.h @@ -31,9 +31,9 @@ namespace Mask { double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth); double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth); std::size_t FindLayerContaining(int Z, int A); - std::size_t GetNumberOfLayers() { return m_layers.size(); } + std::size_t GetNumberOfLayers() const { return m_layers.size(); } void SetName(std::string& n) { m_name = n; } - const Target& GetLayerInfo(int index) { return m_layers[index]; } + const Target& GetLayerInfo(int index) const { return m_layers[index]; } const std::string& GetName() { return m_name; } private: diff --git a/src/Mask/MaskApp.cpp b/src/Mask/MaskApp.cpp index 7e76d14..5806ca1 100644 --- a/src/Mask/MaskApp.cpp +++ b/src/Mask/MaskApp.cpp @@ -1,185 +1,80 @@ #include "MaskApp.h" +#include "ConfigSerializer.h" #include #include #include "TFile.h" #include "TTree.h" +#include "yaml-cpp/yaml.h" + namespace Mask { MaskApp::MaskApp() : - m_system(nullptr), m_resources(nullptr) + m_resources(nullptr) { std::cout<<"----------Monte Carlo Simulation of Kinematics----------"<>junk>>m_outputName; - input>>junk>>m_nthreads; - std::getline(input, junk); - std::getline(input, junk); - input>>junk>>m_nsamples; - std::cout<<"Allocating resources... Asking for " << m_nthreads << " threads..."; - m_resources = std::make_unique>(m_nthreads); - std::cout<<" Complete."< params; - int z, a; - while(input>>junk) + //Initialize the system + //Reaction chain + for(uint32_t i=0; i> junk >> junk; - currentParams.rxnType = StringToRxnType(junk); - if(currentParams.rxnType == RxnType::Reaction) - { - input >> junk; - for(int i=0; i<3; i++) - { - input >> z >> a; - currentParams.Z.push_back(z); - currentParams.A.push_back(a); - } - input >> junk; - input >> junk >> currentParams.meanBeamEnergy; - input >> junk >> currentParams.sigmaBeamEnergy; - input >> junk >> junk; - currentParams.thetaType = StringToRxnThetaType(junk); - input >> junk >> currentParams.thetaMin; - input >> junk >> currentParams.thetaMax; - input >> junk >> currentParams.phiMin; - input >> junk >> currentParams.phiMax; - input >> junk >> currentParams.meanResidualEx; - input >> junk >> currentParams.sigmaResidualEx; - params.push_back(currentParams); - } - else if(currentParams.rxnType == RxnType::Decay) - { - input >> junk; - for(int i=0; i<2; i++) - { - input >> z >> a; - currentParams.Z.push_back(z); - currentParams.A.push_back(a); - } - input >> junk; - input >> junk >> currentParams.phiMin; - input >> junk >> currentParams.phiMax; - input >> junk >> currentParams.meanResidualEx; - input >> junk >> currentParams.sigmaResidualEx; - input >> junk >> currentParams.angularDistFile; - params.push_back(currentParams); - } - else - { - std::cerr << "Invalid reaction information at MaskApp::LoadConfig!" << std::endl; - return false; - } - } - } - - m_system = CreateSystem(params); - if(m_system == nullptr || !m_system->IsValid()) - { - std::cerr<<"Failure to parse reaction system... configuration not loaded."<IsValid()) { std::cerr<<"Failure to parse reaction system... configuration not loaded."< avec, zvec, svec; - int s; - input>>junk>>nlayers; - for(int i=0; i>junk>>junk>>thickness; - avec.clear(); zvec.clear(); svec.clear(); - while(input>>junk) - { - if(junk == "begin_elements") - { - input>>junk>>junk>>junk; - continue; - } - else if (junk == "end_elements") - break; - input>>z>>a>>s; - zvec.push_back(z); avec.push_back(a); svec.push_back(s); - } - target.AddLayer(zvec, avec, svec, thickness); - input>>junk; - } - m_system->SetLayeredTarget(target); - + //Link the target for(auto system : m_systemList) { - system->SetLayeredTarget(target); + system->SetLayeredTarget(m_params.target); } + //Setup threading + m_resources = std::make_unique>(m_params.nThreads); + //Little bit of integer division mangling to make sure we do the total number of samples + uint64_t quotient = m_params.nSamples / m_params.nThreads; + uint64_t remainder = m_params.nSamples % m_params.nThreads; + m_chunkSamples.push_back(quotient + remainder); + for(uint64_t i=1; iGetSystemEquation()<GetSystemEquation() << std::endl; + std::cout << "Number of samples: " << m_params.nSamples << std::endl; + std::cout << "Number of threads: " << m_params.nThreads << std::endl; + std::cout << "Outputing data to file: " << m_params.outputFileName << std::endl; return true; } - - //Not implemented... yet! + bool MaskApp::SaveConfig(const std::string& filename) - { - return true; + { + std::cout << "Writing configuration to " << filename << std::endl; + return ConfigSerializer::SerializeConfig(filename, m_params); } - + void MaskApp::Run() { std::cout<<"Running simulation..."<Branch("nuclei", m_system->GetNuclei()); - - //For progress tracking - uint64_t percent5 = 0.05*m_nsamples; - uint64_t count = 0; - uint64_t npercent = 0; - - for(uint64_t i=0; iRunSystem(); - tree->Fill(); - } - - tree->Write(tree->GetName(), TObject::kOverwrite); - output->Close(); - - std::cout< chainParams; + LayeredTarget target; + }; + class MaskApp { public: @@ -23,17 +32,9 @@ namespace Mask { bool SaveConfig(const std::string& filename); void Run(); - void RunSingleThread(); private: - void RunChunk(ReactionSystem* system); - - ReactionSystem* m_system; - - std::string m_outputName; - RxnType m_rxnType; - uint64_t m_nsamples; - uint32_t m_nthreads; + AppParameters m_params; std::vector m_systemList; //One system for each thread std::vector m_chunkSamples; diff --git a/src/Mask/Nucleus.h b/src/Mask/Nucleus.h index 2769cf6..c820051 100644 --- a/src/Mask/Nucleus.h +++ b/src/Mask/Nucleus.h @@ -12,6 +12,7 @@ #include #include #include "Math/Vector4D.h" +#include "Math/Point3D.h" #include "MassLookup.h" namespace Mask { @@ -48,6 +49,8 @@ namespace Mask { double detectedKE = 0.0; double detectedTheta = 0.0; double detectedPhi = 0.0; + + ROOT::Math::XYZPoint detectedPos = ROOT::Math::XYZPoint(0., 0., 0.); }; Nucleus CreateNucleus(uint32_t z, uint32_t a); diff --git a/src/Mask/Reaction.cpp b/src/Mask/Reaction.cpp index 751dfa7..1a1b871 100644 --- a/src/Mask/Reaction.cpp +++ b/src/Mask/Reaction.cpp @@ -229,7 +229,6 @@ namespace Mask { m_residual->vec4 = m_target->vec4 - m_ejectile->vec4; //energy loss for the *light* break up nucleus - double keorig = 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); 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; + } } diff --git a/src/Mask/Reaction.h b/src/Mask/Reaction.h index 12b9e07..f6ea1e3 100644 --- a/src/Mask/Reaction.h +++ b/src/Mask/Reaction.h @@ -42,6 +42,9 @@ namespace Mask { void SetRxnLayer(std::size_t layer) { m_rxnLayer = layer; }; 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; }; std::size_t GetRxnLayer() const { return m_rxnLayer; }; diff --git a/src/Mask/RxnType.h b/src/Mask/RxnType.h index cb46d59..b868758 100644 --- a/src/Mask/RxnType.h +++ b/src/Mask/RxnType.h @@ -37,6 +37,17 @@ namespace Mask { return RxnType::None; } + static std::string RxnTypeToString(RxnType type) + { + switch(type) + { + case RxnType::Decay: return "Decay"; + case RxnType::Reaction: return "Reaction"; + case RxnType::None: return "None"; + default: return "None"; + } + } + static RxnThetaType StringToRxnThetaType(const std::string& type) { if(type == "CenterOfMass") @@ -52,6 +63,17 @@ namespace Mask { return RxnThetaType::None; } } + + static std::string RxnThetaTypeToString(RxnThetaType type) + { + switch(type) + { + case RxnThetaType::CenterOfMass: return "CenterOfMass"; + case RxnThetaType::Lab: return "Lab"; + case RxnThetaType::None: return "None"; + default: return "None"; + } + } } #endif \ No newline at end of file diff --git a/src/Mask/Target.h b/src/Mask/Target.h index be8b4a4..6f3a8b1 100644 --- a/src/Mask/Target.h +++ b/src/Mask/Target.h @@ -31,12 +31,12 @@ namespace Mask { double GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle); double GetEnergyLossFractionalDepth(int zp, int ap, double startEnergy, double angle, double percent_depth); double GetReverseEnergyLossFractionalDepth(int zp, int ap, double finalEnergy, double angle, double percent_depth); - inline const double& GetThickness() { return m_thickness; } - inline int GetNumberOfElements() { return m_Z.size(); } - inline int GetElementZ(int index) { return m_Z[index]; } - inline int GetElementA(int index) { return m_A[index]; } - inline int GetElementStoich(int index) { return m_stoich[index]; } - + const double& GetThickness() const { return m_thickness; } + int GetNumberOfElements() const { return m_Z.size(); } + int GetElementZ(int index) const { return m_Z[index]; } + int GetElementA(int index) const { return m_A[index]; } + int GetElementStoich(int index) const { return m_stoich[index]; } + private: void Init(const std::vector& z, const std::vector& a, const std::vector& stoich); diff --git a/src/Mask/ThreeStepSystem.cpp b/src/Mask/ThreeStepSystem.cpp index d66710e..4d43acb 100644 --- a/src/Mask/ThreeStepSystem.cpp +++ b/src/Mask/ThreeStepSystem.cpp @@ -2,6 +2,11 @@ #include "RandomGenerator.h" #include "KinematicsExceptions.h" +#include "Math/Boost.h" +#include "Math/Vector3D.h" +#include "Math/RotationZ.h" +#include "Math/RotationY.h" + namespace Mask { ThreeStepSystem::ThreeStepSystem(const std::vector& params) : @@ -56,8 +61,8 @@ namespace Mask { 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[5] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3 - m_nuclei[5] = CreateNucleus(zb4, ab4); //breakup4 + 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])); @@ -114,44 +119,58 @@ namespace Mask { << m_nuclei[7].isotopicSymbol; 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() { - //Sample parameters - std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); - double bke = m_beamDistributions[0](gen); - double rxnTheta = std::acos((m_thetaRanges[0])(gen)); - double rxnPhi = m_phiRanges[0](gen); - double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta(); - double decay1Theta = std::acos(decay1costheta); - double decay1Phi = m_phiRanges[1](gen); - 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)); + ThreeStepParameters 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(rxnDepth); - m_step1.SetBeamKE(bke); - m_step1.SetPolarRxnAngle(rxnTheta); - m_step1.SetAzimRxnAngle(rxnPhi); - m_step1.SetExcitation(residEx); + 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(rxnDepth); - m_step2.SetPolarRxnAngle(decay1Theta); - m_step2.SetAzimRxnAngle(decay1Phi); - m_step2.SetExcitation(decay1Ex); + m_step2.SetReactionDepth(params.rxnDepth); + m_step2.SetPolarRxnAngle(params.decay1Theta); + m_step2.SetAzimRxnAngle(params.decay1Phi); + m_step2.SetExcitation(params.decay1Ex); - m_step3.SetReactionDepth(rxnDepth); - m_step3.SetPolarRxnAngle(decay2Theta); - m_step3.SetAzimRxnAngle(decay2Phi); - m_step3.SetExcitation(decay2Ex); + m_step3.SetReactionDepth(params.rxnDepth); + m_step3.SetPolarRxnAngle(params.decay2Theta); + m_step3.SetAzimRxnAngle(params.decay2Phi); + m_step3.SetResidualEnergyLoss(true); + m_step3.SetExcitation(params.decay2Ex); m_step1.Calculate(); - if(decay1costheta == -10) + 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); @@ -161,15 +180,13 @@ namespace Mask { } m_step2.Calculate(); - if(decay2costheta == -10) + if(params.cosdecay2Theta == -10) { 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_step3.SetResidualEnergyLoss(true); m_step3.Calculate(); - } } \ No newline at end of file diff --git a/src/Mask/ThreeStepSystem.h b/src/Mask/ThreeStepSystem.h index 7b933f4..5238d88 100644 --- a/src/Mask/ThreeStepSystem.h +++ b/src/Mask/ThreeStepSystem.h @@ -6,6 +6,23 @@ 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 { public: @@ -18,6 +35,7 @@ namespace Mask { protected: void Init(const std::vector& params); void SetSystemEquation() override; + ThreeStepParameters SampleParameters(); Reaction m_step1, m_step2, m_step3; diff --git a/src/Mask/TwoStepSystem.cpp b/src/Mask/TwoStepSystem.cpp index 691367f..2b0fcc3 100644 --- a/src/Mask/TwoStepSystem.cpp +++ b/src/Mask/TwoStepSystem.cpp @@ -92,35 +92,59 @@ namespace Mask { << m_nuclei[5].isotopicSymbol; 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() { //Sample parameters - std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); - double bke = (m_beamDistributions[0])(gen); - double rxnTheta = std::acos((m_thetaRanges[0])(gen)); - double rxnPhi = (m_phiRanges[0])(gen); - double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta(); - double decay1Theta = std::acos(decay1costheta); - double decay1Phi = m_phiRanges[1](gen); - double residEx = (m_exDistributions[0])(gen); - double decay2Ex = m_exDistributions[1](gen); - double rxnDepth = (m_rxnDepthDist(gen)); + // std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); + // double bke = (m_beamDistributions[0])(gen); + // double rxnTheta = std::acos((m_thetaRanges[0])(gen)); + // double rxnPhi = (m_phiRanges[0])(gen); + // double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta(); + // double decay1Theta = std::acos(decay1costheta); + // double decay1Phi = m_phiRanges[1](gen); + // double residEx = (m_exDistributions[0])(gen); + // double decay2Ex = m_exDistributions[1](gen); + // double rxnDepth = (m_rxnDepthDist(gen)); + + TwoStepParameters params; + do + { + params = SampleParameters(); + } + while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) + && m_step2.CheckDecayThreshold(params.residEx, params.decay2Ex))); - m_step1.SetReactionDepth(rxnDepth); - m_step1.SetBeamKE(bke); - m_step1.SetPolarRxnAngle(rxnTheta); - m_step1.SetAzimRxnAngle(rxnPhi); - m_step1.SetExcitation(residEx); + 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(rxnDepth); - m_step2.SetPolarRxnAngle(decay1Theta); - m_step2.SetAzimRxnAngle(decay1Phi); - m_step2.SetExcitation(decay2Ex); + m_step2.SetReactionDepth(params.rxnDepth); + m_step2.SetPolarRxnAngle(params.decay1Theta); + m_step2.SetAzimRxnAngle(params.decay1Phi); + m_step2.SetExcitation(params.decay2Ex); m_step1.Calculate(); - if(decay1costheta == -10) + 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); diff --git a/src/Mask/TwoStepSystem.h b/src/Mask/TwoStepSystem.h index 70a1773..16d29df 100644 --- a/src/Mask/TwoStepSystem.h +++ b/src/Mask/TwoStepSystem.h @@ -6,6 +6,19 @@ 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 { public: @@ -18,6 +31,7 @@ namespace Mask { private: void Init(const std::vector& params); void SetSystemEquation() override; + TwoStepParameters SampleParameters(); Reaction m_step1, m_step2; }; diff --git a/src/Plotters/RootPlotter.cpp b/src/Plotters/RootPlotter.cpp index 83809e0..a56b8ae 100644 --- a/src/Plotters/RootPlotter.cpp +++ b/src/Plotters/RootPlotter.cpp @@ -1,7 +1,8 @@ #include "RootPlotter.h" #include #include -#include +#include +#include "Math/Boost.h" #include @@ -10,7 +11,8 @@ static double FullPhi(double 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); } @@ -42,10 +44,16 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam std::cout<<"\rPercent of data processed: "<GetEntry(i); - for(Mask::Nucleus& nuc : *(dataHandle)) + // for(Mask::Nucleus& nuc : *(dataHandle)) + // { + // FillData(nuc); + // } + for(int i=0; isize(); i++) { - FillData(nuc); + FillData(dataHandle->at(i), i); } + //Don't leave this in! + //Correlations(*dataHandle); } std::cout<Close(); @@ -57,40 +65,142 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam output->Close(); } -void RootPlotter::FillData(const Mask::Nucleus& nuc) +void RootPlotter::FillData(const Mask::Nucleus& nuc, int i) { - std::string modifier = ""; - if(nuc.isDetected) - modifier = "detected"; + std::string mod = "_detected"; + std::string num = "_" + std::to_string(i); 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_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 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 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 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_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_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(ex_name, ex_title, 260, -1.0, 25, nuc.GetExcitationEnergy()); 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_ph_name, ke_vs_ph_title, 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(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(ke_vs_th_name_det, ke_vs_th_title_det, nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2); + 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_det, th_vs_ph_title_det, nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2); + MyFill(ex_name_det, ex_title_det, 260, -1.0, 25, nuc.GetExcitationEnergy()); + 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& 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) { auto iter = m_map.find(name); diff --git a/src/Plotters/RootPlotter.h b/src/Plotters/RootPlotter.h index 1cc263e..19e4131 100644 --- a/src/Plotters/RootPlotter.h +++ b/src/Plotters/RootPlotter.h @@ -6,6 +6,7 @@ #include #include "Nucleus.h" +#include "Target.h" #include #include @@ -20,7 +21,8 @@ public: void Run(const std::string& inputname, const std::string& outputname); private: - void FillData(const Mask::Nucleus& nuc); + void FillData(const Mask::Nucleus& nuc, int i); + void Correlations(const std::vector& 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 binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy); @@ -29,6 +31,7 @@ private: std::unordered_map> m_map; static constexpr double s_rad2deg = 180.0/M_PI; + Mask::Target m_target; }; #endif \ No newline at end of file diff --git a/src/vendor/yaml-cpp b/src/vendor/yaml-cpp new file mode 160000 index 0000000..0e6e28d --- /dev/null +++ b/src/vendor/yaml-cpp @@ -0,0 +1 @@ +Subproject commit 0e6e28d1a38224fc8172fae0109ea7f673c096db