docs: def minimize_qes(fitted_rates, rates, weights, combs):
      """Varies QEs to minimize the deviation of the rates from the fitted_rates.
  
      Parameters
      ----------
      fitted_rates: numpy array of fitted rates from fit_angular_distribution
      rates: numpy array of rates of all PMT combinations
      weights: numpy array of weights for the squared sum
      combs: pmt combinations to use for minimization
  
      Returns
      -------
      opt_qes: optimal qe values for all PMTs
  
      """
  
      def make_quality_function(fitted_rates, rates, weights, combs):
          def quality_function(qes):
              sq_sum = 0
              for fitted_rate, comb, rate, weight  \
                      in zip(fitted_rates, combs, rates, weights):
                  sq_sum += ((rate / qes[comb[0]] / qes[comb[1]] - fitted_rate) *
                             weight)**2
              return sq_sum
  
          return quality_function
  
      qfunc = make_quality_function(fitted_rates, rates, weights, combs)
      qes = np.ones(31)
      bounds = [(0.1, 2.)] * 31
      opt_qes = optimize.minimize(qfunc, qes, bounds=bounds)
      return opt_qes