vimos_calib.c

00001 /* $Id: vimos_calib.c,v 1.7 2008/12/10 14:31:39 lbilbao 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: lbilbao $
00023  * $Date: 2008/12/10 14:31:39 $
00024  * $Revision: 1.7 $
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 vimos_calib_create(cpl_plugin *);
00039 static int vimos_calib_exec(cpl_plugin *);
00040 static int vimos_calib_destroy(cpl_plugin *);
00041 static int vimos_calib(cpl_parameterlist *, cpl_frameset *);
00042 
00043 static char vimos_calib_description[] =
00044 "This recipe is used to identify reference lines on MOS 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 vimos_science. The input arc lamp and flat field exposures are\n"
00049 "assumed to be obtained quasi-simultaneously, so that they would be\n" 
00050 "described by exactly the same instrument distortions.\n"
00051 "A line catalog must be specified, containing the wavelengths of the\n"
00052 "reference arc lamp lines used for the wavelength calibration. A grism\n"
00053 "table (typically depending on the instrument mode, and in particular on\n"
00054 "the grism used) may also be specified: this table contains a default\n"
00055 "recipe parameter setting to control the way spectra are extracted for\n"
00056 "a specific instrument mode, as it is used for automatic run of the\n"
00057 "pipeline on Paranal and in Garching. If this table is specified, it\n"
00058 "will modify the default recipe parameter setting, with the exception of\n"
00059 "those parameters which have been explicitly modifyed on the command line.\n"
00060 "If a grism table is not specified, the input recipe parameters values\n"
00061 "will always be read from the command line, or from an esorex configuration\n"
00062 "file if present, or from their generic default values (that are rarely\n"
00063 "meaningful). Finally a master bias frame must be input to this recipe.\n" 
00064 "In the table below the MOS_CURV_COEFF, MOS_CURV_TRACES, MOS_SPATIAL_MAP\n"
00065 "MOS_ARC_SPECTRUM_EXTRACTED, MOS_SPECTRA_DETECTION, MOS_SLIT_MAP, and\n" 
00066 "MOS_SLIT_LOCATION, are never created in case of long-slit-like data.\n" 
00067 "The products MOS_SPECTRA_DETECTION, MOS_SLIT_MAP, and MOS_DISP_RESIDUALS,\n" 
00068 "are just created if the --check parameter is set to true. The product\n"
00069 "GLOBAL_DISTORTION_TABLE is just created if more than 12 separate spectra\n"
00070 "are found in the CCD.\n\n"
00071 "Input files:\n\n"
00072 "  DO category:               Type:       Explanation:         Required:\n"
00073 "  MOS_SCREEN_FLAT            Raw         Flat field exposures    Y\n"
00074 "  MOS_ARC_SPECTRUM           Raw         Arc lamp exposure       Y\n"
00075 "  MASTER_BIAS or BIAS        Calib       Bias frame              Y\n"
00076 "  LINE_CATALOG               Calib       Line catalog            Y\n"
00077 "  GRISM_TABLE                Calib       Grism table             .\n\n"
00078 "Output files:\n\n"
00079 "  DO category:               Data type:  Explanation:\n"
00080 "  MOS_COMBINED_SCREEN_FLAT   FITS image  Combined (sum) flat field\n"
00081 "  MOS_MASTER_SCREEN_FLAT     FITS image  Normalised flat field\n"
00082 "  MOS_ARC_SPECTRUM_EXTRACTED FITS image  Wavelength calibrated arc spectrum\n"
00083 "  MOS_DISP_COEFF             FITS table  Inverse dispersion coefficients\n"
00084 "  MOS_DISP_RESIDUALS         FITS image  Residuals in wavelength calibration\n"
00085 "  MOS_DISP_RESIDUALS_TABLE   FITS table  Residuals in wavelength calibration\n"
00086 "  MOS_DELTA_IMAGE            FITS image  Offset vs linear wavelength calib\n"
00087 "  MOS_WAVELENGTH_MAP         FITS image  Wavelength for each pixel on CCD\n"
00088 "  MOS_SPECTRA_DETECTION      FITS image  Check for preliminary detection\n"
00089 "  MOS_SLIT_MAP               FITS image  Map of central wavelength on CCD\n"
00090 "  MOS_CURV_TRACES            FITS table  Spectral curvature traces\n"
00091 "  MOS_CURV_COEFF             FITS table  Spectral curvature coefficients\n"
00092 "  MOS_SPATIAL_MAP            FITS image  Spatial position along slit on CCD\n"
00093 "  MOS_SPECTRAL_RESOLUTION    FITS table  Resolution at reference arc lines\n"
00094 "  MOS_SLIT_LOCATION          FITS table  Slits on product frames and CCD\n"
00095 "  GLOBAL_DISTORTION_TABLE    FITS table  Global distortions table\n\n";
00096 
00097 #define vimos_calib_exit(message)             \
00098 {                                             \
00099 if (message) cpl_msg_error(recipe, message);  \
00100 cpl_free(instrume);                           \
00101 cpl_free(pipefile);                           \
00102 cpl_free(fiterror);                           \
00103 cpl_free(fitlines);                           \
00104 cpl_image_delete(bias);                       \
00105 cpl_image_delete(master_bias);                \
00106 cpl_image_delete(coordinate);                 \
00107 cpl_image_delete(checkwave);                  \
00108 cpl_image_delete(flat);                       \
00109 cpl_image_delete(master_flat);                \
00110 cpl_image_delete(norm_flat);                  \
00111 cpl_image_delete(rainbow);                    \
00112 cpl_image_delete(rectified);                  \
00113 cpl_image_delete(residual);                   \
00114 cpl_image_delete(smo_flat);                   \
00115 cpl_image_delete(spatial);                    \
00116 cpl_image_delete(spectra);                    \
00117 cpl_image_delete(wavemap);                    \
00118 cpl_image_delete(delta);                      \
00119 cpl_mask_delete(refmask);                     \
00120 cpl_propertylist_delete(header);              \
00121 cpl_propertylist_delete(save_header);         \
00122 cpl_propertylist_delete(qclist);              \
00123 cpl_table_delete(grism_table);                \
00124 cpl_table_delete(idscoeff);                   \
00125 cpl_table_delete(restable);                   \
00126 cpl_table_delete(maskslits);                  \
00127 cpl_table_delete(overscans);                  \
00128 cpl_table_delete(traces);                     \
00129 cpl_table_delete(polytraces);                 \
00130 cpl_table_delete(slits);                      \
00131 cpl_table_delete(restab);                     \
00132 cpl_table_delete(global);                     \
00133 cpl_table_delete(wavelengths);                \
00134 cpl_vector_delete(lines);                     \
00135 cpl_msg_indent_less();                        \
00136 return -1;                                    \
00137 }
00138 
00139 #define vimos_calib_exit_memcheck(message)       \
00140 {                                               \
00141 if (message) cpl_msg_info(recipe, message);     \
00142 printf("free instrume (%p)\n", instrume);       \
00143 cpl_free(instrume);                             \
00144 printf("free pipefile (%p)\n", pipefile);       \
00145 cpl_free(pipefile);                             \
00146 printf("free fiterror (%p)\n", fiterror);       \
00147 cpl_free(fiterror);                             \
00148 printf("free fitlines (%p)\n", fitlines);       \
00149 cpl_free(fitlines);                             \
00150 printf("free bias (%p)\n", bias);               \
00151 cpl_image_delete(bias);                         \
00152 printf("free master_bias (%p)\n", master_bias); \
00153 cpl_image_delete(master_bias);                  \
00154 printf("free coordinate (%p)\n", coordinate);   \
00155 cpl_image_delete(coordinate);                   \
00156 printf("free checkwave (%p)\n", checkwave);     \
00157 cpl_image_delete(checkwave);                    \
00158 printf("free flat (%p)\n", flat);               \
00159 cpl_image_delete(flat);                         \
00160 printf("free master_flat (%p)\n", master_flat); \
00161 cpl_image_delete(master_flat);                  \
00162 printf("free norm_flat (%p)\n", norm_flat);     \
00163 cpl_image_delete(norm_flat);                    \
00164 printf("free rainbow (%p)\n", rainbow);         \
00165 cpl_image_delete(rainbow);                      \
00166 printf("free rectified (%p)\n", rectified);     \
00167 cpl_image_delete(rectified);                    \
00168 printf("free residual (%p)\n", residual);       \
00169 cpl_image_delete(residual);                     \
00170 printf("free smo_flat (%p)\n", smo_flat);       \
00171 cpl_image_delete(smo_flat);                     \
00172 printf("free spatial (%p)\n", spatial);         \
00173 cpl_image_delete(spatial);                      \
00174 printf("free spectra (%p)\n", spectra);         \
00175 cpl_image_delete(spectra);                      \
00176 printf("free wavemap (%p)\n", wavemap);         \
00177 cpl_image_delete(wavemap);                      \
00178 printf("free delta (%p)\n", delta);             \
00179 cpl_image_delete(delta);                        \
00180 printf("free refmask (%p)\n", refmask);         \
00181 cpl_mask_delete(refmask);                       \
00182 printf("free header (%p)\n", header);           \
00183 cpl_propertylist_delete(header);                \
00184 printf("free save_header (%p)\n", save_header); \
00185 cpl_propertylist_delete(save_header);           \
00186 printf("free qclist (%p)\n", qclist);           \
00187 cpl_propertylist_delete(qclist);                \
00188 printf("free grism_table (%p)\n", grism_table); \
00189 cpl_table_delete(grism_table);                  \
00190 printf("free idscoeff (%p)\n", idscoeff);       \
00191 cpl_table_delete(idscoeff);                     \
00192 printf("free restable (%p)\n", restable);       \
00193 cpl_table_delete(restable);                     \
00194 printf("free maskslits (%p)\n", maskslits);     \
00195 cpl_table_delete(maskslits);                    \
00196 printf("free overscans (%p)\n", overscans);     \
00197 cpl_table_delete(overscans);                    \
00198 printf("free traces (%p)\n", traces);           \
00199 cpl_table_delete(traces);                       \
00200 printf("free polytraces (%p)\n", polytraces);   \
00201 cpl_table_delete(polytraces);                   \
00202 printf("free slits (%p)\n", slits);             \
00203 cpl_table_delete(slits);                        \
00204 printf("free restab (%p)\n", restab);           \
00205 cpl_table_delete(restab);                       \
00206 printf("free global (%p)\n", global);           \
00207 cpl_table_delete(global);                       \
00208 printf("free wavelengths (%p)\n", wavelengths); \
00209 cpl_table_delete(wavelengths);                  \
00210 printf("free lines (%p)\n", lines);             \
00211 cpl_vector_delete(lines);                       \
00212 cpl_msg_indent_less();                          \
00213 return 0;                                       \
00214 }
00215 
00216 
00228 int cpl_plugin_get_info(cpl_pluginlist *list)
00229 {
00230     cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe );
00231     cpl_plugin *plugin = &recipe->interface;
00232 
00233     cpl_plugin_init(plugin,
00234                     CPL_PLUGIN_API,
00235                     FORS_BINARY_VERSION,
00236                     CPL_PLUGIN_TYPE_RECIPE,
00237                     "vimos_calib",
00238                     "Determination of the extraction mask",
00239                     vimos_calib_description,
00240                     "Carlo Izzo",
00241                     PACKAGE_BUGREPORT,
00242     "This file is currently part of the FORS Instrument Pipeline\n"
00243     "Copyright (C) 2002-2006 European Southern Observatory\n\n"
00244     "This program is free software; you can redistribute it and/or modify\n"
00245     "it under the terms of the GNU General Public License as published by\n"
00246     "the Free Software Foundation; either version 2 of the License, or\n"
00247     "(at your option) any later version.\n\n"
00248     "This program is distributed in the hope that it will be useful,\n"
00249     "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
00250     "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
00251     "GNU General Public License for more details.\n\n"
00252     "You should have received a copy of the GNU General Public License\n"
00253     "along with this program; if not, write to the Free Software Foundation,\n"
00254     "Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n",
00255                     vimos_calib_create,
00256                     vimos_calib_exec,
00257                     vimos_calib_destroy);
00258 
00259     cpl_pluginlist_append(list, plugin);
00260     
00261     return 0;
00262 }
00263 
00264 
00275 static int vimos_calib_create(cpl_plugin *plugin)
00276 {
00277     cpl_recipe    *recipe;
00278     cpl_parameter *p;
00279 
00280 
00281     /* 
00282      * Check that the plugin is part of a valid recipe 
00283      */
00284 
00285     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00286         recipe = (cpl_recipe *)plugin;
00287     else 
00288         return -1;
00289 
00290     /* 
00291      * Create the parameters list in the cpl_recipe object 
00292      */
00293 
00294     recipe->parameters = cpl_parameterlist_new(); 
00295 
00296 
00297     /*
00298      * Dispersion
00299      */
00300 
00301     p = cpl_parameter_new_value("fors.vimos_calib.dispersion",
00302                                 CPL_TYPE_DOUBLE,
00303                                 "Expected spectral dispersion (Angstrom/pixel)",
00304                                 "fors.vimos_calib",
00305                                 0.0);
00306     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dispersion");
00307     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00308     cpl_parameterlist_append(recipe->parameters, p);
00309 
00310     /*
00311      * Peak detection level
00312      */
00313 
00314     p = cpl_parameter_new_value("fors.vimos_calib.peakdetection",
00315                                 CPL_TYPE_DOUBLE,
00316                                 "Initial peak detection threshold (ADU)",
00317                                 "fors.vimos_calib",
00318                                 0.0);
00319     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "peakdetection");
00320     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00321     cpl_parameterlist_append(recipe->parameters, p);
00322 
00323     /* 
00324      * Degree of wavelength calibration polynomial
00325      */
00326 
00327     p = cpl_parameter_new_value("fors.vimos_calib.wdegree",
00328                                 CPL_TYPE_INT,
00329                                 "Degree of wavelength calibration polynomial",
00330                                 "fors.vimos_calib",
00331                                 0);
00332     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wdegree");
00333     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00334     cpl_parameterlist_append(recipe->parameters, p);
00335 
00336     /*
00337      * Reference lines search radius
00338      */
00339 
00340     p = cpl_parameter_new_value("fors.vimos_calib.wradius",
00341                                 CPL_TYPE_INT,
00342                                 "Search radius if iterating pattern-matching "
00343                                 "with first-guess method",
00344                                 "fors.vimos_calib",
00345                                 4);
00346     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wradius");
00347     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00348     cpl_parameterlist_append(recipe->parameters, p);
00349 
00350     /*
00351      * Rejection threshold in dispersion relation polynomial fitting
00352      */
00353 
00354     p = cpl_parameter_new_value("fors.vimos_calib.wreject",
00355                                 CPL_TYPE_DOUBLE,
00356                                 "Rejection threshold in dispersion "
00357                                 "relation fit (pixel)",
00358                                 "fors.vimos_calib",
00359                                 0.7);
00360     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wreject");
00361     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00362     cpl_parameterlist_append(recipe->parameters, p);
00363 
00364     /*
00365      * Wavelength solution interpolation (for LSS data)
00366      */
00367 
00368     p = cpl_parameter_new_value("fors.vimos_calib.wmode",
00369                                 CPL_TYPE_INT,
00370                                 "Interpolation mode of wavelength solution "
00371                                 "applicable to LSS-like data (0 = no "
00372                                 "interpolation, 1 = fill gaps, 2 = global "
00373                                 "model)",
00374                                 "fors.vimos_calib",
00375                                 2);
00376     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wmode");
00377     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00378     cpl_parameterlist_append(recipe->parameters, p);
00379 
00380     /*
00381      * Line catalog table column containing the reference wavelengths
00382      */
00383 
00384     p = cpl_parameter_new_value("fors.vimos_calib.wcolumn",
00385                                 CPL_TYPE_STRING,
00386                                 "Name of line catalog table column "
00387                                 "with wavelengths",
00388                                 "fors.vimos_calib",
00389                                 "WLEN");
00390     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wcolumn");
00391     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00392     cpl_parameterlist_append(recipe->parameters, p);
00393 
00394     /*
00395      * Degree of spectral curvature polynomial
00396      */
00397 
00398     p = cpl_parameter_new_value("fors.vimos_calib.cdegree",
00399                                 CPL_TYPE_INT,
00400                                 "Degree of spectral curvature polynomial",
00401                                 "fors.vimos_calib",
00402                                 0);
00403     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cdegree");
00404     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00405     cpl_parameterlist_append(recipe->parameters, p);
00406 
00407     /*
00408      * Curvature solution interpolation (for MOS-like data)
00409      */
00410  
00411     p = cpl_parameter_new_value("fors.vimos_calib.cmode",
00412                                 CPL_TYPE_INT,
00413                                 "Interpolation mode of curvature solution "
00414                                 "applicable to MOS-like data (0 = no "
00415                                 "interpolation, 1 = fill gaps, 2 = global "
00416                                 "model)",
00417                                 "fors.vimos_calib",
00418                                 1);
00419     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cmode");
00420     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00421     cpl_parameterlist_append(recipe->parameters, p);
00422 
00423     /*
00424      * Start wavelength for spectral extraction
00425      */
00426 
00427     p = cpl_parameter_new_value("fors.vimos_calib.startwavelength",
00428                                 CPL_TYPE_DOUBLE,
00429                                 "Start wavelength in spectral extraction",
00430                                 "fors.vimos_calib",
00431                                 0.0);
00432     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "startwavelength");
00433     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00434     cpl_parameterlist_append(recipe->parameters, p);
00435 
00436     /*
00437      * End wavelength for spectral extraction
00438      */
00439 
00440     p = cpl_parameter_new_value("fors.vimos_calib.endwavelength",
00441                                 CPL_TYPE_DOUBLE,
00442                                 "End wavelength in spectral extraction",
00443                                 "fors.vimos_calib",
00444                                 0.0);
00445     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "endwavelength");
00446     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00447     cpl_parameterlist_append(recipe->parameters, p);
00448 
00449     /*
00450      * Reference wavelength for wavelength calibration
00451      */
00452 
00453     p = cpl_parameter_new_value("fors.vimos_calib.reference",
00454                                 CPL_TYPE_DOUBLE,
00455                                 "Reference wavelength for calibration",
00456                                 "fors.vimos_calib",
00457                                 0.0);
00458     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "reference");
00459     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00460     cpl_parameterlist_append(recipe->parameters, p);
00461 
00462     /*
00463      * Try slit identification
00464      */
00465 
00466     p = cpl_parameter_new_value("fors.vimos_calib.slit_ident",
00467                                 CPL_TYPE_BOOL,
00468                                 "Attempt slit identification",
00469                                 "fors.vimos_calib",
00470                                 TRUE);
00471     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "slit_ident");
00472     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00473     cpl_parameterlist_append(recipe->parameters, p);
00474 
00475     /*
00476      * Degree of flat field fitting polynomial along spatial direction 
00477      * (used for LSS data)
00478      */
00479 
00480     p = cpl_parameter_new_value("fors.vimos_calib.sdegree",
00481                                 CPL_TYPE_INT,
00482                                 "Degree of flat field fitting polynomial "
00483                                 "along spatial direction (used for LSS-like "
00484                                 "data only)",
00485                                 "fors.vimos_calib",
00486                                 4);
00487     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "sdegree");
00488     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00489     cpl_parameterlist_append(recipe->parameters, p);
00490 
00491     /*
00492      * Degree of flat field fitting polynomial along dispersion direction
00493      * (used for MOS data)
00494      */
00495 
00496     p = cpl_parameter_new_value("fors.vimos_calib.ddegree",
00497                                 CPL_TYPE_INT,
00498                                 "Degree of flat field fitting polynomial "
00499                                 "along dispersion direction (not used for "
00500                                 "long-slit-like data)",
00501                                 "fors.vimos_calib",
00502                                 -1);
00503     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ddegree");
00504     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00505     cpl_parameterlist_append(recipe->parameters, p);
00506 
00507     /*
00508      * Smooth box radius for flat field along dispersion direction
00509      */
00510 
00511     p = cpl_parameter_new_value("fors.vimos_calib.dradius",
00512                                 CPL_TYPE_INT,
00513                                 "Smooth box radius for flat field along "
00514                                 "dispersion direction",
00515                                 "fors.vimos_calib",
00516                                 10);
00517     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dradius");
00518     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00519     cpl_parameterlist_append(recipe->parameters, p);
00520 
00521     /*
00522      * Smooth box radius for flat field along spatial direction
00523      * (used for LSS data only)
00524      */
00525 
00526     p = cpl_parameter_new_value("fors.vimos_calib.sradius",
00527                                 CPL_TYPE_INT,
00528                                 "Smooth box radius for flat field along "
00529                                 "spatial direction",
00530                                 "fors.vimos_calib",
00531                                 10);
00532     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "sradius");
00533     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00534     cpl_parameterlist_append(recipe->parameters, p);
00535 
00536     /*
00537      * Computation of QC1 parameters
00538      */
00539 
00540     p = cpl_parameter_new_value("fors.vimos_calib.qc",
00541                                 CPL_TYPE_BOOL,
00542                                 "Compute QC1 parameters",
00543                                 "fors.vimos_calib",
00544                                 FALSE);
00545     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "qc");
00546     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00547     cpl_parameterlist_append(recipe->parameters, p);
00548 
00549     /*
00550      * Create check products
00551      */
00552 
00553     p = cpl_parameter_new_value("fors.vimos_calib.check",
00554                                 CPL_TYPE_BOOL,
00555                                 "Create intermediate products",
00556                                 "fors.vimos_calib",
00557                                 FALSE);
00558     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "check");
00559     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00560     cpl_parameterlist_append(recipe->parameters, p);
00561 
00562     return 0;
00563 }
00564 
00565 
00574 static int vimos_calib_exec(cpl_plugin *plugin)
00575 {
00576     cpl_recipe *recipe;
00577     
00578     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00579         recipe = (cpl_recipe *)plugin;
00580     else 
00581         return -1;
00582 
00583     return vimos_calib(recipe->parameters, recipe->frames);
00584 }
00585 
00586 
00595 static int vimos_calib_destroy(cpl_plugin *plugin)
00596 {
00597     cpl_recipe *recipe;
00598     
00599     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00600         recipe = (cpl_recipe *)plugin;
00601     else 
00602         return -1;
00603 
00604     cpl_parameterlist_delete(recipe->parameters); 
00605 
00606     return 0;
00607 }
00608 
00609 
00619 static int vimos_calib(cpl_parameterlist *parlist, cpl_frameset *frameset)
00620 {
00621 
00622     const char *recipe = "vimos_calib";
00623 
00624     /*
00625      * Input parameters
00626      */
00627 
00628     double      dispersion;
00629     double      peakdetection;
00630     int         wdegree;
00631     int         wradius;
00632     double      wreject;
00633     int         wmode;
00634     const char *wcolumn;
00635     int         cdegree;
00636     int         cmode;
00637     double      startwavelength;
00638     double      endwavelength;
00639     double      reference;
00640     int         slit_ident;
00641     int         sdegree;
00642     int         ddegree;
00643     int         sradius;
00644     int         dradius;
00645     int         qc;
00646     int         check;
00647 
00648     /*
00649      * CPL objects
00650      */
00651 
00652     cpl_imagelist    *biases      = NULL;
00653     cpl_image        *bias        = NULL;
00654     cpl_image        *master_bias = NULL;
00655     cpl_image        *multi_bias  = NULL;
00656     cpl_image        *flat        = NULL;
00657     cpl_image        *master_flat = NULL;
00658     cpl_image        *smo_flat    = NULL;
00659     cpl_image        *norm_flat   = NULL;
00660     cpl_image        *spectra     = NULL;
00661     cpl_image        *wavemap     = NULL;
00662     cpl_image        *delta       = NULL;
00663     cpl_image        *residual    = NULL;
00664     cpl_image        *checkwave   = NULL;
00665     cpl_image        *rectified   = NULL;
00666     cpl_image        *dummy       = NULL;
00667     cpl_image        *refimage    = NULL;
00668     cpl_image        *coordinate  = NULL;
00669     cpl_image        *rainbow     = NULL;
00670     cpl_image        *spatial     = NULL;
00671 
00672     cpl_mask         *refmask     = NULL;
00673 
00674     cpl_table        *grism_table = NULL;
00675     cpl_table        *overscans   = NULL;
00676     cpl_table        *wavelengths = NULL;
00677     cpl_table        *idscoeff    = NULL;
00678     cpl_table        *restable    = NULL;
00679     cpl_table        *slits       = NULL;
00680     cpl_table        *positions   = NULL;
00681     cpl_table        *maskslits   = NULL;
00682     cpl_table        *traces      = NULL;
00683     cpl_table        *polytraces  = NULL;
00684     cpl_table        *restab      = NULL;
00685     cpl_table        *global      = NULL;
00686 
00687     cpl_vector       *lines       = NULL;
00688 
00689     cpl_propertylist *header      = NULL;
00690     cpl_propertylist *save_header = NULL;
00691     cpl_propertylist *qclist      = NULL;
00692 
00693     /*
00694      * Auxiliary variables
00695      */
00696 
00697     char        version[80];
00698     const char *arc_tag;
00699     const char *flat_tag;
00700     const char *master_screen_flat_tag;
00701     const char *master_norm_flat_tag;
00702     const char *reduced_lamp_tag;
00703     const char *disp_residuals_tag;
00704     const char *disp_coeff_tag;
00705     const char *wavelength_map_tag;
00706     const char *spectra_detection_tag;
00707     const char *spectral_resolution_tag;
00708     const char *slit_map_tag;
00709     const char *curv_traces_tag;
00710     const char *curv_coeff_tag;
00711     const char *spatial_map_tag;
00712     const char *slit_location_tag;
00713     const char *global_distortion_tag = "GLOBAL_DISTORTION_TABLE";
00714     const char *disp_residuals_table_tag;
00715     const char *delta_image_tag;
00716     const char *keyname;
00717     const char *key_gris_name;
00718     const char *key_gris_id;
00719     const char *key_filt_name;
00720     const char *key_filt_id;
00721     int         quadrant;
00722     int         mos;
00723     int         treat_as_lss = 0;
00724     int         nslits;
00725     float      *data;
00726     double     *xpos;
00727     double      mxpos;
00728     double      mean_rms;
00729     double      alltime;
00730     int         nflats;
00731     int         nbias;
00732     int         nlines;
00733     double     *line;
00734     double     *fiterror = NULL;
00735     int        *fitlines = NULL;
00736     int         nx, ny;
00737     double      gain;
00738     int         compute_qc;
00739     int         ccd_xsize, ccd_ysize;
00740     int         rotate = 1;
00741     int         rotate_back = -1;
00742     int         i;
00743 
00744     char       *instrume = NULL;
00745     char       *pipefile = NULL;
00746 
00747 
00748     snprintf(version, 80, "%s-%s", PACKAGE, PACKAGE_VERSION);
00749 
00750     cpl_msg_set_indentation(2);
00751 
00752     if (dfs_files_dont_exist(frameset))
00753         vimos_calib_exit(NULL);
00754 
00755 
00756     /* 
00757      * Get configuration parameters
00758      */
00759 
00760     cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
00761     cpl_msg_indent_more();
00762 
00763     if (cpl_frameset_count_tags(frameset, "GRISM_TABLE") > 1)
00764         vimos_calib_exit("Too many in input: GRISM_TABLE");
00765 
00766     grism_table = dfs_load_table(frameset, "GRISM_TABLE", 1);
00767 
00768     dispersion = dfs_get_parameter_double(parlist, 
00769                     "fors.vimos_calib.dispersion", grism_table);
00770 
00771     if (dispersion <= 0.0)
00772         vimos_calib_exit("Invalid spectral dispersion value");
00773 
00774     peakdetection = dfs_get_parameter_double(parlist, 
00775                     "fors.vimos_calib.peakdetection", grism_table);
00776     if (peakdetection <= 0.0)
00777         vimos_calib_exit("Invalid peak detection level");
00778 
00779     wdegree = dfs_get_parameter_int(parlist, 
00780                     "fors.vimos_calib.wdegree", grism_table);
00781 
00782     if (wdegree < 1)
00783         vimos_calib_exit("Invalid polynomial degree");
00784 
00785     if (wdegree > 5)
00786         vimos_calib_exit("Max allowed polynomial degree is 5");
00787 
00788     wradius = dfs_get_parameter_int(parlist, "fors.vimos_calib.wradius", NULL);
00789 
00790     if (wradius < 0)
00791         vimos_calib_exit("Invalid search radius");
00792 
00793     wreject = dfs_get_parameter_double(parlist, 
00794                                        "fors.vimos_calib.wreject", NULL);
00795 
00796     if (wreject <= 0.0)
00797         vimos_calib_exit("Invalid rejection threshold");
00798 
00799     wmode = dfs_get_parameter_int(parlist, "fors.vimos_calib.wmode", NULL);
00800 
00801     if (wmode < 0 || wmode > 2)
00802         vimos_calib_exit("Invalid wavelength solution interpolation mode");
00803 
00804     wcolumn = dfs_get_parameter_string(parlist, 
00805                                        "fors.vimos_calib.wcolumn", NULL);
00806 
00807     cdegree = dfs_get_parameter_int(parlist, 
00808                     "fors.vimos_calib.cdegree", grism_table);
00809 
00810     if (cdegree < 1)
00811         vimos_calib_exit("Invalid polynomial degree");
00812 
00813     if (cdegree > 5)
00814         vimos_calib_exit("Max allowed polynomial degree is 5");
00815 
00816     cmode = dfs_get_parameter_int(parlist, "fors.vimos_calib.cmode", NULL);
00817 
00818     if (cmode < 0 || cmode > 2)
00819         vimos_calib_exit("Invalid curvature solution interpolation mode");
00820 
00821     startwavelength = dfs_get_parameter_double(parlist, 
00822                     "fors.vimos_calib.startwavelength", grism_table);
00823     if (startwavelength > 1.0)
00824         if (startwavelength < 3000.0 || startwavelength > 13000.0)
00825             vimos_calib_exit("Invalid wavelength");
00826 
00827     endwavelength = dfs_get_parameter_double(parlist, 
00828                     "fors.vimos_calib.endwavelength", grism_table);
00829     if (endwavelength > 1.0) {
00830         if (endwavelength < 3000.0 || endwavelength > 13000.0)
00831             vimos_calib_exit("Invalid wavelength");
00832         if (startwavelength < 1.0)
00833             vimos_calib_exit("Invalid wavelength interval");
00834     }
00835 
00836     if (startwavelength > 1.0)
00837         if (endwavelength - startwavelength <= 0.0)
00838             vimos_calib_exit("Invalid wavelength interval");
00839 
00840     reference = dfs_get_parameter_double(parlist,
00841                 "fors.vimos_calib.reference", grism_table);
00842 
00843     if (reference < startwavelength || reference > endwavelength)
00844         vimos_calib_exit("Invalid reference wavelength");
00845 
00846     slit_ident = dfs_get_parameter_bool(parlist, 
00847                     "fors.vimos_calib.slit_ident", NULL);
00848 
00849     sdegree = dfs_get_parameter_int(parlist, "fors.vimos_calib.sdegree", NULL);
00850     ddegree = dfs_get_parameter_int(parlist, "fors.vimos_calib.ddegree", NULL);
00851     sradius = dfs_get_parameter_int(parlist, "fors.vimos_calib.sradius", NULL);
00852     dradius = dfs_get_parameter_int(parlist, "fors.vimos_calib.dradius", NULL);
00853 
00854     if (sradius < 1 || dradius < 1)
00855         vimos_calib_exit("Invalid smoothing box radius");
00856 
00857     qc = dfs_get_parameter_bool(parlist, "fors.vimos_calib.qc", NULL);
00858 
00859     check = dfs_get_parameter_bool(parlist, "fors.vimos_calib.check", NULL);
00860 
00861     cpl_table_delete(grism_table); grism_table = NULL;
00862 
00863     if (cpl_error_get_code())
00864         vimos_calib_exit("Failure getting the configuration parameters");
00865 
00866 
00867     /* 
00868      * Check input set-of-frames
00869      */
00870 
00871     cpl_msg_indent_less();
00872     cpl_msg_info(recipe, "Check input set-of-frames:");
00873     cpl_msg_indent_more();
00874 
00875     if (!dfs_equal_keyword(frameset, "ESO OCS CON QUAD")) 
00876         vimos_calib_exit("Input frames are not from the same quadrant");
00877 
00878     mos = cpl_frameset_count_tags(frameset, "MOS_ARC_SPECTRUM");
00879 
00880     if (mos == 0)
00881         vimos_calib_exit("Missing input arc lamp frame");
00882 
00883     if (mos > 1)
00884         vimos_calib_exit("Just one input arc lamp frame is allowed"); 
00885 
00886     arc_tag                  = "MOS_ARC_SPECTRUM";
00887     flat_tag                 = "MOS_SCREEN_FLAT";
00888     master_screen_flat_tag   = "MOS_COMBINED_SCREEN_FLAT";
00889     master_norm_flat_tag     = "MOS_MASTER_SCREEN_FLAT";
00890     reduced_lamp_tag         = "MOS_ARC_SPECTRUM_EXTRACTED";
00891     disp_residuals_tag       = "MOS_DISP_RESIDUALS";
00892     disp_coeff_tag           = "MOS_DISP_COEFF";
00893     wavelength_map_tag       = "MOS_WAVELENGTH_MAP";
00894     spectra_detection_tag    = "MOS_SPECTRA_DETECTION";
00895     spectral_resolution_tag  = "MOS_SPECTRAL_RESOLUTION";
00896     slit_map_tag             = "MOS_SLIT_MAP";
00897     curv_traces_tag          = "MOS_CURV_TRACES";
00898     curv_coeff_tag           = "MOS_CURV_COEFF";
00899     spatial_map_tag          = "MOS_SPATIAL_MAP";
00900     slit_location_tag        = "MOS_SLIT_LOCATION";
00901     disp_residuals_table_tag = "MOS_DISP_RESIDUALS_TABLE";
00902     delta_image_tag          = "MOS_DELTA_IMAGE";
00903 
00904     nbias = 0;
00905     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") == 0) {
00906         if (cpl_frameset_count_tags(frameset, "BIAS") == 0)
00907             vimos_calib_exit("Missing required input: MASTER_BIAS or BIAS");
00908         nbias = cpl_frameset_count_tags(frameset, "BIAS");
00909     }
00910 
00911     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") > 1)
00912         vimos_calib_exit("Too many in input: MASTER_BIAS");
00913 
00914     if (cpl_frameset_count_tags(frameset, "LINE_CATALOG") == 0)
00915         vimos_calib_exit("Missing required input: LINE_CATALOG");
00916 
00917     if (cpl_frameset_count_tags(frameset, "LINE_CATALOG") > 1)
00918         vimos_calib_exit("Too many in input: LINE_CATALOG");
00919 
00920     nflats = cpl_frameset_count_tags(frameset, flat_tag);
00921 
00922     if (nflats < 1) {
00923         cpl_msg_error(recipe, "Missing required input: %s", flat_tag);
00924         vimos_calib_exit(NULL);
00925     }
00926 
00927     cpl_msg_indent_less();
00928 
00929     if (nflats > 1)
00930         cpl_msg_info(recipe, "Load %d flat field frames and sum them...",
00931                      nflats);
00932     else
00933         cpl_msg_info(recipe, "Load flat field exposure...");
00934 
00935     cpl_msg_indent_more();
00936 
00937     header = dfs_load_header(frameset, flat_tag, 0);
00938 
00939     if (header == NULL)
00940         vimos_calib_exit("Cannot load flat field frame header");
00941 
00942     alltime = cpl_propertylist_get_double(header, "EXPTIME");
00943 
00944     if (cpl_error_get_code() != CPL_ERROR_NONE)
00945         vimos_calib_exit("Missing keyword EXPTIME in flat field frame header");
00946 
00947     cpl_propertylist_delete(header);
00948 
00949     for (i = 1; i < nflats; i++) {
00950 
00951         header = dfs_load_header(frameset, NULL, 0);
00952 
00953         if (header == NULL)
00954             vimos_calib_exit("Cannot load flat field frame header");
00955 
00956         alltime += cpl_propertylist_get_double(header, "EXPTIME");
00957 
00958         if (cpl_error_get_code() != CPL_ERROR_NONE)
00959             vimos_calib_exit("Missing keyword EXPTIME in flat field "
00960                             "frame header");
00961 
00962         cpl_propertylist_delete(header);
00963 
00964     }
00965 
00966     master_flat = dfs_load_image(frameset, flat_tag, CPL_TYPE_FLOAT, 0, 0);
00967 
00968     if (master_flat == NULL)
00969         vimos_calib_exit("Cannot load flat field");
00970 
00971     for (i = 1; i < nflats; i++) {
00972         flat = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
00973         if (flat) {
00974             cpl_image_add(master_flat, flat);
00975             cpl_image_delete(flat); flat = NULL;
00976         }
00977         else
00978             vimos_calib_exit("Cannot load flat field");
00979     }
00980 
00981 
00982     /*
00983      * Get some info from arc lamp header
00984      */
00985 
00986     header = dfs_load_header(frameset, arc_tag, 0);
00987 
00988     if (header == NULL)
00989         vimos_calib_exit("Cannot load arc lamp header");
00990 
00991     instrume = (char *)cpl_propertylist_get_string(header, "INSTRUME");
00992     if (instrume == NULL)
00993         vimos_calib_exit("Missing keyword INSTRUME in arc lamp header");
00994     instrume = cpl_strdup(instrume);
00995 
00996     quadrant = cpl_propertylist_get_int(header, "ESO OCS CON QUAD");
00997 
00998     switch (quadrant) {
00999     case 1:
01000         key_gris_name = "ESO INS GRIS1 NAME";
01001         key_gris_id = "ESO INS GRIS1 ID";
01002         key_filt_name = "ESO INS FILT1 NAME";
01003         key_filt_id = "ESO INS FILT1 ID";
01004         break;
01005     case 2:
01006         key_gris_name = "ESO INS GRIS2 NAME";
01007         key_gris_id = "ESO INS GRIS2 ID";
01008         key_filt_name = "ESO INS FILT2 NAME";
01009         key_filt_id = "ESO INS FILT2 ID";
01010         break;
01011     case 3:
01012         key_gris_name = "ESO INS GRIS3 NAME";
01013         key_gris_id = "ESO INS GRIS3 ID";
01014         key_filt_name = "ESO INS FILT3 NAME";
01015         key_filt_id = "ESO INS FILT3 ID";
01016         break;
01017     case 4:
01018         key_gris_name = "ESO INS GRIS4 NAME";
01019         key_gris_id = "ESO INS GRIS4 ID";
01020         key_filt_name = "ESO INS FILT4 NAME";
01021         key_filt_id = "ESO INS FILT4 ID";
01022         break;
01023     }
01024 
01025 /*
01026     if (!dfs_equal_keyword(frameset, key_gris_id))
01027         vimos_calib_exit("Input frames are not from the same grism");
01028 
01029     if (!dfs_equal_keyword(frameset, key_filt_id))
01030         vimos_calib_exit("Input frames are not from the same filter");
01031 */
01032 
01033     gain = cpl_propertylist_get_double(header, "ESO DET OUT1 CONAD");
01034 
01035     if (cpl_error_get_code() != CPL_ERROR_NONE)
01036         vimos_calib_exit("Missing keyword ESO DET OUT1 CONAD in arc lamp "
01037                         "frame header");
01038 
01039     cpl_msg_info(recipe, "The gain factor is: %.2f e-/ADU", gain);
01040 
01041     cpl_msg_info(recipe, "Produce mask slit position table...");
01042 
01043     maskslits = mos_load_slits_vimos(header);
01044 
01045     /*
01046      * Check if all slits have the same X offset: in such case, 
01047      * treat the observation as a long-slit one!
01048      */
01049 
01050     mxpos = cpl_table_get_column_median(maskslits, "xtop");
01051     xpos = cpl_table_get_data_double(maskslits, "xtop");
01052     nslits = cpl_table_get_nrow(maskslits);
01053 
01054     treat_as_lss = 1;
01055     for (i = 0; i < nslits; i++) {
01056         if (fabs(mxpos-xpos[i]) > 0.01) {
01057             treat_as_lss = 0;
01058             break;
01059         }
01060     }
01061 
01062     if (treat_as_lss) {
01063         cpl_msg_warning(recipe, "All MOS slits have the same offset: %.2f\n"
01064                         "The long-slit data reduction strategy is applied!", 
01065                         mxpos);
01066         cpl_table_delete(maskslits); maskslits = NULL;
01067     }
01068 
01069     if (slit_ident == 0) {
01070         cpl_table_delete(maskslits); maskslits = NULL;
01071     }
01072 
01073 
01074     /* Leave the header on for the next step... */
01075 
01076 
01077     /*
01078      * Remove the master bias
01079      */
01080 
01081     if (nbias) {
01082 
01083         /*
01084          * Set of raw BIASes in input, need to create master bias!
01085          */
01086 
01087         cpl_msg_info(recipe, "Generate the master from input raw biases...");
01088 
01089         if (nbias > 1) {
01090 
01091             biases = cpl_imagelist_new();
01092 
01093             bias = dfs_load_image(frameset, "BIAS", CPL_TYPE_FLOAT, 0, 0);
01094     
01095             if (bias == NULL)
01096                 vimos_calib_exit("Cannot load bias frame");
01097 
01098             cpl_imagelist_set(biases, bias, 0);
01099     
01100             for (i = 1; i < nbias; i++) {
01101                 bias = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
01102                 if (bias)
01103                     cpl_imagelist_set(biases, bias, i);
01104                 else
01105                     vimos_calib_exit("Cannot load bias frame");
01106             }
01107     
01108             master_bias = cpl_imagelist_collapse_median_create(biases);
01109 
01110             cpl_imagelist_delete(biases);
01111         }
01112         else
01113             master_bias = dfs_load_image(frameset, "BIAS", 
01114                                          CPL_TYPE_FLOAT, 0, 1);
01115 
01116     }
01117     else {
01118         master_bias = dfs_load_image(frameset, "MASTER_BIAS", 
01119                                      CPL_TYPE_FLOAT, 0, 1);
01120         if (master_bias == NULL)
01121             vimos_calib_exit("Cannot load master bias");
01122     }
01123 
01124     cpl_msg_info(recipe, "Remove the master bias...");
01125 
01126     overscans = mos_load_overscans_vimos(header, 1);
01127     cpl_propertylist_delete(header); header = NULL;
01128 
01129     if (nbias) {
01130         int xlow = cpl_table_get_int(overscans, "xlow", 0, NULL);
01131         int ylow = cpl_table_get_int(overscans, "ylow", 0, NULL);
01132         int xhig = cpl_table_get_int(overscans, "xhig", 0, NULL);
01133         int yhig = cpl_table_get_int(overscans, "yhig", 0, NULL);
01134         dummy = cpl_image_extract(master_bias, xlow+1, ylow+1, xhig, yhig);
01135         cpl_image_delete(master_bias); master_bias = dummy;
01136 
01137         if (dfs_save_image(frameset, master_bias, "MASTER_BIAS",
01138                            NULL, parlist, recipe, version))
01139             vimos_calib_exit(NULL);
01140     }
01141 
01142     if (nflats > 1) {
01143         multi_bias = cpl_image_multiply_scalar_create(master_bias, nflats);
01144         dummy = mos_remove_bias(master_flat, multi_bias, overscans);
01145         cpl_image_delete(multi_bias);
01146     }
01147     else {
01148         dummy = mos_remove_bias(master_flat, master_bias, overscans);
01149     }
01150     cpl_image_delete(master_flat);
01151     master_flat = dummy;
01152 
01153     if (master_flat == NULL)
01154         vimos_calib_exit("Cannot remove bias from flat field");
01155 
01156     cpl_msg_indent_less();
01157     cpl_msg_info(recipe, "Load arc lamp exposure...");
01158     cpl_msg_indent_more();
01159 
01160     spectra = dfs_load_image(frameset, arc_tag, CPL_TYPE_FLOAT, 0, 0);
01161 
01162     if (spectra == NULL)
01163         vimos_calib_exit("Cannot load arc lamp exposure");
01164 
01165     cpl_msg_info(recipe, "Remove the master bias...");
01166 
01167     dummy = mos_remove_bias(spectra, master_bias, overscans);
01168     cpl_table_delete(overscans); overscans = NULL;
01169     cpl_image_delete(master_bias); master_bias = NULL;
01170     cpl_image_delete(spectra); spectra = dummy;
01171 
01172     if (spectra == NULL)
01173         vimos_calib_exit("Cannot remove bias from arc lamp exposure");
01174 
01175     cpl_msg_indent_less();
01176     cpl_msg_info(recipe, "Load input line catalog...");
01177     cpl_msg_indent_more();
01178 
01179     wavelengths = dfs_load_table(frameset, "LINE_CATALOG", 1);
01180 
01181     if (wavelengths == NULL)
01182         vimos_calib_exit("Cannot load line catalog");
01183 
01184 
01185     /*
01186      * Cast the wavelengths into a (double precision) CPL vector
01187      */
01188 
01189     nlines = cpl_table_get_nrow(wavelengths);
01190 
01191     if (nlines == 0)
01192         vimos_calib_exit("Empty input line catalog");
01193 
01194     if (cpl_table_has_column(wavelengths, wcolumn) != 1) {
01195         cpl_msg_error(recipe, "Missing column %s in input line catalog table",
01196                       wcolumn);
01197         vimos_calib_exit(NULL);
01198     }
01199 
01200     line = cpl_malloc(nlines * sizeof(double));
01201     
01202     for (i = 0; i < nlines; i++)
01203         line[i] = cpl_table_get(wavelengths, wcolumn, i, NULL);
01204 
01205     cpl_table_delete(wavelengths); wavelengths = NULL;
01206 
01207     lines = cpl_vector_wrap(nlines, line);
01208 
01209 
01210     /*
01211      * Rotate frames horizontally with red to the right
01212      */
01213 
01214     cpl_image_turn(spectra, rotate);
01215     cpl_image_turn(master_flat, rotate);
01216 
01217     if (treat_as_lss) {
01218 
01219         /* FIXME:
01220          * The LSS data calibration is still dirty: it doesn't apply
01221          * any spatial rectification, and only in future an external
01222          * spectral curvature model would be provided in input. Here
01223          * and there temporary solutions are adpted, such as accepting
01224          * the preliminary wavelength calibration.
01225          */
01226 
01227         /*
01228          * Flat field normalisation is done directly on the master flat
01229          * field (without spatial rectification first). The spectral
01230          * curvature model may be provided in input, in future releases.
01231          */
01232 
01233         cpl_msg_indent_less();
01234         cpl_msg_info(recipe, "Perform flat field normalisation...");
01235         cpl_msg_indent_more();
01236 
01237         norm_flat = cpl_image_duplicate(master_flat);
01238 
01239         smo_flat = mos_normalise_longflat(norm_flat, sradius, dradius,
01240                                           sdegree);
01241 
01242         cpl_image_delete(smo_flat); smo_flat = NULL; /* It may be a product */
01243 
01244         save_header = dfs_load_header(frameset, flat_tag, 0);
01245         cpl_propertylist_update_int(save_header, "ESO PRO DATANCOM", nflats);
01246 
01247         cpl_image_turn(master_flat, rotate_back);
01248         if (dfs_save_image(frameset, master_flat, master_screen_flat_tag,
01249                            save_header, parlist, recipe, version))
01250             vimos_calib_exit(NULL);
01251         cpl_image_delete(master_flat); master_flat = NULL;
01252 
01253         cpl_image_turn(norm_flat, rotate_back);
01254         if (dfs_save_image(frameset, norm_flat, master_norm_flat_tag,
01255                            save_header, parlist, recipe, version))
01256             vimos_calib_exit(NULL);
01257 
01258         cpl_propertylist_delete(save_header); save_header = NULL;
01259         cpl_image_delete(norm_flat); norm_flat = NULL;
01260 
01261 
01262         /*
01263          * In the case of LSS data, extract the spectra directly
01264          * on the first attempt. The spectral curvature model may
01265          * be provided in input, in future releases.
01266          */
01267 
01268         cpl_msg_indent_less();
01269         cpl_msg_info(recipe, "Perform wavelength calibration...");
01270         cpl_msg_indent_more();
01271 
01272         nx = cpl_image_get_size_x(spectra);
01273         ny = cpl_image_get_size_y(spectra);
01274 
01275         wavemap = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01276         if (check)
01277             residual = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01278 
01279         fiterror = cpl_calloc(ny, sizeof(double));
01280         fitlines = cpl_calloc(ny, sizeof(int));
01281         idscoeff = cpl_table_new(ny);
01282         restable = cpl_table_new(nlines);
01283 
01284     if (mos_saturation_process(spectra))
01285         vimos_calib_exit("Cannot process saturation");
01286 
01287     if (mos_subtract_background(spectra))
01288         vimos_calib_exit("Cannot subtract the background");
01289 
01290         rectified = mos_wavelength_calibration_raw(spectra, lines, dispersion,
01291                                                    peakdetection, wradius,
01292                                                    wdegree, wreject, reference,
01293                                                    &startwavelength, 
01294                                                    &endwavelength, fitlines, 
01295                                                    fiterror, idscoeff, wavemap,
01296                                                    residual, restable, NULL);
01297 
01298         if (rectified == NULL)
01299             vimos_calib_exit("Wavelength calibration failure.");
01300 
01301         if (dfs_save_table(frameset, restable, disp_residuals_table_tag, NULL,
01302                            parlist, recipe, version))
01303             vimos_calib_exit(NULL);
01304 
01305         cpl_table_delete(restable); restable = NULL;
01306 
01307         if (wmode) {
01308             cpl_image_delete(rectified); rectified = NULL;
01309             cpl_image_delete(wavemap); wavemap = NULL;
01310             mos_interpolate_wavecalib(idscoeff, wavemap, wmode);
01311             wavemap = mos_map_idscoeff(idscoeff, nx, reference,
01312                                        startwavelength, endwavelength);
01313             rectified = mos_wavelength_calibration(spectra, reference,
01314                                                    startwavelength, 
01315                                                    endwavelength, dispersion, 
01316                                                    idscoeff, 0);
01317         }
01318 
01319         cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
01320         cpl_table_set_column_unit(idscoeff, "error", "pixel");
01321         cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;
01322 
01323         for (i = 0; i < ny; i++)
01324             if (!cpl_table_is_valid(idscoeff, "c0", i))
01325                 cpl_table_set_invalid(idscoeff, "error", i);
01326 
01327         delta = mos_map_pixel(idscoeff, reference, startwavelength,
01328                               endwavelength, dispersion, 2);
01329 
01330         header = cpl_propertylist_new();
01331         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01332         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01333         cpl_propertylist_update_double(header, "CRVAL1",
01334                                        startwavelength + dispersion/2);
01335         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01336         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01337         cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01338         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01339         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01340         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01341         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01342         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01343         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01344 
01345         if (dfs_save_image(frameset, delta, delta_image_tag,
01346                            header, parlist, recipe, version))
01347             vimos_calib_exit(NULL);
01348 
01349         cpl_image_delete(delta); delta = NULL;
01350         cpl_propertylist_delete(header); header = NULL;
01351 
01352         cpl_msg_info(recipe, "Valid solutions found: %d out of %d rows", 
01353                      ny - cpl_table_count_invalid(idscoeff, "c0"), ny);
01354 
01355         cpl_image_delete(spectra); spectra = NULL;
01356 
01357         mean_rms = mos_distortions_rms(rectified, lines, startwavelength,
01358                                        dispersion, 6, 0);
01359 
01360         cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);
01361 
01362         mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01363 
01364         cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)",
01365                      mean_rms, mean_rms * dispersion);
01366 
01367         restab = mos_resolution_table(rectified, startwavelength, dispersion,
01368                                       60000, lines);
01369 
01370         if (restab) {
01371             cpl_msg_info(recipe, "Mean spectral resolution: %.2f",
01372                   cpl_table_get_column_mean(restab, "resolution"));
01373             cpl_msg_info(recipe, 
01374                   "Mean reference lines FWHM: %.2f +/- %.2f pixel",
01375                   cpl_table_get_column_mean(restab, "fwhm") / dispersion,
01376                   cpl_table_get_column_mean(restab, "fwhm_rms") / dispersion);
01377 
01378             if (qc) {
01379 
01380                 header = dfs_load_header(frameset, arc_tag, 0);
01381 
01382                 if (header == NULL)
01383                     vimos_calib_exit("Cannot reload arc lamp header");
01384 
01385                 qclist = cpl_propertylist_new();
01386     
01387                 fors_qc_start_group(qclist, "1.1", instrume);
01388 
01389 
01390                 /*
01391                  * QC1 group header
01392                  */
01393 
01394                 if (fors_qc_write_string("PRO.CATG", spectral_resolution_tag,
01395                                         "Product category", instrume))
01396                     vimos_calib_exit("Cannot write product category to "
01397                                     "QC log file");
01398 
01399                 if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL,
01400                                           "DPR type", instrume))
01401                     vimos_calib_exit("Missing keyword DPR TYPE in arc "
01402                                     "lamp header");
01403 
01404                 if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL,
01405                                           "Template", instrume))
01406                     vimos_calib_exit("Missing keyword TPL ID in arc "
01407                                     "lamp header");
01408 
01409                 if (fors_qc_keyword_to_paf(header, key_gris_name, NULL,
01410                                           "Grism name", instrume))
01411                     vimos_calib_exit("Missing keyword INS GRISx NAME in arc "
01412                                     "lamp header");
01413 
01414 
01415                 if (fors_qc_keyword_to_paf(header, key_gris_id, NULL,
01416                                           "Grism identifier", instrume)) {
01417                     cpl_msg_error(recipe, "Missing keyword %s in arc "
01418                                     "lamp header", key_gris_id);
01419                     vimos_calib_exit(NULL);
01420                 }
01421 
01422                 if (cpl_propertylist_has(header, key_filt_name))
01423                     fors_qc_keyword_to_paf(header, key_filt_name, NULL,
01424                                           "Filter name", instrume);
01425 
01426                 if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
01427                                           "Collimator name", instrume))
01428                     vimos_calib_exit("Missing keyword INS COLL NAME in arc "
01429                                     "lamp header");
01430 
01431                 if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
01432                                           "Chip identifier", instrume))
01433                     vimos_calib_exit("Missing keyword DET CHIP1 ID in arc "
01434                                     "lamp header");
01435 
01436                 if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
01437                                           "Archive name of input data",
01438                                           instrume))
01439                     vimos_calib_exit("Missing keyword ARCFILE in arc "
01440                                     "lamp header");
01441 
01442                 cpl_propertylist_delete(header); header = NULL;
01443 
01444                 pipefile = dfs_generate_filename_tfits(spectral_resolution_tag);
01445                 if (fors_qc_write_string("PIPEFILE", pipefile,
01446                                         "Pipeline product name", instrume))
01447                     vimos_calib_exit("Cannot write PIPEFILE to QC log file");
01448                 cpl_free(pipefile); pipefile = NULL;
01449 
01450 
01451                 /*
01452                  * QC1 parameters
01453                  */
01454     
01455                 keyname = "QC.MOS.RESOLUTION";
01456 
01457                 if (fors_qc_write_qc_double(qclist, 
01458                                            cpl_table_get_column_mean(restab,
01459                                                                  "resolution"),
01460                                            keyname, NULL, 
01461                                            "Mean spectral resolution",
01462                                            instrume)) {
01463                     vimos_calib_exit("Cannot write mean spectral resolution to "
01464                                     "QC log file");
01465                 }
01466 
01467                 keyname = "QC.MOS.RESOLUTION.RMS"; 
01468 
01469                 if (fors_qc_write_qc_double(qclist,
01470                                            cpl_table_get_column_stdev(restab,
01471                                                                 "resolution"),
01472                                            keyname, NULL, 
01473                                            "Scatter of spectral resolution",
01474                                            instrume)) {
01475                     vimos_calib_exit("Cannot write spectral resolution scatter "
01476                                     "to QC log file");
01477                 }
01478 
01479                 keyname = "QC.MOS.RESOLUTION.NLINES";
01480 
01481                 if (fors_qc_write_qc_int(qclist, cpl_table_get_nrow(restab) -
01482                                         cpl_table_count_invalid(restab,
01483                                                                 "resolution"),
01484                                         keyname, NULL,
01485                                         "Number of lines for spectral "
01486                                         "resolution computation",
01487                                         instrume)) {
01488                     vimos_calib_exit("Cannot write number of lines used in "
01489                                     "spectral resolution computation "
01490                                     "to QC log file");
01491                 }
01492 
01493                 fors_qc_end_group();
01494 
01495             }  /* End of QC1 computation */
01496 
01497             if (dfs_save_table(frameset, restab, spectral_resolution_tag, 
01498                                qclist, parlist, recipe, version))
01499                 vimos_calib_exit(NULL);
01500 
01501             cpl_table_delete(restab); restab = NULL;
01502             cpl_propertylist_delete(qclist); qclist = NULL;
01503 
01504         }
01505         else
01506             vimos_calib_exit("Cannot compute the spectral resolution table");
01507 
01508         cpl_vector_delete(lines); lines = NULL;
01509 
01510 
01511         /*
01512          * Save rectified arc lamp spectrum to disk
01513          */
01514 
01515         header = cpl_propertylist_new();
01516         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01517         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01518         cpl_propertylist_update_double(header, "CRVAL1", 
01519                                        startwavelength + dispersion/2);
01520         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01521         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01522         cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01523         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01524         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01525         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01526         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01527         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01528         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01529         cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);
01530 
01531         if (dfs_save_image(frameset, rectified, reduced_lamp_tag, header, 
01532                            parlist, recipe, version))
01533             vimos_calib_exit(NULL);
01534 
01535         cpl_image_delete(rectified); rectified = NULL;
01536         cpl_propertylist_delete(header); header = NULL;
01537 
01538         if (dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL, 
01539                            parlist, recipe, version))
01540             vimos_calib_exit(NULL);
01541 
01542         cpl_table_delete(idscoeff); idscoeff = NULL;
01543 
01544         header = dfs_load_header(frameset, arc_tag, 0);
01545 
01546         if (header == NULL)
01547             vimos_calib_exit("Cannot reload arc lamp header");
01548 
01549         if (qc) {
01550 
01551             compute_qc = 0;
01552             if (fabs(mxpos) < 0.05)
01553                 compute_qc = 1;
01554 
01555             if (compute_qc) {
01556 
01557                 fors_qc_start_group(header, "1.1", instrume);
01558 
01559                 /*
01560                  * QC1 group header
01561                  */
01562 
01563                 if (fors_qc_write_string("PRO.CATG", wavelength_map_tag,
01564                                         "Product category", instrume))
01565                     vimos_calib_exit("Cannot write product category to "
01566                                     "QC log file");
01567 
01568                 if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL, 
01569                                           "DPR type", instrume))
01570                     vimos_calib_exit("Missing keyword DPR TYPE in arc "
01571                                     "lamp header");
01572 
01573                 if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL, 
01574                                           "Template", instrume))
01575                     vimos_calib_exit("Missing keyword TPL ID in arc "
01576                                     "lamp header");
01577 
01578                 if (fors_qc_keyword_to_paf(header, key_gris_name, NULL, 
01579                                           "Grism name", instrume)) {
01580                     cpl_msg_error(recipe, "Missing keyword %s in arc "
01581                                     "lamp header", key_gris_name);
01582                     vimos_calib_exit(NULL);
01583                 }
01584 
01585                 if (fors_qc_keyword_to_paf(header, key_gris_id, NULL,
01586                                           "Grism identifier", instrume)) {
01587                     cpl_msg_error(recipe, "Missing keyword %s in arc "
01588                                     "lamp header", key_gris_id);
01589                     vimos_calib_exit(NULL);
01590                 }
01591 
01592                 if (cpl_propertylist_has(header, key_filt_name))
01593                     fors_qc_keyword_to_paf(header, key_filt_name, NULL,
01594                                           "Filter name", instrume);
01595 
01596                 if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
01597                                           "Collimator name", instrume))
01598                     vimos_calib_exit("Missing keyword INS COLL NAME in arc "
01599                                     "lamp header");
01600 
01601                 if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
01602                                           "Chip identifier", instrume))
01603                     vimos_calib_exit("Missing keyword DET CHIP1 ID in arc "
01604                                     "lamp header");
01605 
01606                 if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
01607                                           "Archive name of input data", 
01608                                           instrume))
01609                     vimos_calib_exit("Missing keyword ARCFILE in arc "
01610                                     "lamp header");
01611     
01612                 pipefile = dfs_generate_filename(wavelength_map_tag);
01613                 if (fors_qc_write_string("PIPEFILE", pipefile, 
01614                                         "Pipeline product name", instrume))
01615                     vimos_calib_exit("Cannot write PIPEFILE to QC log file");
01616                 cpl_free(pipefile); pipefile = NULL;
01617 
01618 
01619                 /*
01620                  * QC1 parameters
01621                  */
01622 
01623                 data = cpl_image_get_data(wavemap);
01624     
01625                 if (fors_qc_write_qc_double(header, data[nx/2 + ny*nx/2],
01626                                            "QC.MOS.CENTRAL.WAVELENGTH",
01627                                            "Angstrom", 
01628                                            "Wavelength at CCD center",
01629                                            instrume)) {
01630                     vimos_calib_exit("Cannot write central wavelength to QC "
01631                                     "log file");
01632                 }
01633 
01634                 fors_qc_end_group();
01635 
01636             }
01637 
01638         }
01639 
01640         cpl_image_turn(wavemap, rotate_back);
01641         if (dfs_save_image(frameset, wavemap, wavelength_map_tag, header,
01642                            parlist, recipe, version))
01643             vimos_calib_exit(NULL);
01644 
01645         cpl_image_delete(wavemap); wavemap = NULL;
01646 
01647         if (check) {
01648             cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01649             cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01650             /* cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01651             cpl_propertylist_update_double(header, "CD1_1", 1.0);
01652             cpl_propertylist_update_double(header, "CD1_2", 0.0);
01653             cpl_propertylist_update_double(header, "CD2_1", 0.0);
01654             cpl_propertylist_update_double(header, "CD2_2", 1.0);
01655             cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01656             cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01657 
01658             if (dfs_save_image(frameset, residual, disp_residuals_tag, header,
01659                                parlist, recipe, version))
01660                 vimos_calib_exit(NULL);
01661 
01662             cpl_image_delete(residual); residual = NULL;
01663         }
01664 
01665         cpl_propertylist_delete(header); header = NULL;
01666 
01667         return 0;         /* Successful LSS data reduction */
01668 
01669     }   /* End of LSS data reduction section */
01670 
01671 
01672     /*
01673      * Here the MOS calibration is carried out.
01674      */
01675 
01676     /*
01677      * Detecting spectra on the CCD
01678      */
01679 
01680     cpl_msg_indent_less();
01681     cpl_msg_info(recipe, "Detecting spectra on CCD...");
01682     cpl_msg_indent_more();
01683 
01684     ccd_xsize = nx = cpl_image_get_size_x(spectra);
01685     ccd_ysize = ny = cpl_image_get_size_y(spectra);
01686 
01687     refmask = cpl_mask_new(nx, ny);
01688 
01689     if (mos_saturation_process(spectra))
01690     vimos_calib_exit("Cannot process saturation");
01691 
01692     if (mos_subtract_background(spectra))
01693     vimos_calib_exit("Cannot subtract the background");
01694 
01695     checkwave = mos_wavelength_calibration_raw(spectra, lines, dispersion, 
01696                                                peakdetection, wradius, 
01697                                                wdegree, wreject, reference,
01698                                                &startwavelength, &endwavelength,
01699                                                NULL, NULL, NULL, NULL, NULL, 
01700                                                NULL, refmask);
01701 
01702     if (checkwave == NULL)
01703         vimos_calib_exit("Wavelength calibration failure.");
01704 
01705     /*
01706      * Save check image to disk
01707      */
01708 
01709     header = cpl_propertylist_new();
01710     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01711     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01712     cpl_propertylist_update_double(header, "CRVAL1", 
01713                                    startwavelength + dispersion/2);
01714     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01715     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01716     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01717     cpl_propertylist_update_double(header, "CD1_1", dispersion);
01718     cpl_propertylist_update_double(header, "CD1_2", 0.0);
01719     cpl_propertylist_update_double(header, "CD2_1", 0.0);
01720     cpl_propertylist_update_double(header, "CD2_2", 1.0);
01721     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01722     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01723 
01724     if (check) {
01725         if (dfs_save_image(frameset, checkwave, spectra_detection_tag, header, 
01726                            parlist, recipe, version))
01727             vimos_calib_exit(NULL);
01728     }
01729 
01730     cpl_image_delete(checkwave); checkwave = NULL;
01731     cpl_propertylist_delete(header); header = NULL;
01732 
01733     cpl_msg_info(recipe, "Locate slits at reference wavelength on CCD...");
01734     slits = mos_locate_spectra(refmask);
01735 
01736     if (!slits) {
01737         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
01738         vimos_calib_exit("No slits could be detected!");
01739     }
01740 
01741     refimage = cpl_image_new_from_mask(refmask);
01742     cpl_mask_delete(refmask); refmask = NULL;
01743 
01744     if (check) {
01745         save_header = dfs_load_header(frameset, arc_tag, 0);
01746         cpl_image_turn(refimage, rotate_back);
01747         if (dfs_save_image(frameset, refimage, slit_map_tag, NULL,
01748                            parlist, recipe, version))
01749             vimos_calib_exit(NULL);
01750         cpl_propertylist_delete(save_header); save_header = NULL;
01751     }
01752 
01753     cpl_image_delete(refimage); refimage = NULL;
01754 
01755     if (slit_ident) {
01756 
01757         /*
01758          * Attempt slit identification: this recipe may continue even
01759          * in case of failed identification (i.e., the position table is 
01760          * not produced, but an error is not set). In case of failure,
01761          * the spectra would be still extracted, even if they would not
01762          * be associated to slits on the mask.
01763          * 
01764          * The reason for making the slit identification an user option 
01765          * (via the parameter slit_ident) is to offer the possibility 
01766          * to avoid identifications that are only apparently successful, 
01767          * as it would happen in the case of an incorrect slit description 
01768          * in the data header.
01769          */
01770 
01771         cpl_msg_indent_less();
01772         cpl_msg_info(recipe, "Attempt slit identification (optional)...");
01773         cpl_msg_indent_more();
01774 
01775         mos_rotate_slits(maskslits, -rotate, 0, 0);
01776         positions = mos_identify_slits(slits, maskslits, NULL);
01777 
01778         if (positions) {
01779             cpl_table_delete(slits);
01780             slits = positions;
01781 
01782             /*
01783              * Eliminate slits which are _entirely_ outside the CCD
01784              */
01785 
01786             cpl_table_and_selected_double(slits, 
01787                                           "ybottom", CPL_GREATER_THAN, ny-1);
01788             cpl_table_or_selected_double(slits, 
01789                                           "ytop", CPL_LESS_THAN, 0);
01790             cpl_table_erase_selected(slits);
01791 
01792             nslits = cpl_table_get_nrow(slits);
01793 
01794             if (nslits == 0)
01795                 vimos_calib_exit("No slits found on the CCD");
01796 
01797             cpl_msg_info(recipe, "%d slits are entirely or partially "
01798                          "contained in CCD", nslits);
01799 
01800         }
01801         else {
01802             slit_ident = 0;
01803             cpl_msg_info(recipe, "Global distortion model cannot be computed");
01804             if (cpl_error_get_code() != CPL_ERROR_NONE) {
01805                 vimos_calib_exit(NULL);
01806             }
01807         }
01808     }
01809 
01810 
01811     /*
01812      * Determination of spectral curvature
01813      */
01814 
01815     cpl_msg_indent_less();
01816     cpl_msg_info(recipe, "Determining spectral curvature...");
01817     cpl_msg_indent_more();
01818 
01819     cpl_msg_info(recipe, "Tracing master flat field spectra edges...");
01820     traces = mos_trace_flat(master_flat, slits, reference, 
01821                             startwavelength, endwavelength, dispersion);
01822 
01823     if (!traces)
01824         vimos_calib_exit("Tracing failure");
01825 
01826     cpl_msg_info(recipe, "Fitting flat field spectra edges...");
01827     polytraces = mos_poly_trace(slits, traces, cdegree);
01828 
01829     if (!polytraces)
01830         vimos_calib_exit("Trace fitting failure");
01831 
01832     if (cmode) {
01833         cpl_msg_info(recipe, "Computing global spectral curvature model...");
01834         mos_global_trace(slits, polytraces, cmode);
01835     }
01836 
01837     if (dfs_save_table(frameset, traces, curv_traces_tag, NULL, parlist,
01838                        recipe, version))
01839         vimos_calib_exit(NULL);
01840 
01841     cpl_table_delete(traces); traces = NULL;
01842 
01843     coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01844     spatial = mos_spatial_calibration(spectra, slits, polytraces, reference, 
01845                                       startwavelength, endwavelength, 
01846                                       dispersion, 0, coordinate);
01847 
01848     if (!slit_ident) {
01849         cpl_image_delete(spectra); spectra = NULL;
01850     }
01851 
01852     /*
01853      * Flat field normalisation is done directly on the master flat
01854      * field (without spatial rectification first). The spectral
01855      * curvature model may be provided in input, in future releases.
01856      */
01857 
01858     cpl_msg_indent_less();
01859     cpl_msg_info(recipe, "Perform flat field normalisation...");
01860     cpl_msg_indent_more();
01861 
01862     norm_flat = cpl_image_duplicate(master_flat);
01863 
01864     smo_flat = mos_normalise_flat(norm_flat, coordinate, slits, polytraces, 
01865                                   reference, startwavelength, endwavelength,
01866                                   dispersion, dradius, ddegree);
01867 
01868     cpl_image_delete(smo_flat); smo_flat = NULL;  /* It may be a product */
01869  
01870     save_header = dfs_load_header(frameset, flat_tag, 0);
01871     cpl_propertylist_update_int(save_header, "ESO PRO DATANCOM", nflats);
01872 
01873     cpl_image_turn(master_flat, rotate_back);
01874     if (dfs_save_image(frameset, master_flat, master_screen_flat_tag,
01875                        save_header, parlist, recipe, version))
01876         vimos_calib_exit(NULL);
01877 
01878     cpl_image_delete(master_flat); master_flat = NULL;
01879 
01880     cpl_image_turn(norm_flat, rotate_back);
01881     if (dfs_save_image(frameset, norm_flat, master_norm_flat_tag,
01882                        save_header, parlist, recipe, version))
01883         vimos_calib_exit(NULL);
01884 
01885     cpl_image_delete(norm_flat); norm_flat = NULL;
01886     cpl_propertylist_delete(save_header); save_header = NULL;
01887 
01888 
01889     /*
01890      * Final wavelength calibration of spectra having their curvature
01891      * removed
01892      */
01893 
01894     cpl_msg_indent_less();
01895     cpl_msg_info(recipe, "Perform final wavelength calibration...");
01896     cpl_msg_indent_more();
01897 
01898     nx = cpl_image_get_size_x(spatial);
01899     ny = cpl_image_get_size_y(spatial);
01900 
01901     idscoeff = cpl_table_new(ny);
01902     restable = cpl_table_new(nlines);
01903     rainbow = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01904     if (check)
01905         residual = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01906     fiterror = cpl_calloc(ny, sizeof(double));
01907     fitlines = cpl_calloc(ny, sizeof(int));
01908 
01909     rectified = mos_wavelength_calibration_final(spatial, slits, lines, 
01910                                                  dispersion, peakdetection, 
01911                                                  wradius, wdegree, wreject,
01912                                                  reference, &startwavelength, 
01913                                                  &endwavelength, fitlines, 
01914                                                  fiterror, idscoeff, rainbow, 
01915                                                  residual, restable);
01916 
01917     cpl_image_delete(spatial); spatial = NULL;
01918 
01919     if (rectified == NULL)
01920         vimos_calib_exit("Wavelength calibration failure.");
01921 
01922     if (dfs_save_table(frameset, restable, disp_residuals_table_tag, NULL,
01923                        parlist, recipe, version))
01924         vimos_calib_exit(NULL);
01925 
01926     cpl_table_delete(restable); restable = NULL;
01927 
01928     cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
01929     cpl_table_set_column_unit(idscoeff, "error", "pixel");
01930     cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;
01931 
01932     for (i = 0; i < ny; i++)
01933         if (!cpl_table_is_valid(idscoeff, "c0", i))
01934             cpl_table_set_invalid(idscoeff, "error", i);
01935 
01936     delta = mos_map_pixel(idscoeff, reference, startwavelength,
01937                           endwavelength, dispersion, 2);
01938 
01939     header = cpl_propertylist_new();
01940     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01941     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01942     cpl_propertylist_update_double(header, "CRVAL1",
01943                                    startwavelength + dispersion/2);
01944     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01945     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01946     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01947     cpl_propertylist_update_double(header, "CD1_1", dispersion);
01948     cpl_propertylist_update_double(header, "CD1_2", 0.0);
01949     cpl_propertylist_update_double(header, "CD2_1", 0.0);
01950     cpl_propertylist_update_double(header, "CD2_2", 1.0);
01951     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01952     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01953 
01954     if (dfs_save_image(frameset, delta, delta_image_tag,
01955                        header, parlist, recipe, version))
01956         vimos_calib_exit(NULL);
01957 
01958     cpl_image_delete(delta); delta = NULL;
01959     cpl_propertylist_delete(header); header = NULL;
01960 
01961     mean_rms = mos_distortions_rms(rectified, lines, startwavelength, 
01962                                    dispersion, 6, 0);
01963 
01964     cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);
01965 
01966     mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01967 
01968     cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)", 
01969                  mean_rms, mean_rms * dispersion);
01970 
01971     restab = mos_resolution_table(rectified, startwavelength, dispersion, 
01972                                   60000, lines);
01973 
01974     if (restab) {
01975         cpl_msg_info(recipe, "Mean spectral resolution: %.2f", 
01976                    cpl_table_get_column_mean(restab, "resolution"));
01977         cpl_msg_info(recipe, "Mean reference lines FWHM: %.2f +/- %.2f pixel",
01978                    cpl_table_get_column_mean(restab, "fwhm") / dispersion,
01979                    cpl_table_get_column_mean(restab, "fwhm_rms") / dispersion);
01980 
01981         if (qc) {
01982 
01983             header = dfs_load_header(frameset, arc_tag, 0);
01984 
01985             if (header == NULL)
01986                 vimos_calib_exit("Cannot reload arc lamp header");
01987 
01988             qclist = cpl_propertylist_new();
01989 
01990             fors_qc_start_group(qclist, "1.1", instrume);
01991 
01992 
01993             /*
01994              * QC1 group header
01995              */
01996 
01997             if (fors_qc_write_string("PRO.CATG", spectral_resolution_tag,
01998                                     "Product category", instrume))
01999                 vimos_calib_exit("Cannot write product category to "
02000                                 "QC log file");
02001 
02002             if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL,
02003                                       "DPR type", instrume))
02004                 vimos_calib_exit("Missing keyword DPR TYPE in arc "
02005                                 "lamp header");
02006 
02007             if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL,
02008                                       "Template", instrume))
02009                 vimos_calib_exit("Missing keyword TPL ID in arc "
02010                                 "lamp header");
02011 
02012             if (fors_qc_keyword_to_paf(header, key_gris_name, NULL,
02013                                       "Grism name", instrume)) {
02014                 cpl_msg_error(recipe, "Missing keyword %s in arc "
02015                                 "lamp header", key_gris_name);
02016                 vimos_calib_exit(NULL);
02017             }
02018 
02019             if (fors_qc_keyword_to_paf(header, key_gris_id, NULL,
02020                                       "Grism identifier", instrume)) {
02021                 cpl_msg_error(recipe, "Missing keyword %s in arc "
02022                                 "lamp header", key_gris_id);
02023                 vimos_calib_exit(NULL);
02024             }
02025 
02026             if (cpl_propertylist_has(header, key_filt_name))
02027                 fors_qc_keyword_to_paf(header, key_filt_name, NULL,
02028                                       "Filter name", instrume);
02029 
02030             if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
02031                                       "Collimator name", instrume))
02032                 vimos_calib_exit("Missing keyword INS COLL NAME in arc "
02033                                 "lamp header");
02034 
02035             if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
02036                                       "Chip identifier", instrume))
02037                 vimos_calib_exit("Missing keyword DET CHIP1 ID in arc "
02038                                 "lamp header");
02039 
02040             if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
02041                                       "Archive name of input data", 
02042                                       instrume))
02043                 vimos_calib_exit("Missing keyword ARCFILE in arc "
02044                                 "lamp header");
02045 
02046             cpl_propertylist_delete(header); header = NULL;
02047 
02048             pipefile = dfs_generate_filename_tfits(spectral_resolution_tag);
02049             if (fors_qc_write_string("PIPEFILE", pipefile,
02050                                     "Pipeline product name", instrume))
02051                 vimos_calib_exit("Cannot write PIPEFILE to QC log file");
02052             cpl_free(pipefile); pipefile = NULL;
02053 
02054 
02055             /*
02056              * QC1 parameters
02057              */
02058 
02059             keyname = "QC.MOS.RESOLUTION";
02060 
02061             if (fors_qc_write_qc_double(qclist, 
02062                                        cpl_table_get_column_mean(restab,
02063                                                                  "resolution"),
02064                                        keyname,
02065                                        "Angstrom",
02066                                        "Mean spectral resolution",
02067                                        instrume)) {
02068                 vimos_calib_exit("Cannot write mean spectral resolution to QC "
02069                                 "log file");
02070             }
02071 
02072             keyname = "QC.MOS.RESOLUTION.RMS";
02073 
02074             if (fors_qc_write_qc_double(qclist, 
02075                                        cpl_table_get_column_stdev(restab, 
02076                                                                   "resolution"),
02077                                        "QC.MOS.RESOLUTION.RMS",
02078                                        "Angstrom", 
02079                                        "Scatter of spectral resolution",
02080                                        instrume)) {
02081                 vimos_calib_exit("Cannot write spectral resolution scatter "
02082                                 "to QC log file");
02083             }
02084 
02085             keyname = "QC.MOS.RESOLUTION.NLINES";
02086 
02087             if (fors_qc_write_qc_int(qclist, cpl_table_get_nrow(restab) -
02088                                     cpl_table_count_invalid(restab, 
02089                                                             "resolution"),
02090                                     "QC.MOS.RESOLUTION.NLINES",
02091                                     NULL,
02092                                     "Number of lines for spectral resolution "
02093                                     "computation",
02094                                     instrume)) {
02095                 vimos_calib_exit("Cannot write number of lines used in "
02096                                 "spectral resolution computation "
02097                                 "to QC log file");
02098             }
02099 
02100             fors_qc_end_group();
02101 
02102         }
02103 
02104         if (dfs_save_table(frameset, restab, spectral_resolution_tag, qclist,
02105                            parlist, recipe, version))
02106             vimos_calib_exit(NULL);
02107 
02108         cpl_table_delete(restab); restab = NULL;
02109         cpl_propertylist_delete(qclist); qclist = NULL;
02110 
02111     }
02112     else
02113         vimos_calib_exit("Cannot compute the spectral resolution table");
02114 
02115     cpl_free(instrume); instrume = NULL;
02116     cpl_vector_delete(lines); lines = NULL;
02117 
02118     if (dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL,
02119                        parlist, recipe, version))
02120         vimos_calib_exit(NULL);
02121 
02122     /*
02123      * Global distortion models
02124      */
02125 
02126     if (slit_ident) {
02127 
02128         cpl_msg_info(recipe, "Computing global distortions model");
02129         global = mos_global_distortion(slits, maskslits, idscoeff, 
02130                                        polytraces, reference);
02131 
02132         if (global && 0) {
02133             cpl_table *stest;
02134             cpl_table *ctest;
02135             cpl_table *dtest;
02136             cpl_image *itest;
02137 
02138             stest = mos_build_slit_location(global, maskslits, ccd_ysize);
02139 
02140             ctest = mos_build_curv_coeff(global, maskslits, stest);
02141             if (dfs_save_table(frameset, ctest, "CURVS", NULL,
02142                                parlist, recipe, version))
02143                 vimos_calib_exit(NULL);
02144 
02145             itest = mos_spatial_calibration(spectra, stest, ctest, 
02146                                             reference, startwavelength, 
02147                                             endwavelength, dispersion, 
02148                                             0, NULL);
02149             cpl_table_delete(ctest); ctest = NULL;
02150             cpl_image_delete(itest); itest = NULL;
02151             if (dfs_save_table(frameset, stest, "SLITS", NULL,
02152                                parlist, recipe, version))
02153                 vimos_calib_exit(NULL);
02154 
02155             dtest = mos_build_disp_coeff(global, stest);
02156             if (dfs_save_table(frameset, dtest, "DISPS", NULL,
02157                                parlist, recipe, version))
02158                 vimos_calib_exit(NULL);
02159 
02160             cpl_table_delete(dtest); dtest = NULL;
02161             cpl_table_delete(stest); stest = NULL;
02162         }
02163 
02164         if (global) {
02165             if (dfs_save_table(frameset, global, global_distortion_tag, NULL,
02166                                parlist, recipe, version))
02167                 vimos_calib_exit(NULL);
02168             cpl_table_delete(global); global = NULL;
02169         }
02170 
02171         cpl_image_delete(spectra); spectra = NULL;
02172         cpl_table_delete(maskslits); maskslits = NULL;
02173     }
02174 
02175     cpl_table_delete(idscoeff); idscoeff = NULL;
02176 
02177     header = cpl_propertylist_new();
02178     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
02179     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
02180     cpl_propertylist_update_double(header, "CRVAL1", 
02181                                    startwavelength + dispersion/2);
02182     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
02183     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
02184     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
02185     cpl_propertylist_update_double(header, "CD1_1", dispersion);
02186     cpl_propertylist_update_double(header, "CD1_2", 0.0);
02187     cpl_propertylist_update_double(header, "CD2_1", 0.0);
02188     cpl_propertylist_update_double(header, "CD2_2", 1.0);
02189     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
02190     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
02191     cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);
02192 
02193     if (dfs_save_image(frameset, rectified, reduced_lamp_tag, header,
02194                        parlist, recipe, version))
02195         vimos_calib_exit(NULL);
02196 
02197     cpl_image_delete(rectified); rectified = NULL;
02198     cpl_propertylist_delete(header); header = NULL;
02199 
02200     if (check) {
02201         save_header = dfs_load_header(frameset, arc_tag, 0);
02202 
02203         cpl_propertylist_update_double(save_header, "CRPIX2", 1.0);
02204         cpl_propertylist_update_double(save_header, "CRVAL2", 1.0);
02205         /* cpl_propertylist_update_double(save_header, "CDELT2", 1.0); */
02206         cpl_propertylist_update_double(save_header, "CD1_1", 1.0);
02207         cpl_propertylist_update_double(save_header, "CD1_2", 0.0);
02208         cpl_propertylist_update_double(save_header, "CD2_1", 0.0);
02209         cpl_propertylist_update_double(save_header, "CD2_2", 1.0);
02210         cpl_propertylist_update_string(save_header, "CTYPE1", "LINEAR");
02211         cpl_propertylist_update_string(save_header, "CTYPE2", "PIXEL");
02212 
02213         if (dfs_save_image(frameset, residual, disp_residuals_tag, save_header,
02214                            parlist, recipe, version))
02215             vimos_calib_exit(NULL);
02216 
02217         cpl_image_delete(residual); residual = NULL;
02218         cpl_propertylist_delete(save_header); save_header = NULL;
02219     }
02220 
02221     wavemap = mos_map_wavelengths(coordinate, rainbow, slits, polytraces, 
02222                                   reference, startwavelength, endwavelength, 
02223                                   dispersion);
02224 
02225     cpl_image_delete(rainbow); rainbow = NULL;
02226 
02227     save_header = dfs_load_header(frameset, arc_tag, 0);
02228 
02229     cpl_image_turn(wavemap, rotate_back);
02230     if (dfs_save_image(frameset, wavemap, wavelength_map_tag, save_header,
02231                        parlist, recipe, version))
02232         vimos_calib_exit(NULL);
02233 
02234     cpl_image_delete(wavemap); wavemap = NULL;
02235 
02236     cpl_image_turn(coordinate, rotate_back);
02237     if (dfs_save_image(frameset, coordinate, spatial_map_tag, save_header,
02238                        parlist, recipe, version))
02239         vimos_calib_exit(NULL);
02240 
02241     cpl_image_delete(coordinate); coordinate = NULL;
02242     cpl_propertylist_delete(save_header); save_header = NULL;
02243 
02244     if (dfs_save_table(frameset, polytraces, curv_coeff_tag, NULL,
02245                        parlist, recipe, version))
02246         vimos_calib_exit(NULL);
02247 
02248     cpl_table_delete(polytraces); polytraces = NULL;
02249 
02250     mos_rotate_slits(slits, rotate, ccd_ysize, ccd_xsize);
02251     if (dfs_save_table(frameset, slits, slit_location_tag, NULL,
02252                        parlist, recipe, version))
02253         vimos_calib_exit(NULL);
02254 
02255     cpl_table_delete(slits); slits = NULL;
02256 
02257     if (cpl_error_get_code()) {
02258         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
02259         vimos_calib_exit(NULL);
02260     }
02261 
02262     return 0;
02263 }

Generated on Wed Jun 24 14:11:19 2009 for FORS Pipeline Reference Manual by  doxygen 1.4.7