The equilibrium ensemble approach to disordered systems is used to investigate the critical behaviour of the two dimensional Ising model in presence of quenched random site dilution. The numerical transfer matrix technique in semi-infinite strips of finite width, together with phenomenological renormalization and conformal invariance, is particularly suited to put the equilibrium ensemble approach to work. A new method to extract with great precision the critical temperature of the model is proposed and applied. A more systematic finite-size scaling analysis than in previous numerical studies has been performed. A parallel investigation, along the lines of the two main scenarios currently under discussion, namely the logarithmic corrections scenario (with critical exponents fixed in the Ising universality class) versus the weak universality scenario (critical exponents varying with the degree of disorder), is carried out. In interpreting our data, maximum care is constantly taken to be open in both directions. A critical discussion shows that, still, an unambiguous discrimination between the two scenarios is not possible on the basis of the available finite size data.