code: def create_avg_rsp(rsp_fn, skydir, ltc, event_class, event_types, x,
                     egy, cth_bins, npts=None):
      
      if npts is None:
          npts = int(np.ceil(np.max(cth_bins[1:] - cth_bins[:-1]) / 0.05))
  
      wrsp = np.zeros((len(x), len(egy), len(cth_bins) - 1))
      exps = np.zeros((len(egy), len(cth_bins) - 1))
  
      cth_bins = utils.split_bin_edges(cth_bins, npts)
      cth = edge_to_center(cth_bins)
      ltw = ltc.get_skydir_lthist(skydir, cth_bins)
      ltw = ltw.reshape(-1, npts)
  
      for et in event_types:
          rsp = rsp_fn(event_class, et, x, egy, cth)
          aeff = create_aeff(event_class, et, egy, cth)
          rsp = rsp.reshape(wrsp.shape + (npts,))
          aeff = aeff.reshape(exps.shape + (npts,))
          wrsp += np.sum(rsp * aeff[np.newaxis, :, :, :] *
                         ltw[np.newaxis, np.newaxis, :, :], axis=-1)
          exps += np.sum(aeff * ltw[np.newaxis, :, :], axis=-1)
  
      exps_inv = np.zeros_like(exps)
      exps_inv[exps > 0] = 1./exps[exps>0]        
      wrsp *= exps_inv[np.newaxis, :, :]
      return wrsp
docs: Calculate the weighted response function.