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