Added table generation/reading for reverse energy loss tables. Added input configuration file handling. Testing to come.

input.txt

@@ -0,0 +1,11 @@
+TableType: ElossTable
+MinKE(MeV): 0.2 MaxKE(MeV): 20.0 StepKE(MeV): 0.01
+MinTheta(deg): 0.0 MaxTheta(deg): 75.0 StepTheta(deg): 0.1
+ProjectileZ: 1 ProjectileA: 1
+begin_materials
+    begin_material
+        thickness: 116830.0
+        73 181 1
+    end_material
+end_materials
+OutputFilename: tables/ta_eloss.etab

src/CMakeLists.txt

@@ -8,6 +8,8 @@ target_sources(PunchTable PRIVATE
     MassLookup.cpp
     PunchTable.h
     PunchTable.cpp
+    ElossTable.h
+    ElossTable.cpp
     GenerateTable.h
     GenerateTable.cpp
     main.cpp

src/ElossTable.cpp

@@ -0,0 +1,114 @@
+#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;
+
+		std::getline(input, junk);
+		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")
+			{
+				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 finalEnergy, double theta_incident)
+    {
+        theta_incident /= s_deg2rad;
+		if(!m_isValid)
+		{
+			std::cerr<<"ElossTable not initialized at GetEnergyLoss()"<<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 ElossTable::GetEnergyLoss()"<<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 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/ElossTable.h

@@ -0,0 +1,38 @@
+#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);
+        double GetEnergyLoss(double finalEnergy, double theta_incident);
+
+        inline const bool IsValid() const { return m_isValid; }
+
+    private:
+        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/GenerateTable.cpp

@@ -9,7 +9,62 @@
 
 namespace PunchTable {
 
-    void GenerateTable(const TableParameters& params)
+    TableParameters GetTableParameters(const std::string& filename)
+    {
+        TableParameters result;
+        std::ifstream input(filename);
+        if(!input.is_open())
+        {
+            std::cerr<<"Unable to open table parameter file at GetTableParameters()!"<<std::endl;
+            return result;
+        }
+
+        std::string junk;
+        std::vector<int> z, a, s;
+        double thickness;
+
+        input>>junk>>result.tableType;
+        input>>junk>>result.minKE>>junk>>result.maxKE>>junk>>result.stepKE;
+        input>>junk>>result.minTheta>>junk>>result.maxTheta>>junk>>result.stepTheta;
+        input>>junk>>result.projectileZ>>junk>>result.projectileA;
+        input>>junk;
+        if(junk != "begin_materials")
+        {
+            std::cerr<<"Parsing error in GetTableParameters()! Materials not in right place"<<std::endl;
+        }
+        while(input>>junk)
+        {
+            if(junk == "end_materials")
+                break;
+            else if(junk == "begin_material")
+            {
+                input>>junk>>thickness;
+                while(input>>junk)
+                {
+                    if(junk == "end_material")
+                        break;
+                    z.push_back(std::stoi(junk));
+                    input>>junk;
+                    a.push_back(std::stoi(junk));
+                    input>>junk;
+                    s.push_back(std::stoi(junk));
+                }
+
+                result.targetZ.push_back(z);
+                result.targetA.push_back(a);
+                result.targetS.push_back(s);
+                result.targetThickness.push_back(thickness);
+            }
+            else
+            {
+                std::cerr<<"Parsing error in GetTableParameters()! Material structure error"<<std::endl;
+            }
+        }
+        input>>junk>>result.filename;
+        return result;
+    }
+
+    void GeneratePunchTable(const TableParameters& params)
     {
         static double s_deg2rad = M_PI/180.0; //deg -> radians
         static double s_energyPrecision = 0.001; // keV precision
@@ -23,7 +78,7 @@ namespace PunchTable {
            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;
+            std::cerr<<"ERR --  Invalid target parameters passed to GeneratePunchTable. Mismatching target arrays."<<std::endl;
             return;
         }
 
@@ -46,7 +101,7 @@ namespace PunchTable {
         std::ofstream output(params.filename);
         if(!output.is_open())
         {
-            std::cerr<<"ERR -- Unable to open output file "<<params.filename<<" at GenerateTable"<<std::endl;
+            std::cerr<<"ERR -- Unable to open output file "<<params.filename<<" at GeneratePunchTable"<<std::endl;
             return;
         }
 		output<<std::setprecision(5);
@@ -135,6 +190,110 @@ namespace PunchTable {
 		output.close();
     }
 
+    void GenerateElossTable(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 GenerateElossTable. 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 nlayers = materials.size();
+        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 GenerateElossTable"<<std::endl;
+            return;
+        }
+		output<<std::setprecision(5);
+
+		output<<"Materials: "<<std::endl;;
+		for(size_t i=0; i<nlayers; 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<<"---------------------------------"<<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> energyFinalData;
+        std::vector<double> energyLossData;
+
+        double theta, ke;
+        double eloss;
+
+        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;
+            energyLossData.clear();
+            energyFinalData.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;
+                eloss = 0.0;
+                for(size_t k=0; k<nlayers; k++)
+                {
+                    materials[k].thickness(params.targetThickness[k]*s_ug2g/std::fabs(std::cos(theta)));
+                    eloss += catima::reverse_integrate_energyloss(projectile, materials[k]);
+                }
+
+                energyFinalData.push_back(ke);
+                energyLossData.push_back(eloss);
+            }
+
+            output<<"begin_theta "<<theta<<std::endl;
+			for(size_t j=0; j<energyLossData.size(); j++) 
+			{
+				output<<std::setw(16)<<energyFinalData[j]<<" "<<std::setw(16)<<energyLossData[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

src/GenerateTable.h

@@ -8,6 +8,7 @@ namespace PunchTable {
 
     struct TableParameters
     {
+        std::string tableType; //Either PunchTable or ElossTable
         double minKE; //MeV
         double maxKE; //MeV
         double stepKE; //MeV
@@ -16,14 +17,18 @@ namespace PunchTable {
         double stepTheta; //deg
         int projectileZ;
         int projectileA;
-        std::vector<std::vector<int>> targetZ; //Last element is the deposition layer
+        std::vector<std::vector<int>> targetZ; //Last element is the deposition layer for PunchTable
         std::vector<std::vector<int>> targetA;
         std::vector<std::vector<int>> targetS; //Stoichiometry
         std::vector<double> targetThickness; //ug/cm2
         std::string filename;
     };
 
-    void GenerateTable(const TableParameters& params);
+    TableParameters GetTableParameters(const std::string& filename);
+
+    void GeneratePunchTable(const TableParameters& params);
+
+    void GenerateElossTable(const TableParameters& params);
 
     //For testing (thetaIncident = radians, initialKE = MeV)
     double GetEnergyDeposited(const TableParameters& params, double thetaIncident, double initialKE);

src/main.cpp

@@ -1,5 +1,6 @@
 #include "MassLookup.h"
 #include "PunchTable.h"
+#include "ElossTable.h"
 #include "GenerateTable.h"
 #include "catima/gwm_integrators.h"
 #include <string>
@@ -12,41 +13,32 @@
 int main(int argc, char** argv)
 {
 
-	std::string options = "";
-	if(argc == 2)
+	bool test = true;
+	if(argc != 2)
 	{
-		options = argv[1];
+		std::cerr<<"PunchTable requires an input file!"<<std::endl;
+		return 1;
 	}
 
-	PunchTable::TableParameters params;
+	PunchTable::TableParameters params = PunchTable::GetTableParameters(argv[1]);
 
-	params.projectileA = 1;
-	params.projectileZ = 1;
-
-	params.minKE = 8.23;
-	params.maxKE = 25.0;
-	params.stepKE = 0.01;
-
-	params.minTheta = 0.0;
-	params.maxTheta = 75.0;
-	params.stepTheta = 0.1;
-
-	double thicknessSABRE = 500.0 * 1.0e-4 * 2.3216 * 1.0e6; //500 um thick times density of Si-> ug/cm^2
-
-	params.targetZ = {{14}};
-	params.targetA = {{28}};
-	params.targetS = {{1}};
-	params.targetThickness = {thicknessSABRE};
-
-	params.filename = "tables/test.ptab";
-
-	if(options == "--make-table" || options == "")
+	if(params.tableType == "PunchTable")
 	{
-		PunchTable::GenerateTable(params);
+		PunchTable::GeneratePunchTable(params);
+	}
+	else if(params.tableType == "ElossTable")
+	{
+		PunchTable::GenerateElossTable(params);
+	}
+	else
+	{
+		std::cerr<<"Unrecognized table type "<<params.tableType<<std::endl;
+		return 1;
 	}
 
-	if(options == "--test" || options == "")
+	if(test && params.tableType == "PunchTable")
 	{
+		std::cout<<std::endl;
 		std::cout<<"-------------Testing---------"<<std::endl;
 		PunchTable::PunchTable table(params.filename);
 		std::cout<<"Is the table valid? "<<table.IsValid()<<std::endl;
@@ -59,4 +51,19 @@ int main(int argc, char** argv)
 				 <<" and the deposited energy is "<<ke_dep<<" MeV"<<std::endl;
 		std::cout<<"-----------------------------"<<std::endl;
 	}
+	if(test && params.tableType == "ElossTable")
+	{
+		std::cout<<std::endl;
+		std::cout<<"-------------Testing---------"<<std::endl;
+		PunchTable::ElossTable table(params.filename);
+		std::cout<<"Is the table valid? "<<table.IsValid()<<std::endl;
+		double ke_test = 14.5; //MeV
+		double theta_test = 35.5*M_PI/180.0;
+		double eloss = table.GetEnergyLoss(theta_test, ke_test);
+		std::cout<<"For a "<<PunchTable::MassLookup::GetInstance().FindSymbol(params.projectileZ, params.projectileA)<<" with kinetic energy "<<ke_test<<" MeV "
+				 <<"calculate an energy loss of "<<eloss<<" MeV. Please compare to favorite tool (LISE, SRIM, etc.)"<<std::endl;
+		std::cout<<"-----------------------------"<<std::endl;
+	}
+
+	return 0;
 }