diff --git a/src/Detectors/SabreArray.h b/src/Detectors/SabreArray.h index b32ca6c..24cac47 100644 --- a/src/Detectors/SabreArray.h +++ b/src/Detectors/SabreArray.h @@ -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_degraderThickness = 70.0 * 1.0e-4 * 16.69 * 1e6; //tantalum degrader (70 um thick) - static constexpr double s_energyThreshold = 0.25; //in MeV + static constexpr double s_energyThreshold = 0.2; //in MeV }; diff --git a/src/Mask/DecaySystem.cpp b/src/Mask/DecaySystem.cpp index 7146d9b..1c18157 100644 --- a/src/Mask/DecaySystem.cpp +++ b/src/Mask/DecaySystem.cpp @@ -73,7 +73,9 @@ namespace Mask { double rxnTheta = m_decayAngularDistributions[0].GetRandomCosTheta(); double rxnPhi = m_phiRanges[0](gen); double ex = m_exDistributions[0](gen); + double rxnDepth = m_rxnDepthDist(gen); + m_step1.SetReactionDepth(rxnDepth); m_step1.SetPolarRxnAngle(rxnTheta); m_step1.SetAzimRxnAngle(rxnPhi); m_step1.SetExcitation(ex); diff --git a/src/Mask/LayeredTarget.cpp b/src/Mask/LayeredTarget.cpp index 184f49c..3d4609a 100644 --- a/src/Mask/LayeredTarget.cpp +++ b/src/Mask/LayeredTarget.cpp @@ -16,7 +16,7 @@ Written by G.W. McCann Aug. 2020 namespace Mask { LayeredTarget::LayeredTarget() : - m_name(""), m_fractionalDepthDistribution(0.0, 1.0) + m_name("") { } @@ -33,7 +33,7 @@ namespace Mask { Calculates energy loss assuming that the reaction occurs in the middle of the target layer Note that the layer order can matter! */ - double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, std::size_t rxnLayer, double angle) + double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth) { if(rxnLayer > m_layers.size()) { @@ -48,8 +48,7 @@ namespace Mask { { if(i == rxnLayer) { - frac = m_fractionalDepthDistribution(RandomGenerator::GetInstance().GetGenerator()); - eloss += m_layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, frac); + eloss += m_layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, rxnDepth); newEnergy = startEnergy - eloss; } else @@ -67,7 +66,7 @@ namespace Mask { Calculates energy loss assuming that the reaction occurs in the middle of the target Note that the layer order can matter! */ - double LayeredTarget::GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle) { + double LayeredTarget::GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth) { if(rxnLayer > m_layers.size()) { @@ -76,16 +75,15 @@ namespace Mask { } double eloss = 0.0; - RandomGenerator& gen = RandomGenerator::GetInstance(); if(angle < M_PI/2.0) { double newEnergy = startEnergy; - eloss += m_layers[rxnLayer].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator())); + eloss += m_layers[rxnLayer].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, (1.0 - rxnDepth)); newEnergy = startEnergy - eloss; if(rxnLayer == m_layers.size()) return eloss; - for(std::size_t i=rxnLayer; i m_layers.size()) { @@ -130,7 +128,7 @@ namespace Mask { eloss += m_layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle); newEnergy = startEnergy + eloss; } - eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator())); + eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, rxnDepth); newEnergy = startEnergy + eloss; } else @@ -141,7 +139,7 @@ namespace Mask { eloss += m_layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle); newEnergy = startEnergy + eloss; } - eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator())); + eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, rxnDepth); newEnergy = startEnergy + eloss; } diff --git a/src/Mask/LayeredTarget.h b/src/Mask/LayeredTarget.h index 2aed0a4..cf6dcf7 100644 --- a/src/Mask/LayeredTarget.h +++ b/src/Mask/LayeredTarget.h @@ -27,9 +27,9 @@ namespace Mask { LayeredTarget(); ~LayeredTarget(); void AddLayer(const std::vector& Z, const std::vector& A, const std::vector& stoich, double thickness); - double GetProjectileEnergyLoss(int zp, int ap, double startEnergy, std::size_t rxnLayer, double angle); - double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle); - double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle); + double GetProjectileEnergyLoss(int zp, int ap, 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); std::size_t FindLayerContaining(int Z, int A); std::size_t GetNumberOfLayers() { return m_layers.size(); } void SetName(std::string& n) { m_name = n; } @@ -39,7 +39,6 @@ namespace Mask { private: std::vector m_layers; std::string m_name; - std::uniform_real_distribution m_fractionalDepthDistribution; }; } diff --git a/src/Mask/OneStepSystem.cpp b/src/Mask/OneStepSystem.cpp index 746e838..1aa7ac1 100644 --- a/src/Mask/OneStepSystem.cpp +++ b/src/Mask/OneStepSystem.cpp @@ -80,7 +80,9 @@ namespace Mask { double rxnTheta = std::acos((m_thetaRanges[0])(gen)); double rxnPhi = (m_phiRanges[0])(gen); double residEx = (m_exDistributions[0])(gen); + double rxnDepth = (m_rxnDepthDist(gen)); + m_step1.SetReactionDepth(rxnDepth); m_step1.SetBeamKE(bke); m_step1.SetPolarRxnAngle(rxnTheta); m_step1.SetAzimRxnAngle(rxnPhi); diff --git a/src/Mask/Reaction.cpp b/src/Mask/Reaction.cpp index 04f2809..751dfa7 100644 --- a/src/Mask/Reaction.cpp +++ b/src/Mask/Reaction.cpp @@ -15,13 +15,13 @@ namespace Mask { Reaction::Reaction() : m_target(nullptr), m_projectile(nullptr), m_ejectile(nullptr), m_residual(nullptr), m_layeredTarget(nullptr), - m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_ejectThetaType(RxnThetaType::None), m_isInit(false), m_isResidEloss(false) + m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_rxnDepth(0.0), m_ejectThetaType(RxnThetaType::None), m_isInit(false), m_isResidEloss(false) { } Reaction::Reaction(Nucleus* target, Nucleus* projectile, Nucleus* ejectile, Nucleus* residual) : m_target(nullptr), m_projectile(nullptr), m_ejectile(nullptr), m_residual(nullptr), - m_layeredTarget(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_ejectThetaType(RxnThetaType::None), m_isResidEloss(false) + m_layeredTarget(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), m_rxnLayer(0), m_rxnDepth(0.0), m_ejectThetaType(RxnThetaType::None), m_isResidEloss(false) { BindNuclei(target, projectile, ejectile, residual); } @@ -67,7 +67,7 @@ namespace Mask { if(!m_isInit || m_isDecay) return; - m_bke = bke - m_layeredTarget->GetProjectileEnergyLoss(m_projectile->Z, m_projectile->A, bke, m_rxnLayer, 0); + m_bke = bke - m_layeredTarget->GetProjectileEnergyLoss(m_projectile->Z, m_projectile->A, bke, m_rxnLayer, 0, m_rxnDepth); } void Reaction::SetEjectileThetaType(RxnThetaType type) @@ -113,17 +113,27 @@ namespace Mask { m_residual->vec4 = m_target->vec4 + m_projectile->vec4 - m_ejectile->vec4; - ejectKE -= m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, ejectKE, m_rxnLayer, m_theta); - ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * m_ejectile->groundStateMass)); - ejectE = ejectKE + m_ejectile->groundStateMass; - m_ejectile->SetVec4Spherical(m_theta, m_phi, ejectP, ejectE); + ejectKE -= m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, ejectKE, m_rxnLayer, m_theta, m_rxnDepth); + if(ejectKE > 0.0) + { + double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0*m_ejectile->groundStateMass)); + double ejectE = ejectKE + m_ejectile->groundStateMass; + m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), ejectP, ejectE); + } + else + m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), 0.0, m_ejectile->groundStateMass); if(m_isResidEloss) { double residKE = m_residual->GetKE() - m_layeredTarget->GetEjectileEnergyLoss(m_residual->Z, m_residual->A, m_residual->GetKE(), - m_rxnLayer, m_residual->vec4.Theta()); - double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M())); - double residE = residKE + m_residual->vec4.M(); - m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE); + m_rxnLayer, m_residual->vec4.Theta(), m_rxnDepth); + if(residKE > 0.0) + { + double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M())); + double residE = residKE + m_residual->vec4.M(); + m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE); + } + else + m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), 0.0, m_residual->vec4.M()); } } @@ -158,19 +168,30 @@ namespace Mask { double ejectP = m_ejectile->vec4.P(); double ejectE = m_ejectile->vec4.E(); //energy loss for ejectile (after reaction!) - ejectKE -= m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, ejectKE, m_rxnLayer, m_ejectile->vec4.Theta()); - ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * m_ejectile->groundStateMass)); - ejectE = ejectKE + m_ejectile->groundStateMass; - m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), ejectP, ejectE); + ejectKE -= m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, ejectKE, m_rxnLayer, m_ejectile->vec4.Theta(), m_rxnDepth); + if(ejectKE > 0.0) + { + double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0*m_ejectile->groundStateMass)); + double ejectE = ejectKE + m_ejectile->groundStateMass; + m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), ejectP, ejectE); + } + else + m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), 0.0, m_ejectile->groundStateMass); + //if on, get eloss for residual (after reaction!) if(m_isResidEloss) { double residKE = m_residual->GetKE() - - m_layeredTarget->GetEjectileEnergyLoss(m_residual->Z, m_residual->A, m_residual->GetKE(), m_rxnLayer, m_residual->vec4.Theta()); - double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M())); - double residE = residKE + m_residual->vec4.M(); - m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE); + m_layeredTarget->GetEjectileEnergyLoss(m_residual->Z, m_residual->A, m_residual->GetKE(), m_rxnLayer, m_residual->vec4.Theta(), m_rxnDepth); + if(residKE > 0.0) + { + double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M())); + double residE = residKE + m_residual->vec4.M(); + m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE); + } + else + m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), 0.0, m_residual->vec4.M()); } } @@ -208,20 +229,31 @@ 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()); - double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0*m_ejectile->groundStateMass)); - double ejectE = ejectKE + m_ejectile->groundStateMass; - m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), ejectP, ejectE); - + m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, m_ejectile->GetKE(), m_rxnLayer, m_ejectile->vec4.Theta(), m_rxnDepth); + if(ejectKE > 0.0) + { + double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0*m_ejectile->groundStateMass)); + double ejectE = ejectKE + m_ejectile->groundStateMass; + m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), ejectP, ejectE); + } + else + m_ejectile->SetVec4Spherical(m_ejectile->vec4.Theta(), m_ejectile->vec4.Phi(), 0.0, m_ejectile->groundStateMass); + //if on, get eloss for *heavy* break up nucleus if(m_isResidEloss) { double residKE = m_residual->GetKE() - - m_layeredTarget->GetEjectileEnergyLoss(m_residual->Z, m_residual->A, m_residual->GetKE(), m_rxnLayer, m_residual->vec4.Theta()); - double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M())); - double residE = residKE + m_residual->vec4.M(); - m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE); + m_layeredTarget->GetEjectileEnergyLoss(m_residual->Z, m_residual->A, m_residual->GetKE(), m_rxnLayer, m_residual->vec4.Theta(), m_rxnDepth); + if(residKE > 0.0) + { + double residP = std::sqrt(residKE*(residKE + 2.0*m_residual->vec4.M())); + double residE = residKE + m_residual->vec4.M(); + m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), residP, residE); + } + else + m_residual->SetVec4Spherical(m_residual->vec4.Theta(), m_residual->vec4.Phi(), 0.0, m_residual->vec4.M()); } } diff --git a/src/Mask/Reaction.h b/src/Mask/Reaction.h index 7923139..12b9e07 100644 --- a/src/Mask/Reaction.h +++ b/src/Mask/Reaction.h @@ -32,6 +32,7 @@ namespace Mask { void SetPolarRxnAngle(double theta) { m_theta = theta; }; void SetAzimRxnAngle(double phi) { m_phi = phi; }; void SetExcitation(double ex) { m_ex = ex; }; + void SetReactionDepth(double depth) { m_rxnDepth = depth; } void BindTarget(Nucleus* nuc) { m_target = nuc; }; void BindProjectile(Nucleus* nuc) { m_projectile = nuc; }; @@ -59,7 +60,7 @@ namespace Mask { LayeredTarget* m_layeredTarget; //not owned by Reaction - double m_bke, m_theta, m_phi, m_ex; + double m_bke, m_theta, m_phi, m_ex, m_rxnDepth; int m_rxnLayer; RxnThetaType m_ejectThetaType; diff --git a/src/Mask/ReactionSystem.cpp b/src/Mask/ReactionSystem.cpp index 3a8f86d..327e0e5 100644 --- a/src/Mask/ReactionSystem.cpp +++ b/src/Mask/ReactionSystem.cpp @@ -10,7 +10,7 @@ namespace Mask { ReactionSystem::ReactionSystem() : - m_isTargetSet(false), m_isValid(true), m_rxnLayer(0), m_sysEquation("") + m_isTargetSet(false), m_isValid(true), m_rxnLayer(0), m_sysEquation(""), m_rxnDepthDist(0.0, 1.0) { } diff --git a/src/Mask/ReactionSystem.h b/src/Mask/ReactionSystem.h index c4761be..9a69ca5 100644 --- a/src/Mask/ReactionSystem.h +++ b/src/Mask/ReactionSystem.h @@ -63,6 +63,7 @@ namespace Mask { std::vector> m_beamDistributions, m_exDistributions; std::vector> m_thetaRanges, m_phiRanges; std::vector m_decayAngularDistributions; + std::uniform_real_distribution m_rxnDepthDist; bool m_isTargetSet; bool m_isValid; diff --git a/src/Mask/ThreeStepSystem.cpp b/src/Mask/ThreeStepSystem.cpp index 880de24..d66710e 100644 --- a/src/Mask/ThreeStepSystem.cpp +++ b/src/Mask/ThreeStepSystem.cpp @@ -131,16 +131,20 @@ namespace Mask { 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.SetBeamKE(bke); m_step1.SetPolarRxnAngle(rxnTheta); m_step1.SetAzimRxnAngle(rxnPhi); m_step1.SetExcitation(residEx); + m_step2.SetReactionDepth(rxnDepth); m_step2.SetPolarRxnAngle(decay1Theta); m_step2.SetAzimRxnAngle(decay1Phi); m_step2.SetExcitation(decay1Ex); + m_step3.SetReactionDepth(rxnDepth); m_step3.SetPolarRxnAngle(decay2Theta); m_step3.SetAzimRxnAngle(decay2Phi); m_step3.SetExcitation(decay2Ex); diff --git a/src/Mask/TwoStepSystem.cpp b/src/Mask/TwoStepSystem.cpp index 2145e2e..691367f 100644 --- a/src/Mask/TwoStepSystem.cpp +++ b/src/Mask/TwoStepSystem.cpp @@ -105,12 +105,15 @@ namespace Mask { double decay1Phi = m_phiRanges[1](gen); double residEx = (m_exDistributions[0])(gen); double decay2Ex = m_exDistributions[1](gen); - + double rxnDepth = (m_rxnDepthDist(gen)); + + m_step1.SetReactionDepth(rxnDepth); m_step1.SetBeamKE(bke); m_step1.SetPolarRxnAngle(rxnTheta); m_step1.SetAzimRxnAngle(rxnPhi); m_step1.SetExcitation(residEx); + m_step2.SetReactionDepth(rxnDepth); m_step2.SetPolarRxnAngle(decay1Theta); m_step2.SetAzimRxnAngle(decay1Phi); m_step2.SetExcitation(decay2Ex);