diff --git a/scripts/debug/hda_pulse_records_plotter.py b/scripts/debug/hda_pulse_records_plotter.py
index bae72139b..524871fd4 100755
--- a/scripts/debug/hda_pulse_records_plotter.py
+++ b/scripts/debug/hda_pulse_records_plotter.py
@@ -76,6 +76,9 @@ def read_records(path, sep=','):
     intensity_calc = read_record(os.path.join(
         path, 'intensity_calc.csv'
     ), sep)
+    subray_sim = read_record(os.path.join(
+        path, 'subray_sim.csv'
+    ), sep)
     # Return key-word records
     return {
         # Intensity calculation records
@@ -85,7 +88,13 @@ def read_records(path, sep=','):
         'radius_m': intensity_calc[:, 6],
         'bdrf': intensity_calc[:, 7],
         'cross_section': intensity_calc[:, 8],
-        'received_power': intensity_calc[:, 9]
+        'received_power': intensity_calc[:, 9],
+        # Subray simulation records
+        'subray_hit': subray_sim[:, 0].astype(bool),
+        'radius_step': subray_sim[:, 1],
+        'circle_steps': subray_sim[:, 2],
+        'circle_step': subray_sim[:, 3],
+        'divergence_angle_rad': subray_sim[:, 4]
     }
 
 
