diff --git a/benchmarks/nHcal/basic_distribution/scripts/basic_distribution_analysis.cxx b/benchmarks/nHcal/basic_distribution/scripts/basic_distribution_analysis.cxx
index ccdaf01f073456a433384dbe98825299e851455c..14b5820c2cba882d2d720b1254716324886bc7d1 100644
--- a/benchmarks/nHcal/basic_distribution/scripts/basic_distribution_analysis.cxx
+++ b/benchmarks/nHcal/basic_distribution/scripts/basic_distribution_analysis.cxx
@@ -91,6 +91,7 @@ int basic_distribution_analysis(const string &filename, TString outname)
     vector<double> totalEventHits;
     vector<double> zLayerValues;
     vector<int> nHitsPerLayer;
+    vector<double> nEnergyPerLayer;
     
     for (unsigned ev = 0; ev < nEvents; ev++) 
     {
@@ -109,7 +110,7 @@ int basic_distribution_analysis(const string &filename, TString outname)
             cerr << "HcalEndcapNHits collection is not valid!" << endl;
         }
 
-        map<double, int> hitsPerLayerZ;
+        map<double, pair<int, double>> layerData;
         double totalEnergy = 0;
         int totalHits = 0;
 
@@ -126,16 +127,18 @@ int basic_distribution_analysis(const string &filename, TString outname)
             rPosition.push_back(sqrt(pos.x * pos.x + pos.y * pos.y));
 
             double zBin = round(pos.z);
-            hitsPerLayerZ[zBin]++;
+            layerData[zBin].first++;         
+            layerData[zBin].second += hit.getEnergy();
         }
 
         totalEventEnergy.push_back(totalEnergy);
         totalEventHits.push_back(totalHits);
 
-        for (const auto& [zValue, nHits] : hitsPerLayerZ)
+        for (const auto& [zValue, stats] : layerData)
         {
             zLayerValues.push_back(zValue);
-            nHitsPerLayer.push_back(nHits);
+            nHitsPerLayer.push_back(stats.first);
+            nEnergyPerLayer.push_back(stats.second);
         }
 
     }
@@ -149,24 +152,29 @@ int basic_distribution_analysis(const string &filename, TString outname)
     auto minmax_energy           = minmax_element(energy.begin(), energy.end());
     auto minmax_zLayerValues     = minmax_element(zLayerValues.begin(), zLayerValues.end());
     auto minmax_nHitsPerLayer    = minmax_element(nHitsPerLayer.begin(), nHitsPerLayer.end());
+    auto minmax_nEnergyPerLayer  = minmax_element(nEnergyPerLayer.begin(), nEnergyPerLayer.end());
 
     int nbins = nEvents/1000;
 
     TH1D *h_energyTotal = new TH1D("h_energyTotal", "Total Energy per Event; Energy per Event", 
                                     nbins/2, *minmax_totalEventEnergy.first, *minmax_totalEventEnergy.second);
     TH2D *h_layerEnergy = new TH2D("h_layerEnergy", "Energy in Layers (Z); Z; Energy", 
-                                    nbins/4, *minmax_zPosition.first, *minmax_zPosition.second, 
-                                    nbins/5, *minmax_energy.first, *minmax_energy.second);
+                                    nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second, 
+                                    nbins/4, *minmax_nEnergyPerLayer.first, *minmax_nEnergyPerLayer.second);
+    TProfile *p_layerEnergy = new TProfile("p_layerEnergy", "Energy in Layers (Z); Z; Mean Energy",
+                                    nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second);
     TH1D *h_hitCount    = new TH1D("h_hitCount", "Number of Hits per Event; Hits per Event", 
                                     nbins/2, *minmax_totalEventHits.first, *minmax_totalEventHits.second);
     TH2D *h_layerHits   = new TH2D("h_layerHits", "Number of Hits in Layers (Z); Layer(Z); Hits", 
-                                    nbins/4, *minmax_zLayerValues.first, *minmax_zLayerValues.second, 
+                                    nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second, 
                                     *minmax_nHitsPerLayer.second, *minmax_nHitsPerLayer.first, *minmax_nHitsPerLayer.second);
+    TProfile *p_layerHits = new TProfile("p_layerHits", "Number of Hits in Layers (Z); Z; Mean Hits",
+                                    nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second);
     TH2D *h_XYPos       = new TH2D("h_XYPos", "Hits position X,Y;X [mm];Y [mm]", 
                                     nbins, *minmax_xPosition.first, *minmax_xPosition.second, 
                                     nbins, *minmax_yPosition.first, *minmax_yPosition.second);
     TH2D *h_ZRPos       = new TH2D("h_ZRPos", "Hits position Z,R;Z [mm];R [mm]", 
-                                    nbins/4, *minmax_zPosition.first, *minmax_zPosition.second, 
+                                    nbins/5, *minmax_zPosition.first, *minmax_zPosition.second, 
                                     nbins, *minmax_rPosition.first, *minmax_rPosition.second);
     TH2D *h_XYEnergy    = new TH2D("h_XYEnergy", "Hits energy X,Y;X [mm];Y [mm]", 
                                     nbins, *minmax_xPosition.first, *minmax_xPosition.second, 
@@ -174,7 +182,7 @@ int basic_distribution_analysis(const string &filename, TString outname)
 
     for(int i = 0; i < xPosition.size(); i++)
     {
-        h_layerEnergy->Fill(zPosition[i], energy[i]);
+        
         h_XYPos->Fill(xPosition[i], yPosition[i]);
         h_ZRPos->Fill(zPosition[i], rPosition[i]);
         h_XYEnergy->Fill(xPosition[i], yPosition[i], energy[i]);
@@ -183,6 +191,10 @@ int basic_distribution_analysis(const string &filename, TString outname)
     for (int i = 0; i < zLayerValues.size(); i++)
     {
         h_layerHits->Fill(zLayerValues[i], nHitsPerLayer[i]);
+        p_layerHits->Fill(zLayerValues[i], nHitsPerLayer[i]);
+
+        h_layerEnergy->Fill(zLayerValues[i], nEnergyPerLayer[i]);
+        p_layerEnergy->Fill(zLayerValues[i], nEnergyPerLayer[i]);
     }
 
     for(int i = 0; i < nEvents; i++)
@@ -193,16 +205,18 @@ int basic_distribution_analysis(const string &filename, TString outname)
 
     gStyle->SetPadLeftMargin(0.13);
 
-    TCanvas *canvas = new TCanvas("canvas", "canvas", 1600, 600);
+    TCanvas *canvas = new TCanvas("canvas", "canvas", 1600, 800);
     canvas->Divide(4,2);
     canvas->cd(1);
     h_energyTotal->Draw();
     canvas->cd(2);
     h_layerEnergy->Draw("COLZ");
+    p_layerEnergy->Draw("SAME");
     canvas->cd(3);
     h_hitCount->Draw();
     canvas->cd(4);
     h_layerHits->Draw("COLZ");
+    p_layerHits->Draw("SAME");
     canvas->cd(5);
     h_XYPos->Draw("COLZ");
     canvas->cd(6);