Small modifications for slightly increased performance and bugfixes

.gitignore

@@ -3,6 +3,7 @@ bin/
 lib/
 include/
 .vscode/
+tables/
 
 *.o
 *.so
@@ -15,4 +16,4 @@ include/
 .DS_Store
 
 
-!.gitignore
+!.gitignore

CMakeLists.txt

@@ -2,7 +2,14 @@ cmake_minimum_required(VERSION 3.16)
 
 project(PunchTable)
 
-set(CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 17)
+
+if(NOT CMAKE_BUILD_TYPE  STREQUAL "Debug")
+    set(CMAKE_BUILD_TYPE "Release")
+    MESSAGE(STATUS "Build type Release")
+else()
+    MESSAGE(STATUS "Build type Debug")
+endif()
 
 set(PUNCHTABLE_BINARY_DIR ${CMAKE_CURRENT_SOURCE_DIR}/bin)
 

src/CMakeLists.txt

@@ -10,6 +10,7 @@ target_sources(PunchTable PRIVATE
     PunchTable.cpp
     GenerateTable.h
     GenerateTable.cpp
+    main.cpp
 )
 
 target_link_libraries(PunchTable catima)

src/CubicSpline.cpp

@@ -190,7 +190,7 @@ namespace PunchTable {
 			}
 			else if (i == (m_splines.size() -1))
 			{
-				//std::cerr<<"Error at CubicSpline::Evaluate! Input x value: "<<x<<" is not within the spline range min: "<<splines[0].x1<<" max: "<<splines[splines.size()-1].x2<<std::endl;
+				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;
 			}
 		}

src/GenerateTable.cpp

@@ -1,11 +1,183 @@
 #include "GenerateTable.h"
+#include "MassLookup.h"
 #include <fstream>
+#include <iostream>
 #include <cmath>
 
