The dynamical cluster-decay model (DCM), with deformation effects up to hexadecapole deformations and 'compact' orientations included, is used to calculate the fusion evaporation residue cross sections σ ER for 3-and 4-neutron emission in a hot fusion reaction 48 Ca+ 238 U→ 286 112 * at various incident energies, taking three different proton magic numbers Z = 114, 120 or 126 and N = 184 for the superheavy region. In each case, the shell corrections are obtained from an 'empirical' formula, with the corresponding liquid drop energies adjusted to give the experimental binding energies. This is done for all possible mass (and charge) fragmentations of the compound system. The DCM gives a good description of the measured fusion excitation function, σ ER (=σ 3n +σ 4n ) as a function of the compound nucleus excitation energy E * , within one parameter fitting, the neck length R(E * ). Of all the three choices of magic numbers, the fusion evaporation residue cross section remains the largest for the case of Z = 126, N = 184, and the lowest for the Z = 114, N = 184 case, thereby clearly suggesting that Z = 126, N = 184 are the strongest magic numbers (largest shell corrections), and Z = 114, N = 184 are the weakest (smallest shell corrections), with Z = 120, N = 184 lying in between. Thus, the present study seems to support the island of stability for superheavy nuclei centering around Z = 126, N = 184, rather than around Z = 114, N = 184, with Z = 120, N = 184 as the second best possibility.