/* $Id: fors_calib.cc,v 1.14 2013/10/24 16:46:53 cgarcia Exp $
 *
 * This file is part of the FORS Data Reduction Pipeline
 * Copyright (C) 2002-2010 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/*
 * $Author: cgarcia $
 * $Date: 2013/10/24 16:46:53 $
 * $Revision: 1.14 $
 * $Name:  $
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cmath>
#include <limits>
#include <string>
#include <sstream>
#include <exception>
#include <math.h>
#include <string.h>
#include <cpl.h>
#include <moses.h>
#include <fors_image.h>
#include <fors_flat_normalise.h>
#include <fors_detected_slits.h>
#include <fors_stack.h>
#include <fors_dfs.h>
#include <fors_header.h>
#include <fors_utils.h>
#include "fors_detmodel.h"
#include "fors_ccd_config.h"
#include "fors_subtract_bias.h"
#include "fors_saturation_mos.h"
#include "fors_overscan.h"
#include "flat_combine.h"
#include "fiera_config.h"
#include "fors_grism.h"
#include "fors_bpm.h"

struct fors_calib_config 
{
    /* Input parameters */
    double        dispersion;
    double        peakdetection;
    int           wdegree;
    int           wradius;
    double        wreject;
    int           wmode;
    int           wmosmode;
    int           cdegree;
    int           cmode;
    double        startwavelength;
    double        endwavelength;
    int           slit_ident;
    int           spa_polyorder;
    int           disp_nknots;
    int           sradius;
    int           dradius;
    float         fit_threshold;
    const char   *stack_method;
    int           min_reject;
    int           max_reject;
    double        klow;
    double        khigh;
    int           kiter;
    const char   *ignore_lines;
    const char   *used_linesets;
    double        nonlinear_level;
    double        max_nonlinear_ratio;
};

static int fors_calib_create(cpl_plugin *);
static int fors_calib_exec(cpl_plugin *);
static int fors_calib_destroy(cpl_plugin *);
static int fors_calib(cpl_parameterlist *, cpl_frameset *);
int fors_calib_retrieve_input_param(cpl_parameterlist * parlist, 
                                     cpl_frameset * frameset,
                                     fors_calib_config * config);
cpl_vector * fors_calib_get_reference_lines(cpl_frameset * frameset, 
                                            const char * arctag,
                                            const char * ignore_lines,
                                            const char * used_linesets);
std::auto_ptr<mosca::image> fors_calib_flat_mos_create_master_flat
(fors::calibrated_slits& calibrated_slits, 
 const mosca::wavelength_calibration& wave_cal,
 const mosca::grism_config& grism_cfg,
 fors_image *master_bias, const cpl_frame * bias_frame,
 struct fors_calib_config& config, cpl_frameset * frameset,
 const char * flat_tag, 
 double nonlinear_level, double max_nonlinear_ratio,
 std::vector<std::vector<double> >& slit_sat_ratio,
 std::vector<std::vector<int> >& slit_sat_count,
 cpl_mask **& nonlinear_flat_masks,
 cpl_mask **& saturated_flat_masks,
 std::auto_ptr<mosca::fiera_config>& ccd_config);
int fors_calib_flat_mos_normalise
(std::auto_ptr<mosca::image>& master_flat_d,
 const mosca::wavelength_calibration& wave_cal,
 cpl_table * slits, cpl_table * polytraces, cpl_image * coordinate, 
 struct fors_calib_config& config,
 std::auto_ptr<mosca::image>& norm_flat,
 cpl_image ** wave_profiles,
 std::vector<float>& sed_norm);
int fors_calib_flat_mos_rect_mapped
(std::auto_ptr<mosca::image>& master_flat_d, 
 std::auto_ptr<mosca::image>& norm_flat,
 cpl_table * slits,
 cpl_table *idscoeff, cpl_table * polytraces, 
 double reference, struct fors_calib_config& config,
 cpl_image *& mapped_flat, 
 cpl_image *& mapped_nflat);
cpl_mask * fors_calib_mask_rect_mapped
(cpl_mask * mask,
 cpl_table * slits, cpl_table *idscoeff, cpl_table * polytraces,
 double reference, struct fors_calib_config& config);
int fors_calib_flats_save
(std::auto_ptr<mosca::image>& master_flat_d,
 cpl_image * flat_mask,
 std::auto_ptr<mosca::image>& norm_flat,
 cpl_image * mapped_flat,  cpl_image * mapped_nflat,
 const fors::detected_slits detected_slits,
 const std::vector<std::vector<double> >& slit_sat_ratio,
 const std::vector<std::vector<int> >& slit_sat_count,
 struct fors_calib_config& config,
 cpl_frameset * frameset, const char * flat_tag, 
 const char * master_screen_flat_tag, const char * master_norm_flat_tag, 
 const char * mapped_screen_flat_tag, const char * mapped_norm_flat_tag, 
 cpl_parameterlist * parlist, const cpl_frame * ref_flat_frame,
 const mosca::ccd_config& ccd_config);
void fors_calib_qc_saturation
(cpl_propertylist * header, const fors::detected_slits detected_slits,
 const std::vector<std::vector<double> >& slit_sat_ratio,
 const std::vector<std::vector<int> >& slit_sat_count);



static char fors_calib_description[] =
"This recipe is used to identify reference lines on LSS, MOS and MXU arc lamp\n"
"exposures, and trace the spectral edges on the corresponding flat field\n"
"exposures. This information is used to determine the spectral extraction\n"
"mask to be applied in the scientific data reduction, performed with the\n"
"recipe fors_science.\n"
"This recipe accepts both FORS1 and FORS2 frames. The input arc lamp and\n"
"flat field exposures are assumed to be obtained quasi-simultaneously, so\n"
"that they would be described by exactly the same instrument distortions.\n"
"A line catalog must be specified, containing the wavelengths of the\n"
"reference arc lamp lines used for the wavelength calibration. A grism\n"
"table (typically depending on the instrument mode, and in particular on\n"
"the grism used) may also be specified: this table contains a default\n"
"recipe parameter setting to control the way spectra are extracted for\n"
"a specific instrument mode, as it is used for automatic run of the\n"
"pipeline on Paranal and in Garching. If this table is specified, it\n"
"will modify the default recipe parameter setting, with the exception of\n"
"those parameters which have been explicitly modified on the command line.\n"
"If a grism table is not specified, the input recipe parameters values\n"
"will always be read from the command line, or from an esorex configuration\n"
"file if present, or from their generic default values (that are rarely\n"
"meaningful). Finally a master bias frame must be input to this recipe.\n" 
"In the table below the MXU acronym can be read alternatively as MOS\n"
"and LSS, with the exception of CURV_COEFF_LSS, CURV_TRACES_LSS,\n"
"SPATIAL_MAP_LSS, SPECTRA_DETECTION_LSS, and and SLIT_MAP_LSS, which are\n" 
"never created. The products SPECTRA_DETECTION_MXU, SLIT_MAP_MXU, and\n" 
"DISP_RESIDUALS_MXU, are just created if the --check parameter is set to\n" 
"true. The product GLOBAL_DISTORTION_TABLE is just created if more than 12\n" 
"separate spectra are found in the CCD.\n\n"
"Input files:\n\n"
"  DO category:             Type:       Explanation:         Required:\n"
"  SCREEN_FLAT_MXU          Raw         Flat field exposures    Y\n"
"  LAMP_MXU                 Raw         Arc lamp exposure       Y\n"
"  MASTER_BIAS              Calib       Master Bias frame       Y\n"
"  MASTER_LINECAT           Calib       Line catalog            Y\n"
"  GRISM_TABLE              Calib       Grism table             .\n\n"
"Output files:\n\n"
"  DO category:             Data type:  Explanation:\n"
"  MASTER_SCREEN_FLAT_MXU   FITS image  Combined (sum) flat field\n"
"  MASTER_NORM_FLAT_MXU     FITS image  Normalised flat field\n"
"  MAPPED_SCREEN_FLAT_MXU   FITS image  Wavelength calibrated flat field\n"
"  MAPPED_NORM_FLAT_MXU     FITS image  Wavelength calibrated normalised flat\n"
"  REDUCED_LAMP_MXU         FITS image  Wavelength calibrated arc spectrum\n"
"  DISP_COEFF_MXU           FITS table  Inverse dispersion coefficients\n"
"  DISP_RESIDUALS_MXU       FITS image  Residuals in wavelength calibration\n"
"  DISP_RESIDUALS_TABLE_MXU FITS table  Residuals in wavelength calibration\n"
"  DELTA_IMAGE_MXU          FITS image  Offset vs linear wavelength calib\n"
"  WAVELENGTH_MAP_MXU       FITS image  Wavelength for each pixel on CCD\n"
"  SPECTRA_DETECTION_MXU    FITS image  Check for preliminary detection\n"
"  SLIT_MAP_MXU             FITS image  Map of central wavelength on CCD\n"
"  CURV_TRACES_MXU          FITS table  Spectral curvature traces\n"
"  CURV_COEFF_MXU           FITS table  Spectral curvature coefficients\n"
"  SPATIAL_MAP_MXU          FITS image  Spatial position along slit on CCD\n"
"  SPECTRAL_RESOLUTION_MXU  FITS table  Resolution at reference arc lines\n"
"  DETECTED_LINES_MXU       FITS table  All the lines detected in the arc\n"
"  ARC_RECTIFIED_MXU        FITS iamge  The spatial rectified arc\n"
"  SLIT_LOCATION_MXU        FITS table  Slits on product frames and CCD\n"
"  GLOBAL_DISTORTION_TABLE  FITS table  Global distortions table\n\n";

#define fors_calib_exit(message)                      \
{                                                     \
if ((const char *)message != NULL ) cpl_msg_error(recipe, message);  \
cpl_free(fiterror);                                   \
cpl_free(fitlines);                                   \
fors_image_delete(&master_bias);                        \
cpl_image_delete(coordinate);                         \
cpl_image_delete(checkwave);                          \
cpl_image_delete(flat);                               \
cpl_image_delete(master_flat);                        \
cpl_image_delete(added_flat);                         \
cpl_image_delete(rainbow);                            \
cpl_image_delete(rectified);                          \
cpl_image_delete(residual);                           \
cpl_image_delete(smo_flat);                           \
cpl_image_delete(spatial);                            \
cpl_image_delete(spectra);                            \
cpl_image_delete(wavemap);                            \
cpl_image_delete(delta);                              \
cpl_image_delete(rect_flat);                          \
cpl_image_delete(rect_nflat);                         \
cpl_image_delete(mapped_flat);                        \
cpl_image_delete(mapped_nflat);                       \
cpl_mask_delete(refmask);                             \
cpl_propertylist_delete(header);                      \
cpl_propertylist_delete(save_header);                 \
cpl_propertylist_delete(qclist);                      \
cpl_table_delete(idscoeff);                           \
cpl_table_delete(idscoeff_all);                       \
cpl_table_delete(restable);                           \
cpl_table_delete(maskslits);                          \
cpl_table_delete(overscans);                          \
cpl_table_delete(traces);                             \
cpl_table_delete(polytraces);                         \
cpl_table_delete(slits);                              \
cpl_table_delete(restab);                             \
cpl_table_delete(global);                             \
cpl_table_delete(wavelengths);                        \
cpl_vector_delete(lines);                             \
cpl_msg_indent_less();                                \
return -1;                                            \
}

#define fors_calib_exit_memcheck(message)              \
{                                                      \
if ((const char *)message !=NULL ) cpl_msg_info(recipe, message);    \
printf("free fiterror (%p)\n", fiterror);              \
cpl_free(fiterror);                                    \
printf("free fitlines (%p)\n", fitlines);              \
cpl_free(fitlines);                                    \
printf("free bias (%p)\n", bias);                      \
cpl_image_delete(bias);                                \
printf("free master_bias (%p)\n", master_bias);        \
fors_image_delete(&master_bias);                         \
printf("free coordinate (%p)\n", coordinate);          \
cpl_image_delete(coordinate);                          \
printf("free checkwave (%p)\n", checkwave);            \
cpl_image_delete(checkwave);                           \
printf("free flat (%p)\n", flat);                      \
cpl_image_delete(flat);                                \
printf("free master_flat (%p)\n", master_flat);        \
cpl_image_delete(master_flat);                         \
printf("free norm_flat (%p)\n", norm_flat);            \
cpl_image_delete(norm_flat);                           \
printf("free mapped_flat (%p)\n", mapped_flat);        \
cpl_image_delete(mapped_flat);                         \
printf("free mapped_nflat (%p)\n", mapped_nflat);      \
cpl_image_delete(mapped_nflat);                        \
printf("free rainbow (%p)\n", rainbow);                \
cpl_image_delete(rainbow);                             \
printf("free rectified (%p)\n", rectified);            \
cpl_image_delete(rectified);                           \
printf("free residual (%p)\n", residual);              \
cpl_image_delete(residual);                            \
printf("free smo_flat (%p)\n", smo_flat);              \
cpl_image_delete(smo_flat);                            \
printf("free spatial (%p)\n", spatial);                \
cpl_image_delete(spatial);                             \
printf("free spectra (%p)\n", spectra);                \
cpl_image_delete(spectra);                             \
printf("free wavemap (%p)\n", wavemap);                \
cpl_image_delete(wavemap);                             \
printf("free delta (%p)\n", delta);                    \
cpl_image_delete(delta);                               \
printf("free rect_flat (%p)\n", rect_flat);            \
cpl_image_delete(rect_flat);                           \
printf("free rect_nflat (%p)\n", rect_nflat);          \
cpl_image_delete(rect_nflat);                          \
printf("free refmask (%p)\n", refmask);                \
cpl_mask_delete(refmask);                              \
printf("free header (%p)\n", header);                  \
cpl_propertylist_delete(header);                       \
printf("free save_header (%p)\n", save_header);        \
cpl_propertylist_delete(save_header);                  \
printf("free qclist (%p)\n", qclist);                  \
cpl_propertylist_delete(qclist);                       \
printf("free grism_table (%p)\n", grism_table);        \
cpl_table_delete(grism_table);                         \
printf("free idscoeff (%p)\n", idscoeff);              \
cpl_table_delete(idscoeff);                            \
printf("free idscoeff_all (%p)\n", idscoeff_all);      \
cpl_table_delete(idscoeff_all);                        \
printf("free restable (%p)\n", restable);              \
cpl_table_delete(restable);                            \
printf("free maskslits (%p)\n", maskslits);            \
cpl_table_delete(maskslits);                           \
printf("free overscans (%p)\n", overscans);            \
cpl_table_delete(overscans);                           \
printf("free traces (%p)\n", traces);                  \
cpl_table_delete(traces);                              \
printf("free polytraces (%p)\n", polytraces);          \
cpl_table_delete(polytraces);                          \
printf("free slits (%p)\n", slits);                    \
cpl_table_delete(slits);                               \
printf("free restab (%p)\n", restab);                  \
cpl_table_delete(restab);                              \
printf("free global (%p)\n", global);                  \
cpl_table_delete(global);                              \
printf("free wavelengths (%p)\n", wavelengths);        \
cpl_table_delete(wavelengths);                         \
printf("free lines (%p)\n", lines);                    \
cpl_vector_delete(lines);                              \
return 0;                                              \
}


/**
 * @brief    Build the list of available plugins, for this module. 
 *
 * @param    list    The plugin list
 *
 * @return   0 if everything is ok, -1 otherwise
 *
 * Create the recipe instance and make it available to the application 
 * using the interface. This function is exported.
 */

int cpl_plugin_get_info(cpl_pluginlist *plist)
{
    cpl_recipe *recipe = static_cast<cpl_recipe *>(cpl_calloc(1, sizeof *recipe ));
    cpl_plugin *plugin = &recipe->interface;

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    FORS_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "fors_calib",
                    "Determination of the extraction mask",
                    fors_calib_description,
                    "Carlo Izzo",
                    PACKAGE_BUGREPORT,
    "This file is currently part of the FORS Instrument Pipeline\n"
    "Copyright (C) 2002-2010 European Southern Observatory\n\n"
    "This program is free software; you can redistribute it and/or modify\n"
    "it under the terms of the GNU General Public License as published by\n"
    "the Free Software Foundation; either version 2 of the License, or\n"
    "(at your option) any later version.\n\n"
    "This program is distributed in the hope that it will be useful,\n"
    "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
    "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
    "GNU General Public License for more details.\n\n"
    "You should have received a copy of the GNU General Public License\n"
    "along with this program; if not, write to the Free Software Foundation,\n"
    "Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n",
                    fors_calib_create,
                    fors_calib_exec,
                    fors_calib_destroy);

    cpl_pluginlist_append(plist, plugin);
    
    return 0;
}


/**
 * @brief    Setup the recipe options    
 *
 * @param    plugin  The plugin
 *
 * @return   0 if everything is ok
 *
 * Defining the command-line/configuration parameters for the recipe.
 */

