def METHOD_NAME(paramfile, output, redshift, resolution):
    """Do the work and output the file.
    This reads parameters from the MP-GenIC paramfile, generates a table of patchy reionization redshifts, and saves it.
    This table can be read by MP-Gadget. The core of the method is a correlation between large-scale overdensity and
    redshift of reionization calibrated from a radiative transfer simulation. To realise this, the code needs to know the large
    scale overdensity, which it does using FastPM.
    Arguments:
    - paramfile: genic parameter file to read from
    - output: file to save the reionization table to
    - redshift: redshift for the midpoint of reionization
    - Nmesh: sie of the particle grid to use
    """
    config = configobj.ConfigObj(infile=paramfile, configspec=GenICconfigspec, file_error=True)
    #Input sanitisation
    vtor = validate.Validator()
    config.validate(vtor)
    comm = MPI.COMM_WORLD
    logger = logging
    cm_per_mpc = 3.085678e24
    BoxSize = config["BoxSize"] * config["UnitLength_in_cm"] / cm_per_mpc
    Redshift = redshift
    Nmesh = int(BoxSize / resolution)
    # round it to 8.
    Nmesh -= Nmesh % 8
    # Top-hat filter the density field on one grid scale
    filtersize = resolution
    if comm.rank == 0:
        logger.info("output = %s", output)
        logger.info("BoxSize = %g", BoxSize)
        logger.info("Redshift = %g", Redshift)
        logger.info("Nmesh = %g", Nmesh)
    pm = ParticleMesh([Nmesh, Nmesh, Nmesh], BoxSize, comm=comm)
    mnu = numpy.array([config["MNue"], config["MNum"], config["MNut"]])
    omegacdm = config["Omega0"] - config["OmegaBaryon"] - numpy.sum(mnu)/93.14/config["HubbleParam"]**2
    cosmo = Cosmology(h=config["HubbleParam"],Omega0_cdm=omegacdm,T0_cmb=config["CMBTemperature"])
    state = get_lpt(pm,Redshift, cosmo, config["Seed"])
    real = state.to_mesh()
    logger.info("field painted")
    mean = real.cmean()
    if comm.rank == 0:
        logger.info("mean 2lpt density = %s", mean)
    real[...] /= mean
    real[...] -= 1
    cmplx = real.r2c()
    logger.info("field transformed")
    for k, _, slab in zip(cmplx.slabs.x, cmplx.slabs.i, cmplx.slabs):
        k2 = sum(kd ** 2 for kd in k)
        # tophat
        f = tophat(filtersize, k2 ** 0.5)
        slab[...] *= f
        # zreion
        slab[...] *= Bofk(k2 ** 0.5)
        slab[...] *= (1 + Redshift)
    real = cmplx.c2r()
    real[...] += Redshift
    logger.info("filters applied %d",real.size)
    mean = real.cmean()
    if comm.rank == 0:
        logger.info("zreion.mean = %s", mean)
    buffer = numpy.empty(real.size, real.dtype)
    real.ravel(out=buffer)
    if comm.rank == 0:
        logger.info("sorted for output")
    if os.path.exists(output):
        raise IOError("Refusing to write to existing file: ",output)
    with bigfile.BigFileMPI(comm, output, create=True) as ff:
        with ff.create_from_array("Zreion_Table", buffer) as bb:
            bb.attrs['BoxSize'] = BoxSize
            bb.attrs['Redshift'] = Redshift
            bb.attrs['TopHatFilterSize'] = filtersize
            bb.attrs['Nmesh'] = Nmesh
        #
        # hack: compatible with current MPGadget. This is not really needed
        # we'll remove the bins later, since BoxSize and Nmesh are known.
        with ff.create("XYZ_bins", dtype='f8', size=Nmesh) as bb:
            if comm.rank == 0:
                bins = numpy.linspace(0, BoxSize, Nmesh, dtype='f8')
                bb.write(0, bins)
    if comm.rank == 0:
        logger.info("done. written at %s", output)