Integrated catima for more advanced energy loss.

.gitignore

@@ -1,5 +1,6 @@
 *.o
 *.so
+*.dylib
 *.sublime-project
 *.sublime-workspace
 *.code-workspace

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "src/vendor/catima"]
+	path = src/vendor/catima
+	url = https://github.com/gwm17/catima.git

CMakeLists.txt

@@ -7,4 +7,5 @@ set(SABRERECON_LIBRARY_DIR ${CMAKE_CURRENT_SOURCE_DIR}/lib)
 
 set(CMAKE_CXX_STANDARD 17)
 
+add_subdirectory(src/vendor/catima)
 add_subdirectory(src)

input.txt

@@ -1,6 +1,6 @@
 begin_data
-	input /Volumes/Wyndle/10B3He_May2022/calibrated/run_53_149_cal.root
+	input /Volumes/Wyndle/10B3He_May2022/calibrated/run_53_143_cal.root
-	output /Volumes/Wyndle/10B3He_May2022/histograms/calHistograms_53_149.root
+	output /Volumes/Wyndle/10B3He_May2022/histograms/calHistograms_53_143_newEloss.root
 	beamKE(MeV) 24.0
 end_data
 begin_reconstructor
@@ -8,14 +8,14 @@ begin_reconstructor
 		B 8.759
 		theta 15.0
 		begin_fpcal
-			78.5352
+			78.5946
-			0.038
+			0.0382
 		end_fpcal
 	end_focalplane
 	begin_target
 		thickness 74.0
 		begin_elements
-			5 1
+			5 10 1
 		end_elements
 	end_target
 end_reconstructor

src/CMakeLists.txt