static int fors_calib_create(cpl_plugin *plugin)
{
    cpl_recipe    *recipe;
    cpl_parameter *p;


    /* 
     * Check that the plugin is part of a valid recipe 
     */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
        recipe = (cpl_recipe *)plugin;
    else 
        return -1;

    /* 
     * Create the parameters list in the cpl_recipe object 
     */

    recipe->parameters = cpl_parameterlist_new(); 


    /*
     * Dispersion
     */

    p = cpl_parameter_new_value("fors.fors_calib.dispersion",
                                CPL_TYPE_DOUBLE,
                                "Expected spectral dispersion (Angstrom/pixel)",
                                "fors.fors_calib",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dispersion");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Peak detection level
     */

    p = cpl_parameter_new_value("fors.fors_calib.peakdetection",
                                CPL_TYPE_DOUBLE,
                                "Initial peak detection threshold (ADU)",
                                "fors.fors_calib",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "peakdetection");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /* 
     * Degree of wavelength calibration polynomial
     */

    p = cpl_parameter_new_value("fors.fors_calib.wdegree",
                                CPL_TYPE_INT,
                                "Degree of wavelength calibration polynomial",
                                "fors.fors_calib",
                                0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wdegree");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Reference lines search radius
     */

    p = cpl_parameter_new_value("fors.fors_calib.wradius",
                                CPL_TYPE_INT,
                                "Search radius if iterating pattern-matching "
                                "with first-guess method",
                                "fors.fors_calib",
                                4);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wradius");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Rejection threshold in dispersion relation polynomial fitting
     */

    p = cpl_parameter_new_value("fors.fors_calib.wreject",
                                CPL_TYPE_DOUBLE,
                                "Rejection threshold in dispersion "
                                "relation fit (pixel)",
                                "fors.fors_calib",
                                0.7);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wreject");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Wavelength solution interpolation (for LSS data)
     */

    p = cpl_parameter_new_value("fors.fors_calib.wmode",
                                CPL_TYPE_INT,
                                "Interpolation mode of wavelength solution "
                                "applicable to LSS-like data (0 = no "
                                "interpolation, 1 = fill gaps, 2 = global "
                                "model)",
                                "fors.fors_calib",
                                2);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wmode");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Wavelength solution interpolation (for MOS data)
     */

    p = cpl_parameter_new_value("fors.fors_calib.wmosmode",
                                CPL_TYPE_INT,
                                "Interpolation mode of wavelength solution "
                                "(0 = no interpolation, 1 = local (slit) "
                                "solution, 2 = global model)",
                                "fors.fors_calib",
                                0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wmosmode");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Catalog lines to ignore in wavelength calibration
     */

    p = cpl_parameter_new_value("fors.fors_calib.ignore_lines",
                                CPL_TYPE_STRING,
                                "Catalog lines nearest to wavelengths in this "
                                "list will be ignored for wavelength calibration",
                                "fors.fors_calib",
                                "");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ignore_lines");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Linesets to use
     */

    p = cpl_parameter_new_value("fors.fors_calib.used_linesets",
                                CPL_TYPE_STRING,
                                "Linesets to use. Valid are 'standard' and"
                                "'extended' (see column LINE_SET in the "
                                "line catalogue)",
                                "fors.fors_calib",
                                "standard");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "used_linesets");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Degree of spectral curvature polynomial
     */

    p = cpl_parameter_new_value("fors.fors_calib.cdegree",
                                CPL_TYPE_INT,
                                "Degree of spectral curvature polynomial",
                                "fors.fors_calib",
                                -1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cdegree");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Curvature solution interpolation (for MOS-like data)
     */
 
    p = cpl_parameter_new_value("fors.fors_calib.cmode",
                                CPL_TYPE_INT,
                                "Interpolation mode of curvature solution "
                                "applicable to MOS-like data (0 = no "
                                "interpolation, 1 = fill gaps, 2 = global "
                                "model)",
                                "fors.fors_calib",
                                1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cmode");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Start wavelength for spectral extraction
     */

    p = cpl_parameter_new_value("fors.fors_calib.startwavelength",
                                CPL_TYPE_DOUBLE,
                                "Start wavelength in spectral extraction",
                                "fors.fors_calib",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "startwavelength");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * End wavelength for spectral extraction
     */

    p = cpl_parameter_new_value("fors.fors_calib.endwavelength",
                                CPL_TYPE_DOUBLE,
                                "End wavelength in spectral extraction",
                                "fors.fors_calib",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "endwavelength");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Try slit identification
     */

    p = cpl_parameter_new_value("fors.fors_calib.slit_ident",
                                CPL_TYPE_BOOL,
                                "Attempt slit identification for MOS or MXU",
                                "fors.fors_calib",
                                TRUE);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "slit_ident");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Stacking method
     */
    p = cpl_parameter_new_enum("fors.fors_calib.stack_method",
                                CPL_TYPE_STRING,
                                "Frames combination method",
                                "fors.fors_calib",
                                "sum", 4,
                                "sum", "mean", "median", "ksigma");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "stack_method");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /* 
     * Threshold for the sigma clipping algorithm 
     */
    p = cpl_parameter_new_value("fors.fors_calib.ksigma",
                                CPL_TYPE_STRING,
                                "Low and high threshold in ksigma method",
                                "fors.fors_calib",
                                "-3.0,3.0");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ksigma");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /* 
     * Number of iterations for the sigma clipping algorithm 
     */
    p = cpl_parameter_new_value("fors.fors_calib.kiter",
                                CPL_TYPE_INT,
                                "Max number of iterations in ksigma method",
                                "fors.fors_calib",
                                999);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "kiter");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Number of knots in flat field fitting splines along spatial direction 
     */

    p = cpl_parameter_new_value("fors.fors_calib.s_degree",
                                CPL_TYPE_INT,
                                "Polynomial degree for the flat field fitting "
                                "along spatial direction",
                                "fors.fors_calib",
                                -1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "s_degree");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Smooth box radius for flat field along spatial direction
     * (if s_knots < 0)
     */

    p = cpl_parameter_new_value("fors.fors_calib.sradius",
                                CPL_TYPE_INT,
                                "Smooth box radius for flat field along "
                                "spatial direction (used if s_knots < 0)",
                                "fors.fors_calib",
                                10);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "sradius");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Number of knots in flat field fitting splines along dispersion direction 
     */

    p = cpl_parameter_new_value("fors.fors_calib.d_nknots",
                                CPL_TYPE_INT,
                                "Number of knots in flat field fitting "
                                "splines along dispersion direction",
                                "fors.fors_calib",
                                -1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "d_nknots");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Smooth box radius for flat field along dispersion direction
     */

    p = cpl_parameter_new_value("fors.fors_calib.dradius",
                                CPL_TYPE_INT,
                                "Smooth box radius for flat field along "
                                "dispersion direction (if d_knots < 0)",
                                "fors.fors_calib",
                                10);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dradius");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Threshold percentage for flat spline fitting with respect to the maximum
     */

    p = cpl_parameter_new_value("fors.fors_calib.fit_threshold",
                                CPL_TYPE_DOUBLE,
                                "Threshold percentage for flat spline fitting"
                                "with respect to the maximum",
                                "fors.fors_calib",
                                0.01);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "fit_threshold");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /* 
     * Saturation level of the detector
     */
    p = cpl_parameter_new_value("fors.fors_calib.nonlinear_level",
                                CPL_TYPE_DOUBLE,
                                "Level above which the detector is not linear",
                                "fors.fors_calib",
                                60000.);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "nonlinear_level");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /* 
     * Maximum allowed ratio of saturated pixels per slit 
     */
    p = cpl_parameter_new_value("fors.fors_calib.max_nonlinear_ratio",
                                CPL_TYPE_DOUBLE,
                                "Maximum allowed ratio of non-linear pixels per slit",
                                "fors.fors_calib",
                                0.2);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "max_nonlinear_ratio");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    return 0;
}


/**
 * @brief    Execute the plugin instance given by the interface
 *
 * @param    plugin  the plugin
 *
 * @return   0 if everything is ok
 */

static int fors_calib_exec(cpl_plugin *plugin)
{
    cpl_recipe  *   recipe ;
    int             status = 1;

    /* Get the recipe out of the plugin */
    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
        recipe = (cpl_recipe *)plugin ;
    else return -1 ;

    /* Issue a banner */
    fors_print_banner();

    try
    {
        status = fors_calib(recipe->parameters, recipe->frames);
    }
    catch(std::exception& ex)
    {
        cpl_msg_error(cpl_func, "Recipe error: %s", ex.what());
    }
    catch(...)
    {
        cpl_msg_error(cpl_func, "An uncaught error during recipe execution");
    }

    return status;
}


/**
 * @brief    Destroy what has been created by the 'create' function
 *
 * @param    plugin  The plugin
 *
 * @return   0 if everything is ok
 */

static int fors_calib_destroy(cpl_plugin *plugin)
{
    cpl_recipe *recipe;
    
    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
        recipe = (cpl_recipe *)plugin;
    else 
        return -1;

    cpl_parameterlist_delete(recipe->parameters); 

    return 0;
}


/**
 * @brief    Interpret the command line options and execute the data processing
 *
 * @param    parlist     The parameters list
 * @param    frameset    The set-of-frames
 *
 * @return   0 if everything is ok
 */

