gwm17/Mask
1
0
Fork 0
mirror of https://github.com/gwm17/Mask.git synced 2024-11-22 18:28:51 -05:00
Code Issues Projects Releases Wiki Activity

Compare commits

base: gwm17:8b81dda70e80747338eccf793219909ee5cb94a8
Branches Tags
gwm17:main
gwm17:devel
..
compare: gwm17:db719c01828ddbfd945a44b94c919dd498a227d9
Branches Tags
gwm17:main
gwm17:devel

2 Commits
8b81dda70e ... db719c0182

Author SHA1 Message Date
Gordon McCann db719c0182 Update catima for fix in integrator 2022-12-01 09:40:23 -05:00
Gordon McCann 6b0834de5b Fix bug in forward ejectile energy loss where energy loss in reaction layer was run twice. Also better handle randomization of rxn location within rxn layer 2022-12-01 09:27:03 -05:00
12 changed files with 92 additions and 50 deletions
Show Stats Expand all files Collapse all files

2
src/Detectors/SabreArray.h
View File

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

2
src/Mask/DecaySystem.cpp
View File

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

22
src/Mask/LayeredTarget.cpp
View File

@ -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(); i++)
for(std::size_t i=(rxnLayer+1); i<m_layers.size(); i++)
{
eloss += m_layers[i].GetEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
@ -94,7 +92,7 @@ namespace Mask {
else
{ //Travelling backwards through target
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, rxnDepth);
newEnergy = startEnergy - eloss;
if(rxnLayer == 0)
return eloss;
@ -112,7 +110,7 @@ namespace Mask {
}
/*ReverseEnergyLoss version of GetEjectileEnergyLoss*/
double LayeredTarget::GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle)
double LayeredTarget::GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth)
{
if(rxnLayer > 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;
}

7
src/Mask/LayeredTarget.h
View File

@ -27,9 +27,9 @@ namespace Mask {
LayeredTarget();
~LayeredTarget();
void AddLayer(const std::vector<int>& Z, const std::vector<int>& A, const std::vector<int>& 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<Target> m_layers;
std::string m_name;
std::uniform_real_distribution<double> m_fractionalDepthDistribution;
};
}

2
src/Mask/OneStepSystem.cpp
View File

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

88
src/Mask/Reaction.cpp
View File

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

3
src/Mask/Reaction.h
View File

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

2
src/Mask/ReactionSystem.cpp
View File

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

1
src/Mask/ReactionSystem.h
View File

@ -63,6 +63,7 @@ namespace Mask {
std::vector<std::normal_distribution<double>> m_beamDistributions, m_exDistributions;
std::vector<std::uniform_real_distribution<double>> m_thetaRanges, m_phiRanges;
std::vector<AngularDistribution> m_decayAngularDistributions;
std::uniform_real_distribution<double> m_rxnDepthDist;
bool m_isTargetSet;
bool m_isValid;

6
src/Mask/ThreeStepSystem.cpp
View File

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

5
src/Mask/TwoStepSystem.cpp
View File

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

2
src/vendor/catima vendored

@ -1 +1 @@
Subproject commit 629690a6f90ad9d343e02f72d4fad860fce1ed6d
Subproject commit 671cda3c0f13f2eb38988ff0bfd5f5e9579b5fec