Added catima for better energy loss calculations

.gitmodules

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

CMakeLists.txt

@@ -8,6 +8,7 @@ set(MASK_INCLUDE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/include)
 
 set(CMAKE_CXX_STANDARD 17)
 
+add_subdirectory(src/vendor/catima)
 add_subdirectory(src/Mask)
 add_subdirectory(src/MaskApp)
 add_subdirectory(src/Detectors)

include/MassLookup.h

@@ -22,6 +22,7 @@ namespace Mask {
 	public:
 		~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);
 	
 		static MassLookup& GetInstance() {

include/SabreEfficiency.h

@@ -25,6 +25,7 @@ private:
     
 	Mask::Target deadlayer;
     Mask::Target sabre_eloss;
+    Mask::Target degrader;
     SabreDeadChannelMap dmap;
 
 	//Sabre constants
@@ -41,6 +42,7 @@ private:
     const double DEG2RAD = M_PI/180.0;
     static constexpr double DEADLAYER_THIN = 50 * 1e-7 * 2.3296 * 1e6; // ug/cm^2 (50 nm thick * density)
     static constexpr double SABRE_THICKNESS = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density)
+    static constexpr double DEGRADER_THICKNESS = 70.0 * 1.0e-4 * 16.69 * 1e6; //tantalum degrader (70 um thick) 
 
     const double ENERGY_THRESHOLD = 0.2; //in MeV
 

include/Target.h

@@ -17,6 +17,8 @@ Written by G.W. McCann Aug. 2020
 #include <vector>
 #include <cmath>
 #include "EnergyLoss.h"
+#include "catima/gwm_integrators.h"
+#include "MassLookup.h"
 
 namespace Mask {
 
@@ -39,7 +41,9 @@ namespace Mask {
 	
 	private:
 		EnergyLoss eloss;
+		catima::Material target_material;
 		double thickness;
+		double thickness_gcm2;
 		std::vector<int> Z, A, Stoich;
 	
 	};

include/Vec4.h

@@ -15,7 +15,7 @@ namespace Mask {
 	public:
 		Vec4();
 		Vec4(double px, double py, double pz, double E);
-		virtual ~Vec4();
+		virtual ~Vec4(); 
 		void SetVectorCartesian(double px, double py, double pz, double E);
 		void SetVectorSpherical(double theta, double phi, double p, double E);
 	
@@ -66,4 +66,4 @@ namespace Mask {
 
 }
 
-#endif
+#endif

src/Detectors/AnasenEfficiency.cpp

@@ -1,6 +1,7 @@
 #include "AnasenEfficiency.h"
 #include <fstream>
 #include <iomanip>
+#include <iostream>
 
 AnasenEfficiency::AnasenEfficiency() :
 	DetectorEfficiency(), det_silicon(si_thickness)

src/Detectors/CMakeLists.txt

@@ -14,6 +14,7 @@ target_sources(DetectEff PUBLIC
 
 target_link_libraries(DetectEff
     Mask
+    catima
 )
 
 set_target_properties(DetectEff PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MASK_BINARY_DIR})

src/Mask/CMakeLists.txt

@@ -26,4 +26,5 @@ target_sources(Mask PRIVATE
     Vec4.cpp
 )
 
+target_link_libraries(Mask catima)
 set_target_properties(Mask PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${MASK_LIBRARY_DIR})

src/Mask/Target.cpp

@@ -11,12 +11,14 @@ Written by G.W. McCann Aug. 2020
 
 */
 #include "Target.h"
+#include "catima/nucdata.h"
 
 namespace Mask {
 
 	/*Targets must be of known thickness*/
 	Target::Target(double thick) {
 		thickness = thick;
+		thickness_gcm2 = thickness*1.0e-6;
 	}
 	
 	Target::~Target() {}
@@ -26,8 +28,12 @@ namespace Mask {
 		Z = z;
 		A = a;
 		Stoich = stoich;
-		
+		MassLookup& masses = MassLookup::GetInstance();
 		eloss.SetTargetComponents(Z, A, Stoich);
+		for(size_t i=0; i<Z.size(); i++)
+		{
+			target_material.add_element(masses.FindMassU(Z[i], A[i]), Z[i], Stoich[i]);
+		}
 	}
 	
 	/*Element verification*/
@@ -44,8 +50,11 @@ namespace Mask {
 			return startEnergy;
 		else if (theta > M_PI/2.) 
 			theta = M_PI - theta;
-
+		catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0);
-		return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta)));
+		proj.T = startEnergy/proj.A;
+		target_material.thickness(thickness_gcm2/fabs(cos(theta)));
+		return catima::integrate_energyloss(proj, target_material);
+		//return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta)));
 	}
 
 	/*Calculates the energy loss for traveling some fraction through the target*/
