Add better edge sampling of resonance excitations. Add new CoupledThreeStepSystem for angular correlations as well as optional correlation analysis in RootPlot

src/Detectors/SabreArray.cpp

@@ -20,6 +20,12 @@ SabreArray::SabreArray() :
 	m_degradedDetectors[2] = false;
 	m_degradedDetectors[3] = false;
 	m_degradedDetectors[4] = true;
+    //No degraded detectors
+	// m_degradedDetectors[0] = false;
+	// m_degradedDetectors[1] = false;
+	// m_degradedDetectors[2] = false;
+	// m_degradedDetectors[3] = false;
+	// m_degradedDetectors[4] = false;
 
 	//Choose who to look at right now. Usually switch on or off degraded/non-degraded.
 	m_activeDetectors[0] = false;
@@ -142,4 +148,4 @@ DetectorResult SabreArray::IsDetected(const Mask::Nucleus& nucleus)
 
 	observation.detectFlag = false;
 	return observation;
-}
+}

src/Mask/AngularDistribution.cpp

@@ -79,7 +79,7 @@ namespace Mask {
 		//Renormalize distribution such that total prob is 1.0.
 		//Test branching ratio to see if we "make" a decay particle,
 		//then use re-normalized distribution to pick an angle. 
-		if(m_constants[0] < 0.5)
+		if(m_constants[0] != 0.5)
 		{
 			double norm = 0.5/m_constants[0];
 			for(auto& value : m_constants)
@@ -115,4 +115,14 @@ namespace Mask {
 		return costheta;
 	}
 
+	double AngularDistribution::GetProbability(double cosTheta)
+	{
+		double prob = 0.0;
+		for (std::size_t i=0; i<m_constants.size(); i++)
+		{
+			prob += m_constants[i] * P_l(i*2, cosTheta);
+		}
+		return prob;
+	}
+
 }

src/Mask/AngularDistribution.h

@@ -17,6 +17,7 @@ namespace Mask {
 		double GetRandomCosTheta();
 		int GetL() { return m_L; }
 		double GetBranchingRatio() { return m_branchingRatio; }
+		double GetProbability(double cosTheta);
 	
 	private:
 		bool IsIsotropic() { return m_isIsotropic; }

src/Mask/CMakeLists.txt

@@ -55,6 +55,8 @@ target_sources(Mask PRIVATE
     FileWriter.cpp
     FileReader.h
     FileReader.cpp
+    CoupledThreeStepSystem.h
+    CoupledThreeStepSystem.cpp
 )
 
 set(THREADS_PREFER_PTHREAD_FLAG On)

src/Mask/CoupledThreeStepSystem.cpp

@@ -0,0 +1,210 @@
+#include "CoupledThreeStepSystem.h"
+
+#include "Math/Boost.h"
+#include "Math/Vector3D.h"
+#include "Math/RotationY.h"
+#include "Math/RotationZ.h"
+
+namespace Mask {
+
+    CoupledThreeStepSystem::CoupledThreeStepSystem(const std::vector<StepParameters>& params) :
+		ReactionSystem()
+	{
+		m_nuclei.resize(8);
+		Init(params);
+	}
+	
+	CoupledThreeStepSystem::~CoupledThreeStepSystem() {}
+	
+	void CoupledThreeStepSystem::Init(const std::vector<StepParameters>& params)
+	{
+		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;
+		}
+
+		const StepParameters& step1Params = params[0];
+		const StepParameters& step2Params = params[1];
+		const StepParameters& step3Params = params[2];
+
+		//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] << ")" <<std::endl;
+			return;  
+		}
+		int zb2 = step2Params.Z[0] - step2Params.Z[1];
+		int ab2 = step2Params.A[0] - step2Params.A[1];
+		if(zb2 != step3Params.Z[0] || ab2 != step3Params.A[0])
+		{
+			m_isValid = false;
+			std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step two and step three are not sequential! Step two heavy (Z,A): ("
+					  << zb2 << "," << ab2 << ") Step three target (Z,A): (" << step3Params.Z[0] << "," << step3Params.A[0] << ")" <<std::endl;
+			return;
+		}
+		int zb4 = step3Params.Z[0] - step3Params.Z[1];
+		int ab4 = step3Params.A[0] - step3Params.A[1];
+
+		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_nuclei[4] = CreateNucleus(step2Params.Z[1], step2Params.A[1]); //breakup1
+		m_nuclei[5] = CreateNucleus(zb2, ab2); //breakup2
+		m_nuclei[6] = CreateNucleus(step3Params.Z[1], step3Params.A[1]); //breakup3
+		m_nuclei[7] = CreateNucleus(zb4, ab4); //breakup4
+
+		m_step1.BindNuclei(&(m_nuclei[0]), &(m_nuclei[1]), &(m_nuclei[2]), &(m_nuclei[3]));
+		m_step2.BindNuclei(&(m_nuclei[3]), nullptr, &(m_nuclei[4]), &(m_nuclei[5]));
+		m_step3.BindNuclei(&(m_nuclei[5]), nullptr, &(m_nuclei[6]), &(m_nuclei[7]));
+	
+		SetSystemEquation();
+
+		//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 CoupledThreeStepSystem::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_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
+			throw ReactionLayerException();
+	}
+	
+	void CoupledThreeStepSystem::SetSystemEquation()
+	{
+		std::stringstream stream;
+		stream << m_nuclei[0].isotopicSymbol << "("
+			   << m_nuclei[1].isotopicSymbol << ", "
+			   << m_nuclei[2].isotopicSymbol << ")"
+			   << m_nuclei[3].isotopicSymbol << "->"
+			   << m_nuclei[4].isotopicSymbol << "+"
+			   << m_nuclei[5].isotopicSymbol << "->"
+			   << m_nuclei[6].isotopicSymbol << "+"
+			   << m_nuclei[7].isotopicSymbol;
+		m_sysEquation = stream.str();
+	}
+
+	CoupledThreeStepParameters CoupledThreeStepSystem::SampleParameters()
+	{
+		CoupledThreeStepParameters params;
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		params.beamEnergy = m_beamDistributions[0](gen);
+		params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
+		params.rxnPhi = m_phiRanges[0](gen);
+		params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
+        params.cosRelativeAngle = m_decayAngularDistributions[1].GetRandomCosTheta();
+		params.decay1Theta = std::acos(params.cosdecay1Theta);
+		params.decay1Phi = m_phiRanges[1](gen);
+		params.residEx = m_exDistributions[0](gen);
+		params.decay1Ex = m_exDistributions[1](gen);
+		params.decay2Ex = m_exDistributions[2](gen);
+		params.rxnDepth = (m_rxnDepthDist(gen));
+		return params;
+	}
+
+    //Called after running step2
+    void CoupledThreeStepSystem::SampleCoupling(CoupledThreeStepParameters& params)
+    {
+        ROOT::Math::Boost boost(m_nuclei[5].vec4.BoostToCM());
+		ROOT::Math::PxPyPzEVector a1Vec = (boost * m_nuclei[4].vec4);
+		ROOT::Math::XYZVector a1PVec(a1Vec.Px(), a1Vec.Py(), a1Vec.Pz());
+        ROOT::Math::PxPyPzEVector a2Vec;
+		ROOT::Math::XYZVector a2PVec;
+		ROOT::Math::RotationZ zRot;
+		ROOT::Math::RotationY yRot;
+		a1PVec *= 1.0/a1Vec.P();
+        zRot.SetAngle(a1PVec.Phi());
+		yRot.SetAngle(a1PVec.Theta());
+		double relAngle = std::acos(params.cosRelativeAngle);
+		ROOT::Math::XYZVector a1PVecAligned = yRot.Inverse() * (zRot.Inverse() * a1PVec);
+		a2PVec.SetXYZ(
+			std::sin(relAngle),
+			0.0,
+			std::cos(relAngle)
+		);
+		a2PVec = zRot * (yRot * a2PVec);
+		params.decay2Theta = a2PVec.Theta();
+		params.decay2Phi = a2PVec.Phi();
+    }
+	
+	void CoupledThreeStepSystem::RunSystem()
+	{
+		CoupledThreeStepParameters params;
+		do
+		{
+			params = SampleParameters();
+		}
+		while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) 
+				&& m_step2.CheckDecayThreshold(params.residEx, params.decay1Ex)
+				&& m_step3.CheckDecayThreshold(params.decay1Ex, params.decay2Ex)));
+		
+		m_step1.SetReactionDepth(params.rxnDepth);
+		m_step1.SetBeamKE(params.beamEnergy);
+		m_step1.SetPolarRxnAngle(params.rxnTheta);
+		m_step1.SetAzimRxnAngle(params.rxnPhi);
+		m_step1.SetExcitation(params.residEx);
+	
+		m_step2.SetReactionDepth(params.rxnDepth);
+		m_step2.SetPolarRxnAngle(params.decay1Theta);
+		m_step2.SetAzimRxnAngle(params.decay1Phi);
+		m_step2.SetExcitation(params.decay1Ex);
+	
+		m_step3.SetReactionDepth(params.rxnDepth);
+		
+		m_step3.SetExcitation(params.decay2Ex);
+	
+		m_step1.Calculate();
+	
+		if(params.cosdecay1Theta == -10)
+		{
+			m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
+			m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
+			m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
+			m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
+			return;
+		}
+		m_step2.Calculate();
+	
+		m_step3.SetResidualEnergyLoss(true);
+
+		SampleCoupling(params); //must be called here as frame needs to be set
+
+		m_step3.SetPolarRxnAngle(params.decay2Theta);
+		m_step3.SetAzimRxnAngle(params.decay2Phi);
+		m_step3.Calculate();
+	}
+}

