Numerical simulations of three and four top quark hadroproduction processes are carried out in the SMEFT model framework. The simulated data are used to derive expected theoretical constraints on Wilson coefficients of relevant SMEFT operators of dimension six. Obtained limits for both cases are discussed and compared in terms of processes’ sensitivity to possible BSM contribution. Results show that operator O1ttO1tt is better constrained by the process of four top quark production, whereas other four operators O1QQO1QQ, O1QtO1Qt, O8QtO8Qt and O8QQO8QQ, are similarly constrained in three and four top quark production processes. In all cases, the expected limits taken from the simultaneous analysis of the production of three and four top quarks are strengthened. Analytical expressions for the partial amplitudes of the processes tt→tttt→tt and tˉt→tˉttˉt→tˉt caused by the operators O1ttO1tt, O1QQO1QQ, O1QtO1Qt, O8QtO8Qt, O8QQO8QQ were obtained for the first time. Based on the expressions of the obtained partial amplitudes, graphs of the perturbative unitarity boundary for the listed operators were drawn. The question of how kinematic cuts motivated by partial unitarity affect the resulting constraints on the Willson coefficients is addressed. It is shown that in all cases the limits are getting somewhat worse if such cuts are applied.