def METHOD_NAME():
    ## AcquisitionData.set_storage_scheme('mem')
    # no info printing from the engine, warnings and errors sent to stdout
    # msg_red = MessageRedirector()
    # output goes to files
    msg_red = sirf.STIR.MessageRedirector('info.txt', 'warn.txt', 'errr.txt')
    # raw data to be used as a template for the acquisition model
    acq_template = sirf.STIR.AcquisitionData(raw_data_file)
    # create image with suitable sizes
    image = acq_template.create_uniform_image()
    create_sample_image(image)
    image.write("simulated_image.hv")
    # create attenuation image
    uMap = acq_template.create_uniform_image()
    create_sample_image(uMap, attenuation = True)
    uMap.write("simulated_uMap.hv")
    # show the phantom image
    image_array = image.as_array()
    show_2D_array('Phantom image', image_array[0,:,:])
    # show the attenuation image
    uMap_array = uMap.as_array()
    show_2D_array('Attenuation image', uMap_array[0,:,:])
    # require same number slices and equal z-sampling for projection data & image
    image = image.zoom_image(zooms=(0.5, 1.0, 1.0))
    uMap = uMap.zoom_image(zooms=(0.5, 1.0, 1.0))
    # select acquisition model that implements the geometric
    # forward projection by a ray tracing matrix multiplication
    acq_model_matrix = sirf.STIR.SPECTUBMatrix();
    acq_model_matrix.set_attenuation_image(uMap) # add attenuation
    acq_model_matrix.set_keep_all_views_in_cache(True) # choose whethe to keep views in cache
    acq_model_matrix.set_up(acq_template, image)
    acq_model = sirf.STIR.AcquisitionModelUsingMatrix(acq_model_matrix)
    print('projecting image...')
    # project the image to obtain simulated acquisition data
    # data from raw_data_file is used as a template
    acq_model.set_up(acq_template, image)
    simulated_data = acq_template.get_uniform_copy()
    acq_model.forward(image, 0, 1, simulated_data)
    if output_file is not None:
        simulated_data.write(output_file)
    # show simulated acquisition data
    simulated_data_as_array = simulated_data.as_array()
    show_2D_array('Forward projection', simulated_data_as_array[0, 0,:,:])
    # create noisy data
    noisy_data = simulated_data.clone()
    noisy_data_as_array = np.random.poisson(simulated_data.as_array())
    noisy_data.fill(noisy_data_as_array)
    show_2D_array('Forward projection with added noise', noisy_data_as_array[0, 0,:,:])
    # create objective function
    obj_fun = sirf.STIR.make_Poisson_loglikelihood(noisy_data)
    obj_fun.set_acquisition_model(acq_model)
    # create OSEM reconstructor object
    num_subsets = 21 # number of subsets for OSEM reconstruction
    num_subiters = 42 #number of subiterations (i.e two full iterations)
    OSEM_reconstructor = sirf.STIR.OSMAPOSLReconstructor()
    OSEM_reconstructor.set_objective_function(obj_fun)
    OSEM_reconstructor.set_num_subsets(num_subsets)
    OSEM_reconstructor.set_num_subiterations(num_subiters)
    # create initialisation image and set up reconstructor
    init_image = make_cylindrical_FOV(image.get_uniform_copy(1))
    OSEM_reconstructor.set_up(init_image)
    # Reconstruct and show reconstructed image
    OSEM_reconstructor.reconstruct(init_image)
    out_image = OSEM_reconstructor.get_current_estimate()
    out_image_array = out_image.as_array()
    show_2D_array('Reconstructed image', out_image_array[0,:,:])