@@ -110,6 +119,7 @@ def plot_records(arec, brec, outdir):
     """
     do_incidence_angle_plots(arec, brec, outdir)
     do_by_incidence_angle_plots(arec, brec, outdir)
+    do_subray_hit_plots(arec, brec, outdir)
 
 
 def validate_record(key, rec, recid):
@@ -413,6 +423,202 @@ def _do_by_incidence_angle_plots(
     plt.close(fig)
 
 
+def do_subray_hit_subplot_hist2d(
+    fig, ax, x, y, title=None, xlabel=None, ylabel=None, bins="auto"
+):
+    if title is not None:
+        ax.set_title(title, fontsize=15)
+    if bins == "auto":
+        bins = [len(np.unique(x)), len(np.unique(y))]
+    hist2d = ax.hist2d(
+        x, y, bins=bins, cmap='viridis',
+        weights=100*np.ones_like(x)/len(x),
+        edgecolors='black'
+    )
+    fig.colorbar(hist2d[3])
+    if xlabel is not None:
+        ax.set_xlabel(xlabel, fontsize=14)
+    if ylabel is not None:
+        ax.set_ylabel(ylabel, fontsize=14)
+    ax.tick_params(axis='both', which='both', labelsize=12)
+    ax.grid('both')
+    ax.set_axisbelow(True)
+
+
+def do_subray_hit_subplot_hist(
+    fig, ax, hit, x, title=None, xlabel=None, ylabel=None, bins=7,
+    relative=False
+):
+    # TODO Rethink : Implement
+    if title is not None:
+        ax.set_title(title, fontsize=15)
+    x_hit = x[hit]
+    x_nohit = x[~hit]
+    weights = [
+        100*np.ones_like(x_hit)/len(x_hit),
+        100*np.ones_like(x_nohit)/len(x_nohit)
+    ] if relative else None
+    hist = ax.hist(
+        [x_hit, x_nohit], bins=bins, label=['hit', 'miss'], weights=weights
+    )
+    if xlabel is not None:
+        ax.set_xlabel(xlabel, fontsize=14)
+    if ylabel is not None:
+        ax.set_ylabel(ylabel, fontsize=14)
+    ax.tick_params(axis='both', which='both', labelsize=14)
+    ax.legend(loc='upper right', fontsize=12)
+    ax.grid('both')
+    ax.set_axisbelow(True)
+
+
+def do_subray_hit_plots(arec, brec, outdir):
+    # Validate subray hit data
+    if(not validate_record('subray_hit', arec, 'a') or
+        not validate_record('radius_step', arec, 'a') or
+        not validate_record('circle_steps', arec, 'a') or
+        not validate_record('circle_step', arec, 'a') or
+        not validate_record('divergence_angle_rad', arec, 'a') or
+        not validate_record('subray_hit', brec, 'b') or
+        not validate_record('radius_step', brec, 'b') or
+        not validate_record('circle_steps', brec, 'b') or
+        not validate_record('circle_step', brec, 'b') or
+        not validate_record('divergence_angle_rad', brec, 'b')
+       ):
+        print('Cannot do subray hit plots')
+        return
+
+    # Do the subray hit plots
+    fig = init_figure()  # Initialize figure
+    # CASE A
+    ax = fig.add_subplot(4, 5, 1)  # Initialize hit2Dhist on (radstep,circstep)
+    do_subray_hit_subplot_hist2d(
+        fig, ax,
+        arec['circle_step'][arec['subray_hit']],
+        arec['radius_step'][arec['subray_hit']],
+        title='Hit distribution (100%) (A)',
+    )
+    ax = fig.add_subplot(4, 5, 2)  # Initialize a hist on radius step by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, arec['subray_hit'], arec['radius_step'],
+        ylabel='Absolute'
+    )
+    ax = fig.add_subplot(4, 5, 3)  # Initialize a hist on circle steps by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, arec['subray_hit'], arec['circle_steps'],
+    )
+    ax = fig.add_subplot(4, 5, 4)  # Initialize a hist on circle step by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, arec['subray_hit'], arec['circle_step'],
+    )
+    ax = fig.add_subplot(4, 5, 5)  # Initialize a hist on div. angle by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, arec['subray_hit'],
+        1e03*arec['divergence_angle_rad']*180/np.pi,
+    )
+    ax = fig.add_subplot(4, 5, 6)  # Initialize non-hit 2D hist on (rs, cs)
+    do_subray_hit_subplot_hist2d(
+        fig, ax,
+        arec['circle_step'][~arec['subray_hit']],
+        arec['radius_step'][~arec['subray_hit']],
+        title='No-hit distribution (100%)',
+        xlabel='Circle step',
+        ylabel='Radius step'
+    )
+    ax = fig.add_subplot(4, 5, 7)  # Initialize a hist on radius step by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, arec['subray_hit'], arec['radius_step'],
+        ylabel='Relative ($100\\%$)',
+        relative=True,
+        xlabel='Radius step'
+    )
+    ax = fig.add_subplot(4, 5, 8)  # Initialize a hist on circle steps by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, arec['subray_hit'], arec['circle_steps'],
+        relative=True,
+        xlabel='Circle steps'
+    )
+    ax = fig.add_subplot(4, 5, 9)  # Initialize a hist on circle step by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, arec['subray_hit'], arec['circle_step'],
+        relative=True,
+        xlabel='Circle step'
+    )
+    ax = fig.add_subplot(4, 5, 10)  # Initialize a hist on div. angle by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, arec['subray_hit'],
+        1e03*arec['divergence_angle_rad']*180/np.pi,
+        relative=True,
+        xlabel='Divergence angle (deg $\\times 10^{-3}$)'
+    )
+    # CASE B
+    ax = fig.add_subplot(4, 5, 11)  # Initialize hit2Dhist on (radstep,circstep)
+    do_subray_hit_subplot_hist2d(
+        fig, ax,
+        brec['circle_step'][brec['subray_hit']],
+        brec['radius_step'][brec['subray_hit']],
+        title='Hit distribution (100%) (B)',
+    )
+    ax = fig.add_subplot(4, 5, 12)  # Initialize a hist on radius step by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, brec['subray_hit'], brec['radius_step'],
+        ylabel='Absolute'
+    )
+    ax = fig.add_subplot(4, 5, 13)  # Initialize a hist on circle steps by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, brec['subray_hit'], brec['circle_steps'],
+    )
+    ax = fig.add_subplot(4, 5, 14)  # Initialize a hist on circle step by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, brec['subray_hit'], brec['circle_step'],
+    )
+    ax = fig.add_subplot(4, 5, 15)  # Initialize a hist on div. angle by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, brec['subray_hit'],
+        1e03*brec['divergence_angle_rad']*180/np.pi,
+    )
+    ax = fig.add_subplot(4, 5, 16)  # Initialize non-hit 2D hist on (rs, cs)
+    do_subray_hit_subplot_hist2d(
+        fig, ax,
+        brec['circle_step'][~brec['subray_hit']],
+        brec['radius_step'][~brec['subray_hit']],
+        title='No-hit distribution (100%)',
+        xlabel='Circle step',
+        ylabel='Radius step'
+    )
+    ax = fig.add_subplot(4, 5, 17)  # Initialize a hist on radius step by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, brec['subray_hit'], brec['radius_step'],
+        ylabel='Relative ($100\\%$)',
+        relative=True,
+        xlabel='Radius step'
+    )
+    ax = fig.add_subplot(4, 5, 18)  # Initialize a hist on circle steps by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, brec['subray_hit'], brec['circle_steps'],
+        relative=True,
+        xlabel='Circle steps'
+    )
+    ax = fig.add_subplot(4, 5, 19)  # Initialize a hist on circle step by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, brec['subray_hit'], brec['circle_step'],
+        relative=True,
+        xlabel='Circle step'
+    )
+    ax = fig.add_subplot(4, 5, 20)  # Initialize a hist on div. angle by hit
+    do_subray_hit_subplot_hist(
+        fig, ax, brec['subray_hit'],
+        1e03*brec['divergence_angle_rad']*180/np.pi,
+        relative=True,
+        xlabel='Divergence angle (deg $\\times 10^{-3}$)'
+    )
+    # TODO Rethink : Implement
+    fig.tight_layout()
+    # Save figure to file and remove it from memory
+    fig.savefig(os.path.join(outdir, 'subray_hit_plots.png'))
+    fig.clear()
+    plt.close(fig)
+
+
 # ---   M A I N   --- #
 # ------------------- #
 if __name__ == '__main__':
diff --git a/src/dataanalytics/HDA_GlobalVars.cpp b/src/dataanalytics/HDA_GlobalVars.cpp
new file mode 100644
index 000000000..f3964f8e6
--- /dev/null
+++ b/src/dataanalytics/HDA_GlobalVars.cpp
@@ -0,0 +1,118 @@
+#ifdef DATA_ANALYTICS
+#include <HDA_GlobalVars.h>
+
+namespace helios { namespace analytics{
+
+
+// ***  GLOBAL OBJECT  *** //
+// *********************** //
+HDA_GlobalVars HDA_GV;
+
+// ***  WRITE METHODS  *** //
+// *********************** //
+HDA_GlobalVars & HDA_GlobalVars::incrementGeneratedSubraysCount(){
+    std::unique_lock<std::mutex> lock(generatedSubraysCount_mutex);
+    ++generatedSubraysCount;
+    return *this;
+}
+
+HDA_GlobalVars &
+        HDA_GlobalVars::incrementGeneratedRaysBeforeEarlyAbortCount(){
+    std::unique_lock<std::mutex> lock(
+        generatedRaysBeforeEarlyAbortCount_mutex
+    );
+    ++generatedRaysBeforeEarlyAbortCount;
+    return *this;
+}
+
+HDA_GlobalVars &
+        HDA_GlobalVars::incrementGeneratedRaysAfterEarlyAbortCount(){
+    std::unique_lock<std::mutex> lock(
+        generatedRaysAfterEarlyAbortCount_mutex
+    );
+    ++generatedRaysAfterEarlyAbortCount;
+    return *this;
+}
+
+HDA_GlobalVars &
+        HDA_GlobalVars::incrementIntensityComputationsCount() {
+    std::unique_lock<std::mutex> lock(intensityComputationsCount_mutex);
+    ++intensityComputationsCount;
+    return *this;
+}
+
+HDA_GlobalVars & HDA_GlobalVars::incrementIntersectiveSubraysCount(){
+    std::unique_lock<std::mutex> lock(intersectiveSubraysCount_mutex);
+    ++intersectiveSubraysCount;
+    return *this;
+}
+
+HDA_GlobalVars & HDA_GlobalVars::incrementNonIntersectiveSubraysCount(){
+    std::unique_lock<std::mutex> lock(nonIntersectiveSubraysCount_mutex);
+    ++nonIntersectiveSubraysCount;
+    return *this;
+}
+
+HDA_GlobalVars & HDA_GlobalVars::incrementSubrayIntersectionCount() {
+    std::unique_lock<std::mutex> lock(subrayIntersectionCount_mutex);
+    ++subrayIntersectionCount;
+    return *this;
+}
+
+HDA_GlobalVars & HDA_GlobalVars::incrementSubrayNonIntersectionCount() {
+    std::unique_lock<std::mutex> lock(subrayNonIntersectionCount_mutex);
+    ++subrayNonIntersectionCount;
+    return *this;
+}
+
+// ***  READ METHODS  *** //
+// ********************** //
+std::size_t HDA_GlobalVars::getGeneratedSubraysCount(){
+    std::unique_lock<std::mutex> lock(generatedSubraysCount_mutex);
+    return generatedSubraysCount;
+}
+
+std::size_t HDA_GlobalVars::getGeneratedRaysBeforeEarlyAbortCount(){
+    std::unique_lock<std::mutex> lock(
+        generatedRaysBeforeEarlyAbortCount_mutex
+    );
+    return generatedRaysBeforeEarlyAbortCount;
+}
+
+std::size_t HDA_GlobalVars::getGeneratedRaysAfterEarlyAbortCount(){
+    std::unique_lock<std::mutex> lock(
+        generatedRaysAfterEarlyAbortCount_mutex
+    );
+    return generatedRaysAfterEarlyAbortCount;
+}
+
+std::size_t HDA_GlobalVars::getIntersectiveSubraysCount() {
+    std::unique_lock<std::mutex> lock(intersectiveSubraysCount_mutex);
+    return intersectiveSubraysCount;
+}
+
+std::size_t HDA_GlobalVars::getNonIntersectiveSubraysCount() {
+    std::unique_lock<std::mutex> lock(nonIntersectiveSubraysCount_mutex);
+    return nonIntersectiveSubraysCount;
+}
+
+std::size_t HDA_GlobalVars::getSubrayIntersectionCount() {
+    std::unique_lock<std::mutex> lock(subrayIntersectionCount_mutex);
+    return subrayIntersectionCount;
+}
+
+std::size_t HDA_GlobalVars::getSubrayNonIntersectionCount() {
+    std::unique_lock<std::mutex> lock(subrayNonIntersectionCount_mutex);
+    return subrayNonIntersectionCount;
+}
+
+std::size_t HDA_GlobalVars::getIntensityComputationsCount(){
+    std::unique_lock<std::mutex> lock(intensityComputationsCount_mutex);
+    return intensityComputationsCount;
+}
+
+
+
+}}
+
+#endif
diff --git a/src/dataanalytics/HDA_GlobalVars.h b/src/dataanalytics/HDA_GlobalVars.h
new file mode 100644
index 000000000..024f97e91
--- /dev/null
+++ b/src/dataanalytics/HDA_GlobalVars.h
@@ -0,0 +1,79 @@
+#ifdef DATA_ANALYTICS
+#pragma once
+#include <cstdlib>
+#include <mutex>
+
+namespace helios { namespace analytics{
+
+// ***  GLOBAL OBJECT  *** //
+// *********************** //
+extern class HDA_GlobalVars HDA_GV;
+
+// TODO Rethink : Document
+class HDA_GlobalVars{
+public:
+    // ***  ATTRIBUTES  *** //
+    // ******************** //
+    std::size_t generatedRaysBeforeEarlyAbortCount;
+    std::size_t generatedRaysAfterEarlyAbortCount;
+    std::size_t generatedSubraysCount;
+    std::size_t intersectiveSubraysCount;
+    std::size_t nonIntersectiveSubraysCount;
+    std::size_t subrayIntersectionCount;
+    std::size_t subrayNonIntersectionCount;
+    std::size_t intensityComputationsCount;
+
+protected:
+    // ***  CONCURRENCY HANDLING ATTRIBUTES  *** //
+    // ***************************************** //
+    std::mutex generatedRaysBeforeEarlyAbortCount_mutex;
+    std::mutex generatedRaysAfterEarlyAbortCount_mutex;
+    std::mutex generatedSubraysCount_mutex;
+    std::mutex intersectiveSubraysCount_mutex;
+    std::mutex nonIntersectiveSubraysCount_mutex;
+    std::mutex subrayIntersectionCount_mutex;
+    std::mutex subrayNonIntersectionCount_mutex;
+    std::mutex intensityComputationsCount_mutex;
+
+public:
+    // ***  CONSTRUCTION / DESTRUCTION  *** //
+    // ************************************ //
+    HDA_GlobalVars() :
+        generatedRaysBeforeEarlyAbortCount(0),
+        generatedRaysAfterEarlyAbortCount(0),
+        generatedSubraysCount(0),
+        intersectiveSubraysCount(0),
+        nonIntersectiveSubraysCount(0),
+        subrayIntersectionCount(0),
+        subrayNonIntersectionCount(0),
+        intensityComputationsCount(0)
+    {}
+    virtual ~HDA_GlobalVars() = default;
+
+    // ***  WRITE METHODS  *** //
+    // *********************** //
+    HDA_GlobalVars & incrementGeneratedRaysBeforeEarlyAbortCount();
+    HDA_GlobalVars & incrementGeneratedRaysAfterEarlyAbortCount();
+    HDA_GlobalVars & incrementGeneratedSubraysCount();
+    HDA_GlobalVars & incrementIntersectiveSubraysCount();
+    HDA_GlobalVars & incrementNonIntersectiveSubraysCount();
+    HDA_GlobalVars & incrementSubrayIntersectionCount();
+    HDA_GlobalVars & incrementSubrayNonIntersectionCount();
+    HDA_GlobalVars & incrementIntensityComputationsCount();
+
+    // ***  READ METHODS  *** //
+    // ********************** //
+    std::size_t getGeneratedSubraysCount();
+    std::size_t getGeneratedRaysBeforeEarlyAbortCount();
+    std::size_t getGeneratedRaysAfterEarlyAbortCount();
+    std::size_t getIntersectiveSubraysCount();
+    std::size_t getNonIntersectiveSubraysCount();
+    std::size_t getSubrayIntersectionCount();
+    std::size_t getSubrayNonIntersectionCount();
+    std::size_t getIntensityComputationsCount();
+
+};
+
+}}
+
+#endif
diff --git a/src/dataanalytics/HDA_GlobalVarsReporter.h b/src/dataanalytics/HDA_GlobalVarsReporter.h
new file mode 100644
index 000000000..4f03dc431
--- /dev/null
+++ b/src/dataanalytics/HDA_GlobalVarsReporter.h
@@ -0,0 +1,55 @@
+#ifdef DATA_ANALYTICS
+#pragma once
+#include <HDA_GlobalVars.h>
+#include <util/logger/logging.hpp>
+
+#include <sstream>
+#include <string>
+
+namespace helios { namespace analytics{
+
+// TODO Rethink : Document
+class HDA_GlobalVarsReporter{
+protected:
+    // ***  ATTRIBUTES  *** //
+    // ******************** //
+    HDA_GlobalVars & gv;
+
+public:
+    // ***  CONSTRUCTION / DESTRUCTION  *** //
+    // ************************************ //
+    explicit HDA_GlobalVarsReporter(HDA_GlobalVars &gv) : gv(gv) {}
+    virtual ~HDA_GlobalVarsReporter() = default;
+
+    // ***   PRINT   *** //
+    // ***************** //
+    void print(){
+        // Initialize string stream to build the print
+        std::stringstream ss;
+
+        // Print global variables
+        ss << "HDA GLOBAL VARS REPORT:\n\n";
+        ss  << "Generated rays before early abort: "
+            << gv.getGeneratedRaysBeforeEarlyAbortCount() << "\n";
+        ss  << "Generated rays after early abort: "
+            << gv.getGeneratedRaysAfterEarlyAbortCount() << "\n";
+        ss << "Generated subrays: " << gv.getGeneratedSubraysCount() << "\n";
+        ss  << "Intersective subrays: "
+            << gv.getIntersectiveSubraysCount() << "\n";
+        ss  << "Non-intersective subrays: "
+            << gv.getNonIntersectiveSubraysCount() << "\n";
+        ss  << "Subray intersections: "
+            << gv.getSubrayIntersectionCount() << "\n";
+        ss  << "Subray non-intersections: "
+            << gv.getSubrayNonIntersectionCount() << "\n";
+        ss  << "Number of computed intensities: "
+            << gv.getIntensityComputationsCount() << "\n";
+        // Print through info logging level system
+        std::string text = ss.str();
+        logging::INFO(ss.str());
+    }
+};
+
+}}
+
+#endif
diff --git a/src/dataanalytics/HDA_OfstreamWrapper.h b/src/dataanalytics/HDA_OfstreamWrapper.h
index b9fb094d1..88df13ac2 100644
--- a/src/dataanalytics/HDA_OfstreamWrapper.h
+++ b/src/dataanalytics/HDA_OfstreamWrapper.h
@@ -4,7 +4,6 @@
 #include <string>
 #include <fstream>
 
-// TODO Rethink : Document
 
 namespace helios { namespace analytics{
 
diff --git a/src/dataanalytics/HDA_PulseRecorder.cpp b/src/dataanalytics/HDA_PulseRecorder.cpp
index 715d683a2..e8518898d 100644
--- a/src/dataanalytics/HDA_PulseRecorder.cpp
+++ b/src/dataanalytics/HDA_PulseRecorder.cpp
@@ -9,6 +9,7 @@ using namespace helios::analytics;
 bool HDA_PulseRecorder::isAnyBufferOpen(){
     bool anyOpen = false;
     anyOpen |= intensityCalc->isOpen();
+    anyOpen |= subraySim->isOpen();
     return anyOpen;
 }
 
@@ -22,12 +23,19 @@ void HDA_PulseRecorder::openBuffers(){
         sep,
         true  // vectorial flag
     );
+    subraySim = std::make_shared<HDA_RecordBuffer<std::vector<double>>>(
+        craftOutputPath("subray_sim.csv"),
+        maxSize,
+        sep,
+        true  // vectorial flag
+    );
 }
 
 void HDA_PulseRecorder::closeBuffers(){
     // Close subray buffers
     std::unique_lock<std::mutex> lock(intensityCalcMutex);
     intensityCalc->close();
+    subraySim->close();
 }
 
 
@@ -48,4 +56,19 @@ void HDA_PulseRecorder::recordIntensityCalculation(
     }
 }
 
+void HDA_PulseRecorder::recordSubraySimuilation(
+    std::vector<double> const &record
+){
+    std::unique_lock<std::mutex> lock(subraySimMutex);
+    subraySim->push(record);
+}
+
+void HDA_PulseRecorder::recordSubraySimulation(
+    std::vector<std::vector<double>> const &records
+){
+    std::unique_lock<std::mutex> lock(subraySimMutex);
+    for(std::vector<double> const & record : records){
+        subraySim->push(record);
+    }
+}
 #endif
diff --git a/src/dataanalytics/HDA_PulseRecorder.h b/src/dataanalytics/HDA_PulseRecorder.h
index b96cddbdd..fa862732d 100644
--- a/src/dataanalytics/HDA_PulseRecorder.h
+++ b/src/dataanalytics/HDA_PulseRecorder.h
@@ -23,7 +23,7 @@ class HDA_PulseRecorder : public HDA_Recorder{
     // ***  ATTRIBUTES  *** //
     // ******************** //
     /**
-     * @brief The vector which components are variables involved on a
+     * @brief The vectors which components are variables involved on a
      *  particular intensity calculation for a given subray.
      *
      * [0, 1, 2] -> \f$(x, y, z)\f$
@@ -47,12 +47,32 @@ class HDA_PulseRecorder : public HDA_Recorder{
      * [10] -> 1 if the point was captured, 0 otherwise.
      */
     std::shared_ptr<HDA_RecordBuffer<std::vector<double>>> intensityCalc;
+    /**
+     * @brief The vectors which components are variables involved on the
+     *  subray simulation.
+     *
+     * [0] -> Subray hit (0 does not hit, 1 hit)
+     *
+     * [1] -> Radius step
+     *
+     * [2] -> Circle steps
+     *
+     * [3] -> Circle step
+     *
+     * [4] -> Divergence angle (in rad)
+     */
+    std::shared_ptr<HDA_RecordBuffer<std::vector<double>>> subraySim;
 
     /**
      * @brief The mutex to handle concurrent writes to the buffers related to
      *  intensity calculation.
      */
     std::mutex intensityCalcMutex;
+    /**
+     * @brief The mutex to handle concurrent writes to the buffers related to
+     *  subray simulation.
+     */
+    std::mutex subraySimMutex;
 
 public:
     // ***  CONSTRUCTION / DESTRUCTION  *** //
@@ -106,6 +126,19 @@ class HDA_PulseRecorder : public HDA_Recorder{
     virtual void recordIntensityCalculation(
         std::vector<std::vector<double>> const &records
     );
+    /**
+     * @brief Handle all the records for the current subray simulation.
+     */
+    virtual void recordSubraySimuilation(std::vector<double> const &record);
+    /**
+     * @brief Like
+     *  HDA_PulseRecorder::recordSubraySimulation(std::vector<double>)
+     *  but receiving many records at once.
+     * @see HDA_PulseRecorder::recordSubraySimulation(std::vector<double>)
+     */
+    virtual void recordSubraySimulation(
+        std::vector<std::vector<double>> const &records
+    );
 
 };
 
diff --git a/src/main/Main.cpp b/src/main/Main.cpp
index 7644a18f8..ceeb82479 100644
--- a/src/main/Main.cpp
+++ b/src/main/Main.cpp
@@ -9,6 +9,11 @@
 #include <demo/DemoSelector.h>
 #endif
 
+#ifdef DATA_ANALYTICS
+#include <dataanalytics/HDA_GlobalVars.h>
+#include <dataanalytics/HDA_GlobalVarsReporter.h>
+#endif
+
 #include <gdal_priv.h>
 #include <boost/filesystem.hpp>
 
@@ -161,6 +166,13 @@ int main(int argc, char** argv) {
         );
     }
 
+#ifdef DATA_ANALYTICS
+    helios::analytics::HDA_GlobalVarsReporter reporter(
+        helios::analytics::HDA_GV
+    );
+    reporter.print();
+#endif
+
     // Return successful exit status code (0)
     return EXIT_SUCCESS;
 }
diff --git a/src/main/helios_version.cpp b/src/main/helios_version.cpp
index 481e9f2cb..e8e4a5cc4 100644
--- a/src/main/helios_version.cpp
+++ b/src/main/helios_version.cpp
@@ -4,7 +4,7 @@
 
 const char * HELIOS_VERSION = "1.2.0";
 
-const char * HELIOS_GIT_HASH = "9f5ba7f6";
+const char * HELIOS_GIT_HASH = "28bb2e3c";
 
 const char * getHeliosVersion(){
     return HELIOS_VERSION;
diff --git a/src/scanner/MultiScanner.cpp b/src/scanner/MultiScanner.cpp
index c12dc2e36..e60eaa2a5 100644
--- a/src/scanner/MultiScanner.cpp
+++ b/src/scanner/MultiScanner.cpp
@@ -166,12 +166,18 @@ void MultiScanner::computeSubrays(
         NoiseSource<double> &intersectionHandlingNoiseSource,
         std::map<double, double> &reflections,
         vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+       ,bool &subrayHit
+#endif
     )> handleSubray,
     vector<double> const &tMinMax,
     NoiseSource<double> &intersectionHandlingNoiseSource,
     std::map<double, double> &reflections,
     vector<RaySceneIntersection> &intersects,
     size_t const idx
+#ifdef DATA_ANALYTICS
+   ,std::shared_ptr<HDA_PulseRecorder> pulseRecorder
+#endif
 ){
     scanDevs[idx].computeSubrays(
         handleSubray,
@@ -179,6 +185,9 @@ void MultiScanner::computeSubrays(
         intersectionHandlingNoiseSource,
         reflections,
         intersects
+#ifdef DATA_ANALYTICS
+       ,pulseRecorder
+#endif
     );
 }
 
diff --git a/src/scanner/MultiScanner.h b/src/scanner/MultiScanner.h
index 5c6c578cc..f04c47341 100644
--- a/src/scanner/MultiScanner.h
+++ b/src/scanner/MultiScanner.h
@@ -165,12 +165,18 @@ class MultiScanner : public Scanner{
             NoiseSource<double> &intersectionHandlingNoiseSource,
             std::map<double, double> &reflections,
             vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+            ,bool &subrayHit
+#endif
         )> handleSubray,
         vector<double> const &tMinMax,
         NoiseSource<double> &intersectionHandlingNoiseSource,
         std::map<double, double> &reflections,
         vector<RaySceneIntersection> &intersects,
         size_t const idx
+#ifdef DATA_ANALYTICS
+       ,std::shared_ptr<HDA_PulseRecorder> pulseRecorder
+#endif
     ) override;
     /**
      * @see Scanner::initializeFullWaveform
diff --git a/src/scanner/Scanner.h b/src/scanner/Scanner.h
index 85c0fadcf..78bb769eb 100644
--- a/src/scanner/Scanner.h
+++ b/src/scanner/Scanner.h
@@ -450,12 +450,18 @@ class Scanner : public Asset {
             NoiseSource<double> &intersectionHandlingNoiseSource,
             std::map<double, double> &reflections,
             vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+           ,bool &subrayHit
+#endif
         )> handleSubray,
         vector<double> const &tMinMax,
         NoiseSource<double> &intersectionHandlingNoiseSource,
         std::map<double, double> &reflections,
         vector<RaySceneIntersection> &intersects,
         size_t const idx
+#ifdef DATA_ANALYTICS
+       ,std::shared_ptr<HDA_PulseRecorder> pulseRecorder
+#endif
     ) = 0;
 
     /**
diff --git a/src/scanner/ScanningDevice.cpp b/src/scanner/ScanningDevice.cpp
index fd516fb02..5e61ce83e 100644
--- a/src/scanner/ScanningDevice.cpp
+++ b/src/scanner/ScanningDevice.cpp
@@ -3,6 +3,10 @@
 #include <logging.hpp>
 #include <scanner/detector/AbstractDetector.h>
 #include <maths/EnergyMaths.h>
+#ifdef DATA_ANALYTICS
+#include <dataanalytics/HDA_GlobalVars.h>
+using namespace helios::analytics;
+#endif
 
 // ***  CONSTRUCTION / DESTRUCTION  *** //
 // ************************************ //
@@ -222,12 +226,21 @@ void ScanningDevice::computeSubrays(
         NoiseSource<double> &intersectionHandlingNoiseSource,
         std::map<double, double> &reflections,
         vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+       ,bool &subrayHit
+#endif
     )> handleSubray,
     std::vector<double> const &tMinMax,
     NoiseSource<double> &intersectionHandlingNoiseSource,
     std::map<double, double> &reflections,
     std::vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+   ,std::shared_ptr<HDA_PulseRecorder> pulseRecorder
+#endif
 ){
+#ifdef DATA_ANALYTICS
+    bool subrayHit;
+#endif
 
     int const beamSampleQuality = FWF_settings.beamSampleQuality;
     double const radiusStep_rad = beamDivergence_rad/beamSampleQuality;
@@ -260,7 +273,20 @@ void ScanningDevice::computeSubrays(
                 intersectionHandlingNoiseSource,
                 reflections,
                 intersects
+#ifdef DATA_ANALYTICS
+               ,subrayHit
+#endif
             );
+#ifdef DATA_ANALYTICS
+            HDA_GV.incrementGeneratedSubraysCount();
+            pulseRecorder->recordSubraySimuilation(std::vector<double>({
+                (double)subrayHit,
+                (double) radiusStep,
+                (double) circleSteps,
+                (double) circleStep,
+                subrayDivergenceAngle_rad
+            }));
+#endif
         }
     }
 }
diff --git a/src/scanner/ScanningDevice.h b/src/scanner/ScanningDevice.h
index 5f8c111a7..4e919c1d5 100644
--- a/src/scanner/ScanningDevice.h
+++ b/src/scanner/ScanningDevice.h
@@ -308,11 +308,17 @@ class ScanningDevice : public Asset {
             NoiseSource<double> &intersectionHandlingNoiseSource,
             std::map<double, double> &reflections,
             vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+           ,bool &subrayHit
+#endif
         )> handleSubray,
         std::vector<double> const &tMinMax,
         NoiseSource<double> &intersectionHandlingNoiseSource,
         std::map<double, double> &reflections,
         std::vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+       ,std::shared_ptr<HDA_PulseRecorder> pulseRecorder
+#endif
     );
     /**
      * @see Scanner::initializeFullWaveform
diff --git a/src/scanner/SingleScanner.cpp b/src/scanner/SingleScanner.cpp
index d1372e820..1a7c11b80 100644
--- a/src/scanner/SingleScanner.cpp
+++ b/src/scanner/SingleScanner.cpp
@@ -202,12 +202,18 @@ void SingleScanner::computeSubrays(
         NoiseSource<double> &intersectionHandlingNoiseSource,
         std::map<double, double> &reflections,
         vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+       ,bool &subrayHit
+#endif
     )> handleSubray,
     vector<double> const &tMinMax,
     NoiseSource<double> &intersectionHandlingNoiseSource,
     std::map<double, double> &reflections,
     vector<RaySceneIntersection> &intersects,
     size_t const idx
+#ifdef DATA_ANALYTICS
+   ,std::shared_ptr<HDA_PulseRecorder> pulseRecorder
+#endif
 ){
     scanDev.computeSubrays(
         handleSubray,
@@ -215,6 +221,9 @@ void SingleScanner::computeSubrays(
         intersectionHandlingNoiseSource,
         reflections,
         intersects
+#ifdef DATA_ANALYTICS
+       ,pulseRecorder
+#endif
     );
 }
 
diff --git a/src/scanner/SingleScanner.h b/src/scanner/SingleScanner.h
index 3776c156b..34bb88bcb 100644
--- a/src/scanner/SingleScanner.h
+++ b/src/scanner/SingleScanner.h
@@ -138,12 +138,18 @@ class SingleScanner : public Scanner{
             NoiseSource<double> &intersectionHandlingNoiseSource,
             std::map<double, double> &reflections,
             vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+           ,bool &subrayHit
+#endif
         )> handleSubray,
         vector<double> const &tMinMax,
         NoiseSource<double> &intersectionHandlingNoiseSource,
         std::map<double, double> &reflections,
         vector<RaySceneIntersection> &intersects,
         size_t const idx
+#ifdef DATA_ANALYTICS
+       ,std::shared_ptr<HDA_PulseRecorder> pulseRecorder
+#endif
     ) override;
     /**
      * @see Scanner::initializeFullWaveform
diff --git a/src/scanner/detector/FullWaveformPulseRunnable.cpp b/src/scanner/detector/FullWaveformPulseRunnable.cpp
index b9b706ce6..aac64abba 100644
--- a/src/scanner/detector/FullWaveformPulseRunnable.cpp
+++ b/src/scanner/detector/FullWaveformPulseRunnable.cpp
@@ -15,6 +15,11 @@
 #include <scanner/detector/FullWaveform.h>
 #include <scanner/Scanner.h>
 
+#ifdef DATA_ANALYTICS
+#include <dataanalytics/HDA_GlobalVars.h>
+using helios::analytics::HDA_GV;
+#endif
+
 using namespace std;
 
 // ***  CONSTANTS  *** //
@@ -48,11 +53,17 @@ void FullWaveformPulseRunnable::operator()(
 	    pulse.getOriginRef(),
 	    beamDir
     );
+#ifdef DATA_ANALYTICS
+    HDA_GV.incrementGeneratedRaysBeforeEarlyAbortCount();
+#endif
 	if (tMinMax.empty()) {
 		logging::DEBUG("Early abort - beam does not intersect with the scene");
 		scanner->setLastPulseWasHit(false, pulse.getDeviceIndex());
 		return;
 	}
+#ifdef DATA_ANALYTICS
+    HDA_GV.incrementGeneratedRaysAfterEarlyAbortCount();
+#endif
 
 	// Ray casting (find intersections)
 	map<double, double> reflections;
@@ -66,7 +77,8 @@ void FullWaveformPulseRunnable::operator()(
 	    reflections,
 	    intersects
 #ifdef DATA_ANALYTICS
-       ,calcIntensityRecords
+      ,calcIntensityRecords,
+       pulseRecorder
 #endif
     );
 
@@ -102,7 +114,8 @@ void FullWaveformPulseRunnable::computeSubrays(
     std::map<double, double> &reflections,
     vector<RaySceneIntersection> &intersects
 #ifdef DATA_ANALYTICS
-   ,std::vector<std::vector<double>> &calcIntensityRecords
+   ,std::vector<std::vector<double>> &calcIntensityRecords,
+   std::shared_ptr<HDA_PulseRecorder> pulseRecorder
 #endif
 ){
     scanner->computeSubrays(
@@ -115,6 +128,9 @@ void FullWaveformPulseRunnable::computeSubrays(
             NoiseSource<double> &intersectionHandlingNoiseSource,
             std::map<double, double> &reflections,
             vector<RaySceneIntersection> &intersects
+#ifdef DATA_ANALYTICS
+           ,bool &subrayHit
+#endif
         ) -> void {
             handleSubray(
                 _tMinMax,
@@ -126,7 +142,8 @@ void FullWaveformPulseRunnable::computeSubrays(
                 reflections,
                 intersects
 #ifdef DATA_ANALYTICS
-               ,calcIntensityRecords
+               ,subrayHit,
+                calcIntensityRecords
 #endif
             );
         },
@@ -135,6 +152,9 @@ void FullWaveformPulseRunnable::computeSubrays(
         reflections,
         intersects,
         pulse.getDeviceIndex()
+#ifdef DATA_ANALYTICS
+       ,pulseRecorder
+#endif
     );
 }
 
@@ -148,9 +168,13 @@ void FullWaveformPulseRunnable::handleSubray(
     map<double, double> &reflections,
     vector<RaySceneIntersection> &intersects
 #ifdef DATA_ANALYTICS
-   ,std::vector<std::vector<double>> &calcIntensityRecords
+   ,bool &subrayHit,
+    std::vector<std::vector<double>> &calcIntensityRecords
 #endif
 ){
+#ifdef DATA_ANALYTICS
+    subrayHit = false;
+#endif
     // Rotate around the circle:
     vector<double> tMinMax = _tMinMax;
     Rotation r2 = Rotation(Directions::forward, circleStep_rad * circleStep);
@@ -171,8 +195,12 @@ void FullWaveformPulseRunnable::handleSubray(
         );
 
         if (intersect != nullptr && intersect->prim != nullptr) {
+#ifdef DATA_ANALYTICS
+            HDA_GV.incrementSubrayIntersectionCount();
+            subrayHit = true;
+#endif
             // Incidence angle:
-            if(!scanner->isFixedIncidenceAngle()) {
+            if (!scanner->isFixedIncidenceAngle()) {
                 incidenceAngle =
                     intersect->prim->getIncidenceAngle_rad(
                         pulse.getOriginRef(),
@@ -188,7 +216,7 @@ void FullWaveformPulseRunnable::handleSubray(
             );
 
             // Distance must be inside [rangeMin, rangeMax] interval
-            if(detector->isDistanceNotInRange(distance)) continue;
+            if (detector->isDistanceNotInRange(distance)) continue;
 
             // Distance between beam's center line and intersection point:
             double const radius = sin(divergenceAngle) * distance;
@@ -196,7 +224,7 @@ void FullWaveformPulseRunnable::handleSubray(
                 distance, pulse.getDeviceIndex()
             );
             double intensity = 0.0;
-            if(intersect->prim->canComputeSigmaWithLadLut()){
+            if (intersect->prim->canComputeSigmaWithLadLut()) {
                 // LadLut based intensity computation
                 double sigma = intersect->prim->computeSigmaWithLadLut(
                     subrayDirection
@@ -204,9 +232,11 @@ void FullWaveformPulseRunnable::handleSubray(
                 intensity = scanner->calcIntensity(
                     distance, radius, sigma, pulse.getDeviceIndex()
                 );
-            }
-            else{
+            } else {
                 // Lighting-based intensity computation
+#ifdef DATA_ANALYTICS
+                HDA_GV.incrementIntensityComputationsCount();
+#endif
                 intensity = scanner->calcIntensity(
                     incidenceAngle,
                     distance,
@@ -215,14 +245,14 @@ void FullWaveformPulseRunnable::handleSubray(
                     radius,
                     pulse.getDeviceIndex()
 #ifdef DATA_ANALYTICS
-                   ,calcIntensityRecords
+                    , calcIntensityRecords
 #endif
                 );
             }
 
 
             // Intersection handling
-            if(intersect->prim->canHandleIntersections()) {
+            if (intersect->prim->canHandleIntersections()) {
                 glm::dvec3 outsideIntersectionPoint =
                     RayUtils::obtainPointAfterTraversing(
                         *intersect->prim->getAABB(),
@@ -247,15 +277,14 @@ void FullWaveformPulseRunnable::handleSubray(
                         subrayDirection
                     );
                     rayContinues = true;
-                }
-                else{ // Update distance considering noise
+                } else { // Update distance considering noise
                     distance = glm::distance(
                         ihr.getIntersectionPoint(),
                         pulse.getOriginRef()
                     );
                 }
             }
-            if(!rayContinues) { // If ray is not continuing
+            if (!rayContinues) { // If ray is not continuing
                 // Then register hit by default
                 reflections.insert(
                     pair<double, double>(distance, intensity)
@@ -268,9 +297,16 @@ void FullWaveformPulseRunnable::handleSubray(
             calcIntensityRecord[0] = intersect->point.x;
             calcIntensityRecord[1] = intersect->point.y;
             calcIntensityRecord[2] = intersect->point.z;
+        }
+        else {
+            HDA_GV.incrementSubrayNonIntersectionCount();
 #endif
         }
     }
+#ifdef DATA_ANALYTICS
+    if(subrayHit) HDA_GV.incrementIntersectiveSubraysCount();
+    else HDA_GV.incrementNonIntersectiveSubraysCount();
+#endif
 }
 
 void FullWaveformPulseRunnable::digestIntersections(
diff --git a/src/scanner/detector/FullWaveformPulseRunnable.h b/src/scanner/detector/FullWaveformPulseRunnable.h
index f796b0956..0ca586486 100644
--- a/src/scanner/detector/FullWaveformPulseRunnable.h
+++ b/src/scanner/detector/FullWaveformPulseRunnable.h
@@ -88,7 +88,8 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable {
         std::map<double, double> &reflections,
         vector<RaySceneIntersection> &intersects
 #ifdef DATA_ANALYTICS
-      ,std::vector<std::vector<double>> &calcIntensityRecords
+       ,std::vector<std::vector<double>> &calcIntensityRecords,
+        std::shared_ptr<HDA_PulseRecorder> pulseRecorder
 #endif
     );
     /**
@@ -112,7 +113,8 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable {
         std::map<double, double> &reflections,
         vector<RaySceneIntersection> &intersects
 #ifdef DATA_ANALYTICS
-       ,std::vector<std::vector<double>> &calcIntensityRecords
+       ,bool &subrayHit,
+       std::vector<std::vector<double>> &calcIntensityRecords
 #endif
     );
     /**