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()
Compare the hazard curves coming from two different calculations.