def METHOD_NAME():
    ##    PET.AcquisitionData.set_storage_scheme('memory')
    # no info printing from the engine, warnings and errors sent to stdout
    msg_red = PET.MessageRedirector()
    # Create a template Acquisition Model
    #acq_template = AcquisitionData('Siemens mMR', 1, 0, 1)
    acq_template = PET.AcquisitionData(acq_template_filename)#q.get_uniform_copy()
    # create the attenuation image
    atten_image = PET.ImageData(acq_template)
    image_size = atten_image.dimensions()
    voxel_size = atten_image.voxel_sizes()
    # create a cylindrical water phantom
    water_cyl = PET.EllipticCylinder()
    water_cyl.set_length(image_size[0]*voxel_size[0])
    water_cyl.set_radii((image_size[1]*voxel_size[1]*0.25, \
                     image_size[2]*voxel_size[2]*0.25))
    water_cyl.set_origin((image_size[0]*voxel_size[0]*0.5, 0, 0))
    # add the shape to the image
    atten_image.add_shape(water_cyl, scale = 9.687E-02)
    # z-pixel coordinate of the xy-crossection to show
    z = int(image_size[0]*0.5)
    if interactive:
        # show the phantom image
        atten_image_array = atten_image.as_array()
        show_2D_array('Attenuation image', atten_image_array[z,:,:])
    # Create the activity image
    act_image = atten_image.clone()
    act_image.fill(0.0)
    # create the activity cylinder
    act_cyl = PET.EllipticCylinder()
    act_cyl.set_length(image_size[0] * voxel_size[0])
    act_cyl.set_radii((image_size[1] * voxel_size[1] * 0.125, \
                         image_size[2] * voxel_size[2] * 0.125))
    act_cyl.set_origin((0, image_size[1] * voxel_size[1] * 0.06, \
                          image_size[2] * voxel_size[2] * 0.06))
    # add the shape to the image
    act_image.add_shape(act_cyl, scale=1)
    # z-pixel coordinate of the xy-crossection to show
    z = int(image_size[0] * 0.5)
    if interactive:
        # show the phantom image
        act_image_array = act_image.as_array()
        show_2D_array('Activity image', act_image_array[z, :, :])
    # Create the Single Scatter Simulation model
    sss = PET.SingleScatterSimulator()
    # Set the attenuation image
    sss.set_attenuation_image(atten_image)
    # set-up the scatter simulator
    sss.set_up(acq_template, act_image)
    # Simulate!
    sss_data = sss.forward(act_image)
    if interactive:
        # show simulated scatter data
        simulated_scatter_as_array = sss_data.as_array()
        show_2D_array('scatter simulation', simulated_scatter_as_array[0,0,:,:])
    sss_data.write(output_file)
    ## let's also compute the unscattered counts (at the same low resolution) and compare
    acq_model = PET.AcquisitionModelUsingRayTracingMatrix()
    asm = PET.AcquisitionSensitivityModel(atten_image, acq_model)
    acq_model.set_acquisition_sensitivity(asm)
    acq_model.set_up(acq_template, act_image)
    #unscattered_data = acq_template.get_uniform_copy()
    unscattered_data = acq_model.forward(act_image)
    simulated_unscatter_as_array = unscattered_data.as_array()
    if not interactive:
        return
    show_2D_array('unscattered simulation', simulated_unscatter_as_array[0,0,:,:])
    plt.figure()
    ax = plt.subplot(111)
    plt.plot(simulated_unscatter_as_array[0,4,0,:], label='unscattered')
    plt.plot(simulated_scatter_as_array[0,4,0,:], label='scattered')
    ax.legend()
    plt.show()