+#include "catima/gwm_integrators.h"
+
 namespace PunchTable {
 
     void GenerateTable(const TableParameters& params)
     {
+        static double s_deg2rad = M_PI/180.0; //deg -> radians
+        static double s_energyPrecision = 0.001; // keV precision
+        static double s_ug2g = 1.0e-6; //ug -> g
+        MassLookup& masses = MassLookup::GetInstance();
 
+        catima::Projectile projectile(masses.FindMassU(params.projectileZ, params.projectileA), params.projectileZ, 0.0, 0.0);
+        std::vector<catima::Material> materials;
+
+        if(params.targetZ.size() != params.targetA.size() ||
+           params.targetZ.size() != params.targetS.size() ||
+           params.targetZ.size() != params.targetThickness.size())
+        {
+            std::cerr<<"ERR --  Invalid target parameters passed to GenerateTable. Mismatching target arrays."<<std::endl;
+            return;
+        }
+
+        for(size_t i=0; i<params.targetZ.size(); i++)
+        {
+            auto& tZ = params.targetZ[i];
+            auto& tA = params.targetA[i];
+            auto& tS = params.targetS[i];
+            materials.emplace_back();
+            for(size_t j=0; j<tZ.size(); j++)
+            {
+                materials[i].add_element(masses.FindMassU(tZ[j], tA[j]), tZ[j], tS[j]);
+            }
+        }
+
+        size_t depLayer = params.targetZ.size()-1;
+        size_t thetaBins = (params.maxTheta - params.minTheta)/params.stepTheta;
+        size_t energyBins = (params.maxKE - params.minKE)/params.stepKE;
+
+        std::ofstream output(params.filename);
+        if(!output.is_open())
+        {
+            std::cerr<<"ERR -- Unable to open output file "<<params.filename<<" at GenerateTable"<<std::endl;
+            return;
+        }
+		output<<std::setprecision(5);
+
+		output<<"Materials: "<<std::endl;;
+		for(size_t i=0; i<depLayer; i++)
+        {
+            output<<"\tGoing through: ";
+            for(size_t j=0; j<params.targetZ[i].size(); j++)
+			    output<<masses.FindSymbol(params.targetZ[i][j], params.targetA[i][j])<<params.targetS[i][j];
+            output<<" ("<<params.targetThickness[i]<<" ug/cm2)"<<std::endl;
+        }
+        output<<"\tDeposting into: ";
+        for(size_t i=0; i<params.targetZ[depLayer].size(); i++)
+            output<<masses.FindSymbol(params.targetZ[depLayer][i], params.targetA[depLayer][i])<<params.targetS[depLayer][i];
+        output<<" ("<<params.targetThickness[depLayer]<<" ug/cm2)"<<std::endl;
+		output<<"---------------------------------"<<std::endl;
+
+		output<<"IncidentAngleRange(deg): "<<params.minTheta<<" to "<<params.maxTheta<<" stepSize: "<<params.stepTheta<<std::endl;
+		output<<"---------------------------------"<<std::endl;
+
+		output<<std::setw(16)<<"Energy Deposited"<<"|"<<std::setw(16)<<"Intial Energy"<<std::endl;
+		output<<"---------------------------------"<<std::endl;
+
+        std::vector<double> energyInData;
+        std::vector<double> energyDepositedData;
+
+        double theta, ke;
+        double edep;
+        double totalDep;
+        bool stopped;
+
+        size_t totalBins = thetaBins*energyBins;
+        double flush_percent=0.01;
+        size_t count=0, flush_val = flush_percent*totalBins, flush_count=0;
+        for(size_t i=0; i < thetaBins; i++)
+        {
+            theta = (params.minTheta + i*params.stepTheta)*s_deg2rad;
+            energyInData.clear();
+            energyDepositedData.clear();
+            for(size_t j=0; j<energyBins; j++)
+            {
+                count++;
+                if(count == flush_val)
+                {
+                    count=0;
+                    flush_count++;
+                    std::cout<<"\rPercent of data generated: "<<flush_count*flush_percent*100.0<<"%"<<std::flush;
+                }
+                ke = params.maxKE - j*params.stepKE;
+                projectile.T = ke/projectile.A;
+                totalDep = 0.0;
+                stopped = false;
+                for(size_t k=0; k<depLayer; k++)
+                {
+                    materials[k].thickness(params.targetThickness[k]*s_ug2g/std::fabs(std::cos(theta)));
+                    totalDep += catima::integrate_energyloss(projectile, materials[k]);
+                    if((ke - totalDep) < s_energyPrecision)
+                    {
+                        stopped = true;
+                        break;
+                    }
+                }
+
+                if(stopped)
+                    continue;
+
+                materials[depLayer].thickness(params.targetThickness[depLayer]*s_ug2g/std::fabs(std::cos(theta)));
+                edep = catima::integrate_energyloss(projectile, materials[depLayer]);
+                if(std::abs(ke-edep) > s_energyPrecision)
+                {
+                    energyInData.push_back(ke);
+                    energyDepositedData.push_back(edep);
+                }
+            }
+
+            output<<"begin_theta "<<theta<<std::endl;
+			for(size_t j=0; j<energyDepositedData.size(); j++) 
+			{
+				output<<std::setw(16)<<energyDepositedData[j]<<" "<<std::setw(16)<<energyInData[j]<<std::endl;
+			}
+			output<<"end_theta"<<std::endl;   
+        }
+
+        
+		output.close();
+    }
+
+    double GetEnergyDeposited(const TableParameters& params, double thetaIncident, double initialKE)
+    {
+        static double s_energyPrecision = 0.001; // keV precision
+        static double s_ug2g = 1.0e-6; //ug -> g
+        MassLookup& masses = MassLookup::GetInstance();
+
+
+        catima::Projectile projectile(masses.FindMassU(params.projectileZ, params.projectileA), params.projectileZ, 0.0, 0.0);
+        std::vector<catima::Material> materials;
+
+         if(params.targetZ.size() != params.targetA.size() ||
+           params.targetZ.size() != params.targetS.size() ||
+           params.targetZ.size() != params.targetThickness.size())
+        {
+            std::cerr<<"ERR --  Invalid target parameters passed to GenerateTable. Mismatching target arrays."<<std::endl;
+            return 0.0;
+        }
+
+        for(size_t i=0; i<params.targetZ.size(); i++)
+        {
+            auto& tZ = params.targetZ[i];
+            auto& tA = params.targetA[i];
+            auto& tS = params.targetS[i];
+            materials.emplace_back();
+            for(size_t j=0; j<tZ.size(); j++)
+            {
+                materials[i].add_element(masses.FindMassU(tZ[j], tA[j]), tZ[j], tS[j]);
+            }
+        }
+
+        size_t depLayer = params.targetZ.size() - 1;
+        projectile.T = initialKE/projectile.A;
+        double totalDep = 0.0;
+        for(size_t k=0; k<depLayer; k++)
+        {
+            materials[k].thickness(params.targetThickness[k]*s_ug2g/std::fabs(std::cos(thetaIncident)));
+            totalDep += catima::integrate_energyloss(projectile, materials[k]);
+            if((initialKE - totalDep) < s_energyPrecision)
+            {
+                return 0.0;
+            }
+        }
+
+        materials[depLayer].thickness(params.targetThickness[depLayer]*s_ug2g/std::fabs(std::cos(thetaIncident)));
+        return catima::integrate_energyloss(projectile, materials[depLayer]);
     }
 }

