straggling

Armory/ClassAnasen.h

@@ -59,7 +59,7 @@ private:
 
   const float qqqR1 = 50;
   const float qqqR2 = 100;
-  const float qqqZPos = 23 + 75 + 30; 
+  const float qqqZPos = 100; 
 
   void CalGeometry();
 

Armory/ClassQQQ.h

@@ -49,7 +49,7 @@ private:
   const int numDet = 4;
   const float qqqR1 = 50;
   const float qqqR2 = 100;
-  const float qqqZPos = 23 + 75 + 30;
+  const float qqqZPos = 100.0;
 
   short id; // -1 when no hit
   short chUp;
@@ -73,13 +73,19 @@ private:
 
     //see https://nukephysik101.wordpress.com/2023/12/30/intersect-between-2-line-segments/
     //zero all z-component
+
     TVector3 a0 = p1; a0.SetZ(0);
     TVector3 a1 = p2; a1.SetZ(0);
 
     TVector3 b0 = q1; b0.SetZ(0);
     TVector3 b1 = q2; b1.SetZ(0);
 
-    double h = 0, k = 0; // placeholder values, implementation of intersection logic
+    double A = ((b0-b1).Cross(a0-a1)).Mag();
+
+    double h = ((b0-a0).Cross(b1-a0)).Z()/ A;
+    double k = ((a1-b0).Cross(a0-b0)).Z()/ A;
+
+
     if( verbose ) printf(" ----h, k : %f, %f\n", h, k);
 
     return std::pair<double,double>(h,k);