@@ -30,6 +30,7 @@ target_sources(SabreRecon PUBLIC
 	)
 target_link_libraries(SabreRecon
 	CalDict
+	catima
 	${ROOT_LIBRARIES}
 	)
 set_target_properties(SabreRecon PROPERTIES

src/EnergyLoss/Target.cpp

@@ -23,10 +23,10 @@ namespace SabreRecon {
 	}
 
 	/*Targets must be of known thickness*/
-	Target::Target(const std::vector<int>& z, const std::vector<int>& stoich, double thick) :
+	Target::Target(const std::vector<int>& a, const std::vector<int>& z, const std::vector<int>& stoich, double thick) :
 		m_isValid(false)
 	{
-		SetParameters(z, stoich, thick);
+		SetParameters(a, z, stoich, thick);
 	}
 	
 	Target::~Target()
@@ -34,7 +34,7 @@ namespace SabreRecon {
 	}
 	
 	/*Set target elements of given Z, A, S*/
-	void Target::SetParameters(const std::vector<int>& z, const std::vector<int>& stoich, double thick)
+	void Target::SetParameters(const std::vector<int>& a, const std::vector<int>& z, const std::vector<int>& stoich, double thick)
 	{
 		m_params.ZT = z;
 		double denom = 0;
@@ -43,6 +43,13 @@ namespace SabreRecon {
 		for(auto& s : stoich)
 			m_params.composition.push_back(s/denom);
 		m_totalThickness = thick;
+		m_totalThickness_gcm2 = m_totalThickness*1.0e-6;
+
+		auto& masses = MassLookup::GetInstance();
+		for(size_t i=0; i<z.size(); i++)
+		{
+			m_material.add_element(masses.FindMassU(z[i], a[i]), z[i], stoich[i]);
+		}
 		m_isValid = true;
 	}
 	
@@ -60,7 +67,14 @@ namespace SabreRecon {
 		m_params.energy = startEnergy;
 		m_params.thickness = m_totalThickness/(std::fabs(std::cos(theta)));
 
-		return EnergyLoss::GetEnergyLoss(m_params);
+		m_projectile.A = MassLookup::GetInstance().FindMassU(zp, ap);
+		m_projectile.Z = zp;
+		m_projectile.Q = zp;
+		m_projectile.T = startEnergy/m_projectile.A;
+		m_material.thickness(m_totalThickness_gcm2/(std::fabs(std::cos(theta))));
+
+		return catima::integrate_energyloss(m_projectile, m_material);
+		//return EnergyLoss::GetEnergyLoss(m_params);
 	}
 
 	/*Calculates the energy loss for traveling some fraction through the target*/
@@ -77,7 +91,15 @@ namespace SabreRecon {
 		m_params.energy = startEnergy;
 		m_params.thickness = m_totalThickness*percent_depth/(std::fabs(std::cos(theta)));
 
-		return EnergyLoss::GetEnergyLoss(m_params);
+		m_projectile.A = MassLookup::GetInstance().FindMassU(zp, ap);
+		m_projectile.Z = zp;
+		m_projectile.Q = zp;
+		m_projectile.T = startEnergy/m_projectile.A;
+		m_material.thickness(m_totalThickness_gcm2*percent_depth/(std::fabs(std::cos(theta))));
+
+		return catima::integrate_energyloss(m_projectile, m_material);
+
+		//return EnergyLoss::GetEnergyLoss(m_params);
 	}
 	
 	/*Calculates reverse energy loss for travelling all the way through the target*/
@@ -94,7 +116,15 @@ namespace SabreRecon {
 		m_params.energy = finalEnergy;
 		m_params.thickness = m_totalThickness/(std::fabs(std::cos(theta)));
 
-		return EnergyLoss::GetReverseEnergyLoss(m_params);
+		m_projectile.A = MassLookup::GetInstance().FindMassU(zp, ap);
+		m_projectile.Z = zp;
+		m_projectile.Q = zp;
+		m_projectile.T = finalEnergy/m_projectile.A;
+		m_material.thickness(m_totalThickness_gcm2/(std::fabs(std::cos(theta))));
+
+		return catima::reverse_integrate_energyloss(m_projectile, m_material);
+
+		//return EnergyLoss::GetReverseEnergyLoss(m_params);
 	}
 
 	/*Calculates the reverse energy loss for traveling some fraction through the target*/
@@ -111,7 +141,15 @@ namespace SabreRecon {
 		m_params.energy = finalEnergy;
 		m_params.thickness = m_totalThickness*percent_depth/(std::fabs(std::cos(theta)));
 
-		return EnergyLoss::GetReverseEnergyLoss(m_params);
+		m_projectile.A = MassLookup::GetInstance().FindMassU(zp, ap);
+		m_projectile.Z = zp;
+		m_projectile.Q = zp;
+		m_projectile.T = finalEnergy/m_projectile.A;
+		m_material.thickness(m_totalThickness_gcm2*percent_depth/(std::fabs(std::cos(theta))));
+
+		return catima::reverse_integrate_energyloss(m_projectile, m_material);
+
+		//return EnergyLoss::GetReverseEnergyLoss(m_params);
 	}
 
 }

src/EnergyLoss/Target.h

@@ -16,6 +16,7 @@ Written by G.W. McCann Aug. 2020
 #include <string>
 #include <vector>
 #include "EnergyLoss.h"
+#include "catima/gwm_integrators.h"
 
 namespace SabreRecon {
 
@@ -24,10 +25,10 @@ namespace SabreRecon {
 
 	public:
 		Target();
-	 	Target(const std::vector<int>& z, const std::vector<int>& stoich, double thick);
+	 	Target(const std::vector<int>& a, const std::vector<int>& z, const std::vector<int>& stoich, double thick);
 	 	~Target();
 
-	 	void SetParameters(const std::vector<int>& z, const std::vector<int>& stoich, double thick);
+	 	void SetParameters(const std::vector<int>& a, const std::vector<int>& z, const std::vector<int>& stoich, double thick);
 	 	double GetEnergyLossTotal(int zp, int ap, double startEnergy, double angle);
 	 	double GetReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double angle);
 	 	double GetEnergyLossFractionalDepth(int zp, int ap, double startEnergy, double angle, double percent_depth);
@@ -39,7 +40,10 @@ namespace SabreRecon {
 	
 	private:
 		EnergyLoss::Parameters m_params;
+		catima::Material m_material;
+		catima::Projectile m_projectile;
 		double m_totalThickness;
+		double m_totalThickness_gcm2;
 		bool m_isValid;
 	};
 

src/Histogrammer.cpp

@@ -40,6 +40,7 @@ namespace SabreRecon {
 
 		double thickness;
 		std::vector<int> targ_z;
+		std::vector<int> targ_a;
 		std::vector<int> targ_s;
 
 		std::vector<ReconCut> cuts;
@@ -105,8 +106,10 @@ namespace SabreRecon {
 							{
 								targ_z.push_back(std::stoi(junk));
 								input>>junk;
+								targ_a.push_back(std::stoi(junk));
+								input>>junk;
 								targ_s.push_back(std::stoi(junk));
-								std::cout<<"e"<<targ_s.size()-1<<": ("<<targ_z[targ_z.size()-1]<<","<<targ_s[targ_z.size()-1]<<") ";
+								std::cout<<"e"<<targ_s.size()-1<<": ("<<targ_z[targ_z.size()-1]<<","<<targ_a[targ_z.size()-1]<<","<<targ_s[targ_z.size()-1]<<") ";
 							}
 						}
 						std::cout<<std::endl;
@@ -143,7 +146,7 @@ namespace SabreRecon {
 
 		//init resources
 		std::cout<<"Initializing resources..."<<std::endl;
-		Target target(targ_z, targ_s, thickness);
+		Target target(targ_a, targ_z, targ_s, thickness);
 		m_recon.Init(target, theta, B, fpCal);
 		m_cuts.InitCuts(cuts);
 		m_cuts.InitEvent(m_eventPtr);
@@ -215,11 +218,12 @@ namespace SabreRecon {
 		uint64_t count = 0, flush_count = 0, flush_val = nevents*flush_frac;
 
 		//Analysis results data
-		ReconResult recon5Li, recon7Be, recon8Be, recon14N;
+		ReconResult recon5Li, recon7Be, recon8Be, recon14N, recon9B;
-		TVector3 sabreCoords;
+		TVector3 sabreCoords, b9Coords;
+		double relAngle;
 
 		//Temp
-		TFile* punchCutFile = TFile::Open("/Volumes/Wyndle/10B3He_May2022/cuts/protonPunchGate_strict.root");
+		TFile* punchCutFile = TFile::Open("/Volumes/Wyndle/10B3He_May2022/cuts/punchProtons_relAngle.root");
 		TCutG* protonGate = (TCutG*) punchCutFile->Get("CUTG");
 		protonGate->SetName("protonPunchGate");
 
@@ -241,11 +245,14 @@ namespace SabreRecon {
 				if(!m_eventPtr->sabre.empty() && m_eventPtr->sabre[0].ringE > s_weakSabreThreshold)
 				{
 					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}});
 					recon8Be = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,10},{2,3},{2,4},{1,1}});
 					recon7Be = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{5,10},{2,3},{2,4},{1,2}});
 					recon14N = m_recon.RunSabreExcitation(m_eventPtr->xavg, m_beamKE, biggestSabre, {{8,16},{2,3},{2,4},{1,1}});
 					sabreCoords = m_recon.GetSabreCoordinates(biggestSabre);
