vimos_calib.c

00001 /* $Id: vimos_calib.c,v 1.6 2008/08/05 14:21:35 cizzo Exp $
00002  *
00003  * This file is part of the FORS Data Reduction Pipeline
00004  * Copyright (C) 2006 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2008/08/05 14:21:35 $
00024  * $Revision: 1.6 $
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         rectified = mos_wavelength_calibration_raw(spectra, lines, dispersion,
01285                                                    peakdetection, wradius,
01286                                                    wdegree, wreject, reference,
01287                                                    &startwavelength, 
01288                                                    &endwavelength, fitlines, 
01289                                                    fiterror, idscoeff, wavemap,
01290                                                    residual, restable, NULL);
01291 
01292         if (rectified == NULL)
01293             vimos_calib_exit("Wavelength calibration failure.");
01294 
01295         if (dfs_save_table(frameset, restable, disp_residuals_table_tag, NULL,
01296                            parlist, recipe, version))
01297             vimos_calib_exit(NULL);
01298 
01299         cpl_table_delete(restable); restable = NULL;
01300 
01301         if (wmode) {
01302             cpl_image_delete(rectified); rectified = NULL;
01303             cpl_image_delete(wavemap); wavemap = NULL;
01304             mos_interpolate_wavecalib(idscoeff, wavemap, wmode);
01305             wavemap = mos_map_idscoeff(idscoeff, nx, reference,
01306                                        startwavelength, endwavelength);
01307             rectified = mos_wavelength_calibration(spectra, reference,
01308                                                    startwavelength, 
01309                                                    endwavelength, dispersion, 
01310                                                    idscoeff, 0);
01311         }
01312 
01313         cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
01314         cpl_table_set_column_unit(idscoeff, "error", "pixel");
01315         cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;
01316 
01317         for (i = 0; i < ny; i++)
01318             if (!cpl_table_is_valid(idscoeff, "c0", i))
01319                 cpl_table_set_invalid(idscoeff, "error", i);
01320 
01321         delta = mos_map_pixel(idscoeff, reference, startwavelength,
01322                               endwavelength, dispersion, 2);
01323 
01324         header = cpl_propertylist_new();
01325         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01326         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01327         cpl_propertylist_update_double(header, "CRVAL1",
01328                                        startwavelength + dispersion/2);
01329         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01330         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01331         cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01332         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01333         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01334         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01335         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01336         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01337         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01338 
01339         if (dfs_save_image(frameset, delta, delta_image_tag,
01340                            header, parlist, recipe, version))
01341             vimos_calib_exit(NULL);
01342 
01343         cpl_image_delete(delta); delta = NULL;
01344         cpl_propertylist_delete(header); header = NULL;
01345 
01346         cpl_msg_info(recipe, "Valid solutions found: %d out of %d rows", 
01347                      ny - cpl_table_count_invalid(idscoeff, "c0"), ny);
01348 
01349         cpl_image_delete(spectra); spectra = NULL;
01350 
01351         mean_rms = mos_distortions_rms(rectified, lines, startwavelength,
01352                                        dispersion, 6, 0);
01353 
01354         cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);
01355 
01356         mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01357 
01358         cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)",
01359                      mean_rms, mean_rms * dispersion);
01360 
01361         restab = mos_resolution_table(rectified, startwavelength, dispersion,
01362                                       60000, lines);
01363 
01364         if (restab) {
01365             cpl_msg_info(recipe, "Mean spectral resolution: %.2f",
01366                   cpl_table_get_column_mean(restab, "resolution"));
01367             cpl_msg_info(recipe, 
01368                   "Mean reference lines FWHM: %.2f +/- %.2f pixel",
01369                   cpl_table_get_column_mean(restab, "fwhm") / dispersion,
01370                   cpl_table_get_column_mean(restab, "fwhm_rms") / dispersion);
01371 
01372             if (qc) {
01373 
01374                 header = dfs_load_header(frameset, arc_tag, 0);
01375 
01376                 if (header == NULL)
01377                     vimos_calib_exit("Cannot reload arc lamp header");
01378 
01379                 qclist = cpl_propertylist_new();
01380     
01381                 fors_qc_start_group(qclist, "1.1", instrume);
01382 
01383 
01384                 /*
01385                  * QC1 group header
01386                  */
01387 
01388                 if (fors_qc_write_string("PRO.CATG", spectral_resolution_tag,
01389                                         "Product category", instrume))
01390                     vimos_calib_exit("Cannot write product category to "
01391                                     "QC log file");
01392 
01393                 if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL,
01394                                           "DPR type", instrume))
01395                     vimos_calib_exit("Missing keyword DPR TYPE in arc "
01396                                     "lamp header");
01397 
01398                 if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL,
01399                                           "Template", instrume))
01400                     vimos_calib_exit("Missing keyword TPL ID in arc "
01401                                     "lamp header");
01402 
01403                 if (fors_qc_keyword_to_paf(header, key_gris_name, NULL,
01404                                           "Grism name", instrume))
01405                     vimos_calib_exit("Missing keyword INS GRISx NAME in arc "
01406                                     "lamp header");
01407 
01408 
01409                 if (fors_qc_keyword_to_paf(header, key_gris_id, NULL,
01410                                           "Grism identifier", instrume)) {
01411                     cpl_msg_error(recipe, "Missing keyword %s in arc "
01412                                     "lamp header", key_gris_id);
01413                     vimos_calib_exit(NULL);
01414                 }
01415 
01416                 if (cpl_propertylist_has(header, key_filt_name))
01417                     fors_qc_keyword_to_paf(header, key_filt_name, NULL,
01418                                           "Filter name", instrume);
01419 
01420                 if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
01421                                           "Collimator name", instrume))
01422                     vimos_calib_exit("Missing keyword INS COLL NAME in arc "
01423                                     "lamp header");
01424 
01425                 if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
01426                                           "Chip identifier", instrume))
01427                     vimos_calib_exit("Missing keyword DET CHIP1 ID in arc "
01428                                     "lamp header");
01429 
01430                 if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
01431                                           "Archive name of input data",
01432                                           instrume))
01433                     vimos_calib_exit("Missing keyword ARCFILE in arc "
01434                                     "lamp header");
01435 
01436                 cpl_propertylist_delete(header); header = NULL;
01437 
01438                 pipefile = dfs_generate_filename_tfits(spectral_resolution_tag);
01439                 if (fors_qc_write_string("PIPEFILE", pipefile,
01440                                         "Pipeline product name", instrume))
01441                     vimos_calib_exit("Cannot write PIPEFILE to QC log file");
01442                 cpl_free(pipefile); pipefile = NULL;
01443 
01444 
01445                 /*
01446                  * QC1 parameters
01447                  */
01448     
01449                 keyname = "QC.MOS.RESOLUTION";
01450 
01451                 if (fors_qc_write_qc_double(qclist, 
01452                                            cpl_table_get_column_mean(restab,
01453                                                                  "resolution"),
01454                                            keyname, NULL, 
01455                                            "Mean spectral resolution",
01456                                            instrume)) {
01457                     vimos_calib_exit("Cannot write mean spectral resolution to "
01458                                     "QC log file");
01459                 }
01460 
01461                 keyname = "QC.MOS.RESOLUTION.RMS"; 
01462 
01463                 if (fors_qc_write_qc_double(qclist,
01464                                            cpl_table_get_column_stdev(restab,
01465                                                                 "resolution"),
01466                                            keyname, NULL, 
01467                                            "Scatter of spectral resolution",
01468                                            instrume)) {
01469                     vimos_calib_exit("Cannot write spectral resolution scatter "
01470                                     "to QC log file");
01471                 }
01472 
01473                 keyname = "QC.MOS.RESOLUTION.NLINES";
01474 
01475                 if (fors_qc_write_qc_int(qclist, cpl_table_get_nrow(restab) -
01476                                         cpl_table_count_invalid(restab,
01477                                                                 "resolution"),
01478                                         keyname, NULL,
01479                                         "Number of lines for spectral "
01480                                         "resolution computation",
01481                                         instrume)) {
01482                     vimos_calib_exit("Cannot write number of lines used in "
01483                                     "spectral resolution computation "
01484                                     "to QC log file");
01485                 }
01486 
01487                 fors_qc_end_group();
01488 
01489             }  /* End of QC1 computation */
01490 
01491             if (dfs_save_table(frameset, restab, spectral_resolution_tag, 
01492                                qclist, parlist, recipe, version))
01493                 vimos_calib_exit(NULL);
01494 
01495             cpl_table_delete(restab); restab = NULL;
01496             cpl_propertylist_delete(qclist); qclist = NULL;
01497 
01498         }
01499         else
01500             vimos_calib_exit("Cannot compute the spectral resolution table");
01501 
01502         cpl_vector_delete(lines); lines = NULL;
01503 
01504 
01505         /*
01506          * Save rectified arc lamp spectrum to disk
01507          */
01508 
01509         header = cpl_propertylist_new();
01510         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01511         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01512         cpl_propertylist_update_double(header, "CRVAL1", 
01513                                        startwavelength + dispersion/2);
01514         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01515         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01516         cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01517         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01518         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01519         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01520         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01521         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01522         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01523         cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);
01524 
01525         if (dfs_save_image(frameset, rectified, reduced_lamp_tag, header, 
01526                            parlist, recipe, version))
01527             vimos_calib_exit(NULL);
01528 
01529         cpl_image_delete(rectified); rectified = NULL;
01530         cpl_propertylist_delete(header); header = NULL;
01531 
01532         if (dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL, 
01533                            parlist, recipe, version))
01534             vimos_calib_exit(NULL);
01535 
01536         cpl_table_delete(idscoeff); idscoeff = NULL;
01537 
01538         header = dfs_load_header(frameset, arc_tag, 0);
01539 
01540         if (header == NULL)
01541             vimos_calib_exit("Cannot reload arc lamp header");
01542 
01543         if (qc) {
01544 
01545             compute_qc = 0;
01546             if (fabs(mxpos) < 0.05)
01547                 compute_qc = 1;
01548 
01549             if (compute_qc) {
01550 
01551                 fors_qc_start_group(header, "1.1", instrume);
01552 
01553                 /*
01554                  * QC1 group header
01555                  */
01556 
01557                 if (fors_qc_write_string("PRO.CATG", wavelength_map_tag,
01558                                         "Product category", instrume))
01559                     vimos_calib_exit("Cannot write product category to "
01560                                     "QC log file");
01561 
01562                 if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL, 
01563                                           "DPR type", instrume))
01564                     vimos_calib_exit("Missing keyword DPR TYPE in arc "
01565                                     "lamp header");
01566 
01567                 if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL, 
01568                                           "Template", instrume))
01569                     vimos_calib_exit("Missing keyword TPL ID in arc "
01570                                     "lamp header");
01571 
01572                 if (fors_qc_keyword_to_paf(header, key_gris_name, NULL, 
01573                                           "Grism name", instrume)) {
01574                     cpl_msg_error(recipe, "Missing keyword %s in arc "
01575                                     "lamp header", key_gris_name);
01576                     vimos_calib_exit(NULL);
01577                 }
01578 
01579                 if (fors_qc_keyword_to_paf(header, key_gris_id, NULL,
01580                                           "Grism identifier", instrume)) {
01581                     cpl_msg_error(recipe, "Missing keyword %s in arc "
01582                                     "lamp header", key_gris_id);
01583                     vimos_calib_exit(NULL);
01584                 }
01585 
01586                 if (cpl_propertylist_has(header, key_filt_name))
01587                     fors_qc_keyword_to_paf(header, key_filt_name, NULL,
01588                                           "Filter name", instrume);
01589 
01590                 if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
01591                                           "Collimator name", instrume))
01592                     vimos_calib_exit("Missing keyword INS COLL NAME in arc "
01593                                     "lamp header");
01594 
01595                 if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
01596                                           "Chip identifier", instrume))
01597                     vimos_calib_exit("Missing keyword DET CHIP1 ID in arc "
01598                                     "lamp header");
01599 
01600                 if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
01601                                           "Archive name of input data", 
01602                                           instrume))
01603                     vimos_calib_exit("Missing keyword ARCFILE in arc "
01604                                     "lamp header");
01605     
01606                 pipefile = dfs_generate_filename(wavelength_map_tag);
01607                 if (fors_qc_write_string("PIPEFILE", pipefile, 
01608                                         "Pipeline product name", instrume))
01609                     vimos_calib_exit("Cannot write PIPEFILE to QC log file");
01610                 cpl_free(pipefile); pipefile = NULL;
01611 
01612 
01613                 /*
01614                  * QC1 parameters
01615                  */
01616 
01617                 data = cpl_image_get_data(wavemap);
01618     
01619                 if (fors_qc_write_qc_double(header, data[nx/2 + ny*nx/2],
01620                                            "QC.MOS.CENTRAL.WAVELENGTH",
01621                                            "Angstrom", 
01622                                            "Wavelength at CCD center",
01623                                            instrume)) {
01624                     vimos_calib_exit("Cannot write central wavelength to QC "
01625                                     "log file");
01626                 }
01627 
01628                 fors_qc_end_group();
01629 
01630             }
01631 
01632         }
01633 
01634         cpl_image_turn(wavemap, rotate_back);
01635         if (dfs_save_image(frameset, wavemap, wavelength_map_tag, header,
01636                            parlist, recipe, version))
01637             vimos_calib_exit(NULL);
01638 
01639         cpl_image_delete(wavemap); wavemap = NULL;
01640 
01641         if (check) {
01642             cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01643             cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01644             /* cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01645             cpl_propertylist_update_double(header, "CD1_1", 1.0);
01646             cpl_propertylist_update_double(header, "CD1_2", 0.0);
01647             cpl_propertylist_update_double(header, "CD2_1", 0.0);
01648             cpl_propertylist_update_double(header, "CD2_2", 1.0);
01649             cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01650             cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01651 
01652             if (dfs_save_image(frameset, residual, disp_residuals_tag, header,
01653                                parlist, recipe, version))
01654                 vimos_calib_exit(NULL);
01655 
01656             cpl_image_delete(residual); residual = NULL;
01657         }
01658 
01659         cpl_propertylist_delete(header); header = NULL;
01660 
01661         return 0;         /* Successful LSS data reduction */
01662 
01663     }   /* End of LSS data reduction section */
01664 
01665 
01666     /*
01667      * Here the MOS calibration is carried out.
01668      */
01669 
01670     /*
01671      * Detecting spectra on the CCD
01672      */
01673 
01674     cpl_msg_indent_less();
01675     cpl_msg_info(recipe, "Detecting spectra on CCD...");
01676     cpl_msg_indent_more();
01677 
01678     ccd_xsize = nx = cpl_image_get_size_x(spectra);
01679     ccd_ysize = ny = cpl_image_get_size_y(spectra);
01680 
01681     refmask = cpl_mask_new(nx, ny);
01682 
01683     checkwave = mos_wavelength_calibration_raw(spectra, lines, dispersion, 
01684                                                peakdetection, wradius, 
01685                                                wdegree, wreject, reference,
01686                                                &startwavelength, &endwavelength,
01687                                                NULL, NULL, NULL, NULL, NULL, 
01688                                                NULL, refmask);
01689 
01690     if (checkwave == NULL)
01691         vimos_calib_exit("Wavelength calibration failure.");
01692 
01693     /*
01694      * Save check image to disk
01695      */
01696 
01697     header = cpl_propertylist_new();
01698     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01699     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01700     cpl_propertylist_update_double(header, "CRVAL1", 
01701                                    startwavelength + dispersion/2);
01702     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01703     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01704     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01705     cpl_propertylist_update_double(header, "CD1_1", dispersion);
01706     cpl_propertylist_update_double(header, "CD1_2", 0.0);
01707     cpl_propertylist_update_double(header, "CD2_1", 0.0);
01708     cpl_propertylist_update_double(header, "CD2_2", 1.0);
01709     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01710     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01711 
01712     if (check) {
01713         if (dfs_save_image(frameset, checkwave, spectra_detection_tag, header, 
01714                            parlist, recipe, version))
01715             vimos_calib_exit(NULL);
01716     }
01717 
01718     cpl_image_delete(checkwave); checkwave = NULL;
01719     cpl_propertylist_delete(header); header = NULL;
01720 
01721     cpl_msg_info(recipe, "Locate slits at reference wavelength on CCD...");
01722     slits = mos_locate_spectra(refmask);
01723 
01724     if (!slits) {
01725         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
01726         vimos_calib_exit("No slits could be detected!");
01727     }
01728 
01729     refimage = cpl_image_new_from_mask(refmask);
01730     cpl_mask_delete(refmask); refmask = NULL;
01731 
01732     if (check) {
01733         save_header = dfs_load_header(frameset, arc_tag, 0);
01734         cpl_image_turn(refimage, rotate_back);
01735         if (dfs_save_image(frameset, refimage, slit_map_tag, NULL,
01736                            parlist, recipe, version))
01737             vimos_calib_exit(NULL);
01738         cpl_propertylist_delete(save_header); save_header = NULL;
01739     }
01740 
01741     cpl_image_delete(refimage); refimage = NULL;
01742 
01743     if (slit_ident) {
01744 
01745         /*
01746          * Attempt slit identification: this recipe may continue even
01747          * in case of failed identification (i.e., the position table is 
01748          * not produced, but an error is not set). In case of failure,
01749          * the spectra would be still extracted, even if they would not
01750          * be associated to slits on the mask.
01751          * 
01752          * The reason for making the slit identification an user option 
01753          * (via the parameter slit_ident) is to offer the possibility 
01754          * to avoid identifications that are only apparently successful, 
01755          * as it would happen in the case of an incorrect slit description 
01756          * in the data header.
01757          */
01758 
01759         cpl_msg_indent_less();
01760         cpl_msg_info(recipe, "Attempt slit identification (optional)...");
01761         cpl_msg_indent_more();
01762 
01763         mos_rotate_slits(maskslits, -rotate, 0, 0);
01764         positions = mos_identify_slits(slits, maskslits, NULL);
01765 
01766         if (positions) {
01767             cpl_table_delete(slits);
01768             slits = positions;
01769 
01770             /*
01771              * Eliminate slits which are _entirely_ outside the CCD
01772              */
01773 
01774             cpl_table_and_selected_double(slits, 
01775                                           "ybottom", CPL_GREATER_THAN, ny-1);
01776             cpl_table_or_selected_double(slits, 
01777                                           "ytop", CPL_LESS_THAN, 0);
01778             cpl_table_erase_selected(slits);
01779 
01780             nslits = cpl_table_get_nrow(slits);
01781 
01782             if (nslits == 0)
01783                 vimos_calib_exit("No slits found on the CCD");
01784 
01785             cpl_msg_info(recipe, "%d slits are entirely or partially "
01786                          "contained in CCD", nslits);
01787 
01788         }
01789         else {
01790             slit_ident = 0;
01791             cpl_msg_info(recipe, "Global distortion model cannot be computed");
01792             if (cpl_error_get_code() != CPL_ERROR_NONE) {
01793                 vimos_calib_exit(NULL);
01794             }
01795         }
01796     }
01797 
01798 
01799     /*
01800      * Determination of spectral curvature
01801      */
01802 
01803     cpl_msg_indent_less();
01804     cpl_msg_info(recipe, "Determining spectral curvature...");
01805     cpl_msg_indent_more();
01806 
01807     cpl_msg_info(recipe, "Tracing master flat field spectra edges...");
01808     traces = mos_trace_flat(master_flat, slits, reference, 
01809                             startwavelength, endwavelength, dispersion);
01810 
01811     if (!traces)
01812         vimos_calib_exit("Tracing failure");
01813 
01814     cpl_msg_info(recipe, "Fitting flat field spectra edges...");
01815     polytraces = mos_poly_trace(slits, traces, cdegree);
01816 
01817     if (!polytraces)
01818         vimos_calib_exit("Trace fitting failure");
01819 
01820     if (cmode) {
01821         cpl_msg_info(recipe, "Computing global spectral curvature model...");
01822         mos_global_trace(slits, polytraces, cmode);
01823     }
01824 
01825     if (dfs_save_table(frameset, traces, curv_traces_tag, NULL, parlist,
01826                        recipe, version))
01827         vimos_calib_exit(NULL);
01828 
01829     cpl_table_delete(traces); traces = NULL;
01830 
01831     coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01832     spatial = mos_spatial_calibration(spectra, slits, polytraces, reference, 
01833                                       startwavelength, endwavelength, 
01834                                       dispersion, 0, coordinate);
01835 
01836     if (!slit_ident) {
01837         cpl_image_delete(spectra); spectra = NULL;
01838     }
01839 
01840     /*
01841      * Flat field normalisation is done directly on the master flat
01842      * field (without spatial rectification first). The spectral
01843      * curvature model may be provided in input, in future releases.
01844      */
01845 
01846     cpl_msg_indent_less();
01847     cpl_msg_info(recipe, "Perform flat field normalisation...");
01848     cpl_msg_indent_more();
01849 
01850     norm_flat = cpl_image_duplicate(master_flat);
01851 
01852     smo_flat = mos_normalise_flat(norm_flat, coordinate, slits, polytraces, 
01853                                   reference, startwavelength, endwavelength,
01854                                   dispersion, dradius, ddegree);
01855 
01856     cpl_image_delete(smo_flat); smo_flat = NULL;  /* It may be a product */
01857  
01858     save_header = dfs_load_header(frameset, flat_tag, 0);
01859     cpl_propertylist_update_int(save_header, "ESO PRO DATANCOM", nflats);
01860 
01861     cpl_image_turn(master_flat, rotate_back);
01862     if (dfs_save_image(frameset, master_flat, master_screen_flat_tag,
01863                        save_header, parlist, recipe, version))
01864         vimos_calib_exit(NULL);
01865 
01866     cpl_image_delete(master_flat); master_flat = NULL;
01867 
01868     cpl_image_turn(norm_flat, rotate_back);
01869     if (dfs_save_image(frameset, norm_flat, master_norm_flat_tag,
01870                        save_header, parlist, recipe, version))
01871         vimos_calib_exit(NULL);
01872 
01873     cpl_image_delete(norm_flat); norm_flat = NULL;
01874     cpl_propertylist_delete(save_header); save_header = NULL;
01875 
01876 
01877     /*
01878      * Final wavelength calibration of spectra having their curvature
01879      * removed
01880      */
01881 
01882     cpl_msg_indent_less();
01883     cpl_msg_info(recipe, "Perform final wavelength calibration...");
01884     cpl_msg_indent_more();
01885 
01886     nx = cpl_image_get_size_x(spatial);
01887     ny = cpl_image_get_size_y(spatial);
01888 
01889     idscoeff = cpl_table_new(ny);
01890     restable = cpl_table_new(nlines);
01891     rainbow = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01892     if (check)
01893         residual = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01894     fiterror = cpl_calloc(ny, sizeof(double));
01895     fitlines = cpl_calloc(ny, sizeof(int));
01896 
01897     rectified = mos_wavelength_calibration_final(spatial, slits, lines, 
01898                                                  dispersion, peakdetection, 
01899                                                  wradius, wdegree, wreject,
01900                                                  reference, &startwavelength, 
01901                                                  &endwavelength, fitlines, 
01902                                                  fiterror, idscoeff, rainbow, 
01903                                                  residual, restable);
01904 
01905     cpl_image_delete(spatial); spatial = NULL;
01906 
01907     if (rectified == NULL)
01908         vimos_calib_exit("Wavelength calibration failure.");
01909 
01910     if (dfs_save_table(frameset, restable, disp_residuals_table_tag, NULL,
01911                        parlist, recipe, version))
01912         vimos_calib_exit(NULL);
01913 
01914     cpl_table_delete(restable); restable = NULL;
01915 
01916     cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
01917     cpl_table_set_column_unit(idscoeff, "error", "pixel");
01918     cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;
01919 
01920     for (i = 0; i < ny; i++)
01921         if (!cpl_table_is_valid(idscoeff, "c0", i))
01922             cpl_table_set_invalid(idscoeff, "error", i);
01923 
01924     delta = mos_map_pixel(idscoeff, reference, startwavelength,
01925                           endwavelength, dispersion, 2);
01926 
01927     header = cpl_propertylist_new();
01928     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01929     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01930     cpl_propertylist_update_double(header, "CRVAL1",
01931                                    startwavelength + dispersion/2);
01932     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01933     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01934     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01935     cpl_propertylist_update_double(header, "CD1_1", dispersion);
01936     cpl_propertylist_update_double(header, "CD1_2", 0.0);
01937     cpl_propertylist_update_double(header, "CD2_1", 0.0);
01938     cpl_propertylist_update_double(header, "CD2_2", 1.0);
01939     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01940     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01941 
01942     if (dfs_save_image(frameset, delta, delta_image_tag,
01943                        header, parlist, recipe, version))
01944         vimos_calib_exit(NULL);
01945 
01946     cpl_image_delete(delta); delta = NULL;
01947     cpl_propertylist_delete(header); header = NULL;
01948 
01949     mean_rms = mos_distortions_rms(rectified, lines, startwavelength, 
01950                                    dispersion, 6, 0);
01951 
01952     cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);
01953 
01954     mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01955 
01956     cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)", 
01957                  mean_rms, mean_rms * dispersion);
01958 
01959     restab = mos_resolution_table(rectified, startwavelength, dispersion, 
01960                                   60000, lines);
01961 
01962     if (restab) {
01963         cpl_msg_info(recipe, "Mean spectral resolution: %.2f", 
01964                    cpl_table_get_column_mean(restab, "resolution"));
01965         cpl_msg_info(recipe, "Mean reference lines FWHM: %.2f +/- %.2f pixel",
01966                    cpl_table_get_column_mean(restab, "fwhm") / dispersion,
01967                    cpl_table_get_column_mean(restab, "fwhm_rms") / dispersion);
01968 
01969         if (qc) {
01970 
01971             header = dfs_load_header(frameset, arc_tag, 0);
01972 
01973             if (header == NULL)
01974                 vimos_calib_exit("Cannot reload arc lamp header");
01975 
01976             qclist = cpl_propertylist_new();
01977 
01978             fors_qc_start_group(qclist, "1.1", instrume);
01979 
01980 
01981             /*
01982              * QC1 group header
01983              */
01984 
01985             if (fors_qc_write_string("PRO.CATG", spectral_resolution_tag,
01986                                     "Product category", instrume))
01987                 vimos_calib_exit("Cannot write product category to "
01988                                 "QC log file");
01989 
01990             if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL,
01991                                       "DPR type", instrume))
01992                 vimos_calib_exit("Missing keyword DPR TYPE in arc "
01993                                 "lamp header");
01994 
01995             if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL,
01996                                       "Template", instrume))
01997                 vimos_calib_exit("Missing keyword TPL ID in arc "
01998                                 "lamp header");
01999 
02000             if (fors_qc_keyword_to_paf(header, key_gris_name, NULL,
02001                                       "Grism name", instrume)) {
02002                 cpl_msg_error(recipe, "Missing keyword %s in arc "
02003                                 "lamp header", key_gris_name);
02004                 vimos_calib_exit(NULL);
02005             }
02006 
02007             if (fors_qc_keyword_to_paf(header, key_gris_id, NULL,
02008                                       "Grism identifier", instrume)) {
02009                 cpl_msg_error(recipe, "Missing keyword %s in arc "
02010                                 "lamp header", key_gris_id);
02011                 vimos_calib_exit(NULL);
02012             }
02013 
02014             if (cpl_propertylist_has(header, key_filt_name))
02015                 fors_qc_keyword_to_paf(header, key_filt_name, NULL,
02016                                       "Filter name", instrume);
02017 
02018             if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
02019                                       "Collimator name", instrume))
02020                 vimos_calib_exit("Missing keyword INS COLL NAME in arc "
02021                                 "lamp header");
02022 
02023             if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
02024                                       "Chip identifier", instrume))
02025                 vimos_calib_exit("Missing keyword DET CHIP1 ID in arc "
02026                                 "lamp header");
02027 
02028             if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
02029                                       "Archive name of input data", 
02030                                       instrume))
02031                 vimos_calib_exit("Missing keyword ARCFILE in arc "
02032                                 "lamp header");
02033 
02034             cpl_propertylist_delete(header); header = NULL;
02035 
02036             pipefile = dfs_generate_filename_tfits(spectral_resolution_tag);
02037             if (fors_qc_write_string("PIPEFILE", pipefile,
02038                                     "Pipeline product name", instrume))
02039                 vimos_calib_exit("Cannot write PIPEFILE to QC log file");
02040             cpl_free(pipefile); pipefile = NULL;
02041 
02042 
02043             /*
02044              * QC1 parameters
02045              */
02046 
02047             keyname = "QC.MOS.RESOLUTION";
02048 
02049             if (fors_qc_write_qc_double(qclist, 
02050                                        cpl_table_get_column_mean(restab,
02051                                                                  "resolution"),
02052                                        keyname,
02053                                        "Angstrom",
02054                                        "Mean spectral resolution",
02055                                        instrume)) {
02056                 vimos_calib_exit("Cannot write mean spectral resolution to QC "
02057                                 "log file");
02058             }
02059 
02060             keyname = "QC.MOS.RESOLUTION.RMS";
02061 
02062             if (fors_qc_write_qc_double(qclist, 
02063                                        cpl_table_get_column_stdev(restab, 
02064                                                                   "resolution"),
02065                                        "QC.MOS.RESOLUTION.RMS",
02066                                        "Angstrom", 
02067                                        "Scatter of spectral resolution",
02068                                        instrume)) {
02069                 vimos_calib_exit("Cannot write spectral resolution scatter "
02070                                 "to QC log file");
02071             }
02072 
02073             keyname = "QC.MOS.RESOLUTION.NLINES";
02074 
02075             if (fors_qc_write_qc_int(qclist, cpl_table_get_nrow(restab) -
02076                                     cpl_table_count_invalid(restab, 
02077                                                             "resolution"),
02078                                     "QC.MOS.RESOLUTION.NLINES",
02079                                     NULL,
02080                                     "Number of lines for spectral resolution "
02081                                     "computation",
02082                                     instrume)) {
02083                 vimos_calib_exit("Cannot write number of lines used in "
02084                                 "spectral resolution computation "
02085                                 "to QC log file");
02086             }
02087 
02088             fors_qc_end_group();
02089 
02090         }
02091 
02092         if (dfs_save_table(frameset, restab, spectral_resolution_tag, qclist,
02093                            parlist, recipe, version))
02094             vimos_calib_exit(NULL);
02095 
02096         cpl_table_delete(restab); restab = NULL;
02097         cpl_propertylist_delete(qclist); qclist = NULL;
02098 
02099     }
02100     else
02101         vimos_calib_exit("Cannot compute the spectral resolution table");
02102 
02103     cpl_free(instrume); instrume = NULL;
02104     cpl_vector_delete(lines); lines = NULL;
02105 
02106     if (dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL,
02107                        parlist, recipe, version))
02108         vimos_calib_exit(NULL);
02109 
02110     /*
02111      * Global distortion models
02112      */
02113 
02114     if (slit_ident) {
02115 
02116         cpl_msg_info(recipe, "Computing global distortions model");
02117         global = mos_global_distortion(slits, maskslits, idscoeff, 
02118                                        polytraces, reference);
02119 
02120         if (global && 0) {
02121             cpl_table *stest;
02122             cpl_table *ctest;
02123             cpl_table *dtest;
02124             cpl_image *itest;
02125 
02126             stest = mos_build_slit_location(global, maskslits, ccd_ysize);
02127 
02128             ctest = mos_build_curv_coeff(global, maskslits, stest);
02129             if (dfs_save_table(frameset, ctest, "CURVS", NULL,
02130                                parlist, recipe, version))
02131                 vimos_calib_exit(NULL);
02132 
02133             itest = mos_spatial_calibration(spectra, stest, ctest, 
02134                                             reference, startwavelength, 
02135                                             endwavelength, dispersion, 
02136                                             0, NULL);
02137             cpl_table_delete(ctest); ctest = NULL;
02138             cpl_image_delete(itest); itest = NULL;
02139             if (dfs_save_table(frameset, stest, "SLITS", NULL,
02140                                parlist, recipe, version))
02141                 vimos_calib_exit(NULL);
02142 
02143             dtest = mos_build_disp_coeff(global, stest);
02144             if (dfs_save_table(frameset, dtest, "DISPS", NULL,
02145                                parlist, recipe, version))
02146                 vimos_calib_exit(NULL);
02147 
02148             cpl_table_delete(dtest); dtest = NULL;
02149             cpl_table_delete(stest); stest = NULL;
02150         }
02151 
02152         if (global) {
02153             if (dfs_save_table(frameset, global, global_distortion_tag, NULL,
02154                                parlist, recipe, version))
02155                 vimos_calib_exit(NULL);
02156             cpl_table_delete(global); global = NULL;
02157         }
02158 
02159         cpl_image_delete(spectra); spectra = NULL;
02160         cpl_table_delete(maskslits); maskslits = NULL;
02161     }
02162 
02163     cpl_table_delete(idscoeff); idscoeff = NULL;
02164 
02165     header = cpl_propertylist_new();
02166     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
02167     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
02168     cpl_propertylist_update_double(header, "CRVAL1", 
02169                                    startwavelength + dispersion/2);
02170     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
02171     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
02172     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
02173     cpl_propertylist_update_double(header, "CD1_1", dispersion);
02174     cpl_propertylist_update_double(header, "CD1_2", 0.0);
02175     cpl_propertylist_update_double(header, "CD2_1", 0.0);
02176     cpl_propertylist_update_double(header, "CD2_2", 1.0);
02177     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
02178     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
02179     cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);
02180 
02181     if (dfs_save_image(frameset, rectified, reduced_lamp_tag, header,
02182                        parlist, recipe, version))
02183         vimos_calib_exit(NULL);
02184 
02185     cpl_image_delete(rectified); rectified = NULL;
02186     cpl_propertylist_delete(header); header = NULL;
02187 
02188     if (check) {
02189         save_header = dfs_load_header(frameset, arc_tag, 0);
02190 
02191         cpl_propertylist_update_double(save_header, "CRPIX2", 1.0);
02192         cpl_propertylist_update_double(save_header, "CRVAL2", 1.0);
02193         /* cpl_propertylist_update_double(save_header, "CDELT2", 1.0); */
02194         cpl_propertylist_update_double(save_header, "CD1_1", 1.0);
02195         cpl_propertylist_update_double(save_header, "CD1_2", 0.0);
02196         cpl_propertylist_update_double(save_header, "CD2_1", 0.0);
02197         cpl_propertylist_update_double(save_header, "CD2_2", 1.0);
02198         cpl_propertylist_update_string(save_header, "CTYPE1", "LINEAR");
02199         cpl_propertylist_update_string(save_header, "CTYPE2", "PIXEL");
02200 
02201         if (dfs_save_image(frameset, residual, disp_residuals_tag, save_header,
02202                            parlist, recipe, version))
02203             vimos_calib_exit(NULL);
02204 
02205         cpl_image_delete(residual); residual = NULL;
02206         cpl_propertylist_delete(save_header); save_header = NULL;
02207     }
02208 
02209     wavemap = mos_map_wavelengths(coordinate, rainbow, slits, polytraces, 
02210                                   reference, startwavelength, endwavelength, 
02211                                   dispersion);
02212 
02213     cpl_image_delete(rainbow); rainbow = NULL;
02214 
02215     save_header = dfs_load_header(frameset, arc_tag, 0);
02216 
02217     cpl_image_turn(wavemap, rotate_back);
02218     if (dfs_save_image(frameset, wavemap, wavelength_map_tag, save_header,
02219                        parlist, recipe, version))
02220         vimos_calib_exit(NULL);
02221 
02222     cpl_image_delete(wavemap); wavemap = NULL;
02223 
02224     cpl_image_turn(coordinate, rotate_back);
02225     if (dfs_save_image(frameset, coordinate, spatial_map_tag, save_header,
02226                        parlist, recipe, version))
02227         vimos_calib_exit(NULL);
02228 
02229     cpl_image_delete(coordinate); coordinate = NULL;
02230     cpl_propertylist_delete(save_header); save_header = NULL;
02231 
02232     if (dfs_save_table(frameset, polytraces, curv_coeff_tag, NULL,
02233                        parlist, recipe, version))
02234         vimos_calib_exit(NULL);
02235 
02236     cpl_table_delete(polytraces); polytraces = NULL;
02237 
02238     mos_rotate_slits(slits, rotate, ccd_ysize, ccd_xsize);
02239     if (dfs_save_table(frameset, slits, slit_location_tag, NULL,
02240                        parlist, recipe, version))
02241         vimos_calib_exit(NULL);
02242 
02243     cpl_table_delete(slits); slits = NULL;
02244 
02245     if (cpl_error_get_code()) {
02246         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
02247         vimos_calib_exit(NULL);
02248     }
02249 
02250     return 0;
02251 }

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