src/Mask/CoupledThreeStepSystem.h

@@ -0,0 +1,46 @@
+#ifndef COUPLEDTHREESTEPSYSTEM_H
+#define COUPLEDTHREESTEPSYSTEM_H
+
+#include "ReactionSystem.h"
+#include "AngularDistribution.h"
+
+namespace Mask {
+
+	struct CoupledThreeStepParameters
+	{
+		double beamEnergy = 0.;
+		double rxnTheta = 0.;
+		double rxnPhi = 0.;
+		double cosdecay1Theta = 0.;
+		double cosRelativeAngle = 0.;
+		double decay1Theta = 0.;
+		double decay1Phi = 0.;
+		double decay2Theta = 0.;
+		double decay2Phi = 0.;
+		double residEx = 0.;
+		double decay1Ex = 0.;
+		double decay2Ex = 0.;
+		double rxnDepth = 0.;
+	};
+
+	class CoupledThreeStepSystem : public ReactionSystem
+	{
+	public:
+		CoupledThreeStepSystem(const std::vector<StepParameters>& params);
+		~CoupledThreeStepSystem();
+
+		virtual void SetLayeredTarget(const LayeredTarget& target) override;
+		virtual void RunSystem() override;
+		
+	protected:
+		void Init(const std::vector<StepParameters>& params);
+		void SetSystemEquation() override;
+		CoupledThreeStepParameters SampleParameters();
+        void SampleCoupling(CoupledThreeStepParameters& params);
+	
+		Reaction m_step1, m_step2, m_step3;
+	};
+
+}
+
+#endif

