From 76ae86d5957fb7c40b5b7540fc602ce5973d5247 Mon Sep 17 00:00:00 2001 From: gwm17 Date: Thu, 1 Sep 2022 09:52:09 -0400 Subject: [PATCH] Overhaul input and reaction system construction. Much more streamlined, now can have more sampling parameters (chained decay to excited states, etc.) --- input.txt | 47 +++++++---- src/Mask/AngularDistribution.cpp | 2 +- src/Mask/DecaySystem.cpp | 61 ++++++++------- src/Mask/DecaySystem.h | 21 ++--- src/Mask/MaskApp.cpp | 121 ++++++++++++++++++----------- src/Mask/OneStepSystem.cpp | 71 ++++++++--------- src/Mask/OneStepSystem.h | 14 ++-- src/Mask/Reaction.cpp | 5 +- src/Mask/ReactionSystem.cpp | 96 +++++++++++++++++------ src/Mask/ReactionSystem.h | 82 +++++++++----------- src/Mask/RxnType.h | 46 ++--------- src/Mask/ThreeStepSystem.cpp | 129 +++++++++++++++++++------------ src/Mask/ThreeStepSystem.h | 32 +++----- src/Mask/TwoStepSystem.cpp | 95 +++++++++++++---------- src/Mask/TwoStepSystem.h | 21 +---- src/Plotters/RootPlotter.cpp | 20 ++--- 16 files changed, 469 insertions(+), 394 deletions(-) diff --git a/input.txt b/input.txt index 5ca223b..7922dc8 100644 --- a/input.txt +++ b/input.txt @@ -1,12 +1,38 @@ ----------Data Information---------- OutputFile: /media/data/gwm17/mask_tests/10B3Hea_16800keV_5Lia_74B.root ----------Reaction Information---------- -begin_nuclei(Z,A) - 5 10 - 2 3 - 2 4 - 2 4 -end_nuclei(Z,A) +NumberOfSamples: 100000 +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 @@ -15,13 +41,4 @@ begin_layer element 5 10 1 end_elements end_layer -----------Sampling Information---------- -NumberOfSamples: 1000000 -BeamMeanEnergy(MeV): 24.0 BeamEnergySigma(MeV): 0.0 -EjectileThetaType(0=Lab,1=CM): 0 -EjectileThetaMin(deg): 15.0 EjectileThetaMax(deg): 15.0 -EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 0.0 -ResidualExMean(MeV): 16.8 ResidualExSigma(MeV): 0.023 -Decay1_DistributionFile: ./etc/isotropic_dist.txt -Decay2_DistributionFile: ./etc/isotropic_dist.txt -------------------------------------- diff --git a/src/Mask/AngularDistribution.cpp b/src/Mask/AngularDistribution.cpp index 0884e80..8f07b2b 100644 --- a/src/Mask/AngularDistribution.cpp +++ b/src/Mask/AngularDistribution.cpp @@ -13,7 +13,7 @@ namespace Mask { } AngularDistribution::AngularDistribution(const std::string& file) : - m_branchingRatio(1.0), m_L(0), m_isIsotropic(true) + m_uniformCosineDist(-1.0, 1.0), m_branchingRatio(1.0), m_L(0), m_isIsotropic(true) { ReadDistributionFile(file); } diff --git a/src/Mask/DecaySystem.cpp b/src/Mask/DecaySystem.cpp index 11a562e..7146d9b 100644 --- a/src/Mask/DecaySystem.cpp +++ b/src/Mask/DecaySystem.cpp @@ -5,49 +5,52 @@ namespace Mask { - DecaySystem::DecaySystem() : + DecaySystem::DecaySystem(const std::vector& params) : ReactionSystem() { m_nuclei.resize(3); - } - - DecaySystem::DecaySystem(const std::vector& z, const std::vector& a) : - ReactionSystem() - { - m_nuclei.resize(3); - SetNuclei(z, a); + Init(params); } DecaySystem::~DecaySystem() {} - bool DecaySystem::SetNuclei(const std::vector& z, const std::vector& a) + void DecaySystem::Init(const std::vector& params) { - if(z.size() != a.size() || z.size() != 2) - return false; + if(params.size() != 1 || params[0].rxnType != RxnType::Decay || + params[0].Z.size() != 2 || params[0].A.size() != 2) + { + m_isValid = false; + std::cerr << "Invalid parameters at DecaySystem::Init(), does not match Decay signature!" << std::endl; + return; + } - int zr = z[0] - z[1]; - int ar = a[0] - a[1]; + const StepParameters& step1Params = params[0]; - m_nuclei[0] = CreateNucleus(z[0], a[0]); //target - m_nuclei[1] = CreateNucleus(z[1], a[1]); //breakup1 + int zr = step1Params.Z[0] - step1Params.Z[1]; + int ar = step1Params.A[0] - step1Params.A[1]; + + m_nuclei[0] = CreateNucleus(step1Params.Z[0], step1Params.A[0]); //target + m_nuclei[1] = CreateNucleus(step1Params.Z[1], step1Params.A[1]); //breakup1 m_nuclei[2] = CreateNucleus(zr, ar); //breakup2 m_step1.BindNuclei(&(m_nuclei[0]), nullptr, &(m_nuclei[1]), &(m_nuclei[2])); SetSystemEquation(); - return true; - } - std::vector* DecaySystem::GetNuclei() - { - return &m_nuclei; + + AddPhiRange(step1Params.phiMin, step1Params.phiMax); + AddDecayAngularDistribution(step1Params.angularDistFile); + AddExcitationDistribution(step1Params.meanResidualEx, step1Params.sigmaResidualEx); + + return; } - void DecaySystem::LinkTarget() { - m_step1.SetLayeredTarget(&m_target); - + void DecaySystem::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_step1.SetRxnLayer(m_rxnLayer); m_isTargetSet = true; } @@ -66,14 +69,14 @@ namespace Mask { void DecaySystem::RunSystem() { - //Link up the target if it hasn't been done yet - if(!m_isTargetSet) - LinkTarget(); - - double rxnTheta = std::acos(m_step1Distribution.GetRandomCosTheta()); - double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator()); + std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); + double rxnTheta = m_decayAngularDistributions[0].GetRandomCosTheta(); + double rxnPhi = m_phiRanges[0](gen); + double ex = m_exDistributions[0](gen); + m_step1.SetPolarRxnAngle(rxnTheta); m_step1.SetAzimRxnAngle(rxnPhi); + m_step1.SetExcitation(ex); m_step1.SetResidualEnergyLoss(true); m_step1.Calculate(); } diff --git a/src/Mask/DecaySystem.h b/src/Mask/DecaySystem.h index c0a1122..23b6649 100644 --- a/src/Mask/DecaySystem.h +++ b/src/Mask/DecaySystem.h @@ -6,29 +6,20 @@ namespace Mask { - class DecaySystem: public ReactionSystem { + class DecaySystem: public ReactionSystem + { public: - DecaySystem(); - DecaySystem(const std::vector& z, const std::vector& a); + DecaySystem(const std::vector& params); ~DecaySystem(); - bool SetNuclei(const std::vector& z, const std::vector& a) override; - void RunSystem() override; - std::vector*GetNuclei() override; - - virtual void SetDecay1Distribution(const std::string& filename) override { m_step1Distribution.ReadDistributionFile(filename); } - - int GetDecay1AngularMomentum() { return m_step1Distribution.GetL(); } - double GetDecay1BranchingRatio() { return m_step1Distribution.GetBranchingRatio(); } + virtual void SetLayeredTarget(const LayeredTarget& target) override; + virtual void RunSystem() override; private: - void LinkTarget() override; + void Init(const std::vector& params); void SetSystemEquation() override; Reaction m_step1; - - AngularDistribution m_step1Distribution; - }; } diff --git a/src/Mask/MaskApp.cpp b/src/Mask/MaskApp.cpp index 07a92e8..bc27127 100644 --- a/src/Mask/MaskApp.cpp +++ b/src/Mask/MaskApp.cpp @@ -10,7 +10,7 @@ namespace Mask { MaskApp::MaskApp() : m_system(nullptr) { - std::cout<<"----------GWM Kinematics Simulation----------"<>junk>>m_outputName; - - std::vector avec, zvec, svec; - int z, a, s; - getline(input, junk); - getline(input, junk); + std::getline(input, junk); + std::getline(input, junk); + input>>junk>>m_nsamples; + std::vector params; + int z, a; while(input>>junk) { - if(junk == "begin_nuclei(Z,A)") + if(junk == "begin_chain") continue; - else if (junk == "end_nuclei(Z,A)") + else if (junk == "end_chain") break; - else + else if(junk == "begin_step") { - z = std::stoi(junk); - input>>a; - zvec.push_back(z); - avec.push_back(a); + 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(zvec, avec); - if(m_system == nullptr) + m_system = CreateSystem(params); + if(m_system == nullptr || !m_system->IsValid()) { + std::cerr<<"Param size: "< avec, zvec, svec; + int s; input>>junk>>nlayers; for(int i=0; i>z>>a>>s; zvec.push_back(z); avec.push_back(a); svec.push_back(s); } - m_system->AddTargetLayer(zvec, avec, svec, thickness); + target.AddLayer(zvec, avec, svec, thickness); input>>junk; } + m_system->SetLayeredTarget(target); + + std::cout<<"Outputing data to file: "<>junk>>m_nsamples; - input>>junk>>par1>>junk>>par2; - m_system->SetBeamDistro(par1, par2); - input>>junk>>thetaTypeString; - m_system->SetReactionThetaType(StringToRxnThetaType(thetaTypeString)); - - input>>junk>>par1>>junk>>par2; - m_system->SetTheta1Range(par1, par2); - input>>junk>>par1>>junk>>par2; - m_system->SetPhi1Range(par1, par2); - input>>junk>>par1>>junk>>par2; - m_system->SetExcitationDistro(par1, par2); - input>>junk>>dfile1; - input>>junk>>dfile2; - m_system->SetDecay1Distribution(dfile1); - m_system->SetDecay2Distribution(dfile2); - std::cout<<"Number of samples: "<& params) : ReactionSystem() { m_nuclei.resize(4); - } - - OneStepSystem::OneStepSystem(const std::vector& z, const std::vector& a) : - ReactionSystem() - { - m_nuclei.resize(4); - SetNuclei(z, a); + Init(params); } OneStepSystem::~OneStepSystem() {} - bool OneStepSystem::SetNuclei(const std::vector& z, const std::vector& a) + void OneStepSystem::Init(const std::vector& params) { - if(z.size() != a.size() || z.size() != 3) - return false; + if(params.size() != 1 || params[0].rxnType != RxnType::Reaction || + params[0].Z.size() != 3 || params[0].A.size() != 3) + { + m_isValid = false; + std::cerr << "Invalid parameters at OneStepSystem::Init(), does not match OneStep signature!" << std::endl; + return; + } - int zr = z[0] + z[1] - z[2]; - int ar = a[0] + a[1] - a[2]; + const StepParameters& step1Params = params[0]; - m_nuclei[0] = CreateNucleus(z[0], a[0]); //target - m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile - m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile + //Set nuclei + + int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2]; + int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2]; + + 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_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3])); SetSystemEquation(); - return true; - } - std::vector* OneStepSystem::GetNuclei() - { - return &m_nuclei; + //Set 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); } - void OneStepSystem::LinkTarget() + void OneStepSystem::SetLayeredTarget(const LayeredTarget& target) { - m_step1.SetLayeredTarget(&m_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_step1.SetRxnLayer(m_rxnLayer); m_isTargetSet = true; } @@ -67,18 +72,14 @@ namespace Mask { m_sysEquation = stream.str(); } - void OneStepSystem::RunSystem() { - //Link up the target if it hasn't been done yet - if(!m_isTargetSet) - { - LinkTarget(); - } - + void OneStepSystem::RunSystem() + { //Sample parameters - double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator()); - double rxnTheta = std::acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator())); - double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator()); - double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator()); + 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 residEx = (m_exDistributions[0])(gen); m_step1.SetBeamKE(bke); m_step1.SetPolarRxnAngle(rxnTheta); diff --git a/src/Mask/OneStepSystem.h b/src/Mask/OneStepSystem.h index fe2248e..4062fe8 100644 --- a/src/Mask/OneStepSystem.h +++ b/src/Mask/OneStepSystem.h @@ -5,24 +5,20 @@ namespace Mask { - class OneStepSystem: public ReactionSystem { + class OneStepSystem: public ReactionSystem + { public: - OneStepSystem(); - OneStepSystem(const std::vector& z, const std::vector& a); + OneStepSystem(const std::vector& params); ~OneStepSystem(); - bool SetNuclei(const std::vector& z, const std::vector& a) override; + virtual void SetLayeredTarget(const LayeredTarget& target) override; void RunSystem() override; - std::vector* GetNuclei() override; - - virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); }; private: - void LinkTarget() override; + void Init(const std::vector& params); void SetSystemEquation() override; Reaction m_step1; - }; } diff --git a/src/Mask/Reaction.cpp b/src/Mask/Reaction.cpp index fbc1d07..04f2809 100644 --- a/src/Mask/Reaction.cpp +++ b/src/Mask/Reaction.cpp @@ -187,14 +187,15 @@ namespace Mask { //Calculate in CM, where decay is isotropic void Reaction::CalculateDecay() { - double Q = m_target->vec4.M() - m_ejectile->groundStateMass - m_residual->groundStateMass; + double residualMass = m_residual->groundStateMass + m_ex; + double Q = m_target->vec4.M() - m_ejectile->groundStateMass - residualMass; if(Q < 0) throw QValueException(); ROOT::Math::Boost boost(m_target->vec4.BoostToCM()); m_target->vec4 = boost*m_target->vec4; double ejectE_cm = (m_ejectile->groundStateMass*m_ejectile->groundStateMass - - m_residual->groundStateMass*m_residual->groundStateMass + m_target->vec4.E()*m_target->vec4.E()) / + residualMass*residualMass + m_target->vec4.E()*m_target->vec4.E()) / (2.0*m_target->vec4.E()); double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - m_ejectile->groundStateMass*m_ejectile->groundStateMass); diff --git a/src/Mask/ReactionSystem.cpp b/src/Mask/ReactionSystem.cpp index f8e9cb9..3a8f86d 100644 --- a/src/Mask/ReactionSystem.cpp +++ b/src/Mask/ReactionSystem.cpp @@ -10,41 +10,91 @@ namespace Mask { ReactionSystem::ReactionSystem() : - m_beamDist(nullptr), m_theta1Range(nullptr), m_phi1Range(nullptr), m_exDist(nullptr), m_isTargetSet(false), - m_rxnLayer(0), m_sysEquation("") + m_isTargetSet(false), m_isValid(true), m_rxnLayer(0), m_sysEquation("") { } ReactionSystem::~ReactionSystem() { - delete m_beamDist; - delete m_theta1Range; - delete m_phi1Range; - delete m_exDist; - } - - void ReactionSystem::AddTargetLayer(const std::vector& zt, const std::vector& at, const std::vector& stoich, double thickness) - { - m_target.AddLayer(zt, at, stoich, thickness); } - ReactionSystem* CreateSystem(const std::vector& z, const std::vector& a) + ReactionSystem* CreateSystem(const std::vector& params) { - if(z.size() != a.size()) + switch(params.size()) { - std::cerr<<"Size of Z list does not equal size of A list!"< mean: " << mean + << " sigma: " << sigma << std::endl; + return; + } + m_beamDistributions.emplace_back(mean, sigma); + } + + //Again only really one of these, but same issue as beam. This is the sampling of the primary reaction/decay (step 1) + void ReactionSystem::AddThetaRange(double min, double max) + { + if(min == -1.0 || max == -1.0 || min > max) + { + m_isValid = false; + std::cerr << "Invalid parameters at ReactionSystem::AddThetaRange(), range is invalid -> min: " << min + << " max: " << max << std::endl; + return; + } + m_thetaRanges.emplace_back(std::cos(min*s_deg2rad), std::cos(max*s_deg2rad)); + } + + void ReactionSystem::AddPhiRange(double min, double max) + { + if(min == -1.0 || max == -1.0 || min > max) + { + m_isValid = false; + std::cerr << "Invalid parameters at ReactionSystem::AddPhiRange(), range is invalid -> min: " << min + << " max: " << max << std::endl; + return; + } + m_phiRanges.emplace_back(min*s_deg2rad, max*s_deg2rad); + } + + //Each reaction step can generate an excited nucleus (for a reaction can make an excited residual, decay can make an excited + //breakup2 or "heavy") + void ReactionSystem::AddExcitationDistribution(double mean, double sigma) + { + if(mean == -1.0 || sigma == -1.0) + { + m_isValid = false; + std::cerr << "Invalid parameters at ReactionSystem::AddExcitationDistribution(), distribution is invalid -> mean: " << mean + << " sigma: " << sigma << std::endl; + return; + } + m_exDistributions.emplace_back(mean, sigma); + } + + //Apply angular distributions to decay products. Here apply them to the breakup1 or "light" fragment + void ReactionSystem::AddDecayAngularDistribution(const std::string& filename) + { + if(filename.empty()) + m_decayAngularDistributions.emplace_back(); + else + m_decayAngularDistributions.emplace_back(filename); + } } \ No newline at end of file diff --git a/src/Mask/ReactionSystem.h b/src/Mask/ReactionSystem.h index 02f9fbf..c4761be 100644 --- a/src/Mask/ReactionSystem.h +++ b/src/Mask/ReactionSystem.h @@ -12,75 +12,69 @@ #include "Reaction.h" #include "KinematicsExceptions.h" #include "RxnType.h" +#include "AngularDistribution.h" #include #include namespace Mask { + struct StepParameters + { + RxnType rxnType = RxnType::None; + std::vector Z; + std::vector A; + double meanBeamEnergy = -1.0; + double sigmaBeamEnergy = -1.0; + RxnThetaType thetaType = RxnThetaType::None; + double thetaMin = -1.0; + double thetaMax = -1.0; + double phiMin = -1.0; + double phiMax = -1.0; + double meanResidualEx = -1.0; + double sigmaResidualEx = -1.0; + std::string angularDistFile; + }; + class ReactionSystem { public: ReactionSystem(); virtual ~ReactionSystem(); - - virtual bool SetNuclei(const std::vector& z, const std::vector& a) = 0; - virtual void RunSystem() = 0; - virtual std::vector* GetNuclei() = 0; - virtual void SetReactionThetaType(RxnThetaType type) {} - virtual void SetDecay1Distribution(const std::string& filename) {} - virtual void SetDecay2Distribution(const std::string& filename) {} - - void AddTargetLayer(const std::vector& zt, const std::vector& at, const std::vector& stoich, double thickness); - - /*Set sampling parameters*/ - void SetBeamDistro(double mean, double sigma) - { - if(m_beamDist) - delete m_beamDist; - m_beamDist = new std::normal_distribution(mean, sigma); - } - - void SetTheta1Range(double min, double max) - { - if(m_theta1Range) - delete m_theta1Range; - m_theta1Range = new std::uniform_real_distribution(std::cos(min*s_deg2rad), std::cos(max*s_deg2rad)); - } - - void SetPhi1Range(double min, double max) - { - if(m_phi1Range) - delete m_phi1Range; - m_phi1Range = new std::uniform_real_distribution(min*s_deg2rad, max*s_deg2rad); - } - - void SetExcitationDistro(double mean, double sigma) - { - if(m_exDist) - delete m_exDist; - m_exDist = new std::normal_distribution(mean, sigma); - } + virtual void SetLayeredTarget(const LayeredTarget& target) = 0; + virtual void RunSystem() = 0; + + std::vector* GetNuclei() { return &m_nuclei; } const std::string& GetSystemEquation() const { return m_sysEquation; } - + bool IsValid() const { return m_isValid; } + protected: - virtual void LinkTarget() = 0; virtual void SetSystemEquation() = 0; + + void AddBeamDistribution(double mean, double sigma); + void AddThetaRange(double min, double max); + void AddPhiRange(double min, double max); + void AddExcitationDistribution(double mean, double sigma); + void AddDecayAngularDistribution(const std::string& filename); LayeredTarget m_target; //Sampling information - std::normal_distribution *m_beamDist, *m_exDist; - std::uniform_real_distribution *m_theta1Range, *m_phi1Range; - + std::vector> m_beamDistributions, m_exDistributions; + std::vector> m_thetaRanges, m_phiRanges; + std::vector m_decayAngularDistributions; + bool m_isTargetSet; + bool m_isValid; + std::size_t m_rxnLayer; std::string m_sysEquation; std::vector m_nuclei; + static constexpr double s_deg2rad = M_PI/180.0; }; - ReactionSystem* CreateSystem(const std::vector& z, const std::vector& a); + ReactionSystem* CreateSystem(const std::vector& params); } #endif \ No newline at end of file diff --git a/src/Mask/RxnType.h b/src/Mask/RxnType.h index 9b4f841..cb46d59 100644 --- a/src/Mask/RxnType.h +++ b/src/Mask/RxnType.h @@ -7,11 +7,9 @@ namespace Mask { enum class RxnType { - PureDecay=0, - OneStepRxn=1, - TwoStepRxn=2, - ThreeStepRxn=3, - None=4 + Decay, + Reaction, + None }; enum class RxnThetaType @@ -29,44 +27,16 @@ namespace Mask { ThreeStepSize = 5 }; - static RxnType GetRxnTypeFromString(const std::string& type_str) + static RxnType StringToRxnType(const std::string& type) { - if (type_str == "PureDecay") - return RxnType::PureDecay; - else if (type_str == "OneStepRxn") - return RxnType::OneStepRxn; - else if (type_str == "TwoStepRxn") - return RxnType::TwoStepRxn; - else if (type_str == "ThreeStepRxn") - return RxnType::ThreeStepRxn; + if (type == "Decay") + return RxnType::Decay; + else if (type == "Reaction") + return RxnType::Reaction; else return RxnType::None; } - static std::string GetStringFromRxnType(RxnType type) - { - switch(type) - { - case RxnType::PureDecay: return "PureDecay"; - case RxnType::OneStepRxn: return "OneStepRxn"; - case RxnType::TwoStepRxn: return "TwoStepRxn"; - case RxnType::ThreeStepRxn: return "ThreeStepRxn"; - case RxnType::None : return "None"; - } - - return "None"; - } - - static RxnType GetRxnTypeFromInt(uint32_t value) - { - return static_cast(value); - } - - static uint32_t GetIntFromRxnType(RxnType type) - { - return static_cast(type); - } - static RxnThetaType StringToRxnThetaType(const std::string& type) { if(type == "CenterOfMass") diff --git a/src/Mask/ThreeStepSystem.cpp b/src/Mask/ThreeStepSystem.cpp index 800ae3e..880de24 100644 --- a/src/Mask/ThreeStepSystem.cpp +++ b/src/Mask/ThreeStepSystem.cpp @@ -3,41 +3,60 @@ #include "KinematicsExceptions.h" namespace Mask { - - ThreeStepSystem::ThreeStepSystem() : - ReactionSystem(), m_phi2Range(0, 2.0*M_PI) - { - m_nuclei.resize(8); - } - ThreeStepSystem::ThreeStepSystem(const std::vector& z, const std::vector& a) : - ReactionSystem(), m_phi2Range(0, 2.0*M_PI) + ThreeStepSystem::ThreeStepSystem(const std::vector& params) : + ReactionSystem() { m_nuclei.resize(8); - SetNuclei(z, a); + Init(params); } ThreeStepSystem::~ThreeStepSystem() {} - bool ThreeStepSystem::SetNuclei(const std::vector& z, const std::vector& a) + void ThreeStepSystem::Init(const std::vector& params) { - if(z.size() != a.size() || z.size() < 5) - return false; + 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; + } - int zr = z[0] + z[1] - z[2]; - int zb2 = zr - z[3]; - int zb4 = zb2 - z[4]; - int ar = a[0] + a[1] - a[2]; - int ab2 = ar - a[3]; - int ab4 = ab2 - a[4]; + const StepParameters& step1Params = params[0]; + const StepParameters& step2Params = params[1]; + const StepParameters& step3Params = params[2]; - m_nuclei[0] = CreateNucleus(z[0], a[0]); //target - m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile - m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile + //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] << ")" <* ThreeStepSystem::GetNuclei() - { - return &m_nuclei; + + //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 ThreeStepSystem::LinkTarget() + void ThreeStepSystem::SetLayeredTarget(const LayeredTarget& target) { - m_step1.SetLayeredTarget(&m_target); - m_step2.SetLayeredTarget(&m_target); - m_step3.SetLayeredTarget(&m_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 + } + else throw ReactionLayerException(); } @@ -84,22 +115,22 @@ namespace Mask { m_sysEquation = stream.str(); } - void ThreeStepSystem::RunSystem() { - //Link up the target if it hasn't been done yet - if(!m_isTargetSet) - LinkTarget(); - + void ThreeStepSystem::RunSystem() + { //Sample parameters - double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator()); - double rxnTheta = acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator())); - double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator()); - double decay1costheta = m_step2Distribution.GetRandomCosTheta(); + 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_phi2Range(RandomGenerator::GetInstance().GetGenerator()); - double decay2costheta = m_step3Distribution.GetRandomCosTheta(); + double decay1Phi = m_phiRanges[1](gen); + double decay2costheta = m_decayAngularDistributions[1].GetRandomCosTheta(); double decay2Theta = std::acos(decay2costheta); - double decay2Phi = m_phi2Range(RandomGenerator::GetInstance().GetGenerator()); - double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator()); + double decay2Phi = m_phiRanges[2](gen); + double residEx = m_exDistributions[0](gen); + double decay1Ex = m_exDistributions[1](gen); + double decay2Ex = m_exDistributions[2](gen); m_step1.SetBeamKE(bke); m_step1.SetPolarRxnAngle(rxnTheta); @@ -108,9 +139,11 @@ namespace Mask { m_step2.SetPolarRxnAngle(decay1Theta); m_step2.SetAzimRxnAngle(decay1Phi); + m_step2.SetExcitation(decay1Ex); m_step3.SetPolarRxnAngle(decay2Theta); m_step3.SetAzimRxnAngle(decay2Phi); + m_step3.SetExcitation(decay2Ex); m_step1.Calculate(); diff --git a/src/Mask/ThreeStepSystem.h b/src/Mask/ThreeStepSystem.h index a79017e..7b933f4 100644 --- a/src/Mask/ThreeStepSystem.h +++ b/src/Mask/ThreeStepSystem.h @@ -6,35 +6,21 @@ namespace Mask { - class ThreeStepSystem : public ReactionSystem { + class ThreeStepSystem : public ReactionSystem + { public: - ThreeStepSystem(); - ThreeStepSystem(const std::vector& z, const std::vector& a); + ThreeStepSystem(const std::vector& params); ~ThreeStepSystem(); - bool SetNuclei(const std::vector& z, const std::vector& a) override; - void RunSystem() override; - std::vector* GetNuclei() override; - - virtual void SetDecay1Distribution(const std::string& filename) override { m_step2Distribution.ReadDistributionFile(filename); }; - virtual void SetDecay2Distribution(const std::string& filename) override { m_step3Distribution.ReadDistributionFile(filename); }; - - virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); }; - - int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); }; - int GetDecay2AngularMomentum(){ return m_step3Distribution.GetL(); }; - double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); }; - double GetDecay2BranchingRatio(){ return m_step3Distribution.GetBranchingRatio(); }; - - protected: - void LinkTarget() override; - void SetSystemEquation() override; - std::uniform_real_distribution m_phi2Range; + virtual void SetLayeredTarget(const LayeredTarget& target) override; + virtual void RunSystem() override; + + protected: + void Init(const std::vector& params); + void SetSystemEquation() override; Reaction m_step1, m_step2, m_step3; - AngularDistribution m_step2Distribution, m_step3Distribution; - }; } diff --git a/src/Mask/TwoStepSystem.cpp b/src/Mask/TwoStepSystem.cpp index 9750fd6..2145e2e 100644 --- a/src/Mask/TwoStepSystem.cpp +++ b/src/Mask/TwoStepSystem.cpp @@ -5,58 +5,74 @@ #include namespace Mask { - - TwoStepSystem::TwoStepSystem() : - ReactionSystem(), m_phi2Range(0, 2.0*M_PI) - { - m_nuclei.resize(6); - } - TwoStepSystem::TwoStepSystem(const std::vector& z, const std::vector& a) : - ReactionSystem(), m_phi2Range(0, 2.0*M_PI) + TwoStepSystem::TwoStepSystem(const std::vector& params) : + ReactionSystem() { m_nuclei.resize(6); - SetNuclei(z, a); + Init(params); } TwoStepSystem::~TwoStepSystem() {} - bool TwoStepSystem::SetNuclei(const std::vector&z, const std::vector& a) + void TwoStepSystem::Init(const std::vector& params) { - if(z.size() != a.size() || z.size() != 4) - return false; + if(params.size() != 2 || params[0].rxnType != RxnType::Reaction || params[1].rxnType != RxnType::Decay || + params[0].Z.size() != 3 || params[0].A.size() != 3 || params[1].Z.size() != 2 || params[1].A.size() != 2) + { + m_isValid = false; + std::cerr << "Invalid parameters at TwoStepSystem::Init(), does not match TwoStep signature!" << std::endl; + return; + } - int zr = z[0] + z[1] - z[2]; - int ar = a[0] + a[1] - a[2]; - int zb = zr - z[3]; - int ab = ar - a[3]; + const StepParameters& step1Params = params[0]; + const StepParameters& step2Params = params[1]; - m_nuclei[0] = CreateNucleus(z[0], a[0]); //target - m_nuclei[1] = CreateNucleus(z[1], a[1]); //projectile - m_nuclei[2] = CreateNucleus(z[2], a[2]); //ejectile + //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 TwoStepSystem::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] << ")" <* TwoStepSystem::GetNuclei() - { - return &m_nuclei; + //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); } - void TwoStepSystem::LinkTarget() + void TwoStepSystem::SetLayeredTarget(const LayeredTarget& target) { - m_step1.SetLayeredTarget(&m_target); - m_step2.SetLayeredTarget(&m_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_step1.SetRxnLayer(m_rxnLayer); m_step2.SetRxnLayer(m_rxnLayer); m_isTargetSet = true; @@ -79,18 +95,16 @@ namespace Mask { void TwoStepSystem::RunSystem() { - //Link up the target if it hasn't been done yet - if(!m_isTargetSet) - LinkTarget(); - //Sample parameters - double bke = (*m_beamDist)(RandomGenerator::GetInstance().GetGenerator()); - double rxnTheta = acos((*m_theta1Range)(RandomGenerator::GetInstance().GetGenerator())); - double rxnPhi = (*m_phi1Range)(RandomGenerator::GetInstance().GetGenerator()); - double decay1costheta = m_step2Distribution.GetRandomCosTheta(); + 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_phi2Range(RandomGenerator::GetInstance().GetGenerator()); - double residEx = (*m_exDist)(RandomGenerator::GetInstance().GetGenerator()); + double decay1Phi = m_phiRanges[1](gen); + double residEx = (m_exDistributions[0])(gen); + double decay2Ex = m_exDistributions[1](gen); m_step1.SetBeamKE(bke); m_step1.SetPolarRxnAngle(rxnTheta); @@ -99,6 +113,7 @@ namespace Mask { m_step2.SetPolarRxnAngle(decay1Theta); m_step2.SetAzimRxnAngle(decay1Phi); + m_step2.SetExcitation(decay2Ex); m_step1.Calculate(); diff --git a/src/Mask/TwoStepSystem.h b/src/Mask/TwoStepSystem.h index 9b1249b..70a1773 100644 --- a/src/Mask/TwoStepSystem.h +++ b/src/Mask/TwoStepSystem.h @@ -9,30 +9,17 @@ namespace Mask { class TwoStepSystem : public ReactionSystem { public: - TwoStepSystem(); - TwoStepSystem(const std::vector& z, const std::vector& a); + TwoStepSystem(const std::vector& params); ~TwoStepSystem(); - bool SetNuclei(const std::vector& z, const std::vector& a) override; - void RunSystem() override; - std::vector* GetNuclei() override; - - virtual void SetDecay1Distribution(const std::string& filename) override { m_step2Distribution.ReadDistributionFile(filename); }; - - virtual void SetReactionThetaType(RxnThetaType type) override { m_step1.SetEjectileThetaType(type); }; - int GetDecay1AngularMomentum() { return m_step2Distribution.GetL(); }; - double GetDecay1BranchingRatio() { return m_step2Distribution.GetBranchingRatio(); }; + virtual void SetLayeredTarget(const LayeredTarget& target) override; + virtual void RunSystem() override; private: - void LinkTarget() override; + void Init(const std::vector& params); void SetSystemEquation() override; - std::uniform_real_distribution m_phi2Range; - Reaction m_step1, m_step2; - - AngularDistribution m_step2Distribution; - }; } diff --git a/src/Plotters/RootPlotter.cpp b/src/Plotters/RootPlotter.cpp index 9a07dc9..83809e0 100644 --- a/src/Plotters/RootPlotter.cpp +++ b/src/Plotters/RootPlotter.cpp @@ -75,18 +75,18 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc) std::string angdist_name = sym + modifier +"_angDist"; std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts"; - MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2); - MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.GetKE(), 4); - MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2); - MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy()); - MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM)); + 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(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.vec4.Theta()*s_rad2deg, nuc.detectedKE, 2); - MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), FullPhi(nuc.vec4.Phi())*s_rad2deg, nuc.detectedKE, 4); - MyFill(th_vs_ph_name.c_str(), th_vs_ph_title.c_str(), nuc.vec4.Theta()*s_rad2deg, FullPhi(nuc.vec4.Phi())*s_rad2deg, 2); - MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy()); - MyFill(angdist_name.c_str(), angdist_title.c_str(),20,-1.0,1.0,std::cos(nuc.thetaCM)); + 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)); } }