@@ -56,7 +65,11 @@ namespace Mask {
 		else if (theta > M_PI/2.)
 			theta = M_PI-theta;
 
-		return eloss.GetEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta))));
+		catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0);
+		proj.T = finalEnergy/proj.A;
+		target_material.thickness(thickness_gcm2*percent_depth/fabs(cos(theta)));
+		return catima::integrate_energyloss(proj, target_material);
+		//return eloss.GetEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta))));
 	}
 	
 	/*Calculates reverse energy loss for travelling all the way through the target*/
@@ -66,7 +79,11 @@ namespace Mask {
 		else if (theta > M_PI/2.) 
 			theta = M_PI - theta;
 
-		return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta)));
+		catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0);
+		proj.T = finalEnergy/proj.A;
+		target_material.thickness(thickness_gcm2/fabs(cos(theta)));
+		return catima::reverse_integrate_energyloss(proj, target_material);
+		//return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta)));
 	}
 
 	/*Calculates the reverse energy loss for traveling some fraction through the target*/
@@ -76,8 +93,11 @@ namespace Mask {
 			return finalEnergy;
 		else if (theta > M_PI/2.)
 			theta = M_PI-theta;
-
+		catima::Projectile proj(MassLookup::GetInstance().FindMassU(zp, ap), zp, 0.0, 0.0);
-		return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta))));
+		proj.T = finalEnergy/proj.A;
+		target_material.thickness(thickness_gcm2*percent_depth/fabs(cos(theta)));
+		return catima::reverse_integrate_energyloss(proj, target_material);
+		//return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness*percent_depth/(std::fabs(std::cos(theta))));
 	}
 
 }

src/Plotters/ROOT/RootPlotter.cpp

@@ -1,5 +1,6 @@
 #include "RootPlotter.h"
 #include <TFile.h>
+#include <TVector3.h>
 
 #include <iostream>
 
@@ -57,12 +58,19 @@ void RootPlotter::FillCorrelations(const Mask::MaskFileData& data, Mask::RxnType
 
 	if(type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn)
 	{
+		TVector3 p1, p2;
+		p1.SetMagThetaPhi(1.0, data.theta[3], data.phi[3]);
+		p2.SetMagThetaPhi(1.0, data.theta[4], data.phi[4]);
+		double theta_resid_break1 = std::acos(p1.Dot(p2));
 		std::string theta_break1_theta_break2_name = "theta_break1_theta_break2_cor";
 		std::string theta_break1_theta_break2_title = theta_break1_theta_break2_name + ";#theta_{lab} Breakup1 (deg);#theta_{lab} Breakup2 (deg)";
 		MyFill(theta_break1_theta_break2_name, theta_break1_theta_break2_title, data.theta[4]*rad2deg, data.theta[5]*rad2deg, 4);
 		std::string theta_resid_theta_break1_name = "theta_resid_theta_break1_cor";
 		std::string theta_resid_theta_break1_title = theta_resid_theta_break1_name + ";#theta_{lab} Residual (deg);#theta_{lab} Breakup1 (deg)";
 		MyFill(theta_resid_theta_break1_name, theta_resid_theta_break1_title, data.theta[3]*rad2deg, data.theta[4]*rad2deg, 4);
+		std::string ke_break1_theta_rel_name = "ke_break1_theta_rel";
+		std::string ke_break1_theta_rel_title = ke_break1_theta_rel_name + ";#theta_{resid-break1};KE_{break1} (MeV)";
+		MyFill(ke_break1_theta_rel_name, ke_break1_theta_rel_title, theta_resid_break1*rad2deg, data.KE[4], 4);
 	}
 	if(type == Mask::RxnType::ThreeStepRxn)
 	{
@@ -86,9 +94,16 @@ void RootPlotter::FillCorrelationsDetected(const Mask::MaskFileData& data, Mask:
 	}
 	if((type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) && data.detect_flag[4])
 	{
+		TVector3 p1, p2;
+		p1.SetMagThetaPhi(1.0, data.theta[3], data.phi[3]);
+		p2.SetMagThetaPhi(1.0, data.theta[4], data.phi[4]);
+		double theta_resid_break1 = std::acos(p1.Dot(p2));
 		std::string theta_resid_theta_break1_name = "theta_resid_theta_break1_cor_detected";
 		std::string theta_resid_theta_break1_title = theta_resid_theta_break1_name + ";#theta_{lab} Residual (deg);#theta_{lab} Breakup1 (deg)";
 		MyFill(theta_resid_theta_break1_name, theta_resid_theta_break1_title, data.theta[3]*rad2deg, data.theta[4]*rad2deg, 4);
+		std::string ke_break1_theta_rel_name = "ke_break1_theta_rel_detected";
+		std::string ke_break1_theta_rel_title = ke_break1_theta_rel_name + ";#theta_{resid-break1};KE_{break1} (MeV)";
+		MyFill(ke_break1_theta_rel_name, ke_break1_theta_rel_title, theta_resid_break1*rad2deg, data.KE[4], 4);
 	}
 
 	if((type == Mask::RxnType::TwoStepRxn || type == Mask::RxnType::ThreeStepRxn) && data.detect_flag[4] && data.detect_flag[5])

src/vendor/catima

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