src/Mask/Reaction.cpp

@@ -256,6 +256,27 @@ namespace Mask {
 		}
 	}
 
+	bool Reaction::CheckReactionThreshold(double beamEnergy, double residEx)
+	{
+		double Q = m_target->groundStateMass + m_projectile->groundStateMass -
+				   (m_ejectile->groundStateMass + m_residual->groundStateMass + residEx);
+		double Ethresh = -Q*(m_ejectile->groundStateMass+m_residual->groundStateMass) / 
+						 (m_ejectile->groundStateMass + m_residual->groundStateMass - m_projectile->groundStateMass);
+		if (beamEnergy < Ethresh)
+			return false;
+		else
+			return true;
+	}
+
+	bool Reaction::CheckDecayThreshold(double targetEx, double residEx)
+	{
+		double Q = m_target->groundStateMass + targetEx - 
+				  (m_ejectile->groundStateMass + m_residual->groundStateMass + residEx);
+		if ( Q < 0.0)
+			return false;
+		else
+			return true;
+	}
 }
 
 

src/Mask/Reaction.h

@@ -42,6 +42,9 @@ namespace Mask {
 		void SetRxnLayer(std::size_t layer) { m_rxnLayer = layer; };
 		void SetResidualEnergyLoss(bool isEloss) { m_isResidEloss = isEloss; };
 
+		bool CheckReactionThreshold(double beamEnergy, double residualEx);
+		bool CheckDecayThreshold(double targetEx, double residualEx);
+
 		bool IsDecay() const { return m_isDecay; };
 
 		std::size_t GetRxnLayer() const { return m_rxnLayer; };

src/Mask/ThreeStepSystem.cpp

@@ -2,6 +2,11 @@
 #include "RandomGenerator.h"
 #include "KinematicsExceptions.h"
 
+#include "Math/Boost.h"
+#include "Math/Vector3D.h"
+#include "Math/RotationZ.h"
+#include "Math/RotationY.h"
+
 namespace Mask {
 	
 	ThreeStepSystem::ThreeStepSystem(const std::vector<StepParameters>& params) :
@@ -114,44 +119,58 @@ namespace Mask {
 			   << m_nuclei[7].isotopicSymbol;
 		m_sysEquation = stream.str();
 	}
+
+	ThreeStepParameters ThreeStepSystem::SampleParameters()
+	{
+		ThreeStepParameters params;
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		params.beamEnergy = m_beamDistributions[0](gen);
+		params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
+		params.rxnPhi = m_phiRanges[0](gen);
+		params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
+		params.decay1Theta = std::acos(params.cosdecay1Theta);
+		params.decay1Phi = m_phiRanges[1](gen);
+		params.cosdecay2Theta = m_decayAngularDistributions[1].GetRandomCosTheta();
+		params.decay2Theta = std::acos(params.cosdecay2Theta);
+		params.decay2Phi = m_phiRanges[1](gen);
+		params.residEx = m_exDistributions[0](gen);
+		params.decay1Ex = m_exDistributions[1](gen);
+		params.decay2Ex = m_exDistributions[2](gen);
+		params.rxnDepth = (m_rxnDepthDist(gen));
+		return params;
+	}
 	
 	void ThreeStepSystem::RunSystem()
 	{
-		//Sample parameters
+		ThreeStepParameters params;
-		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		do
-		double bke = m_beamDistributions[0](gen);
+		{
-		double rxnTheta = std::acos((m_thetaRanges[0])(gen));
+			params = SampleParameters();
-		double rxnPhi = m_phiRanges[0](gen);
+		}
-		double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
+		while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) 
-		double decay1Theta = std::acos(decay1costheta);
+				&& m_step2.CheckDecayThreshold(params.residEx, params.decay1Ex)
-		double decay1Phi = m_phiRanges[1](gen);
+				&& m_step3.CheckDecayThreshold(params.decay1Ex, params.decay2Ex)));
-		double decay2costheta = m_decayAngularDistributions[1].GetRandomCosTheta();
-		double decay2Theta = std::acos(decay2costheta);
-		double decay2Phi = m_phiRanges[2](gen);
-		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.SetReactionDepth(params.rxnDepth);
-		m_step1.SetBeamKE(bke);
+		m_step1.SetBeamKE(params.beamEnergy);
-		m_step1.SetPolarRxnAngle(rxnTheta);
+		m_step1.SetPolarRxnAngle(params.rxnTheta);
-		m_step1.SetAzimRxnAngle(rxnPhi);
+		m_step1.SetAzimRxnAngle(params.rxnPhi);
-		m_step1.SetExcitation(residEx);
+		m_step1.SetExcitation(params.residEx);
 	
