docs: def compare_mean_curves(calc_ref, calc, nsigma=3):
      """
      Compare the hazard curves coming from two different calculations.
      """
      dstore_ref = datastore.read(calc_ref)
      dstore = datastore.read(calc)
      imtls = dstore_ref['oqparam'].imtls
      if dstore['oqparam'].imtls != imtls:
          raise RuntimeError('The IMTs and levels are different between '
                             'calculation %d and %d' % (calc_ref, calc))
      sitecol_ref = dstore_ref['sitecol']
      sitecol = dstore['sitecol']
      site_id_ref = {(lon, lat): sid for sid, lon, lat in zip(
          sitecol_ref.sids, sitecol_ref.lons, sitecol_ref.lats)}
      site_id = {(lon, lat): sid for sid, lon, lat in zip(
          sitecol.sids, sitecol.lons, sitecol.lats)}
      common = set(site_id_ref) & set(site_id)
      if not common:
          raise RuntimeError('There are no common sites between calculation '
                             '%d and %d' % (calc_ref, calc))
      pmap_ref = PmapGetter(dstore_ref, sids=[site_id_ref[lonlat]
                                              for lonlat in common]).get_mean()
      pmap = PmapGetter(dstore, sids=[site_id[lonlat]
                                      for lonlat in common]).get_mean()
      for lonlat in common:
          mean, std = pmap[site_id[lonlat]].array.T  # shape (2, N)
          mean_ref, std_ref = pmap_ref[site_id_ref[lonlat]].array.T
          err = numpy.sqrt(std**2 + std_ref**2)
          for imt in imtls:
              sl = imtls(imt)
              ok = (numpy.abs(mean[sl] - mean_ref[sl]) < nsigma * err[sl]).all()
              if not ok:
                  md = (numpy.abs(mean[sl] - mean_ref[sl])).max()
                  plt.title('point=%s, imt=%s, maxdiff=%.2e' % (lonlat, imt, md))
                  plt.loglog(imtls[imt], mean_ref[sl] + std_ref[sl],
                             label=str(calc_ref), color='black')
                  plt.loglog(imtls[imt], mean_ref[sl] - std_ref[sl],
                             color='black')
                  plt.loglog(imtls[imt], mean[sl] + std[sl], label=str(calc),
                             color='red')
                  plt.loglog(imtls[imt], mean[sl] - std[sl], color='red')
                  plt.legend()
                  plt.show()