Added in punch through reconstruction. Everything seems functional

src/CMakeLists.txt

@@ -26,6 +26,12 @@ target_sources(SabreRecon PUBLIC
 	EnergyLoss/EnergyLoss.cpp
 	EnergyLoss/Target.h
 	EnergyLoss/Target.cpp
+	EnergyLoss/CubicSpline.h
+	EnergyLoss/CubicSpline.cpp
+	EnergyLoss/PunchTable.h
+	EnergyLoss/PunchTable.cpp
+	EnergyLoss/ElossTable.h
+	EnergyLoss/ElossTable.cpp
 	main.cpp
 	)
 target_link_libraries(SabreRecon

src/EnergyLoss/CubicSpline.cpp

@@ -0,0 +1,235 @@
+/*
+CubicSpline.cpp
+Class for generating cubic splines of data in tables or held in memory. Cubic splines are a form of interpolation,
+somewhat more advanced than linear interpolation, but not so complicated that it significantly slows down a calculation.
+For more information see Wikipedia or Num. Rec. in C.
+
+Gordon M. May 2021
+*/
+#include "CubicSpline.h"
+#include <fstream>
+#include <iostream>
+#include <cmath>
+
+namespace PunchTable {
+
+	CubicSpline::CubicSpline() :
+		m_validFlag(false)
+	{
+	}
+
+	CubicSpline::CubicSpline(const std::string& filename) :
+		m_validFlag(false)
+	{
+		ReadFile(filename);
+	}
+
+	CubicSpline::CubicSpline(const std::vector<double>& x, const std::vector<double>& y) :
+		m_validFlag(false)
+	{
+		ReadData(x, y);
+	}
+
+	CubicSpline::~CubicSpline() {}
+
+	/*
+	Expected file format is a naked (no header) single space separated table:
+	x y\n
+	*/
+	void CubicSpline::ReadFile(const std::string& filename)
+	{
+		std::ifstream input(filename);
+		if(!input.is_open()) {
+			std::cerr<<"Unable to open input data at CubicSpline::ReadFile from filename: "<<filename<<std::endl;
+			m_validFlag = false;
+			return;
+		}
+
+		std::string junk;
+		//std::getline(input, junk);
+
+		double x, y;
+		while(input>>x)
+		{
+			input>>y;
+			m_dataX.push_back(x);
+			m_dataY.push_back(y);
+		}
+
+		if(m_dataX.size() != m_dataY.size())
+		{
+			std::cerr<<"Error in CubicSpline::ReadFile! Number of x points not equal to number of y points!"<<std::endl;
+			m_validFlag = false;
+			return;
+		}
+
+		m_validFlag = true;
+
+		input.close();
+
+		MakeSplines();
+	}
+
+	/*
+	After data is read in splines can be solved. Each data point is referred to as a knot. Endpoint conditions are
+	derivatives of neighbors must be equal and second derivatives must be zero (natural cubic splines). Solved using gaussian elimination.
+	*/
+	void CubicSpline::MakeSplines()
+	{
+		if(!m_validFlag)
+		{
+			std::cerr<<"Error at CubicSpline::MakeSplines! Unable to generate splines without first initializing data."<<std::endl;
+			return;
+		}
+
+		int knots = m_dataX.size();
+
+		//Matrix and vector data init
+		double** a = new double*[knots];
+		for(int i=0; i<knots; i++) 
+			a[i] = new double[knots];
+		double* b = new double[knots];
+		double* k = new double[knots];
+
+		Spline s;
+		m_splines.clear();
+		double x0, x1, x2;
+		double y0, y1, y2;
+
+		//Setup matrix eqn.
+		for(int i=0; i<knots; i++)
+		{
+			if(i == 0)
+			{
+				x1 = m_dataX[i];
+				x2 = m_dataX[i+1];
+				y1 = m_dataY[i];
+				y2 = m_dataY[i+1];
+				a[i][i] = 2.0/(x2-x1);
+				a[i][i+1] = 1.0/(x2-x1);
+				b[i] = 3.0*(y2-y1)/(std::pow((x2-x1), 2.0));
+	      		s.x1 = x1; s.x2 = x2;
+	      		s.y1 = y1; s.y2 = y2;
+	      		m_splines.push_back(s);
+			}
+			else if(i == (knots-1))
+			{
+				x0 = m_dataX[i-1];
+				x1 = m_dataX[i];
+	      		y0 = m_dataY[i-1];
+	      		y1 = m_dataY[i];
+	      		a[i][i-1] = 1.0/(x1-x0);
+	      		a[i][i] = 2.0/(x1-x0);
+	      		b[i] = 3.0*(y1-y0)/(std::pow((x1-x0), 2.0));
+			}
+			else
+			{
+		    	x0 = m_dataX[i-1]; 
+		    	x1 = m_dataX[i]; 
+		    	x2 = m_dataX[i+1];
+	      		y0 = m_dataY[i-1];
+	      		y1 = m_dataY[i];
+	      		y2 = m_dataY[i+1];
+	      		a[i][i-1] = 1.0/(x1-x0);
+	      		a[i][i] = 2.0/(x1-x0)+2.0/(x2-x1);
+	      		a[i][i+1] = 1.0/(x2-x1);
+	      		b[i] = 3.0*(y1-y0)/(std::pow((x1-x0), 2.0))+3.0*(y2-y1)/(std::pow((x2-x1), 2.0));
+	      		s.x1 = x1; s.x2 = x2;
+	      		s.y1 = y1; s.y2 = y2;
+	      		m_splines.push_back(s);
+			}
+		}
+
+		//solve for curvature vector k using gaussian elimination
+		a[0][1] /= a[0][0];
+		b[0] /= a[0][0];
+		for(int i=1; i<(knots-1); i++)
+		{
+			a[i][i+1] /= a[i][i]-a[i][i-1]*a[i-1][i];
+	   	 	b[i] = (b[i] - a[i][i-1]*b[i-1])/(a[i][i]-a[i][i-1]*a[i-1][i]);
+		}
+		int g1 = knots-1;
+	  	int g2 = knots-2;
+	  	b[g1] = (b[g1]-a[g1][g2]*b[g2])/(a[g1][g1]-a[g1][g2]*a[g2][g1]);
+
+	  	k[g1] = b[g1];
+	  	for(int i=(knots-2); i>=0; i--)
+	    	k[i] = b[i] - a[i][i+1]*k[i+1];
+
+	  	//Fill the spline data
+	  	for(size_t i=0; i<m_splines.size(); i++)
+	  	{
+	  		m_splines[i].k1 = k[i];
+	  		m_splines[i].k2 = k[i+1];
+	  	}
+
+	  	//deallocate
+	  	delete[] b;
+	  	delete[] k;
+	  	for(int i=0; i<knots; i++)
+	  		delete[] a[i];
+	  	delete[] a;
+	}
+
+	double CubicSpline::Evaluate(double x)
+	{
+		if(!m_validFlag)
+		{
+			std::cerr<<"Error at CubicSpline::Evaluate! Unable to evaluate without first generating splines."<<std::endl;
+			return 0.0;
+		}
+
+		Spline s;
+		for(size_t i=0; i<m_splines.size(); i++)
+		{
+			auto& spline = m_splines[i];
+			if(x >= spline.x1 && x <= spline.x2)
+			{
+				s = spline;
+				break;
+			}
+			else if (i == (m_splines.size() -1))
+			{
+				//std::cerr<<"Error at CubicSpline::Evaluate! Input x value: "<<x<<" is not within the spline range min: "<<m_splines[0].x1<<" max: "<<m_splines[m_splines.size()-1].x2<<std::endl;
+				return 0.0;
+			}
+		}
+
+		double t = (x-s.x1)/(s.x2-s.x1);
+	    double a = s.k1*(s.x2-s.x1)-(s.y2-s.y1);
+	    double b = -s.k2*(s.x2-s.x1)+(s.y2-s.y1);
+	    return (1.0-t)*s.y1+t*s.y2+t*(1.0-t)*((1.0-t)*a+t*b);
+	}
+
+	//Purely for plotting in ROOT, do not use for caluculations. 
+	double CubicSpline::EvaluateROOT(double* x, double* p)
+	{
+		if(!m_validFlag)
+		{
+			std::cerr<<"Error at CubicSpline::EvaluateROOT! Unable to evaluate without first generating splines."<<std::endl;
+			return 0.0;
+		}
+
+		double xval = x[0];
+		Spline s;
+		for(size_t i=0; i<m_splines.size(); i++)
+		{
+			auto& spline = m_splines[i];
+			if(xval >= spline.x1 && xval <= spline.x2)
+			{
+				s = spline;
+				break;
+			}
+			else if (i == (m_splines.size() -1))
+			{
+				return 0.0;
+			}
+		}
+
+		double t = (xval-s.x1)/(s.x2-s.x1);
+	    double a = s.k1*(s.x2-s.x1)-(s.y2-s.y1);
+	    double b = -s.k2*(s.x2-s.x1)+(s.y2-s.y1);
+	    return (1.0-t)*s.y1+t*s.y2+t*(1.0-t)*((1.0-t)*a+t*b);
+	}
+
+}

src/EnergyLoss/CubicSpline.h

@@ -0,0 +1,55 @@
+/*
+CubicSpline.h 
+Class for generating cubic splines of data in tables or held in memory. Cubic splines are a form of interpolation,
+somewhat more advanced than linear interpolation, but not so complicated that it significantly slows down a calculation.
+For more information see Wikipedia or Num. Rec. in C.
+
+Gordon M. May 2021
+*/
+#ifndef CUBICSPLINE_H
+#define CUBICSPLINE_H
+
+#include <vector>
+#include <string>
+
+namespace PunchTable {
+
+	//struct holding final spline info
+	struct Spline
+	{
+		double y1=0, y2=0;
+		double x1=0, x2=0;
+		double k1=0, k2=0;
+	};
+
+	class CubicSpline
+	{
+	public:
+		CubicSpline();
+		CubicSpline(const std::string& filename);
+		CubicSpline(const std::vector<double>& x, const std::vector<double>& y);
+		~CubicSpline();
+		void ReadFile(const std::string& filename);
+		inline void ReadData(const std::vector<double>& x, const std::vector<double>& y)
+		{ 
+			m_dataX = x; 
+			m_dataY = y; 
+			m_validFlag=true;
+			MakeSplines();
+		}
+		bool IsValid() { return m_validFlag; }
+		double Evaluate(double x);
+		double EvaluateROOT(double* x, double* p); //for plotting as ROOT function
+
+	private:
+		void MakeSplines();
+
+		std::vector<double> m_dataX;
+		std::vector<double> m_dataY;
+		std::vector<Spline> m_splines;
+
+		bool m_validFlag;
+	};
+
+}
+#endif

src/EnergyLoss/ElossTable.cpp

@@ -0,0 +1,120 @@
+#include "ElossTable.h"
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+
+namespace PunchTable {
+
+    ElossTable::ElossTable() :
+        m_isValid(false)
+    {
+    }
+
+    ElossTable::ElossTable(const std::string& filename) :
+        m_isValid(false)
+    {
+        ReadFile(filename);
+    }
+
+    ElossTable::~ElossTable() {}
+
+    void ElossTable::ReadFile(const std::string& filename)
+    {
+        std::ifstream input(filename);
+		if(!input.is_open())
+		{
+			std::cerr<<"Unable to open table file named "<<filename<<"! Exiting."<<std::endl;
+			m_isValid = false;
+			return;
+		}
+
+		std::string junk;
+		double thickness;
+		double theta;
+		double value;
+		std::vector<double> energyLoss, energyFinal;
+
+		input>>junk>>junk>>m_projectileString;
+		input>>junk>>junk;
+
+		while(input>>junk)
+		{
+			if(junk == "---------------------------------")
+				break;
+			input>>junk;
+			m_materialString += junk;
+			input>>junk;
+		}
+		input>>junk>>m_thetaMin>>junk>>m_thetaMax>>junk>>m_thetaStep;
+		std::getline(input, junk);
+		std::getline(input, junk);
+		std::getline(input, junk);
+		std::getline(input, junk);
+		while(input>>junk)
+		{
+			if(junk == "begin_theta")
+			{
+				energyLoss.clear();
+				energyFinal.clear();
+				input>>theta;
+				while(input>>junk)
+				{
+					if(junk == "end_theta")
+						break;
+					energyFinal.push_back(std::stod(junk));
+					input>>value;
+					energyLoss.push_back(value);
+				}
+
+				if(!energyFinal.empty())
+					m_splines.emplace_back(energyFinal, energyLoss);
+				else
+					m_splines.emplace_back();
+			}
+			else
+			{
+				std::cerr<<"Unexpected expression found when reading punch table: "<<junk<<std::endl;
+				m_isValid = false;
+				return;
+			}
+		}
+
+		m_isValid = true;
+    }
+
+    double ElossTable::GetEnergyLoss(double thetaIncident, double finalEnergy)
+    {
+        thetaIncident /= s_deg2rad;
+		if(!m_isValid)
+		{
+			std::cerr<<"ElossTable not initialized at GetEnergyLoss()"<<std::endl;
+			return 0.0;
+		}
+		else if(thetaIncident < m_thetaMin || thetaIncident > m_thetaMax)
+		{
+			//std::cerr<<"Theta incident "<<thetaIncident<<" outside of range of calculated values for ElossTable::GetEnergyLoss()"<<std::endl;
+			return 0.0;
+		}
+
+		int theta_bin = (thetaIncident - m_thetaMin)/m_thetaStep;
+
+		//std::cout<<"theta bin: "<<theta_bin<<" theta_inc: "<<thetaIncident<<std::endl;
+
+		if(m_splines[theta_bin].IsValid())
+		{
+			double eloss =  m_splines[theta_bin].Evaluate(finalEnergy);
+			if(eloss == 0.0) //Not in the spline
+			{
+				//std::cout<<"here"<<std::endl;
+				return 0.0;
+			}
+			else
+				return eloss;
+		}
+		else
+        {
+            std::cerr<<"Spline error at ElossTable::GetEnergyLoss()!"<<std::endl;
+			return 0.0;
+        }
+    }
+}

src/EnergyLoss/ElossTable.h

@@ -0,0 +1,43 @@
+#ifndef ELOSS_TABLE_H
+#define ELOSS_TABLE_H
+
+#include "CubicSpline.h"
+#include <vector>
+#include <string>
+#include <cmath>
+
+namespace PunchTable {
+
+    /*
+        Class for interpolating a reverse energy loss file. Note that
+        this is for REVERSE energy loss, which is typically what is interesting
+        for experimental data.
+    */
+    class ElossTable
+    {
+    public:
+        ElossTable();
+        ElossTable(const std::string& filename);
+        ~ElossTable();
+
+        void ReadFile(const std::string& filename);
+        std::string GetProjectile() { return m_projectileString; }
+        std::string GetMaterial() { return m_materialString; }
+
+        double GetEnergyLoss(double thetaIncident, double finalEnergy);
+        
+        inline const bool IsValid() const { return m_isValid; }
+
+    private:
+        std::string m_projectileString;
+        std::string m_materialString;
+        std::vector<CubicSpline> m_splines;
+        double m_thetaStep, m_thetaMin, m_thetaMax;
+
+        bool m_isValid;
+
+        static constexpr double s_deg2rad = M_PI/180.0;
+    };
+}
+
+#endif

src/EnergyLoss/PunchTable.cpp

@@ -0,0 +1,119 @@
+#include "PunchTable.h"
+#include <fstream>
+#include <iostream>
+
+namespace PunchTable {
+
+	PunchTable::PunchTable() :
+		m_validFlag(false)
+	{
+	}
+
+	PunchTable::PunchTable(const std::string& filename) :
+		m_validFlag(false)
+	{
+		ReadFile(filename);
+	}
+
+	PunchTable::~PunchTable() {}
+
+	void PunchTable::ReadFile(const std::string& filename)
+	{
+		std::ifstream input(filename);
+		if(!input.is_open())
+		{
+			std::cerr<<"Unable to open table file named "<<filename<<"! Exiting."<<std::endl;
+			m_validFlag = false;
+			return;
+		}
+
+		std::string junk;
+		double thickness;
+		double theta;
+		double value;
+		std::vector<double> energyIn, energyDep;
+
+		input>>junk>>junk>>m_projectileString;
+		input>>junk>>junk;
+
+		while(input>>junk)
+		{
+			if(junk == "---------------------------------")
+				break;
+			input>>junk;
+			m_materialString += junk;
+			input>>junk;
+		}
+		
+		input>>junk>>m_thetaMin>>junk>>m_thetaMax>>junk>>m_thetaStep;
+		std::getline(input, junk);
+		std::getline(input, junk);
+		std::getline(input, junk);
+		std::getline(input, junk);
+		while(input>>junk)
+		{
+			if(junk == "begin_theta")
+			{
+				energyIn.clear();
+				energyDep.clear();
+				input>>theta;
+				while(input>>junk)
+				{
+					if(junk == "end_theta")
+						break;
+					energyDep.push_back(std::stod(junk));
+					input>>value;
+					energyIn.push_back(value);
+				}
+
+				if(!energyDep.empty())
+					m_splines.emplace_back(energyDep, energyIn);
+				else
+					m_splines.emplace_back();
+			}
+			else
+			{
+				std::cerr<<"Unexpected expression found when reading punch table: "<<junk<<std::endl;
+				m_validFlag = false;
+				return;
+			}
+		}
+
+		m_validFlag = true;
+
+	}
+
+	double PunchTable::GetInitialKineticEnergy(double theta_incident, double e_deposited)
+	{
+		theta_incident /= s_deg2rad;
+		if(!m_validFlag)
+		{
+			std::cerr<<"PunchTable not initialized at GetInitialKineticEnergy()"<<std::endl;
+			return 0.0;
+		}
+		else if(theta_incident < m_thetaMin || theta_incident > m_thetaMax)
+		{
+			//std::cerr<<"Theta incident outside of range of calculated values for PunchTable::GetInitialKineticEnergy"<<std::endl;
+			return 0.0;
+		}
+
+		int theta_bin = (theta_incident - m_thetaMin)/m_thetaStep;
+
+		//std::cout<<"theta bin: "<<theta_bin<<" theta_inc: "<<theta_incident<<std::endl;
+
+		if(m_splines[theta_bin].IsValid())
+		{
+			double initialE =  m_splines[theta_bin].Evaluate(e_deposited);
+			if(initialE == 0.0) //Not in the spline, stopped completely
+			{
+				//std::cout<<"here"<<std::endl;
+				return e_deposited;
+			}
+			else
+				return initialE;
+		}
+		else
+			return e_deposited;
+	}
+
+}

src/EnergyLoss/PunchTable.h

@@ -0,0 +1,36 @@
+#ifndef PUNCHTABLE_H
+#define PUNCHTABLE_H
+
+#include "CubicSpline.h"
+#include <vector>
+#include <string>
+#include <cmath>
+
+namespace PunchTable {
+
+	class PunchTable {
+	public:
+		PunchTable();
+		PunchTable(const std::string& filename);
+		~PunchTable();
+
+		void ReadFile(const std::string& filename);
+		std::string GetProjectile() { return m_projectileString; }
+        std::string GetMaterial() { return m_materialString; }
+
+		double GetInitialKineticEnergy(double theta_incident, double e_deposited); //radians, MeV
+		inline bool IsValid() const { return m_validFlag; }
+
+	private:
+		std::string m_projectileString;
+        std::string m_materialString;
+		std::vector<CubicSpline> m_splines;
+		double m_thetaStep, m_thetaMin, m_thetaMax;
+
+		bool m_validFlag;
+
+		static constexpr double s_deg2rad = M_PI/180.0;
+	};
+
+}
+#endif	

src/Histogrammer.cpp

@@ -20,6 +20,12 @@ namespace SabreRecon {
 	{
 		TH1::AddDirectory(kFALSE);
 		ParseConfig(input);
+		//As of Ubuntu 22.04: linker aggresively strips (particularly for release) to the point where non-header symbols MUST be called AND used
+		//within the executable. Use EnforceDictionaryLinked as shown.
+		if(EnforceDictionaryLinked())
+		{
+			std::cout<<"Enforcing dictionary linking"<<std::endl;
+		}
 	}
 
 	Histogrammer::~Histogrammer() {}
@@ -46,6 +52,9 @@ namespace SabreRecon {
 		std::vector<ReconCut> cuts;
 		ReconCut this_cut;
 
+		std::vector<std::string> ptables;
+		std::vector<std::string> etables;
+
 		input>>junk;
 		if(junk == "begin_data")
 		{
@@ -115,6 +124,28 @@ namespace SabreRecon {
 						std::cout<<std::endl;
 					}
 				}
+				else if(junk == "begin_punchtables")
+				{
+					std::cout<<"Looking for PunchTables..."<<std::endl;
+					while(input>>junk)
+					{
+						if(junk == "end_punchtables")
+							break;
+						ptables.push_back(junk);
+						std::cout<<"Adding PunchTable: "<<junk<<std::endl;
+					}
+				}
+				else if(junk == "begin_elosstables")
+				{
+					std::cout<<"Looking for ElossTables..."<<std::endl;
+					while(input>>junk)
+					{
+						if(junk == "end_elosstables")
+							break;
+						etables.push_back(junk);
+						std::cout<<"Adding ElossTable: "<<junk<<std::endl;
+					}
+				}
 				else if(junk == "end_focalplane")
 					continue;
 				else if(junk == "end_target")
@@ -148,6 +179,10 @@ namespace SabreRecon {
 		std::cout<<"Initializing resources..."<<std::endl;
 		Target target(targ_a, targ_z, targ_s, thickness);
 		m_recon.Init(target, theta, B, fpCal);
+		for(auto& table : ptables)
+			m_recon.AddPunchThruTable(table);
+		for(auto& table : etables)
+			m_recon.AddEnergyLossTable(table);
 		m_cuts.InitCuts(cuts);
 		m_cuts.InitEvent(m_eventPtr);
 
@@ -217,16 +252,6 @@ namespace SabreRecon {
 		float flush_frac = 0.01f;
 		uint64_t count = 0, flush_count = 0, flush_val = nevents*flush_frac;
 
-		//Analysis results data
-		ReconResult recon5Li, recon7Be, recon8Be, recon14N, recon9B;
-		TVector3 sabreCoords, b9Coords;
-		double relAngle;
-
-		//Temp
-		TFile* punchCutFile = TFile::Open("/Volumes/Wyndle/10B3He_May2022/cuts/punchProtons_relAngle.root");
-		TCutG* protonGate = (TCutG*) punchCutFile->Get("CUTG");
-		protonGate->SetName("protonPunchGate");
-
 		for(uint64_t i=0; i<nevents; i++)
 		{
 			tree->GetEntry(i);
@@ -241,9 +266,36 @@ namespace SabreRecon {
 			//Only analyze data that passes cuts, has sabre, and passes a weak threshold requirement
 			if(m_cuts.IsInside())
 			{
+				
 				FillHistogram1D({"xavg_gated","xavg_gated;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
+				
 				if(!m_eventPtr->sabre.empty() && m_eventPtr->sabre[0].ringE > s_weakSabreThreshold)
 				{
+					if(m_eventPtr->sabre[0].detID == 0 || m_eventPtr->sabre[0].detID == 1 || m_eventPtr->sabre[0].detID == 4)
+					{
+						RunDegradedSabre();
+					}
+					else
+					{
+						RunSabre();
+					}
+				}
+			}	
+		}
+		std::cout<<std::endl;
+		input->Close();
+		output->cd();
+		for(auto& gram : m_histoMap)
+			gram.second->Write(gram.second->GetName(), TObject::kOverwrite);
+		output->Close();
+	}
+
+	void Histogrammer::RunSabre()
+	{
+		static ReconResult recon5Li, recon7Be, recon8Be, recon14N, recon9B;
+		static TVector3 sabreCoords, b9Coords;
+		static double relAngle;
+
 		auto& biggestSabre = m_eventPtr->sabre[0];
 		recon9B = m_recon.RunFPResidExcitation(m_eventPtr->xavg, m_beamKE, {{5,10},{2,3},{3,4}});
 		recon5Li = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,10},{2,3},{2,4},{2,4}});
@@ -255,35 +307,27 @@ namespace SabreRecon {
 		relAngle = std::acos(b9Coords.Dot(sabreCoords)/(sabreCoords.Mag()*b9Coords.Mag()));
 
 		FillHistogram1D({"xavg_gated_sabre","xavg_gated_sabre;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
-					FillHistogram2D({"scintE_cathodeE","scintE_cathodeE;scintE;cathodeE",512,0,4096,512,0,4096},
+		FillHistogram2D({"scintE_cathodeE","scintE_cathodeE;scintE;cathodeE",512,0,4096,512,0,4096}, m_eventPtr->scintE, m_eventPtr->cathodeE);
-									m_eventPtr->scintE, m_eventPtr->cathodeE);
 		FillHistogram2D({"xavg_theta","xavg_theta;xavg;theta",600,-300.0,300.0,500,0.0,1.5}, m_eventPtr->xavg, m_eventPtr->theta);
 
 		FillHistogram1D({"ex_5Li", "ex_5Li;E_x(MeV);counts",3000,-5.0,25.0}, recon5Li.excitation);
 		FillHistogram1D({"ex_7Be", "ex_7Be;E_x(MeV);counts",3000,-20.0,10.0}, recon7Be.excitation);
 		FillHistogram1D({"ex_8Be", "ex_8Be;E_x(MeV);counts",3000,-5.0,25.0}, recon8Be.excitation);
 		FillHistogram1D({"ex_14N", "ex_14N;E_x(MeV);counts",3000,-20.0,10.0}, recon14N.excitation);
-					FillHistogram2D({"ex_14N_7Be","ex_14N_7Be;E_x 14N;E_x 7Be",500,-10.0,10.0,500,-10.,10.0}, recon14N.excitation,
+		FillHistogram2D({"ex_14N_7Be","ex_14N_7Be;E_x 14N;E_x 7Be",500,-10.0,10.0,500,-10.,10.0}, recon14N.excitation, recon7Be.excitation);
-									recon7Be.excitation);
+		FillHistogram2D({"sabreTheta_sabreE","sabreTheta_sabreE;#theta (deg); E(MeV)",180,0,180,400,0,20.0},sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE);
-					FillHistogram2D({"sabreTheta_sabreE","sabreTheta_sabreE;#theta (deg); E(MeV)",180,0,180,400,0,20.0},
+		FillHistogram2D({"xavg_sabreE","xavg_sabreE;xavg; E(MeV)",600,-300.0,300.0,400,0,20.0},m_eventPtr->xavg, biggestSabre.ringE);
-									sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE);
+		FillHistogram2D({"9Btheta_sabreTheta","9Btheta_sabreTheta;#theta_{9B};#theta_{SABRE}",180,0.0,180.0,180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg, sabreCoords.Theta()*s_rad2deg);
-					FillHistogram2D({"xavg_sabreE","xavg_sabreE;xavg; E(MeV)",600,-300.0,300.0,400,0,20.0},
+		FillHistogram2D({"sabreE_relAngle","sabreE_relAngle;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg,biggestSabre.ringE);
-									m_eventPtr->xavg, biggestSabre.ringE);
-					FillHistogram2D({"9Btheta_sabreTheta","9Btheta_sabreTheta;#theta_{9B};#theta_{SABRE}",180,0.0,180.0,
-									180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg, sabreCoords.Theta()*s_rad2deg);
-					FillHistogram2D({"sabreE_relAngle","sabreE_relAngle;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg,
-									biggestSabre.ringE);
 
 		if(m_eventPtr->xavg > -186.0 && m_eventPtr->xavg < -178.0 && (biggestSabre.detID == 2 || biggestSabre.detID == 3))
 		{
-						FillHistogram2D({"sabreE_sabreTheta_nub","sabreE_sabreTheta_nub;#theta (deg);E(MeV)",
+			FillHistogram2D({"sabreE_sabreTheta_nub","sabreE_sabreTheta_nub;#theta (deg);E(MeV)",180,0.0,180.0,400,0.0,20.0},sabreCoords.Theta()*s_rad2deg,biggestSabre.ringE);
-										180,0.0,180.0,400,0.0,20.0},sabreCoords.Theta()*s_rad2deg,biggestSabre.ringE);
+			FillHistogram2D({"sabreE_sabrePhi_nub","sabreE_sabreTheta_nub;#phi (deg);E(MeV)",360,0.0,360.0,400,0.0,20.0},Phi360(sabreCoords.Phi())*s_rad2deg,biggestSabre.ringE);
-						FillHistogram2D({"sabreE_sabrePhi_nub","sabreE_sabreTheta_nub;#phi (deg);E(MeV)",
-										360,0.0,360.0,400,0.0,20.0},Phi360(sabreCoords.Phi())*s_rad2deg,biggestSabre.ringE);
 		}
 
 		//Gate on reconstr. excitation structures; overlaping cases are possible!
-					if(recon5Li.excitation > -1.5 && recon5Li.excitation < 1.5)
+		if(recon5Li.excitation > -2.0 && recon5Li.excitation < 2.0)
 		{
 			FillHistogram1D({"xavg_gated5Ligs", "xavg_gated5Ligs;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
 		}
@@ -294,68 +338,70 @@ namespace SabreRecon {
 		if(recon7Be.excitation > -0.1 && recon7Be.excitation < 0.15)
 		{
 			FillHistogram1D({"xavg_gated7Begs", "xavg_gated7Begs;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
-						FillHistogram2D({"xavg_sabreE_7Begs","xavg_sabreE_7Begs;xavg;E(MeV)",600,-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg,
+			FillHistogram2D({"xavg_sabreE_7Begs","xavg_sabreE_7Begs;xavg;E(MeV)",600,-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE);
-										biggestSabre.ringE);
+			if(!(recon14N.excitation > -0.1 && recon14N.excitation < 2.0))
+				FillHistogram1D({"xavg_gated7Begs_reject14Ngs", "xavg_gated7Begs_reject14Ngs;xavg;counts",600, -300.0, 300.0}, m_eventPtr->xavg);
 			if(m_eventPtr->xavg > -186.0 && m_eventPtr->xavg < -178.0)
 			{
-							FillHistogram2D({"sabreE_sabreTheta_7begs_nub","sabreE_sabreTheta_7begs_nub;#theta (deg);E(MeV)",
+				FillHistogram2D({"sabreE_sabreTheta_7begs_nub","sabreE_sabreTheta_7begs_nub;#theta (deg);E(MeV)",180,0.0,180.0,400,0.0,20.0},sabreCoords.Theta()*s_rad2deg,biggestSabre.ringE);
-											180,0.0,180.0,400,0.0,20.0},sabreCoords.Theta()*s_rad2deg,biggestSabre.ringE);
+				FillHistogram2D({"sabreE_sabrePhi_7begs_nub","sabreE_sabreTheta_7begs_nub;#phi (deg);E(MeV)",360,0.0,360.0,400,0.0,20.0},Phi360(sabreCoords.Phi())*s_rad2deg,biggestSabre.ringE);
-							FillHistogram2D({"sabreE_sabrePhi_7begs_nub","sabreE_sabreTheta_7begs_nub;#phi (deg);E(MeV)",
-											360,0.0,360.0,400,0.0,20.0},Phi360(sabreCoords.Phi())*s_rad2deg,biggestSabre.ringE);
 			}
 		}
-					if(recon14N.excitation > -0.1 && recon14N.excitation < 0.1)
+		if(recon14N.excitation > -0.1 && recon14N.excitation < 0.2)
 		{
 			FillHistogram1D({"xavg_gated14Ngs", "xavg_gated14Ngs;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
 		}
-					
+		if(!(recon14N.excitation > -0.1 && recon14N.excitation < 0.2) && !(recon7Be.excitation > -0.1 && recon7Be.excitation < 0.15)
-					if(!(recon14N.excitation > -0.1 && recon14N.excitation < 0.1) && !(recon7Be.excitation > -0.1 && recon7Be.excitation < 0.15)
+			&& !(recon8Be.excitation > -0.1 && recon8Be.excitation < 0.1) && !(recon5Li.excitation > -2.0 && recon5Li.excitation < 2.0))
-						&& !(recon8Be.excitation > -0.1 && recon8Be.excitation < 0.1) && !(recon5Li.excitation > -1.5 && recon5Li.excitation < 1.5))
 		{
 			FillHistogram1D({"xavg_notGatedAllChannels", "xavg_notGatedAllChannels;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
 		}
+	}
 
-					//Degrader analysis... SABRE detectors 0, 1, 4 are covered with tantalum
+	void Histogrammer::RunDegradedSabre()
-					if(biggestSabre.detID == 0 || biggestSabre.detID == 1 || biggestSabre.detID == 4)
 	{
-						FillHistogram2D({"sabreE_sabreTheta_degraded", "sabreE_sabreTheta_degraded;#theta (deg);E(MeV)",
+		static ReconResult recon8Be, recon8BePunch, recon9Be, recon9BePunch, recon9B;
-										180,0.0,180.0,400,0.0,20.0}, sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE);
+		static TVector3 sabreCoords, b9Coords;
-						FillHistogram2D({"xavg_sabreE_degraded", "xavg_sabreE_degraded;xavg;E(MeV)",
+		static double relAngle;
-										600,0.-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE);
+
-						FillHistogram2D({"9Btheta_sabreTheta_degraded","9Btheta_sabreTheta_degraded;#theta_{9B};#theta_{SABRE}",180,0.0,180.0,
+		auto& biggestSabre = m_eventPtr->sabre[0];
-										180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg, sabreCoords.Theta()*s_rad2deg);
+
-						FillHistogram2D({"sabreE_relAngle_degraded","sabreE_relAngle_degraded;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg,
+		recon8Be = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,10},{2,3},{2,4},{1,1}});
-										biggestSabre.ringE);
+		recon8BePunch = m_recon.RunSabreExcitationPunchDegraded(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,10},{2,3},{2,4},{1,1}});
+		recon9Be = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,11},{2,3},{2,4},{1,1}});
+		recon9BePunch = m_recon.RunSabreExcitationPunchDegraded(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,11},{2,3},{2,4},{1,1}});
+		recon9B = m_recon.RunFPResidExcitation(m_eventPtr->xavg, m_beamKE, {{5,10},{2,3},{3,4}});
+		sabreCoords = m_recon.GetSabreCoordinates(biggestSabre);
+		b9Coords.SetMagThetaPhi(1.0, recon9B.residThetaLab, recon9B.residPhiLab);
+		relAngle = std::acos(b9Coords.Dot(sabreCoords)/(sabreCoords.Mag()*b9Coords.Mag()));
+
+		FillHistogram1D({"ex_8be_degraded","ex_8be_degraded; E_x(MeV); counts",3000,-5.0,25.0}, recon8Be.excitation);
+		FillHistogram2D({"xavg_ex8be_degraded","xavg_ex8be_degraded;xavg;E_x(MeV)",600,-300.0,300.0,300,-5.0,25.0}, m_eventPtr->xavg, recon8Be.excitation);
+		FillHistogram1D({"ex_8be_degradedPunched","ex_8be_degradedPunched; E_x(MeV); counts",3000,-15.0,15.0}, recon8BePunch.excitation);
+		FillHistogram2D({"xavg_ex8be_degradedPunched","xavg_ex8be_degradedPunched;xavg;E_x(MeV)",600,-300.0,300.0,300,-15.0,15.0}, m_eventPtr->xavg, recon8BePunch.excitation);
+		FillHistogram1D({"ex_9be_degraded","ex_9be_degraded; E_x(MeV); counts",3000,-5.0,25.0}, recon9Be.excitation);
+		FillHistogram2D({"xavg_ex9be_degraded","xavg_ex9be_degraded;xavg;E_x(MeV)",600,-300.0,300.0,300,-5.0,25.0}, m_eventPtr->xavg, recon9Be.excitation);
+		FillHistogram1D({"ex_9be_degradedPunched","ex_9be_degradedPunched; E_x(MeV); counts",3000,-15.0,15.0}, recon9BePunch.excitation);
+		FillHistogram2D({"xavg_ex9be_degradedPunched","xavg_ex9be_degradedPunched;xavg;E_x(MeV)",600,-300.0,300.0,300,-15.0,15.0}, m_eventPtr->xavg, recon9BePunch.excitation);
+		FillHistogram2D({"relAngle_recovSabreKE_8bePunchRecon","relAngle_recovSabreKe;#theta_{rel}(deg);Recovered KE (MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg,recon8BePunch.sabreRxnKE);
+		FillHistogram2D({"relAngle_recovSabreKE_9bePunchRecon","relAngle_recovSabreKe;#theta_{rel}(deg);Recovered KE (MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg,recon9BePunch.sabreRxnKE);
+		if(recon8Be.excitation > 2.2 && recon8Be.excitation < 3.8)
+		{
+			FillHistogram1D({"xavg_gated_8be1ex_degraded", "xavg_gated_8be1ex_degraded;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
+		}
+
+		FillHistogram2D({"sabreE_sabreTheta_degraded", "sabreE_sabreTheta_degraded;#theta (deg);E(MeV)",180,0.0,180.0,400,0.0,20.0}, sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE);
+		FillHistogram2D({"xavg_sabreE_degraded", "xavg_sabreE_degraded;xavg;E(MeV)", 600,0.-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE);
+		FillHistogram2D({"9Btheta_sabreTheta_degraded","9Btheta_sabreTheta_degraded;#theta_{9B};#theta_{SABRE}",180,0.0,180.0,180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg, sabreCoords.Theta()*s_rad2deg);
+		FillHistogram2D({"sabreE_relAngle_degraded","sabreE_relAngle_degraded;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg, biggestSabre.ringE);
 		if(biggestSabre.local_wedge != 0 && biggestSabre.local_wedge != 7 && biggestSabre.local_ring != 15 && biggestSabre.local_ring != 0) //Edges might not be degraded right
 		{
-							FillHistogram2D({"sabreE_sabreTheta_degraded_rejectEdge", "sabreE_sabreTheta_degraded;#theta (deg);E(MeV)",
+			FillHistogram2D({"sabreE_sabreTheta_degraded_rejectEdge", "sabreE_sabreTheta_degraded;#theta (deg);E(MeV)",180,0.0,180.0,400,0.0,20.0}, sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE);
-										180,0.0,180.0,400,0.0,20.0}, sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE);
+			FillHistogram2D({"xavg_sabreE_degraded_rejectEdge", "xavg_sabreE_degraded;xavg;E(MeV)",600,0.-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE);
-							FillHistogram2D({"xavg_sabreE_degraded_rejectEdge", "xavg_sabreE_degraded;xavg;E(MeV)",
-										600,0.-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE);
 			FillHistogram1D({"xavg_degraded_rejectEdge","xavg_degraded_rejectEdge;xavg",600,-300.0,300.0}, m_eventPtr->xavg);
-							FillHistogram2D({"9Btheta_sabreTheta_degraded_rejectEdge","9Btheta_sabreTheta_degraded_rejectEdge;#theta_{9B};#theta_{SABRE}",180,0.0,180.0,
+			FillHistogram2D({"9Btheta_sabreTheta_degraded_rejectEdge","9Btheta_sabreTheta_degraded_rejectEdge;#theta_{9B};#theta_{SABRE}",180,0.0,180.0,180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg,
-										180,0.0,180.0}, recon9B.residThetaLab*s_rad2deg, sabreCoords.Theta()*s_rad2deg);
+							  sabreCoords.Theta()*s_rad2deg);
-							FillHistogram2D({"sabreE_relAngle_degraded_rejectEdge","sabreE_relAngle_degraded_rejectEdge;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg,
+			FillHistogram2D({"sabreE_relAngle_degraded_rejectEdge","sabreE_relAngle_degraded_rejectEdge;#theta_{rel};E(MeV)",180,0.0,180.0,400,0.0,20.0},relAngle*s_rad2deg, biggestSabre.ringE);
-											biggestSabre.ringE);
-							
-							if(protonGate->IsInside(relAngle*s_rad2deg, biggestSabre.ringE))
-							{
-								FillHistogram2D({"xavg_sabreE_degraded_rejectEdge_pGate", "xavg_sabreE_degraded;xavg;E(MeV)",
-												600,0.-300.0,300.0,400,0.0,20.0}, m_eventPtr->xavg, biggestSabre.ringE);
-								FillHistogram1D({"xavg_degraded_rejectEdge_pGate","xavg_degraded_rejectEdge;xavg",600,-300.0,300.0}, m_eventPtr->xavg);
 		}
 	}
-					}
-				}
-			}	
-		}
-		std::cout<<std::endl;
-		input->Close();
-		output->cd();
-		for(auto& gram : m_histoMap)
-			gram.second->Write(gram.second->GetName(), TObject::kOverwrite);
-		protonGate->Write();
-		output->Close();
-		punchCutFile->Close();
-	}
 }

src/Histogrammer.h

@@ -41,6 +41,9 @@ namespace SabreRecon {
 		void Run();
 
 	private:
+		void RunSabre();
+		void RunDegradedSabre();
+
 		void ParseConfig(const std::string& name);
 		void FillHistogram1D(const Histogram1DParams& params, double value);
 		void FillHistogram2D(const Histogram2DParams& params, double valueX, double valueY);

src/Reconstructor.cpp

@@ -23,45 +23,162 @@ namespace SabreRecon {
 		m_focalPlane.Init({spsB, spsTheta, spsCal});
 		for(int i=0; i<5; i++)
 			m_sabreArray.emplace_back(SabreDetector::Parameters(s_phiDet[i], s_tiltAngle, s_zOffset, false, i));
+
+		//Setup intermediate energy loss layers
+		m_sabreDeadLayer.SetParameters({28}, {14}, {1}, s_sabreDeadlayerThickness);
+	}
+
+	void Reconstructor::AddEnergyLossTable(const std::string& filename)
+	{
+		m_elossTables.emplace_back(filename);
+	}
+
+	void Reconstructor::AddPunchThruTable(const std::string& filename)
+	{
+		m_punchTables.emplace_back(filename);
+	}
+
+	PunchTable::ElossTable* Reconstructor::GetElossTable(const NucID& projectile, const NucID& material)
+	{
+		MassLookup& masses = MassLookup::GetInstance();
+		std::string projString = masses.FindSymbol(projectile.Z, projectile.A);
+		std::string matString = masses.FindSymbol(material.Z, material.A) + "1"; //temp
+
+		for(auto& table : m_elossTables)
+		{
+			if(table.GetProjectile() == projString && table.GetMaterial() == matString)
+				return &table;
+		}
+
+		return nullptr;
+	}
+
+	PunchTable::PunchTable* Reconstructor::GetPunchThruTable(const NucID& projectile, const NucID& material)
+	{
+		MassLookup& masses = MassLookup::GetInstance();
+		std::string projString = masses.FindSymbol(projectile.Z, projectile.A);
+		std::string matString = masses.FindSymbol(material.Z, material.A) + "1"; //temp
+
+		for(auto& table : m_punchTables)
+		{
+			if(table.GetProjectile() == projString && table.GetMaterial() == matString)
+				return &table;
+		}
+
+		return nullptr;
 	}
 
 	TLorentzVector Reconstructor::GetSabre4Vector(const SabrePair& pair, double mass)
 	{
 		TVector3 coords;
 		TLorentzVector result;
-		double p, E;
+		double p, E, theta, phi;
-		if(pair.detID == 4)
+		//if(pair.detID == 4)
-			coords = m_sabreArray[4].GetHitCoordinates(15-pair.local_ring, pair.local_wedge);
+		//	coords = m_sabreArray[4].GetHitCoordinates(15-pair.local_ring, pair.local_wedge);
-		else
+		//else
 			coords = m_sabreArray[pair.detID].GetHitCoordinates(pair.local_ring, pair.local_wedge);
 		p = std::sqrt(pair.ringE*(pair.ringE + 2.0*mass));
 		E = pair.ringE + mass;
-		result.SetPxPyPzE(p*std::sin(coords.Theta())*std::cos(coords.Phi()),
+		theta = coords.Theta();
-						  p*std::sin(coords.Theta())*std::sin(coords.Phi()),
+		phi = coords.Phi();
-						  p*std::cos(coords.Theta()),
+		result.SetPxPyPzE(p*std::sin(theta)*std::cos(phi),
+						  p*std::sin(theta)*std::sin(phi),
+						  p*std::cos(theta),
 						  E);
 		return result;
 	}
 
 	TLorentzVector Reconstructor::GetSabre4VectorEloss(const SabrePair& pair, double mass, const NucID& id)
 	{
-		TVector3 coords;
+		TVector3 coords, sabreNorm;
 		TLorentzVector result;
-		double p, E, rxnKE;
+		double incidentAngle, p, E, rxnKE, theta, phi;
-		if(pair.detID == 4)
+
-			coords = m_sabreArray[4].GetHitCoordinates(15-pair.local_ring, pair.local_wedge);
+		//if(pair.detID == 4)
-		else
+		//	coords = m_sabreArray[4].GetHitCoordinates(15-pair.local_ring, pair.local_wedge);
+		//else
 			coords = m_sabreArray[pair.detID].GetHitCoordinates(pair.local_ring, pair.local_wedge);
-		rxnKE = pair.ringE + m_target.GetReverseEnergyLossFractionalDepth(id.Z, id.A, pair.ringE, coords.Theta(), 0.5);
+		sabreNorm = m_sabreArray[pair.detID].GetNormTilted();
-		p = std::sqrt(pair.ringE*(pair.ringE + 2.0*mass));
+		incidentAngle = std::acos(sabreNorm.Dot(coords)/(sabreNorm.Mag()*coords.Mag()));
-		E = pair.ringE + mass;
+
-		result.SetPxPyPzE(p*std::sin(coords.Theta())*std::cos(coords.Phi()),
+		rxnKE = pair.ringE + m_sabreDeadLayer.GetReverseEnergyLossTotal(id.Z, id.A, pair.ringE, incidentAngle);
-						  p*std::sin(coords.Theta())*std::sin(coords.Phi()),
+		rxnKE += m_target.GetReverseEnergyLossFractionalDepth(id.Z, id.A, rxnKE, coords.Theta(), 0.5);
-						  p*std::cos(coords.Theta()),
+		p = std::sqrt(rxnKE*(rxnKE + 2.0*mass));
+		E = rxnKE + mass;
+		theta = coords.Theta();
+		phi = coords.Phi();
+		result.SetPxPyPzE(p*std::sin(theta)*std::cos(phi),
+						  p*std::sin(theta)*std::sin(phi),
+						  p*std::cos(theta),
 						  E);
 		return result;
 	}
 
+	TLorentzVector Reconstructor::GetSabre4VectorElossPunchThru(const SabrePair& pair, double mass, const NucID& id)
+	{
+		TVector3 coords, sabreNorm;
+		TLorentzVector result;
+		double incidentAngle, p, E, rxnKE, theta, phi;
+
+		PunchTable::PunchTable* table = GetPunchThruTable(id, {14, 28});
+		if(table == nullptr)
+			return result;
+
+		//if(pair.detID == 4)
+		//	coords = m_sabreArray[4].GetHitCoordinates(15-pair.local_ring, pair.local_wedge);
+		//else
+			coords = m_sabreArray[pair.detID].GetHitCoordinates(pair.local_ring, pair.local_wedge);
+		sabreNorm = m_sabreArray[pair.detID].GetNormTilted();
+		incidentAngle = std::acos(sabreNorm.Dot(coords)/(sabreNorm.Mag()*coords.Mag()));
+		if(incidentAngle > M_PI/2.0)
+			incidentAngle = M_PI - incidentAngle;
+		rxnKE = table->GetInitialKineticEnergy(incidentAngle, pair.ringE);
+		if(rxnKE == 0.0)
+			return result;
+		rxnKE += m_target.GetReverseEnergyLossFractionalDepth(id.Z, id.A, rxnKE, coords.Theta(), 0.5);
+		p = std::sqrt(rxnKE*(rxnKE + 2.0*mass));
+		E = rxnKE + mass;
+		theta = coords.Theta();
+		phi = coords.Phi();
+		result.SetPxPyPzE(p*std::sin(theta)*std::cos(phi), p*std::sin(theta)*std::sin(phi), p*std::cos(theta), E);
+		return result;
+	}
+
+	TLorentzVector Reconstructor::GetSabre4VectorElossPunchThruDegraded(const SabrePair& pair, double mass, const NucID& id)
+	{
+		TVector3 coords, sabreNorm;
+		TLorentzVector result;
+		double incidentAngle, p, E, rxnKE, theta, phi;
+
+		PunchTable::PunchTable* ptable = GetPunchThruTable(id, {14, 28});
+		PunchTable::ElossTable* etable = GetElossTable(id, {73, 181});
+		if(ptable == nullptr || etable == nullptr)
+			return result;
+
+		//if(pair.detID == 4)
+		//	coords = m_sabreArray[4].GetHitCoordinates(15-pair.local_ring, pair.local_wedge);
+		//else
+			coords = m_sabreArray[pair.detID].GetHitCoordinates(pair.local_ring, pair.local_wedge);
+		sabreNorm = m_sabreArray[pair.detID].GetNormTilted();
+		incidentAngle = std::acos(sabreNorm.Dot(coords)/(sabreNorm.Mag()*coords.Mag()));
+		if(incidentAngle > M_PI/2.0)
+			incidentAngle = M_PI - incidentAngle;
+
+		rxnKE = ptable->GetInitialKineticEnergy(incidentAngle, pair.ringE);
+		if(rxnKE == pair.ringE)
+			return result;
+		rxnKE += etable->GetEnergyLoss(incidentAngle, rxnKE);
+		if(rxnKE == 0.0)
+			return result;
+		rxnKE += m_target.GetReverseEnergyLossFractionalDepth(id.Z, id.A, rxnKE, coords.Theta(), 0.5);
+		p = std::sqrt(rxnKE*(rxnKE + 2.0*mass));
+		E = rxnKE + mass;
+		theta = coords.Theta();
+		phi = coords.Phi();
+		result.SetPxPyPzE(p*std::sin(theta)*std::cos(phi), p*std::sin(theta)*std::sin(phi), p*std::cos(theta), E);
+		return result;
+	}
+
 	TLorentzVector Reconstructor::GetFP4VectorEloss(double xavg, double mass, const NucID& id)
 	{
 		TLorentzVector result;
@@ -236,6 +353,7 @@ namespace SabreRecon {
 		result.excitation = resid_vec.M() - massResid;
 
 		auto parent_vec = targ_vec + proj_vec;
+		result.sabreRxnKE = eject_vec.E() - massEject;
 		auto boost = parent_vec.BoostVector();
 		eject_vec.Boost(-1.0*boost);
 		result.ejectThetaCM = eject_vec.Theta();
@@ -280,6 +398,7 @@ namespace SabreRecon {
 		TLorentzVector decayFrag_vec = resid_vec - decayBreak_vec;
 
 		result.excitation = decayFrag_vec.M() - massDecayFrag;
+		result.sabreRxnKE = decayBreak_vec.E() - massDecayBreak;
 		auto boost = resid_vec.BoostVector();
 		decayBreak_vec.Boost(-1.0*boost);
 		result.ejectThetaCM = decayBreak_vec.Theta();
@@ -324,6 +443,105 @@ namespace SabreRecon {
 		TLorentzVector decayFrag_vec = resid_vec - decayBreak_vec;
 
 		result.excitation = decayFrag_vec.M() - massDecayFrag;
+		result.sabreRxnKE = eject_vec.E() - massEject;
+		auto boost = resid_vec.BoostVector();
+		decayBreak_vec.Boost(-1.0*boost);
+		result.ejectThetaCM = decayBreak_vec.Theta();
+		result.ejectPhiCM = decayBreak_vec.Phi();
+
+		return result;
+	}
+
+	ReconResult Reconstructor::RunSabreExcitationPunch(double xavg, double beamKE, const SabrePair& sabre, const std::vector<NucID>& nuclei)
+	{
+		ReconResult result;
+
+		NucID decayFrag;
+		decayFrag.Z = nuclei[0].Z + nuclei[1].Z - nuclei[2].Z - nuclei[3].Z;
+		decayFrag.A = nuclei[0].A + nuclei[1].A - nuclei[2].A - nuclei[3].A;
+
+		if(decayFrag.Z > decayFrag.A || decayFrag.A <= 0 || decayFrag.Z < 0)
+		{
+			std::cerr<<"Invalid reisdual nucleus at Reconstructor::RunSabreExcitation with Z: "<<decayFrag.Z<<" A: "<<decayFrag.A<<std::endl;
+			return result;
+		}
+
+		MassLookup& masses = MassLookup::GetInstance();
+		double massTarg = masses.FindMass(nuclei[0].Z, nuclei[0].A);
+		double massProj = masses.FindMass(nuclei[1].Z, nuclei[1].A);
+		double massEject = masses.FindMass(nuclei[2].Z, nuclei[2].A);
+		double massDecayBreak = masses.FindMass(nuclei[3].Z, nuclei[3].A);
+		double massDecayFrag = masses.FindMass(decayFrag.Z, decayFrag.A);
+
+		if(massTarg == 0.0 || massProj == 0.0 || massEject == 0.0 || massDecayBreak == 0.0 || massDecayFrag == 0.0)
+		{
+			std::cerr<<"Invalid nuclei at Reconstructor::RunSabreExcitation by mass!"<<std::endl;
+			return result;
+		}
+
+		TLorentzVector targ_vec;
+		targ_vec.SetPxPyPzE(0.0, 0.0, 0.0, massTarg);
+		auto proj_vec = GetProj4VectorEloss(beamKE, massProj, nuclei[1]);
+		auto eject_vec = GetFP4VectorEloss(xavg, massEject, nuclei[2]);
+		auto decayBreak_vec = GetSabre4VectorElossPunchThru(sabre, massDecayBreak, nuclei[3]);
+		if(decayBreak_vec.E() == 0.0)
+		{
+			return result;
+		}
+		TLorentzVector resid_vec = targ_vec + proj_vec - eject_vec;
+		TLorentzVector decayFrag_vec = resid_vec - decayBreak_vec;
+
+		result.excitation = decayFrag_vec.M() - massDecayFrag;
+		result.sabreRxnKE = decayBreak_vec.E() - massDecayBreak;
+		auto boost = resid_vec.BoostVector();
+		decayBreak_vec.Boost(-1.0*boost);
+		result.ejectThetaCM = decayBreak_vec.Theta();
+		result.ejectPhiCM = decayBreak_vec.Phi();
+
+		return result;
+	}
+
+	ReconResult Reconstructor::RunSabreExcitationPunchDegraded(double xavg, double beamKE, const SabrePair& sabre, const std::vector<NucID>& nuclei)
+	{
+		ReconResult result;
+
+		NucID decayFrag;
+		decayFrag.Z = nuclei[0].Z + nuclei[1].Z - nuclei[2].Z - nuclei[3].Z;
+		decayFrag.A = nuclei[0].A + nuclei[1].A - nuclei[2].A - nuclei[3].A;
+
+		if(decayFrag.Z > decayFrag.A || decayFrag.A <= 0 || decayFrag.Z < 0)
+		{
+			std::cerr<<"Invalid reisdual nucleus at Reconstructor::RunSabreExcitation with Z: "<<decayFrag.Z<<" A: "<<decayFrag.A<<std::endl;
+			return result;
+		}
+
+		MassLookup& masses = MassLookup::GetInstance();
+		double massTarg = masses.FindMass(nuclei[0].Z, nuclei[0].A);
+		double massProj = masses.FindMass(nuclei[1].Z, nuclei[1].A);
+		double massEject = masses.FindMass(nuclei[2].Z, nuclei[2].A);
+		double massDecayBreak = masses.FindMass(nuclei[3].Z, nuclei[3].A);
+		double massDecayFrag = masses.FindMass(decayFrag.Z, decayFrag.A);
+
+		if(massTarg == 0.0 || massProj == 0.0 || massEject == 0.0 || massDecayBreak == 0.0 || massDecayFrag == 0.0)
+		{
+			std::cerr<<"Invalid nuclei at Reconstructor::RunSabreExcitation by mass!"<<std::endl;
+			return result;
+		}
+
+		TLorentzVector targ_vec;
+		targ_vec.SetPxPyPzE(0.0, 0.0, 0.0, massTarg);
+		auto proj_vec = GetProj4VectorEloss(beamKE, massProj, nuclei[1]);
+		auto eject_vec = GetFP4VectorEloss(xavg, massEject, nuclei[2]);
+		auto decayBreak_vec = GetSabre4VectorElossPunchThruDegraded(sabre, massDecayBreak, nuclei[3]);
+		if(decayBreak_vec.E() == 0.0)
+		{
+			return result;
+		}
+		TLorentzVector resid_vec = targ_vec + proj_vec - eject_vec;
+		TLorentzVector decayFrag_vec = resid_vec - decayBreak_vec;
+
+		result.excitation = decayFrag_vec.M() - massDecayFrag;
+		result.sabreRxnKE = decayBreak_vec.E() - massDecayBreak;
 		auto boost = resid_vec.BoostVector();
 		decayBreak_vec.Boost(-1.0*boost);
 		result.ejectThetaCM = decayBreak_vec.Theta();
@@ -334,9 +552,9 @@ namespace SabreRecon {
 
 	TVector3 Reconstructor::GetSabreCoordinates(const SabrePair& pair)
 	{
-		if(pair.detID == 4)
+		//if(pair.detID == 4)
-			return m_sabreArray[4].GetHitCoordinates(15-pair.local_ring, pair.local_wedge);
+		//	return m_sabreArray[4].GetHitCoordinates(15-pair.local_ring, pair.local_wedge);
-		else
+		//else
 			return m_sabreArray[pair.detID].GetHitCoordinates(pair.local_ring, pair.local_wedge);
 	}
 }

src/Reconstructor.h

@@ -4,6 +4,8 @@
 #include <string>
 #include <vector>
 #include "EnergyLoss/Target.h"
+#include "EnergyLoss/ElossTable.h"
+#include "EnergyLoss/PunchTable.h"
 #include "CalDict/DataStructs.h"
 #include "TLorentzVector.h"
 #include "Detectors/SabreDetector.h"
@@ -13,13 +15,14 @@ namespace SabreRecon {
 
 	struct ReconResult
 	{
-		double excitation;
+		double excitation = -100.0;
-		double ejectThetaCM;
+		double sabreRxnKE = -100.0;
-		double ejectPhiCM;
+		double ejectThetaCM = -100.0;
-		double residThetaLab;
+		double ejectPhiCM = -100.0;
-		double residPhiLab;
+		double residThetaLab = -100.0;
-		double residThetaCM;
+		double residPhiLab = -100.0;
-		double residPhiCM;
+		double residThetaCM = -100.0;
+		double residPhiCM = -100.0;
 	};
 
 	struct NucID
@@ -46,6 +49,9 @@ namespace SabreRecon {
 
 		void Init(const Target& target, double spsTheta, double spsB, const std::vector<double>& spsCal);
 
+		void AddEnergyLossTable(const std::string& filename);
+		void AddPunchThruTable(const std::string& filename);
+
 		ReconResult RunThreeParticleExcitation(const SabrePair& p1, const SabrePair& p2, const SabrePair& p3, const std::vector<NucID>& nuclei);
 		ReconResult RunTwoParticleExcitation(const SabrePair& p1, const SabrePair& p2, const std::vector<NucID>& nuclei);
 		//nuclei: target, projectile, ejectile
@@ -57,23 +63,36 @@ namespace SabreRecon {
     	//nuclei: target, projectile, ejectile, decayFP
     	ReconResult RunSabreExcitationDetEject(double xavg, double beamKE, const SabrePair& sabre, const std::vector<NucID>& nuclei);
 
+		ReconResult RunSabreExcitationPunch(double xavg, double beamKE, const SabrePair& sabre, const std::vector<NucID>& nuclei);
+		ReconResult RunSabreExcitationPunchDegraded(double xavg, double beamKE, const SabrePair& sabre, const std::vector<NucID>& nuclei);
+
 		TVector3 GetSabreCoordinates(const SabrePair& pair);
     	
 
 	private:
 		TLorentzVector GetSabre4Vector(const SabrePair& pair, double mass);
     	TLorentzVector GetSabre4VectorEloss(const SabrePair& pair, double mass, const NucID& id);
+    	TLorentzVector GetSabre4VectorElossPunchThru(const SabrePair& pair, double mass, const NucID& id);
+    	TLorentzVector GetSabre4VectorElossPunchThruDegraded(const SabrePair& pair, double mass, const NucID& id);
     	TLorentzVector GetFP4VectorEloss(double xavg, double mass, const NucID& id);
     	TLorentzVector GetProj4VectorEloss(double beamKE, double mass, const NucID& id);
 
+		PunchTable::PunchTable* GetPunchThruTable(const NucID& projectile, const NucID& material);
+		PunchTable::ElossTable* GetElossTable(const NucID& projectile, const NucID& material);
+
     	std::vector<SabreDetector> m_sabreArray;
     	FocalPlaneDetector m_focalPlane;
     	Target m_target;
+		Target m_sabreDeadLayer;
+
+		std::vector<PunchTable::PunchTable> m_punchTables;
+		std::vector<PunchTable::ElossTable> m_elossTables;
 
     	//SABRE constants
     	static constexpr double s_phiDet[5] = { 306.0, 18.0, 234.0, 162.0, 90.0 };
     	static constexpr double s_tiltAngle = 40.0;
     	static constexpr double s_zOffset = -0.1245;
+		static constexpr double s_sabreDeadlayerThickness = 50.0 * 1.0e-7 * 2.3296 * 1.0e6; // 50 nm deadlayer -> ug/cm^2
 
     	//Kinematics constants
     	static constexpr double s_deg2rad = M_PI/180.0; //rad/deg