fors_calib.c

00001 /* $Id: fors_calib.c,v 1.11 2008/09/03 12:41:43 cizzo Exp $
00002  *
00003  * This file is part of the FORS Data Reduction Pipeline
00004  * Copyright (C) 2006 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2008/09/03 12:41:43 $
00024  * $Revision: 1.11 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <math.h>
00033 #include <cpl.h>
00034 #include <moses.h>
00035 #include <fors_dfs.h>
00036 #include <fors_qc.h>
00037 
00038 static int fors_calib_create(cpl_plugin *);
00039 static int fors_calib_exec(cpl_plugin *);
00040 static int fors_calib_destroy(cpl_plugin *);
00041 static int fors_calib(cpl_parameterlist *, cpl_frameset *);
00042 
00043 static char fors_calib_description[] =
00044 "This recipe is used to identify reference lines on LSS, MOS and MXU arc lamp\n"
00045 "exposures, and trace the spectral edges on the corresponding flat field\n"
00046 "exposures. This information is used to determine the spectral extraction\n"
00047 "mask to be applied in the scientific data reduction, performed with the\n"
00048 "recipe fors_science.\n"
00049 "This recipe accepts both FORS1 and FORS2 frames. The input arc lamp and\n"
00050 "flat field exposures are assumed to be obtained quasi-simultaneously, so\n"
00051 "that they would be described by exactly the same instrument distortions.\n"
00052 "A line catalog must be specified, containing the wavelengths of the\n"
00053 "reference arc lamp lines used for the wavelength calibration. A grism\n"
00054 "table (typically depending on the instrument mode, and in particular on\n"
00055 "the grism used) may also be specified: this table contains a default\n"
00056 "recipe parameter setting to control the way spectra are extracted for\n"
00057 "a specific instrument mode, as it is used for automatic run of the\n"
00058 "pipeline on Paranal and in Garching. If this table is specified, it\n"
00059 "will modify the default recipe parameter setting, with the exception of\n"
00060 "those parameters which have been explicitly modifyed on the command line.\n"
00061 "If a grism table is not specified, the input recipe parameters values\n"
00062 "will always be read from the command line, or from an esorex configuration\n"
00063 "file if present, or from their generic default values (that are rarely\n"
00064 "meaningful). Finally a master bias frame must be input to this recipe.\n" 
00065 "In the table below the MXU acronym can be read alternatively as MOS\n"
00066 "and LSS, with the exception of CURV_COEFF_LSS, CURV_TRACES_LSS,\n"
00067 "SPATIAL_MAP_LSS, SPECTRA_DETECTION_LSS, and and SLIT_MAP_LSS, which are\n" 
00068 "never created. The products SPECTRA_DETECTION_MXU, SLIT_MAP_MXU, and\n" 
00069 "DISP_RESIDUALS_MXU, are just created if the --check parameter is set to\n" 
00070 "true. The product GLOBAL_DISTORTION_TABLE is just created if more than 12\n" 
00071 "separate spectra are found in the CCD.\n\n"
00072 "Input files:\n\n"
00073 "  DO category:             Type:       Explanation:         Required:\n"
00074 "  SCREEN_FLAT_MXU          Raw         Flat field exposures    Y\n"
00075 "  LAMP_MXU                 Raw         Arc lamp exposure       Y\n"
00076 "  MASTER_BIAS or BIAS      Calib       Bias frame              Y\n"
00077 "  MASTER_LINECAT           Calib       Line catalog            Y\n"
00078 "  GRISM_TABLE              Calib       Grism table             .\n\n"
00079 "Output files:\n\n"
00080 "  DO category:             Data type:  Explanation:\n"
00081 "  MASTER_SCREEN_FLAT_MXU   FITS image  Combined (sum) flat field\n"
00082 "  MASTER_NORM_FLAT_MXU     FITS image  Normalised flat field\n"
00083 "  MAPPED_SCREEN_FLAT_MXU   FITS image  Wavelength calibrated flat field\n"
00084 "  MAPPED_NORM_FLAT_MXU     FITS image  Wavelength calibrated normalised flat\n"
00085 "  REDUCED_LAMP_MXU         FITS image  Wavelength calibrated arc spectrum\n"
00086 "  DISP_COEFF_MXU           FITS table  Inverse dispersion coefficients\n"
00087 "  DISP_RESIDUALS_MXU       FITS image  Residuals in wavelength calibration\n"
00088 "  DISP_RESIDUALS_TABLE_MXU FITS table  Residuals in wavelength calibration\n"
00089 "  DELTA_IMAGE_MXU          FITS image  Offset vs linear wavelength calib\n"
00090 "  WAVELENGTH_MAP_MXU       FITS image  Wavelength for each pixel on CCD\n"
00091 "  SPECTRA_DETECTION_MXU    FITS image  Check for preliminary detection\n"
00092 "  SLIT_MAP_MXU             FITS image  Map of central wavelength on CCD\n"
00093 "  CURV_TRACES_MXU          FITS table  Spectral curvature traces\n"
00094 "  CURV_COEFF_MXU           FITS table  Spectral curvature coefficients\n"
00095 "  SPATIAL_MAP_MXU          FITS image  Spatial position along slit on CCD\n"
00096 "  SPECTRAL_RESOLUTION_MXU  FITS table  Resolution at reference arc lines\n"
00097 "  SLIT_LOCATION_MXU        FITS table  Slits on product frames and CCD\n"
00098 "  GLOBAL_DISTORTION_TABLE  FITS table  Global distortions table\n\n";
00099 
00100 #define fors_calib_exit(message)              \
00101 {                                             \
00102 if (message) cpl_msg_error(recipe, message);  \
00103 cpl_free(instrume);                           \
00104 cpl_free(pipefile);                           \
00105 cpl_free(fiterror);                           \
00106 cpl_free(fitlines);                           \
00107 cpl_image_delete(bias);                       \
00108 cpl_image_delete(master_bias);                \
00109 cpl_image_delete(coordinate);                 \
00110 cpl_image_delete(checkwave);                  \
00111 cpl_image_delete(flat);                       \
00112 cpl_image_delete(master_flat);                \
00113 cpl_image_delete(norm_flat);                  \
00114 cpl_image_delete(mapped_flat);                \
00115 cpl_image_delete(mapped_nflat);               \
00116 cpl_image_delete(rainbow);                    \
00117 cpl_image_delete(rectified);                  \
00118 cpl_image_delete(residual);                   \
00119 cpl_image_delete(smo_flat);                   \
00120 cpl_image_delete(spatial);                    \
00121 cpl_image_delete(spectra);                    \
00122 cpl_image_delete(wavemap);                    \
00123 cpl_image_delete(delta);                      \
00124 cpl_image_delete(rect_flat);                  \
00125 cpl_image_delete(rect_nflat);                 \
00126 cpl_image_delete(mapped_flat);                \
00127 cpl_image_delete(mapped_nflat);               \
00128 cpl_mask_delete(refmask);                     \
00129 cpl_propertylist_delete(header);              \
00130 cpl_propertylist_delete(save_header);         \
00131 cpl_propertylist_delete(qclist);              \
00132 cpl_table_delete(grism_table);                \
00133 cpl_table_delete(idscoeff);                   \
00134 cpl_table_delete(idscoeff_all);               \
00135 cpl_table_delete(restable);                   \
00136 cpl_table_delete(maskslits);                  \
00137 cpl_table_delete(overscans);                  \
00138 cpl_table_delete(traces);                     \
00139 cpl_table_delete(polytraces);                 \
00140 cpl_table_delete(slits);                      \
00141 cpl_table_delete(restab);                     \
00142 cpl_table_delete(global);                     \
00143 cpl_table_delete(wavelengths);                \
00144 cpl_vector_delete(lines);                     \
00145 cpl_msg_indent_less();                        \
00146 return -1;                                    \
00147 }
00148 
00149 #define fors_calib_exit_memcheck(message)       \
00150 {                                               \
00151 if (message) cpl_msg_info(recipe, message);     \
00152 printf("free instrume (%p)\n", instrume);       \
00153 cpl_free(instrume);                             \
00154 printf("free pipefile (%p)\n", pipefile);       \
00155 cpl_free(pipefile);                             \
00156 printf("free fiterror (%p)\n", fiterror);       \
00157 cpl_free(fiterror);                             \
00158 printf("free fitlines (%p)\n", fitlines);       \
00159 cpl_free(fitlines);                             \
00160 printf("free bias (%p)\n", bias);               \
00161 cpl_image_delete(bias);                         \
00162 printf("free master_bias (%p)\n", master_bias); \
00163 cpl_image_delete(master_bias);                  \
00164 printf("free coordinate (%p)\n", coordinate);   \
00165 cpl_image_delete(coordinate);                   \
00166 printf("free checkwave (%p)\n", checkwave);     \
00167 cpl_image_delete(checkwave);                    \
00168 printf("free flat (%p)\n", flat);               \
00169 cpl_image_delete(flat);                         \
00170 printf("free master_flat (%p)\n", master_flat); \
00171 cpl_image_delete(master_flat);                  \
00172 printf("free norm_flat (%p)\n", norm_flat);     \
00173 cpl_image_delete(norm_flat);                    \
00174 printf("free mapped_flat (%p)\n", mapped_flat); \
00175 cpl_image_delete(mapped_flat);                  \
00176 printf("free mapped_nflat (%p)\n", mapped_nflat); \
00177 cpl_image_delete(mapped_nflat);                 \
00178 printf("free rainbow (%p)\n", rainbow);         \
00179 cpl_image_delete(rainbow);                      \
00180 printf("free rectified (%p)\n", rectified);     \
00181 cpl_image_delete(rectified);                    \
00182 printf("free residual (%p)\n", residual);       \
00183 cpl_image_delete(residual);                     \
00184 printf("free smo_flat (%p)\n", smo_flat);       \
00185 cpl_image_delete(smo_flat);                     \
00186 printf("free spatial (%p)\n", spatial);         \
00187 cpl_image_delete(spatial);                      \
00188 printf("free spectra (%p)\n", spectra);         \
00189 cpl_image_delete(spectra);                      \
00190 printf("free wavemap (%p)\n", wavemap);         \
00191 cpl_image_delete(wavemap);                      \
00192 printf("free delta (%p)\n", delta);             \
00193 cpl_image_delete(delta);                        \
00194 printf("free rect_flat (%p)\n", rect_flat);     \
00195 cpl_image_delete(rect_flat);                    \
00196 printf("free rect_nflat (%p)\n", rect_nflat);   \
00197 cpl_image_delete(rect_nflat);                   \
00198 printf("free refmask (%p)\n", refmask);         \
00199 cpl_mask_delete(refmask);                       \
00200 printf("free header (%p)\n", header);           \
00201 cpl_propertylist_delete(header);                \
00202 printf("free save_header (%p)\n", save_header); \
00203 cpl_propertylist_delete(save_header);           \
00204 printf("free qclist (%p)\n", qclist);           \
00205 cpl_propertylist_delete(qclist);                \
00206 printf("free grism_table (%p)\n", grism_table); \
00207 cpl_table_delete(grism_table);                  \
00208 printf("free idscoeff (%p)\n", idscoeff);       \
00209 cpl_table_delete(idscoeff);                     \
00210 printf("free idscoeff_all (%p)\n", idscoeff_all);  \
00211 cpl_table_delete(idscoeff_all);                 \
00212 printf("free restable (%p)\n", restable);       \
00213 cpl_table_delete(restable);                     \
00214 printf("free maskslits (%p)\n", maskslits);     \
00215 cpl_table_delete(maskslits);                    \
00216 printf("free overscans (%p)\n", overscans);     \
00217 cpl_table_delete(overscans);                    \
00218 printf("free traces (%p)\n", traces);           \
00219 cpl_table_delete(traces);                       \
00220 printf("free polytraces (%p)\n", polytraces);   \
00221 cpl_table_delete(polytraces);                   \
00222 printf("free slits (%p)\n", slits);             \
00223 cpl_table_delete(slits);                        \
00224 printf("free restab (%p)\n", restab);           \
00225 cpl_table_delete(restab);                       \
00226 printf("free global (%p)\n", global);           \
00227 cpl_table_delete(global);                       \
00228 printf("free wavelengths (%p)\n", wavelengths); \
00229 cpl_table_delete(wavelengths);                  \
00230 printf("free lines (%p)\n", lines);             \
00231 cpl_vector_delete(lines);                       \
00232 cpl_msg_indent_less();                          \
00233 return 0;                                       \
00234 }
00235 
00236 
00248 int cpl_plugin_get_info(cpl_pluginlist *list)
00249 {
00250     cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe );
00251     cpl_plugin *plugin = &recipe->interface;
00252 
00253     cpl_plugin_init(plugin,
00254                     CPL_PLUGIN_API,
00255                     FORS_BINARY_VERSION,
00256                     CPL_PLUGIN_TYPE_RECIPE,
00257                     "fors_calib",
00258                     "Determination of the extraction mask",
00259                     fors_calib_description,
00260                     "Carlo Izzo",
00261                     PACKAGE_BUGREPORT,
00262     "This file is currently part of the FORS Instrument Pipeline\n"
00263     "Copyright (C) 2002-2006 European Southern Observatory\n\n"
00264     "This program is free software; you can redistribute it and/or modify\n"
00265     "it under the terms of the GNU General Public License as published by\n"
00266     "the Free Software Foundation; either version 2 of the License, or\n"
00267     "(at your option) any later version.\n\n"
00268     "This program is distributed in the hope that it will be useful,\n"
00269     "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
00270     "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
00271     "GNU General Public License for more details.\n\n"
00272     "You should have received a copy of the GNU General Public License\n"
00273     "along with this program; if not, write to the Free Software Foundation,\n"
00274     "Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n",
00275                     fors_calib_create,
00276                     fors_calib_exec,
00277                     fors_calib_destroy);
00278 
00279     cpl_pluginlist_append(list, plugin);
00280     
00281     return 0;
00282 }
00283 
00284 
00295 static int fors_calib_create(cpl_plugin *plugin)
00296 {
00297     cpl_recipe    *recipe;
00298     cpl_parameter *p;
00299 
00300 
00301     /* 
00302      * Check that the plugin is part of a valid recipe 
00303      */
00304 
00305     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00306         recipe = (cpl_recipe *)plugin;
00307     else 
00308         return -1;
00309 
00310     /* 
00311      * Create the parameters list in the cpl_recipe object 
00312      */
00313 
00314     recipe->parameters = cpl_parameterlist_new(); 
00315 
00316 
00317     /*
00318      * Dispersion
00319      */
00320 
00321     p = cpl_parameter_new_value("fors.fors_calib.dispersion",
00322                                 CPL_TYPE_DOUBLE,
00323                                 "Expected spectral dispersion (Angstrom/pixel)",
00324                                 "fors.fors_calib",
00325                                 0.0);
00326     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dispersion");
00327     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00328     cpl_parameterlist_append(recipe->parameters, p);
00329 
00330     /*
00331      * Peak detection level
00332      */
00333 
00334     p = cpl_parameter_new_value("fors.fors_calib.peakdetection",
00335                                 CPL_TYPE_DOUBLE,
00336                                 "Initial peak detection threshold (ADU)",
00337                                 "fors.fors_calib",
00338                                 0.0);
00339     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "peakdetection");
00340     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00341     cpl_parameterlist_append(recipe->parameters, p);
00342 
00343     /* 
00344      * Degree of wavelength calibration polynomial
00345      */
00346 
00347     p = cpl_parameter_new_value("fors.fors_calib.wdegree",
00348                                 CPL_TYPE_INT,
00349                                 "Degree of wavelength calibration polynomial",
00350                                 "fors.fors_calib",
00351                                 0);
00352     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wdegree");
00353     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00354     cpl_parameterlist_append(recipe->parameters, p);
00355 
00356     /*
00357      * Reference lines search radius
00358      */
00359 
00360     p = cpl_parameter_new_value("fors.fors_calib.wradius",
00361                                 CPL_TYPE_INT,
00362                                 "Search radius if iterating pattern-matching "
00363                                 "with first-guess method",
00364                                 "fors.fors_calib",
00365                                 4);
00366     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wradius");
00367     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00368     cpl_parameterlist_append(recipe->parameters, p);
00369 
00370     /*
00371      * Rejection threshold in dispersion relation polynomial fitting
00372      */
00373 
00374     p = cpl_parameter_new_value("fors.fors_calib.wreject",
00375                                 CPL_TYPE_DOUBLE,
00376                                 "Rejection threshold in dispersion "
00377                                 "relation fit (pixel)",
00378                                 "fors.fors_calib",
00379                                 0.7);
00380     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wreject");
00381     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00382     cpl_parameterlist_append(recipe->parameters, p);
00383 
00384     /*
00385      * Wavelength solution interpolation (for LSS data)
00386      */
00387 
00388     p = cpl_parameter_new_value("fors.fors_calib.wmode",
00389                                 CPL_TYPE_INT,
00390                                 "Interpolation mode of wavelength solution "
00391                                 "applicable to LSS-like data (0 = no "
00392                                 "interpolation, 1 = fill gaps, 2 = global "
00393                                 "model)",
00394                                 "fors.fors_calib",
00395                                 2);
00396     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wmode");
00397     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00398     cpl_parameterlist_append(recipe->parameters, p);
00399 
00400     /*
00401      * Line catalog table column containing the reference wavelengths
00402      */
00403 
00404     p = cpl_parameter_new_value("fors.fors_calib.wcolumn",
00405                                 CPL_TYPE_STRING,
00406                                 "Name of line catalog table column "
00407                                 "with wavelengths",
00408                                 "fors.fors_calib",
00409                                 "WLEN");
00410     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wcolumn");
00411     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00412     cpl_parameterlist_append(recipe->parameters, p);
00413 
00414     /*
00415      * Degree of spectral curvature polynomial
00416      */
00417 
00418     p = cpl_parameter_new_value("fors.fors_calib.cdegree",
00419                                 CPL_TYPE_INT,
00420                                 "Degree of spectral curvature polynomial",
00421                                 "fors.fors_calib",
00422                                 0);
00423     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cdegree");
00424     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00425     cpl_parameterlist_append(recipe->parameters, p);
00426 
00427     /*
00428      * Curvature solution interpolation (for MOS-like data)
00429      */
00430  
00431     p = cpl_parameter_new_value("fors.fors_calib.cmode",
00432                                 CPL_TYPE_INT,
00433                                 "Interpolation mode of curvature solution "
00434                                 "applicable to MOS-like data (0 = no "
00435                                 "interpolation, 1 = fill gaps, 2 = global "
00436                                 "model)",
00437                                 "fors.fors_calib",
00438                                 1);
00439     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cmode");
00440     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00441     cpl_parameterlist_append(recipe->parameters, p);
00442 
00443     /*
00444      * Start wavelength for spectral extraction
00445      */
00446 
00447     p = cpl_parameter_new_value("fors.fors_calib.startwavelength",
00448                                 CPL_TYPE_DOUBLE,
00449                                 "Start wavelength in spectral extraction",
00450                                 "fors.fors_calib",
00451                                 0.0);
00452     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "startwavelength");
00453     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00454     cpl_parameterlist_append(recipe->parameters, p);
00455 
00456     /*
00457      * End wavelength for spectral extraction
00458      */
00459 
00460     p = cpl_parameter_new_value("fors.fors_calib.endwavelength",
00461                                 CPL_TYPE_DOUBLE,
00462                                 "End wavelength in spectral extraction",
00463                                 "fors.fors_calib",
00464                                 0.0);
00465     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "endwavelength");
00466     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00467     cpl_parameterlist_append(recipe->parameters, p);
00468 
00469     /*
00470      * Try slit identification
00471      */
00472 
00473     p = cpl_parameter_new_value("fors.fors_calib.slit_ident",
00474                                 CPL_TYPE_BOOL,
00475                                 "Attempt slit identification for MOS or MXU",
00476                                 "fors.fors_calib",
00477                                 TRUE);
00478     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "slit_ident");
00479     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00480     cpl_parameterlist_append(recipe->parameters, p);
00481 
00482     /*
00483      * Degree of flat field fitting polynomial along spatial direction 
00484      * (used for LSS data)
00485      */
00486 
00487     p = cpl_parameter_new_value("fors.fors_calib.sdegree",
00488                                 CPL_TYPE_INT,
00489                                 "Degree of flat field fitting polynomial "
00490                                 "along spatial direction (used for LSS "
00491                                 "data only)",
00492                                 "fors.fors_calib",
00493                                 4);
00494     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "sdegree");
00495     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00496     cpl_parameterlist_append(recipe->parameters, p);
00497 
00498     /*
00499      * Degree of flat field fitting polynomial along dispersion direction
00500      * (used for MOS and MXU data)
00501      */
00502 
00503     p = cpl_parameter_new_value("fors.fors_calib.ddegree",
00504                                 CPL_TYPE_INT,
00505                                 "Degree of flat field fitting polynomial "
00506                                 "along dispersion direction (used for MOS "
00507                                 "and MXU data only)",
00508                                 "fors.fors_calib",
00509                                 -1);
00510     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ddegree");
00511     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00512     cpl_parameterlist_append(recipe->parameters, p);
00513 
00514     /*
00515      * Smooth box radius for flat field along dispersion direction
00516      */
00517 
00518     p = cpl_parameter_new_value("fors.fors_calib.dradius",
00519                                 CPL_TYPE_INT,
00520                                 "Smooth box radius for flat field along "
00521                                 "dispersion direction",
00522                                 "fors.fors_calib",
00523                                 10);
00524     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dradius");
00525     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00526     cpl_parameterlist_append(recipe->parameters, p);
00527 
00528     /*
00529      * Smooth box radius for flat field along spatial direction
00530      * (used for LSS data only)
00531      */
00532 
00533     p = cpl_parameter_new_value("fors.fors_calib.sradius",
00534                                 CPL_TYPE_INT,
00535                                 "Smooth box radius for flat field along "
00536                                 "spatial direction",
00537                                 "fors.fors_calib",
00538                                 10);
00539     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "sradius");
00540     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00541     cpl_parameterlist_append(recipe->parameters, p);
00542 
00543     /*
00544      * Computation of QC1 parameters
00545      */
00546 
00547     p = cpl_parameter_new_value("fors.fors_calib.qc",
00548                                 CPL_TYPE_BOOL,
00549                                 "Compute QC1 parameters",
00550                                 "fors.fors_calib",
00551                                 TRUE);
00552     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "qc");
00553     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00554     cpl_parameterlist_append(recipe->parameters, p);
00555 
00556     /*
00557      * Create check products
00558      */
00559 
00560     p = cpl_parameter_new_value("fors.fors_calib.check",
00561                                 CPL_TYPE_BOOL,
00562                                 "Create intermediate products",
00563                                 "fors.fors_calib",
00564                                 FALSE);
00565     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "check");
00566     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00567     cpl_parameterlist_append(recipe->parameters, p);
00568 
00569     return 0;
00570 }
00571 
00572 
00581 static int fors_calib_exec(cpl_plugin *plugin)
00582 {
00583     cpl_recipe *recipe;
00584     
00585     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00586         recipe = (cpl_recipe *)plugin;
00587     else 
00588         return -1;
00589 
00590     return fors_calib(recipe->parameters, recipe->frames);
00591 }
00592 
00593 
00602 static int fors_calib_destroy(cpl_plugin *plugin)
00603 {
00604     cpl_recipe *recipe;
00605     
00606     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00607         recipe = (cpl_recipe *)plugin;
00608     else 
00609         return -1;
00610 
00611     cpl_parameterlist_delete(recipe->parameters); 
00612 
00613     return 0;
00614 }
00615 
00616 
00626 static int fors_calib(cpl_parameterlist *parlist, cpl_frameset *frameset)
00627 {
00628 
00629     const char *recipe = "fors_calib";
00630 
00631 
00632     /*
00633      * Input parameters
00634      */
00635 
00636     double      dispersion;
00637     double      peakdetection;
00638     int         wdegree;
00639     int         wradius;
00640     double      wreject;
00641     int         wmode;
00642     const char *wcolumn;
00643     int         cdegree;
00644     int         cmode;
00645     double      startwavelength;
00646     double      endwavelength;
00647     int         slit_ident;
00648     int         sdegree;
00649     int         ddegree;
00650     int         sradius;
00651     int         dradius;
00652     int         qc;
00653     int         check;
00654 
00655     /*
00656      * CPL objects
00657      */
00658 
00659     cpl_imagelist    *biases       = NULL;
00660     cpl_image        *bias         = NULL;
00661     cpl_image        *master_bias  = NULL;
00662     cpl_image        *multi_bias   = NULL;
00663     cpl_image        *flat         = NULL;
00664     cpl_image        *master_flat  = NULL;
00665     cpl_image        *smo_flat     = NULL;
00666     cpl_image        *norm_flat    = NULL;
00667     cpl_image        *spectra      = NULL;
00668     cpl_image        *wavemap      = NULL;
00669     cpl_image        *delta        = NULL;
00670     cpl_image        *residual     = NULL;
00671     cpl_image        *checkwave    = NULL;
00672     cpl_image        *rectified    = NULL;
00673     cpl_image        *dummy        = NULL;
00674     cpl_image        *refimage     = NULL;
00675     cpl_image        *coordinate   = NULL;
00676     cpl_image        *rainbow      = NULL;
00677     cpl_image        *spatial      = NULL;
00678     cpl_image        *rect_flat    = NULL;
00679     cpl_image        *rect_nflat   = NULL;
00680     cpl_image        *mapped_flat  = NULL;
00681     cpl_image        *mapped_nflat = NULL;
00682 
00683     cpl_mask         *refmask      = NULL;
00684 
00685     cpl_table        *grism_table  = NULL;
00686     cpl_table        *overscans    = NULL;
00687     cpl_table        *wavelengths  = NULL;
00688     cpl_table        *idscoeff     = NULL;
00689     cpl_table        *idscoeff_all = NULL;
00690     cpl_table        *restable     = NULL;
00691     cpl_table        *slits        = NULL;
00692     cpl_table        *positions    = NULL;
00693     cpl_table        *maskslits    = NULL;
00694     cpl_table        *traces       = NULL;
00695     cpl_table        *polytraces   = NULL;
00696     cpl_table        *restab       = NULL;
00697     cpl_table        *global       = NULL;
00698 
00699     cpl_vector       *lines        = NULL;
00700 
00701     cpl_propertylist *header       = NULL;
00702     cpl_propertylist *save_header  = NULL;
00703     cpl_propertylist *qclist       = NULL;
00704 
00705     /*
00706      * Auxiliary variables
00707      */
00708 
00709     char    version[80];
00710     const char   *arc_tag;
00711     const char   *flat_tag;
00712     const char   *master_screen_flat_tag;
00713     const char   *master_norm_flat_tag;
00714     const char   *reduced_lamp_tag;
00715     const char   *disp_residuals_tag;
00716     const char   *disp_coeff_tag;
00717     const char   *wavelength_map_tag;
00718     const char   *spectra_detection_tag;
00719     const char   *spectral_resolution_tag;
00720     const char   *slit_map_tag;
00721     const char   *curv_traces_tag;
00722     const char   *curv_coeff_tag;
00723     const char   *spatial_map_tag;
00724     const char   *slit_location_tag;
00725     const char   *global_distortion_tag = "GLOBAL_DISTORTION_TABLE";
00726     const char   *disp_residuals_table_tag;
00727     const char   *delta_image_tag;
00728     const char   *mapped_screen_flat_tag;
00729     const char   *mapped_norm_flat_tag;
00730     const char   *keyname;
00731     int     mxu, mos, lss;
00732     int     treat_as_lss = 0;
00733     int     nslits;
00734     float  *data;
00735     double *xpos;
00736     double  mxpos;
00737     double  mean_rms;
00738     double  mean_rms_err;
00739     double  alltime;
00740     int     nflats;
00741     int     nbias;
00742     int     nlines;
00743     int     rebin;
00744     double *line;
00745     double *fiterror = NULL;
00746     int    *fitlines = NULL;
00747     int     nx, ny;
00748     double  reference;
00749     double  gain;
00750     int     compute_central_wave;
00751     int     ccd_xsize, ccd_ysize;
00752     int     i;
00753 
00754     char   *instrume = NULL;
00755     char   *pipefile = NULL;
00756 
00757 
00758     snprintf(version, 80, "%s-%s", PACKAGE, PACKAGE_VERSION);
00759 
00760     cpl_msg_set_indentation(2);
00761 
00762     if (dfs_files_dont_exist(frameset))
00763         fors_calib_exit(NULL);
00764 
00765     /* 
00766      * Get configuration parameters
00767      */
00768 
00769     cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
00770     cpl_msg_indent_more();
00771 
00772     if (cpl_frameset_count_tags(frameset, "GRISM_TABLE") > 1)
00773         fors_calib_exit("Too many in input: GRISM_TABLE");
00774 
00775     grism_table = dfs_load_table(frameset, "GRISM_TABLE", 1);
00776 
00777     dispersion = dfs_get_parameter_double(parlist, 
00778                     "fors.fors_calib.dispersion", grism_table);
00779 
00780     if (dispersion <= 0.0)
00781         fors_calib_exit("Invalid spectral dispersion value");
00782 
00783     peakdetection = dfs_get_parameter_double(parlist, 
00784                     "fors.fors_calib.peakdetection", grism_table);
00785     if (peakdetection <= 0.0)
00786         fors_calib_exit("Invalid peak detection level");
00787 
00788     wdegree = dfs_get_parameter_int(parlist, 
00789                     "fors.fors_calib.wdegree", grism_table);
00790 
00791     if (wdegree < 1)
00792         fors_calib_exit("Invalid polynomial degree");
00793 
00794     if (wdegree > 5)
00795         fors_calib_exit("Max allowed polynomial degree is 5");
00796 
00797     wradius = dfs_get_parameter_int(parlist, "fors.fors_calib.wradius", NULL);
00798 
00799     if (wradius < 0)
00800         fors_calib_exit("Invalid search radius");
00801 
00802     wreject = dfs_get_parameter_double(parlist, "fors.fors_calib.wreject", NULL);
00803 
00804     if (wreject <= 0.0)
00805         fors_calib_exit("Invalid rejection threshold");
00806 
00807     wmode = dfs_get_parameter_int(parlist, "fors.fors_calib.wmode", NULL);
00808 
00809     if (wmode < 0 || wmode > 2)
00810         fors_calib_exit("Invalid wavelength solution interpolation mode");
00811 
00812     wcolumn = dfs_get_parameter_string(parlist, "fors.fors_calib.wcolumn", NULL);
00813 
00814     cdegree = dfs_get_parameter_int(parlist, 
00815                     "fors.fors_calib.cdegree", grism_table);
00816 
00817     if (cdegree < 1)
00818         fors_calib_exit("Invalid polynomial degree");
00819 
00820     if (cdegree > 5)
00821         fors_calib_exit("Max allowed polynomial degree is 5");
00822 
00823     cmode = dfs_get_parameter_int(parlist, "fors.fors_calib.cmode", NULL);
00824 
00825     if (cmode < 0 || cmode > 2)
00826         fors_calib_exit("Invalid curvature solution interpolation mode");
00827 
00828     startwavelength = dfs_get_parameter_double(parlist, 
00829                     "fors.fors_calib.startwavelength", grism_table);
00830     if (startwavelength > 1.0)
00831         if (startwavelength < 3000.0 || startwavelength > 13000.0)
00832             fors_calib_exit("Invalid wavelength");
00833 
00834     endwavelength = dfs_get_parameter_double(parlist, 
00835                     "fors.fors_calib.endwavelength", grism_table);
00836     if (endwavelength > 1.0) {
00837         if (endwavelength < 3000.0 || endwavelength > 13000.0)
00838             fors_calib_exit("Invalid wavelength");
00839         if (startwavelength < 1.0)
00840             fors_calib_exit("Invalid wavelength interval");
00841     }
00842 
00843     if (startwavelength > 1.0)
00844         if (endwavelength - startwavelength <= 0.0)
00845             fors_calib_exit("Invalid wavelength interval");
00846 
00847     slit_ident = dfs_get_parameter_bool(parlist, 
00848                     "fors.fors_calib.slit_ident", NULL);
00849 
00850     sdegree = dfs_get_parameter_int(parlist, "fors.fors_calib.sdegree", NULL);
00851     ddegree = dfs_get_parameter_int(parlist, "fors.fors_calib.ddegree", NULL);
00852     sradius = dfs_get_parameter_int(parlist, "fors.fors_calib.sradius", NULL);
00853     dradius = dfs_get_parameter_int(parlist, "fors.fors_calib.dradius", NULL);
00854 
00855     if (sradius < 1 || dradius < 1)
00856         fors_calib_exit("Invalid smoothing box radius");
00857 
00858     qc = dfs_get_parameter_bool(parlist, "fors.fors_calib.qc", NULL);
00859 
00860     check = dfs_get_parameter_bool(parlist, "fors.fors_calib.check", NULL);
00861 
00862     cpl_table_delete(grism_table); grism_table = NULL;
00863 
00864     if (cpl_error_get_code())
00865         fors_calib_exit("Failure getting the configuration parameters");
00866 
00867 
00868     /* 
00869      * Check input set-of-frames
00870      */
00871 
00872     cpl_msg_indent_less();
00873     cpl_msg_info(recipe, "Check input set-of-frames:");
00874     cpl_msg_indent_more();
00875 
00876     if (!dfs_equal_keyword(frameset, "ESO INS GRIS1 ID")) 
00877         fors_calib_exit("Input frames are not from the same grism");
00878 
00879     if (!dfs_equal_keyword(frameset, "ESO INS FILT1 ID")) 
00880         fors_calib_exit("Input frames are not from the same filter");
00881 
00882     if (!dfs_equal_keyword(frameset, "ESO DET CHIP1 ID")) 
00883         fors_calib_exit("Input frames are not from the same chip");
00884 
00885     mxu = cpl_frameset_count_tags(frameset, "LAMP_MXU");
00886     mos = cpl_frameset_count_tags(frameset, "LAMP_MOS");
00887     lss = cpl_frameset_count_tags(frameset, "LAMP_LSS");
00888 
00889     if (mxu + mos + lss == 0)
00890         fors_calib_exit("Missing input arc lamp frame");
00891 
00892     if (mxu + mos + lss > 1)
00893         fors_calib_exit("Just one input arc lamp frame is allowed"); 
00894 
00895     if (mxu) {
00896         cpl_msg_info(recipe, "MXU data found");
00897         arc_tag                  = "LAMP_MXU";
00898         flat_tag                 = "SCREEN_FLAT_MXU";
00899         master_screen_flat_tag   = "MASTER_SCREEN_FLAT_MXU";
00900         master_norm_flat_tag     = "MASTER_NORM_FLAT_MXU";
00901         reduced_lamp_tag         = "REDUCED_LAMP_MXU";
00902         disp_residuals_tag       = "DISP_RESIDUALS_MXU";
00903         disp_coeff_tag           = "DISP_COEFF_MXU";
00904         wavelength_map_tag       = "WAVELENGTH_MAP_MXU";
00905         spectra_detection_tag    = "SPECTRA_DETECTION_MXU";
00906         spectral_resolution_tag  = "SPECTRAL_RESOLUTION_MXU";
00907         slit_map_tag             = "SLIT_MAP_MXU";
00908         curv_traces_tag          = "CURV_TRACES_MXU";
00909         curv_coeff_tag           = "CURV_COEFF_MXU";
00910         spatial_map_tag          = "SPATIAL_MAP_MXU";
00911         slit_location_tag        = "SLIT_LOCATION_MXU";
00912         disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_MXU";
00913         delta_image_tag          = "DELTA_IMAGE_MXU";
00914         mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_MXU";
00915         mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_MXU";
00916     }
00917 
00918     if (lss) {
00919         cpl_msg_info(recipe, "LSS data found");
00920         arc_tag                  = "LAMP_LSS";
00921         flat_tag                 = "SCREEN_FLAT_LSS";
00922         master_screen_flat_tag   = "MASTER_SCREEN_FLAT_LSS";
00923         master_norm_flat_tag     = "MASTER_NORM_FLAT_LSS";
00924         reduced_lamp_tag         = "REDUCED_LAMP_LSS";
00925         spectral_resolution_tag  = "SPECTRAL_RESOLUTION_LSS";
00926         disp_residuals_tag       = "DISP_RESIDUALS_LSS";
00927         disp_coeff_tag           = "DISP_COEFF_LSS";
00928         slit_location_tag        = "SLIT_LOCATION_LSS";
00929         wavelength_map_tag       = "WAVELENGTH_MAP_LSS";
00930         slit_map_tag             = "SLIT_MAP_LSS";
00931         disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_LSS";
00932         delta_image_tag          = "DELTA_IMAGE_LSS";
00933         mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_LSS";
00934         mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_LSS";
00935     }
00936 
00937     if (mos) {
00938         cpl_msg_info(recipe, "MOS data found");
00939         arc_tag                  = "LAMP_MOS";
00940         flat_tag                 = "SCREEN_FLAT_MOS";
00941         master_screen_flat_tag   = "MASTER_SCREEN_FLAT_MOS";
00942         master_norm_flat_tag     = "MASTER_NORM_FLAT_MOS";
00943         reduced_lamp_tag         = "REDUCED_LAMP_MOS";
00944         disp_residuals_tag       = "DISP_RESIDUALS_MOS";
00945         disp_coeff_tag           = "DISP_COEFF_MOS";
00946         wavelength_map_tag       = "WAVELENGTH_MAP_MOS";
00947         spectra_detection_tag    = "SPECTRA_DETECTION_MOS";
00948         spectral_resolution_tag  = "SPECTRAL_RESOLUTION_MOS";
00949         slit_map_tag             = "SLIT_MAP_MOS";
00950         curv_traces_tag          = "CURV_TRACES_MOS";
00951         curv_coeff_tag           = "CURV_COEFF_MOS";
00952         spatial_map_tag          = "SPATIAL_MAP_MOS";
00953         slit_location_tag        = "SLIT_LOCATION_MOS";
00954         disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_MOS";
00955         delta_image_tag          = "DELTA_IMAGE_MOS";
00956         mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_MOS";
00957         mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_MOS";
00958     }
00959 
00960     nbias = 0;
00961     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") == 0) {
00962         if (cpl_frameset_count_tags(frameset, "BIAS") == 0)
00963             fors_calib_exit("Missing required input: MASTER_BIAS or BIAS");
00964         nbias = cpl_frameset_count_tags(frameset, "BIAS");
00965     }
00966 
00967     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") > 1)
00968         fors_calib_exit("Too many in input: MASTER_BIAS");
00969 
00970     if (cpl_frameset_count_tags(frameset, "MASTER_LINECAT") == 0)
00971         fors_calib_exit("Missing required input: MASTER_LINECAT");
00972 
00973     if (cpl_frameset_count_tags(frameset, "MASTER_LINECAT") > 1)
00974         fors_calib_exit("Too many in input: MASTER_LINECAT");
00975 
00976     nflats = cpl_frameset_count_tags(frameset, flat_tag);
00977 
00978     if (nflats < 1) {
00979         cpl_msg_error(recipe, "Missing required input: %s", flat_tag);
00980         fors_calib_exit(NULL);
00981     }
00982 
00983     cpl_msg_indent_less();
00984 
00985     if (nflats > 1)
00986         cpl_msg_info(recipe, "Load %d flat field frames and sum them...",
00987                      nflats);
00988     else
00989         cpl_msg_info(recipe, "Load flat field exposure...");
00990 
00991     cpl_msg_indent_more();
00992 
00993     header = dfs_load_header(frameset, flat_tag, 0);
00994 
00995     if (header == NULL)
00996         fors_calib_exit("Cannot load flat field frame header");
00997 
00998     alltime = cpl_propertylist_get_double(header, "EXPTIME");
00999 
01000     if (cpl_error_get_code() != CPL_ERROR_NONE)
01001         fors_calib_exit("Missing keyword EXPTIME in flat field frame header");
01002 
01003     cpl_propertylist_delete(header);
01004 
01005     for (i = 1; i < nflats; i++) {
01006 
01007         header = dfs_load_header(frameset, NULL, 0);
01008 
01009         if (header == NULL)
01010             fors_calib_exit("Cannot load flat field frame header");
01011 
01012         alltime += cpl_propertylist_get_double(header, "EXPTIME");
01013 
01014         if (cpl_error_get_code() != CPL_ERROR_NONE)
01015             fors_calib_exit("Missing keyword EXPTIME in flat field "
01016                             "frame header");
01017 
01018         cpl_propertylist_delete(header);
01019 
01020     }
01021 
01022     master_flat = dfs_load_image(frameset, flat_tag, CPL_TYPE_FLOAT, 0, 0);
01023 
01024     if (master_flat == NULL)
01025         fors_calib_exit("Cannot load flat field");
01026 
01027     for (i = 1; i < nflats; i++) {
01028         flat = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
01029         if (flat) {
01030             cpl_image_add(master_flat, flat);
01031             cpl_image_delete(flat); flat = NULL;
01032         }
01033         else
01034             fors_calib_exit("Cannot load flat field");
01035     }
01036 
01037 /***
01038     if (nflats > 1)
01039         cpl_image_divide_scalar(master_flat, nflats);
01040 ***/
01041 
01042 
01043     /*
01044      * Get the reference wavelength and the rebin factor along the
01045      * dispersion direction from the arc lamp exposure
01046      */
01047 
01048     header = dfs_load_header(frameset, arc_tag, 0);
01049 
01050     if (header == NULL)
01051         fors_calib_exit("Cannot load arc lamp header");
01052 
01053     instrume = (char *)cpl_propertylist_get_string(header, "INSTRUME");
01054     if (instrume == NULL)
01055         fors_calib_exit("Missing keyword INSTRUME in arc lamp header");
01056     instrume = cpl_strdup(instrume);
01057 
01058     if (instrume[4] == '1')
01059         snprintf(version, 80, "%s/%s", "fors1", VERSION);
01060     if (instrume[4] == '2')
01061         snprintf(version, 80, "%s/%s", "fors2", VERSION);
01062 
01063     reference = cpl_propertylist_get_double(header, "ESO INS GRIS1 WLEN");
01064 
01065     if (cpl_error_get_code() != CPL_ERROR_NONE)
01066         fors_calib_exit("Missing keyword ESO INS GRIS1 WLEN in arc lamp "
01067                         "frame header");
01068 
01069     if (reference < 3000.0)   /* Perhaps in nanometers... */
01070         reference *= 10;
01071 
01072     if (reference < 3000.0 || reference > 13000.0) {
01073         cpl_msg_error(recipe, "Invalid central wavelength %.2f read from "
01074                       "keyword ESO INS GRIS1 WLEN in arc lamp frame header",
01075                       reference);
01076         fors_calib_exit(NULL);
01077     }
01078 
01079     cpl_msg_info(recipe, "The central wavelength is: %.2f", reference);
01080 
01081     rebin = cpl_propertylist_get_int(header, "ESO DET WIN1 BINX");
01082 
01083     if (cpl_error_get_code() != CPL_ERROR_NONE)
01084         fors_calib_exit("Missing keyword ESO DET WIN1 BINX in arc lamp "
01085                         "frame header");
01086 
01087     if (rebin != 1) {
01088         dispersion *= rebin;
01089         cpl_msg_warning(recipe, "The rebin factor is %d, and therefore the "
01090                         "working dispersion used is %f A/pixel", rebin, 
01091                         dispersion);
01092     }
01093 
01094     gain = cpl_propertylist_get_double(header, "ESO DET OUT1 CONAD");
01095 
01096     if (cpl_error_get_code() != CPL_ERROR_NONE)
01097         fors_calib_exit("Missing keyword ESO DET OUT1 CONAD in arc lamp "
01098                         "frame header");
01099 
01100     cpl_msg_info(recipe, "The gain factor is: %.2f e-/ADU", gain);
01101 
01102     if (mos || mxu) {
01103 /* goto skip; */
01104         cpl_msg_info(recipe, "Produce mask slit position table...");
01105         if (mos)
01106             maskslits = mos_load_slits_fors_mos(header);
01107         else
01108             maskslits = mos_load_slits_fors_mxu(header);
01109 
01110         /*
01111          * Check if all slits have the same X offset: in such case, 
01112          * treat the observation as a long-slit one!
01113          */
01114 
01115         mxpos = cpl_table_get_column_median(maskslits, "xtop");
01116         xpos = cpl_table_get_data_double(maskslits, "xtop");
01117         nslits = cpl_table_get_nrow(maskslits);
01118 
01119         treat_as_lss = 1;
01120         for (i = 0; i < nslits; i++) {
01121             if (fabs(mxpos-xpos[i]) > 0.01) {
01122                 treat_as_lss = 0;
01123                 break;
01124             }
01125         }
01126 /* skip: */
01127 
01128         if (treat_as_lss) {
01129             cpl_msg_warning(recipe, "All MOS slits have the same offset: %.2f\n"
01130                             "The LSS data reduction strategy is applied!", 
01131                             mxpos);
01132             cpl_table_delete(maskslits); maskslits = NULL;
01133             if (mos) {
01134                 master_screen_flat_tag   = "MASTER_SCREEN_FLAT_LONG_MOS";
01135                 master_norm_flat_tag     = "MASTER_NORM_FLAT_LONG_MOS";
01136                 disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_LONG_MOS";
01137                 delta_image_tag          = "DELTA_IMAGE_LONG_MOS";
01138                 spectral_resolution_tag  = "SPECTRAL_RESOLUTION_LONG_MOS";
01139                 reduced_lamp_tag         = "REDUCED_LAMP_LONG_MOS";
01140                 disp_coeff_tag           = "DISP_COEFF_LONG_MOS";
01141                 wavelength_map_tag       = "WAVELENGTH_MAP_LONG_MOS";
01142                 slit_location_tag        = "SLIT_LOCATION_LONG_MOS";
01143                 mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_LONG_MOS";
01144                 mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_LONG_MOS";
01145             }
01146             else {
01147                 master_screen_flat_tag   = "MASTER_SCREEN_FLAT_LONG_MXU";
01148                 master_norm_flat_tag     = "MASTER_NORM_FLAT_LONG_MXU";
01149                 disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_LONG_MXU";
01150                 delta_image_tag          = "DELTA_IMAGE_LONG_MXU";
01151                 spectral_resolution_tag  = "SPECTRAL_RESOLUTION_LONG_MXU";
01152                 reduced_lamp_tag         = "REDUCED_LAMP_LONG_MXU";
01153                 disp_coeff_tag           = "DISP_COEFF_LONG_MXU";
01154                 wavelength_map_tag       = "WAVELENGTH_MAP_LONG_MXU";
01155                 slit_location_tag        = "SLIT_LOCATION_LONG_MXU";
01156                 mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_LONG_MXU";
01157                 mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_LONG_MXU";
01158             }
01159         }
01160     }
01161 
01162     if (slit_ident == 0) {
01163         cpl_table_delete(maskslits); maskslits = NULL;
01164     }
01165 
01166 
01167     /* Leave the header on for the next step... */
01168 
01169 
01170     /*
01171      * Remove the master bias
01172      */
01173 
01174     if (nbias) {
01175 
01176         /*
01177          * Set of raw BIASes in input, need to create master bias!
01178          */
01179 
01180         cpl_msg_info(recipe, "Generate the master from input raw biases...");
01181 
01182         if (nbias > 1) {
01183 
01184             biases = cpl_imagelist_new();
01185 
01186             bias = dfs_load_image(frameset, "BIAS", CPL_TYPE_FLOAT, 0, 0);
01187     
01188             if (bias == NULL)
01189                 fors_calib_exit("Cannot load bias frame");
01190 
01191             cpl_imagelist_set(biases, bias, 0);
01192     
01193             for (i = 1; i < nbias; i++) {
01194                 bias = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
01195                 if (bias)
01196                     cpl_imagelist_set(biases, bias, i);
01197                 else
01198                     fors_calib_exit("Cannot load bias frame");
01199             }
01200     
01201             master_bias = cpl_imagelist_collapse_median_create(biases);
01202 
01203             cpl_imagelist_delete(biases);
01204         }
01205         else {
01206             master_bias = dfs_load_image(frameset, "BIAS", 
01207                                          CPL_TYPE_FLOAT, 0, 1);
01208             if (master_bias == NULL)
01209                 fors_calib_exit("Cannot load bias");
01210         }
01211 
01212     }
01213     else {
01214         master_bias = dfs_load_image(frameset, "MASTER_BIAS", 
01215                                      CPL_TYPE_FLOAT, 0, 1);
01216         if (master_bias == NULL)
01217             fors_calib_exit("Cannot load master bias");
01218     }
01219 
01220     cpl_msg_info(recipe, "Remove the master bias...");
01221 
01222     overscans = mos_load_overscans_vimos(header, 1);
01223     cpl_propertylist_delete(header); header = NULL;
01224 
01225     if (nbias) {
01226         int xlow = cpl_table_get_int(overscans, "xlow", 0, NULL);
01227         int ylow = cpl_table_get_int(overscans, "ylow", 0, NULL);
01228         int xhig = cpl_table_get_int(overscans, "xhig", 0, NULL);
01229         int yhig = cpl_table_get_int(overscans, "yhig", 0, NULL);
01230         dummy = cpl_image_extract(master_bias, xlow+1, ylow+1, xhig, yhig);
01231         cpl_image_delete(master_bias); master_bias = dummy;
01232 
01233         if (dfs_save_image(frameset, master_bias, "MASTER_BIAS",
01234                            NULL, parlist, recipe, version))
01235             fors_calib_exit(NULL);
01236     }
01237 
01238     if (nflats > 1) {
01239         multi_bias = cpl_image_multiply_scalar_create(master_bias, nflats);
01240         dummy = mos_remove_bias(master_flat, multi_bias, overscans);
01241         cpl_image_delete(multi_bias);
01242     }
01243     else {
01244         dummy = mos_remove_bias(master_flat, master_bias, overscans);
01245     }
01246     cpl_image_delete(master_flat);
01247     master_flat = dummy;
01248 
01249     if (master_flat == NULL)
01250         fors_calib_exit("Cannot remove bias from flat field");
01251 
01252     cpl_msg_indent_less();
01253     cpl_msg_info(recipe, "Load arc lamp exposure...");
01254     cpl_msg_indent_more();
01255 
01256     spectra = dfs_load_image(frameset, arc_tag, CPL_TYPE_FLOAT, 0, 0);
01257 
01258     if (spectra == NULL)
01259         fors_calib_exit("Cannot load arc lamp exposure");
01260 
01261     cpl_msg_info(recipe, "Remove the master bias...");
01262 
01263     dummy = mos_remove_bias(spectra, master_bias, overscans);
01264     cpl_table_delete(overscans); overscans = NULL;
01265     cpl_image_delete(master_bias); master_bias = NULL;
01266     cpl_image_delete(spectra); spectra = dummy;
01267 
01268     if (spectra == NULL)
01269         fors_calib_exit("Cannot remove bias from arc lamp exposure");
01270 
01271     cpl_msg_indent_less();
01272     cpl_msg_info(recipe, "Load input line catalog...");
01273     cpl_msg_indent_more();
01274 
01275     wavelengths = dfs_load_table(frameset, "MASTER_LINECAT", 1);
01276 
01277     if (wavelengths == NULL)
01278         fors_calib_exit("Cannot load line catalog");
01279 
01280 
01281     /*
01282      * Cast the wavelengths into a (double precision) CPL vector
01283      */
01284 
01285     nlines = cpl_table_get_nrow(wavelengths);
01286 
01287     if (nlines == 0)
01288         fors_calib_exit("Empty input line catalog");
01289 
01290     if (cpl_table_has_column(wavelengths, wcolumn) != 1) {
01291         cpl_msg_error(recipe, "Missing column %s in input line catalog table",
01292                       wcolumn);
01293         fors_calib_exit(NULL);
01294     }
01295 
01296     line = cpl_malloc(nlines * sizeof(double));
01297     
01298     for (i = 0; i < nlines; i++)
01299         line[i] = cpl_table_get(wavelengths, wcolumn, i, NULL);
01300 
01301     cpl_table_delete(wavelengths); wavelengths = NULL;
01302 
01303     lines = cpl_vector_wrap(nlines, line);
01304 
01305 
01306     if (lss || treat_as_lss) {
01307 
01308         int first_row, last_row;
01309         int ylow, yhig;
01310 
01311         /* FIXME:
01312          * The LSS data calibration is still dirty: it doesn't apply
01313          * any spatial rectification, and only in future an external
01314          * spectral curvature model would be provided in input. Here
01315          * and there temporary solutions are adpted, such as accepting
01316          * the preliminary wavelength calibration.
01317          */
01318 
01319         /*
01320          * Flat field normalisation is done directly on the master flat
01321          * field (without spatial rectification first). The spectral
01322          * curvature model may be provided in input, in future releases.
01323          */
01324 
01325         cpl_msg_indent_less();
01326         cpl_msg_info(recipe, "Perform flat field normalisation...");
01327         cpl_msg_indent_more();
01328 
01329         norm_flat = cpl_image_duplicate(master_flat);
01330 
01331         smo_flat = mos_normalise_longflat(norm_flat, sradius, dradius,
01332                                           sdegree);
01333 
01334         cpl_image_delete(smo_flat); smo_flat = NULL; /* It may be a product */
01335 
01336 /*        if (!qc) {         enable when if (1) disabled *******/
01337 
01338         if (1) {
01339             save_header = dfs_load_header(frameset, flat_tag, 0);
01340             cpl_propertylist_update_int(save_header, 
01341                                         "ESO PRO DATANCOM", nflats);
01342 
01343             if (dfs_save_image(frameset, master_flat, master_screen_flat_tag,
01344                                save_header, parlist, recipe, version))
01345                 fors_calib_exit(NULL);
01346 //%%%%%            cpl_image_delete(master_flat); master_flat = NULL;
01347         }
01348 
01349         if (dfs_save_image(frameset, norm_flat, master_norm_flat_tag,
01350                            save_header, parlist, recipe, version))
01351             fors_calib_exit(NULL);
01352 
01353         cpl_propertylist_delete(save_header); save_header = NULL;
01354 //%%%%%        cpl_image_delete(norm_flat); norm_flat = NULL;
01355 
01356 
01357         /*
01358          * In the case of LSS data, extract the spectra directly
01359          * on the first attempt. The spectral curvature model may
01360          * be provided in input, in future releases.
01361          */
01362 
01363         cpl_msg_indent_less();
01364         cpl_msg_info(recipe, "Perform wavelength calibration...");
01365         cpl_msg_indent_more();
01366 
01367         nx = cpl_image_get_size_x(spectra);
01368         ny = cpl_image_get_size_y(spectra);
01369 
01370         wavemap = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01371         idscoeff_all = cpl_table_new(ny);
01372 
01373         rectified = mos_wavelength_calibration_raw(spectra, lines, dispersion,
01374                                                    peakdetection, wradius,
01375                                                    wdegree, wreject, reference,
01376                                                    &startwavelength,
01377                                                    &endwavelength, NULL,
01378                                                    NULL, idscoeff_all, wavemap,
01379                                                    NULL, NULL, NULL);
01380 
01381         if (rectified == NULL)
01382             fors_calib_exit("Wavelength calibration failure.");
01383 
01384         cpl_image_delete(rectified); rectified = NULL;
01385 
01386         first_row = 0;
01387         while (!cpl_table_is_valid(idscoeff_all, "c0", first_row))
01388             first_row++;
01389 
01390         last_row = ny - 1;
01391         while (!cpl_table_is_valid(idscoeff_all, "c0", last_row))
01392             last_row--;
01393 
01394         ylow = first_row + 1;
01395         yhig = last_row + 1;
01396 
01397         if (ylow >= yhig) {
01398             cpl_error_reset();
01399             fors_calib_exit("No spectra could be detected.");
01400         }
01401 
01402         cpl_msg_info(recipe, 
01403                      "Spectral pattern was detected on %d out of %d CCD rows", 
01404                      yhig - ylow, ny);
01405 
01406         dummy = cpl_image_extract(spectra, 1, ylow, nx, yhig);
01407         cpl_image_delete(spectra); spectra = dummy;
01408 
01409         ccd_ysize = ny;
01410         ny = cpl_image_get_size_y(spectra);
01411 
01412         if (check)
01413             residual = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01414 
01415         fiterror = cpl_calloc(ny, sizeof(double));
01416         fitlines = cpl_calloc(ny, sizeof(int));
01417         idscoeff = cpl_table_new(ny);
01418         restable = cpl_table_new(nlines);
01419 
01420         rectified = mos_wavelength_calibration_raw(spectra, lines, dispersion,
01421                                                    peakdetection, wradius,
01422                                                    wdegree, wreject, reference,
01423                                                    &startwavelength, 
01424                                                    &endwavelength, fitlines, 
01425                                                    fiterror, idscoeff, NULL,
01426                                                    residual, restable, NULL);
01427 
01428         if (rectified == NULL)
01429             fors_calib_exit("Wavelength calibration failure.");
01430 
01431         /*
01432          * A dummy slit locations table
01433          */
01434 
01435         slits = cpl_table_new(1);
01436         cpl_table_new_column(slits, "slit_id", CPL_TYPE_INT);
01437         cpl_table_new_column(slits, "xtop", CPL_TYPE_DOUBLE);
01438         cpl_table_new_column(slits, "ytop", CPL_TYPE_DOUBLE);
01439         cpl_table_new_column(slits, "xbottom", CPL_TYPE_DOUBLE);
01440         cpl_table_new_column(slits, "ybottom", CPL_TYPE_DOUBLE);
01441         cpl_table_new_column(slits, "position", CPL_TYPE_INT);
01442         cpl_table_new_column(slits, "length", CPL_TYPE_INT);
01443         cpl_table_set_column_unit(slits, "xtop", "pixel");
01444         cpl_table_set_column_unit(slits, "ytop", "pixel");
01445         cpl_table_set_column_unit(slits, "xbottom", "pixel");
01446         cpl_table_set_column_unit(slits, "ybottom", "pixel");
01447         cpl_table_set_column_unit(slits, "position", "pixel");
01448         cpl_table_set_column_unit(slits, "length", "pixel");
01449         cpl_table_set_int(slits, "slit_id", 0, 0);
01450         cpl_table_set_double(slits, "xtop", 0, 0);
01451         cpl_table_set_double(slits, "ytop", 0, last_row);
01452         cpl_table_set_double(slits, "xbottom", 0, 0);
01453         cpl_table_set_double(slits, "ybottom", 0, first_row);
01454         cpl_table_set_int(slits, "position", 0, 0);
01455         cpl_table_set_int(slits, "length", 0, ny);
01456 
01457         if (dfs_save_table(frameset, slits, slit_location_tag, NULL,
01458                            parlist, recipe, version))
01459             fors_calib_exit(NULL);
01460 
01461         cpl_table_delete(slits); slits = NULL;
01462 
01463         if (dfs_save_table(frameset, restable, disp_residuals_table_tag, NULL,
01464                            parlist, recipe, version))
01465             fors_calib_exit(NULL);
01466 
01467         cpl_table_delete(restable); restable = NULL;
01468 
01469         if (wmode) {
01470             cpl_image_delete(rectified); rectified = NULL;
01471             cpl_image_delete(wavemap); wavemap = NULL;
01472             mos_interpolate_wavecalib(idscoeff, wavemap, wmode);
01473             mos_interpolate_wavecalib(idscoeff_all, wavemap, wmode);
01474             wavemap = mos_map_idscoeff(idscoeff_all, nx, reference,
01475                                        startwavelength, endwavelength);
01476             rectified = mos_wavelength_calibration(spectra, reference,
01477                                                    startwavelength, 
01478                                                    endwavelength, dispersion, 
01479                                                    idscoeff, 0);
01480         }
01481 
01482         cpl_table_delete(idscoeff_all); idscoeff_all = NULL;
01483 
01484         cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
01485         cpl_table_set_column_unit(idscoeff, "error", "pixel");
01486         cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;
01487 
01488         for (i = 0; i < ny; i++)
01489             if (!cpl_table_is_valid(idscoeff, "c0", i))
01490                 cpl_table_set_invalid(idscoeff, "error", i);
01491 
01492         delta = mos_map_pixel(idscoeff, reference, startwavelength,
01493                               endwavelength, dispersion, 2);
01494 
01495 //%%%%%
01496         dummy = cpl_image_extract(master_flat, 1, ylow, nx, yhig);
01497         cpl_image_delete(master_flat); master_flat = dummy;
01498 
01499         mapped_flat = mos_wavelength_calibration(master_flat, reference,
01500                                       startwavelength, endwavelength,
01501                                       dispersion, idscoeff, 0);
01502 
01503         cpl_image_delete(master_flat); master_flat = NULL;
01504 
01505         dummy = cpl_image_extract(norm_flat, 1, ylow, nx, yhig);
01506         cpl_image_delete(norm_flat); norm_flat = dummy;
01507 
01508         mapped_nflat = mos_wavelength_calibration(norm_flat, reference,
01509                                       startwavelength, endwavelength,
01510                                       dispersion, idscoeff, 0);
01511 
01512         cpl_image_delete(norm_flat); norm_flat = NULL;
01513 
01514         header = cpl_propertylist_new();
01515         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01516         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01517         cpl_propertylist_update_double(header, "CRVAL1",
01518                                        startwavelength + dispersion/2);
01519         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01520         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01521         cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01522         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01523         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01524         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01525         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01526         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01527         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01528 
01529         if (dfs_save_image(frameset, delta, delta_image_tag,
01530                            header, parlist, recipe, version))
01531             fors_calib_exit(NULL);
01532 
01533         cpl_image_delete(delta); delta = NULL;
01534 
01535         if (dfs_save_image(frameset, mapped_flat, mapped_screen_flat_tag,
01536                            header, parlist, recipe, version))
01537             fors_calib_exit(NULL);
01538 
01539         cpl_image_delete(mapped_flat); mapped_flat = NULL;
01540 
01541         if (dfs_save_image(frameset, mapped_nflat, mapped_norm_flat_tag,
01542                            header, parlist, recipe, version))
01543             fors_calib_exit(NULL);
01544 
01545         cpl_image_delete(mapped_nflat); mapped_nflat = NULL;
01546 
01547         cpl_propertylist_delete(header); header = NULL;
01548 
01549         cpl_msg_info(recipe, "Valid solutions found: %d out of %d rows", 
01550                      ny - cpl_table_count_invalid(idscoeff, "c0"), ny);
01551 
01552         cpl_image_delete(spectra); spectra = NULL;
01553 
01554         mean_rms = mos_distortions_rms(rectified, lines, startwavelength,
01555                                        dispersion, 6, 0);
01556 
01557         cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);
01558 
01559         mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01560         mean_rms_err = cpl_table_get_column_stdev(idscoeff, "error");
01561 
01562         cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)",
01563                      mean_rms, mean_rms * dispersion);
01564 
01565         restab = mos_resolution_table(rectified, startwavelength, dispersion,
01566                                       60000, lines);
01567 
01568         if (restab) {
01569             cpl_msg_info(recipe, "Mean spectral resolution: %.2f",
01570                   cpl_table_get_column_mean(restab, "resolution"));
01571             cpl_msg_info(recipe, 
01572                   "Mean reference lines FWHM: %.2f +/- %.2f pixel",
01573                   cpl_table_get_column_mean(restab, "fwhm") / dispersion,
01574                   cpl_table_get_column_mean(restab, "fwhm_rms") / dispersion);
01575 
01576             if (qc) {
01577 
01578                 header = dfs_load_header(frameset, arc_tag, 0);
01579 
01580                 if (header == NULL)
01581                     fors_calib_exit("Cannot reload arc lamp header");
01582 
01583                 qclist = cpl_propertylist_new();
01584 
01585                 fors_qc_start_group(qclist, "2.0", instrume);
01586 
01587 
01588                 /*
01589                  * QC1 group header
01590                  */
01591 
01592                 if (fors_qc_write_string("PRO.CATG", spectral_resolution_tag,
01593                                         "Product category", instrume))
01594                     fors_calib_exit("Cannot write product category to "
01595                                     "QC log file");
01596 
01597                 if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL,
01598                                           "DPR type", instrume))
01599                     fors_calib_exit("Missing keyword DPR TYPE in arc "
01600                                     "lamp header");
01601 
01602                 if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL,
01603                                           "Template", instrume))
01604                     fors_calib_exit("Missing keyword TPL ID in arc "
01605                                     "lamp header");
01606 
01607                 if (fors_qc_keyword_to_paf(header, "ESO INS GRIS1 NAME", NULL,
01608                                           "Grism name", instrume))
01609                     fors_calib_exit("Missing keyword INS GRIS1 NAME in arc "
01610                                     "lamp header");
01611 
01612                 if (fors_qc_keyword_to_paf(header, "ESO INS GRIS1 ID", NULL,
01613                                           "Grism identifier", instrume))
01614                     fors_calib_exit("Missing keyword INS GRIS1 ID in arc "
01615                                     "lamp header");
01616 
01617                 if (cpl_propertylist_has(header, "ESO INS FILT1 NAME"))
01618                     fors_qc_keyword_to_paf(header, "ESO INS FILT1 NAME", NULL,
01619                                           "Filter name", instrume);
01620 
01621                 if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
01622                                           "Collimator name", instrume))
01623                     fors_calib_exit("Missing keyword INS COLL NAME in arc "
01624                                     "lamp header");
01625 
01626                 if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
01627                                           "Chip identifier", instrume))
01628                     fors_calib_exit("Missing keyword DET CHIP1 ID in arc "
01629                                     "lamp header");
01630 
01631                 if (mos) {
01632                     if (fors_qc_keyword_to_paf(header, "ESO INS MOS10 WID",
01633                                               "arcsec", "Slit width", instrume))
01634                         fors_calib_exit("Missing keyword ESO INS MOS1 WID in "
01635                                         "arc lamp header");
01636                 }
01637                 else {
01638                     if (fors_qc_keyword_to_paf(header, "ESO INS SLIT WID",
01639                                               "arcsec", "Slit width", instrume))
01640                         fors_calib_exit("Missing keyword ESO INS SLIT WID in "
01641                                         "arc lamp header");
01642                 }
01643 
01644                 if (fors_qc_keyword_to_paf(header, "ESO DET WIN1 BINX", NULL,    
01645                                           "Binning factor along X", instrume))
01646                     fors_calib_exit("Missing keyword ESO DET WIN1 BINX "
01647                                       "in frame header");
01648 
01649                 if (fors_qc_keyword_to_paf(header, "ESO DET WIN1 BINY", NULL,
01650                                           "Binning factor along Y", instrume))
01651                     fors_calib_exit("Missing keyword ESO DET WIN1 BINY "
01652                                       "in frame header");
01653 
01654                 if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
01655                                           "Archive name of input data",
01656                                           instrume))
01657                     fors_calib_exit("Missing keyword ARCFILE in arc "
01658                                     "lamp header");
01659 
01660                 cpl_propertylist_delete(header); header = NULL;
01661 
01662                 pipefile = dfs_generate_filename_tfits(spectral_resolution_tag);
01663                 if (fors_qc_write_string("PIPEFILE", pipefile,
01664                                         "Pipeline product name", instrume))
01665                     fors_calib_exit("Cannot write PIPEFILE to QC log file");
01666                 cpl_free(pipefile); pipefile = NULL;
01667 
01668 
01669                 /*
01670                  * QC1 parameters
01671                  */
01672    
01673                 if (mos)
01674                     keyname = "QC.MOS.RESOLUTION";
01675                 else
01676                     keyname = "QC.LSS.RESOLUTION";
01677 
01678                 if (fors_qc_write_qc_double(qclist, 
01679                                            cpl_table_get_column_mean(restab,
01680                                                                  "resolution"),
01681                                            keyname, NULL, 
01682                                            "Mean spectral resolution",
01683                                            instrume)) {
01684                     fors_calib_exit("Cannot write mean spectral resolution to "
01685                                     "QC log file");
01686                 }
01687 
01688                 if (mos)
01689                     keyname = "QC.MOS.RESOLUTION.RMS"; 
01690                 else
01691                     keyname = "QC.LSS.RESOLUTION.RMS";
01692 
01693                 if (fors_qc_write_qc_double(qclist,
01694                                            cpl_table_get_column_stdev(restab,
01695                                                                 "resolution"),
01696                                            keyname, NULL, 
01697                                            "Scatter of spectral resolution",
01698                                            instrume)) {
01699                     fors_calib_exit("Cannot write spectral resolution scatter "
01700                                     "to QC log file");
01701                 }
01702 
01703                 if (mos)
01704                     keyname = "QC.MOS.RESOLUTION.NWAVE";
01705                 else
01706                     keyname = "QC.LSS.RESOLUTION.NWAVE";
01707 
01708                 if (fors_qc_write_qc_int(qclist, cpl_table_get_nrow(restab) -
01709                                         cpl_table_count_invalid(restab,
01710                                                                 "resolution"),
01711                                         keyname, NULL,
01712                                         "Number of examined wavelengths "
01713                                         "for resolution computation",
01714                                         instrume)) {
01715                     fors_calib_exit("Cannot write number of lines used in "
01716                                     "spectral resolution computation "
01717                                     "to QC log file");
01718                 }
01719 
01720                 if (mos)
01721                     keyname = "QC.MOS.RESOLUTION.MEANRMS";
01722                 else
01723                     keyname = "QC.LSS.RESOLUTION.MEANRMS";
01724 
01725                 if (fors_qc_write_qc_double(qclist, 
01726                                            cpl_table_get_column_mean(restab,
01727                                                          "resolution_rms"),
01728                                         keyname, NULL,
01729                                         "Mean error on spectral "
01730                                         "resolution computation",
01731                                         instrume)) {
01732                     fors_calib_exit("Cannot write mean error in "
01733                                     "spectral resolution computation "
01734                                     "to QC log file");
01735                 }
01736 
01737                 if (mos)
01738                     keyname = "QC.MOS.RESOLUTION.NLINES";
01739                 else
01740                     keyname = "QC.LSS.RESOLUTION.NLINES";
01741 
01742                 if (fors_qc_write_qc_int(qclist, 
01743                          cpl_table_get_column_mean(restab, "nlines") *
01744                                            cpl_table_get_nrow(restab),
01745                                         keyname, NULL,
01746                                         "Number of lines for spectral "
01747                                         "resolution computation",
01748                                         instrume)) {
01749                     fors_calib_exit("Cannot write number of examined "
01750                          "wavelengths in spectral resolution computation "
01751                          "to QC log file");
01752                 }
01753 
01754                 fors_qc_end_group();
01755 
01756             }  /* End of QC1 computation */
01757 
01758             if (dfs_save_table(frameset, restab, spectral_resolution_tag, 
01759                                qclist, parlist, recipe, version))
01760                 fors_calib_exit(NULL);
01761 
01762             cpl_table_delete(restab); restab = NULL;
01763             cpl_propertylist_delete(qclist); qclist = NULL;
01764 
01765         }
01766         else
01767             fors_calib_exit("Cannot compute the spectral resolution table");
01768 
01769         cpl_vector_delete(lines); lines = NULL;
01770 
01771 
01772         /*
01773          * Save rectified arc lamp spectrum to disk
01774          */
01775 
01776         header = cpl_propertylist_new();
01777         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01778         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01779         cpl_propertylist_update_double(header, "CRVAL1", 
01780                                        startwavelength + dispersion/2);
01781         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01782         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01783         cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01784         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01785         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01786         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01787         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01788         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01789         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01790         cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);
01791 
01792         if (dfs_save_image(frameset, rectified, reduced_lamp_tag, header, 
01793                            parlist, recipe, version))
01794             fors_calib_exit(NULL);
01795 
01796         cpl_image_delete(rectified); rectified = NULL;
01797         cpl_propertylist_delete(header); header = NULL;
01798 
01799         if (dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL, 
01800                            parlist, recipe, version))
01801             fors_calib_exit(NULL);
01802 
01803         cpl_table_delete(idscoeff); idscoeff = NULL;
01804 
01805         header = dfs_load_header(frameset, arc_tag, 0);
01806 
01807         if (header == NULL)
01808             fors_calib_exit("Cannot reload arc lamp header");
01809 
01810         if (qc) {
01811 
01812             compute_central_wave = 0;
01813             if (lss) {
01814 /***
01815                 if (fabs(1.0 - cpl_propertylist_get_double(header,
01816                                              "ESO INS SLIT WID")) < 0.05)
01817  ***/
01818                 compute_central_wave = 1;
01819             }
01820             else {
01821                 if (fabs(mxpos) < 0.05)
01822                     compute_central_wave = 1;
01823             }
01824 
01825             fors_qc_start_group(header, "2.0", instrume);
01826 
01827             /*
01828              * QC1 group header
01829              */
01830 
01831             if (fors_qc_write_string("PRO.CATG", wavelength_map_tag,
01832                                     "Product category", instrume))
01833                 fors_calib_exit("Cannot write product category to "
01834                                 "QC log file");
01835 
01836             if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL, 
01837                                       "DPR type", instrume))
01838                 fors_calib_exit("Missing keyword DPR TYPE in arc "
01839                                 "lamp header");
01840 
01841             if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL, 
01842                                       "Template", instrume))
01843                 fors_calib_exit("Missing keyword TPL ID in arc "
01844                                 "lamp header");
01845 
01846             if (fors_qc_keyword_to_paf(header, "ESO INS GRIS1 NAME", NULL, 
01847                                       "Grism name", instrume))
01848                 fors_calib_exit("Missing keyword INS GRIS1 NAME in arc "
01849                                 "lamp header");
01850 
01851             if (fors_qc_keyword_to_paf(header, "ESO INS GRIS1 ID", NULL,
01852                                       "Grism identifier", instrume))
01853                 fors_calib_exit("Missing keyword INS GRIS1 ID in arc "
01854                                 "lamp header");
01855 
01856             if (cpl_propertylist_has(header, "ESO INS FILT1 NAME"))
01857                 fors_qc_keyword_to_paf(header, "ESO INS FILT1 NAME", NULL,
01858                                       "Filter name", instrume);
01859 
01860             if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
01861                                       "Collimator name", instrume))
01862                 fors_calib_exit("Missing keyword INS COLL NAME in arc "
01863                                 "lamp header");
01864 
01865             if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
01866                                       "Chip identifier", instrume))
01867                 fors_calib_exit("Missing keyword DET CHIP1 ID in arc "
01868                                 "lamp header");
01869 
01870             if (mos) {
01871                 if (fors_qc_keyword_to_paf(header, "ESO INS MOS10 WID",
01872                                           "arcsec", "Slit width", instrume))
01873                     fors_calib_exit("Missing keyword ESO INS MOS1 WID in "
01874                                     "arc lamp header");
01875             }
01876             else {
01877                 if (fors_qc_keyword_to_paf(header, "ESO INS SLIT WID",
01878                                           "arcsec", "Slit width", instrume))
01879                     fors_calib_exit("Missing keyword ESO INS SLIT WID in "
01880                                     "arc lamp header");
01881             }
01882 
01883             if (fors_qc_keyword_to_paf(header, "ESO DET WIN1 BINX", NULL,    
01884                                       "Binning factor along X", instrume))
01885                 fors_calib_exit("Missing keyword ESO DET WIN1 BINX "
01886                                   "in frame header");
01887 
01888             if (fors_qc_keyword_to_paf(header, "ESO DET WIN1 BINY", NULL,
01889                                       "Binning factor along Y", instrume))
01890                 fors_calib_exit("Missing keyword ESO DET WIN1 BINY "
01891                                   "in frame header");
01892 
01893             if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
01894                                       "Archive name of input data", 
01895                                       instrume))
01896                 fors_calib_exit("Missing keyword ARCFILE in arc "
01897                                 "lamp header");
01898 
01899             pipefile = dfs_generate_filename(wavelength_map_tag);
01900             if (fors_qc_write_string("PIPEFILE", pipefile, 
01901                                     "Pipeline product name", instrume))
01902                 fors_calib_exit("Cannot write PIPEFILE to QC log file");
01903             cpl_free(pipefile); pipefile = NULL;
01904 
01905 
01906             /*
01907              * QC1 parameters
01908              */
01909 
01910             if (fors_qc_write_qc_double(header,
01911                                        mean_rms,
01912                                        "QC.WAVE.ACCURACY",
01913                                        "pixel",
01914                                        "Mean accuracy of wavecalib model",
01915                                        instrume)) {
01916                 fors_calib_exit("Cannot write mean wavelength calibration "
01917                                 "accuracy to QC log file");
01918             }
01919 
01920             if (fors_qc_write_qc_double(header,
01921                                        mean_rms_err,
01922                                        "QC.WAVE.ACCURACY.ERROR",
01923                                        "pixel",
01924                                        "Error on accuracy of wavecalib model",
01925                                        instrume)) {
01926                 fors_calib_exit("Cannot write error on wavelength calibration "
01927                                 "accuracy to QC log file");
01928             }
01929 
01930             if (compute_central_wave) {
01931 
01932                 data = cpl_image_get_data(wavemap);
01933     
01934                 if (lss) {
01935                     if (fors_qc_write_qc_double(header, 
01936                                                data[nx/2 + ccd_ysize*nx/2],
01937                                                "QC.LSS.CENTRAL.WAVELENGTH",
01938                                                "Angstrom", 
01939                                                "Wavelength at CCD center",
01940                                                instrume)) {
01941                         fors_calib_exit("Cannot write central wavelength to QC "
01942                                         "log file");
01943                     }
01944                 }
01945                 else {
01946                     if (fors_qc_write_qc_double(header, 
01947                                                data[nx/2 + ccd_ysize*nx/2],
01948                                                "QC.MOS.CENTRAL.WAVELENGTH",
01949                                                "Angstrom", 
01950                                                "Wavelength at CCD center",
01951                                                instrume)) {
01952                         fors_calib_exit("Cannot write central wavelength to QC "
01953                                         "log file");
01954                     }
01955                 }
01956             }
01957 
01958             fors_qc_end_group();
01959 
01960         }
01961 
01962         if (dfs_save_image(frameset, wavemap, wavelength_map_tag, header,
01963                            parlist, recipe, version))
01964             fors_calib_exit(NULL);
01965 
01966         cpl_image_delete(wavemap); wavemap = NULL;
01967 
01968         cpl_propertylist_erase_regexp(header, "^ESO QC ", 0);
01969 
01970         if (check) {
01971             cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01972             cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01973             /* cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01974             cpl_propertylist_update_double(header, "CD1_1", 1.0);
01975             cpl_propertylist_update_double(header, "CD1_2", 0.0);
01976             cpl_propertylist_update_double(header, "CD2_1", 0.0);
01977             cpl_propertylist_update_double(header, "CD2_2", 1.0);
01978             cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01979             cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01980 
01981             if (dfs_save_image(frameset, residual, disp_residuals_tag, header,
01982                                parlist, recipe, version))
01983                 fors_calib_exit(NULL);
01984 
01985             cpl_image_delete(residual); residual = NULL;
01986         }
01987 
01988         cpl_propertylist_delete(header); header = NULL;
01989         cpl_free(instrume); instrume = NULL;
01990 
01991         return 0;         /* Successful LSS data reduction */
01992 
01993     }   /* End of LSS data reduction section */
01994 
01995 
01996     /*
01997      * Here the MOS and MXU calibration is carried out.
01998      */
01999 
02000     /*
02001      * Detecting spectra on the CCD
02002      */
02003 
02004     cpl_msg_indent_less();
02005     cpl_msg_info(recipe, "Detecting spectra on CCD...");
02006     cpl_msg_indent_more();
02007 
02008     ccd_xsize = nx = cpl_image_get_size_x(spectra);
02009     ccd_ysize = ny = cpl_image_get_size_y(spectra);
02010 
02011     refmask = cpl_mask_new(nx, ny);
02012 
02013     checkwave = mos_wavelength_calibration_raw(spectra, lines, dispersion, 
02014                                                peakdetection, wradius, 
02015                                                wdegree, wreject, reference,
02016                                                &startwavelength, &endwavelength,
02017                                                NULL, NULL, NULL, NULL, NULL, 
02018                                                NULL, refmask);
02019 
02020     if (checkwave == NULL)
02021         fors_calib_exit("Wavelength calibration failure.");
02022 
02023     /*
02024      * Save check image to disk
02025      */
02026 
02027     header = cpl_propertylist_new();
02028     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
02029     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
02030     cpl_propertylist_update_double(header, "CRVAL1", 
02031                                    startwavelength + dispersion/2);
02032     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
02033     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
02034     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
02035     cpl_propertylist_update_double(header, "CD1_1", dispersion);
02036     cpl_propertylist_update_double(header, "CD1_2", 0.0);
02037     cpl_propertylist_update_double(header, "CD2_1", 0.0);
02038     cpl_propertylist_update_double(header, "CD2_2", 1.0);
02039     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
02040     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
02041 
02042     if (check) {
02043         if (dfs_save_image(frameset, checkwave, spectra_detection_tag, header, 
02044                            parlist, recipe, version))
02045             fors_calib_exit(NULL);
02046     }
02047 
02048     cpl_image_delete(checkwave); checkwave = NULL;
02049     cpl_propertylist_delete(header); header = NULL;
02050 
02051     cpl_msg_info(recipe, "Locate slits at reference wavelength on CCD...");
02052     slits = mos_locate_spectra(refmask);
02053 
02054     if (!slits) {
02055         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
02056         fors_calib_exit("No slits could be detected!");
02057     }
02058 
02059     refimage = cpl_image_new_from_mask(refmask);
02060     cpl_mask_delete(refmask); refmask = NULL;
02061 
02062     if (check) {
02063         save_header = dfs_load_header(frameset, arc_tag, 0);
02064         if (dfs_save_image(frameset, refimage, slit_map_tag, save_header,
02065                            parlist, recipe, version))
02066             fors_calib_exit(NULL);
02067         cpl_propertylist_delete(save_header); save_header = NULL;
02068     }
02069 
02070     cpl_image_delete(refimage); refimage = NULL;
02071 
02072     if (slit_ident) {
02073 
02074         /*
02075          * Attempt slit identification: this recipe may continue even
02076          * in case of failed identification (i.e., the position table is 
02077          * not produced, but an error is not set). In case of failure,
02078          * the spectra would be still extracted, even if they would not
02079          * be associated to slits on the mask.
02080          * 
02081          * The reason for making the slit identification an user option 
02082          * (via the parameter slit_ident) is to offer the possibility 
02083          * to avoid identifications that are only apparently successful, 
02084          * as it would happen in the case of an incorrect slit description 
02085          * in the data header.
02086          */
02087 
02088         cpl_msg_indent_less();
02089         cpl_msg_info(recipe, "Attempt slit identification (optional)...");
02090         cpl_msg_indent_more();
02091 
02092         positions = mos_identify_slits(slits, maskslits, NULL);
02093 
02094         if (positions) {
02095             cpl_table_delete(slits);
02096             slits = positions;
02097 
02098             /*
02099              * Eliminate slits which are _entirely_ outside the CCD
02100              */
02101 
02102             cpl_table_and_selected_double(slits, 
02103                                           "ybottom", CPL_GREATER_THAN, ny-1);
02104             cpl_table_or_selected_double(slits, 
02105                                           "ytop", CPL_LESS_THAN, 0);
02106             cpl_table_erase_selected(slits);
02107 
02108             nslits = cpl_table_get_nrow(slits);
02109 
02110             if (nslits == 0)
02111                 fors_calib_exit("No slits found on the CCD");
02112 
02113             cpl_msg_info(recipe, "%d slits are entirely or partially "
02114                          "contained in CCD", nslits);
02115 
02116         }
02117         else {
02118             slit_ident = 0;
02119             cpl_msg_info(recipe, "Global distortion model cannot be computed");
02120             if (cpl_error_get_code() != CPL_ERROR_NONE) {
02121                 fors_calib_exit(NULL);
02122             }
02123         }
02124     }
02125 
02126 
02127     /*
02128      * Determination of spectral curvature
02129      */
02130 
02131     cpl_msg_indent_less();
02132     cpl_msg_info(recipe, "Determining spectral curvature...");
02133     cpl_msg_indent_more();
02134 
02135     cpl_msg_info(recipe, "Tracing master flat field spectra edges...");
02136     traces = mos_trace_flat(master_flat, slits, reference, 
02137                             startwavelength, endwavelength, dispersion);
02138 
02139     if (!traces)
02140         fors_calib_exit("Tracing failure");
02141 
02142     cpl_msg_info(recipe, "Fitting flat field spectra edges...");
02143     polytraces = mos_poly_trace(slits, traces, cdegree);
02144 
02145     if (!polytraces)
02146         fors_calib_exit("Trace fitting failure");
02147 
02148     if (cmode) {
02149         cpl_msg_info(recipe, "Computing global spectral curvature model...");
02150         mos_global_trace(slits, polytraces, cmode);
02151     }
02152 
02153     if (dfs_save_table(frameset, traces, curv_traces_tag, NULL, parlist,
02154                        recipe, version))
02155         fors_calib_exit(NULL);
02156 
02157     cpl_table_delete(traces); traces = NULL;
02158 
02159     coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
02160     spatial = mos_spatial_calibration(spectra, slits, polytraces, reference, 
02161                                       startwavelength, endwavelength, 
02162                                       dispersion, 0, coordinate);
02163 
02164     if (!slit_ident) {
02165         cpl_image_delete(spectra); spectra = NULL;
02166     }
02167 
02168     /*
02169      * Flat field normalisation is done directly on the master flat
02170      * field (without spatial rectification first). The spectral
02171      * curvature model may be provided in input, in future releases.
02172      */
02173 
02174     cpl_msg_indent_less();
02175     cpl_msg_info(recipe, "Perform flat field normalisation...");
02176     cpl_msg_indent_more();
02177 
02178     norm_flat = cpl_image_duplicate(master_flat);
02179 
02180     smo_flat = mos_normalise_flat(norm_flat, coordinate, slits, polytraces, 
02181                                   reference, startwavelength, endwavelength,
02182                                   dispersion, dradius, ddegree);
02183 
02184     cpl_image_delete(smo_flat); smo_flat = NULL;  /* It may be a product */
02185 
02186  
02187     save_header = dfs_load_header(frameset, flat_tag, 0);
02188     cpl_propertylist_update_int(save_header, "ESO PRO DATANCOM", nflats);
02189 //%%%
02190     rect_flat = mos_spatial_calibration(master_flat, slits, polytraces, 
02191                                         reference, startwavelength, 
02192                                         endwavelength, dispersion, 0, NULL);
02193     rect_nflat = mos_spatial_calibration(norm_flat, slits, polytraces, 
02194                                         reference, startwavelength, 
02195                                         endwavelength, dispersion, 0, NULL);
02196 //%%%
02197 
02198     if (dfs_save_image(frameset, master_flat, master_screen_flat_tag,
02199                        save_header, parlist, recipe, version))
02200         fors_calib_exit(NULL);
02201 
02202     cpl_image_delete(master_flat); master_flat = NULL;
02203 
02204     if (dfs_save_image(frameset, norm_flat, master_norm_flat_tag,
02205                        save_header, parlist, recipe, version))
02206         fors_calib_exit(NULL);
02207 
02208     cpl_image_delete(norm_flat); norm_flat = NULL;
02209     cpl_propertylist_delete(save_header); save_header = NULL;
02210 
02211 
02212     /*
02213      * Final wavelength calibration of spectra having their curvature
02214      * removed
02215      */
02216 
02217     cpl_msg_indent_less();
02218     cpl_msg_info(recipe, "Perform final wavelength calibration...");
02219     cpl_msg_indent_more();
02220 
02221     nx = cpl_image_get_size_x(spatial);
02222     ny = cpl_image_get_size_y(spatial);
02223 
02224     idscoeff = cpl_table_new(ny);
02225     restable = cpl_table_new(nlines);
02226     rainbow = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
02227     if (check)
02228         residual = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
02229     fiterror = cpl_calloc(ny, sizeof(double));
02230     fitlines = cpl_calloc(ny, sizeof(int));
02231 
02232     rectified = mos_wavelength_calibration_final(spatial, slits, lines, 
02233                                                  dispersion, peakdetection, 
02234                                                  wradius, wdegree, wreject,
02235                                                  reference, &startwavelength, 
02236                                                  &endwavelength, fitlines, 
02237                                                  fiterror, idscoeff, rainbow, 
02238                                                  residual, restable);
02239 
02240 /*
02241 dfs_save_image(frameset, rainbow, "rainbow_calib", NULL, parlist, recipe, version);
02242 */
02243 
02244     cpl_image_delete(spatial); spatial = NULL;
02245 
02246     if (rectified == NULL)
02247         fors_calib_exit("Wavelength calibration failure.");
02248 
02249     if (dfs_save_table(frameset, restable, disp_residuals_table_tag, NULL,
02250                        parlist, recipe, version))
02251         fors_calib_exit(NULL);
02252 
02253     cpl_table_delete(restable); restable = NULL;
02254 
02255     cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
02256     cpl_table_set_column_unit(idscoeff, "error", "pixel");
02257     cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;
02258 
02259     for (i = 0; i < ny; i++)
02260         if (!cpl_table_is_valid(idscoeff, "c0", i))
02261             cpl_table_set_invalid(idscoeff, "error", i);
02262 
02263     delta = mos_map_pixel(idscoeff, reference, startwavelength,
02264                           endwavelength, dispersion, 2);
02265 
02266     header = cpl_propertylist_new();
02267     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
02268     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
02269     cpl_propertylist_update_double(header, "CRVAL1",
02270                                    startwavelength + dispersion/2);
02271     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
02272     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
02273     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
02274     cpl_propertylist_update_double(header, "CD1_1", dispersion);
02275     cpl_propertylist_update_double(header, "CD1_2", 0.0);
02276     cpl_propertylist_update_double(header, "CD2_1", 0.0);
02277     cpl_propertylist_update_double(header, "CD2_2", 1.0);
02278     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
02279     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
02280 
02281     if (dfs_save_image(frameset, delta, delta_image_tag,
02282                        header, parlist, recipe, version))
02283         fors_calib_exit(NULL);
02284 
02285     cpl_image_delete(delta); delta = NULL;
02286     cpl_propertylist_delete(header); header = NULL;
02287 
02288     mean_rms = mos_distortions_rms(rectified, lines, startwavelength, 
02289                                    dispersion, 6, 0);
02290 
02291     cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);
02292 
02293     mean_rms = cpl_table_get_column_mean(idscoeff, "error");
02294     mean_rms_err = cpl_table_get_column_stdev(idscoeff, "error");
02295 
02296     cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)", 
02297                  mean_rms, mean_rms * dispersion);
02298 
02299     restab = mos_resolution_table(rectified, startwavelength, dispersion, 
02300                                   60000, lines);
02301 
02302     if (restab) {
02303         cpl_msg_info(recipe, "Mean spectral resolution: %.2f", 
02304                    cpl_table_get_column_mean(restab, "resolution"));
02305         cpl_msg_info(recipe, "Mean reference lines FWHM: %.2f +/- %.2f pixel",
02306                    cpl_table_get_column_mean(restab, "fwhm") / dispersion,
02307                    cpl_table_get_column_mean(restab, "fwhm_rms") / dispersion);
02308 
02309         if (qc) {
02310 
02311             header = dfs_load_header(frameset, arc_tag, 0);
02312 
02313             if (header == NULL)
02314                 fors_calib_exit("Cannot reload arc lamp header");
02315 
02316             qclist = cpl_propertylist_new();
02317 
02318             fors_qc_start_group(qclist, "2.0", instrume);
02319 
02320 
02321             /*
02322              * QC1 group header
02323              */
02324 
02325             if (fors_qc_write_string("PRO.CATG", spectral_resolution_tag,
02326                                     "Product category", instrume))
02327                 fors_calib_exit("Cannot write product category to "
02328                                 "QC log file");
02329 
02330             if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL,
02331                                       "DPR type", instrume))
02332                 fors_calib_exit("Missing keyword DPR TYPE in arc "
02333                                 "lamp header");
02334 
02335             if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL,
02336                                       "Template", instrume))
02337                 fors_calib_exit("Missing keyword TPL ID in arc "
02338                                 "lamp header");
02339 
02340             if (fors_qc_keyword_to_paf(header, "ESO INS GRIS1 NAME", NULL,
02341                                       "Grism name", instrume))
02342                 fors_calib_exit("Missing keyword INS GRIS1 NAME in arc "
02343                                 "lamp header");
02344 
02345             if (fors_qc_keyword_to_paf(header, "ESO INS GRIS1 ID", NULL,
02346                                       "Grism identifier", instrume))
02347                 fors_calib_exit("Missing keyword INS GRIS1 ID in arc "
02348                                 "lamp header");
02349 
02350             if (cpl_propertylist_has(header, "ESO INS FILT1 NAME"))
02351                 fors_qc_keyword_to_paf(header, "ESO INS FILT1 NAME", NULL,
02352                                       "Filter name", instrume);
02353 
02354             if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
02355                                       "Collimator name", instrume))
02356                 fors_calib_exit("Missing keyword INS COLL NAME in arc "
02357                                 "lamp header");
02358 
02359             if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
02360                                       "Chip identifier", instrume))
02361                 fors_calib_exit("Missing keyword DET CHIP1 ID in arc "
02362                                 "lamp header");
02363 
02364             if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
02365                                       "Archive name of input data", 
02366                                       instrume))
02367                 fors_calib_exit("Missing keyword ARCFILE in arc "
02368                                 "lamp header");
02369 
02370             cpl_propertylist_delete(header); header = NULL;
02371 
02372             pipefile = dfs_generate_filename_tfits(spectral_resolution_tag);
02373             if (fors_qc_write_string("PIPEFILE", pipefile,
02374                                     "Pipeline product name", instrume))
02375                 fors_calib_exit("Cannot write PIPEFILE to QC log file");
02376             cpl_free(pipefile); pipefile = NULL;
02377 
02378 
02379             /*
02380              * QC1 parameters
02381              */
02382 
02383             if (mos)
02384                 keyname = "QC.MOS.RESOLUTION";
02385             else
02386                 keyname = "QC.MXU.RESOLUTION";
02387 
02388             if (fors_qc_write_qc_double(qclist, 
02389                                        cpl_table_get_column_mean(restab,
02390                                                                  "resolution"),
02391                                        keyname,
02392                                        "Angstrom",
02393                                        "Mean spectral resolution",
02394                                        instrume)) {
02395                 fors_calib_exit("Cannot write mean spectral resolution to QC "
02396                                 "log file");
02397             }
02398 
02399             if (mos)
02400                 keyname = "QC.MOS.RESOLUTION.RMS";
02401             else
02402                 keyname = "QC.MXU.RESOLUTION.RMS";
02403 
02404             if (fors_qc_write_qc_double(qclist, 
02405                                        cpl_table_get_column_stdev(restab, 
02406                                                                   "resolution"),
02407                                        keyname,
02408                                        "Angstrom", 
02409                                        "Scatter of spectral resolution",
02410                                        instrume)) {
02411                 fors_calib_exit("Cannot write spectral resolution scatter "
02412                                 "to QC log file");
02413             }
02414 
02415             if (mos)
02416                 keyname = "QC.MOS.RESOLUTION.NWAVE";
02417             else
02418                 keyname = "QC.MXU.RESOLUTION.NWAVE";
02419 
02420             if (fors_qc_write_qc_int(qclist, cpl_table_get_nrow(restab) -
02421                                     cpl_table_count_invalid(restab, 
02422                                                             "resolution"),
02423                                     keyname,
02424                                     NULL,
02425                                     "Number of examined wavelengths "
02426                                     "for resolution computation",
02427                                     instrume)) {
02428                 fors_calib_exit("Cannot write number of lines used in "
02429                                 "spectral resolution computation "
02430                                 "to QC log file");
02431             }
02432 
02433             if (mos)
02434                 keyname = "QC.MOS.RESOLUTION.MEANRMS";
02435             else
02436                 keyname = "QC.MXU.RESOLUTION.MEANRMS";
02437 
02438             if (fors_qc_write_qc_double(qclist,
02439                                        cpl_table_get_column_mean(restab,
02440                                                      "resolution_rms"),
02441                                     keyname, NULL,
02442                                     "Mean error on spectral "
02443                                     "resolution computation",
02444                                     instrume)) {
02445                 fors_calib_exit("Cannot write mean error in "
02446                                 "spectral resolution computation "
02447                                 "to QC log file");
02448             }
02449 
02450             if (mos)
02451                 keyname = "QC.MOS.RESOLUTION.NLINES";
02452             else
02453                 keyname = "QC.MXU.RESOLUTION.NLINES";
02454 
02455             if (fors_qc_write_qc_int(qclist,
02456                      cpl_table_get_column_mean(restab, "nlines") *
02457                                        cpl_table_get_nrow(restab),
02458                                     keyname, NULL,
02459                                     "Number of lines for spectral "
02460                                     "resolution computation",
02461                                     instrume)) {
02462                 fors_calib_exit("Cannot write number of examined "
02463                      "wavelengths in spectral resolution computation "
02464                      "to QC log file");
02465             }
02466 
02467             fors_qc_end_group();
02468 
02469         }
02470 
02471         if (dfs_save_table(frameset, restab, spectral_resolution_tag, qclist,
02472                            parlist, recipe, version))
02473             fors_calib_exit(NULL);
02474 
02475         cpl_table_delete(restab); restab = NULL;
02476         cpl_propertylist_delete(qclist); qclist = NULL;
02477 
02478     }
02479     else
02480         fors_calib_exit("Cannot compute the spectral resolution table");
02481 
02482     cpl_vector_delete(lines); lines = NULL;
02483 
02484     if (dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL,
02485                        parlist, recipe, version))
02486         fors_calib_exit(NULL);
02487 //%%%
02488 
02489     mapped_flat = mos_wavelength_calibration(rect_flat, reference,
02490                                       startwavelength, endwavelength,
02491                                       dispersion, idscoeff, 0);
02492 
02493     mapped_nflat = mos_wavelength_calibration(rect_nflat, reference,
02494                                       startwavelength, endwavelength,
02495                                       dispersion, idscoeff, 0);
02496 
02497     cpl_image_delete(rect_flat); rect_flat = NULL;
02498     cpl_image_delete(rect_nflat); rect_nflat = NULL;
02499 //%%%
02500 
02501     /*
02502      * Global distortion models
02503      */
02504 
02505     if (slit_ident) {
02506 
02507         cpl_msg_info(recipe, "Computing global distortions model");
02508         global = mos_global_distortion(slits, maskslits, idscoeff, 
02509                                        polytraces, reference);
02510 
02511         if (global && 0) {
02512             cpl_table *stest;
02513             cpl_table *ctest;
02514             cpl_table *dtest;
02515             cpl_image *itest;
02516 
02517             stest = mos_build_slit_location(global, maskslits, ccd_ysize);
02518 
02519             ctest = mos_build_curv_coeff(global, maskslits, stest);
02520             if (dfs_save_table(frameset, ctest, "CURVS", NULL,
02521                                parlist, recipe, version))
02522                 fors_calib_exit(NULL);
02523 
02524             itest = mos_spatial_calibration(spectra, stest, ctest, 
02525                                             reference, startwavelength, 
02526                                             endwavelength, dispersion, 
02527                                             0, NULL);
02528             cpl_table_delete(ctest); ctest = NULL;
02529             cpl_image_delete(itest); itest = NULL;
02530             if (dfs_save_table(frameset, stest, "SLITS", NULL,
02531                                parlist, recipe, version))
02532                 fors_calib_exit(NULL);
02533 
02534             dtest = mos_build_disp_coeff(global, stest);
02535             if (dfs_save_table(frameset, dtest, "DISPS", NULL,
02536                                parlist, recipe, version))
02537                 fors_calib_exit(NULL);
02538 
02539             cpl_table_delete(dtest); dtest = NULL;
02540             cpl_table_delete(stest); stest = NULL;
02541         }
02542 
02543         if (global) {
02544             if (dfs_save_table(frameset, global, global_distortion_tag, NULL,
02545                                parlist, recipe, version))
02546                 fors_calib_exit(NULL);
02547             cpl_table_delete(global); global = NULL;
02548         }
02549 
02550         cpl_image_delete(spectra); spectra = NULL;
02551         cpl_table_delete(maskslits); maskslits = NULL;
02552     }
02553 
02554     cpl_table_delete(idscoeff); idscoeff = NULL;
02555 
02556     header = cpl_propertylist_new();
02557     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
02558     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
02559     cpl_propertylist_update_double(header, "CRVAL1", 
02560                                    startwavelength + dispersion/2);
02561     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
02562     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
02563     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
02564     cpl_propertylist_update_double(header, "CD1_1", dispersion);
02565     cpl_propertylist_update_double(header, "CD1_2", 0.0);
02566     cpl_propertylist_update_double(header, "CD2_1", 0.0);
02567     cpl_propertylist_update_double(header, "CD2_2", 1.0);
02568     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
02569     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
02570     cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);
02571 
02572     if (dfs_save_image(frameset, rectified, reduced_lamp_tag, header,
02573                        parlist, recipe, version))
02574         fors_calib_exit(NULL);
02575 
02576     cpl_image_delete(rectified); rectified = NULL;
02577 //%%%
02578     cpl_propertylist_update_int(header, "ESO PRO DATANCOM", nflats);
02579 
02580     if (dfs_save_image(frameset, mapped_flat, mapped_screen_flat_tag, header,
02581                        parlist, recipe, version))
02582         fors_calib_exit(NULL);
02583 
02584     cpl_image_delete(mapped_flat); mapped_flat = NULL;
02585 
02586     if (dfs_save_image(frameset, mapped_nflat, mapped_norm_flat_tag, header,
02587                        parlist, recipe, version))
02588         fors_calib_exit(NULL);
02589 
02590     cpl_image_delete(mapped_nflat); mapped_nflat = NULL;
02591 
02592     cpl_propertylist_delete(header); header = NULL;
02593 
02594     if (check) {
02595         save_header = dfs_load_header(frameset, arc_tag, 0);
02596 
02597         cpl_propertylist_update_double(save_header, "CRPIX2", 1.0);
02598         cpl_propertylist_update_double(save_header, "CRVAL2", 1.0);
02599         /* cpl_propertylist_update_double(save_header, "CDELT2", 1.0); */
02600         cpl_propertylist_update_double(save_header, "CD1_1", 1.0);
02601         cpl_propertylist_update_double(save_header, "CD1_2", 0.0);
02602         cpl_propertylist_update_double(save_header, "CD2_1", 0.0);
02603         cpl_propertylist_update_double(save_header, "CD2_2", 1.0);
02604         cpl_propertylist_update_string(save_header, "CTYPE1", "LINEAR");
02605         cpl_propertylist_update_string(save_header, "CTYPE2", "PIXEL");
02606 
02607         if (dfs_save_image(frameset, residual, disp_residuals_tag, save_header,
02608                            parlist, recipe, version))
02609             fors_calib_exit(NULL);
02610 
02611         cpl_image_delete(residual); residual = NULL;
02612         cpl_propertylist_delete(save_header); save_header = NULL;
02613     }
02614 
02615     wavemap = mos_map_wavelengths(coordinate, rainbow, slits, polytraces, 
02616                                   reference, startwavelength, endwavelength, 
02617                                   dispersion);
02618 
02619     cpl_image_delete(rainbow); rainbow = NULL;
02620 
02621     save_header = dfs_load_header(frameset, arc_tag, 0);
02622 
02623     if (qc) {
02624         fors_qc_start_group(save_header, "2.0", instrume);
02625 
02626         /*
02627          * QC1 group header
02628          */
02629 
02630         if (fors_qc_write_string("PRO.CATG", wavelength_map_tag,
02631                                 "Product category", instrume))
02632             fors_calib_exit("Cannot write product category to "
02633                             "QC log file");
02634 
02635         if (fors_qc_keyword_to_paf(save_header, "ESO DPR TYPE", NULL,
02636                                   "DPR type", instrume))
02637             fors_calib_exit("Missing keyword DPR TYPE in arc "
02638                             "lamp header");
02639 
02640         if (fors_qc_keyword_to_paf(save_header, "ESO TPL ID", NULL,
02641                                   "Template", instrume))
02642             fors_calib_exit("Missing keyword TPL ID in arc "
02643                             "lamp header");
02644 
02645         if (fors_qc_keyword_to_paf(save_header, "ESO INS GRIS1 NAME", NULL,
02646                                   "Grism name", instrume))
02647             fors_calib_exit("Missing keyword INS GRIS1 NAME in arc "
02648                             "lamp header");
02649 
02650         if (fors_qc_keyword_to_paf(save_header, "ESO INS GRIS1 ID", NULL,
02651                                   "Grism identifier", instrume))
02652             fors_calib_exit("Missing keyword INS GRIS1 ID in arc "
02653                             "lamp header");
02654 
02655         if (cpl_propertylist_has(save_header, "ESO INS FILT1 NAME"))
02656             fors_qc_keyword_to_paf(save_header, "ESO INS FILT1 NAME", NULL,
02657                                   "Filter name", instrume);
02658 
02659         if (fors_qc_keyword_to_paf(save_header, "ESO INS COLL NAME", NULL,
02660                                   "Collimator name", instrume))
02661             fors_calib_exit("Missing keyword INS COLL NAME in arc "
02662                             "lamp header");
02663 
02664         if (fors_qc_keyword_to_paf(save_header, "ESO DET CHIP1 ID", NULL,
02665                                   "Chip identifier", instrume))
02666             fors_calib_exit("Missing keyword DET CHIP1 ID in arc "
02667                             "lamp header");
02668 
02669         if (fors_qc_keyword_to_paf(save_header, "ESO DET WIN1 BINX", NULL,
02670                                   "Binning factor along X", instrume))
02671             fors_calib_exit("Missing keyword ESO DET WIN1 BINX "
02672                               "in frame header");
02673 
02674         if (fors_qc_keyword_to_paf(save_header, "ESO DET WIN1 BINY", NULL,
02675                                   "Binning factor along Y", instrume))
02676             fors_calib_exit("Missing keyword ESO DET WIN1 BINY "
02677                               "in frame header");
02678 
02679         if (fors_qc_keyword_to_paf(save_header, "ARCFILE", NULL,
02680                                   "Archive name of input data",
02681                                   instrume))
02682             fors_calib_exit("Missing keyword ARCFILE in arc "
02683                             "lamp header");
02684 
02685         pipefile = dfs_generate_filename(wavelength_map_tag);
02686         if (fors_qc_write_string("PIPEFILE", pipefile,
02687                                 "Pipeline product name", instrume))
02688             fors_calib_exit("Cannot write PIPEFILE to QC log file");
02689         cpl_free(pipefile); pipefile = NULL;
02690 
02691 
02692         /*
02693          * QC1 parameters
02694          */
02695 
02696         if (fors_qc_write_qc_double(save_header,
02697                                    mean_rms,
02698                                    "QC.WAVE.ACCURACY",
02699                                    "pixel",
02700                                    "Mean accuracy of wavecalib model",
02701                                    instrume)) {
02702             fors_calib_exit("Cannot write mean wavelength calibration "
02703                             "accuracy to QC log file");
02704         }
02705 
02706 
02707         if (fors_qc_write_qc_double(save_header,
02708                                    mean_rms_err,
02709                                    "QC.WAVE.ACCURACY.ERROR",
02710                                    "pixel",
02711                                    "Error on accuracy of wavecalib model",
02712                                    instrume)) {
02713             fors_calib_exit("Cannot write error on wavelength calibration "
02714                             "accuracy to QC log file");
02715         }
02716 
02717         fors_qc_end_group();
02718 
02719     }
02720 
02721     cpl_free(instrume); instrume = NULL;
02722 
02723     if (dfs_save_image(frameset, wavemap, wavelength_map_tag, save_header,
02724                        parlist, recipe, version))
02725         fors_calib_exit(NULL);
02726 
02727     cpl_image_delete(wavemap); wavemap = NULL;
02728 
02729     cpl_propertylist_erase_regexp(save_header, "^ESO QC ", 0);
02730 
02731     if (dfs_save_image(frameset, coordinate, spatial_map_tag, save_header,
02732                        parlist, recipe, version))
02733         fors_calib_exit(NULL);
02734 
02735     cpl_image_delete(coordinate); coordinate = NULL;
02736     cpl_propertylist_delete(save_header); save_header = NULL;
02737 
02738     if (dfs_save_table(frameset, polytraces, curv_coeff_tag, NULL,
02739                        parlist, recipe, version))
02740         fors_calib_exit(NULL);
02741 
02742     cpl_table_delete(polytraces); polytraces = NULL;
02743 
02744     if (dfs_save_table(frameset, slits, slit_location_tag, NULL,
02745                        parlist, recipe, version))
02746         fors_calib_exit(NULL);
02747 
02748     cpl_table_delete(slits); slits = NULL;
02749 
02750     if (cpl_error_get_code()) {
02751         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
02752         fors_calib_exit(NULL);
02753     }
02754 
02755     return 0;
02756 }

Generated on Wed Sep 10 07:31:51 2008 for FORS Pipeline Reference Manual by  doxygen 1.4.6