/*                                                                              *
 *   This file is part of the ESO UVES Pipeline                                 *
 *   Copyright (C) 2004,2005 European Southern Observatory                      *
 *                                                                              *
 *   This library 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, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA       *
 *                                                                              */
 
/*
 * $Author: amodigli $
 * $Date: 2013-08-08 13:36:46 $
 * $Revision: 1.85 $
 * $Name: not supported by cvs2svn $
 *
 */
#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

/*----------------------------------------------------------------------------*/
/**
 * @defgroup uves_response  Recipe: Response
 *
 * This recipe calculates the response function and the quantum detection
 * efficiency.
 * See man-page for details.
 */
/*----------------------------------------------------------------------------*/

/*-----------------------------------------------------------------------------
                                Includes
 -----------------------------------------------------------------------------*/

/* Definitions */
#include <uves.h>

/* Macro steps */
#include <uves_reduce.h>
#include <uves_reduce_utils.h>
#include <uves_response_efficiency.h>
#include <uves_response_utils.h>

/* Utility functions */
#include <uves_extract.h>
#include <uves_plot.h>
#include <uves_dfs.h>
#include <uves_pfits.h>
#include <uves_parameters.h>
#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_utils_cpl.h>
#include <uves_qclog.h>
#include <uves_recipe.h>
#include <uves_error.h>
#include <uves_msg.h>
#include <uves_merge.h>

/* Library */
#include <cpl.h>
#include <stdbool.h>
#include <string.h>

/*-----------------------------------------------------------------------------
                            Functions prototypes
 -----------------------------------------------------------------------------*/
static void uves_efficiency_qclog(cpl_table* table,
                  uves_propertylist* raw_header, 
                  enum uves_chip chip,
                  cpl_table* qclog,
                  const char *ref_obj_name);

static int
uves_response_define_parameters(cpl_parameterlist *parameters);

/*-----------------------------------------------------------------------------
                            Recipe standard code
 -----------------------------------------------------------------------------*/
#define cpl_plugin_get_info uves_response_get_info
UVES_RECIPE_DEFINE(
    UVES_RESPONSE_ID, UVES_RESPONSE_DOM, uves_response_define_parameters,
    "Jonas M. Larsen", "cpl@eso.org",
    "Determines response function and quantum efficiency",
"This recipe reduces a standard star frame (STANDARD_xxx or STANDARD_xxx,\n"
"where xxx = BLUE, RED) using a combination (depending on recipe parameters\n"
"and provided input frames) of the steps:\n"
"  - bias subtraction,\n"
"  - dark subtraction,\n"
"  - background subtraction,\n"
"  - extraction/cosmic ray removal,\n"
"  - flat-field correction,\n"
"  - wavelength rebinning,\n"
"  - sky subtraction,\n"
"  - order merging.\n"
"\n"
" Expected input for this recipe is an raw std star frame, STANDARD_xxx or \n"
"order table(s) for each chip, ORDER_TABLE_xxxx (where xxxx=BLUE, REDL, REDU),\n"
"line table(s) for each chip, LINE_TABLE_xxxx, a master bias frame,\n"
"MASTER_BIAS_xxxx, a master flat, MASTER_FLAT_xxxx, a reference standard star\n"
"flux table, FLUX_STD_TABLE, a table describing the atmospheric extintion,\n"
"EXTCOEFF_TABLE. \n"

"Two reductions are performed, the first using optimal extraction (used to\n"
"compute the instrument response function), the second using linear extraction\n"
"(used to get the Quantum Detection Efficiency)\n"
"\n"
"For each chip (xxxx = BLUE, REDL, REDU) the recipe produces\n"
"  INSTR_RESPONSE_xxxx          Response curve\n"
"  WCALIB_FF_RESPONSE_xxxx      Response curve in (lambda,order) space before\n"
"                               correcting for exposure time, gain, binning and\n"
"                               atmospheric absorption\n"
"  RED_STD_xxxx                 Reduced spectrum\n"
"  EFFICIENCY_TABLE_xxxx        Efficiency table\n"
    "  BKG_STD_xxxx                 The subtracted background\n");

/**@{*/
/*----------------------------------------------------------------------------*/
/**
  @brief    Setup the recipe options    
  @param    parameters        the parameterlist to fill
  @return   0 if everything is ok
 */