-		m_step2.SetReactionDepth(rxnDepth);
+		m_step2.SetReactionDepth(params.rxnDepth);
-		m_step2.SetPolarRxnAngle(decay1Theta);
+		m_step2.SetPolarRxnAngle(params.decay1Theta);
-		m_step2.SetAzimRxnAngle(decay1Phi);
+		m_step2.SetAzimRxnAngle(params.decay1Phi);
-		m_step2.SetExcitation(decay1Ex);
+		m_step2.SetExcitation(params.decay1Ex);
 	
-		m_step3.SetReactionDepth(rxnDepth);
+		m_step3.SetReactionDepth(params.rxnDepth);
-		m_step3.SetPolarRxnAngle(decay2Theta);
+		m_step3.SetPolarRxnAngle(params.decay2Theta);
-		m_step3.SetAzimRxnAngle(decay2Phi);
+		m_step3.SetAzimRxnAngle(params.decay2Phi);
-		m_step3.SetExcitation(decay2Ex);
+		m_step3.SetResidualEnergyLoss(true);
+		m_step3.SetExcitation(params.decay2Ex);
 	
 		m_step1.Calculate();
 	
-		if(decay1costheta == -10)
+		if(params.cosdecay1Theta == -10)
 		{
 			m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
 			m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);
@@ -161,15 +180,13 @@ namespace Mask {
 		}
 		m_step2.Calculate();
 	
