## 5 Discussion
The method relies solely on the earthquake locations and is particularly suited for exploring data such as tectono-volcanic earthquakes, which are difficult to analyze using other methods due to their typically 3D distribution, low magnitude and obscured alignments. As shown in Figure 4a, the classical method requires first identifying the fault (Figure 4b, blue points), then selecting a subset of data where alignment may be evident, and finally, estimating the probable fault structure. In contrast, our method eliminates the need to select a data subset and quickly arrives at the solution. Notably, the computation time for 5000 events and 50 pivots is only tens of minutes when using an I7 PC with 16 cores. Assuming that most of the earthquakes are clustered on fault planes, the algorithm finds the principal fault using a density-driven search. It works as a Monte Carlo method applied to a subspace constrained by the pivot hypocenters. The strength of the proposed method is based on the use of a very robust norm such that when earthquakes belong to the explorer parallelepiped, they are assigned a weight 1; if not, the weight is zero, and this avoids the well-known problem of outliers in the statistical analyses. Conversely, the capability to explore the whole parameter space avoids the trapping of the solution into a local minimum. The algorithm's ease of use, lies in the use of a distance which is, due to rotation, derived from a single variable from the plane. We tested our algorithm on a synthetic dataset and on the real case of the L'Aquila 2009 seismic sequence. The obtained results for both study cases show that HYDE is a very suitable and reliable tool that is able to identify the main fault plane even in the presence of many fault segment substructures. In particular, for the real L'Aquila case, we identified, using our algorithm, that the geometry of the main fault plane is in good agreement with the independently calculated focal mechanism [32]. This validates our results, paving the way to several applications of the proposed method where focal mechanisms are not available or not well constrained to reduce uncertainties. Moreover, regarding the L'Aquila earthquake sequence, we performed a search using 1000 pivots to obtain stable fault parameters, but a search with only 200 pivots already gave us the same solution. This depends on the nature of the clustering of the earthquakes on the fault plane, which attracts the explorers to the right subspace, considerably reducing the algorithm run time. The results of the application of our algorithm on both the synthetic and the real L'Aquila case show that HYDE is able to quickly individuate the seismogenicfault planes even for complex fault geometry. Moreover, when refined seismic catalogs are not available, it is still possible to successfully use our algorithm, calibrating the thickness of the explorer parallelepiped according to the localization error magnitude. Note that the program includes a search for substructures of any dimension, simply extracting data in a local subspace. Future improvements to the HYDE algorithm may include time constraints to study the hypocenter's temporal evolution in order to explore the kinematics of the fault plains. Once the thickness and dimensions of the seismogenic fault have been estimated, it is possible to extract from the entire earthquake dataset those events associated with the fault identified by the explorer parallelepiped in order to track the temporal evolution of the ruptures. The obtained results and their validation, by means of comparison with independent studies, pave the way to several applications. The HYDE algorithm may be particularly useful where the earthquake intensity is not so high that the radiation pattern can be easily determined, as often happens for most aftershock sequences or in volcanic areas, where the percentage of earthquakes for which the radiation pattern can be defined is very low.