src/GenerateTable.h

@@ -8,22 +8,25 @@ namespace PunchTable {
 
     struct TableParameters
     {
-        double minKE;
+        double minKE; //MeV
-        double maxKE;
+        double maxKE; //MeV
-        double stepKE;
+        double stepKE; //MeV
-        double minTheta;
+        double minTheta; //deg
-        double maxTheta;
+        double maxTheta; //deg
-        double stepTheta;
+        double stepTheta; //deg
         int projectileZ;
         int projectileA;
-        std::vector<std::vector<int>> targetZ;
+        std::vector<std::vector<int>> targetZ; //Last element is the deposition layer
         std::vector<std::vector<int>> targetA;
-        std::vector<std::vector<int>> targetS;
+        std::vector<std::vector<int>> targetS; //Stoichiometry
-        double targetThickness;
+        std::vector<double> targetThickness; //ug/cm2
         std::string filename;
     };
 
     void GenerateTable(const TableParameters& params);
+
+    //For testing (thetaIncident = radians, initialKE = MeV)
+    double GetEnergyDeposited(const TableParameters& params, double thetaIncident, double initialKE);
 }
 
 #endif

src/PunchTable.cpp

@@ -34,12 +34,15 @@ namespace PunchTable {
 		std::vector<double> energyIn, energyDep;
 
 		std::getline(input, junk);
-		input>>junk>>thickness;
+		while(junk != "---------------------------------")
+		{
+			std::getline(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")
@@ -87,7 +90,6 @@ namespace PunchTable {
 			return 0.0;
 		}
 
-		float thetaf_bin = (theta_incident - m_thetaMin)/m_thetaStep;
 		int theta_bin = (theta_incident - m_thetaMin)/m_thetaStep;
 
 		std::cout<<"theta bin: "<<theta_bin<<" theta_inc: "<<theta_incident<<std::endl;
@@ -97,6 +99,7 @@ namespace PunchTable {
 			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

src/PunchTable.h

@@ -14,7 +14,7 @@ namespace PunchTable {
 		PunchTable(const std::string& filename);
 		~PunchTable();
 		void ReadFile(const std::string& filename);
-		double GetInitialKineticEnergy(double theta_incident, double e_deposited);
+		double GetInitialKineticEnergy(double theta_incident, double e_deposited); //radians, MeV
 		inline bool IsValid() const { return m_validFlag; }
 
 	private:

src/main.cpp

@@ -1,5 +1,7 @@
 #include "MassLookup.h"
 #include "PunchTable.h"
+#include "GenerateTable.h"
+#include "catima/gwm_integrators.h"
 #include <string>
 #include <vector>
 #include <cmath>
@@ -16,94 +18,43 @@ int main(int argc, char** argv)
 		options = argv[1];
 	}
 
-	int Zp = 1;
+	PunchTable::TableParameters params;
-	int Ap = 1;
 
-	float keStep = 0.01;
+	params.projectileA = 1;
-	float thickness = 500.0 * 1e-4 * 2.3216 * 1e6; //1000 um thick times density of Si-> ug/cm^2
+	params.projectileZ = 1;
-	float deg2rad = M_PI/180.0;
-	float thetaStep = 0.1*deg2rad;
-	float precision = 1.0e-6;
 
-	std::vector<int> ZT(1), AT(1), ST(1);
+	params.minKE = 8.23;
-	ZT[0] = 14;
+	params.maxKE = 25.0;
-	AT[0] = 28;
+	params.stepKE = 0.01;
-	ST[0] = 1;
 
-	EnergyLoss energyLoss;
+	params.minTheta = 0.0;
-	energyLoss.SetTargetComponents(ZT, AT, ST);
+	params.maxTheta = 75.0;
+	params.stepTheta = 0.1;
 
-	float keMin = 8.23f;
+	double thicknessSABRE = 500.0 * 1.0e-4 * 2.3216 * 1.0e6; //500 um thick times density of Si-> ug/cm^2
-	float keMax = 25.0f;
 
-	float thetaMin = 0.0f*deg2rad;
+	params.targetZ = {{14}};
-	float thetaMax = 85.0f*deg2rad;
+	params.targetA = {{28}};
+	params.targetS = {{1}};
+	params.targetThickness = {thicknessSABRE};
 
-	int nebins = int((keMax-keMin)/keStep);
+	params.filename = "tables/test.ptab";
-	int ntbins = int ((thetaMax-thetaMin)/thetaStep);
 
 	if(options == "--make-table" || options == "")
 	{
-
+		PunchTable::GenerateTable(params);
-		std::cout<<"nebins: "<<nebins<<" ntbins: "<<ntbins<<std::endl;
-	
-		float keCur;
-		float thetaCur;
-		float energy;
-	
-		std::vector<std::vector<float>> energy_dep(ntbins);
-		std::vector<std::vector<float>> energy_in(ntbins);
-	
-		for(int i=0; i<ntbins; i++)
-		{
-			std::cout<<"\rWorking on theta bin "<<i+1<<" of "<<ntbins<<std::flush;
-			thetaCur = thetaMin + i*thetaStep;
-			for(int j=0; j<nebins; j++)
-			{
-				keCur = keMax - j*keStep;
-				energy = energyLoss.GetEnergyLoss(Zp, Ap, keCur, thickness / std::fabs(std::cos(thetaCur)));
-				if(std::fabs(energy-keCur) > precision)
-				{
-					energy_dep[i].push_back(energy);
-					energy_in[i].push_back(keCur);
-				}
-			}
-		}
-		std::cout<<std::endl;
-	
-		std::ofstream output("tables/test.ptab");
-		output<<std::setprecision(5);
-		output<<"Material: ";
-		for(size_t i=0; i<ZT.size(); i++)
-			output<<MassLookup::GetInstance()->FindSymbol(ZT[i], AT[i]);
-		output<<std::endl;
-		output<<"Thickness: "<<thickness<<std::endl;
-		output<<"IncidentAngleRange(deg): "<<thetaMin/deg2rad<<" to "<<thetaMax/deg2rad<<" stepSize: "<<thetaStep/deg2rad<<std::endl;
-		output<<std::setw(16)<<"Energy Deposited"<<"|"<<std::setw(16)<<"Intial Energy"<<std::endl;
-		output<<"---------------------------------"<<std::endl;
-		for(int i=0; i<ntbins; i++)
-		{
-			thetaCur = thetaMin + i*thetaStep;
-			output<<"begin_theta "<<thetaCur/deg2rad<<std::endl;
-			for(unsigned int j=0; j<energy_dep[i].size(); j++) 
-			{
-				output<<std::setw(16)<<energy_dep[i][j]<<" "<<std::setw(16)<<energy_in[i][j]<<std::endl;
-			}
-			output<<"end_theta"<<std::endl;
-		}
-		output.close();
 	}
 
 	if(options == "--test" || options == "")
 	{
 		std::cout<<"-------------Testing---------"<<std::endl;
-		PunchTable table("tables/test.ptab");
+		PunchTable::PunchTable table(params.filename);
 		std::cout<<"Is the table valid? "<<table.IsValid()<<std::endl;
 		double ke_test = 14.5; //MeV
-		double theta_test = 35.0*deg2rad;
+		double theta_test = 35.5*M_PI/180.0;
-		double ke_dep = energyLoss.GetEnergyLoss(Zp, Ap, ke_test, thickness/std::fabs(std::cos(theta_test)));
+		double ke_dep = PunchTable::GetEnergyDeposited(params, theta_test, ke_test);
 		double recov_ke = table.GetInitialKineticEnergy(theta_test, ke_dep);
-		std::cout<<"For a "<<MassLookup::GetInstance()->FindSymbol(Zp, Ap)<<" with kinetic energy "<<ke_test<<" MeV "
+		std::cout<<"For a "<<PunchTable::MassLookup::GetInstance().FindSymbol(params.projectileZ, params.projectileA)<<" with kinetic energy "<<ke_test<<" MeV "
 				 <<" the percent error on recovery is "<<(ke_test-recov_ke)/ke_test*100.0<<" where the recovered energy is "<<recov_ke<<" MeV"
 				 <<" and the deposited energy is "<<ke_dep<<" MeV"<<std::endl;
 		std::cout<<"-----------------------------"<<std::endl;