-		if(decay2costheta == -10)
+		if(params.cosdecay2Theta == -10)
 		{
 			m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass);
 			m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass);
 			return;
 		}
-		m_step3.SetResidualEnergyLoss(true);
 		m_step3.Calculate();
-	
 	}
 
 }

src/Mask/ThreeStepSystem.h

@@ -6,6 +6,23 @@
 
 namespace Mask {
 
+	struct ThreeStepParameters
+	{
+		double beamEnergy = 0.;
+		double rxnTheta = 0.;
+		double rxnPhi = 0.;
+		double cosdecay1Theta = 0.;
+		double cosdecay2Theta = 0.;
+		double decay1Theta = 0.;
+		double decay1Phi = 0.;
+		double decay2Theta = 0.;
+		double decay2Phi = 0.;
+		double residEx = 0.;
+		double decay1Ex = 0.;
+		double decay2Ex = 0.;
+		double rxnDepth = 0.;
+	};
+
 	class ThreeStepSystem : public ReactionSystem
 	{
 	public:
@@ -18,6 +35,7 @@ namespace Mask {
 	protected:
 		void Init(const std::vector<StepParameters>& params);
 		void SetSystemEquation() override;
+		ThreeStepParameters SampleParameters();
 	
 		Reaction m_step1, m_step2, m_step3;
 	

src/Mask/TwoStepSystem.cpp

@@ -92,35 +92,59 @@ namespace Mask {
 			   << m_nuclei[5].isotopicSymbol;
 		m_sysEquation = stream.str();
 	}
+
+	TwoStepParameters TwoStepSystem::SampleParameters()
+	{
+		TwoStepParameters params;
+		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		params.beamEnergy = (m_beamDistributions[0])(gen);
+		params.rxnTheta = std::acos((m_thetaRanges[0])(gen));
+		params.rxnPhi = (m_phiRanges[0])(gen);
+		params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta();
+		params.decay1Theta = std::acos(params.cosdecay1Theta);
+		params.decay1Phi = m_phiRanges[1](gen);
+		params.residEx = (m_exDistributions[0])(gen);
+		params.decay2Ex = m_exDistributions[1](gen);
+		params.rxnDepth = (m_rxnDepthDist(gen));
+		return params;
+	}
 	
 	void TwoStepSystem::RunSystem()
 	{
 		//Sample parameters
-		std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
+		// std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator();
-		double bke = (m_beamDistributions[0])(gen);
+		// double bke = (m_beamDistributions[0])(gen);
-		double rxnTheta = std::acos((m_thetaRanges[0])(gen));
+		// double rxnTheta = std::acos((m_thetaRanges[0])(gen));
-		double rxnPhi = (m_phiRanges[0])(gen);
+		// double rxnPhi = (m_phiRanges[0])(gen);
-		double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
+		// double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta();
-		double decay1Theta = std::acos(decay1costheta);
+		// double decay1Theta = std::acos(decay1costheta);
-		double decay1Phi = m_phiRanges[1](gen);
+		// double decay1Phi = m_phiRanges[1](gen);
-		double residEx = (m_exDistributions[0])(gen);
+		// double residEx = (m_exDistributions[0])(gen);
-		double decay2Ex = m_exDistributions[1](gen);
+		// double decay2Ex = m_exDistributions[1](gen);
-		double rxnDepth = (m_rxnDepthDist(gen));
+		// double rxnDepth = (m_rxnDepthDist(gen));
+
+		TwoStepParameters params;
+		do
+		{
+			params = SampleParameters();
+		}
+		while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) 
+				&& m_step2.CheckDecayThreshold(params.residEx, params.decay2Ex)));
 		
-		m_step1.SetReactionDepth(rxnDepth);
+		m_step1.SetReactionDepth(params.rxnDepth);
-		m_step1.SetBeamKE(bke);
+		m_step1.SetBeamKE(params.beamEnergy);
-		m_step1.SetPolarRxnAngle(rxnTheta);
+		m_step1.SetPolarRxnAngle(params.rxnTheta);
-		m_step1.SetAzimRxnAngle(rxnPhi);
+		m_step1.SetAzimRxnAngle(params.rxnPhi);
-		m_step1.SetExcitation(residEx);
+		m_step1.SetExcitation(params.residEx);
 	
-		m_step2.SetReactionDepth(rxnDepth);
+		m_step2.SetReactionDepth(params.rxnDepth);
-		m_step2.SetPolarRxnAngle(decay1Theta);
+		m_step2.SetPolarRxnAngle(params.decay1Theta);
-		m_step2.SetAzimRxnAngle(decay1Phi);
+		m_step2.SetAzimRxnAngle(params.decay1Phi);
-		m_step2.SetExcitation(decay2Ex);
+		m_step2.SetExcitation(params.decay2Ex);
 		
 		m_step1.Calculate();
 	
-		if(decay1costheta == -10)
+		if(params.cosdecay1Theta == -10)
 		{
 			m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass);
 			m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass);