+					b9Coords.SetMagThetaPhi(1.0, recon9B.residThetaLab, recon9B.residPhiLab);
+					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},
@@ -262,6 +269,10 @@ namespace SabreRecon {
 									sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE);
 					FillHistogram2D({"xavg_sabreE","xavg_sabreE;xavg; E(MeV)",600,-300.0,300.0,400,0,20.0},
 									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))
 					{
@@ -297,9 +308,11 @@ namespace SabreRecon {
 					{
 						FillHistogram1D({"xavg_gated14Ngs", "xavg_gated14Ngs;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
 					}
-					else
+					
+					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 > -1.5 && recon5Li.excitation < 1.5))
 					{
-						FillHistogram1D({"xavg_notGated14Ngs", "xavg_notGated14Ngs;xavg;counts",600,-300.0,300.0}, m_eventPtr->xavg);
+						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
@@ -309,6 +322,10 @@ namespace SabreRecon {
 										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)",
@@ -316,7 +333,12 @@ namespace SabreRecon {
 							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);
-							if(protonGate->IsInside(sabreCoords.Theta()*s_rad2deg, biggestSabre.ringE))
+							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, 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,
+											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);

src/MassLookup.h

@@ -24,6 +24,7 @@ namespace SabreRecon {
 		MassLookup();
 		~MassLookup();
 		double FindMass(int Z, int A);
+		double FindMassU(int Z, int A) { return FindMass(Z, A)/u_to_mev; }
 		std::string FindSymbol(int Z, int A);
 	
 		inline static MassLookup& GetInstance() { return *s_instance; }

src/Reconstructor.cpp

@@ -187,12 +187,15 @@ namespace SabreRecon {
 		auto resid_vec = targ_vec + proj_vec - eject_vec;
 
 		result.excitation = resid_vec.M() - massResid;
+		result.residThetaLab = resid_vec.Theta();
+		result.residPhiLab = resid_vec.Phi();
 
 		auto parent_vec = targ_vec + proj_vec;
 		auto boost = parent_vec.BoostVector();
 		eject_vec.Boost(-1.0*boost);
-		result.theta_cm = eject_vec.Theta();
+
-		result.phi_cm = eject_vec.Phi();
+		result.ejectThetaCM = eject_vec.Theta();
+		result.ejectPhiCM = eject_vec.Phi();
 
 		return result;
 	}
@@ -235,8 +238,8 @@ namespace SabreRecon {
 		auto parent_vec = targ_vec + proj_vec;
 		auto boost = parent_vec.BoostVector();
 		eject_vec.Boost(-1.0*boost);
-		result.theta_cm = eject_vec.Theta();
+		result.ejectThetaCM = eject_vec.Theta();
-		result.phi_cm = eject_vec.Phi();
+		result.ejectPhiCM = eject_vec.Phi();
 
 		return result;
 	}
@@ -279,8 +282,8 @@ namespace SabreRecon {
 		result.excitation = decayFrag_vec.M() - massDecayFrag;
 		auto boost = resid_vec.BoostVector();
 		decayBreak_vec.Boost(-1.0*boost);
-		result.theta_cm = decayBreak_vec.Theta();
+		result.ejectThetaCM = decayBreak_vec.Theta();
-		result.phi_cm = decayBreak_vec.Phi();
+		result.ejectPhiCM = decayBreak_vec.Phi();
 
 		return result;
 	}
@@ -323,8 +326,8 @@ namespace SabreRecon {
 		result.excitation = decayFrag_vec.M() - massDecayFrag;
 		auto boost = resid_vec.BoostVector();
 		decayBreak_vec.Boost(-1.0*boost);
-		result.theta_cm = decayBreak_vec.Theta();
+		result.ejectThetaCM = decayBreak_vec.Theta();
-		result.phi_cm = decayBreak_vec.Phi();
+		result.ejectPhiCM = decayBreak_vec.Phi();
 
 		return result;
 	}

src/Reconstructor.h

@@ -14,8 +14,12 @@ namespace SabreRecon {
 	struct ReconResult
 	{
 		double excitation;
-		double theta_cm;
+		double ejectThetaCM;
-		double phi_cm;
+		double ejectPhiCM;
+		double residThetaLab;
+		double residPhiLab;
+		double residThetaCM;
+		double residPhiCM;
 	};
 
 	struct NucID

src/vendor/catima

@@ -0,0 +1 @@
+Subproject commit 629690a6f90ad9d343e02f72d4fad860fce1ed6d