diff --git a/benchmarks/rich/include/CherenkovPIDAnalysis.h b/benchmarks/rich/include/CherenkovPIDAnalysis.h
index 76a310e73d028856d5068c5984059deea75863dc..857d693bd4bccf5e595dcf7d3b0160d34b508ca9 100644
--- a/benchmarks/rich/include/CherenkovPIDAnalysis.h
+++ b/benchmarks/rich/include/CherenkovPIDAnalysis.h
@@ -52,16 +52,22 @@ namespace benchmarks {
       TH1D *m_mcWavelength_dist;
       TH1D *m_mcRindex_dist;
       TH1D *m_highestWeight_dist;
+      TH1D *m_speTheta_dist;
+      TH1D *m_speResid_dist;
       TH2D *m_photonTheta_vs_photonPhi;
       // - momentum scans
       TH2D *m_npe_vs_p;
       TH2D *m_theta_vs_p;
+      TH2D *m_spetheta_vs_p;
       TH2D *m_thetaResid_vs_p;
+      TH2D *m_spethetaResid_vs_p;
       TH2D *m_highestWeight_vs_p;
       // - pseudorapidity scans
       TH2D *m_npe_vs_eta;
       TH2D *m_theta_vs_eta;
+      TH2D *m_spetheta_vs_eta;
       TH2D *m_thetaResid_vs_eta;
+      TH2D *m_spehetaResid_vs_eta;
       TH2D *m_highestWeight_vs_eta;
 
       // logger
diff --git a/benchmarks/rich/src/CherenkovPIDAnalysis.cc b/benchmarks/rich/src/CherenkovPIDAnalysis.cc
index 795dbe814d9a7640b253fddec5f7570063a4ef47..d57d7c7e796115d0e499f0ca3b14848dc10307bb 100644
--- a/benchmarks/rich/src/CherenkovPIDAnalysis.cc
+++ b/benchmarks/rich/src/CherenkovPIDAnalysis.cc
@@ -31,20 +31,32 @@ namespace benchmarks {
         );
     m_theta_dist = new TH1D(
         "theta_dist_"+m_rad_name_trun,
-        "Estimated Cherenkov Angle for "+m_rad_title+";#theta [mrad]",
+        "Reconstructed Cherenkov Angle for "+m_rad_title+";#theta [mrad]",
         Tools::theta_bins, 0, Tools::theta_max
         );
     m_thetaResid_dist = new TH1D(
         "thetaResid_dist_"+m_rad_name_trun,
-        "Estimated Cherenkov Angle Residual for "+m_rad_title+";#Delta#theta [mrad]",
+        "Reconstructed Cherenkov Angle Residual for "+m_rad_title+";#Delta#theta [mrad]",
         Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
         );
     m_photonTheta_vs_photonPhi = new TH2D(
         "photonTheta_vs_photonPhi_"+m_rad_name_trun,
-        "Estimated Photon #theta vs #phi for "+m_rad_title+";#phi [rad];#theta [mrad]",
+        "Reconstructed Photon #theta vs #phi for "+m_rad_title+";#phi [rad];#theta [mrad]",
         Tools::phi_bins, -TMath::Pi(), TMath::Pi(),
         Tools::theta_bins, 0, Tools::theta_max
         );
+    m_speTheta_dist = new TH1D(
+      "speTheta_dist_"+m_rad_name_trun,
+      "Reconstructed Photon SPE #theta for" +m_rad_title+";#theta[mrad],
+      Tools::theta_bins, 0, Tools::theta_max
+    );
+     
+    m_speResid_dist = new TH1D(
+      "speResid_dist_"+m_rad_name_trun,
+      "Reconstructed SPE Cherenkov Angle Residual for "+m_rad_title+";#Delta#theta [mrad]",
+      Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
+    );
+
 
     // truth distributions
     m_mcWavelength_dist = new TH1D(
@@ -74,13 +86,25 @@ namespace benchmarks {
         );
     m_theta_vs_p = new TH2D(
         "theta_vs_p_"+m_rad_name_trun,
-        "Estimated Cherenkov Angle vs. Particle Momentum for "+m_rad_title+";p [GeV];#theta [mrad]",
+        "Reconstructed Cherenkov Angle vs. Particle Momentum for "+m_rad_title+";p [GeV];#theta [mrad]",
         Tools::momentum_bins, 0, Tools::momentum_max,
         Tools::theta_bins, 0, Tools::theta_max
         );
     m_thetaResid_vs_p = new TH2D(
         "thetaResid_vs_p_"+m_rad_name_trun,
-        "Estimated Cherenkov Angle Residual vs. Particle Momentum for "+m_rad_title+";p [GeV];#Delta#theta [mrad]",
+        "Reconstructed Cherenkov Angle Residual vs. Particle Momentum for "+m_rad_title+";p [GeV];#Delta#theta [mrad]",
+        Tools::momentum_bins, 0, Tools::momentum_max,
+        Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
+        );
+    m_spetheta_vs_p = new TH2D(
+        "theta_vs_p_"+m_rad_name_trun,
+        "Reconstructed Cherenkov Angle vs. Particle Momentum for "+m_rad_title+";p [GeV];#theta [mrad]",
+        Tools::momentum_bins, 0, Tools::momentum_max,
+        Tools::theta_bins, 0, Tools::theta_max
+        );
+    m_spethetaResid_vs_p = new TH2D(
+        "thetaResid_vs_p_"+m_rad_name_trun,
+        "Reconstructed Cherenkov Angle Residual vs. Particle Momentum for "+m_rad_title+";p [GeV];#Delta#theta [mrad]",
         Tools::momentum_bins, 0, Tools::momentum_max,
         Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
         );
@@ -100,16 +124,29 @@ namespace benchmarks {
         );
     m_theta_vs_eta = new TH2D(
         "theta_vs_eta_"+m_rad_name_trun,
-        "Estimated Cherenkov Angle vs. Pseudorapidity for "+m_rad_title+";#eta;#theta [mrad]",
+        "Reconstructed Cherenkov Angle vs. Pseudorapidity for "+m_rad_title+";#eta;#theta [mrad]",
         Tools::eta_bins, Tools::eta_min, Tools::eta_max,
         Tools::theta_bins, 0, Tools::theta_max
         );
     m_thetaResid_vs_eta = new TH2D(
         "thetaResid_vs_eta_"+m_rad_name_trun,
-        "Estimated Cherenkov Angle Residual vs. Pseudorapidity for "+m_rad_title+";#eta;#Delta#theta [mrad]",
+        "Reconstructed Cherenkov Angle Residual vs. Pseudorapidity for "+m_rad_title+";#eta;#Delta#theta [mrad]",
         Tools::eta_bins, Tools::eta_min, Tools::eta_max,
         Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
         );
+    m_spetheta_vs_eta = new TH2D(
+        "theta_vs_eta_"+m_rad_name_trun,
+        "Reconstructed SPE Cherenkov Angle vs. Pseudorapidity for "+m_rad_title+";#eta;#theta [mrad]",
+        Tools::eta_bins, Tools::eta_min, Tools::eta_max,
+        Tools::theta_bins, 0, Tools::theta_max
+        );
+    m_spethetaResid_vs_eta = new TH2D(
+        "thetaResid_vs_eta_"+m_rad_name_trun,
+        "Reconstructed SPE Cherenkov Angle Residual vs. Pseudorapidity for "+m_rad_title+";#eta;#Delta#theta [mrad]",
+        Tools::eta_bins, Tools::eta_min, Tools::eta_max,
+        Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
+        );
+
     m_highestWeight_vs_eta = new TH2D(
         "highestWeight_vs_eta_"+m_rad_name_trun,
         "Highest PDG Weight vs. Pseudorapidity for "+m_rad_title+";#eta",
@@ -200,8 +237,15 @@ namespace benchmarks {
       m_thetaResid_dist->Fill(   theta_resid_mrad);
       m_thetaResid_vs_p->Fill(   charged_particle_momentum, theta_resid_mrad);
       m_thetaResid_vs_eta->Fill( charged_particle_eta,      theta_resid_mrad);
-      for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons())
-        m_photonTheta_vs_photonPhi->Fill( phi, theta*1e3); // [rad] -> [mrad]
+      for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons()){
+        m_photonTheta_vs_photonPhi->Fill( phi, theta*1e3);
+        m_speTheta_dist->Fill(theta*1e3);
+        m_speResid_dist->Fill((theta-theta_exp)*1e3);
+        m_spetheta_vs_p->Fill(charged_particle_momentum,theta*1e3);
+        m_spethetaResid_vs_p->Fill(charged_particle_momentum,(theta-theta_exp)*1e3);
+        m_spetheta_vs_eta->Fill(charged_particle_eta,theta*1e3);
+        m_spethetaResid_vs_eta->Fill(charged_particle_eta,(theta-theta_exp)*1e3);
+      } // [rad] -> [mrad]
       m_mcWavelength_dist->Fill(  Tools::HC / cherenkov_pid.getPhotonEnergy() ); // energy [GeV] -> wavelength [nm]
       m_mcRindex_dist->Fill( mc_rindex);