src/Mask/TwoStepSystem.h

@@ -6,6 +6,19 @@
 
 namespace Mask {
 
+	struct TwoStepParameters
+	{
+		double beamEnergy = 0.;
+		double rxnTheta = 0.;
+		double rxnPhi = 0.;
+		double cosdecay1Theta = 0.;
+		double decay1Theta = 0.;
+		double decay1Phi = 0.;
+		double residEx = 0.;
+		double decay2Ex = 0.;
+		double rxnDepth = 0.;
+	};
+
 	class TwoStepSystem : public ReactionSystem
 	{
 	public:
@@ -18,6 +31,7 @@ namespace Mask {
 	private:
 		void Init(const std::vector<StepParameters>& params);
 		void SetSystemEquation() override;
+		TwoStepParameters SampleParameters();
 
 		Reaction m_step1, m_step2;
 	};

src/Plotters/RootPlotter.cpp

@@ -1,7 +1,8 @@
 #include "RootPlotter.h"
 #include <TFile.h>
 #include <TTree.h>
-#include <TVector3.h>
+#include <Math/Vector3D.h>
+#include "Math/Boost.h"
 
 #include <iostream>
 
@@ -10,7 +11,8 @@ static double FullPhi(double phi)
 	return phi < 0.0 ? 2.0*M_PI + phi : phi; 
 }
 
-RootPlotter::RootPlotter()
+RootPlotter::RootPlotter() :
+	m_target({5},{9},{1}, 74.0)
 {
 	TH1::AddDirectory(kFALSE);
 }
@@ -42,10 +44,16 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam
 			std::cout<<"\rPercent of data processed: "<<flushCount*flushFrac*100<<"%"<<std::flush;
 		}
 		tree->GetEntry(i);
-		for(Mask::Nucleus& nuc : *(dataHandle))
+		// for(Mask::Nucleus& nuc : *(dataHandle))
+		// {
+		// 	FillData(nuc);
+		// }
+		for(int i=0; i<dataHandle->size(); i++)
 		{
-			FillData(nuc);
+			FillData(dataHandle->at(i), i);
 		}
+		//Don't leave this in!
+		//Correlations(*dataHandle);
 	}
 	std::cout<<std::endl;
 	input->Close();
@@ -57,34 +65,35 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam
 	output->Close();
 }
 
-void RootPlotter::FillData(const Mask::Nucleus& nuc)
+void RootPlotter::FillData(const Mask::Nucleus& nuc, int i)
 {
-	std::string mod = "detected";
+	std::string mod = "_detected";
+	std::string num = "_" + std::to_string(i);
 
 	std::string sym = nuc.isotopicSymbol;
-	std::string ke_vs_th_name = sym + "_ke_vs_theta";
+	std::string ke_vs_th_name = sym + num + "_ke_vs_theta";
 	std::string ke_vs_th_title = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
-	std::string ke_vs_th_name_det = sym + mod + "_ke_vs_theta";
+	std::string ke_vs_th_name_det = sym + num + mod + "_ke_vs_theta";
 	std::string ke_vs_th_title_det = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)";
