00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00037 static int fors_science_create(cpl_plugin *);
00038 static int fors_science_exec(cpl_plugin *);
00039 static int fors_science_destroy(cpl_plugin *);
00040 static int fors_science(cpl_parameterlist *, cpl_frameset *);
00041
00042 static char fors_science_description[] =
00043 "This recipe is used to reduce scientific spectra using the extraction\n"
00044 "mask and the products created by the recipe fors_calib. The spectra are\n"
00045 "bias subtracted, flat fielded (if a normalised flat field is specified)\n"
00046 "and remapped eliminating the optical distortions. The wavelength calibration\n"
00047 "can be optionally upgraded using a number of sky lines: if no sky lines\n"
00048 "catalog of wavelengths is specified, an internal one is used instead.\n"
00049 "If the alignment to the sky lines is performed, the input dispersion\n"
00050 "coefficients table is upgraded and saved to disk, and a new CCD wavelengths\n"
00051 "map is created.\n"
00052 "This recipe accepts both FORS1 and FORS2 frames. A grism table (typically\n"
00053 "depending on the instrument mode, and in particular on the grism used)\n"
00054 "may also be specified: this table contains a default recipe parameter\n"
00055 "setting to control the way spectra are extracted for a specific instrument\n"
00056 "mode, as it is used for automatic run of the pipeline on Paranal and in\n"
00057 "Garching. If this table is specified, it will modify the default recipe\n"
00058 "parameter setting, with the exception of those parameters which have been\n"
00059 "explicitly modifyed on the command line. If a grism table is not specified,\n"
00060 "the input recipe parameters values will always be read from the command\n"
00061 "line, or from an esorex configuration file if present, or from their\n"
00062 "generic default values (that are rarely meaningful).\n"
00063 "In the table below the MXU acronym can be read alternatively as MOS\n"
00064 "and LSS, depending on the instrument mode of the input data. The acronym\n"
00065 "SCI on products should be read STD in case of standard stars observations\n"
00066 "Either a scientific or a standard star exposure can be specified in input.\n"
00067 "A CURV_COEFF table is not (yet) expected for LSS data.\n\n"
00068 "Input files:\n\n"
00069 " DO category: Type: Explanation: Required:\n"
00070 " SCIENCE_MXU Raw Scientific exposure Y\n"
00071 " or STANDARD_MXU Raw Standard star exposure Y\n"
00072 " MASTER_BIAS Calib Master bias Y\n"
00073 " GRISM_TABLE Calib Grism table .\n"
00074 " MASTER_SKYLINECAT Calib Sky lines catalog .\n"
00075 "\n"
00076 " MASTER_NORM_FLAT_MXU Calib Normalised flat field .\n"
00077 " DISP_COEFF_MXU Calib Inverse dispersion Y\n"
00078 " CURV_COEFF_MXU Calib Spectral curvature Y\n"
00079 " SLIT_LOCATION_MXU Calib Slits positions table Y\n"
00080 "\n"
00081 " or, in case of LSS-like MOS/MXU data,\n"
00082 "\n"
00083 " MASTER_NORM_FLAT_LONG_MXU Calib Normalised flat field .\n"
00084 " DISP_COEFF_LONG_MXU Calib Inverse dispersion Y\n\n"
00085 "Output files:\n\n"
00086 " DO category: Data type: Explanation:\n"
00087 " REDUCED_SCI_MXU FITS image Extracted scientific spectra\n"
00088 " REDUCED_SKY_SCI_MXU FITS image Extracted sky spectra\n"
00089 " REDUCED_ERROR_SCI_MXU FITS image Errors on extracted spectra\n"
00090 " UNMAPPED_SCI_MXU FITS image Sky subtracted scientific spectra\n"
00091 " MAPPED_SCI_MXU FITS image Rectified scientific spectra\n"
00092 " MAPPED_ALL_SCI_MXU FITS image Rectified science spectra with sky\n"
00093 " MAPPED_SKY_SCI_MXU FITS image Rectified sky spectra\n"
00094 " UNMAPPED_SKY_SCI_MXU FITS image Sky on CCD\n"
00095 " GLOBAL_SKY_SPECTRUM_MXU FITS table Global sky spectrum\n"
00096 " OBJECT_TABLE_SCI_MXU FITS table Positions of detected objects\n"
00097 "\n"
00098 " Only if the sky-alignment of the wavelength solution is requested:\n"
00099 " SKY_SHIFTS_LONG_SCI_MXU FITS table Sky lines offsets (LSS-like data)\n"
00100 " or SKY_SHIFTS_SLIT_SCI_MXU FITS table Sky lines offsets (MOS-like data)\n"
00101 " DISP_COEFF_SCI_MXU FITS table Upgraded dispersion coefficients\n"
00102 " WAVELENGTH_MAP_SCI_MXU FITS image Upgraded wavelength map\n\n";
00103
00104 #define fors_science_exit(message) \
00105 { \
00106 if (message) cpl_msg_error(recipe, message); \
00107 cpl_free(exptime); \
00108 cpl_free(instrume); \
00109 cpl_image_delete(dummy); \
00110 cpl_image_delete(mapped); \
00111 cpl_image_delete(mapped_sky); \
00112 cpl_image_delete(mapped_cleaned); \
00113 cpl_image_delete(skylocalmap); \
00114 cpl_image_delete(skymap); \
00115 cpl_image_delete(smapped); \
00116 cpl_table_delete(offsets); \
00117 cpl_table_delete(sky); \
00118 cpl_image_delete(bias); \
00119 cpl_image_delete(spectra); \
00120 cpl_image_delete(coordinate); \
00121 cpl_image_delete(norm_flat); \
00122 cpl_image_delete(rainbow); \
00123 cpl_image_delete(rectified); \
00124 cpl_image_delete(wavemap); \
00125 cpl_image_delete(wavemaplss); \
00126 cpl_propertylist_delete(header); \
00127 cpl_propertylist_delete(save_header); \
00128 cpl_table_delete(grism_table); \
00129 cpl_table_delete(idscoeff); \
00130 cpl_table_delete(maskslits); \
00131 cpl_table_delete(overscans); \
00132 cpl_table_delete(polytraces); \
00133 cpl_table_delete(slits); \
00134 cpl_table_delete(wavelengths); \
00135 cpl_vector_delete(lines); \
00136 cpl_msg_indent_less(); \
00137 return -1; \
00138 }
00139
00140
00141 #define fors_science_exit_memcheck(message) \
00142 { \
00143 if (message) cpl_msg_info(recipe, message); \
00144 cpl_free(exptime); \
00145 cpl_free(instrume); \
00146 cpl_image_delete(dummy); \
00147 cpl_image_delete(mapped); \
00148 cpl_image_delete(mapped_cleaned); \
00149 cpl_image_delete(mapped_sky); \
00150 cpl_image_delete(skylocalmap); \
00151 cpl_image_delete(skymap); \
00152 cpl_image_delete(smapped); \
00153 cpl_table_delete(offsets); \
00154 cpl_table_delete(sky); \
00155 cpl_image_delete(bias); \
00156 cpl_image_delete(spectra); \
00157 cpl_image_delete(coordinate); \
00158 cpl_image_delete(norm_flat); \
00159 cpl_image_delete(rainbow); \
00160 cpl_image_delete(rectified); \
00161 cpl_image_delete(wavemap); \
00162 cpl_image_delete(wavemaplss); \
00163 cpl_propertylist_delete(header); \
00164 cpl_propertylist_delete(save_header); \
00165 cpl_table_delete(grism_table); \
00166 cpl_table_delete(idscoeff); \
00167 cpl_table_delete(maskslits); \
00168 cpl_table_delete(overscans); \
00169 cpl_table_delete(polytraces); \
00170 cpl_table_delete(slits); \
00171 cpl_table_delete(wavelengths); \
00172 cpl_vector_delete(lines); \
00173 cpl_msg_indent_less(); \
00174 return 0; \
00175 }
00176
00177
00189 int cpl_plugin_get_info(cpl_pluginlist *list)
00190 {
00191 cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe );
00192 cpl_plugin *plugin = &recipe->interface;
00193
00194 cpl_plugin_init(plugin,
00195 CPL_PLUGIN_API,
00196 FORS_BINARY_VERSION,
00197 CPL_PLUGIN_TYPE_RECIPE,
00198 "fors_science",
00199 "Extraction of scientific spectra",
00200 fors_science_description,
00201 "Carlo Izzo",
00202 PACKAGE_BUGREPORT,
00203 "This file is currently part of the FORS Instrument Pipeline\n"
00204 "Copyright (C) 2002-2006 European Southern Observatory\n\n"
00205 "This program is free software; you can redistribute it and/or modify\n"
00206 "it under the terms of the GNU General Public License as published by\n"
00207 "the Free Software Foundation; either version 2 of the License, or\n"
00208 "(at your option) any later version.\n\n"
00209 "This program is distributed in the hope that it will be useful,\n"
00210 "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
00211 "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
00212 "GNU General Public License for more details.\n\n"
00213 "You should have received a copy of the GNU General Public License\n"
00214 "along with this program; if not, write to the Free Software Foundation,\n"
00215 "Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA\n",
00216 fors_science_create,
00217 fors_science_exec,
00218 fors_science_destroy);
00219
00220 cpl_pluginlist_append(list, plugin);
00221
00222 return 0;
00223 }
00224
00225
00236 static int fors_science_create(cpl_plugin *plugin)
00237 {
00238 cpl_recipe *recipe;
00239 cpl_parameter *p;
00240
00241
00242
00243
00244
00245
00246 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00247 recipe = (cpl_recipe *)plugin;
00248 else
00249 return -1;
00250
00251
00252
00253
00254
00255 recipe->parameters = cpl_parameterlist_new();
00256
00257
00258
00259
00260
00261
00262 p = cpl_parameter_new_value("fors.fors_science.dispersion",
00263 CPL_TYPE_DOUBLE,
00264 "Resampling step (Angstrom/pixel)",
00265 "fors.fors_science",
00266 0.0);
00267 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dispersion");
00268 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00269 cpl_parameterlist_append(recipe->parameters, p);
00270
00271
00272
00273
00274
00275 p = cpl_parameter_new_value("fors.fors_science.skyalign",
00276 CPL_TYPE_INT,
00277 "Polynomial order for sky lines alignment, "
00278 "or -1 to avoid alignment",
00279 "fors.fors_science",
00280 0);
00281 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skyalign");
00282 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00283 cpl_parameterlist_append(recipe->parameters, p);
00284
00285
00286
00287
00288
00289 p = cpl_parameter_new_value("fors.fors_science.wcolumn",
00290 CPL_TYPE_STRING,
00291 "Name of sky line catalog table column "
00292 "with wavelengths",
00293 "fors.fors_science",
00294 "WLEN");
00295 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wcolumn");
00296 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00297 cpl_parameterlist_append(recipe->parameters, p);
00298
00299
00300
00301
00302
00303 p = cpl_parameter_new_value("fors.fors_science.startwavelength",
00304 CPL_TYPE_DOUBLE,
00305 "Start wavelength in spectral extraction",
00306 "fors.fors_science",
00307 0.0);
00308 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "startwavelength");
00309 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00310 cpl_parameterlist_append(recipe->parameters, p);
00311
00312
00313
00314
00315
00316 p = cpl_parameter_new_value("fors.fors_science.endwavelength",
00317 CPL_TYPE_DOUBLE,
00318 "End wavelength in spectral extraction",
00319 "fors.fors_science",
00320 0.0);
00321 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "endwavelength");
00322 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00323 cpl_parameterlist_append(recipe->parameters, p);
00324
00325
00326
00327
00328
00329 p = cpl_parameter_new_value("fors.fors_science.flux",
00330 CPL_TYPE_BOOL,
00331 "Apply flux conservation",
00332 "fors.fors_science",
00333 TRUE);
00334 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "flux");
00335 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00336 cpl_parameterlist_append(recipe->parameters, p);
00337
00338
00339
00340
00341
00342 p = cpl_parameter_new_value("fors.fors_science.flatfield",
00343 CPL_TYPE_BOOL,
00344 "Apply flat field",
00345 "fors.fors_science",
00346 TRUE);
00347 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "flatfield");
00348 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00349 cpl_parameterlist_append(recipe->parameters, p);
00350
00351
00352
00353
00354
00355 p = cpl_parameter_new_value("fors.fors_science.skyglobal",
00356 CPL_TYPE_BOOL,
00357 "Subtract global sky spectrum from CCD",
00358 "fors.fors_science",
00359 FALSE);
00360 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skyglobal");
00361 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00362 cpl_parameterlist_append(recipe->parameters, p);
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381 p = cpl_parameter_new_value("fors.fors_science.skymedian",
00382 CPL_TYPE_BOOL,
00383 "Sky subtraction from extracted slit spectra",
00384 "fors.fors_science",
00385 FALSE);
00386 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skymedian");
00387 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00388 cpl_parameterlist_append(recipe->parameters, p);
00389
00390
00391
00392
00393
00394 p = cpl_parameter_new_value("fors.fors_science.skylocal",
00395 CPL_TYPE_BOOL,
00396 "Sky subtraction from CCD slit spectra",
00397 "fors.fors_science",
00398 TRUE);
00399 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skylocal");
00400 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00401 cpl_parameterlist_append(recipe->parameters, p);
00402
00403
00404
00405
00406
00407 p = cpl_parameter_new_value("fors.fors_science.cosmics",
00408 CPL_TYPE_BOOL,
00409 "Eliminate cosmic rays hits (only if global "
00410 "sky subtraction is also requested)",
00411 "fors.fors_science",
00412 FALSE);
00413 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cosmics");
00414 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00415 cpl_parameterlist_append(recipe->parameters, p);
00416
00417
00418
00419
00420
00421 p = cpl_parameter_new_value("fors.fors_science.slit_margin",
00422 CPL_TYPE_INT,
00423 "Number of pixels to exclude at each slit "
00424 "in object detection and extraction",
00425 "fors.fors_science",
00426 3);
00427 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "slit_margin");
00428 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00429 cpl_parameterlist_append(recipe->parameters, p);
00430
00431
00432
00433
00434
00435 p = cpl_parameter_new_value("fors.fors_science.ext_radius",
00436 CPL_TYPE_INT,
00437 "Maximum extraction radius for detected "
00438 "objects (pixel)",
00439 "fors.fors_science",
00440 6);
00441 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ext_radius");
00442 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00443 cpl_parameterlist_append(recipe->parameters, p);
00444
00445
00446
00447
00448
00449 p = cpl_parameter_new_value("fors.fors_science.cont_radius",
00450 CPL_TYPE_INT,
00451 "Minimum distance at which two objects "
00452 "of equal luminosity do not contaminate "
00453 "each other (pixel)",
00454 "fors.fors_science",
00455 0);
00456 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cont_radius");
00457 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00458 cpl_parameterlist_append(recipe->parameters, p);
00459
00460
00461
00462
00463
00464 p = cpl_parameter_new_value("fors.fors_science.ext_mode",
00465 CPL_TYPE_INT,
00466 "Object extraction method: 0 = aperture, "
00467 "1 = Horne optimal extraction",
00468 "fors.fors_science",
00469 1);
00470 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ext_mode");
00471 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00472 cpl_parameterlist_append(recipe->parameters, p);
00473
00474
00475
00476
00477
00478 p = cpl_parameter_new_value("fors.fors_science.time_normalise",
00479 CPL_TYPE_BOOL,
00480 "Normalise output spectra by the exposure time",
00481 "fors.fors_science",
00482 TRUE);
00483 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "time_normalise");
00484 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00485 cpl_parameterlist_append(recipe->parameters, p);
00486
00487 return 0;
00488 }
00489
00490
00499 static int fors_science_exec(cpl_plugin *plugin)
00500 {
00501 cpl_recipe *recipe;
00502
00503 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00504 recipe = (cpl_recipe *)plugin;
00505 else
00506 return -1;
00507
00508 return fors_science(recipe->parameters, recipe->frames);
00509 }
00510
00511
00520 static int fors_science_destroy(cpl_plugin *plugin)
00521 {
00522 cpl_recipe *recipe;
00523
00524 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
00525 recipe = (cpl_recipe *)plugin;
00526 else
00527 return -1;
00528
00529 cpl_parameterlist_delete(recipe->parameters);
00530
00531 return 0;
00532 }
00533
00534
00544 static int fors_science(cpl_parameterlist *parlist, cpl_frameset *frameset)
00545 {
00546
00547 const char *recipe = "fors_science";
00548
00549
00550
00551
00552
00553
00554 double dispersion;
00555 int skyalign;
00556 const char *wcolumn;
00557 double startwavelength;
00558 double endwavelength;
00559 int flux;
00560 int flatfield;
00561 int skyglobal;
00562 int skylocal;
00563 int skymedian;
00564 int cosmics;
00565 int slit_margin;
00566 int ext_radius;
00567 int cont_radius;
00568 int ext_mode;
00569 int time_normalise;
00570
00571
00572
00573
00574
00575 cpl_imagelist *all_science;
00576 cpl_image **images;
00577
00578 cpl_image *bias = NULL;
00579 cpl_image *norm_flat = NULL;
00580 cpl_image *spectra = NULL;
00581 cpl_image *rectified = NULL;
00582 cpl_image *coordinate = NULL;
00583 cpl_image *rainbow = NULL;
00584 cpl_image *mapped = NULL;
00585 cpl_image *mapped_sky = NULL;
00586 cpl_image *mapped_cleaned = NULL;
00587 cpl_image *smapped = NULL;
00588 cpl_image *wavemap = NULL;
00589 cpl_image *wavemaplss = NULL;
00590 cpl_image *skymap = NULL;
00591 cpl_image *skylocalmap = NULL;
00592 cpl_image *dummy = NULL;
00593
00594 cpl_table *grism_table = NULL;
00595 cpl_table *overscans = NULL;
00596 cpl_table *wavelengths = NULL;
00597 cpl_table *idscoeff = NULL;
00598 cpl_table *slits = NULL;
00599 cpl_table *maskslits = NULL;
00600 cpl_table *polytraces = NULL;
00601 cpl_table *offsets = NULL;
00602 cpl_table *sky = NULL;
00603
00604 cpl_vector *lines = NULL;
00605
00606 cpl_propertylist *header = NULL;
00607 cpl_propertylist *save_header = NULL;
00608
00609
00610
00611
00612
00613 char version[80];
00614 char *instrume = NULL;
00615 const char *science_tag;
00616 const char *master_norm_flat_tag;
00617 const char *disp_coeff_tag;
00618 const char *disp_coeff_sky_tag;
00619 const char *wavelength_map_sky_tag;
00620 const char *curv_coeff_tag;
00621 const char *slit_location_tag;
00622 const char *reduced_science_tag;
00623 const char *reduced_sky_tag;
00624 const char *reduced_error_tag;
00625 const char *mapped_science_tag;
00626 const char *unmapped_science_tag;
00627 const char *mapped_science_sky_tag;
00628 const char *mapped_sky_tag;
00629 const char *unmapped_sky_tag;
00630 const char *global_sky_spectrum_tag;
00631 const char *object_table_tag;
00632 const char *skylines_offsets_tag;
00633 int mxu, mos, lss;
00634 int treat_as_lss = 0;
00635 int nslits;
00636 int nscience;
00637 double *xpos;
00638 double *exptime = NULL;
00639 double alltime;
00640 double mxpos;
00641 double mean_rms;
00642 int nlines;
00643 int rebin;
00644 double *line;
00645 int nx, ny;
00646 int ccd_xsize, ccd_ysize;
00647 double reference;
00648 double gain;
00649 double ron;
00650 int standard;
00651 int highres;
00652 int i;
00653
00654
00655 snprintf(version, 80, "%s-%s", PACKAGE, PACKAGE_VERSION);
00656
00657 cpl_msg_set_indentation(2);
00658
00659 if (dfs_files_dont_exist(frameset))
00660 fors_science_exit(NULL);
00661
00662
00663
00664
00665
00666
00667 cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
00668 cpl_msg_indent_more();
00669
00670 if (cpl_frameset_count_tags(frameset, "GRISM_TABLE") > 1)
00671 fors_science_exit("Too many in input: GRISM_TABLE");
00672
00673 grism_table = dfs_load_table(frameset, "GRISM_TABLE", 1);
00674
00675 dispersion = dfs_get_parameter_double(parlist,
00676 "fors.fors_science.dispersion", grism_table);
00677
00678 if (dispersion <= 0.0)
00679 fors_science_exit("Invalid resampling step");
00680
00681 skyalign = dfs_get_parameter_int(parlist,
00682 "fors.fors_science.skyalign", NULL);
00683
00684 if (skyalign > 2)
00685 fors_science_exit("Max polynomial degree for sky alignment is 2");
00686
00687 wcolumn = dfs_get_parameter_string(parlist,
00688 "fors.fors_science.wcolumn", NULL);
00689
00690 startwavelength = dfs_get_parameter_double(parlist,
00691 "fors.fors_science.startwavelength", grism_table);
00692 if (startwavelength < 3000.0 || startwavelength > 13000.0)
00693 fors_science_exit("Invalid wavelength");
00694
00695 endwavelength = dfs_get_parameter_double(parlist,
00696 "fors.fors_science.endwavelength", grism_table);
00697 if (endwavelength < 3000.0 || endwavelength > 13000.0)
00698 fors_science_exit("Invalid wavelength");
00699
00700 if (endwavelength - startwavelength <= 0.0)
00701 fors_science_exit("Invalid wavelength interval");
00702
00703 flux = dfs_get_parameter_bool(parlist, "fors.fors_science.flux", NULL);
00704
00705 flatfield = dfs_get_parameter_bool(parlist, "fors.fors_science.flatfield",
00706 NULL);
00707
00708 skyglobal = dfs_get_parameter_bool(parlist, "fors.fors_science.skyglobal",
00709 NULL);
00710 skylocal = dfs_get_parameter_bool(parlist, "fors.fors_science.skylocal",
00711 NULL);
00712 skymedian = dfs_get_parameter_bool(parlist, "fors.fors_science.skymedian",
00713 NULL);
00714
00715
00716
00717
00718
00719 if (skylocal && skyglobal)
00720 fors_science_exit("Cannot apply both local and global sky subtraction");
00721
00722 if (skylocal && skymedian)
00723 fors_science_exit("Cannot apply sky subtraction both on extracted "
00724 "and non-extracted spectra");
00725
00726 cosmics = dfs_get_parameter_bool(parlist,
00727 "fors.fors_science.cosmics", NULL);
00728
00729 if (cosmics)
00730 if (!(skyglobal || skylocal))
00731 fors_science_exit("Cosmic rays correction requires "
00732 "either skylocal=true or skyglobal=true");
00733
00734 slit_margin = dfs_get_parameter_int(parlist,
00735 "fors.fors_science.slit_margin",
00736 NULL);
00737 if (slit_margin < 0)
00738 fors_science_exit("Value must be zero or positive");
00739
00740 ext_radius = dfs_get_parameter_int(parlist, "fors.fors_science.ext_radius",
00741 NULL);
00742 if (ext_radius < 0)
00743 fors_science_exit("Value must be zero or positive");
00744
00745 cont_radius = dfs_get_parameter_int(parlist,
00746 "fors.fors_science.cont_radius",
00747 NULL);
00748 if (cont_radius < 0)
00749 fors_science_exit("Value must be zero or positive");
00750
00751 ext_mode = dfs_get_parameter_int(parlist, "fors.fors_science.ext_mode",
00752 NULL);
00753 if (ext_mode < 0 || ext_mode > 1)
00754 fors_science_exit("Invalid object extraction mode");
00755
00756 time_normalise = dfs_get_parameter_bool(parlist,
00757 "fors.fors_science.time_normalise", NULL);
00758
00759 cpl_table_delete(grism_table); grism_table = NULL;
00760
00761 if (cpl_error_get_code())
00762 fors_science_exit("Failure getting the configuration parameters");
00763
00764
00765
00766
00767
00768
00769 cpl_msg_indent_less();
00770 cpl_msg_info(recipe, "Check input set-of-frames:");
00771 cpl_msg_indent_more();
00772
00773 if (!dfs_equal_keyword(frameset, "ESO INS GRIS1 ID"))
00774 fors_science_exit("Input frames are not from the same grism");
00775
00776 if (!dfs_equal_keyword(frameset, "ESO INS FILT1 ID"))
00777 fors_science_exit("Input frames are not from the same filter");
00778
00779 if (!dfs_equal_keyword(frameset, "ESO DET CHIP1 ID"))
00780 fors_science_exit("Input frames are not from the same chip");
00781
00782 mxu = cpl_frameset_count_tags(frameset, "SCIENCE_MXU");
00783 mos = cpl_frameset_count_tags(frameset, "SCIENCE_MOS");
00784 lss = cpl_frameset_count_tags(frameset, "SCIENCE_LSS");
00785 standard = 0;
00786
00787 if (mxu + mos + lss == 0) {
00788 mxu = cpl_frameset_count_tags(frameset, "STANDARD_MXU");
00789 mos = cpl_frameset_count_tags(frameset, "STANDARD_MOS");
00790 lss = cpl_frameset_count_tags(frameset, "STANDARD_LSS");
00791 standard = 1;
00792 }
00793
00794 if (mxu + mos + lss == 0)
00795 fors_science_exit("Missing input scientific frame");
00796
00797 nscience = mxu + mos + lss;
00798
00799 if (mxu && mxu < nscience)
00800 fors_science_exit("Input scientific frames must be of the same type");
00801
00802 if (mos && mos < nscience)
00803 fors_science_exit("Input scientific frames must be of the same type");
00804
00805 if (lss && lss < nscience)
00806 fors_science_exit("Input scientific frames must be of the same type");
00807
00808 if (mxu) {
00809 if (standard) {
00810 cpl_msg_info(recipe, "MXU data found");
00811 science_tag = "STANDARD_MXU";
00812 reduced_science_tag = "REDUCED_STD_MXU";
00813 unmapped_science_tag = "UNMAPPED_STD_MXU";
00814 mapped_science_tag = "MAPPED_STD_MXU";
00815 mapped_science_sky_tag = "MAPPED_ALL_STD_MXU";
00816 skylines_offsets_tag = "SKY_SHIFTS_SLIT_STD_MXU";
00817 wavelength_map_sky_tag = "WAVELENGTH_MAP_STD_MXU";
00818 disp_coeff_sky_tag = "DISP_COEFF_STD_MXU";
00819 mapped_sky_tag = "MAPPED_SKY_STD_MXU";
00820 unmapped_sky_tag = "UNMAPPED_SKY_STD_MXU";
00821 object_table_tag = "OBJECT_TABLE_STD_MXU";
00822 reduced_sky_tag = "REDUCED_SKY_STD_MXU";
00823 reduced_error_tag = "REDUCED_ERROR_STD_MXU";
00824 }
00825 else {
00826 cpl_msg_info(recipe, "MXU data found");
00827 science_tag = "SCIENCE_MXU";
00828 reduced_science_tag = "REDUCED_SCI_MXU";
00829 unmapped_science_tag = "UNMAPPED_SCI_MXU";
00830 mapped_science_tag = "MAPPED_SCI_MXU";
00831 mapped_science_sky_tag = "MAPPED_ALL_SCI_MXU";
00832 skylines_offsets_tag = "SKY_SHIFTS_SLIT_SCI_MXU";
00833 wavelength_map_sky_tag = "WAVELENGTH_MAP_SCI_MXU";
00834 disp_coeff_sky_tag = "DISP_COEFF_SCI_MXU";
00835 mapped_sky_tag = "MAPPED_SKY_SCI_MXU";
00836 unmapped_sky_tag = "UNMAPPED_SKY_SCI_MXU";
00837 object_table_tag = "OBJECT_TABLE_SCI_MXU";
00838 reduced_sky_tag = "REDUCED_SKY_SCI_MXU";
00839 reduced_error_tag = "REDUCED_ERROR_SCI_MXU";
00840 }
00841
00842 master_norm_flat_tag = "MASTER_NORM_FLAT_MXU";
00843 disp_coeff_tag = "DISP_COEFF_MXU";
00844 curv_coeff_tag = "CURV_COEFF_MXU";
00845 slit_location_tag = "SLIT_LOCATION_MXU";
00846 global_sky_spectrum_tag = "GLOBAL_SKY_SPECTRUM_MXU";
00847
00848 if (!cpl_frameset_count_tags(frameset, master_norm_flat_tag)) {
00849 master_norm_flat_tag = "MASTER_NORM_FLAT_LONG_MXU";
00850 disp_coeff_tag = "DISP_COEFF_LONG_MXU";
00851 slit_location_tag = "SLIT_LOCATION_LONG_MXU";
00852 }
00853 }
00854
00855 if (lss) {
00856
00857 if (cosmics && !skyglobal)
00858 fors_science_exit("Cosmic rays correction for LSS "
00859 "data requires --skyglobal=true");
00860
00861 if (standard) {
00862 cpl_msg_info(recipe, "LSS data found");
00863 science_tag = "STANDARD_LSS";
00864 reduced_science_tag = "REDUCED_STD_LSS";
00865 unmapped_science_tag = "UNMAPPED_STD_LSS";
00866 mapped_science_tag = "MAPPED_STD_LSS";
00867 mapped_science_sky_tag = "MAPPED_ALL_STD_LSS";
00868 skylines_offsets_tag = "SKY_SHIFTS_LONG_STD_LSS";
00869 wavelength_map_sky_tag = "WAVELENGTH_MAP_STD_LSS";
00870 disp_coeff_sky_tag = "DISP_COEFF_STD_LSS";
00871 mapped_sky_tag = "MAPPED_SKY_STD_LSS";
00872 unmapped_sky_tag = "UNMAPPED_SKY_STD_LSS";
00873 object_table_tag = "OBJECT_TABLE_STD_LSS";
00874 reduced_sky_tag = "REDUCED_SKY_STD_LSS";
00875 reduced_error_tag = "REDUCED_ERROR_STD_LSS";
00876 }
00877 else {
00878 cpl_msg_info(recipe, "LSS data found");
00879 science_tag = "SCIENCE_LSS";
00880 reduced_science_tag = "REDUCED_SCI_LSS";
00881 unmapped_science_tag = "UNMAPPED_SCI_LSS";
00882 mapped_science_tag = "MAPPED_SCI_LSS";
00883 mapped_science_sky_tag = "MAPPED_ALL_SCI_LSS";
00884 skylines_offsets_tag = "SKY_SHIFTS_LONG_SCI_LSS";
00885 wavelength_map_sky_tag = "WAVELENGTH_MAP_SCI_LSS";
00886 disp_coeff_sky_tag = "DISP_COEFF_SCI_LSS";
00887 mapped_sky_tag = "MAPPED_SKY_SCI_LSS";
00888 unmapped_sky_tag = "UNMAPPED_SKY_SCI_LSS";
00889 object_table_tag = "OBJECT_TABLE_SCI_LSS";
00890 reduced_sky_tag = "REDUCED_SKY_SCI_LSS";
00891 reduced_error_tag = "REDUCED_ERROR_SCI_LSS";
00892 }
00893
00894 master_norm_flat_tag = "MASTER_NORM_FLAT_LSS";
00895 disp_coeff_tag = "DISP_COEFF_LSS";
00896 slit_location_tag = "SLIT_LOCATION_LSS";
00897 global_sky_spectrum_tag = "GLOBAL_SKY_SPECTRUM_LSS";
00898 }
00899
00900 if (mos) {
00901 cpl_msg_info(recipe, "MOS data found");
00902 if (standard) {
00903 science_tag = "STANDARD_MOS";
00904 reduced_science_tag = "REDUCED_STD_MOS";
00905 unmapped_science_tag = "UNMAPPED_STD_MOS";
00906 mapped_science_tag = "MAPPED_STD_MOS";
00907 mapped_science_sky_tag = "MAPPED_ALL_STD_MOS";
00908 skylines_offsets_tag = "SKY_SHIFTS_SLIT_STD_MOS";
00909 wavelength_map_sky_tag = "WAVELENGTH_MAP_STD_MOS";
00910 disp_coeff_sky_tag = "DISP_COEFF_STD_MOS";
00911 mapped_sky_tag = "MAPPED_SKY_STD_MOS";
00912 unmapped_sky_tag = "UNMAPPED_SKY_STD_MOS";
00913 object_table_tag = "OBJECT_TABLE_STD_MOS";
00914 reduced_sky_tag = "REDUCED_SKY_STD_MOS";
00915 reduced_error_tag = "REDUCED_ERROR_STD_MOS";
00916 }
00917 else {
00918 science_tag = "SCIENCE_MOS";
00919 reduced_science_tag = "REDUCED_SCI_MOS";
00920 unmapped_science_tag = "UNMAPPED_SCI_MOS";
00921 mapped_science_tag = "MAPPED_SCI_MOS";
00922 mapped_science_sky_tag = "MAPPED_ALL_SCI_MOS";
00923 skylines_offsets_tag = "SKY_SHIFTS_SLIT_SCI_MOS";
00924 wavelength_map_sky_tag = "WAVELENGTH_MAP_SCI_MOS";
00925 disp_coeff_sky_tag = "DISP_COEFF_SCI_MOS";
00926 mapped_sky_tag = "MAPPED_SKY_SCI_MOS";
00927 unmapped_sky_tag = "UNMAPPED_SKY_SCI_MOS";
00928 object_table_tag = "OBJECT_TABLE_SCI_MOS";
00929 reduced_sky_tag = "REDUCED_SKY_SCI_MOS";
00930 reduced_error_tag = "REDUCED_ERROR_SCI_MOS";
00931 }
00932
00933 master_norm_flat_tag = "MASTER_NORM_FLAT_MOS";
00934 disp_coeff_tag = "DISP_COEFF_MOS";
00935 curv_coeff_tag = "CURV_COEFF_MOS";
00936 slit_location_tag = "SLIT_LOCATION_MOS";
00937 global_sky_spectrum_tag = "GLOBAL_SKY_SPECTRUM_MOS";
00938
00939 if (!cpl_frameset_count_tags(frameset, master_norm_flat_tag)) {
00940 master_norm_flat_tag = "MASTER_NORM_FLAT_LONG_MOS";
00941 disp_coeff_tag = "DISP_COEFF_LONG_MOS";
00942 slit_location_tag = "SLIT_LOCATION_LONG_MOS";
00943 }
00944 }
00945
00946 if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") == 0)
00947 fors_science_exit("Missing required input: MASTER_BIAS");
00948
00949 if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") > 1)
00950 fors_science_exit("Too many in input: MASTER_BIAS");
00951
00952 if (skyalign >= 0)
00953 if (cpl_frameset_count_tags(frameset, "MASTER_SKYLINECAT") > 1)
00954 fors_science_exit("Too many in input: MASTER_SKYLINECAT");
00955
00956 if (cpl_frameset_count_tags(frameset, disp_coeff_tag) == 0) {
00957 cpl_msg_error(recipe, "Missing required input: %s", disp_coeff_tag);
00958 fors_science_exit(NULL);
00959 }
00960
00961 if (cpl_frameset_count_tags(frameset, disp_coeff_tag) > 1) {
00962 cpl_msg_error(recipe, "Too many in input: %s", disp_coeff_tag);
00963 fors_science_exit(NULL);
00964 }
00965
00966 if (cpl_frameset_count_tags(frameset, slit_location_tag) == 0) {
00967 cpl_msg_error(recipe, "Missing required input: %s",
00968 slit_location_tag);
00969 fors_science_exit(NULL);
00970 }
00971
00972 if (cpl_frameset_count_tags(frameset, slit_location_tag) > 1) {
00973 cpl_msg_error(recipe, "Too many in input: %s", slit_location_tag);
00974 fors_science_exit(NULL);
00975 }
00976
00977 if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) > 1) {
00978 if (flatfield) {
00979 cpl_msg_error(recipe, "Too many in input: %s",
00980 master_norm_flat_tag);
00981 fors_science_exit(NULL);
00982 }
00983 else {
00984 cpl_msg_warning(recipe, "%s in input are ignored, "
00985 "since flat field correction was not requested",
00986 master_norm_flat_tag);
00987 }
00988 }
00989
00990 if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) == 1) {
00991 if (!flatfield) {
00992 cpl_msg_warning(recipe, "%s in input is ignored, "
00993 "since flat field correction was not requested",
00994 master_norm_flat_tag);
00995 }
00996 }
00997
00998 if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) == 0) {
00999 if (flatfield) {
01000 cpl_msg_error(recipe, "Flat field correction was requested, "
01001 "but no %s are found in input",
01002 master_norm_flat_tag);
01003 fors_science_exit(NULL);
01004 }
01005 }
01006
01007 cpl_msg_indent_less();
01008
01009
01010
01011
01012
01013
01014 exptime = cpl_calloc(nscience, sizeof(double));
01015
01016 if (nscience > 1) {
01017
01018 cpl_msg_info(recipe, "Load %d scientific frames and median them...",
01019 nscience);
01020 cpl_msg_indent_more();
01021
01022 all_science = cpl_imagelist_new();
01023
01024 header = dfs_load_header(frameset, science_tag, 0);
01025
01026 if (header == NULL)
01027 fors_science_exit("Cannot load scientific frame header");
01028
01029 alltime = exptime[0] = cpl_propertylist_get_double(header, "EXPTIME");
01030
01031 if (cpl_error_get_code() != CPL_ERROR_NONE)
01032 fors_science_exit("Missing keyword EXPTIME in scientific "
01033 "frame header");
01034
01035 cpl_propertylist_delete(header); header = NULL;
01036
01037 cpl_msg_info(recipe, "Scientific frame 1 exposure time: %.2f s",
01038 exptime[0]);
01039
01040 for (i = 1; i < nscience; i++) {
01041
01042 header = dfs_load_header(frameset, NULL, 0);
01043
01044 if (header == NULL)
01045 fors_science_exit("Cannot load scientific frame header");
01046
01047 exptime[i] = cpl_propertylist_get_double(header, "EXPTIME");
01048
01049 alltime += exptime[i];
01050
01051 if (cpl_error_get_code() != CPL_ERROR_NONE)
01052 fors_science_exit("Missing keyword EXPTIME in scientific "
01053 "frame header");
01054
01055 cpl_propertylist_delete(header); header = NULL;
01056
01057 cpl_msg_info(recipe, "Scientific frame %d exposure time: %.2f s",
01058 i+1, exptime[i]);
01059 }
01060
01061 spectra = dfs_load_image(frameset, science_tag, CPL_TYPE_FLOAT, 0, 0);
01062
01063 if (spectra == NULL)
01064 fors_science_exit("Cannot load scientific frame");
01065
01066 cpl_image_divide_scalar(spectra, exptime[0]);
01067 cpl_imagelist_set(all_science, spectra, 0);
01068
01069 for (i = 1; i < nscience; i++) {
01070
01071 spectra = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
01072
01073 if (spectra) {
01074 cpl_image_divide_scalar(spectra, exptime[i]);
01075 cpl_imagelist_set(all_science, spectra, i);
01076 }
01077 else
01078 fors_science_exit("Cannot load scientific frame");
01079
01080 }
01081
01082 spectra = cpl_imagelist_collapse_median_create(all_science);
01083 cpl_image_multiply_scalar(spectra, alltime);
01084
01085 cpl_imagelist_delete(all_science);
01086 }
01087 else {
01088 cpl_msg_info(recipe, "Load scientific exposure...");
01089 cpl_msg_indent_more();
01090
01091 header = dfs_load_header(frameset, science_tag, 0);
01092
01093 if (header == NULL)
01094 fors_science_exit("Cannot load scientific frame header");
01095
01096 alltime = exptime[0] = cpl_propertylist_get_double(header, "EXPTIME");
01097
01098 if (cpl_error_get_code() != CPL_ERROR_NONE)
01099 fors_science_exit("Missing keyword EXPTIME in scientific "
01100 "frame header");
01101
01102 cpl_propertylist_delete(header); header = NULL;
01103
01104 cpl_msg_info(recipe, "Scientific frame exposure time: %.2f s",
01105 exptime[0]);
01106
01107 spectra = dfs_load_image(frameset, science_tag, CPL_TYPE_FLOAT, 0, 0);
01108 }
01109
01110 if (spectra == NULL)
01111 fors_science_exit("Cannot load scientific frame");
01112
01113 cpl_free(exptime); exptime = NULL;
01114
01115 cpl_msg_indent_less();
01116
01117
01118
01119
01120
01121
01122
01123 header = dfs_load_header(frameset, science_tag, 0);
01124
01125 if (header == NULL)
01126 fors_science_exit("Cannot load scientific frame header");
01127
01128 instrume = (char *)cpl_propertylist_get_string(header, "INSTRUME");
01129 if (instrume == NULL)
01130 fors_science_exit("Missing keyword INSTRUME in scientific header");
01131 instrume = cpl_strdup(instrume);
01132
01133 if (instrume[4] == '1')
01134 snprintf(version, 80, "%s/%s", "fors1", VERSION);
01135 if (instrume[4] == '2')
01136 snprintf(version, 80, "%s/%s", "fors2", VERSION);
01137
01138 cpl_free(instrume); instrume = NULL;
01139
01140 reference = cpl_propertylist_get_double(header, "ESO INS GRIS1 WLEN");
01141
01142 if (cpl_error_get_code() != CPL_ERROR_NONE)
01143 fors_science_exit("Missing keyword ESO INS GRIS1 WLEN in scientific "
01144 "frame header");
01145
01146 if (reference < 3000.0)
01147 reference *= 10;
01148
01149 if (reference < 3000.0 || reference > 13000.0) {
01150 cpl_msg_error(recipe, "Invalid central wavelength %.2f read from "
01151 "keyword ESO INS GRIS1 WLEN in scientific frame header",
01152 reference);
01153 fors_science_exit(NULL);
01154 }
01155
01156 cpl_msg_info(recipe, "The central wavelength is: %.2f", reference);
01157
01158 rebin = cpl_propertylist_get_int(header, "ESO DET WIN1 BINX");
01159
01160 if (cpl_error_get_code() != CPL_ERROR_NONE)
01161 fors_science_exit("Missing keyword ESO DET WIN1 BINX in scientific "
01162 "frame header");
01163
01164 if (rebin != 1) {
01165 dispersion *= rebin;
01166 cpl_msg_warning(recipe, "The rebin factor is %d, and therefore the "
01167 "resampling step used is %f A/pixel", rebin,
01168 dispersion);
01169 }
01170
01171 gain = cpl_propertylist_get_double(header, "ESO DET OUT1 CONAD");
01172
01173 if (cpl_error_get_code() != CPL_ERROR_NONE)
01174 fors_science_exit("Missing keyword ESO DET OUT1 CONAD in scientific "
01175 "frame header");
01176
01177 cpl_msg_info(recipe, "The gain factor is: %.2f e-/ADU", gain);
01178
01179 ron = cpl_propertylist_get_double(header, "ESO DET OUT1 RON");
01180
01181 if (cpl_error_get_code() != CPL_ERROR_NONE)
01182 fors_science_exit("Missing keyword ESO DET OUT1 RON in scientific "
01183 "frame header");
01184
01185 ron /= gain;
01186
01187 cpl_msg_info(recipe, "The read-out-noise is: %.2f ADU", ron);
01188
01189 if (mos || mxu) {
01190
01191 if (mos)
01192 maskslits = mos_load_slits_fors_mos(header);
01193 else
01194 maskslits = mos_load_slits_fors_mxu(header);
01195
01196
01197
01198
01199
01200
01201 mxpos = cpl_table_get_column_median(maskslits, "xtop");
01202 xpos = cpl_table_get_data_double(maskslits, "xtop");
01203 nslits = cpl_table_get_nrow(maskslits);
01204
01205 treat_as_lss = 1;
01206 for (i = 0; i < nslits; i++) {
01207 if (fabs(mxpos-xpos[i]) > 0.01) {
01208 treat_as_lss = 0;
01209 break;
01210 }
01211 }
01212
01213 cpl_table_delete(maskslits); maskslits = NULL;
01214
01215
01216 if (treat_as_lss) {
01217 cpl_msg_warning(recipe, "All MOS slits have the same offset: %.2f\n"
01218 "The LSS data reduction strategy is applied!",
01219 mxpos);
01220 if (mos) {
01221 if (standard) {
01222 skylines_offsets_tag = "SKY_SHIFTS_LONG_STD_MOS";
01223 }
01224 else {
01225 skylines_offsets_tag = "SKY_SHIFTS_LONG_SCI_MOS";
01226 }
01227 }
01228 else {
01229 if (standard) {
01230 skylines_offsets_tag = "SKY_SHIFTS_LONG_STD_MXU";
01231 }
01232 else {
01233 skylines_offsets_tag = "SKY_SHIFTS_LONG_SCI_MXU";
01234 }
01235 }
01236 }
01237 }
01238
01239 if (lss || treat_as_lss) {
01240 if (skylocal) {
01241 if (cosmics)
01242 fors_science_exit("Cosmic rays correction for LSS or LSS-like "
01243 "data requires --skyglobal=true");
01244 skymedian = skylocal;
01245 skylocal = 0;
01246 }
01247 }
01248 else {
01249 if (cpl_frameset_count_tags(frameset, curv_coeff_tag) == 0) {
01250 cpl_msg_error(recipe, "Missing required input: %s", curv_coeff_tag);
01251 fors_science_exit(NULL);
01252 }
01253
01254 if (cpl_frameset_count_tags(frameset, curv_coeff_tag) > 1) {
01255 cpl_msg_error(recipe, "Too many in input: %s", curv_coeff_tag);
01256 fors_science_exit(NULL);
01257 }
01258 }
01259
01260
01261
01262
01263
01264
01265
01266
01267 cpl_msg_info(recipe, "Remove the master bias...");
01268
01269 bias = dfs_load_image(frameset, "MASTER_BIAS", CPL_TYPE_FLOAT, 0, 1);
01270
01271 if (bias == NULL)
01272 fors_science_exit("Cannot load master bias");
01273
01274 overscans = mos_load_overscans_vimos(header, 1);
01275 cpl_propertylist_delete(header); header = NULL;
01276 dummy = mos_remove_bias(spectra, bias, overscans);
01277 cpl_image_delete(spectra); spectra = dummy; dummy = NULL;
01278 cpl_image_delete(bias); bias = NULL;
01279 cpl_table_delete(overscans); overscans = NULL;
01280
01281 if (spectra == NULL)
01282 fors_science_exit("Cannot remove bias from scientific frame");
01283
01284 ccd_xsize = nx = cpl_image_get_size_x(spectra);
01285 ccd_ysize = ny = cpl_image_get_size_y(spectra);
01286
01287 cpl_msg_indent_less();
01288 cpl_msg_info(recipe, "Load normalised flat field (if present)...");
01289 cpl_msg_indent_more();
01290
01291 if (flatfield) {
01292
01293 norm_flat = dfs_load_image(frameset, master_norm_flat_tag,
01294 CPL_TYPE_FLOAT, 0, 1);
01295
01296 if (norm_flat) {
01297 cpl_msg_info(recipe, "Apply flat field correction...");
01298 if (cpl_image_divide(spectra, norm_flat) != CPL_ERROR_NONE) {
01299 cpl_msg_error(recipe, "Failure of flat field correction: %s",
01300 cpl_error_get_message());
01301 fors_science_exit(NULL);
01302 }
01303 cpl_image_delete(norm_flat); norm_flat = NULL;
01304 }
01305 else {
01306 cpl_msg_error(recipe, "Cannot load input %s for flat field "
01307 "correction", master_norm_flat_tag);
01308 fors_science_exit(NULL);
01309 }
01310
01311 }
01312
01313
01314 if (skyalign >= 0) {
01315 cpl_msg_indent_less();
01316 cpl_msg_info(recipe, "Load input sky line catalog...");
01317 cpl_msg_indent_more();
01318
01319 wavelengths = dfs_load_table(frameset, "MASTER_SKYLINECAT", 1);
01320
01321 if (wavelengths) {
01322
01323
01324
01325
01326
01327 nlines = cpl_table_get_nrow(wavelengths);
01328
01329 if (nlines == 0)
01330 fors_science_exit("Empty input sky line catalog");
01331
01332 if (cpl_table_has_column(wavelengths, wcolumn) != 1) {
01333 cpl_msg_error(recipe, "Missing column %s in input line "
01334 "catalog table", wcolumn);
01335 fors_science_exit(NULL);
01336 }
01337
01338 line = cpl_malloc(nlines * sizeof(double));
01339
01340 for (i = 0; i < nlines; i++)
01341 line[i] = cpl_table_get(wavelengths, wcolumn, i, NULL);
01342
01343 cpl_table_delete(wavelengths); wavelengths = NULL;
01344
01345 lines = cpl_vector_wrap(nlines, line);
01346 }
01347 else {
01348 cpl_msg_info(recipe, "No sky line catalog found in input - fine!");
01349 }
01350 }
01351
01352
01353
01354
01355
01356
01357
01358 if (!(lss || treat_as_lss)) {
01359 polytraces = dfs_load_table(frameset, curv_coeff_tag, 1);
01360 if (polytraces == NULL)
01361 fors_science_exit("Cannot load spectral curvature table");
01362 }
01363
01364
01365
01366
01367
01368
01369 slits = dfs_load_table(frameset, slit_location_tag, 1);
01370 if (slits == NULL)
01371 fors_science_exit("Cannot load slits location table");
01372
01373 if (lss || treat_as_lss) {
01374 int first_row = cpl_table_get_double(slits, "ybottom", 0, NULL);
01375 int last_row = cpl_table_get_double(slits, "ytop", 0, NULL);
01376 int ylow, yhig;
01377
01378 ylow = first_row + 1;
01379 yhig = last_row + 1;
01380
01381 dummy = cpl_image_extract(spectra, 1, ylow, nx, yhig);
01382 cpl_image_delete(spectra); spectra = dummy; dummy = NULL;
01383 ny = cpl_image_get_size_y(spectra);
01384 }
01385
01386
01387
01388
01389
01390
01391 idscoeff = dfs_load_table(frameset, disp_coeff_tag, 1);
01392 if (idscoeff == NULL)
01393 fors_science_exit("Cannot load wavelength calibration table");
01394
01395 cpl_msg_indent_less();
01396 cpl_msg_info(recipe, "Processing scientific spectra...");
01397 cpl_msg_indent_more();
01398
01399
01400
01401
01402
01403
01404 if (lss || treat_as_lss) {
01405 smapped = cpl_image_duplicate(spectra);
01406 }
01407 else {
01408 coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01409
01410 smapped = mos_spatial_calibration(spectra, slits, polytraces, reference,
01411 startwavelength, endwavelength,
01412 dispersion, flux, coordinate);
01413 }
01414
01415
01416
01417
01418
01419
01420
01421 rainbow = mos_map_idscoeff(idscoeff, nx, reference, startwavelength,
01422 endwavelength);
01423
01424 if (dispersion > 1.0)
01425 highres = 0;
01426 else
01427 highres = 1;
01428
01429 if (skyalign >= 0) {
01430 if (skyalign) {
01431 cpl_msg_info(recipe, "Align wavelength solution to reference "
01432 "skylines applying %d order residual fit...", skyalign);
01433 }
01434 else {
01435 cpl_msg_info(recipe, "Align wavelength solution to reference "
01436 "skylines applying median offset...");
01437 }
01438
01439 if (lss || treat_as_lss) {
01440 offsets = mos_wavelength_align_lss(smapped, reference,
01441 startwavelength, endwavelength,
01442 idscoeff, lines, highres,
01443 skyalign, rainbow, 4);
01444 }
01445 else {
01446 offsets = mos_wavelength_align(smapped, slits, reference,
01447 startwavelength, endwavelength,
01448 idscoeff, lines, highres, skyalign,
01449 rainbow, 4);
01450 }
01451
01452 if (offsets) {
01453 if (standard)
01454 cpl_msg_warning(recipe, "Alignment of the wavelength solution "
01455 "to reference sky lines may be unreliable in "
01456 "this case!");
01457
01458 if (dfs_save_table(frameset, offsets, skylines_offsets_tag, NULL,
01459 parlist, recipe, version))
01460 fors_science_exit(NULL);
01461
01462 cpl_table_delete(offsets); offsets = NULL;
01463 }
01464 else {
01465 cpl_msg_warning(recipe, "Alignment of the wavelength solution "
01466 "to reference sky lines could not be done!");
01467 skyalign = -1;
01468 }
01469
01470 }
01471
01472 if (lss || treat_as_lss) {
01473 int first_row = cpl_table_get_double(slits, "ybottom", 0, NULL);
01474 int last_row = cpl_table_get_double(slits, "ytop", 0, NULL);
01475 int ylow, yhig;
01476
01477 ylow = first_row + 1;
01478 yhig = last_row + 1;
01479
01480 wavemap = cpl_image_new(ccd_xsize, ccd_ysize, CPL_TYPE_FLOAT);
01481 cpl_image_copy(wavemap, rainbow, 1, ylow);
01482
01483 wavemaplss = cpl_image_extract(wavemap, 1, ylow, nx, yhig);
01484 }
01485 else {
01486 wavemap = mos_map_wavelengths(coordinate, rainbow, slits,
01487 polytraces, reference,
01488 startwavelength, endwavelength,
01489 dispersion);
01490 }
01491
01492 cpl_image_delete(rainbow); rainbow = NULL;
01493 cpl_image_delete(coordinate); coordinate = NULL;
01494
01495
01496
01497
01498
01499
01500 mapped_sky = mos_wavelength_calibration(smapped, reference,
01501 startwavelength, endwavelength,
01502 dispersion, idscoeff, flux);
01503
01504 cpl_msg_indent_less();
01505 cpl_msg_info(recipe, "Check applied wavelength against skylines...");
01506 cpl_msg_indent_more();
01507
01508 mean_rms = mos_distortions_rms(mapped_sky, lines, startwavelength,
01509 dispersion, 6, highres);
01510
01511 cpl_vector_delete(lines); lines = NULL;
01512
01513 cpl_msg_info(recipe, "Mean residual: %f", mean_rms);
01514
01515 mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01516
01517 cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)",
01518 mean_rms, mean_rms * dispersion);
01519
01520 header = cpl_propertylist_new();
01521 cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01522 cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01523 cpl_propertylist_update_double(header, "CRVAL1",
01524 startwavelength + dispersion/2);
01525 cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01526
01527
01528 cpl_propertylist_update_double(header, "CD1_1", dispersion);
01529 cpl_propertylist_update_double(header, "CD1_2", 0.0);
01530 cpl_propertylist_update_double(header, "CD2_1", 0.0);
01531 cpl_propertylist_update_double(header, "CD2_2", 1.0);
01532 cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01533 cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01534
01535 if (time_normalise) {
01536 dummy = cpl_image_divide_scalar_create(mapped_sky, alltime);
01537 if (dfs_save_image(frameset, dummy, mapped_science_sky_tag, header,
01538 parlist, recipe, version))
01539 fors_science_exit(NULL);
01540 cpl_image_delete(dummy); dummy = NULL;
01541 }
01542 else {
01543 if (dfs_save_image(frameset, mapped_sky, mapped_science_sky_tag,
01544 header, parlist, recipe, version))
01545 fors_science_exit(NULL);
01546 }
01547
01548
01549 if (skyglobal == 0 && skymedian == 0 && skylocal == 0) {
01550 cpl_image_delete(mapped_sky); mapped_sky = NULL;
01551 }
01552
01553 if (skyglobal || skylocal) {
01554
01555 cpl_msg_indent_less();
01556
01557 if (skyglobal) {
01558 cpl_msg_info(recipe, "Global sky determination...");
01559 cpl_msg_indent_more();
01560 skymap = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01561 if (lss || treat_as_lss) {
01562 sky = mos_sky_map_super(spectra, wavemaplss, dispersion,
01563 2.0, 50, skymap);
01564 cpl_image_delete(wavemaplss); wavemaplss = NULL;
01565 }
01566 else {
01567 sky = mos_sky_map_super(spectra, wavemap, dispersion,
01568 2.0, 50, skymap);
01569 }
01570 if (sky) {
01571 cpl_image_subtract(spectra, skymap);
01572 }
01573 else {
01574 cpl_image_delete(skymap); skymap = NULL;
01575 }
01576 }
01577 else {
01578 cpl_msg_info(recipe, "Local sky determination...");
01579 cpl_msg_indent_more();
01580 skymap = mos_subtract_sky(spectra, slits, polytraces, reference,
01581 startwavelength, endwavelength, dispersion);
01582 }
01583
01584 if (skymap) {
01585 if (skyglobal) {
01586 if (time_normalise)
01587 cpl_table_divide_scalar(sky, "sky", alltime);
01588 if (dfs_save_table(frameset, sky, global_sky_spectrum_tag,
01589 NULL, parlist, recipe, version))
01590 fors_science_exit(NULL);
01591
01592 cpl_table_delete(sky); sky = NULL;
01593 }
01594
01595 save_header = dfs_load_header(frameset, science_tag, 0);
01596
01597 if (time_normalise)
01598 cpl_image_divide_scalar(skymap, alltime);
01599 if (dfs_save_image(frameset, skymap, unmapped_sky_tag,
01600 save_header, parlist, recipe, version))
01601 fors_science_exit(NULL);
01602
01603 cpl_image_delete(skymap); skymap = NULL;
01604
01605 if (dfs_save_image(frameset, spectra, unmapped_science_tag,
01606 save_header, parlist, recipe, version))
01607 fors_science_exit(NULL);
01608
01609 cpl_propertylist_delete(save_header); save_header = NULL;
01610
01611 if (cosmics) {
01612 cpl_msg_info(recipe, "Removing cosmic rays...");
01613 mos_clean_cosmics(spectra, gain, -1., -1.);
01614 }
01615
01616
01617
01618
01619
01620
01621 cpl_image_delete(smapped); smapped = NULL;
01622
01623 if (lss || treat_as_lss) {
01624 smapped = cpl_image_duplicate(spectra);
01625 }
01626 else {
01627 smapped = mos_spatial_calibration(spectra, slits, polytraces,
01628 reference, startwavelength,
01629 endwavelength, dispersion,
01630 flux, NULL);
01631 }
01632 }
01633 else {
01634 cpl_msg_warning(recipe, "Sky subtraction failure");
01635 if (cosmics)
01636 cpl_msg_warning(recipe, "Cosmic rays removal not performed!");
01637 cosmics = skylocal = skyglobal = 0;
01638 }
01639 }
01640
01641 cpl_image_delete(spectra); spectra = NULL;
01642 cpl_table_delete(polytraces); polytraces = NULL;
01643
01644 if (skyalign >= 0) {
01645 save_header = dfs_load_header(frameset, science_tag, 0);
01646 if (dfs_save_image(frameset, wavemap, wavelength_map_sky_tag,
01647 save_header, parlist, recipe, version))
01648 fors_science_exit(NULL);
01649 cpl_propertylist_delete(save_header); save_header = NULL;
01650 }
01651
01652 cpl_image_delete(wavemap); wavemap = NULL;
01653
01654 mapped = mos_wavelength_calibration(smapped, reference,
01655 startwavelength, endwavelength,
01656 dispersion, idscoeff, flux);
01657
01658 cpl_image_delete(smapped); smapped = NULL;
01659
01660
01661 if (skymedian) {
01662 cpl_msg_indent_less();
01663 cpl_msg_info(recipe, "Local sky determination...");
01664 cpl_msg_indent_more();
01665
01666
01667
01668 skylocalmap = mos_sky_local_old(mapped, slits);
01669 cpl_image_subtract(mapped, skylocalmap);
01670
01671
01672
01673
01674
01675 cpl_image_delete(skylocalmap); skylocalmap = NULL;
01676 }
01677
01678
01679 if (skyglobal || skymedian || skylocal) {
01680
01681 skylocalmap = cpl_image_subtract_create(mapped_sky, mapped);
01682
01683 cpl_image_delete(mapped_sky); mapped_sky = NULL;
01684
01685 if (time_normalise) {
01686 dummy = cpl_image_divide_scalar_create(skylocalmap, alltime);
01687 if (dfs_save_image(frameset, dummy, mapped_sky_tag, header,
01688 parlist, recipe, version))
01689 fors_science_exit(NULL);
01690 cpl_image_delete(dummy); dummy = NULL;
01691 }
01692 else {
01693 if (dfs_save_image(frameset, skylocalmap, mapped_sky_tag, header,
01694 parlist, recipe, version))
01695 fors_science_exit(NULL);
01696 }
01697
01698 cpl_msg_indent_less();
01699 cpl_msg_info(recipe, "Object detection...");
01700 cpl_msg_indent_more();
01701
01702 if (cosmics || nscience > 1) {
01703 dummy = mos_detect_objects(mapped, slits, slit_margin, ext_radius,
01704 cont_radius);
01705 }
01706 else {
01707 mapped_cleaned = cpl_image_duplicate(mapped);
01708 mos_clean_cosmics(mapped_cleaned, gain, -1., -1.);
01709 dummy = mos_detect_objects(mapped_cleaned, slits, slit_margin,
01710 ext_radius, cont_radius);
01711
01712 cpl_image_delete(mapped_cleaned); mapped_cleaned = NULL;
01713 }
01714
01715 cpl_image_delete(dummy); dummy = NULL;
01716
01717 if (dfs_save_table(frameset, slits, object_table_tag, NULL, parlist,
01718 recipe, version))
01719 fors_science_exit(NULL);
01720
01721 cpl_msg_indent_less();
01722 cpl_msg_info(recipe, "Object extraction...");
01723 cpl_msg_indent_more();
01724
01725 images = mos_extract_objects(mapped, skylocalmap, slits,
01726 ext_mode, ron, gain, 1);
01727
01728 cpl_image_delete(skylocalmap); skylocalmap = NULL;
01729
01730 if (images) {
01731 if (time_normalise)
01732 cpl_image_divide_scalar(images[0], alltime);
01733 if (dfs_save_image(frameset, images[0], reduced_science_tag, header,
01734 parlist, recipe, version))
01735 fors_science_exit(NULL);
01736 cpl_image_delete(images[0]);
01737
01738 if (time_normalise)
01739 cpl_image_divide_scalar(images[1], alltime);
01740 if (dfs_save_image(frameset, images[1], reduced_sky_tag, header,
01741 parlist, recipe, version))
01742 fors_science_exit(NULL);
01743 cpl_image_delete(images[1]);
01744
01745 if (time_normalise)
01746 cpl_image_divide_scalar(images[2], alltime);
01747 if (dfs_save_image(frameset, images[2], reduced_error_tag, header,
01748 parlist, recipe, version))
01749 fors_science_exit(NULL);
01750 cpl_image_delete(images[2]);
01751
01752 cpl_free(images);
01753 }
01754 else {
01755 cpl_msg_warning(recipe, "No objects found: the products "
01756 "%s, %s, and %s are not created",
01757 reduced_science_tag, reduced_sky_tag,
01758 reduced_error_tag);
01759 }
01760
01761 }
01762
01763 cpl_table_delete(slits); slits = NULL;
01764
01765 if (skyalign >= 0) {
01766 if (dfs_save_table(frameset, idscoeff, disp_coeff_sky_tag, NULL,
01767 parlist, recipe, version))
01768 fors_science_exit(NULL);
01769 }
01770
01771 cpl_table_delete(idscoeff); idscoeff = NULL;
01772
01773
01774 if (skyglobal || skymedian || skylocal) {
01775 if (time_normalise)
01776 cpl_image_divide_scalar(mapped, alltime);
01777 if (dfs_save_image(frameset, mapped, mapped_science_tag, header,
01778 parlist, recipe, version))
01779 fors_science_exit(NULL);
01780 }
01781
01782 cpl_image_delete(mapped); mapped = NULL;
01783 cpl_propertylist_delete(header); header = NULL;
01784
01785 if (cpl_error_get_code()) {
01786 cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
01787 fors_science_exit(NULL);
01788 }
01789 else
01790 return 0;
01791 }