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

Compare commits

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

No commits in common. "db719c01828ddbfd945a44b94c919dd498a227d9" and "8b81dda70e80747338eccf793219909ee5cb94a8" have entirely different histories.
db719c0182 ... 8b81dda70e

12 changed files with 50 additions and 92 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.2; //in MeV
static constexpr double s_energyThreshold = 0.25; //in MeV
};

2
src/Mask/DecaySystem.cpp
View File

@ -73,9 +73,7 @@ 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_name(""), m_fractionalDepthDistribution(0.0, 1.0)
{
}
@ -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 rxnDepth)
double LayeredTarget::GetProjectileEnergyLoss(int zp, int ap, double startEnergy, std::size_t rxnLayer, double angle)
{
if(rxnLayer > m_layers.size())
{
@ -48,7 +48,8 @@ namespace Mask {
{
if(i == rxnLayer)
{
eloss += m_layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, rxnDepth);
frac = m_fractionalDepthDistribution(RandomGenerator::GetInstance().GetGenerator());
eloss += m_layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, frac);
newEnergy = startEnergy - eloss;
}
else
@ -66,7 +67,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 rxnDepth) {
double LayeredTarget::GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle) {
if(rxnLayer > m_layers.size())
{
@ -75,15 +76,16 @@ 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, (1.0 - rxnDepth));
eloss += m_layers[rxnLayer].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator()));
newEnergy = startEnergy - eloss;
if(rxnLayer == m_layers.size())
return eloss;
for(std::size_t i=(rxnLayer+1); i<m_layers.size(); i++)
for(std::size_t i=rxnLayer; i<m_layers.size(); i++)
{
eloss += m_layers[i].GetEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy - eloss;
@ -92,7 +94,7 @@ namespace Mask {
else
{ //Travelling backwards through target
double newEnergy = startEnergy;
eloss += m_layers[rxnLayer].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, rxnDepth);
eloss += m_layers[rxnLayer].GetEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator()));
newEnergy = startEnergy - eloss;
if(rxnLayer == 0)
return eloss;
@ -110,7 +112,7 @@ namespace Mask {
}
/*ReverseEnergyLoss version of GetEjectileEnergyLoss*/
double LayeredTarget::GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth)
double LayeredTarget::GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle)
{
if(rxnLayer > m_layers.size())
{
@ -128,7 +130,7 @@ namespace Mask {
eloss += m_layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy + eloss;
}
eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, rxnDepth);
eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator()));
newEnergy = startEnergy + eloss;
}
else
@ -139,7 +141,7 @@ namespace Mask {
eloss += m_layers[i].GetReverseEnergyLossTotal(ze, ae, newEnergy, angle);
newEnergy = startEnergy + eloss;
}
eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, rxnDepth);
eloss += m_layers[rxnLayer].GetReverseEnergyLossFractionalDepth(ze, ae, newEnergy, angle, m_fractionalDepthDistribution(gen.GetGenerator()));
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 rxnDepth);
double GetEjectileEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth);
double GetEjectileReverseEnergyLoss(int ze, int ae, double startEnergy, std::size_t rxnLayer, double angle, double rxnDepth);
double 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);
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,6 +39,7 @@ 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,9 +80,7 @@ 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_rxnDepth(0.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_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_rxnDepth(0.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_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_rxnDepth);
m_bke = bke - m_layeredTarget->GetProjectileEnergyLoss(m_projectile->Z, m_projectile->A, bke, m_rxnLayer, 0);
}
void Reaction::SetEjectileThetaType(RxnThetaType type)
@ -113,27 +113,17 @@ 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, 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);
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);
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(), 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());
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);
}
}
@ -168,30 +158,19 @@ 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(), 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);
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);
//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(), 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());
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);
}
}
@ -229,31 +208,20 @@ 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)
{
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);
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);
//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(), 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());
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);
}
}

3
src/Mask/Reaction.h
View File

@ -32,7 +32,6 @@ 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; };
@ -60,7 +59,7 @@ namespace Mask {
LayeredTarget* m_layeredTarget; //not owned by Reaction
double m_bke, m_theta, m_phi, m_ex, m_rxnDepth;
double m_bke, m_theta, m_phi, m_ex;
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_rxnDepthDist(0.0, 1.0)
m_isTargetSet(false), m_isValid(true), m_rxnLayer(0), m_sysEquation("")
{
}

1
src/Mask/ReactionSystem.h
View File

@ -63,7 +63,6 @@ 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,20 +131,16 @@ 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,15 +105,12 @@ 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 671cda3c0f13f2eb38988ff0bfd5f5e9579b5fec
Subproject commit 629690a6f90ad9d343e02f72d4fad860fce1ed6d