-	std::string ke_vs_ph_name = sym + "_ke_vs_phi";
+	std::string ke_vs_ph_name = sym + num + "_ke_vs_phi";
 	std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
-	std::string ke_vs_ph_name_det = sym + mod + "_ke_vs_phi";
+	std::string ke_vs_ph_name_det = sym + num + mod + "_ke_vs_phi";
 	std::string ke_vs_ph_title_det = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)";
-	std::string th_vs_ph_name = sym + "_theta_vs_phi";
+	std::string th_vs_ph_name = sym + num + "_theta_vs_phi";
 	std::string th_vs_ph_title = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
-	std::string th_vs_ph_name_det = sym + mod + "_theta_vs_phi";
+	std::string th_vs_ph_name_det = sym + num + mod + "_theta_vs_phi";
 	std::string th_vs_ph_title_det = th_vs_ph_name + ";#theta_{lab};#phi_{lab}";
-	std::string ex_name = sym + "_ex";
+	std::string ex_name = sym + num + "_ex";
 	std::string ex_title = ex_name + ";E_{ex} (MeV);counts";
-	std::string ex_name_det = sym + mod + "_ex";
+	std::string ex_name_det = sym + num + mod + "_ex";
 	std::string ex_title_det = ex_name + ";E_{ex} (MeV);counts";
-	std::string angdist_name = sym +"_angDist";
+	std::string angdist_name = sym + num + "_angDist";
 	std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts";
-	std::string angdist_name_det = sym + mod +"_angDist";
+	std::string angdist_name_det = sym + num +  mod +"_angDist";
 	std::string angdist_title_det = angdist_name+";cos#right(#theta_{CM}#left);counts";
-	std::string hist_ke_th_name = sym + "_hist_ke_vs_theta";
+	std::string hist_ke_th_name = sym + num + "_hist_ke_vs_theta";
 	std::string hist_ke_th_title = hist_ke_th_name + ";#theta_{lab};Kinetic Energy (MeV)";
-	std::string hist_ke_th_name_det = sym + mod + "_hist_ke_vs_theta";
+	std::string hist_ke_th_name_det = sym + num + mod + "_hist_ke_vs_theta";
 	std::string hist_ke_th_title_det = hist_ke_th_name + ";#theta_{lab};Kinetic Energy (MeV)";
 	
 	MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2);
@@ -105,6 +114,93 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc)
 	
 }
 