/*----------------------------------------------------------------------------*/
static int
uves_response_define_parameters(cpl_parameterlist *parameters)
{
    /*****************
     *    General    *
     *****************/

    if (uves_define_global_parameters(parameters) != CPL_ERROR_NONE)
    {
        return -1;
    }
    
    /*******************
     *    Reduce       *
     ******************/
    /* Get reduce parameters for the top level and also for 'efficiency' substep */
    if (uves_propagate_parameters_step(
        UVES_REDUCE_ID, parameters, make_str(UVES_RESPONSE_ID), NULL) != 0)
    {
        return -1;
    }


    check(uves_define_efficiency_parameters(parameters),
          "Defining efficiency parameters");


  cleanup:
    return (cpl_error_get_code() != CPL_ERROR_NONE);
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Reduce one chip of a UVES std frame
  @param    raw_image           The raw image
  @param    raw_header          FITS header of raw image
  @param    rotated_header      Header describing the geometry of the raw image after
                                rotation and removal of pre- and overscan areas
  @param    master_bias         The master bias image for this chip, or NULL
  @param    master_dark         The master dark image for this chip, or NULL
  @param    mflat_header        FITS header of master flat frame
  @param    master_flat         The master flat image for this chip, or NULL
  @param    mdark_header        FITS header of master dark frame
  @param    ordertable          Order table describing the order locations on the raw image
  @param    order_locations     The polynomial describing the order positions
  @param    linetable           Length 3 array of linetable for sky, object, sky.
  @param    linetable_header    Length 3 array of linetable headers for sky, object, sky.
  @param    dispersion_relation Length 3 array of dispersion relations for sky, object, sky.
  @param    flux_table          Catalogue standard star fluxes
  @param    atm_extinction      Atmospheric extinction coefficients
  @param    chip                CCD chip
  @param    debug_mode               If set to true, intermediate results are saved to 
                                the current directory
  @param    parameters          The recipe parameter list containing parameters
                                for background subtraction, flat-fielding, extraction,
                rebinning
  @param    calc_response       Also calculate response curve (or just do the reduction)?
  @param    PACCURACY           Pointing accuracy (in arcseconds) used to identify object
  @param    ref_obj_id          (out) Reference object ID
  @param    reduced_spectrum    (out) The extracted, flat-fielded, rebinned, merged, 
                                sky-subtracted (but not flux calibrated) science spectrum
  @param    reduced_header      (out) Header describing the geometry of the @em reduced_spectrum
  @param    background          (out) The background image that was subtracted
  @param    response_orders     (out) The response curve per order (2d image), which was
                                not corrected for exposure time, gain, binning, atm. absorption.
  @param    response_header_2d  (out) Header describing the geometry of @em response_orders
  @param    response_curve      (out) The (merged and therefore 1d) response curve
  @param    response_header     (out) Header describing the geometry of @em response_curve
  @param    blaze_efficiency    (out) Efficiency at blaze function maximum, for each order
  @param    efficiency          (out) The quantum detection efficiency table
  @param    extraction_slit     (out) extraction slit length used
  @return   The reduced spectrum

  This function reduces the provided @em raw_image (see @c uves_reduce) and
  computes the response function as

  @em R ( @em lambda ) = @em flux_cat / @em flux_std

  where @em flux_cat is the tabulated standard star flux, and @em flux_std is
  the observed flux (corrected for exposure time, gain, binning and moved 
  to the top of the atmosphere)

  Additionally the instrument Quantum Detection Efficiency is computed.
  (see @c uves_response_efficiency)
**/
/*----------------------------------------------------------------------------*/

static cpl_error_code
uves_response_process_chip(const cpl_image *raw_image, 
                           const uves_propertylist *raw_header, 
                           const uves_propertylist *rotated_header,
               const cpl_image *master_bias,
                           const uves_propertylist *mbias_header,
               const cpl_image *master_dark, 
                           const uves_propertylist *mdark_header, 
               const cpl_image *master_flat, 
                           const uves_propertylist *mflat_header,
               const cpl_table *ordertable, 
                           const polynomial *order_locations,
               const cpl_table *linetable[3], 
                           const uves_propertylist *linetable_header[3], 
                           const polynomial *dispersion_relation[3],
               const cpl_table *flux_table,
               const cpl_table *atm_extinction,
               enum uves_chip chip,
               /* General */
               bool   debug_mode,
               /* Backsub */
               /* Flat fielding */
               /* Extraction */
               /* Rebinning  */
               const cpl_parameterlist *parameters,
               bool calc_response,
               /* Identification */
               double PACCURACY,
               /* Output */
               char **ref_obj_id,
               cpl_image **reduced_spectrum, 
                           uves_propertylist **reduced_header, 
                           cpl_image **background,
               cpl_image **response_orders, 
                           uves_propertylist **response_header_2d,
               cpl_image **response_curve,  
                           uves_propertylist **response_header,
               cpl_table **efficiency, 
                           cpl_table** blaze_efficiency,
                           cpl_table** info_tbl,
               double *extraction_slit)
{
    cpl_image        *rebinned_spectrum = NULL;
    cpl_image        *rebinned_noise    = NULL;
    uves_propertylist *rebinned_header   = NULL;
    cpl_image        *reduced_rebinned = NULL;
    cpl_image        *reduced_rebinned_noise = NULL;

    cpl_image        *reduced_noise    = NULL;

    cpl_table *cosmic_mask     = NULL;   /* Cosmic ray table  (not a product of this recipe) */
    cpl_table *order_trace     = NULL;   /* Order trace table (not a product of this recipe) */
    cpl_image *merged_spectrum = NULL;   /* Not sky-subtracted (if simple extraction)        */
    cpl_image *wave_map      = NULL;
    cpl_image *merged_sky      = NULL;
    cpl_image *merged_noise    = NULL;
    cpl_image *reduced_scaled  = NULL;   /* Merged, sky-subtracted and normalized 
                        (exposure time, gain...) */

    cpl_table *catalogue_flux = NULL;    /* Std star catalogue flux */
 
    /* Do the science reduction. Produces wave.cal. spectra. */
    uves_msg("Reducing standard star");

    check( uves_reduce(raw_image, 
                       raw_header, 
                       rotated_header,
               master_bias,
                       mbias_header,
               master_dark, 
                       mdark_header,
               master_flat, 
                       mflat_header,
               ordertable, 
                       order_locations,
               linetable, 
                       linetable_header, 
                       dispersion_relation,
               chip,
               debug_mode, 
               parameters, 
                       make_str(UVES_RESPONSE_ID),
		       "",
               /* Output */
               NULL, 
                       NULL, 
                       NULL,                   /* 2d products */
               &cosmic_mask,
               &wave_map,
               background,
               NULL, 
                       NULL,                          /* Variance of flat-fielded */
               NULL, 
                       NULL,                         /* Don't need these 
                                  intermediate products */
               &merged_sky,
               &rebinned_spectrum, 
                       &rebinned_noise, 
                       &rebinned_header,
               &merged_spectrum,   
                       &merged_noise, 
                       reduced_header,
               &reduced_rebinned,  
                       &reduced_rebinned_noise,
               reduced_spectrum,   
                       &reduced_noise,
                       info_tbl,
               extraction_slit,
                       &order_trace),
       "Could not reduce frame");


    check( uves_plot_image_rows(*reduced_spectrum, 1, 1, 1,
                "Wavelength (arbitrary units)", NULL,
                "Reduced spectrum (%s chip)", 
                uves_chip_tostring_upper(chip)),
       "Plotting failed");

    if (calc_response)
    {
        /* Calculate 2d response curve  (but don't scale to unit exposure time, ...) */
        uves_msg("Filtering rebinned spectrum");
        check( uves_filter_image_median(&reduced_rebinned, 10, 0, false), 
           "Could not smooth spectrum");
        check( uves_filter_image_average(reduced_rebinned, 10, 0), 
           "Could not smooth spectrum");    
        
        uves_msg("Calculating 2d response curve");
        
        check( *response_orders = uves_calculate_response(reduced_rebinned, 
                                  rebinned_header,
                                  flux_table,
                                  raw_header, 
                                  PACCURACY,
                                  true,     
                                  /*  std_flux / flux  */
                                  ref_obj_id),
           "Could not calculate response curve");
        
        check( *response_header_2d = uves_propertylist_duplicate(rebinned_header),
           "Error creating FITS header for 2d response curve");
        
        check( uves_pfits_set_bunit(*response_header_2d, "FLUX_CAT / FLUX_STD"),
           "Error writing BUNIT keyword");
        
        /*
         *  Calculate 1d response curve
         */
        
        uves_msg("Normalizing reduced spectrum");
        
        {
        int n_traces = cpl_image_get_size_y(*reduced_spectrum);
        assure( n_traces == 1, CPL_ERROR_ILLEGAL_INPUT,
            "2d extraction/reduction not supported");
        
        check( reduced_scaled = uves_normalize_spectrum(*reduced_spectrum, 
                                reduced_noise, 
                                *reduced_header,
                                raw_header,
                                n_traces,
                                chip,
                                atm_extinction,
                                true,  /* Yes, divide by binning */
                                NULL), /* Noise spectrum         */
               "Error normalizing reduced spectrum");
        }
        
        uves_msg("Filtering reduced spectrum");
        check( uves_filter_image_median(&reduced_scaled, 10, 0, false), 
           "Could not smooth spectrum");
        check( uves_filter_image_average(reduced_scaled, 10, 0), 
           "Could not smooth spectrum");
        
        uves_msg("Calculating response curve from scaled spectrum");
        
        cpl_free(*ref_obj_id); *ref_obj_id = NULL;
        check( *response_curve = uves_calculate_response(reduced_scaled, 
                                 *reduced_header,
                                 flux_table,
                                 raw_header, 
                                 PACCURACY, 
                                 true, /*  catalogue_flux / flux  */
                                 ref_obj_id),
           "Could not calculate response curve");
        
        check( *response_header = uves_propertylist_duplicate(*reduced_header),
           "Error creating FITS header for response curve");
        
        check( uves_pfits_set_bunit(*response_header, "FLUX_CAT / FLUX_STD"),
           "Error writing BUNIT keyword");
        
        if (debug_mode)
        {
            check( uves_save_image_local("Pre-smoothed response curve", "raw_response", 
                         *response_curve, chip, -1, -1, *response_header, true), 
               "Error saving image");
        }
        
        check( uves_plot_image_rows(*response_curve, 1, 1, 1, "Wavelength (arbitrary units)", 
                    NULL,
                    "Raw response (%s chip)", uves_chip_tostring_upper(chip)),
           "Plotting failed");
        
        
        /* Rebin the response curve to 50 wlu:
           1) smooth it using a radius of    25 wlu,
           2) then extract every n'th pixel where n = 50/step */
        uves_msg("Rebinning response curve to step size = 50 w.l.u.");
        {
        double dlambda, lambda_start;
        int n, bin, newbin;
    
        check( lambda_start = uves_pfits_get_crval1(*response_header),
               "Error reading start wavelength from header");
        check( dlambda = uves_pfits_get_cdelt1(*response_header),
               "Error reading wavelength step from header");
        
        n = uves_round_double(50.0/dlambda);
        
        assure( n >= 1, CPL_ERROR_ILLEGAL_OUTPUT,
            "Cannot rebin to 50 w.l.u. Current step is only %f w.l.u.", dlambda);
        
        /* Filter radius = 25 wlu, 0    (It's a 1d image) */
        check( uves_filter_image_average(*response_curve, n/2, 0),
               "Error filtering response curve");
        
        newbin = 1;
        for (bin = 1+n/2; bin <= cpl_image_get_size_x(*response_curve); bin += n)
            {
            int pis_rejected;
            
            /* Write to the same image buffer */
            cpl_image_set(*response_curve, 
                      newbin, 1,
                      cpl_image_get(*response_curve, bin, 1, &pis_rejected)
                );
            newbin++;
            }
        
        
        /* Extract image, update start+step wavelengths */
        uves_crop_image(response_curve, 1, 1, newbin-1, 1);
        
        lambda_start = lambda_start + dlambda * ((1+n/2) - 1);  /* Center of first bin */
        dlambda      = n * dlambda;
        
        check( uves_pfits_set_crval1(*response_header, lambda_start),
               "Error updating start wavelength");
        check( uves_pfits_set_cdelt1(*response_header, dlambda),
               "Error updating wavelength step");
        
        }
        
        check( uves_plot_image_rows(*response_curve, 1, 1, 1, "Wavelength (arbitrary units)", 
                    NULL,
                    "Response curve (%s chip)", 
                    uves_chip_tostring_upper(chip)),
           "Plotting failed");
        
        /* Calculate efficiency table */
        uves_msg("Calculating efficiency curve");
        
        check( uves_response_efficiency(raw_image, 
                        raw_header, 
                        rotated_header,
                        master_bias,
                        mbias_header,
                        master_dark, 
                        mdark_header, 
                        ordertable, 
                        order_locations,
                        linetable, 
                        linetable_header, 
                        dispersion_relation,
                        flux_table,
                        atm_extinction,
                        chip,
                        debug_mode,
                        parameters,
                        PACCURACY,
                        efficiency, 
                        blaze_efficiency),
           "Efficiency calculation failed");
        
        check( uves_plot_table(*efficiency, "Wave", "Eff",
                   "Detection Quantum Efficiency (%s chip)", 
                   uves_chip_tostring_upper(chip)), 
           "Plotting failed");
        
        /* Save blaze function efficiency (efficiency at center of order) */
        if (debug_mode) check( uves_save_table_local("Blaze efficiency table",
                            "blaze_efficiency", 
                            *blaze_efficiency, chip, -1, -1, rotated_header, NULL),
                  "Error saving blaze efficiency table");
    }
    else
    {
        uves_msg("Skipping response/efficiency computation");
    }

  cleanup:
    uves_free_propertylist(&rebinned_header);
    uves_free_image(&rebinned_noise);
    uves_free_image(&rebinned_spectrum);
    uves_free_table(&cosmic_mask);
    uves_free_table(&order_trace);
    uves_free_image(&merged_spectrum);
    uves_free_image(&merged_noise);
    uves_free_image(&merged_sky);
    uves_free_image(&reduced_rebinned);
    uves_free_image(&reduced_rebinned_noise);
    uves_free_image(&reduced_noise);
    uves_free_image(&reduced_scaled);
    uves_free_table(&catalogue_flux);

    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        /* Output */
        uves_free_image(reduced_spectrum);
        uves_free_image(background);
        uves_free_image(response_orders);
        uves_free_image(response_curve);
        uves_free_propertylist(reduced_header);
        uves_free_propertylist(response_header);
        uves_free_propertylist(response_header_2d);
        uves_free_table(efficiency);
    }
    
    return cpl_error_get_code();
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Get the command line options and execute the data reduction
  @param    parameters  the parameters list
  @param    frames      the frames list
  @return   CPL_ERROR_NONE if everything is ok
 */
/*----------------------------------------------------------------------------*/
static void
UVES_CONCAT2X(UVES_RESPONSE_ID,exe)(cpl_frameset *frames,
          const cpl_parameterlist *parameters,
          const char *starttime)
{
    /*
     * Variables containg the values of recipe parameters 
     */

    /* General */
    bool debug_mode;

    /* Background subtraction */
    /* Implicitly passed */

    /* Flat-fielding */

    /* Extraction */
    /* Implicitly passed */
    
    /* Rebinning */
    /* Implicitly passed */

    /* Efficiency */
    double PACCURACY;

    /* CPL objects */
    /* Input, raw */
    cpl_image        *raw_image[2]      = {NULL, NULL};
    uves_propertylist *raw_header[2]     = {NULL, NULL};
    uves_propertylist *rotated_header[2] = {NULL, NULL};

    /* Input, calib */
    cpl_image        *master_bias        = NULL;
    uves_propertylist *master_bias_header = NULL;

    cpl_image        *master_dark        = NULL;
    uves_propertylist *master_dark_header = NULL;

    cpl_image        *master_flat        = NULL;
    uves_propertylist *master_flat_header = NULL;

    cpl_table        *ordertable       = NULL;
    uves_propertylist *ordertable_header= NULL;
    polynomial       *order_locations  = NULL;
    cpl_table        *traces           = NULL;

    cpl_table        *flux_table       = NULL;

    cpl_table        *atm_extinction   = NULL;
    uves_propertylist *table_header= NULL;
    /* Line tables for sky, object, sky (UVES specific) */
    const cpl_table        *linetable[3]           = {NULL, NULL, NULL};
    const uves_propertylist *linetable_header[3]    = {NULL, NULL, NULL};
    const polynomial       *dispersion_relation[3] = {NULL, NULL, NULL};

    /* Output */
    cpl_image        *background         = NULL;
    cpl_image        *reduced_spectrum   = NULL;
    uves_propertylist *spectrum_header    = NULL;
    cpl_image        *response_orders    = NULL;
    uves_propertylist *response_header_2d = NULL;
    cpl_image        *response_curve     = NULL;
    uves_propertylist *response_header    = NULL;
    cpl_table        *efficiency         = NULL;
    cpl_table        *blaze_efficiency   = NULL;
    uves_propertylist *efficiency_header  = NULL;
    cpl_table* info_tbl=NULL;

    /* Local variables */
    cpl_table *qclog[2] = {NULL, NULL};
    cpl_table *qclog_optext[2] = {NULL, NULL};
    cpl_table *catalogue_flux = NULL;
    const char *raw_filename = "";        /* Static */
    const char *flux_table_filename = ""; /* Static */
    const char *atm_ext_filename = "";    /* Static */
    char *product_filename = NULL;        /* Dynamically allocated */
    char *prod_catg = NULL; 
    char *ref_obj_name = NULL;            /* Reference object id */
    bool calc_response = false;           /* Calculate instr response? */
    double extraction_slit;

    bool blue  = false;
    enum uves_chip chip;
    int binx = 0;
    int biny = 0;
    const char* PROCESS_CHIP=NULL;

    const char *ordertable_filename = "";
    const char *linetable_filename = "";
    const char *master_bias_filename = "";
    const char *master_dark_filename = "";
    const char *master_flat_filename = "";
    const char *chip_name = "";
    int tracerow=0;                  /* Index of table row */
        
    int raw_index = 0;

    char extname[80];
    double trace_offset=0;
    int trace_number=0;
    int trace_enabled=0;
    int window=0;          /* window number */
    merge_method m_method;

    /* Read recipe parameters */
    {
    /* General */
    check( uves_get_parameter(parameters, NULL, "uves", "debug", 
                  CPL_TYPE_BOOL  , &debug_mode      ), "Could not read parameter");

    check( uves_get_parameter(parameters, NULL, "uves", "process_chip", CPL_TYPE_STRING, &PROCESS_CHIP),
               "Could not read parameter");
    uves_string_toupper((char*)PROCESS_CHIP);

    /* Background subtraction, Flat-fielding, Rebinning */
    /* The input parameter list is passed */
    
    /* For both response curve and efficiency step:
       Allow only extraction methods average/linear/optimal */
    {
        extract_method em;

        /* Validate uves_response.reduce.extract.method */
        check( em = uves_get_extract_method(
               parameters, NULL,
               make_str(UVES_RESPONSE_ID) "." UVES_REDUCE_ID "." UVES_EXTRACT_ID),
           "Could not read extraction method");
        
        assure( em == EXTRACT_LINEAR || em == EXTRACT_AVERAGE || em == EXTRACT_OPTIMAL,
            CPL_ERROR_UNSUPPORTED_MODE, 
            "Use linear/average/optimal extraction method to calculate response curve");
        
        /* Validate uves_response.efficiency.reduce.extract.method */
        check( em = uves_get_extract_method(
               parameters, NULL,
               make_str(UVES_RESPONSE_ID) ".efficiency." UVES_REDUCE_ID "." UVES_EXTRACT_ID),
           "Could not read extraction method");
        
        assure( em == EXTRACT_LINEAR || em == EXTRACT_AVERAGE || em == EXTRACT_OPTIMAL,
            CPL_ERROR_UNSUPPORTED_MODE, 
            "Use linear/average/optimal extraction "
            "method to calculate quantum efficiency");
    }
    
    /* Identification */
    check( uves_get_parameter(parameters, NULL, 
                  make_str(UVES_RESPONSE_ID) ".efficiency", "paccuracy", 
                  CPL_TYPE_DOUBLE, &PACCURACY), 
           "Could not read parameter");
    }
    
    /* Load raw image and header, and identify input frame as red or blue */
    check( uves_load_standard(frames,
                  &raw_filename, raw_image, raw_header, rotated_header, 
                  &blue), 
       "Error loading raw frame");

    /* Load flux table */
    check( uves_load_flux_table(frames, &flux_table_filename, &flux_table),
       "Error loading standard flux table");

    uves_msg_low("Using standard star flux table in '%s'", flux_table_filename);

    /* Before doing the reduction, find out if the standard star is in the table.
     * If not, still do the science reduction, but skip the instr.response part
     */
    catalogue_flux = uves_align(raw_header[0], flux_table, PACCURACY, &ref_obj_name);

    calc_response = true;
    if (cpl_error_get_code() == CPL_ERROR_INCOMPATIBLE_INPUT)
    {
        uves_error_reset();

        uves_msg_warning("No catalogue object found within %.2f arcsecs. "
                 "Instrument response curve will not be computed",
                 PACCURACY);

        calc_response = false;
    }

    /* Load atmospheric extinction table */
    check( uves_load_atmo_ext(frames, &atm_ext_filename, &atm_extinction), 
       "Error loading extinction coefficients");
    
    uves_msg_low("Using atmospheric extinction table in '%s'", atm_ext_filename);

    //cpl_parameterlist_dump(parameters,stdout);
    //uves_msg("recipe=%s",UVES_RESPONSE_ID);
   check( m_method = uves_get_merge_method(parameters, "uves_cal_response", "reduce"),
           "Could not get merging method");
    /* Adjust parameters according to binning 
     */
    check (binx = uves_pfits_get_binx(raw_header[0]), 
       "Could not read x binning factor from input header");
    check (biny = uves_pfits_get_biny(raw_header[0]),
       "Could not read y binning factor from input header");

    /* Loop over one or two chips, over traces and
       over extraction windows */
    for (chip = uves_chip_get_first(blue);
     chip != UVES_CHIP_INVALID;
     chip = uves_chip_get_next(chip))
    {
      table_header = uves_propertylist_new();


      if(strcmp(PROCESS_CHIP,"REDU") == 0) {
	chip = uves_chip_get_next(chip);
      }


        /* const char *drs_filename        = "";    not used */
        raw_index = uves_chip_get_index(chip);
        uves_msg("Processing %s chip in '%s'",
             uves_chip_tostring_upper(chip), raw_filename);
        
        check_nomsg( chip_name = uves_pfits_get_chipid(raw_header[raw_index], chip));

        uves_msg_debug("Binning = %dx%d", binx, biny);
        
        /* Load master bias, set pointer to NULL if not present */
        uves_free_image(&master_bias);
        uves_free_propertylist(&master_bias_header);
        if (cpl_frameset_find(frames, UVES_MASTER_BIAS(chip)) != NULL)
        {
            check( uves_load_mbias(frames, chip_name, 
                       &master_bias_filename,
                       &master_bias,
                       &master_bias_header, chip), 
               "Error loading master bias");
            
            uves_msg_low("Using master bias in '%s'", master_bias_filename);
        }
        else
        {
            uves_msg_low("No master bias in SOF. Bias subtraction not done");
        }
        
        
        /* Load master dark, set pointer to NULL if not present */
        uves_free_image(&master_dark);
        uves_free_propertylist(&master_dark_header);
        if (cpl_frameset_find(frames, UVES_MASTER_DARK(chip)) != NULL)
        {
            check( uves_load_mdark(frames, chip_name,
                       &master_dark_filename,
                       &master_dark,
                       &master_dark_header, chip),
               "Error loading master dark");
            
            uves_msg_low("Using master dark in '%s'", master_dark_filename);
        }
        else
        {
            uves_msg_low("No master dark in SOF. Dark subtraction not done");
        }
        
        /* Load master flat */
        uves_free_image(&master_flat);
        uves_free_propertylist(&master_flat_header);
        check( uves_load_mflat_const(frames, chip_name, 
                     &master_flat_filename,
                     &master_flat, 
                     &master_flat_header, 
                     chip, NULL), "Error loading master flat");
        
        uves_msg_low("Using master flat in '%s'", master_flat_filename);
        
        
        /* Load the order table for this chip */
        uves_free_table       (&ordertable);
        uves_free_propertylist(&ordertable_header);
        uves_polynomial_delete(&order_locations);
        uves_free_table       (&traces);

        check( uves_load_ordertable(frames, 
                    false,  /* FLAMES? */
                    chip_name,
                    &ordertable_filename, 
                    &ordertable, 
                    &ordertable_header,
                                        NULL,
                    &order_locations, &traces, 
                    NULL, NULL,
                                       NULL, NULL, /* fibre_pos,fibre_mask */
                    chip, false),
           "Could not load order table");
        uves_msg_low("Using order table in '%s'", ordertable_filename);
        
        
        /* Loop over all traces (1 trace for UVES) */
        for(tracerow = 0; tracerow < cpl_table_get_nrow(traces); tracerow++)
        {

            
            trace_offset  = cpl_table_get_double(traces, "Offset"    , tracerow, NULL);
            trace_number  = cpl_table_get_int   (traces, "TraceID"   , tracerow, NULL);
            trace_enabled = cpl_table_get_int   (traces, "Tracemask" , tracerow, NULL);
            
            if (trace_enabled != 0)
            {
                
                            if (cpl_table_get_nrow(traces) > 1) {
                                uves_msg("Processing trace %d", trace_number);
                            }
                
                /* This is UVES specific. Load linetable for the 
                   two sky windows (number 1, 3) and for the object
                   window (number 2) */
                
                for (window = 1; window <= 3; window ++)
                {
                    uves_free_table_const ( &(linetable[window-1]) );
                    uves_free_propertylist_const( &(linetable_header[window-1]) );
                    uves_polynomial_delete_const( &(dispersion_relation[window-1]) );
                    check( uves_load_linetable_const(
                           frames, 
                           false,  /* FLAMES? */
                           chip_name,
                           order_locations,
                           cpl_table_get_column_min(ordertable, "Order"),
                           cpl_table_get_column_max(ordertable, "Order"),
                           &linetable_filename,
                           &(linetable          [window-1]),
                           &(linetable_header   [window-1]),
                           &(dispersion_relation[window-1]),
                           NULL,
                           chip,
                           trace_number,
                           window),
                       "Could not load line table, window #%d", window);
                }
                
                uves_msg_low("Using line tables in '%s'", linetable_filename);
                /* end, UVES specific */
                
                /* Do the reduction + response calculation */
                cpl_free(ref_obj_name); ref_obj_name = NULL;
                uves_free_image(&reduced_spectrum);
                uves_free_image(&background);
                uves_free_image(&response_orders);
                uves_free_image(&response_curve);
                uves_free_propertylist(&response_header);
                uves_free_propertylist(&spectrum_header);
                uves_free_propertylist(&response_header_2d);
                uves_free_table(&efficiency);
                uves_free_table(&blaze_efficiency);
                uves_free_table(&info_tbl);
                
                check( uves_response_process_chip(
                       raw_image[raw_index],   /* Raw         */
                                       raw_header[raw_index],  
                       rotated_header[raw_index],
                       master_bias,            /* Calibration */
                                       master_bias_header,
                       master_dark, 
                                       master_dark_header,
                       master_flat, 
                                       master_flat_header,
                       ordertable, 
                                       order_locations,
                       linetable, 
                                       linetable_header, 
                                       dispersion_relation,
                       flux_table,
                       atm_extinction,
                       chip,                  /* Parameters  */
                       debug_mode, 
                       parameters,
                       calc_response,
                       PACCURACY,             /* Identification */
                       &ref_obj_name,
                       &reduced_spectrum, 
                                       &spectrum_header, 
                                       &background,
                       &response_orders, 
                                       &response_header_2d,
                       &response_curve,  
                                       &response_header,
                       &efficiency, 
                                       &blaze_efficiency,
                                       &info_tbl,
                       &extraction_slit),
                   "Response computation failed");
                
                uves_msg("Saving products...");
                
                /* Calculate QC (two tables) */
                
                if (calc_response)
                {
                    uves_qclog_delete(&qclog[0]);
                    qclog[0] = uves_qclog_init(raw_header[raw_index], chip);

                    check( uves_efficiency_qclog(blaze_efficiency,
                                 raw_header[raw_index],
                                 chip,
                                 qclog[0],
                                 ref_obj_name), 
                       "Error generating efficiency QC log");
                }
                
                uves_qclog_delete(&qclog_optext[0]);
                        qclog_optext[0] = cpl_table_new(0);
                /* Do not:  
                   qclog_optext[0] = uves_qclog_init(raw_header[raw_index], chip);
                   because we don't want QC.DID for this
                */
                cpl_table_new_column(qclog_optext[0],"key_name", CPL_TYPE_STRING);
                cpl_table_new_column(qclog_optext[0],"key_type", CPL_TYPE_STRING);
                cpl_table_new_column(qclog_optext[0],"key_value",CPL_TYPE_STRING);
                cpl_table_new_column(qclog_optext[0],"key_help", CPL_TYPE_STRING);
                
                check( uves_qclog_add_sci(qclog_optext[0],
                              raw_header[raw_index],
                              raw_image[raw_index],
                              extraction_slit,
                              info_tbl),
                   "Error generating extraction QC log");

                if (calc_response)
                {
                    /* Save response curve */
                    cpl_free(product_filename);
                    check( product_filename = uves_response_curve_filename(chip),
                       "Error getting filename");
                    uves_pfits_set_extname(response_header,"instrument response");

                    check( uves_frameset_insert(
                           frames,
                           response_curve,
                           CPL_FRAME_GROUP_PRODUCT,
                           CPL_FRAME_TYPE_IMAGE,
                           CPL_FRAME_LEVEL_INTERMEDIATE,
                           product_filename,
                           UVES_INSTR_RESPONSE(chip),
                           raw_header[raw_index],
                           response_header,
                           NULL,
                           parameters,
                           make_str(UVES_RESPONSE_ID),
                           PACKAGE "/" PACKAGE_VERSION, 
                           qclog_optext,
                           starttime, false, 
                           UVES_ALL_STATS),
                       "Could not add response curve '%s' (%s) to frameset", 
                       product_filename, UVES_INSTR_RESPONSE(chip));
                    
                    uves_msg("Response curve '%s' (%s) added to frameset",
                         product_filename, UVES_INSTR_RESPONSE(chip));
                    
                    /* Save response curve (2d) */
                    cpl_free(product_filename);
                    check( product_filename = 
                       uves_response_curve_2d_filename(chip), 
                       "Error getting filename");
                    uves_pfits_set_extname(response_header_2d,"2D instrument response");
                    check( uves_frameset_insert(
                           frames,
                           response_orders,
                           CPL_FRAME_GROUP_PRODUCT,
                           CPL_FRAME_TYPE_IMAGE,
                           CPL_FRAME_LEVEL_INTERMEDIATE,
                           product_filename,
                           UVES_WCALIB_FF_RESPONSE(chip),
                           raw_header[raw_index],
                           response_header_2d,
                           NULL,
                           parameters,
                           make_str(UVES_RESPONSE_ID),
                           PACKAGE "/" PACKAGE_VERSION,
                           qclog_optext,
                           starttime, false, 
                           UVES_ALL_STATS),
                       "Could not add response curve (2d) "
                       "'%s' (%s) to frameset",
                       product_filename, UVES_WCALIB_FF_RESPONSE(chip));
            
                    uves_msg("Response curve (2d) '%s' (%s) added to frameset", 
                         product_filename, UVES_WCALIB_FF_RESPONSE(chip));
                }

                /* Save reduced spectrum */
                cpl_free(product_filename);
                if (m_method == MERGE_NOAPPEND) {
               check( product_filename = uves_response_red_noappend_standard_filename(chip),
                   "Error getting filename");    
                   prod_catg=UVES_RED_NOAPPEND_STD(chip);
                } else {

                check( product_filename = uves_response_red_standard_filename(chip),
                   "Error getting filename");    
                   prod_catg=UVES_RED_STD(chip);

                }
                uves_pfits_set_extname(spectrum_header,"reduced spectrum");
                check( uves_frameset_insert(frames,
                            reduced_spectrum,
                            CPL_FRAME_GROUP_PRODUCT,
                            CPL_FRAME_TYPE_IMAGE,
                            CPL_FRAME_LEVEL_INTERMEDIATE,
                            product_filename,
                            prod_catg,
                            raw_header[raw_index],
                            spectrum_header,
                            NULL,
                            parameters,
                            make_str(UVES_RESPONSE_ID),
                            PACKAGE "/" PACKAGE_VERSION,
                            qclog_optext,
                            starttime, false, 
                            UVES_ALL_STATS),
                   "Could not add reduced spectrum '%s' (%s) to frameset",
                   product_filename, UVES_RED_STD(chip));
                
                uves_msg("Reduced spectrum '%s' (%s) added to frameset", 
                     product_filename, UVES_RED_STD(chip));

                if (calc_response)
                {
                    /* Save efficiency table */
                    uves_free_propertylist(&efficiency_header);
                    efficiency_header = uves_propertylist_new();
                    
                    cpl_free(product_filename);
                    check( product_filename = 
                       uves_response_efficiency_filename(chip),
                       "Error getting filename");
                    uves_pfits_set_extname(efficiency_header,"efficiency");
                    sprintf(extname,"EFFICIENCY");
                    uves_pfits_set_extname(table_header,extname);
                    check( uves_frameset_insert(
                           frames,
                           efficiency,
                           CPL_FRAME_GROUP_PRODUCT,
                           CPL_FRAME_TYPE_TABLE,
                           CPL_FRAME_LEVEL_INTERMEDIATE,
                           product_filename,
                           UVES_EFFICIENCY_TABLE(chip),
                           raw_header[raw_index],
                           efficiency_header,
                           table_header,
                           parameters,
                           make_str(UVES_RESPONSE_ID),
                           PACKAGE "/" PACKAGE_VERSION,
                           qclog,
                           starttime, true, 0),
                       "Could not add background image '%s' (%s) to frameset",
                       product_filename, UVES_EFFICIENCY_TABLE(chip));
                    
                    uves_msg("Efficiency table '%s' (%s) added to frameset", 
                         product_filename, UVES_EFFICIENCY_TABLE(chip));
                } /* end if calc_response */
                
                /* Save background image */
                cpl_free(product_filename);
                check( product_filename = 
                   uves_response_bkg_standard_filename(chip), 
                   "Error getting filename");
                uves_pfits_set_extname(rotated_header[raw_index],"background");
                check( uves_frameset_insert(frames,
                            background,
                            CPL_FRAME_GROUP_PRODUCT,
                            CPL_FRAME_TYPE_IMAGE,
                            CPL_FRAME_LEVEL_INTERMEDIATE,
                            product_filename,
                            UVES_BKG_STD(chip),
                            raw_header[raw_index],
                            rotated_header[raw_index],
                            NULL,
                            parameters,
                            make_str(UVES_RESPONSE_ID),
                            PACKAGE "/" PACKAGE_VERSION, NULL,
                            starttime, false, 
                            CPL_STATS_MIN | CPL_STATS_MAX),
                   "Could not add background image '%s' (%s) to frameset",
                   product_filename, UVES_BKG_STD(chip));

                uves_msg("Background image '%s' (%s) added to frameset",
                     product_filename, UVES_BKG_STD(chip));
                

                /* Save info table */
               cpl_free(product_filename);
                check( product_filename = 
                   uves_order_extract_qc_standard_filename(chip), 
                   "Error getting filename");
                uves_pfits_set_extname(rotated_header[raw_index],"quality control table");
                sprintf(extname,"QC_INFO");
                uves_pfits_set_extname(table_header,extname);

                check( uves_frameset_insert(frames,
                            info_tbl,
                            CPL_FRAME_GROUP_PRODUCT,
                            CPL_FRAME_TYPE_TABLE,
                            CPL_FRAME_LEVEL_INTERMEDIATE,
                            product_filename,
                            UVES_ORDER_EXTRACT_QC(chip),
                            raw_header[raw_index],
                            rotated_header[raw_index],
                            table_header,
                            parameters,
                            make_str(UVES_RESPONSE_ID),
                            PACKAGE "/" PACKAGE_VERSION, NULL,
                            starttime, true, 
                            0),
                   "Could not add extraction quality table '%s' (%s) to frameset",
                   product_filename, UVES_ORDER_EXTRACT_QC(chip));

                uves_msg("Extraction quality table '%s' (%s) added to frameset",
                     product_filename, UVES_ORDER_EXTRACT_QC(chip));
                


            }/* ... if trace is enabled */
            else
            {
                uves_msg_low("Skipping trace number %d", trace_number);
            }
        
        }/* For each trace */
    

      if(strcmp(PROCESS_CHIP,"REDL") == 0) {
	chip = uves_chip_get_next(chip);
      }

      uves_free_propertylist(&table_header);



    } /* For each chip */
    
  cleanup:
    /* Input */
    uves_free_image(&raw_image[0]);
    uves_free_image(&raw_image[1]);
    uves_free_propertylist(&raw_header[0]);
    uves_free_propertylist(&raw_header[1]);
    uves_free_propertylist(&rotated_header[0]);
    uves_free_propertylist(&rotated_header[1]);
    
    /* Input, calib */
    uves_free_image(&master_bias);
    uves_free_propertylist(&master_bias_header);
    
    uves_free_image(&master_dark);
    uves_free_propertylist(&master_dark_header);

    uves_free_image(&master_flat);
    uves_free_propertylist(&master_flat_header);
    
    uves_free_table(&ordertable);
    uves_free_propertylist(&ordertable_header);
    uves_polynomial_delete(&order_locations);
    uves_free_table(&traces);
    uves_free_propertylist(&table_header);
    uves_free_table_const( &(linetable[0]) );
    uves_free_table_const( &(linetable[1]) );
    uves_free_table_const( &(linetable[2]) );
    uves_free_propertylist_const( &(linetable_header[0]) );
    uves_free_propertylist_const( &(linetable_header[1]) );
    uves_free_propertylist_const( &(linetable_header[2]) );
    uves_polynomial_delete_const( &(dispersion_relation[0]) );
    uves_polynomial_delete_const( &(dispersion_relation[1]) );
    uves_polynomial_delete_const( &(dispersion_relation[2]) );
    uves_free_table( &flux_table );
    uves_free_table( &atm_extinction );
    
    /* Output */
    uves_qclog_delete(&qclog[0]);
    uves_qclog_delete(&qclog_optext[0]);
    uves_free_image(&background);
    uves_free_image(&reduced_spectrum);
    uves_free_propertylist(&spectrum_header);
    uves_free_propertylist(&response_header_2d);
    uves_free_propertylist(&response_header);
    uves_free_propertylist(&efficiency_header);
    uves_free_image(&response_orders);
    uves_free_image(&response_curve);
    uves_free_table(&efficiency);
    uves_free_table(&blaze_efficiency);
    uves_free_table(&info_tbl);
    
    /* Local */
    uves_free_table(&catalogue_flux);    
    cpl_free(product_filename);
    cpl_free(ref_obj_name);
    
    return;
}

/**
@@ brief QC log generation
@param table              efficiency table
@param info_tbl           table with info on object location during extraction
@param raw_header         FITS header
@param chip               chip ID
@param qclog              ouput QC log table
@param ref_obj_name       reference object name
*/
static void uves_efficiency_qclog(cpl_table* table,
                  uves_propertylist* raw_header, 
                  enum uves_chip chip,
                  cpl_table* qclog,
                  const char *ref_obj_name)
{
  int i=0;
  bool new_format;
  
  check( new_format = uves_format_is_new(raw_header),
     "Error determining FITS header format");

  check_nomsg(uves_qclog_add_string(qclog,
                        "QC TEST1 ID",
                        "Efficiency-Test-Results",
                        "Name of QC test",
                    "%s"));

  check_nomsg(uves_qclog_add_string(qclog,
                        uves_remove_string_prefix(UVES_INSPATH,"ESO "),
                         uves_pfits_get_inspath(raw_header),
                        "Optical path used.",
                   "%s"));

  check_nomsg(uves_qclog_add_string(qclog,
                        uves_remove_string_prefix(UVES_INSMODE,"ESO "),
                         uves_pfits_get_insmode(raw_header),
                        "Instrument mode used.",
                   "%s"));

  check_nomsg(uves_qclog_add_string(qclog,
                    uves_remove_string_prefix(UVES_GRATNAME(chip),"ESO "),
                    uves_pfits_get_gratname(raw_header,chip),
                    "Cross disperser ID",
                    "%s"));

  check_nomsg(uves_qclog_add_common_wave(raw_header,
                                         chip, qclog));

  check_nomsg(
      uves_qclog_add_string(qclog,
                uves_remove_string_prefix(UVES_CHIP_NAME(chip),"ESO "),
                /* UVES_QC_CHIP_VAL(chip), */
                uves_pfits_get_chip_name(raw_header, chip),
                "Detector chip name",
                "%s"));
  
  check_nomsg(
      uves_qclog_add_double(qclog,
                uves_remove_string_prefix(UVES_GRATWLEN(chip),"ESO "),
                uves_pfits_get_gratwlen(raw_header,chip),
                "Grating central wavelength [nm] (hs).",
                "%.1f"));
  
  check_nomsg(
      uves_qclog_add_double(qclog,
                uves_remove_string_prefix(UVES_CONAD(new_format, chip),"ESO "),
                uves_pfits_get_conad(raw_header, chip),
                "Conversion from ADUs to electrons",
                "%8.2f"));
  
 
  check_nomsg(
      uves_qclog_add_double(qclog,
                uves_remove_string_prefix(UVES_QC_UIT(new_format, chip), "ESO "),
                uves_pfits_get_uit(raw_header),
                "user defined subintegration time",
                "%8.0f"));
  
  check_nomsg(
      uves_qclog_add_double(qclog,
                "AIRMASS",
                (uves_pfits_get_airmass_start(raw_header) +
                 uves_pfits_get_airmass_end(raw_header))/2.0,
                "Averaged airmass",
                "%8.4f"));

  check_nomsg(
      uves_qclog_add_string(qclog,
                uves_remove_string_prefix(UVES_TARG_NAME, "ESO "),
                ref_obj_name,
                "OB target name",
                "%s"));
  
 
  for(i = 0; i < cpl_table_get_nrow(table); i++) 
      {
      char key_name[25];

      int order = cpl_table_get_int(table, "Order", i, NULL);

      sprintf(key_name,"QC BLAZEFF%d", order);

      check_nomsg(uves_qclog_add_double(qclog,
                        key_name,
                        cpl_table_get_double(table, "Eff", i, NULL),
                        "Blaze Efficiency",
                        "%13.6f"));
      /*
        uves_msg("QC-LOG: Wlen =%g Eff=%g", 
        cpl_table_get_double(table,"Wave",i,NULL),
        cpl_table_get_double(table,"Eff",i,NULL));
      */

      sprintf(key_name,"QC BLAZWLEN%d", order);
      check_nomsg(uves_qclog_add_double(qclog,
                        key_name,
                        cpl_table_get_double(table, "Wave", i, NULL)/10.,
                        "Blaze wavelength",  /* nm */
                        "%13.6f"));
      }

  /* Are these QC parameters needed anywhere? */
#if 0
  for(i = 0; i < cpl_table_get_nrow(info_tbl); i++) 
      {
      char key_name[25];

      int order = cpl_table_get_int(info_tbl, "Order", i, NULL);

      sprintf(key_name,"QC ORDER NUM%d", order);
      check_nomsg(uves_qclog_add_int(qclog,
                     key_name,
                     order,
                     "Order Number",
                     "%d"));
      /*
        uves_msg("QC-LOG: Order =%d S/N=%g", 
        cpl_table_get_int(info_tbl,"Order",i,NULL),
        cpl_table_get_double(info_tbl,"S/N",i,NULL));
      */

      sprintf(key_name,"QC OBJ SN%d", order);
      check_nomsg(uves_qclog_add_double(qclog,
                        key_name,
                        cpl_table_get_double(info_tbl,"S/N",i,NULL),
                        "Order S/N",
                        "%f13.6"));
      }
#endif

 cleanup:
  return;
}

/**@}*/