@@ -151,38 +157,38 @@ inline void QQQ::FindQQQPos(TVector3 pos,
 
   for(int det = 0; det < 4; det++){
 
-    double phiMin = det*90.0 + 5.0;
+    double phiMin = det*90.0;
-    double phiMax = phiMin + 85.0;
+    double phiMax = phiMin + 87.0;
 
     if( phi >= phiMin && phi <= phiMax ){
       id = det;
       break;
     }
   }
-
   if( id < 0 ) return;
+  //std::cout << id << std::endl;
 
   const double ringWidth =
-      (qqqR2 - qqqR1)/32.0;
+      (qqqR2 - qqqR1)/16.0;
 
   int ring =
       (int)((r - qqqR1)/ringWidth);
 
   if( ring < 0 ) ring = 0;
-  if( ring > 31 ) ring = 31;
+  if( ring > 15 ) ring = 15;
 
   //--------------------------------------------
   // Sector number (4 strips)
   //--------------------------------------------
 
   double localPhi =
-      phi - (id*90.0 + 5.0);
+      phi - (id*90.0);
 
   int sector =
-      (int)(localPhi/(85.0/4.0));
+      (int)(localPhi/(87.0/16.0));
 
   if( sector < 0 ) sector = 0;
-  if( sector > 3 ) sector = 3;
+  if( sector > 15 ) sector = 15;
 
   chBk = ring;
   chDn = sector;
@@ -246,20 +252,20 @@ inline void QQQ::CalQQQPos(unsigned short ID,
   hitPos.Clear();
 
   if( ID > 3 ) return;
-  if( chBack > 31 ) return;
+  if( chBack > 15 ) return;
-  if( chDown > 3 ) return;
+  if( chDown > 15 ) return;
 
   const double ringWidth =
-      (qqqR2 - qqqR1)/32.0;
+      (qqqR2 - qqqR1)/16.0;
 
   double r =
       qqqR1 + (chBack + 0.5)*ringWidth;
 
   const double sectorWidth =
-      85.0/4.0;
+      87.0/16;
 
   double phiDeg =
-      ID*90.0 + 5.0 +
+      ID*90.0 +
       (chDown + 0.5)*sectorWidth;
 
   double phi =
@@ -274,7 +280,7 @@ inline void QQQ::CalQQQPos(unsigned short ID,
   id   = ID;
   chBk = chBack;
   chDn = chDown;
-  //chUp = chUp;
+  chUp = chUp;
 }
 
 #endif 

Armory/Eloss_and_straggle.py

@@ -0,0 +1,120 @@
+import pycatima as catima
+import numpy as np, sys
+import matplotlib.pyplot as plt
+
+# --- 1. Constants ---
+P_TORR = 250
+P_CO2 = 3
+if(len(sys.argv)==3):
+    P_TORR=int(sys.argv[1])
+    P_CO2 = int(sys.argv[2])
+
+TEMP_K = 293.15
+R = 8.3144
+MEV2U = 1.0 / 931.494
+
+# Gas Density Calculations
+p_pa = P_TORR * 133.322
+molar_density = p_pa / (R * TEMP_K)
+m_he, m_c, m_o= 4.0026, 12.0000, 15.9949
+m_mix_avg = ((1 - P_CO2 / 100) * m_he) + (P_CO2 / 100 * (m_c + 2*m_o))
+rho_g_cm3 = (molar_density * m_mix_avg) / 1e6
+print(f"Gas density at {P_TORR} Torr: {rho_g_cm3:.6e} g/cm^3")
+
+# --- 2. Material & Step Setup ---
+material_def = [(m_he, 2, (1 - P_CO2 / 100)), (m_c, 6, P_CO2 / 100), (m_o, 8, 2*P_CO2 / 100)]
+gas_mix = catima.Material(material_def)
+gas_mix.density(rho_g_cm3)
+
+# Thickness step settings
+step_mg_cm2 = 0.001           # 1 ug/cm2 steps as per your example -- kept fine for
+                               # numerical accuracy of the dedx integration itself.
+step_g_cm2 = step_mg_cm2 / 1000.0
+max_steps = 1000000000             # Adjust based on how far you want to track
+
+coarse_step_cm = 0.2   # row spacing over most of the track
+fine_step_cm = 0.03     # row spacing near the Bragg peak
+fine_zone_frac = 0.085   # fraction of the *total* range treated as "near the peak"
+
+# Set relative integration tolerance (lower = higher precision, slower calculation)
+catima.Config.epsrel = 1e-6
+# Set absolute integration tolerance
+catima.Config.epsabs = 1e-9
+
+def generate_lookup(z, mass_u, e_start_mev, label):
+    filename = f"Eloss/E_vs_x_{label}.dat"
+    header = f"Energy(MeV) \tmg/cm2 \tcm\tEin\tsigE\tsigA\tsigR\tsigX\tcov\ttof\tsp\nStarting Energy: {e_start_mev} MeV"
+
+    projectile = catima.Projectile(mass_u, z)
+    e_u_init = e_start_mev / mass_u
+
+    # 1. Get exact analytical range directly from CATIMA
+    projectile.T(e_u_init)
+    total_range_g_cm2 = catima.range(projectile, gas_mix)
+    total_range_cm = total_range_g_cm2 / rho_g_cm3
+    fine_zone_start_cm = total_range_cm * (1.0 - fine_zone_frac)
+
+    output = []
+    current_dist_cm = 0.0
+
+    # 2. Step directly at checkpoint resolution
+    while current_dist_cm <= total_range_cm * 1.05: # Track slightly past range
+        thickness_g_cm2 = current_dist_cm * rho_g_cm3
+
+        # Configure target thickness for this checkpoint
+        gas_mix.density(rho_g_cm3).thickness(thickness_g_cm2)
+        projectile.T(e_u_init)
+        result = catima.calculate(projectile, gas_mix)
+
+        # Save checkpoint row
+        e_total_out = result.Eout * mass_u
+        output.append([
+            e_total_out, thickness_g_cm2 * 1000.0, current_dist_cm,
+            result.Ein, result.sigma_E, result.sigma_a, result.sigma_r,
+            result.sigma_x, result.cov, result.tof, result.sp
+        ])
+
+        if result.Eout == 0:
+            break
+
+        # Adaptive spatial step
+        step = fine_step_cm if current_dist_cm >= fine_zone_start_cm else coarse_step_cm
+        current_dist_cm += step
+
+    np.savetxt(filename, output, fmt='%.6f', delimiter='\t', header=header)
+    print(f"Lookup table created: {filename} ({len(output)} rows)")
+
+    data = np.array(output)
+    energy_mev, dist_cm, sigma_e, tof = data[:, 0], data[:, 2], data[:, 4], data[:, 9]
+
+    fig, (ax_e, ax_sig, ax_tof) = plt.subplots(3, 1, figsize=(8, 12), sharex=True)
+
+    ax_e.plot(dist_cm, energy_mev)
+    ax_e.fill_between(dist_cm, energy_mev - sigma_e, energy_mev + sigma_e, alpha=0.3)
+    ax_e.set_ylabel("Energy (MeV)")
+    ax_e.set_title(f"Energy Loss Curve {label.capitalize()}")
+    ax_e.grid(True)
+
+    ax_sig.plot(dist_cm, sigma_e)
+    ax_sig.set_ylabel("Energy straggle $\\sigma_E$ (MeV)")
+    ax_sig.grid(True)
+
+    ax_tof.plot(dist_cm, tof)
+    ax_tof.set_xlabel("Distance (cm)")
+    ax_tof.set_ylabel("Time of flight (ns)")
+    ax_tof.grid(True)
+
+    plt.tight_layout()
+    plt.show()
+    
+
+
+# --- 3. Run ---
+# Format: generate_lookup(Z, mass_u, E_start_MeV, label)
+generate_lookup(1, 1.0078, 30, "proton")
+generate_lookup(1, 2.01355, 30, "deuteron")
+generate_lookup(2, 4.0026, 50, "alpha")
+generate_lookup(13,26.9815, 80, "aluminum")
+#generate_lookup(9,17.0021, 70, "fluorine")
+#generate_lookup(8,15.9949, 70, "oxygen")
+

Armory/Makefile

@@ -69,6 +69,6 @@ EventBuilder: EventBuilder.cpp ClassData.h fsuReader.h Hit.h
 	@echo "--------- making EventBuilder"
 	$(CXX) $(CXXFLAGS) EventBuilder.cpp -o EventBuilder $(LDFLAGS)
 
-AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h EnergyLoss.h AutoHist2D.h
+AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassQQQ.h ClassPW.h ClassAnasen.h EnergyLoss.h AutoHist2D.h anasen_anode_cathode_hyperboloids.h
 	@echo "--------- making ANASEN Monte Carlo"
 	$(CXX) $(CXXFLAGS) anasenMS.cpp -o AnasenMS $(LDFLAGS)

Armory/anasenMS.cpp

@@ -8,6 +8,7 @@
 #include "TBenchmark.h"   // timing measurement
 #include "TGraph.h"       // for energy loss interpolation
 #include <cstring>
+#include <algorithm>
 #include "TApplication.h" // ROOT app loop
 #include "ClassTransfer.h" // Reaction kinematics and MC event generation
 #include "ClassAnasen.h"   // ANASEN detector model classes (SX3, PW, etc.)
@@ -42,6 +43,13 @@ TGraph* LoadELoss(const std::string& filename) {
   return g;
 }
 
+// Loads column 2 (Energy_MeV) vs column 4 (Sigma_x_cm) from an E_vs_x_*.dat table
+TGraph* LoadSigmaXVsEnergy(const std::string& filename) {
+  TGraph* g = new TGraph(filename.c_str(), "%*lg %lg %*lg %lg");
+  g->Sort(); // TGraph::Eval requires ascending x (Energy_MeV)
+  return g;
+}
+
 bool IsDeadAnode(int id){
     static std::set<int> dead = {}; // add dead anode IDs here, 0-23
     return dead.count(id);
@@ -53,7 +61,7 @@ bool IsDeadCathode(int id){
 }
 
 bool IsDeadSX3(int id){ 
-    static std::set<int> dead = {0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3
+    static std::set<int> dead = {};//{0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3
     return dead.count(id);
 }
 
@@ -73,7 +81,7 @@ bool IsDeadSX3FrontDnChannel(int sx3ID, int chDn){
 
 bool IsDeadSX3BackChannel(int sx3ID, int chBk){
   static std::set<std::pair<int, int>> dead = {
-    {1, 10}
+    //{1, 10}
     // {sx3ID, back-channel} 
   };
   return dead.count({sx3ID, chBk});
@@ -186,6 +194,7 @@ int main(int argc, char **argv){
   const double beamEntranceZ = -280 - 174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
   
   TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = path length (cm), y = beam energy (MeV)
+  TGraph* sigmaXBeam = LoadSigmaXVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_x (cm)
   // Build a temporary inverse (energy -> path) to locate the path at beamE.
   TGraph* elossBeamInverseRaw = new TGraph(elossBeam->GetN());
   for( int p = 0; p < elossBeam->GetN(); p++ ){
@@ -383,15 +392,24 @@ int main(int argc, char **argv){
   tree1->Branch("cDist",  cathodeDist, "cathodeDist/D");
 
   // SX3 channel assignment and Z fraction (depth) information
-  int sx3ID, sx3Up, sx3Dn, sx3Bk, qqqID;
+  int sx3ID, sx3Up, sx3Dn, sx3Bk, qqqID, qqqUp, qqqBk;
   double sx3ZFrac;
   tree1->Branch("sx3ID",   &sx3ID, "sx3ID/I");
-  tree1->Branch("qqqID",   &qqqID, "qqqID/I");
   tree1->Branch("sx3Up",   &sx3Up, "sx3Up/I");
   tree1->Branch("sx3Dn",   &sx3Dn, "sx3Dn/I");
   tree1->Branch("sx3Bk",   &sx3Bk, "sx3Bk/I");
   tree1->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
 
+  tree1->Branch("qqqID",   &qqqID, "qqqID/I");
+  tree1->Branch("qqqUp",   &qqqUp, "qqqUp/I");
+  tree1->Branch("qqqBk",   &qqqBk, "qqqBk/I");
+
+  double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN//, Ex_recon=NAN
+    ;
+  tree1->Branch("EBeam_Kin", &EBeam_Kin_gs, "EBeam_Kin/D");
+  tree1->Branch("EBeam_Kin_2.2", &EBeam_Kin_2_2, "EBeam_Kin_2.2/D");
+  tree1->Branch("EBeam_Kin_3.4", &EBeam_Kin_3_4, "EBeam_Kin_3.4/D");
+
   // reconstructed angles from PW track fit, method 1 and 2
   double reTheta, rePhi;
   tree1->Branch("reTheta", &reTheta, "reconstucted_theta/D");
@@ -464,13 +482,18 @@ int main(int argc, char **argv){
     transfer.CalReactionConstant();
 
     // vertex position in target volume
-    vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() +  vertexXRange[0];
-    vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() +  vertexYRange[0];
     beamEnergy = gRandom->Uniform(0, beamE); // MeV, sample beam energy at vertex from uniform distribution between 0 and initial beam energy
     KEA = beamEnergy / beamA;
     beamPath_cm = elossBeamInverse->Eval(beamEnergy); // beamE maps to x=0 after path shift
     vertexZ = beamEntranceZ + beamPath_cm * 10.0; // cm -> mm
 
+    // transverse sampling range from the beam's distance straggle at this energy
+    const double sigmaX_mm = std::max(0.0, sigmaXBeam->Eval(beamEnergy)) * 10.0; // cm -> mm
+    vertexX = 2.0 * sigmaX_mm * gRandom->Rndm() - sigmaX_mm;
+    vertexY = 2.0 * sigmaX_mm * gRandom->Rndm() - sigmaX_mm;
+
+    //vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() +  vertexXRange[0];
+    //vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() +  vertexYRange[0];
     //vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() +  vertexZRange[0];
 
     TVector3 vertex(vertexX, vertexY, vertexZ);
@@ -514,6 +537,8 @@ int main(int argc, char **argv){
     dir.SetTheta(thetab * TMath::DegToRad());
     dir.SetPhi(phib * TMath::DegToRad());
 
+    qqq->Clear();
+    sx3->Clear();
     // run detector response models for PW and SX3
     pw->FindWireID(vertex, dir, false);
     sx3->FindSX3Pos(vertex, dir, false);
@@ -529,14 +554,14 @@ int main(int argc, char **argv){
     anodeDist[1]   = hitInfo.nextNearestDist.first; // distance to next nearest anode wire
     cathodeDist[1] = hitInfo.nextNearestDist.second; // distance to next nearest cathode wire
 
-    //if(IsDeadAnode(anodeID[0])) continue;
+    if(IsDeadAnode(anodeID[0])) anodeID[0] = -1; // mark as no hit if anode is dead
-    //if(IsDeadCathode(cathodeID[0])) continue;
+    if(IsDeadCathode(cathodeID[0])) cathodeID[0] = -1; // mark as no hit if cathode is dead
 
     // SX3 hit channel info and depth fraction
     sx3ID = sx3->GetID();
     qqqID = qqq->GetID();
 
-    //if(IsDeadSX3(sx3ID)) continue;
+    if(IsDeadSX3(sx3ID)) sx3ID = -1; // mark as no hit if SX3 is dead
 
     anodeDist[0]   = hitInfo.nearestDist.first; // distance to nearest anode wire
     cathodeDist[0] = hitInfo.nearestDist.second; // distance to nearest cathode wire
@@ -547,7 +572,7 @@ int main(int argc, char **argv){
       sx3Up    = sx3->GetChUp();
       sx3Dn    = sx3->GetChDn();
       sx3Bk    = sx3->GetChBk();
-      //if(IsDeadSX3ChannelCombo(sx3ID, sx3Up, sx3Dn, sx3Bk)) continue;
+      if(IsDeadSX3ChannelCombo(sx3ID, sx3Up, sx3Dn, sx3Bk)) sx3Up = -1, sx3Dn = -1, sx3Bk = -1; // mark as no hit if any SX3 channel is dead
       sx3ZFrac = sx3->GetZFrac();
       
       // apply intrinsic detector resolution to true SX3 hit position
@@ -638,6 +663,15 @@ int main(int argc, char **argv){
       AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)");
       AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(thetab * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)");
       //tree1->Fill();
+      /*
+      //Kinematics aakin_27Al(26.981538408,4.00260325413,4.0026035413,26.981538408,beam_energy_at_vertex/26.981538408); //m3 is alpha
+      Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,beamEnergy/26.981538408); //m3 is proton
+      //Kinematics apkin_27Al(1.00782503224,4.00260325413,4.00260325413,1.00782503224,beamEnergy/1.00782503224); //m3 is proton
+      Ex_recon = apkin_27Al.getExc(Tb, thetab);
+      EBeam_Kin_gs = apkin_27Al.getEbeam_givenQ(Tb, 0.0, thetab);
+      EBeam_Kin_2_2 = apkin_27Al.getEbeam_givenQ(Tb, 2.2, thetab);
+      EBeam_Kin_3_4 = apkin_27Al.getEbeam_givenQ(Tb, 3.4, thetab);
+      //std::cout << EBeam_Kin << std::endl;*/
 
     }else if (qqqID >= 0){
 

Armory/anasen_anode_cathode_hyperboloids.h

@@ -30,8 +30,12 @@ inline std::tuple<TVector3,TVector3,double> find_PC_PathLength(const TVector3& x
 	// Cathode a, c values are found by scaling up the anode waist by 43/37, the ratio of the outermost radii
 	TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895);
 	TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895);
+	TVector3 gw_intersect = onesheet_hyperboloid_intersect(27.712,301.895);
 	if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321)
-		return std::tuple(cathode_intersect,anode_intersect,(cathode_intersect-anode_intersect).Mag()*0.1);
+		return std::tuple(cathode_intersect,gw_intersect,(cathode_intersect-gw_intersect).Mag()*0.1);
 	else
 		return std::tuple(TVector3(0,0,0), TVector3(0,0,0), 54321);
 }
+
+
+

ELoss/EnergyLoss.py

@@ -14,7 +14,6 @@ from scipy.interpolate import interp1d
 import uproot
 import pycatima as catima
 from scipy.integrate import cumulative_trapezoid
-#matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 import cmd
 import shlex
@@ -168,9 +167,22 @@ def make_E_vs_x(
     x = x[::-1]
     x = -x
 
+    # Energy and range straggling at each cumulative depth, starting from emax_mev
+    e_u_init = emax_mev / mass_u
+    sigma_E = np.zeros_like(E)
+    sigma_x = np.zeros_like(E)
+    for i, depth_cm in enumerate(x):
+        gas.density(rho_g_cm3).thickness(depth_cm * rho_g_cm3)
+        projectile.T(e_u_init)
+        result = catima.calculate(projectile, gas)
+        sigma_E[i] = result.sigma_E
+        sigma_x[i] = result.sigma_x
+
     df = pd.DataFrame({
         "Distance_cm": x,
-        "Energy_MeV": E
+        "Energy_MeV": E,
+        "Sigma_E_MeV": sigma_E,
+        "Sigma_x_cm": sigma_x
     })
 
     outfile = get_loss_table_path(medium, label)
@@ -188,10 +200,10 @@ def make_E_vs_x(
 def load_table(filename):
     """
     Load table with columns:
-        x(cm)   E(MeV)   Sigma_E(MeV) [optional]
+        x(cm)   E(MeV)   Sigma_E(MeV) [optional]   Sigma_x(cm) [optional]
 
     Returns:
-        x_array, E_array, sigma_E_array (None if column absent)
+        x_array, E_array, sigma_E_array, sigma_x_array (None if column absent)
     """
 
     data = pd.read_csv(
@@ -205,8 +217,9 @@ def load_table(filename):
     x = data.iloc[:, 0].to_numpy()
     E = data.iloc[:, 1].to_numpy()
     sigma_E = data.iloc[:, 2].to_numpy() if data.shape[1] > 2 else None
+    sigma_x = data.iloc[:, 3].to_numpy() if data.shape[1] > 3 else None
 
-    return x, E, sigma_E
+    return x, E, sigma_E, sigma_x
 
 def get_interpolators(particle, medium):
 
@@ -243,7 +256,7 @@ def get_interpolators(particle, medium):
     if filename is None:
         filename = get_loss_table_path(medium, canonical_particle)
 
-    x, E, sigma_E = load_table(filename)
+    x, E, sigma_E, sigma_x = load_table(filename)
 
     E_of_x = interp1d(
         x,
@@ -269,13 +282,23 @@ def get_interpolators(particle, medium):
     else:
         sigma_of_x = None
 
-    interp_cache[cache_key] = (E_of_x, x_of_E, sigma_of_x)
+    if sigma_x is not None:
+        sigma_x_of_x = interp1d(
+            x,
+            sigma_x,
+            bounds_error=False,
+            fill_value="extrapolate"
+        )
+    else:
+        sigma_x_of_x = None
 
-    return E_of_x, x_of_E, sigma_of_x
+    interp_cache[cache_key] = (E_of_x, x_of_E, sigma_of_x, sigma_x_of_x)
+
+    return E_of_x, x_of_E, sigma_of_x, sigma_x_of_x
 
 def energy_loss(particle, medium, Ei, dl):
 
-    E_of_x, x_of_E, _ = get_interpolators(particle, medium)
+    E_of_x, x_of_E, _, _ = get_interpolators(particle, medium)
 
     xi = x_of_E(Ei)
 
@@ -291,7 +314,7 @@ def energy_loss(particle, medium, Ei, dl):
     
 def energy_reconstruction(particle, medium, Ef, dl):
 
-    E_of_x, x_of_E, _ = get_interpolators(particle, medium)
+    E_of_x, x_of_E, _, _ = get_interpolators(particle, medium)
 
     xf = x_of_E(Ef)
 
@@ -303,7 +326,7 @@ def energy_reconstruction(particle, medium, Ef, dl):
 
 def energy_distance(particle, medium, Ei, Ef):
 
-    _, x_of_E, _ = get_interpolators(particle, medium)
+    _, x_of_E, _, _ = get_interpolators(particle, medium)
 
     xi = x_of_E(Ei)