+//Correlation analysis for coupled step reactions, will need to be customized
+//for a given reaction. By default it is off.
+void RootPlotter::Correlations(const std::vector<Mask::Nucleus>& event)
+{
+	static constexpr double lowerBoundEx = -2.7;
+    static constexpr double upperBoundEx = 2.7;
+	static constexpr double targetDepth = 0.5;
+	const Mask::Nucleus& secondary = event[6];
+	const Mask::Nucleus& primary = event[4];
+	const Mask::Nucleus& parent = event[3];
+    const Mask::Nucleus& li = event[5];
+	if(secondary.vec4.P() == 0.0 || primary.vec4.P() == 0.0)
+    {
+        return;
+    }
+	ROOT::Math::Boost boostParent(parent.vec4.BoostToCM());
+	ROOT::Math::Boost boostIntermediate(li.vec4.BoostToCM());
+	ROOT::Math::PxPyPzEVector a2Vec = (secondary.vec4);
+    ROOT::Math::PxPyPzEVector a1Vec = (primary.vec4);
+
+	double a2KE = secondary.GetKE() + m_target.GetReverseEnergyLossFractionalDepth(secondary.Z, secondary.A, secondary.GetKE(), a2Vec.Theta(), targetDepth);
+	double a2P = std::sqrt(a2KE * (a2KE + 2.0 * secondary.groundStateMass));
+	a2Vec.SetPxPyPzE(
+		a2P * std::sin(a2Vec.Theta()) * std::cos(a2Vec.Phi()),
+		a2P * std::sin(a2Vec.Theta()) * std::sin(a2Vec.Phi()),
+		a2P * std::cos(a2Vec.Theta()),
+		a2KE + secondary.groundStateMass
+	);
+
+	double a1KE = primary.GetKE() + m_target.GetReverseEnergyLossFractionalDepth(primary.Z, primary.A, primary.GetKE(), a1Vec.Theta(), targetDepth);
+	double a1P = std::sqrt(a1KE * (a1KE + 2.0 * primary.groundStateMass));
+	a1Vec.SetPxPyPzE(
+		a1P * std::sin(a1Vec.Theta()) * std::cos(a1Vec.Phi()),
+		a1P * std::sin(a1Vec.Theta()) * std::sin(a1Vec.Phi()),
+		a1P * std::cos(a1Vec.Theta()),
+		a1KE + primary.groundStateMass
+	);
+
+	ROOT::Math::XYZVector a1PVec;
+	ROOT::Math::XYZVector a2PVec;
+
+	ROOT::Math::PxPyPzEVector boostedIntA1Vec = boostIntermediate * a1Vec;
+	ROOT::Math::PxPyPzEVector boostedIntA2Vec = boostIntermediate * a2Vec;
+	ROOT::Math::PxPyPzEVector boostedParA2Vec = boostParent * a2Vec;
+    a1PVec.SetXYZ(boostedIntA1Vec.Px(), boostedIntA1Vec.Py(), boostedIntA1Vec.Pz());
+    a2PVec.SetXYZ(boostedIntA2Vec.Px(), boostedIntA2Vec.Py(), boostedIntA2Vec.Pz());
+    double relCosThetaCM = a1PVec.Dot(a2PVec) / (a1PVec.R() * a2PVec.R());
+
+    ROOT::Math::PxPyPzEVector badLi = parent.vec4 - a2Vec;
+	double secondaryExcitation = badLi.M() - li.groundStateMass;
+
+	MyFill("seoncdaryEx","secondaryEx",1000,-5.0,5.0,secondaryExcitation);
+	MyFill("primaryEx","primaryEx",1000,-5.0,5.0,li.GetExcitationEnergy());
+	MyFill("relativeDist_full","relativeDist_full",20,-1.0,1.0, relCosThetaCM);
+	MyFill("primaryDist_full","primaryDist_full",20,-1.0,1.0,std::cos(primary.thetaCM));
+	MyFill("seondaryDist_full","seconaryDist_full",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
+	MyFill("secondaryDist_fullOwn","seondaryDist_fullOwn",20,-1.0,1.0,std::cos(secondary.thetaCM));
+
+    bool isPrimary = false;
+    bool isSecondary = false;
+    if(primary.isDetected)
+    {
+		MyFill("primaryExMatched","primaryExMatched",1000,-5.0,5.0,li.GetExcitationEnergy());
+		MyFill("primaryDist","primaryDist",20,-1.0,1.0,std::cos(primary.thetaCM));
+        isPrimary = true;
+    }
+    if( secondary.isDetected && secondaryExcitation > lowerBoundEx && secondaryExcitation < upperBoundEx )
+    {
+		MyFill("secondaryDist","secondaryDist",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
+		MyFill("secondaryExMatched","secondaryExMatched",1000,-5.0,5.0,secondaryExcitation);
+        isSecondary = true;
+    }
+
+    if(isSecondary && isPrimary)
+    {
+		return;
+    }
+    else if (isSecondary)
+    {
+        MyFill("completeDist","completeDist",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta()));
+    }
+    else if (isPrimary)
+    {
+        MyFill("completeDist","completeDist",20,-1.0,1.0,std::cos(primary.thetaCM));
+    }
+}
+
 void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val)
 {
 	auto iter = m_map.find(name);

src/Plotters/RootPlotter.h

@@ -6,6 +6,7 @@
 #include <unordered_map>
 
 #include "Nucleus.h"
+#include "Target.h"
 
 #include <TH1.h>
 #include <TH2.h>
@@ -20,7 +21,8 @@ public:
 	void Run(const std::string& inputname, const std::string& outputname);
 	
 private:
-	void FillData(const Mask::Nucleus& nuc);
+	void FillData(const Mask::Nucleus& nuc, int i);
+	void Correlations(const std::vector<Mask::Nucleus>& event);
 	void MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val);
 	void MyFill(const std::string& name, const std::string& title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, 
 				double valx, double valy);
@@ -29,6 +31,7 @@ private:
 	std::unordered_map<std::string, std::shared_ptr<TObject>> m_map;
 
 	static constexpr double s_rad2deg = 180.0/M_PI;
+	Mask::Target m_target;
 };
 
 #endif