Using pure numpy instead of multiple nested for loops that include a calculation with changing values

The following function uses nested for loops that include a calculation used to populate a numpy array:

def foo():
  # Initializing variables
  xrange = 2
  yrange = 2
  t_array = np.random.randint(1, 5, size=(3, 5))
  const_var = 1 
  output_array = np.zeros((xrange, yrange, t_array.shape[1]))
  dist = np.zeros(t_array.shape[1])                                                                                                                                                                                  

  # Populating output_array
  for x in range(xrange):                                                                                                                                                                                                            
              for y in range(yrange):                                                                                                                                                                                                        
                  for x1 in range(const_var):                                                                                                                                                                                               
                      for y1 in range(const_var):                                                                                                                                                                                           
                          for i in range(t_array.shape[1]):                                                                                                                                                                              
                              dist[i] = (x+x1)*(y+y1)+const_var                                                                                                                               
                          output_array[x,y,np.argmin(dist)] = output_array[x,y,np.argmin(dist)] + 1
  return output_array

Is there a general way to do this in pure numpy without using for loops?

Additional explanation:

The first four for loops calculate their Cartesian product. The fifth for loop populates dist with the values from the iterators of the first four for loops. Once dist is populated, the output array gets updated at the respective indices with the values of the first four iterators and the index of the smallest value from dist. The final result should be a populated output_array.

I would like to do this in pure numpy as it would be much faster. My main problem is that the values used in the calculation for dist (x,x1,y,y1) change according to the iterations of the first four for loops and I do not know how to account for this using pure numpy.