The study of QCD with two light dynamical fermions is of fundamental importance to understand the mechanism of color confinement. We present results of a numerical investigation on the order of the chiral phase transition with $N_f = 2$ by use of a novel strategy in finite size scaling analysis. We compare the critical behaviour of the specific heat, of the chiral susceptibility and of the equation of state with the possible critical behaviours. A second order transition in the O(4) and O(2) universality classes are excluded by our data and substantial evidence emerges for a first order transition. Like in most of previous works we have used the standard staggered action with $L_t = 4$: possible scaling violations and the need for further studies are discussed.