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

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

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

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, rxnDepth);
-				eloss += m_layers[i].GetEnergyLossFractionalDepth(zp, ap, newEnergy, angle, frac);
 				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;
 		}
 	

src/Mask/LayeredTarget.h

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

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

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,18 +113,28 @@ 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);
+		ejectKE -= m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, ejectKE, m_rxnLayer, m_theta, m_rxnDepth);
-		ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * m_ejectile->groundStateMass));
+		if(ejectKE > 0.0)
-		ejectE = ejectKE + m_ejectile->groundStateMass;
+		{
-		m_ejectile->SetVec4Spherical(m_theta, m_phi, ejectP, ejectE);
+			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());
+																						  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());
+		}
 	}
 	
 	//Methods from original ANASEN. Gives proper distribution for inverse kinematics.
@@ -158,20 +168,31 @@ 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());
+		ejectKE -= m_layeredTarget->GetEjectileEnergyLoss(m_ejectile->Z, m_ejectile->A, ejectKE, m_rxnLayer, m_ejectile->vec4.Theta(), m_rxnDepth);
-		ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * m_ejectile->groundStateMass));
+		if(ejectKE > 0.0)
-		ejectE = ejectKE + m_ejectile->groundStateMass;
+		{
+			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());
+					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());
+		}
 	}
 	
 	void Reaction::CalculateReaction()
@@ -208,21 +229,32 @@ 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_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());
+					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());
+		}
 	}
 
 }

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; 

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

src/Mask/ReactionSystem.h

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

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

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