static int fors_calib(cpl_parameterlist *parlist, cpl_frameset *frameset)
{

    const char *recipe = "fors_calib";

    /*
     * Input parameters
     */
    struct fors_calib_config config;

    /*
     * CPL objects
     */

    fors_image       *master_bias  = NULL;
    cpl_image        *multi_bias   = NULL;
    cpl_image        *flat         = NULL;
    cpl_image        *master_flat  = NULL;
    cpl_image        *added_flat   = NULL;
    cpl_image        *trace_flat   = NULL;
    cpl_image        *smo_flat     = NULL;
    cpl_image        *spectra      = NULL;
    cpl_image        *wavemap      = NULL;
    cpl_image        *delta        = NULL;
    cpl_image        *residual     = NULL;
    cpl_image        *checkwave    = NULL;
    cpl_image        *rectified    = NULL;
    cpl_image        *dummy        = NULL;
    cpl_image        *add_dummy    = NULL;
    cpl_image        *refimage     = NULL;
    cpl_image        *coordinate   = NULL;
    cpl_image        *rainbow      = NULL;
    cpl_image        *spatial      = NULL;
    cpl_image        *rect_flat    = NULL;
    cpl_image        *rect_nflat   = NULL;
    cpl_image        *mapped_flat  = NULL;
    cpl_image        *mapped_nflat = NULL;

    cpl_mask         *refmask      = NULL;

    cpl_table        *overscans    = NULL;
    cpl_table        *wavelengths  = NULL;
    cpl_table        *idscoeff     = NULL;
    cpl_table        *idscoeff_all = NULL;
    cpl_table        *restable     = NULL;
    cpl_table        *slits        = NULL;
    cpl_table        *positions    = NULL;
    cpl_table        *maskslits    = NULL;
    cpl_table        *traces       = NULL;
    cpl_table        *polytraces   = NULL;
    cpl_table        *restab       = NULL;
    cpl_table        *global       = NULL;

    cpl_vector       *lines        = NULL;

    cpl_propertylist *header       = NULL;
    cpl_propertylist *save_header  = NULL;
    cpl_propertylist *qclist       = NULL;

    const cpl_frame  *ref_flat_frame = NULL;
    const cpl_frame  *ref_arc_frame  = NULL;
    
    /*
     * Auxiliary variables
     */

    const char *arc_tag;
    const char *flat_tag;
    const char *master_screen_flat_tag;
    const char *master_norm_flat_tag;
    const char *reduced_lamp_tag;
    const char *disp_residuals_tag;
    const char *disp_coeff_tag;
    const char *wavelength_map_tag;
    const char *spectra_detection_tag;
    const char *spectral_resolution_tag;
    const char *slit_map_tag;
    const char *curv_traces_tag;
    const char *curv_coeff_tag;
    const char *spatial_map_tag;
    const char *slit_location_tag;
    const char *global_distortion_tag = "GLOBAL_DISTORTION_TABLE";
    const char *disp_residuals_table_tag;
    const char *detected_lines_tag;
    const char *arc_rectified_tag;
    const char *delta_image_tag;
    const char *mapped_screen_flat_tag;
    const char *mapped_norm_flat_tag;
    const char *flat_disp_profile_tag;
    const char *keyname;
    int         mxu, mos, lss;
    int         treat_as_lss = 0;
    int         nslits;
    float      *data;
    double      mxpos;
    double      mean_rms;
    double      mean_rms_err;
    double      alltime;
    int         nflats;
    int         nlines;
    int         rebin;
    double     *fiterror = NULL;
    int        *fitlines = NULL;
    cpl_size    nx, ny;
    cpl_size    size_spec;
    double      ref_wave;
    double      gain;
    int         compute_central_wave;
    int         ccd_xsize, ccd_ysize;
    int         i;

    cpl_errorstate   error_prevstate = cpl_errorstate_get();


    cpl_msg_set_indentation(2);

    /* 
     * Get configuration parameters
     */
    if (cpl_frameset_count_tags(frameset, "GRISM_TABLE") > 1)
        fors_calib_exit("Too many in input: GRISM_TABLE");

    if(fors_calib_retrieve_input_param(parlist, frameset, &config) != 0)
        fors_calib_exit("Failed to read input parameters");;
    
    /* Check input parameters */

    if (config.dispersion <= 0.0)
        fors_calib_exit("Invalid spectral dispersion value");

    if (config.peakdetection <= 0.0)
        fors_calib_exit("Invalid peak detection level");
    
    if (config.wdegree < 1)
        fors_calib_exit("Invalid polynomial degree");

    if (config.wdegree > 5)
        fors_calib_exit("Max allowed polynomial degree is 5");
    

    if (config.wradius < 0)
        fors_calib_exit("Invalid search radius");

    if (config.wreject <= 0.0)
        fors_calib_exit("Invalid rejection threshold");

    if (config.wmode < 0 || config.wmode > 2)
        fors_calib_exit("Invalid wavelength solution interpolation mode");

    if (config.wmosmode < 0 || config.wmosmode > 2)
        fors_calib_exit("Invalid wavelength solution interpolation mode");

    if (config.cdegree < 1)
        fors_calib_exit("Invalid polynomial degree");

    if (config.cdegree > 5)
        fors_calib_exit("Max allowed polynomial degree is 5");

    if (config.cmode < 0 || config.cmode > 2)
        fors_calib_exit("Invalid curvature solution interpolation mode");

    if (config.startwavelength > 1.0)
        if (config.startwavelength < 3000.0 || config.startwavelength > 13000.0)
            fors_calib_exit("Invalid wavelength");

    if (config.endwavelength > 1.0) {
        if (config.endwavelength < 3000.0 || config.endwavelength > 13000.0)
            fors_calib_exit("Invalid wavelength");
        if (config.startwavelength < 1.0)
            fors_calib_exit("Invalid wavelength interval");
    }

    if (config.startwavelength > 1.0)
        if (config.endwavelength - config.startwavelength <= 0.0)
            fors_calib_exit("Invalid wavelength interval");

    std::string stack_method_str = config.stack_method; 
    if(stack_method_str != "mean" && stack_method_str != "median" && 
       stack_method_str != "ksigma" && stack_method_str != "sum")
        throw std::invalid_argument(stack_method_str+" stacking algorithm invalid");

    if (strcmp(config.stack_method, "minmax") == 0) {
        if (config.min_reject < 0)
            fors_calib_exit("Invalid number of lower rejections");
        if (config.max_reject < 0)
            fors_calib_exit("Invalid number of upper rejections");
    }

    if (strcmp(config.stack_method, "ksigma") == 0) {
        if (config.klow < 0.1)
            fors_calib_exit("Invalid lower K-sigma");
        if (config.khigh < 0.1)
            fors_calib_exit("Invalid lower K-sigma");
        if (config.kiter < 1)
            fors_calib_exit("Invalid number of iterations");
    }

    if (cpl_error_get_code())
        fors_calib_exit("Failure getting the configuration parameters");


    /* 
     * Check input set-of-frames
     */

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Check input set-of-frames:");
    cpl_msg_indent_more();

    if (!dfs_equal_keyword(frameset, "ESO INS GRIS1 ID")) 
        cpl_msg_warning(cpl_func,"Input frames are not from the same grism");

    if (!dfs_equal_keyword(frameset, "ESO INS FILT1 ID")) 
        cpl_msg_warning(cpl_func,"Input frames are not from the same filter");

    if (!dfs_equal_keyword(frameset, "ESO DET CHIP1 ID")) 
        cpl_msg_warning(cpl_func,"Input frames are not from the same chip");

    mxu = cpl_frameset_count_tags(frameset, "LAMP_MXU");
    mos = cpl_frameset_count_tags(frameset, "LAMP_MOS");
    lss = cpl_frameset_count_tags(frameset, "LAMP_LSS");

    if (mxu + mos + lss == 0)
        fors_calib_exit("Missing input arc lamp frame");

    if (mxu + mos + lss > 1)
        fors_calib_exit("Just one input arc lamp frame is allowed"); 

    if (mxu) {
        cpl_msg_info(recipe, "MXU data found");
        arc_tag                  = "LAMP_MXU";
        flat_tag                 = "SCREEN_FLAT_MXU";
        master_screen_flat_tag   = "MASTER_SCREEN_FLAT_MXU";
        master_norm_flat_tag     = "MASTER_NORM_FLAT_MXU";
        reduced_lamp_tag         = "REDUCED_LAMP_MXU";
        disp_residuals_tag       = "DISP_RESIDUALS_MXU";
        disp_coeff_tag           = "DISP_COEFF_MXU";
        wavelength_map_tag       = "WAVELENGTH_MAP_MXU";
        spectra_detection_tag    = "SPECTRA_DETECTION_MXU";
        spectral_resolution_tag  = "SPECTRAL_RESOLUTION_MXU";
        slit_map_tag             = "SLIT_MAP_MXU";
        curv_traces_tag          = "CURV_TRACES_MXU";
        curv_coeff_tag           = "CURV_COEFF_MXU";
        spatial_map_tag          = "SPATIAL_MAP_MXU";
        slit_location_tag        = "SLIT_LOCATION_MXU";
        disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_MXU";
        detected_lines_tag       = "DETECTED_LINES_MXU";
        arc_rectified_tag        = "ARC_RECTIFIED_MXU";
        delta_image_tag          = "DELTA_IMAGE_MXU";
        mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_MXU";
        mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_MXU";
        flat_disp_profile_tag    = "FLAT_SED_MXU";
    }

    if (lss) {
        cpl_msg_info(recipe, "LSS data found");
        arc_tag                  = "LAMP_LSS";
        flat_tag                 = "SCREEN_FLAT_LSS";
        master_screen_flat_tag   = "MASTER_SCREEN_FLAT_LSS";
        master_norm_flat_tag     = "MASTER_NORM_FLAT_LSS";
        reduced_lamp_tag         = "REDUCED_LAMP_LSS";
        spectral_resolution_tag  = "SPECTRAL_RESOLUTION_LSS";
        disp_residuals_tag       = "DISP_RESIDUALS_LSS";
        disp_coeff_tag           = "DISP_COEFF_LSS";
        slit_location_tag        = "SLIT_LOCATION_LSS";
        wavelength_map_tag       = "WAVELENGTH_MAP_LSS";
        slit_map_tag             = "SLIT_MAP_LSS";
        disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_LSS";
        detected_lines_tag       = "DETECTED_LINES_LSS";
        arc_rectified_tag        = "ARC_RECTIFIED_LSS";
        delta_image_tag          = "DELTA_IMAGE_LSS";
        mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_LSS";
        mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_LSS";
        flat_disp_profile_tag    = "FLAT_SED_LSS";
    }

    if (mos) {
        cpl_msg_info(recipe, "MOS data found");
        arc_tag                  = "LAMP_MOS";
        flat_tag                 = "SCREEN_FLAT_MOS";
        master_screen_flat_tag   = "MASTER_SCREEN_FLAT_MOS";
        master_norm_flat_tag     = "MASTER_NORM_FLAT_MOS";
        reduced_lamp_tag         = "REDUCED_LAMP_MOS";
        disp_residuals_tag       = "DISP_RESIDUALS_MOS";
        disp_coeff_tag           = "DISP_COEFF_MOS";
        wavelength_map_tag       = "WAVELENGTH_MAP_MOS";
        spectra_detection_tag    = "SPECTRA_DETECTION_MOS";
        spectral_resolution_tag  = "SPECTRAL_RESOLUTION_MOS";
        slit_map_tag             = "SLIT_MAP_MOS";
        curv_traces_tag          = "CURV_TRACES_MOS";
        curv_coeff_tag           = "CURV_COEFF_MOS";
        spatial_map_tag          = "SPATIAL_MAP_MOS";
        slit_location_tag        = "SLIT_LOCATION_MOS";
        disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_MOS";
        detected_lines_tag       = "DETECTED_LINES_MOS";
        arc_rectified_tag        = "ARC_RECTIFIED_MOS";
        delta_image_tag          = "DELTA_IMAGE_MOS";
        mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_MOS";
        mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_MOS";
        flat_disp_profile_tag    = "FLAT_SED_MOS";
    }

    if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") == 0)
            fors_calib_exit("Missing required input: MASTER_BIAS");

    if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") > 1)
        fors_calib_exit("Too many in input: MASTER_BIAS");

    if (cpl_frameset_count_tags(frameset, "MASTER_LINECAT") == 0)
        fors_calib_exit("Missing required input: MASTER_LINECAT");

    if (cpl_frameset_count_tags(frameset, "MASTER_LINECAT") > 1)
        fors_calib_exit("Too many in input: MASTER_LINECAT");

    nflats = cpl_frameset_count_tags(frameset, flat_tag);

    if (nflats < 1) {
        cpl_msg_error(recipe, "Missing required input: %s", flat_tag);
        fors_calib_exit(NULL);
    }

    /* 
     * Get the reference frames used to inherit all the saved products 
     */
    ref_flat_frame = cpl_frameset_find_const(frameset, flat_tag); 
    ref_arc_frame = cpl_frameset_find_const(frameset, arc_tag); 

    
    cpl_msg_indent_less();

    if (nflats > 1)
        cpl_msg_info(recipe, "Load %d flat field frames and stack them "
                     "with method \"%s\"", nflats, config.stack_method);
    else
        cpl_msg_info(recipe, "Load flat field exposure...");

    cpl_msg_indent_more();

    header = dfs_load_header(frameset, flat_tag, 0);

    if (header == NULL)
        fors_calib_exit("Cannot load flat field frame header");

    /*
     * Insert here a check on supported filters:
     */

    std::string wheel4 = 
         std::string(cpl_propertylist_get_string(header, "ESO INS OPTI9 TYPE"));
    if (cpl_error_get_code() != CPL_ERROR_NONE) {
        fors_calib_exit("Missing keyword ESO INS OPTI9 TYPE in flat header");
    }

    if (wheel4 == "FILT") {
        wheel4 = std::string(cpl_propertylist_get_string(header, 
                                                     "ESO INS OPTI9 NAME"));
        cpl_msg_error(recipe, "Unsupported filter: %s", wheel4.c_str());
        fors_calib_exit(NULL);
    }

    alltime = cpl_propertylist_get_double(header, "EXPTIME");

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit("Missing keyword EXPTIME in flat field frame header");

    cpl_propertylist_delete(header);

    for (i = 1; i < nflats; i++) {

        header = dfs_load_header(frameset, NULL, 0);

        if (header == NULL)
            fors_calib_exit("Cannot load flat field frame header");

        alltime += cpl_propertylist_get_double(header, "EXPTIME");

        if (cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit("Missing keyword EXPTIME in flat field "
                            "frame header");

        cpl_propertylist_delete(header);

    }

    /* 
     * Creating master flat
     */ 
    master_flat = dfs_load_image(frameset, flat_tag, CPL_TYPE_FLOAT, 0, 0);

    if (master_flat == NULL)
        fors_calib_exit("Cannot load flat field");

    if (nflats > 1) {
        if (strcmp(config.stack_method, "sum") == 0) {
            for (i = 1; i < nflats; i++) {
                flat = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
                if (flat) {
                    cpl_image_add(master_flat, flat);
                    cpl_image_delete(flat); flat = NULL;
                }
                else
                    fors_calib_exit("Cannot load flat field");
            }

        /***
            if (nflats > 1)
                cpl_image_divide_scalar(master_flat, nflats);
        ***/

        }
        else {
            cpl_imagelist *flatlist = NULL;
            double rflux, flux;

            /*
             * added_flat is needed for tracing (masters obtained with
             * rejections are not suitable for tracing)
             */

            added_flat = cpl_image_duplicate(master_flat);

            flatlist = cpl_imagelist_new();
            cpl_imagelist_set(flatlist, master_flat, 
                              cpl_imagelist_get_size(flatlist));
            master_flat = NULL;

            /*
             * Stacking with rejection requires normalization
             * at the same flux. We normalise according to mean
             * flux. This is equivalent to determining the
             * flux ration for each image as the average of the
             * flux ratio of all pixels weighted on the actual
             * flux of each pixel.
             */

            rflux = cpl_image_get_mean(added_flat);

            for (i = 1; i < nflats; i++) {
                flat = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
                if (flat) {
                    cpl_image_add(added_flat, flat);
                    flux = cpl_image_get_mean(flat);
                    cpl_image_multiply_scalar(flat, rflux / flux);
                    cpl_imagelist_set(flatlist, flat, 
                                      cpl_imagelist_get_size(flatlist));
                    flat = NULL;
                }
                else {
                    fors_calib_exit("Cannot load flat field");
                }
            }

            if (strcmp(config.stack_method, "mean") == 0) {
                master_flat = cpl_imagelist_collapse_create(flatlist);
            }

            if (strcmp(config.stack_method, "median") == 0) {
                master_flat = cpl_imagelist_collapse_median_create(flatlist);
            }

            if (strcmp(config.stack_method, "minmax") == 0) {
                master_flat = cpl_imagelist_collapse_minmax_create(flatlist, 
                                                                   config.min_reject,
                                                                   config.max_reject);
            }

            if (strcmp(config.stack_method, "ksigma") == 0) {
                master_flat = mos_ksigma_stack(flatlist, 
                                               config.klow, config.khigh, config.kiter, NULL);
            }

            cpl_imagelist_delete(flatlist);
        }
    }

    /*
     * Get the reference wavelength and the rebin factor along the
     * dispersion direction from the arc lamp exposure
     */

    header = dfs_load_header(frameset, arc_tag, 0);

    if (header == NULL)
        fors_calib_exit("Cannot load arc lamp header");

    ref_wave = cpl_propertylist_get_double(header, "ESO INS GRIS1 WLEN");

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit("Missing keyword ESO INS GRIS1 WLEN in arc lamp "
                        "frame header");

    if (ref_wave < 3000.0)   /* Perhaps in nanometers... */
        ref_wave *= 10;

    if (ref_wave < 3000.0 || ref_wave > 13000.0) {
        cpl_msg_error(recipe, "Invalid central wavelength %.2f read from "
                      "keyword ESO INS GRIS1 WLEN in arc lamp frame header",
                      ref_wave);
        fors_calib_exit(NULL);
    }

    cpl_msg_info(recipe, "The central wavelength is: %.2f", ref_wave);

    rebin = cpl_propertylist_get_int(header, "ESO DET WIN1 BINX");

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit("Missing keyword ESO DET WIN1 BINX in arc lamp "
                        "frame header");

    if (rebin != 1) {
        config.dispersion *= rebin;
        cpl_msg_warning(recipe, "The rebin factor is %d, and therefore the "
                        "working dispersion used is %f A/pixel", rebin, 
                        config.dispersion);
    }

    gain = cpl_propertylist_get_double(header, "ESO DET OUT1 CONAD");

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit("Missing keyword ESO DET OUT1 CONAD in arc lamp "
                        "frame header");

    cpl_msg_info(recipe, "The gain factor is: %.2f e-/ADU", gain);

    if (mos || mxu) {
        int nslits_out_det = 0;


        cpl_msg_info(recipe, "Produce mask slit position table...");
        if (mos)
            maskslits = mos_load_slits_fors_mos(header, &nslits_out_det);
        else
            maskslits = mos_load_slits_fors_mxu(header);

        /*
         * Check if all slits have the same X offset: in such case, 
         * treat the observation as a long-slit one!
         */

        mxpos = cpl_table_get_column_median(maskslits, "xtop");
        nslits = cpl_table_get_nrow(maskslits);

        treat_as_lss = fors_mos_is_lss_like(maskslits,  nslits_out_det);

        if (treat_as_lss) {
            cpl_msg_warning(recipe, "All MOS slits have the same offset: %.2f\n"
                            "The LSS data reduction strategy is applied!", 
                            mxpos);
            cpl_table_delete(maskslits); maskslits = NULL;
            if (mos) {
                master_screen_flat_tag   = "MASTER_SCREEN_FLAT_LONG_MOS";
                master_norm_flat_tag     = "MASTER_NORM_FLAT_LONG_MOS";
                disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_LONG_MOS";
                delta_image_tag          = "DELTA_IMAGE_LONG_MOS";
                spectral_resolution_tag  = "SPECTRAL_RESOLUTION_LONG_MOS";
                reduced_lamp_tag         = "REDUCED_LAMP_LONG_MOS";
                disp_coeff_tag           = "DISP_COEFF_LONG_MOS";
                detected_lines_tag       = "DETECTED_LINES_LONG_MOS";
                wavelength_map_tag       = "WAVELENGTH_MAP_LONG_MOS";
                slit_location_tag        = "SLIT_LOCATION_LONG_MOS";
                mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_LONG_MOS";
                mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_LONG_MOS";
                flat_disp_profile_tag    = "FLAT_SED_LONG_MOS";
            }
            else {
                master_screen_flat_tag   = "MASTER_SCREEN_FLAT_LONG_MXU";
                master_norm_flat_tag     = "MASTER_NORM_FLAT_LONG_MXU";
                disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_LONG_MXU";
                delta_image_tag          = "DELTA_IMAGE_LONG_MXU";
                spectral_resolution_tag  = "SPECTRAL_RESOLUTION_LONG_MXU";
                reduced_lamp_tag         = "REDUCED_LAMP_LONG_MXU";
                disp_coeff_tag           = "DISP_COEFF_LONG_MXU";
                detected_lines_tag       = "DETECTED_LINES_LONG_MXU";
                wavelength_map_tag       = "WAVELENGTH_MAP_LONG_MXU";
                slit_location_tag        = "SLIT_LOCATION_LONG_MXU";
                mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_LONG_MXU";
                mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_LONG_MXU";
                flat_disp_profile_tag    = "FLAT_SED_LONG_MXU";
            }
        }
    }

    if (config.slit_ident == 0) {
        cpl_table_delete(maskslits); maskslits = NULL;
    }


    /* Leave the header on for the next step... */


    /*
     * Remove the master bias
     */

    //TODO: Remove all this. Create a trace_flat from scratch adding all flats 
    //and removing the overscan.
    const cpl_frame * bias_frame = 
            cpl_frameset_find_const(frameset, "MASTER_BIAS");
    master_bias = fors_image_load(bias_frame);
    if (master_bias == NULL)
        fors_calib_exit("Cannot load master bias");

    cpl_msg_info(recipe, "Remove the master bias...");

    overscans = mos_load_overscans_vimos(header, 1);

    if (nflats > 1) {
        multi_bias = cpl_image_multiply_scalar_create(master_bias->data, nflats);
        dummy = mos_remove_bias(master_flat, multi_bias, overscans);
        if (added_flat)
            add_dummy = mos_remove_bias(added_flat, multi_bias, overscans);
        cpl_image_delete(multi_bias);
    }
    else {
        dummy = mos_remove_bias(master_flat, master_bias->data, overscans);
    }
    cpl_table_delete(overscans); overscans = NULL;
    cpl_image_delete(master_flat);
    master_flat = dummy; //This master flat is only used for tracing, I think

    if (master_flat == NULL)
        fors_calib_exit("Cannot remove bias from flat field");

    if (added_flat) {
        cpl_image_delete(added_flat);
        added_flat = add_dummy;

        if (added_flat == NULL)
            fors_calib_exit("Cannot remove bias from added flat field");

        trace_flat = added_flat;
    }
    else
        trace_flat = master_flat;

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Load arc lamp exposure...");
    cpl_msg_indent_more();

    /* Load arc */
    mosca::fiera_config arc_ccd_config(header);
    cpl_propertylist_delete(header); header = NULL;
    fors_image * arc_raw = fors_image_load
                (cpl_frameset_find_const(frameset, arc_tag));

    /* Update RON estimation from bias */
    cpl_propertylist * master_bias_header =
       cpl_propertylist_load(cpl_frame_get_filename(bias_frame), 0);
    fors_update_ccd_ron(arc_ccd_config, master_bias_header);
    cpl_propertylist_delete(master_bias_header);
    if(!cpl_errorstate_is_equal(error_prevstate))
        fors_calib_exit("Could not get RON from master bias "
                "(missing QC DET OUT? RON keywords)");
    
    /* Create arc variances map */
    fors_image_variance_from_detmodel(arc_raw, arc_ccd_config);
    if(!cpl_errorstate_is_equal(error_prevstate))
        fors_calib_exit("Cannot create variances map");
    
    /* Get the non-linear pixels */
    cpl_mask * nonlinear_arc_mask = 
            cpl_mask_threshold_image_create(arc_raw->data,
                         config.nonlinear_level, std::numeric_limits<double>::max());
    /* Get the A/D saturated pixels */
    cpl_mask *  saturated_arc_mask = 
            cpl_mask_threshold_image_create(arc_raw->data,
                            65535., std::numeric_limits<double>::max());
    cpl_mask * saturated_0 = 
            cpl_mask_threshold_image_create(arc_raw->data,
                            -std::numeric_limits<double>::max(), 
                            std::numeric_limits<double>::min());
    cpl_mask_or(saturated_arc_mask, saturated_0);
    cpl_mask_delete(saturated_0);

    /* Subtract overscan */
    fors_image * arc_red = fors_subtract_prescan(arc_raw, arc_ccd_config);
    if(!cpl_errorstate_is_equal(error_prevstate))
        fors_calib_exit("Cannot subtract pre/overscan");

    /* Trimm pre/overscan */
    fors_trimm_preoverscan(arc_red, arc_ccd_config);
    fors_trimm_preoverscan(nonlinear_arc_mask, arc_ccd_config);
    fors_trimm_preoverscan(saturated_arc_mask, arc_ccd_config);
    fors_image_delete(&arc_raw);
    if(!cpl_errorstate_is_equal(error_prevstate))
        fors_calib_exit("Cannot trimm pre/overscan");

    /* Subtract bias */
    cpl_msg_info(recipe, "Remove the master bias from arc...");
    fors_subtract_bias(arc_red, master_bias);
    if(!cpl_errorstate_is_equal(error_prevstate))
        fors_calib_exit("Cannot remove bias from arc lamp exposure");
    
    /* Assigning the data buffer to the spectra */
    spectra = arc_red->data;
    size_spec = cpl_image_get_size_x(spectra);

    
    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Load input line catalog...");
    cpl_msg_indent_more();

    /*
     * Get the reference lines
     */
    
    lines = fors_calib_get_reference_lines(frameset, arc_tag, 
                                           config.ignore_lines,
                                           config.used_linesets);
    if(lines == NULL)
        fors_calib_exit("Cannot get reference lines");
    nlines = cpl_vector_get_size(lines);
    

    /*
     * Start actual calibration
     */
    
    if (lss || treat_as_lss) {

        cpl_size first_row, last_row;
        cpl_size ylow, yhig;
        cpl_propertylist * wcs_header; 

        /* FIXME:
         * The LSS data calibration is still dirty: it doesn't apply
         * any spatial rectification, and only in future an external
         * spectral curvature model would be provided in input. Here
         * and there temporary solutions are adopted, such as accepting
         * the preliminary wavelength calibration.
         */


        /*
         * In the case of LSS data, extract the spectra directly
         * on the first attempt. The spectral curvature model may
         * be provided in input, in future releases.
         */

        cpl_msg_indent_less();
        cpl_msg_info(recipe, "Perform wavelength calibration...");
        cpl_msg_indent_more();

        nx = cpl_image_get_size_x(spectra);
        ny = cpl_image_get_size_y(spectra);

        wavemap = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
        idscoeff_all = cpl_table_new(ny);

        if (mos_subtract_background(spectra))
            fors_calib_exit("Cannot subtract the background");

        rectified = mos_wavelength_calibration_raw(spectra, lines, config.dispersion,
                                                   config.peakdetection, config.wradius,
                                                   config.wdegree, config.wreject, ref_wave,
                                                   &config.startwavelength,
                                                   &config.endwavelength, NULL,
                                                   NULL, idscoeff_all, wavemap, 
                                                   NULL, NULL, NULL, NULL);

        if (rectified == NULL)
            fors_calib_exit("Wavelength calibration failure.");

        if (!cpl_table_has_valid(idscoeff_all, "c0"))
            fors_calib_exit("Wavelength calibration failure.");

        cpl_image_delete(rectified); rectified = NULL;

        first_row = 0;
        while (!cpl_table_is_valid(idscoeff_all, "c0", first_row))
            first_row++;

        last_row = ny - 1;
        while (!cpl_table_is_valid(idscoeff_all, "c0", last_row))
            last_row--;

        ylow = first_row + 1;
        yhig = last_row + 1;

        if (ylow >= yhig) {
            cpl_error_reset();
            fors_calib_exit("No spectra could be detected.");
        }

        cpl_msg_info(recipe, 
                     "Spectral pattern was detected on %"CPL_SIZE_FORMAT
                     " out of %"CPL_SIZE_FORMAT" CCD rows", 
                     yhig - ylow, ny);

        dummy = cpl_image_extract(spectra, 1, ylow, nx, yhig);
        cpl_image_delete(spectra); spectra = dummy;

        ccd_ysize = (int)ny;
        ny = cpl_image_get_size_y(spectra);

        residual = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);

        fiterror = static_cast<double *>(cpl_calloc(ny, sizeof(double)));
        fitlines = static_cast<int *>(cpl_calloc(ny, sizeof(int)));
        idscoeff = cpl_table_new(ny);
        restable = cpl_table_new(nlines);

       //Table with positions of the detected lines used for wavelength calibration
       cpl_table * detected_lines = cpl_table_new(1);

       rectified = mos_wavelength_calibration_raw(spectra, lines, config.dispersion,
                                                   config.peakdetection, config.wradius,
                                                   config.wdegree, config.wreject, ref_wave,
                                                   &config.startwavelength, 
                                                   &config.endwavelength, fitlines, 
                                                   fiterror, idscoeff, NULL,
                                                   residual, restable, NULL,
                                                   detected_lines);
       fors_dfs_save_table(frameset, detected_lines, detected_lines_tag, NULL,
                           parlist, recipe, ref_arc_frame);
       cpl_table_delete(detected_lines);
       if(cpl_error_get_code() != CPL_ERROR_NONE)
           fors_calib_exit(NULL);

        if (rectified == NULL)
            fors_calib_exit("Wavelength calibration failure.");

        if (!cpl_table_has_valid(idscoeff, "c0"))
            fors_calib_exit("Wavelength calibration failure.");

        /*
         * A dummy slit locations table
         */

        slits = cpl_table_new(1);
        const char *clab[6] = {"c0", "c1", "c2", "c3", "c4", "c5"};
        cpl_table_new_column(slits, "slit_id", CPL_TYPE_INT);
        cpl_table_new_column(slits, "xtop", CPL_TYPE_DOUBLE);
        cpl_table_new_column(slits, "ytop", CPL_TYPE_DOUBLE);
        cpl_table_new_column(slits, "xbottom", CPL_TYPE_DOUBLE);
        cpl_table_new_column(slits, "ybottom", CPL_TYPE_DOUBLE);
        cpl_table_new_column(slits, "position", CPL_TYPE_INT);
        cpl_table_new_column(slits, "length", CPL_TYPE_INT);
        cpl_table_set_column_unit(slits, "xtop", "pixel");
        cpl_table_set_column_unit(slits, "ytop", "pixel");
        cpl_table_set_column_unit(slits, "xbottom", "pixel");
        cpl_table_set_column_unit(slits, "ybottom", "pixel");
        cpl_table_set_column_unit(slits, "position", "pixel");
        cpl_table_set_column_unit(slits, "length", "pixel");
        cpl_table_set_int(slits, "slit_id", 0, 0);
        int cnull; //It is not checked because the first and last row should contain valid calibrations... I think
        cpl_polynomial * top_ids = cpl_polynomial_new(1);
        for (cpl_size k = 0; k <= config.wdegree; k++) {
            double c = cpl_table_get_double(idscoeff, clab[k], 0, &cnull);
            cpl_polynomial_set_coeff(top_ids, &k, c);
        }
        double xtop = cpl_polynomial_eval_1d(top_ids, 0.0, NULL) + 0.5;
        cpl_table_set_double(slits, "xtop", 0, xtop);
        cpl_table_set_double(slits, "ytop", 0, (double)last_row+1);
        cpl_polynomial * botom_ids = cpl_polynomial_new(1);
        for (cpl_size k = 0; k <= config.wdegree; k++) {
            double c = cpl_table_get_double(idscoeff, clab[k], cpl_table_get_nrow(idscoeff) -1, &cnull);
            cpl_polynomial_set_coeff(botom_ids, &k, c);
        }
        double xbottom = cpl_polynomial_eval_1d(botom_ids, 0.0, NULL) + 0.5;
        cpl_table_set_double(slits, "xbottom", 0, xbottom);
        cpl_table_set_double(slits, "ybottom", 0, (double)first_row+1);
        cpl_table_set_int(slits, "position", 0, 0);
        cpl_table_set_int(slits, "length", 0, (int)ny);

        fors_dfs_save_table(frameset, slits, slit_location_tag, NULL,
                            parlist, recipe, ref_flat_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit(NULL);

        /*
         * A dummy tracing table
         */
        polytraces = cpl_table_new(2);
        cpl_table_new_column(polytraces, "slit_id", CPL_TYPE_INT);
        cpl_table_new_column(polytraces, "c0", CPL_TYPE_DOUBLE);
        cpl_table_new_column(polytraces, "c1", CPL_TYPE_DOUBLE);
        cpl_table_new_column(polytraces, "c2", CPL_TYPE_DOUBLE);
        
        cpl_table_set_int(polytraces, "slit_id", 0, 0);
        cpl_table_set_int(polytraces, "slit_id", 1, 0);
        cpl_table_set_double(polytraces, "c0", 0, (double)last_row);
        cpl_table_set_double(polytraces, "c0", 1, (double)first_row+1);
        cpl_table_set_double(polytraces, "c1", 0, 0.);
        cpl_table_set_double(polytraces, "c1", 1, 0.);
        cpl_table_set_double(polytraces, "c2", 0, 0.);
        cpl_table_set_double(polytraces, "c2", 1, 0.);

        fors_dfs_save_table(frameset, restable, disp_residuals_table_tag, NULL,
                            parlist, recipe, ref_arc_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit(NULL);

        cpl_table_delete(restable); restable = NULL;

        if (config.wmode) {
            cpl_image_delete(rectified); rectified = NULL;
            cpl_image_delete(wavemap); wavemap = NULL;
            mos_interpolate_wavecalib(idscoeff, wavemap, config.wmode, 2);
            mos_interpolate_wavecalib(idscoeff_all, wavemap, config.wmode, 2);
            wavemap = mos_map_idscoeff(idscoeff_all, nx, ref_wave,
                                       config.startwavelength, config.endwavelength);
            rectified = mos_wavelength_calibration(spectra, ref_wave,
                                                   config.startwavelength, 
                                                   config.endwavelength, config.dispersion, 
                                                   idscoeff, 0);
        }

        cpl_table_delete(idscoeff_all); idscoeff_all = NULL;

        cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
        cpl_table_set_column_unit(idscoeff, "error", "pixel");
        cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;

        for (i = 0; i < ny; i++)
            if (!cpl_table_is_valid(idscoeff, "c0", i))
                cpl_table_set_invalid(idscoeff, "error", i);

        delta = mos_map_pixel(idscoeff, ref_wave, config.startwavelength,
                              config.endwavelength, config.dispersion, 2);

        
        /* Get the mosca wave calib */
        mosca::wavelength_calibration wave_cal(idscoeff, ref_wave);
        
        /* Check that the wavelength solution is monotonically increasing */
        for(size_t spa_row = 0 ; spa_row < (size_t)ny; spa_row++)
            if(wave_cal.has_valid_cal((double)spa_row))
                if(!wave_cal.is_monotonical(spa_row, config.startwavelength,
                                            config.endwavelength, 
                                            config.dispersion))
                {
                    std::stringstream error_msg;
                    error_msg <<"The wavelength solution at row "<<spa_row<<
                         " does not increase monotonically, "
                         "which is physically impossible. Try with new parameters.";
                    throw std::range_error(error_msg.str());
                }

        
//%%%%%
        wcs_header = cpl_propertylist_new();
        cpl_propertylist_update_double(wcs_header, "CRPIX1", 1.0);
        cpl_propertylist_update_double(wcs_header, "CRPIX2", 1.0);
        cpl_propertylist_update_double(wcs_header, "CRVAL1",
                                       config.startwavelength + config.dispersion/2);
        cpl_propertylist_update_double(wcs_header, "CRVAL2", 1.0);
        /* cpl_propertylist_update_double(header, "CDELT1", config.dispersion);
        cpl_propertylist_update_double(header, "CDELT2", 1.0); */
        cpl_propertylist_update_double(wcs_header, "CD1_1", config.dispersion);
        cpl_propertylist_update_double(wcs_header, "CD1_2", 0.0);
        cpl_propertylist_update_double(wcs_header, "CD2_1", 0.0);
        cpl_propertylist_update_double(wcs_header, "CD2_2", 1.0);
        cpl_propertylist_update_string(wcs_header, "CTYPE1", "LINEAR");
        cpl_propertylist_update_string(wcs_header, "CTYPE2", "PIXEL");

        fors_dfs_save_image(frameset, delta, delta_image_tag,
                            wcs_header, parlist, recipe, ref_arc_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit(NULL);

        cpl_image_delete(delta); delta = NULL;

        cpl_msg_info(recipe, "Valid solutions found: %"CPL_SIZE_FORMAT
                     " out of %"CPL_SIZE_FORMAT" rows", 
                     ny - cpl_table_count_invalid(idscoeff, "c0"), ny);

        cpl_image_delete(spectra); spectra = NULL;

        mean_rms = mos_distortions_rms(rectified, lines, config.startwavelength,
                                       config.dispersion, 6, 0);

        cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);

        mean_rms = cpl_table_get_column_mean(idscoeff, "error");
        mean_rms_err = cpl_table_get_column_stdev(idscoeff, "error");

        cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)",
                     mean_rms, mean_rms * config.dispersion);

        restab = mos_resolution_table(rectified, config.startwavelength, config.dispersion,
                                      60000, lines);

        if (restab) {
            cpl_msg_info(recipe, "Mean spectral resolution: %.2f",
                  cpl_table_get_column_mean(restab, "resolution"));
            cpl_msg_info(recipe, 
                  "Mean reference lines FWHM: %.2f +/- %.2f pixel",
                  cpl_table_get_column_mean(restab, "fwhm") / config.dispersion,
                  cpl_table_get_column_mean(restab, "fwhm_rms") / config.dispersion);

            qclist = cpl_propertylist_new();


            /*
             * QC1 parameters
             */
            keyname = "QC.DID";

            if (fors_header_write_string(qclist,
                    keyname,
                    "2.0",
                    "QC1 dictionary")) {
                fors_calib_exit("Cannot write dictionary version "
                        "to QC log file");
            }

            if (mos)
                keyname = "QC.MOS.RESOLUTION";
            else
                keyname = "QC.LSS.RESOLUTION";

            if (fors_header_write_double(qclist, 
                    cpl_table_get_column_mean(restab,
                            "resolution"),
                            keyname, NULL, 
                            "Mean spectral resolution")) {
                fors_calib_exit("Cannot write mean spectral resolution to "
                        "QC log file");
            }

            if (mos)
                keyname = "QC.MOS.RESOLUTION.RMS"; 
            else
                keyname = "QC.LSS.RESOLUTION.RMS";

            if (fors_header_write_double(qclist,
                    cpl_table_get_column_stdev(restab,
                            "resolution"),
                            keyname, NULL, 
                            "Scatter of spectral resolution")) {
                fors_calib_exit("Cannot write spectral resolution scatter "
                        "to QC log file");
            }

            if (mos)
                keyname = "QC.MOS.RESOLUTION.NWAVE";
            else
                keyname = "QC.LSS.RESOLUTION.NWAVE";

            if (fors_header_write_int(qclist, cpl_table_get_nrow(restab) -
                    cpl_table_count_invalid(restab,
                            "resolution"),
                            keyname, NULL,
                            "Number of examined wavelengths "
                            "for resolution computation")) {
                fors_calib_exit("Cannot write number of lines used in "
                        "spectral resolution computation "
                        "to QC log file");
            }

            if (mos)
                keyname = "QC.MOS.RESOLUTION.MEANRMS";
            else
                keyname = "QC.LSS.RESOLUTION.MEANRMS";

            if (fors_header_write_double(qclist, 
                    cpl_table_get_column_mean(restab,
                            "resolution_rms"),
                            keyname, NULL,
                            "Mean error on spectral "
                            "resolution computation")) {
                fors_calib_exit("Cannot write mean error in "
                        "spectral resolution computation "
                        "to QC log file");
            }

            if (mos)
                keyname = "QC.MOS.RESOLUTION.NLINES";
            else
                keyname = "QC.LSS.RESOLUTION.NLINES";

            if (fors_header_write_int(qclist, 
                    cpl_table_get_column_mean(restab, "nlines") *
                    cpl_table_get_nrow(restab),
                    keyname, NULL,
                    "Number of lines for spectral "
                    "resolution computation")) {
                fors_calib_exit("Cannot write number of examined "
                        "wavelengths in spectral resolution computation "
                        "to QC log file");
            }

            fors_dfs_save_table(frameset, restab, spectral_resolution_tag, 
                                qclist, parlist, recipe, ref_arc_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_calib_exit(NULL);

            cpl_table_delete(restab); restab = NULL;
            cpl_propertylist_delete(qclist); qclist = NULL;

        }
        else
            fors_calib_exit("Cannot compute the spectral resolution table");

        cpl_vector_delete(lines); lines = NULL;
        cpl_msg_indent_less();

        /*
         * Save rectified arc lamp spectrum to disk
         */

        header = cpl_propertylist_new();
        cpl_propertylist_update_double(header, "CRPIX1", 1.0);
        cpl_propertylist_update_double(header, "CRPIX2", 1.0);
        cpl_propertylist_update_double(header, "CRVAL1", 
                                       config.startwavelength + config.dispersion/2);
        cpl_propertylist_update_double(header, "CRVAL2", 1.0);
        /* cpl_propertylist_update_double(header, "CDELT1", config.dispersion);
        cpl_propertylist_update_double(header, "CDELT2", 1.0); */
        cpl_propertylist_update_double(header, "CD1_1", config.dispersion);
        cpl_propertylist_update_double(header, "CD1_2", 0.0);
        cpl_propertylist_update_double(header, "CD2_1", 0.0);
        cpl_propertylist_update_double(header, "CD2_2", 1.0);
        cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
        cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
        cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);

        /* Rectify the arc masks */
        cpl_mask * nonlinear_arc_mask_rect_mapped =
                fors_calib_mask_rect_mapped(nonlinear_arc_mask, slits, idscoeff,
                                            polytraces,ref_wave, config);

        cpl_mask * saturated_arc_mask_rect_mapped =
                fors_calib_mask_rect_mapped(saturated_arc_mask, slits, idscoeff,
                                            polytraces,ref_wave, config);

        cpl_image * combined_arc_mask = 
                fors_bpm_create_combined_bpm(nonlinear_arc_mask_rect_mapped,
                                             saturated_arc_mask_rect_mapped);
        
        /* Save the arc */
        fors_dfs_save_image_mask(frameset, rectified, combined_arc_mask,
                                 reduced_lamp_tag, header, 
                                 parlist, recipe, ref_arc_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit(NULL);

        cpl_mask_delete(nonlinear_arc_mask);
        cpl_mask_delete(nonlinear_arc_mask_rect_mapped);
        cpl_mask_delete(saturated_arc_mask);
        cpl_mask_delete(saturated_arc_mask_rect_mapped);
        cpl_image_delete(combined_arc_mask);
        cpl_image_delete(rectified); rectified = NULL;
        cpl_propertylist_delete(header); header = NULL;

        fors_dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL, 
                            parlist, recipe, ref_arc_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit(NULL);

        header = dfs_load_header(frameset, arc_tag, 0);

        if (header == NULL)
            fors_calib_exit("Cannot reload arc lamp header");

        compute_central_wave = 0;
        if (lss) {
            /***
                if (fabs(1.0 - cpl_propertylist_get_double(header,
                                             "ESO INS SLIT WID")) < 0.05)
             ***/
            compute_central_wave = 1;
        }
        else {
            if (fabs(mxpos) < 0.05)
                compute_central_wave = 1;
        }

        /*
         * QC1 parameters
         */
        keyname = "QC.DID";

        if (fors_header_write_string(header,
                keyname,
                "2.0",
                "QC1 dictionary")) {
            fors_calib_exit("Cannot write dictionary version "
                    "to QC log file");
        }

        if (fors_header_write_double(header,
                mean_rms,
                "QC.WAVE.ACCURACY",
                "pixel",
                "Mean accuracy of wavecalib model")) {
            fors_calib_exit("Cannot write mean wavelength calibration "
                    "accuracy to QC log file");
        }

        if (fors_header_write_double(header,
                mean_rms_err,
                "QC.WAVE.ACCURACY.ERROR",
                "pixel",
                "Error on accuracy of wavecalib model")) {
            fors_calib_exit("Cannot write error on wavelength calibration "
                    "accuracy to QC log file");
        }

        if (compute_central_wave) {

            data = cpl_image_get_data_float(wavemap);

            if (lss) {
                if (fors_header_write_double(header, 
                        data[nx/2 + ccd_ysize*nx/2],
                        "QC.LSS.CENTRAL.WAVELENGTH",
                        "Angstrom", 
                        "Wavelength at CCD center")) {
                    fors_calib_exit("Cannot write central wavelength to QC "
                            "log file");
                }
            }
            else {
                if (fors_header_write_double(header, 
                        data[nx/2 + ccd_ysize*nx/2],
                        "QC.MOS.CENTRAL.WAVELENGTH",
                        "Angstrom", 
                        "Wavelength at CCD center")) {
                    fors_calib_exit("Cannot write central wavelength to QC "
                            "log file");
                }
            }
        }

        fors_dfs_save_image(frameset, wavemap, wavelength_map_tag, header,
                            parlist, recipe, ref_arc_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit(NULL);

        cpl_image_delete(wavemap); wavemap = NULL;

        cpl_propertylist_erase_regexp(header, "^ESO QC ", 0);

        cpl_propertylist_update_double(header, "CRPIX2", 1.0);
        cpl_propertylist_update_double(header, "CRVAL2", 1.0);
        /* cpl_propertylist_update_double(header, "CDELT2", 1.0); */
        cpl_propertylist_update_double(header, "CD1_1", 1.0);
        cpl_propertylist_update_double(header, "CD1_2", 0.0);
        cpl_propertylist_update_double(header, "CD2_1", 0.0);
        cpl_propertylist_update_double(header, "CD2_2", 1.0);
        cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
        cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");

        fors_dfs_save_image(frameset, residual, disp_residuals_tag, header,
                            parlist, recipe, ref_arc_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit(NULL);

        cpl_image_delete(residual); residual = NULL;

        cpl_propertylist_delete(header); header = NULL;
        
        /* Read grism configuration */
        //TODO: Add the waveref to the grism_tables
        cpl_propertylist * flat_header = dfs_load_header(frameset, flat_tag, 0);
        if (flat_header == NULL) {
            cpl_msg_error(recipe, "Cannot load header of %s frame", flat_tag);
            return 0;
        }
        double wave_ref =  cpl_propertylist_get_double(flat_header, "ESO INS GRIS1 WLEN");
        if (wave_ref < 3000.0)   /* Perhaps in nanometers... */
            wave_ref *= 10;
        cpl_propertylist_delete(flat_header);
        cpl_frameset * grism_frame = fors_frameset_extract(frameset, "GRISM_TABLE");
        std::auto_ptr<mosca::grism_config> grism_cfg = 
                fors_grism_config_from_frame(cpl_frameset_get_position(grism_frame, 0), wave_ref);
        cpl_frameset_delete(grism_frame);

        /* Get the detected slit locations */
        fors::detected_slits det_slits = 
            fors::detected_slits_from_tables(slits, polytraces, size_spec);
        
        /* Get the calibrated slits */
        fors::calibrated_slits calib_slits(det_slits, wave_cal, *grism_cfg,
                                           (size_t)nx, (size_t)ccd_ysize); 
        for(std::vector<mosca::calibrated_slit>::const_iterator 
                slit_it = calib_slits.begin();
            slit_it != calib_slits.end() ; slit_it++)
        {
            if(!slit_it->has_valid_wavecal())
                cpl_msg_warning(cpl_func, "Slit %d does not contain valid "
                        "wavelength calibration. Skipping it for master flat", 
                        slit_it->slit_id());
        }

        /* Flat field creation */
        cpl_msg_info(cpl_func, "Perform flat field combination");
        cpl_msg_indent_more();
        std::auto_ptr<mosca::image> master_flat_d;
        std::vector<std::vector<double> > slit_sat_ratio;
        std::vector<std::vector<int> > slit_sat_count;
        cpl_mask ** nonlinear_flat_masks;
        cpl_mask ** saturated_flat_masks;
        std::auto_ptr<mosca::fiera_config> ccd_config;
        master_flat_d = fors_calib_flat_mos_create_master_flat(calib_slits,
                wave_cal, *grism_cfg, master_bias, bias_frame,
                config, frameset, flat_tag, 
                config.nonlinear_level, config.max_nonlinear_ratio, 
                slit_sat_ratio, slit_sat_count,
                nonlinear_flat_masks, 
                saturated_flat_masks, ccd_config);
        if(master_flat_d.get() == 0)
            fors_calib_exit("Cannot combine flat frames");

        /* Getting combined mask */
        cpl_image * flat_mask = fors_bpm_create_combined_bpm(nonlinear_flat_masks,
                                                             saturated_flat_masks,
                                                             nflats);
        /*
         * Flat field normalisation is done directly on the master flat
         * field (without spatial rectification first). The spectral
         * curvature model may be provided in input, in future releases.
         */

        cpl_msg_indent_less();
        cpl_msg_info(recipe, "Perform LSS flat field normalisation...");
        cpl_msg_indent_more();
        std::auto_ptr<mosca::image> norm_flat(NULL);
        norm_flat.reset(new mosca::image(cpl_image_cast(master_flat_d->get_cpl_image(), CPL_TYPE_FLOAT),
                                         cpl_image_cast(master_flat_d->get_cpl_image_err(), CPL_TYPE_FLOAT), true));
        
        fors::flat_normaliser normaliser;
        normaliser.lss_normalise(*norm_flat, wave_cal, 
                config.sradius, config.dradius, 
                config.spa_polyorder, config.disp_nknots,
                config.fit_threshold);

        cpl_msg_indent_less();

        /* Getting slit dispersion profile */
        cpl_image * wave_profiles;
        wave_profiles = normaliser.get_wave_profiles_im();
        
        /* Getting the normalisation factors used in the SED */
        cpl_propertylist * sed_header = cpl_propertylist_new();
        const std::vector<float>& sed_norm = normaliser.get_wave_profiles_norm();
        for(size_t ised = 0 ; ised < sed_norm.size();ised++)
        {
            std::ostringstream norm_key;
            norm_key<< "ESO QC FLAT SED"<<ised+1<<" NORM ";
            cpl_propertylist_append_float(sed_header, norm_key.str().c_str(), 
                                          sed_norm[ised]);
        }

        /* Saving slit dispersion profiles */
        std::ostringstream prof_filename_oss;
        prof_filename_oss << flat_disp_profile_tag << ".fits";
        std::string prof_filename = prof_filename_oss.str();
        std::transform(prof_filename.begin(), prof_filename.end(), prof_filename.begin(), ::tolower);
        cpl_propertylist_append_string(sed_header, CPL_DFS_PRO_TYPE, "REDUCED");
        cpl_propertylist_append_string(sed_header, CPL_DFS_PRO_CATG, flat_disp_profile_tag);
        cpl_dfs_save_image(frameset, NULL, parlist, frameset, 
                           ref_flat_frame, wave_profiles, CPL_BPP_IEEE_FLOAT,
                           recipe, sed_header, NULL, PACKAGE "/" PACKAGE_VERSION,  
                           prof_filename.c_str());
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit("Cannot save flats");
        cpl_propertylist_delete(sed_header);
        cpl_image_delete(wave_profiles);


        /* Computing mapped flats */
        cpl_msg_info(recipe, "Perform mapping of flats...");
        cpl_image * master_flat_f = cpl_image_cast(master_flat_d->get_cpl_image(), 
                                                   CPL_TYPE_FLOAT);
        cpl_image * master_flat_rect = cpl_image_extract(master_flat_f, 1, ylow, nx, yhig);

        mapped_flat = mos_wavelength_calibration(master_flat_rect , ref_wave,
                                      config.startwavelength, config.endwavelength,
                                      config.dispersion, idscoeff, 0);


        if(norm_flat.get() != NULL)
        {
            cpl_image * norm_flat_f = cpl_image_cast(norm_flat->get_cpl_image(), 
                                                     CPL_TYPE_FLOAT);
            cpl_image * norm_flat_rect = cpl_image_extract(norm_flat_f, 1, ylow, nx, yhig);

            mapped_nflat = mos_wavelength_calibration(norm_flat_rect, ref_wave,
                    config.startwavelength, config.endwavelength,
                    config.dispersion, idscoeff, 0);
            cpl_image_delete(norm_flat_rect); norm_flat_rect = NULL;
            cpl_image_delete(norm_flat_f); norm_flat_f = NULL;
        }

        /* Saving flats */
        if(fors_calib_flats_save(master_flat_d, flat_mask,
                norm_flat, mapped_flat, mapped_nflat, det_slits,
                slit_sat_ratio, slit_sat_count, config,
                frameset, flat_tag, master_screen_flat_tag, master_norm_flat_tag,
                mapped_screen_flat_tag, mapped_norm_flat_tag,
                parlist, ref_flat_frame, *ccd_config) != 0)
            fors_calib_exit("Cannot save flats");

        cpl_image_delete(master_flat_f); master_flat_f = NULL;
        cpl_image_delete(master_flat_rect); master_flat_rect = NULL;
        cpl_image_delete(mapped_flat); mapped_flat = NULL;
        cpl_image_delete(mapped_nflat); mapped_nflat = NULL;
        cpl_propertylist_delete(wcs_header); wcs_header = NULL;
        cpl_image_delete(master_flat); master_flat = NULL;
        cpl_image_delete(flat_mask); 
        cpl_table_delete(idscoeff); idscoeff = NULL;
        cpl_propertylist_delete(save_header); save_header = NULL;
        cpl_table_delete(slits); slits = NULL;
        cpl_table_delete(polytraces); polytraces = NULL;
        fors_image_delete(&master_bias); master_bias = NULL;
        for (size_t i_flat = 0; i_flat < nflats; i_flat++)
        {
            cpl_mask_delete(nonlinear_flat_masks[i_flat]);
            cpl_mask_delete(saturated_flat_masks[i_flat]);
        }
        cpl_free(nonlinear_flat_masks);
        cpl_free(saturated_flat_masks);

        return 0;         /* Successful LSS data reduction */

    }   /* End of LSS data reduction section */


    /*
     * Here the MOS and MXU calibration is carried out.
     */

    /*
     * Detecting spectra on the CCD
     */

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Detecting spectra on CCD...");
    cpl_msg_indent_more();

    ccd_xsize = nx = cpl_image_get_size_x(spectra);
    ccd_ysize = ny = cpl_image_get_size_y(spectra);

    refmask = cpl_mask_new(nx, ny);

    if (mos_subtract_background(spectra))
	fors_calib_exit("Cannot subtract the background");
 
    checkwave = mos_wavelength_calibration_raw(spectra, lines, 
                      config.dispersion, config.peakdetection, config.wradius, 
                      config.wdegree, config.wreject, ref_wave,
                      &config.startwavelength, &config.endwavelength,
                      NULL, NULL, NULL, NULL, NULL, NULL, refmask, NULL);

    if (checkwave == NULL)
        fors_calib_exit("Wavelength calibration failure.");

    /*
     * Save check image to disk
     */

    header = cpl_propertylist_new();
    cpl_propertylist_update_double(header, "CRPIX1", 1.0);
    cpl_propertylist_update_double(header, "CRPIX2", 1.0);
    cpl_propertylist_update_double(header, "CRVAL1", 
                                   config.startwavelength+config.dispersion/2);
    cpl_propertylist_update_double(header, "CRVAL2", 1.0);
    /* cpl_propertylist_update_double(header, "CDELT1", config.dispersion);
    cpl_propertylist_update_double(header, "CDELT2", 1.0); */
    cpl_propertylist_update_double(header, "CD1_1", config.dispersion);
    cpl_propertylist_update_double(header, "CD1_2", 0.0);
    cpl_propertylist_update_double(header, "CD2_1", 0.0);
    cpl_propertylist_update_double(header, "CD2_2", 1.0);
    cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
    cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");

    fors_dfs_save_image(frameset, checkwave, spectra_detection_tag, header, 
                        parlist, recipe, ref_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_image_delete(checkwave); checkwave = NULL;
    cpl_propertylist_delete(header); header = NULL;

    cpl_msg_info(recipe, "Locate slits at reference wavelength on CCD...");
    slits = mos_locate_spectra(refmask);

    if (!slits) {
        cpl_msg_error(cpl_func, "Error found in %s: %s",
                      cpl_error_get_where(), cpl_error_get_message());
        fors_calib_exit("No slits could be detected!");
    }

    refimage = cpl_image_new_from_mask(refmask);
    cpl_mask_delete(refmask); refmask = NULL;

    save_header = dfs_load_header(frameset, arc_tag, 0);
    fors_dfs_save_image(frameset, refimage, slit_map_tag, save_header,
                        parlist, recipe, ref_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_propertylist_delete(save_header); save_header = NULL;

    cpl_image_delete(refimage); refimage = NULL;

    if (config.slit_ident) {

        /*
         * Attempt slit identification: this recipe may continue even
         * in case of failed identification (i.e., the position table is 
         * not produced, but an error is not set). In case of failure,
         * the spectra would be still extracted, even if they would not
         * be associated to slits on the mask.
         * 
         * The reason for making the slit identification an user option 
         * (via the parameter slit_ident) is to offer the possibility 
         * to avoid identifications that are only apparently successful, 
         * as it would happen in the case of an incorrect slit description 
         * in the data header.
         */

        cpl_msg_indent_less();
        cpl_msg_info(recipe, "Attempt slit identification (optional)...");
        cpl_msg_indent_more();

        positions = mos_identify_slits(slits, maskslits, NULL);

        if (positions) {
            cpl_table_delete(slits);
            slits = positions;

            /*
             * Eliminate slits which are _entirely_ outside the CCD
             */

            cpl_table_and_selected_double(slits, 
                                          "ybottom", CPL_GREATER_THAN, ny-1);
            cpl_table_or_selected_double(slits, 
                                          "ytop", CPL_LESS_THAN, 0);
            cpl_table_erase_selected(slits);

            nslits = cpl_table_get_nrow(slits);

            if (nslits == 0)
                fors_calib_exit("No slits found on the CCD");

            cpl_msg_info(recipe, "%d slits are entirely or partially "
                         "contained in CCD", nslits);

        }
        else {
            config.slit_ident = 0;
            cpl_table_delete(maskslits); maskslits = NULL;
            cpl_msg_info(recipe, "Global distortion model cannot be computed");
            if (cpl_error_get_code() != CPL_ERROR_NONE) {
                fors_calib_exit(NULL);
            }
        }
    }


    /*
     * Determination of spectral curvature
     */

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Determining spectral curvature...");
    cpl_msg_indent_more();

    cpl_msg_info(recipe, "Tracing master flat field spectra edges...");
    traces = mos_trace_flat(trace_flat, slits, ref_wave, 
                            config.startwavelength, config.endwavelength, 
                            config.dispersion);

    if (!traces)
        fors_calib_exit("Tracing failure");

    cpl_image_delete(added_flat); added_flat = NULL;

    cpl_msg_info(recipe, "Fitting flat field spectra edges...");
    polytraces = mos_poly_trace(slits, traces, config.cdegree);

    if (!polytraces)
        fors_calib_exit("Trace fitting failure");

    if (config.cmode) {
        cpl_msg_info(recipe, "Computing global spectral curvature model...");
        mos_global_trace(slits, polytraces, config.cmode);
    }

    fors_dfs_save_table(frameset, traces, curv_traces_tag, NULL, parlist,
                        recipe, ref_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_table_delete(traces); traces = NULL;

    coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
    spatial = mos_spatial_calibration(spectra, slits, polytraces, ref_wave, 
                  config.startwavelength, config.endwavelength, 
                  config.dispersion, 0, coordinate);

    if (!config.slit_ident) {
        cpl_image_delete(spectra); spectra = NULL;
    }


    /*
     * Final wavelength calibration of spectra having their curvature
     * removed
     */

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Perform final wavelength calibration...");
    cpl_msg_indent_more();

    nx = cpl_image_get_size_x(spatial);
    ny = cpl_image_get_size_y(spatial);

    idscoeff = cpl_table_new(ny);
    restable = cpl_table_new(nlines);
    rainbow = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
    residual = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
    fiterror = static_cast<double*>(cpl_calloc(ny, sizeof(double)));
    fitlines = static_cast<int*>(cpl_calloc(ny, sizeof(int)));

    //Table with positions of the detected lines used for wavelength calibration
    cpl_table * detected_lines = cpl_table_new(1);
    
    rectified = mos_wavelength_calibration_final(spatial, slits, lines, 
                     config.dispersion, config.peakdetection, config.wradius, 
                     config.wdegree, config.wreject, ref_wave, 
                     &config.startwavelength, &config.endwavelength, fitlines, 
                     fiterror, idscoeff, rainbow, residual, restable, 
                     detected_lines);

        fors_dfs_save_image(frameset, spatial, arc_rectified_tag, header,
                        parlist, recipe, ref_arc_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    if (rectified == NULL)
        fors_calib_exit("Wavelength calibration failure.");

    fors_dfs_save_table(frameset, restable, disp_residuals_table_tag, NULL,
                        parlist, recipe, ref_arc_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    fors_dfs_save_table(frameset, detected_lines, detected_lines_tag, NULL,
                        parlist, recipe, ref_arc_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_table_delete(restable); restable = NULL;
    cpl_table_delete(detected_lines); detected_lines = NULL;

    cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
    cpl_table_set_column_unit(idscoeff, "error", "pixel");
    cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;

    for (i = 0; i < ny; i++)
        if (!cpl_table_is_valid(idscoeff, "c0", i))
            cpl_table_set_invalid(idscoeff, "error", i);

    if (config.wmosmode > 0) {
        mos_interpolate_wavecalib_slit(idscoeff, slits, 1, config.wmosmode - 1);

        cpl_image_delete(rectified);

        rectified = mos_wavelength_calibration(spatial, ref_wave,
                         config.startwavelength, config.endwavelength,
                         config.dispersion, idscoeff, 0);
    }

    cpl_image_delete(spatial); spatial = NULL;

    delta = mos_map_pixel(idscoeff, ref_wave, config.startwavelength,
                          config.endwavelength, config.dispersion, 2);

    /* Get the mosca wave calib */
    mosca::wavelength_calibration wave_cal(idscoeff, ref_wave);
    
    /* Check that the wavelength solution is monotonically increasing */
    for(size_t spa_row = 0 ; spa_row < (size_t)ny; spa_row++)
        if(wave_cal.has_valid_cal((double)spa_row))
            if(!wave_cal.is_monotonical(spa_row, config.startwavelength,
                                        config.endwavelength, 
                                        config.dispersion))
            {
                std::stringstream error_msg;
                error_msg <<"The wavelength solution at row "<<spa_row<<
                     " does not increase monotonically, "
                     "which is physically impossible. Try with new parameters.";
                throw std::range_error(error_msg.str());
            }

    header = cpl_propertylist_new();
    cpl_propertylist_update_double(header, "CRPIX1", 1.0);
    cpl_propertylist_update_double(header, "CRPIX2", 1.0);
    cpl_propertylist_update_double(header, "CRVAL1",
                                   config.startwavelength+config.dispersion/2);
    cpl_propertylist_update_double(header, "CRVAL2", 1.0);
    /* cpl_propertylist_update_double(header, "CDELT1", config.dispersion);
    cpl_propertylist_update_double(header, "CDELT2", 1.0); */
    cpl_propertylist_update_double(header, "CD1_1", config.dispersion);
    cpl_propertylist_update_double(header, "CD1_2", 0.0);
    cpl_propertylist_update_double(header, "CD2_1", 0.0);
    cpl_propertylist_update_double(header, "CD2_2", 1.0);
    cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
    cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");

    fors_dfs_save_image(frameset, delta, delta_image_tag,
                        header, parlist, recipe, ref_arc_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_image_delete(delta); delta = NULL;
    cpl_propertylist_delete(header); header = NULL;

    mean_rms = mos_distortions_rms(rectified, lines, config.startwavelength, 
                                   config.dispersion, 6, 0);

    cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);

    mean_rms = cpl_table_get_column_mean(idscoeff, "error");
    mean_rms_err = cpl_table_get_column_stdev(idscoeff, "error");

    cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)", 
                 mean_rms, mean_rms * config.dispersion);

    restab = mos_resolution_table(rectified, config.startwavelength, 
                                  config.dispersion, 60000, lines);

    if (restab) {
        cpl_msg_info(recipe, "Mean spectral resolution: %.2f", 
                   cpl_table_get_column_mean(restab, "resolution"));
        cpl_msg_info(recipe, "Mean reference lines FWHM: %.2f +/- %.2f pixel",
             cpl_table_get_column_mean(restab, "fwhm") / config.dispersion,
             cpl_table_get_column_mean(restab, "fwhm_rms") / config.dispersion);

        qclist = cpl_propertylist_new();

        /*
         * QC1 parameters
         */
        keyname = "QC.DID";

        if (fors_header_write_string(qclist,
                keyname,
                "2.0",
                "QC1 dictionary")) {
            fors_calib_exit("Cannot write dictionary version "
                    "to QC log file");
        }

        if (mos)
            keyname = "QC.MOS.RESOLUTION";
        else
            keyname = "QC.MXU.RESOLUTION";

        if (fors_header_write_double(qclist, 
                cpl_table_get_column_mean(restab,
                        "resolution"),
                        keyname,
                        "Angstrom",
                        "Mean spectral resolution")) {
            fors_calib_exit("Cannot write mean spectral resolution to QC "
                    "log file");
        }

        if (mos)
            keyname = "QC.MOS.RESOLUTION.RMS";
        else
            keyname = "QC.MXU.RESOLUTION.RMS";

        if (fors_header_write_double(qclist, 
                cpl_table_get_column_stdev(restab, 
                        "resolution"),
                        keyname,
                        "Angstrom", 
                        "Scatter of spectral resolution")) {
            fors_calib_exit("Cannot write spectral resolution scatter "
                    "to QC log file");
        }

        if (mos)
            keyname = "QC.MOS.RESOLUTION.NWAVE";
        else
            keyname = "QC.MXU.RESOLUTION.NWAVE";

        if (fors_header_write_int(qclist, cpl_table_get_nrow(restab) -
                cpl_table_count_invalid(restab, 
                        "resolution"),
                        keyname,
                        NULL,
                        "Number of examined wavelengths "
                        "for resolution computation")) {
            fors_calib_exit("Cannot write number of lines used in "
                    "spectral resolution computation "
                    "to QC log file");
        }

        if (mos)
            keyname = "QC.MOS.RESOLUTION.MEANRMS";
        else
            keyname = "QC.MXU.RESOLUTION.MEANRMS";

        if (fors_header_write_double(qclist,
                cpl_table_get_column_mean(restab,
                        "resolution_rms"),
                        keyname, NULL,
                        "Mean error on spectral "
                        "resolution computation")) {
            fors_calib_exit("Cannot write mean error in "
                    "spectral resolution computation "
                    "to QC log file");
        }

        if (mos)
            keyname = "QC.MOS.RESOLUTION.NLINES";
        else
            keyname = "QC.MXU.RESOLUTION.NLINES";

        if (fors_header_write_int(qclist,
                cpl_table_get_column_mean(restab, "nlines") *
                cpl_table_get_nrow(restab),
                keyname, NULL,
                "Number of lines for spectral "
                "resolution computation")) {
            fors_calib_exit("Cannot write number of examined "
                    "wavelengths in spectral resolution computation "
                    "to QC log file");
        }

        fors_dfs_save_table(frameset, restab, spectral_resolution_tag, qclist,
                            parlist, recipe, ref_arc_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_calib_exit(NULL);

        cpl_table_delete(restab); restab = NULL;
        cpl_propertylist_delete(qclist); qclist = NULL;

    }
    else
        fors_calib_exit("Cannot compute the spectral resolution table");

    cpl_vector_delete(lines); lines = NULL;

    fors_dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL,
                        parlist, recipe, ref_arc_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);
//%%%

//%%%

    /*
     * Global distortion models
     */

    if (config.slit_ident) {

        cpl_msg_info(recipe, "Computing global distortions model");
        global = mos_global_distortion(slits, maskslits, idscoeff, 
                                       polytraces, ref_wave);

        if (global && 0) {
            cpl_table *stest;
            cpl_table *ctest;
            cpl_table *dtest;
            cpl_image *itest;

            stest = mos_build_slit_location(global, maskslits, ccd_ysize);

            ctest = mos_build_curv_coeff(global, maskslits, stest);
            fors_dfs_save_table(frameset, ctest, "CURVS", NULL,
                                parlist, recipe, ref_flat_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_calib_exit(NULL);

            itest = mos_spatial_calibration(spectra, stest, ctest, 
                                            ref_wave, config.startwavelength, 
                                            config.endwavelength, 
                                            config.dispersion, 0, NULL);
            cpl_table_delete(ctest); ctest = NULL;
            cpl_image_delete(itest); itest = NULL;
            fors_dfs_save_table(frameset, stest, "SLITS", NULL,
                                parlist, recipe, ref_flat_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_calib_exit(NULL);

            dtest = mos_build_disp_coeff(global, stest);
            fors_dfs_save_table(frameset, dtest, "DISPS", NULL,
                                parlist, recipe, ref_flat_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_calib_exit(NULL);

            cpl_table_delete(dtest); dtest = NULL;
            cpl_table_delete(stest); stest = NULL;
        }

        if (global) {
            fors_dfs_save_table(frameset, global, global_distortion_tag, NULL,
                    parlist, recipe, ref_arc_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_calib_exit(NULL);

            cpl_table_delete(global); global = NULL;
        }

        cpl_image_delete(spectra); spectra = NULL;
        cpl_table_delete(maskslits); maskslits = NULL;
    }

    /* Create header for wavelength calibrated images */
    header = cpl_propertylist_new();
    cpl_propertylist_update_double(header, "CRPIX1", 1.0);
    cpl_propertylist_update_double(header, "CRPIX2", 1.0);
    cpl_propertylist_update_double(header, "CRVAL1", 
                                   config.startwavelength+config.dispersion/2);
    cpl_propertylist_update_double(header, "CRVAL2", 1.0);
    /* cpl_propertylist_update_double(header, "CDELT1", config.dispersion);
    cpl_propertylist_update_double(header, "CDELT2", 1.0); */
    cpl_propertylist_update_double(header, "CD1_1", config.dispersion);
    cpl_propertylist_update_double(header, "CD1_2", 0.0);
    cpl_propertylist_update_double(header, "CD2_1", 0.0);
    cpl_propertylist_update_double(header, "CD2_2", 1.0);
    cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
    cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
    cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);

    /* Rectify the arc masks */
    cpl_mask * nonlinear_arc_mask_rect_mapped =
            fors_calib_mask_rect_mapped(nonlinear_arc_mask, slits, idscoeff,
                                        polytraces,ref_wave, config);

    cpl_mask * saturated_arc_mask_rect_mapped =
            fors_calib_mask_rect_mapped(saturated_arc_mask, slits, idscoeff,
                                        polytraces,ref_wave, config);

    cpl_image * combined_arc_mask = 
            fors_bpm_create_combined_bpm(nonlinear_arc_mask_rect_mapped,
                                         saturated_arc_mask_rect_mapped);
    
    /* Save the arc */
    fors_dfs_save_image_mask(frameset, rectified, combined_arc_mask,
                             reduced_lamp_tag, header, 
                             parlist, recipe, ref_arc_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_mask_delete(nonlinear_arc_mask);
    cpl_mask_delete(nonlinear_arc_mask_rect_mapped);
    cpl_mask_delete(saturated_arc_mask);
    cpl_mask_delete(saturated_arc_mask_rect_mapped);
    cpl_image_delete(combined_arc_mask);
    cpl_image_delete(rectified); rectified = NULL;
    cpl_propertylist_delete(header); header = NULL;
    
    save_header = dfs_load_header(frameset, arc_tag, 0);

    cpl_propertylist_update_double(save_header, "CRPIX2", 1.0);
    cpl_propertylist_update_double(save_header, "CRVAL2", 1.0);
    /* cpl_propertylist_update_double(save_header, "CDELT2", 1.0); */
    cpl_propertylist_update_double(save_header, "CD1_1", 1.0);
    cpl_propertylist_update_double(save_header, "CD1_2", 0.0);
    cpl_propertylist_update_double(save_header, "CD2_1", 0.0);
    cpl_propertylist_update_double(save_header, "CD2_2", 1.0);
    cpl_propertylist_update_string(save_header, "CTYPE1", "LINEAR");
    cpl_propertylist_update_string(save_header, "CTYPE2", "PIXEL");

    fors_dfs_save_image(frameset, residual, disp_residuals_tag, save_header,
                        parlist, recipe, ref_arc_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_image_delete(residual); residual = NULL;
    cpl_propertylist_delete(save_header); save_header = NULL;

    wavemap = mos_map_wavelengths(coordinate, rainbow, slits, polytraces, 
                                  ref_wave, config.startwavelength, 
                                  config.endwavelength, config.dispersion);

    cpl_image_delete(rainbow); rainbow = NULL;

    save_header = dfs_load_header(frameset, arc_tag, 0);

    /*
     * QC1 parameters
     */
    keyname = "QC.DID";

    if (fors_header_write_string(save_header,
            keyname,
            "2.0",
            "QC1 dictionary")) {
        fors_calib_exit("Cannot write dictionary version "
                "to QC log file");
    }

    if (fors_header_write_double(save_header,
            mean_rms,
            "QC.WAVE.ACCURACY",
            "pixel",
            "Mean accuracy of wavecalib model")) {
        fors_calib_exit("Cannot write mean wavelength calibration "
                "accuracy to QC log file");
    }


    if (fors_header_write_double(save_header,
            mean_rms_err,
            "QC.WAVE.ACCURACY.ERROR",
            "pixel",
            "Error on accuracy of wavecalib model")) {
        fors_calib_exit("Cannot write error on wavelength calibration "
                "accuracy to QC log file");
    }

    fors_dfs_save_image(frameset, wavemap, wavelength_map_tag, save_header,
            parlist, recipe, ref_arc_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_image_delete(wavemap); wavemap = NULL;

    cpl_propertylist_erase_regexp(save_header, "^ESO QC ", 0);

    fors_dfs_save_image(frameset, coordinate, spatial_map_tag, save_header,
                        parlist, recipe, ref_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_propertylist_delete(save_header); save_header = NULL;

    header = NULL;    /* To be really, really, REALLY sure... */

    /*
     * QC1 parameters
     */
    double maxpos, maxneg, maxcurve, maxslope;

    header = dfs_load_header(frameset, arc_tag, 0);

    keyname = "QC.DID";

    if (fors_header_write_string(header,
            keyname,
            "2.0",
            "QC1 dictionary")) {
        fors_calib_exit("Cannot write dictionary version "
                "to QC log file");
    }

    maxpos = fabs(cpl_table_get_column_max(polytraces, "c2"));
    maxneg = fabs(cpl_table_get_column_min(polytraces, "c2"));
    maxcurve = maxpos > maxneg ? maxpos : maxneg;
    if (fors_header_write_double(header,
            maxcurve,
            "QC.TRACE.MAX.CURVATURE",
            "Y pixel / X pixel ^2",
            "Max observed curvature in "
            "spectral tracing")) {
        fors_calib_exit("Cannot write max observed curvature in spectral "
                "tracing to QC log file");
    }

    maxpos = fabs(cpl_table_get_column_max(polytraces, "c1"));
    maxneg = fabs(cpl_table_get_column_min(polytraces, "c1"));
    maxslope = maxpos > maxneg ? maxpos : maxneg;
    if (fors_header_write_double(header,
            maxslope,
            "QC.TRACE.MAX.SLOPE",
            "Y pixel / X pixel",
            "Max observed slope in spectral tracing")) {
        fors_calib_exit("Cannot write max observed slope in spectral "
                "tracing to QC log file");
    }
    
    /* Saving slits and polytraces */
    
    fors_dfs_save_table(frameset, polytraces, curv_coeff_tag, header,
                        parlist, recipe, ref_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);

    cpl_propertylist_delete(header); header = NULL;

    fors_dfs_save_table(frameset, slits, slit_location_tag, NULL,
                        parlist, recipe, ref_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit(NULL);
    
    /* Read grism configuration */
    //TODO: Add the waveref to the grism_tables
    cpl_frameset * grism_frame = fors_frameset_extract(frameset, "GRISM_TABLE");
    std::auto_ptr<mosca::grism_config> grism_cfg = 
            fors_grism_config_from_frame(cpl_frameset_get_position(grism_frame, 0), ref_wave);
    cpl_frameset_delete(grism_frame);

    /* Get the detected slit locations */
    fors::detected_slits det_slits = 
        fors::detected_slits_from_tables(slits, polytraces, size_spec);
    
    /* Get the calibrated slits */
    fors::calibrated_slits calib_slits(det_slits, wave_cal, *grism_cfg,
                                       ccd_xsize, ccd_ysize); 
    for(std::vector<mosca::calibrated_slit>::const_iterator 
            slit_it = calib_slits.begin();
        slit_it != calib_slits.end() ; slit_it++)
    {
        if(!slit_it->has_valid_wavecal())
            cpl_msg_warning(cpl_func, "Slit %d does not contain valid "
                    "wavelength calibration. Skipping it for master flat", 
                    slit_it->slit_id());
    }
    
    /* Compute master flat.
     * TODO: master flat has already been computed above using the old method
     * Here we use the new method and is the one saved. The other is not yet
     * deleted in case it is used for something else.
     */
    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Perform flat field combination...");
    
    cpl_image_delete(master_flat);
    std::auto_ptr<mosca::image> master_flat_d;
    std::vector<std::vector<double> > slit_sat_ratio;
    std::vector<std::vector<int> > slit_sat_count;
    cpl_mask ** nonlinear_flat_masks;
    cpl_mask ** saturated_flat_masks;
    std::auto_ptr<mosca::fiera_config> ccd_config;
    master_flat_d = fors_calib_flat_mos_create_master_flat(calib_slits,
            wave_cal, *grism_cfg, master_bias, bias_frame,
            config, frameset, flat_tag, 
            config.nonlinear_level, config.max_nonlinear_ratio, 
            slit_sat_ratio, slit_sat_count,
            nonlinear_flat_masks, 
            saturated_flat_masks, ccd_config);
    if(master_flat_d.get() == 0)
        fors_calib_exit("Cannot combine flat frames");

    /* Create flat bad pixel mask */
    cpl_image * flat_mask = fors_bpm_create_combined_bpm(nonlinear_flat_masks,
                                                         saturated_flat_masks,
                                                         nflats);
    
    /*
     * Flat field normalisation is done directly on the master flat
     * field (without spatial rectification first). The spectral
     * curvature model may be provided in input, in future releases.
     */
    cpl_msg_info(recipe, "Performing flat field normalisation");
    std::auto_ptr<mosca::image> norm_flat(NULL);
    cpl_image * wave_profiles;
    std::vector<float> sed_norm;
    if(fors_calib_flat_mos_normalise(master_flat_d, wave_cal, slits, 
            polytraces,  coordinate, config, 
            norm_flat, &wave_profiles, sed_norm) != 0)
        fors_calib_exit("Cannot normalise flat");

    cpl_msg_info(recipe, "Performing flat field distortion correction");
    if(fors_calib_flat_mos_rect_mapped(master_flat_d, norm_flat, 
            slits, idscoeff, polytraces, ref_wave, config, 
            mapped_flat, mapped_nflat) != 0)
        fors_calib_exit("Cannot correct flat field from distortion");

     /* Getting the normalisation factors used in the SED */
     cpl_propertylist * sed_header = cpl_propertylist_new();
     for(size_t ised = 0 ; ised < sed_norm.size();ised++)
     {
         std::ostringstream norm_key;
         norm_key<< "ESO QC FLAT SED"<<ised+1<<" NORM ";
         cpl_propertylist_append_float(sed_header, norm_key.str().c_str(),
                                       sed_norm[ised]);
     }

    /* Saving all flats */
    cpl_msg_info(recipe, "Saving flats");
    if(fors_calib_flats_save(master_flat_d, flat_mask,
            norm_flat, mapped_flat, mapped_nflat, det_slits,
            slit_sat_ratio, slit_sat_count, config,
            frameset, flat_tag, master_screen_flat_tag, master_norm_flat_tag,
            mapped_screen_flat_tag, mapped_norm_flat_tag,
            parlist, ref_flat_frame, *ccd_config) != 0)
        fors_calib_exit("Cannot save flats");

    /* Saving slit dispersion profiles */
    std::ostringstream prof_filename_oss;
    prof_filename_oss << flat_disp_profile_tag << ".fits";
    std::string prof_filename = prof_filename_oss.str();
    std::transform(prof_filename.begin(), prof_filename.end(), prof_filename.begin(), ::tolower);
    cpl_propertylist_append_string(sed_header, CPL_DFS_PRO_TYPE, "REDUCED");
    cpl_propertylist_append_string(sed_header, CPL_DFS_PRO_CATG, flat_disp_profile_tag);
    cpl_dfs_save_image(frameset, NULL, parlist, frameset, 
                       ref_flat_frame, wave_profiles, CPL_BPP_IEEE_FLOAT,
                       recipe, sed_header, NULL, PACKAGE "/" PACKAGE_VERSION,  
                       prof_filename.c_str());
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_calib_exit("Cannot save flats");
    cpl_propertylist_delete(sed_header); sed_header = NULL;

    cpl_table_delete(polytraces); polytraces = NULL;
    cpl_table_delete(slits); slits = NULL;
    cpl_table_delete(idscoeff); idscoeff = NULL;
    cpl_image_delete(coordinate); coordinate = NULL;
    fors_image_delete(&master_bias); master_bias = NULL;
    cpl_image_delete(mapped_flat);
    cpl_image_delete(mapped_nflat);
    cpl_image_delete(wave_profiles);
    cpl_image_delete(flat_mask);
    for (size_t i_flat = 0; i_flat < nflats; i_flat++)
    {
        cpl_mask_delete(nonlinear_flat_masks[i_flat]);
        cpl_mask_delete(saturated_flat_masks[i_flat]);
    }
    cpl_free(nonlinear_flat_masks);
    cpl_free(saturated_flat_masks);

    if (cpl_error_get_code()) {
        cpl_msg_error(cpl_func, "Error found in %s: %s",
                      cpl_error_get_where(), cpl_error_get_message());
        fors_calib_exit(NULL);
    }

    return 0;
}

int fors_calib_retrieve_input_param(cpl_parameterlist * parlist, 
                                     cpl_frameset * frameset,
                                     fors_calib_config * config)
{
    const char *recipe = "fors_calib";

    cpl_table        *grism_table  = NULL;

    cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
    cpl_msg_indent_more();
    
    grism_table = dfs_load_table(frameset, "GRISM_TABLE", 1);

    config->dispersion = dfs_get_parameter_double(parlist, 
                    "fors.fors_calib.dispersion", grism_table);


    config->peakdetection = dfs_get_parameter_double(parlist, 
                    "fors.fors_calib.peakdetection", grism_table);

    config->wdegree = dfs_get_parameter_int(parlist, 
                    "fors.fors_calib.wdegree", grism_table);

    config->wradius = dfs_get_parameter_int(parlist, "fors.fors_calib.wradius", NULL);

    config->wreject = dfs_get_parameter_double(parlist, 
                                       "fors.fors_calib.wreject", NULL);

    config->wmode = dfs_get_parameter_int(parlist, "fors.fors_calib.wmode", NULL);

    config->wmosmode = dfs_get_parameter_int(parlist,
                                     "fors.fors_calib.wmosmode", NULL);

    config->cdegree = dfs_get_parameter_int(parlist, "fors.fors_calib.cdegree", 
                                    grism_table);

    config->cmode = dfs_get_parameter_int(parlist, "fors.fors_calib.cmode", NULL);

    config->startwavelength = dfs_get_parameter_double(parlist, 
                    "fors.fors_calib.startwavelength", grism_table);

    config->endwavelength = dfs_get_parameter_double(parlist, 
                    "fors.fors_calib.endwavelength", grism_table);

    config->slit_ident = dfs_get_parameter_bool(parlist, 
                    "fors.fors_calib.slit_ident", NULL);

    config->stack_method = dfs_get_parameter_string(parlist, 
                                        "fors.fors_calib.stack_method", NULL);

    if (strcmp(config->stack_method, "ksigma") == 0) {
        std::string ksigma = dfs_get_parameter_string(parlist,
                                         "fors.fors_calib.ksigma", NULL);
        std::string::size_type comma = ksigma.find(',');
        std::istringstream klow(ksigma.substr(0, comma));
        std::istringstream khigh(ksigma.substr(comma+1));
        bool invalid_klow = !(klow >> config->klow) || !(klow.eof()  || (klow >> std::ws && klow.eof())); 
        bool invalid_khigh = !(khigh >> config->khigh) || !(khigh.eof()  || (khigh >> std::ws && khigh.eof()));
        if (comma == std::string::npos || invalid_klow || invalid_khigh)
            throw std::invalid_argument("ksigma must contain two "
                                        "comma-separated numbers");

        config->kiter = dfs_get_parameter_int(parlist, 
                                         "fors.fors_calib.kiter", NULL);
    }

    config->spa_polyorder = dfs_get_parameter_int(parlist, "fors.fors_calib.s_degree", NULL);

    config->disp_nknots = dfs_get_parameter_int(parlist, "fors.fors_calib.d_nknots", NULL);

    config->sradius = dfs_get_parameter_int(parlist, "fors.fors_calib.sradius", NULL);

    config->dradius = dfs_get_parameter_int(parlist, "fors.fors_calib.dradius", NULL);

    config->fit_threshold = dfs_get_parameter_double(parlist, 
            "fors.fors_calib.fit_threshold", NULL);
    
    config->ignore_lines= dfs_get_parameter_string(parlist, 
            "fors.fors_calib.ignore_lines", NULL);

    config->used_linesets= dfs_get_parameter_string(parlist, 
            "fors.fors_calib.used_linesets", NULL);

    config->nonlinear_level = dfs_get_parameter_double(parlist, 
            "fors.fors_calib.nonlinear_level", NULL);

    config->max_nonlinear_ratio = dfs_get_parameter_double(parlist, 
            "fors.fors_calib.max_nonlinear_ratio", NULL);
    
    cpl_table_delete(grism_table); grism_table = NULL;

    return 0; 
}

//functor to search substrings with find_if. C++11 lambdas would be shorter
struct contains_subtring : std::unary_function<std::string ,bool>
{
    contains_subtring(std::string& substr) : m_substr(substr) {};
    bool operator()(const std::string& string) 
    {
        return (string.find(m_substr) != std::string::npos); 
    }
    std::string m_substr;
};

//Get the list of lines after proper filtering
cpl_vector * fors_calib_get_reference_lines(cpl_frameset * frameset, 
                                            const char * arctag,
                                            const char * ignore_lines,
                                            const char * used_linesets)
{
    cpl_table        *wavelengths  = NULL;
    cpl_propertylist *archeader  = NULL;
    cpl_size          nlines_all;
    cpl_size          n_selected = 0;
    cpl_size          i;
    cpl_vector       *lines;
    double            lambda;
    int               null;
    const char *      wcolumn = "WLEN";
    const char *      ioncolumn = "CHEMICAL_ION";
    const char *      linesetcolumn = "LINE_SET";

    /*
     * Read the wavelengths table 
     */
    wavelengths = dfs_load_table(frameset, "MASTER_LINECAT", 1);
    archeader   = dfs_load_header(frameset, arctag, 0);

    if (wavelengths == NULL)
    {
        cpl_msg_error(cpl_func, "Cannot load line catalog");
        return NULL;
    }


    nlines_all = cpl_table_get_nrow(wavelengths);

    if (nlines_all == 0)
    {
        cpl_msg_error(cpl_func, "Empty input line catalog");
        cpl_table_delete(wavelengths);
        return NULL;
    }

    if (cpl_table_has_column(wavelengths, wcolumn) != 1 ||
        cpl_table_has_column(wavelengths, ioncolumn) != 1 || 
        cpl_table_has_column(wavelengths, linesetcolumn) != 1) 
    {
        cpl_msg_error(cpl_func, "Missing columns %s %s %s in input line catalog",
                      wcolumn, ioncolumn, linesetcolumn);
        cpl_table_delete(wavelengths);
        return NULL;
    }

    /*
     * Deselect lines which are not present in the lamps
     */
    std::vector<std::string> observed_lamps;
    for(size_t ikeylamp = 1; ikeylamp < 9; ++ikeylamp)
    {
        std::ostringstream oss;
        oss<<"ESO INS LAMP"<<ikeylamp<<" NAME";
        if(cpl_propertylist_has(archeader, oss.str().c_str()))
            observed_lamps.push_back
                (cpl_propertylist_get_string(archeader, oss.str().c_str()));
    }
    //Check line by line if the chemical element was part of any of the lamps
    for (cpl_size iline = 0; iline < nlines_all; iline++)
    {
        if(cpl_table_get_string(wavelengths, ioncolumn, iline) == NULL)
        {
            cpl_table_unselect_row(wavelengths, iline);
            break;
        }
        std::string chem_elem = 
            cpl_table_get_string(wavelengths, ioncolumn, iline);
        //Get up to the first space (i. e. from 'He I', get just 'He')
        std::string::size_type space = chem_elem.find(' ');
        if (space != std::string::npos)
            chem_elem = chem_elem.substr(0, space);
        //Check if the chemical element is contained in any of the lamps names
        if (std::find_if(observed_lamps.begin(), observed_lamps.end(), 
            contains_subtring(chem_elem)) == observed_lamps.end())
            cpl_table_unselect_row(wavelengths, iline);
    }

    /*
     * Deselect lines which are not in the lineset 
     */
    //Parse command line option used_linesets
    std::string used_linesets_str(used_linesets);
    std::vector<std::string> linesets;
    while(used_linesets_str.length() > 0)
    {
        //Parsing used_linesets (values are separated by comma)
        std::string::size_type found = used_linesets_str.find(',');
        std::string lineset;
        if(found != std::string::npos)
        {
            lineset = used_linesets_str.substr(0, found);
            used_linesets_str = used_linesets_str.substr(found+1);
        }
        else
        {
            lineset = used_linesets_str;
            used_linesets_str = "";
        }
        linesets.push_back(lineset);
    }
    for(cpl_size iline = 0; iline < nlines_all; iline++)
    {
        if(cpl_table_get_string(wavelengths, linesetcolumn, iline) == NULL)
        {
            cpl_table_unselect_row(wavelengths, iline);
            break;
        }
        std::string table_lineset = 
                cpl_table_get_string(wavelengths, linesetcolumn, iline);
        if(std::find(linesets.begin(), linesets.end(), table_lineset) == linesets.end())
            cpl_table_unselect_row(wavelengths, iline);
    }

    /*
     * Deselect lines which are present in ignore_lines 
     */
    std::string ignore_lines_str(ignore_lines);
    while(ignore_lines_str.length() > 0)
    {
        //Parsing ignore_lines (values are separated by comma)
        std::string::size_type found = ignore_lines_str.find(',');
        std::string lambda_str;
        if(found != std::string::npos)
        {
            lambda_str = ignore_lines_str.substr(0, found);
            ignore_lines_str = ignore_lines_str.substr(found+1);
        }
        else
        {
            lambda_str = ignore_lines_str;
            ignore_lines_str = "";
        }
        std::istringstream iss(lambda_str);
        if ( !(iss >> lambda) || !(iss.eof() || (iss >> std::ws && iss.eof())) )
        {
            cpl_msg_error(cpl_func, "Cannot interpret number in ignored_lines");
            cpl_table_delete(wavelengths);
            return NULL;
        }

        //Search for closest line in catalog. The line is unselected but
        //it will be checked again against the next ignored line. In this way,
        //if a value appears many times in the ignored_lines, only one line
        //will be removed
        cpl_size i_ignore = 0;
        double min_lambda_dif = 
             std::fabs(lambda - cpl_table_get(wavelengths, wcolumn, 0, &null));
        for (i = 1; i < nlines_all; i++)
        {
            double lambda_dif = 
              std::fabs(lambda - cpl_table_get(wavelengths, wcolumn, i, &null));
            if(lambda_dif < min_lambda_dif)
            {
                min_lambda_dif = lambda_dif;
                i_ignore = i;
            }
         }
        cpl_table_unselect_row(wavelengths, i_ignore);
    } 
    
    /* Create the final list of reference lines */
    n_selected = cpl_table_count_selected(wavelengths);
    lines = cpl_vector_new(n_selected);
    cpl_size i_line = 0;
    for (i = 0; i < nlines_all; i++)
    {
        lambda = cpl_table_get(wavelengths, wcolumn, i, &null);
        if(cpl_table_is_selected(wavelengths, i))
        {
            cpl_vector_set(lines, i_line, lambda);
            i_line++;
        }
    }

    cpl_table_delete(wavelengths);
    cpl_propertylist_delete(archeader);

    return lines;
}

int fors_calib_flat_mos_normalise
(std::auto_ptr<mosca::image>& master_flat_d,
 const mosca::wavelength_calibration& wave_cal,
 cpl_table * slits, cpl_table * polytraces, cpl_image * coordinate, 
 struct fors_calib_config& config,
 std::auto_ptr<mosca::image>& norm_flat, 
 cpl_image ** wave_profiles,
 std::vector<float>& sed_norm)
{
    cpl_msg_indent_more();

    norm_flat.reset(new mosca::image(cpl_image_cast(master_flat_d->get_cpl_image(), 
                                                    CPL_TYPE_FLOAT),
                                     cpl_image_cast(master_flat_d->get_cpl_image_err(), 
                                                    CPL_TYPE_FLOAT), true));
    
    /* Flat normalisation */
    fors::flat_normaliser normaliser;
    normaliser.mos_normalise(*norm_flat, wave_cal,
                             coordinate, slits, polytraces,
                             config.startwavelength, config.endwavelength,
                             config.dispersion, config.sradius, config.dradius,
                             config.spa_polyorder, config.disp_nknots, config.fit_threshold);

    /* Get the spectral shape of the slits and save it */
    *wave_profiles = normaliser.get_wave_profiles_im();

    /* Get the normalisation factors used */
    sed_norm = normaliser.get_wave_profiles_norm();

    cpl_msg_indent_less();

    return 0;
}

int fors_calib_flat_mos_rect_mapped
(std::auto_ptr<mosca::image>& master_flat_d, 
 std::auto_ptr<mosca::image>& norm_flat,
 cpl_table * slits,
 cpl_table *idscoeff, cpl_table * polytraces,
 double reference, struct fors_calib_config& config,
 cpl_image *& mapped_flat, 
 cpl_image *& mapped_nflat)
{
    cpl_image * rect_flat;
    cpl_image * rect_nflat = NULL;
    
    cpl_msg_indent_more();

    /* mos_spatial-calibration cannot accept doubles 
     * At the end I changed the master flat calibration for float, but the 
     * output is still double (see TODO comment on flat_combine) */
    cpl_image * master_flat = cpl_image_cast(master_flat_d->get_cpl_image(),
                                             CPL_TYPE_FLOAT);
    
    /* Flat spatial distortion correction */ 
    rect_flat = mos_spatial_calibration(master_flat, slits, polytraces, 
                                        reference, config.startwavelength, 
                                        config.endwavelength, config.dispersion,
                                        0, NULL);
    if(norm_flat.get() != NULL)
    {
        cpl_image * norm_flat_f = cpl_image_cast(norm_flat->get_cpl_image(),
                                                 CPL_TYPE_FLOAT);
        rect_nflat = mos_spatial_calibration(norm_flat_f, slits, polytraces, 
                                             reference, config.startwavelength, 
                                             config.endwavelength, 
                                             config.dispersion, 0, NULL);
        cpl_image_delete(norm_flat_f);
    }

    /* Flat wavelength calibration */
    mapped_flat = mos_wavelength_calibration(rect_flat, reference,
                                             config.startwavelength, 
                                             config.endwavelength,
                                             config.dispersion, idscoeff, 0);

    if(norm_flat.get() != NULL)
        mapped_nflat = mos_wavelength_calibration(rect_nflat, reference,
                                                  config.startwavelength, 
                                                  config.endwavelength, 
                                                  config.dispersion, idscoeff, 
                                                  0);

    cpl_image_delete(master_flat);
    cpl_image_delete(rect_flat);
    if(norm_flat.get() != NULL)
        cpl_image_delete(rect_nflat);
    cpl_msg_indent_less();

    return 0;
}

cpl_mask * fors_calib_mask_rect_mapped
(cpl_mask * mask,
 cpl_table * slits, cpl_table *idscoeff, cpl_table * polytraces,
 double reference, struct fors_calib_config& config)
{
    cpl_image * target_image = cpl_image_new_from_mask(mask);
        
    /* Spatial distortion correction */ 
    cpl_image * rect_image = mos_spatial_calibration(target_image, slits, 
                                                     polytraces, 
                                                     reference, 
                                                     config.startwavelength, 
                                                     config.endwavelength, 
                                                     config.dispersion,
                                                     0, NULL);

    /* Wavelength calibration */
    cpl_image * mapped_image = mos_wavelength_calibration(rect_image, reference,
                                                          config.startwavelength, 
                                                          config.endwavelength,
                                                          config.dispersion, 
                                                          idscoeff, 0);


    cpl_mask * rect_mapped_mask = 
            cpl_mask_threshold_image_create(mapped_image,
                            0, std::numeric_limits<double>::max());

    cpl_image_delete(rect_image);
    cpl_image_delete(mapped_image);
    cpl_image_delete(target_image);
    
    return rect_mapped_mask;
}

int fors_calib_flats_save
(std::auto_ptr<mosca::image>& master_flat_d, 
 cpl_image * flat_mask,
 std::auto_ptr<mosca::image>& norm_flat,
 cpl_image * mapped_flat,  cpl_image * mapped_nflat,
 const fors::detected_slits detected_slits,
 const std::vector<std::vector<double> >& slit_sat_ratio,
 const std::vector<std::vector<int> >& slit_sat_count,
 struct fors_calib_config& config,
 cpl_frameset * frameset, const char * flat_tag, 
 const char * master_screen_flat_tag, const char * master_norm_flat_tag, 
 const char * mapped_screen_flat_tag, const char * mapped_norm_flat_tag, 
 cpl_parameterlist * parlist, const cpl_frame * ref_flat_frame, 
 const mosca::ccd_config& ccd_config)
{
    cpl_propertylist * save_header;
    cpl_propertylist * wave_header;
    const char *recipe_name = "fors_calib";
    
    cpl_msg_indent_more();

    size_t nflats = cpl_frameset_get_size(frameset);
    save_header = dfs_load_header(frameset, flat_tag, 0);
    cpl_propertylist_update_int(save_header, "ESO PRO DATANCOM", nflats);

    /* Computing QC for saturation */
    fors_calib_qc_saturation(save_header, detected_slits,
                             slit_sat_ratio, slit_sat_count);
    
    /* Adding the trimming keywords */
    fors_trimm_fill_info(save_header, ccd_config);

    /* Saving regular flat */
    fors_image * fors_master_flat = fors_image_new(
            cpl_image_duplicate(master_flat_d->get_cpl_image()),
            cpl_image_power_create(master_flat_d->get_cpl_image_err(), 2.));
    fors_dfs_save_image_err_mask(frameset, fors_master_flat, flat_mask,
                                 master_screen_flat_tag, save_header, parlist,
                                 recipe_name, ref_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
    {
        cpl_propertylist_delete(save_header);
        return -1;
    }

    /* Saving normalised flats */
    if(norm_flat.get() != NULL)
    {
        fors_image * fors_norm_flat = fors_image_new(
                cpl_image_duplicate(norm_flat->get_cpl_image()),
                cpl_image_power_create(norm_flat->get_cpl_image_err(), 2.));
        fors_dfs_save_image_err_mask(frameset, fors_norm_flat, flat_mask,
                                     master_norm_flat_tag, save_header, parlist,
                                     recipe_name, ref_flat_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
        {
            cpl_propertylist_delete(save_header);
            return -1;
        }
        fors_image_delete(&fors_norm_flat);
    }

    /* Create header for wavelength calibrated images */
    wave_header = cpl_propertylist_new();
    cpl_propertylist_update_double(wave_header, "CRPIX1", 1.0);
    cpl_propertylist_update_double(wave_header, "CRPIX2", 1.0);
    cpl_propertylist_update_double(wave_header, "CRVAL1", 
                                   config.startwavelength + config.dispersion/2);
    cpl_propertylist_update_double(wave_header, "CRVAL2", 1.0);
    /* cpl_propertylist_update_double(header, "CDELT1", config.dispersion);
    cpl_propertylist_update_double(header, "CDELT2", 1.0); */
    cpl_propertylist_update_double(wave_header, "CD1_1", config.dispersion);
    cpl_propertylist_update_double(wave_header, "CD1_2", 0.0);
    cpl_propertylist_update_double(wave_header, "CD2_1", 0.0);
    cpl_propertylist_update_double(wave_header, "CD2_2", 1.0);
    cpl_propertylist_update_string(wave_header, "CTYPE1", "LINEAR");
    cpl_propertylist_update_string(wave_header, "CTYPE2", "PIXEL");

    cpl_propertylist_update_int(wave_header, "ESO PRO DATANCOM", nflats);

    /* Saving mapped flat */
    fors_dfs_save_image(frameset, mapped_flat, mapped_screen_flat_tag,
                        wave_header, parlist, recipe_name, ref_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
    {
        cpl_propertylist_delete(wave_header);
        cpl_propertylist_delete(save_header);
        return -1;
    }

    /* Saving normalised mapped flat */
    if(mapped_nflat != NULL)
    {
        fors_dfs_save_image(frameset, mapped_nflat, mapped_norm_flat_tag, 
                wave_header, parlist, recipe_name, ref_flat_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
        {
            cpl_propertylist_delete(wave_header);
            cpl_propertylist_delete(save_header);
            return -1;
        }
    }

    cpl_propertylist_delete(wave_header);
    cpl_propertylist_delete(save_header);
    fors_image_delete(&fors_master_flat);

    cpl_msg_indent_less();

    return 0;
}

std::auto_ptr<mosca::image> fors_calib_flat_mos_create_master_flat
(fors::calibrated_slits& calibrated_slits, 
 const mosca::wavelength_calibration& wave_cal,
 const mosca::grism_config& grism_cfg,
 fors_image *master_bias, const cpl_frame * bias_frame,
 struct fors_calib_config& config, cpl_frameset * frameset,
 const char * flat_tag, 
 double nonlinear_level, double max_nonlinear_ratio,
 std::vector<std::vector<double> >& slit_sat_ratio,
 std::vector<std::vector<int> >& slit_sat_count,
 cpl_mask **& nonlinear_flat_masks,
 cpl_mask **& saturated_flat_masks,
 std::auto_ptr<mosca::fiera_config>& ccd_config)
{
    const char     * recipe_name = "fors_calib";
    cpl_errorstate   error_prevstate = cpl_errorstate_get();
    std::auto_ptr<mosca::image> master_flat;

    cpl_msg_indent_more();

    ccd_config = 
           fors_ccd_config_read(cpl_frameset_find_const(frameset, flat_tag), 
                                bias_frame);

    if(ccd_config.get() == 0)
    {
        cpl_msg_error(recipe_name, "Cannot get CCD configuration from header "
                                   "or RON from master bias"
                                   "(missing QC DET OUT? RON keywords)");
        return master_flat;
    }

    /* Get the flat frames */
    cpl_frameset * flatframes = fors_frameset_extract(frameset, flat_tag);
    size_t nflats = cpl_frameset_get_size(flatframes);

    /* Allocate the bad pixel masks*/
    nonlinear_flat_masks = 
            (cpl_mask **)cpl_malloc(nflats* sizeof(cpl_mask*));; 
    saturated_flat_masks = 
            (cpl_mask **)cpl_malloc(nflats* sizeof(cpl_mask*));; 

    /* Reading individual raw flats */
    //TODO: This has copy overhead. Substitute with shared_ptr
    std::vector<mosca::image> basiccal_flats;
    for (size_t i_flat = 0; i_flat < nflats; i_flat++)
    {
        cpl_frame * flatframe = cpl_frameset_get_position(flatframes, i_flat);
        fors_image * flat_raw = fors_image_load(flatframe);

        if (!flat_raw)
            return master_flat;
        
        /* Create variances map */
        fors_image_variance_from_detmodel(flat_raw, *ccd_config);
        if(!cpl_errorstate_is_equal(error_prevstate))
            return master_flat;

        /* Get the non-linear pixels */
        nonlinear_flat_masks[i_flat] = 
                cpl_mask_threshold_image_create(flat_raw->data,
                                nonlinear_level, std::numeric_limits<double>::max());
        /* Get the A/D saturated pixels */
        saturated_flat_masks[i_flat] = 
                cpl_mask_threshold_image_create(flat_raw->data,
                                65535., std::numeric_limits<double>::max());
        cpl_mask * saturated_0 = 
                cpl_mask_threshold_image_create(flat_raw->data,
                                -std::numeric_limits<double>::max(), 
                                std::numeric_limits<double>::min());
        cpl_mask_or(saturated_flat_masks[i_flat], saturated_0);
        cpl_mask_delete(saturated_0);

        /* Subtract overscan */
        fors_image * flat = fors_subtract_prescan(flat_raw, *ccd_config);
        if(!cpl_errorstate_is_equal(error_prevstate))
            return master_flat;

        /* Trimm pre/overscan */
        fors_trimm_preoverscan(flat, *ccd_config);
        fors_trimm_preoverscan(nonlinear_flat_masks[i_flat], *ccd_config);
        fors_trimm_preoverscan(saturated_flat_masks[i_flat], *ccd_config);
        fors_image_delete(&flat_raw);
        if(!cpl_errorstate_is_equal(error_prevstate))
            return master_flat;

        /* Subtract master bias */
        fors_subtract_bias(flat, master_bias);
        if(!cpl_errorstate_is_equal(error_prevstate))
            return master_flat;

        /* Transforming into mosca::image, which takes ownership */
        cpl_image * flat_data = flat->data; 
        cpl_image * flat_err = flat->variance;
        cpl_image_power(flat_err, 0.5);

        mosca::image new_flat(flat_data, flat_err, true, mosca::X_AXIS);
        basiccal_flats.push_back(new_flat);
        //Only the structure is freed, the images are taken over by new_flat
        cpl_free(flat); 
    }

    if(!cpl_errorstate_is_equal(error_prevstate))
    {
        cpl_msg_error(recipe_name, "Could not read the flats");
        return master_flat;
    }
    
    /* Reject slits that have too many non-linear pixels */
    cpl_msg_info(cpl_func, "Computing saturation of flats");
    cpl_msg_indent_more();
    fors_saturation_reject_sat_slits(basiccal_flats, calibrated_slits, 
                                     nonlinear_flat_masks, 
                                     saturated_flat_masks, max_nonlinear_ratio,
                                     slit_sat_ratio, slit_sat_count);
    cpl_msg_indent_less();

    /* Computing master flat */
    cpl_msg_info(cpl_func, "Computing master flat");
    std::string stacking_method(config.stack_method);
    if(stacking_method == "mean" || stacking_method == "sum")
    {
        //TODO: Hardcoded value!! 
        int smooth_size = 10; 
        mosca::reduce_mean reduce_method;
        master_flat = mosca::flat_combine<float, mosca::reduce_mean>
            (basiccal_flats, calibrated_slits, wave_cal, grism_cfg, smooth_size, reduce_method);
        if(stacking_method == "sum")
        {
            cpl_image_multiply_scalar(master_flat->get_cpl_image(), nflats);
            cpl_image_multiply_scalar(master_flat->get_cpl_image_err(), nflats);
        }
    }
    else if(stacking_method == "median")
    {
        //TODO: Hardcoded value!! 
        int smooth_size = 10; 
        mosca::reduce_median reduce_method;
        master_flat = mosca::flat_combine<float, mosca::reduce_median>
            (basiccal_flats, calibrated_slits, wave_cal, grism_cfg, smooth_size, reduce_method);        
    }
    else if(stacking_method == "ksigma")
    {
        //TODO: Hardcoded value!! 
        int smooth_size = 10; 
        mosca::reduce_sigma_clipping reduce_method(config.khigh, config.klow, config.kiter);
        master_flat = mosca::flat_combine<float, mosca::reduce_sigma_clipping>
            (basiccal_flats, calibrated_slits, wave_cal, grism_cfg, smooth_size, reduce_method);    
    }

    //Cleanup
    cpl_frameset_delete(flatframes);

    cpl_msg_indent_less();
    return master_flat;
}

void fors_calib_qc_saturation
(cpl_propertylist * header, const fors::detected_slits detected_slits,
 const std::vector<std::vector<double> >& slit_sat_ratio,
 const std::vector<std::vector<int> >& slit_sat_count)
{
    size_t n_slits = slit_sat_ratio.size();
    size_t n_flats = slit_sat_ratio[0].size();
    std::vector<double> flat_sat_total_sat(n_flats);
    for(size_t i_slit = 0; i_slit < n_slits; i_slit++)
    {
        int slit_id = detected_slits[i_slit].slit_id();
        for(size_t i_flat = 0; i_flat < n_flats; i_flat++)
        {
            flat_sat_total_sat[i_flat] += slit_sat_count[i_slit][i_flat];
            char * keyname;
            keyname = cpl_sprintf("ESO QC FLAT%02zd SLIT%02d SAT RATIO",
                        i_flat + 1, slit_id);
            cpl_propertylist_append_double(header, keyname, 
                                           slit_sat_ratio[i_slit][i_flat]);
            cpl_free(keyname);
            keyname = cpl_sprintf("ESO QC FLAT%02zd SLIT%02d SAT COUNT",
                        i_flat + 1, slit_id);
            cpl_propertylist_append_double(header, keyname, 
                                           slit_sat_count[i_slit][i_flat]);
            cpl_free(keyname);
        }
    }
    
    for(size_t i_flat = 0; i_flat < n_flats; i_flat++)
    {
        char * keyname;
        keyname = cpl_sprintf("ESO QC FLAT%02zd SAT COUNT", i_flat + 1);
        cpl_propertylist_append_double(header, keyname, 
                                       flat_sat_total_sat[i_flat]);
        cpl_free(keyname);        
    }
}