File: flames_obs_scired_interact.py

package info (click to toggle)
cpl-plugin-uves 6.1.3+dfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 23,128 kB
  • sloc: ansic: 171,056; sh: 4,359; python: 3,002; makefile: 1,322
file content (464 lines) | stat: -rwxr-xr-x 22,808 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
from __future__ import absolute_import
from __future__ import print_function
def log_invocation(fn='/tmp/reflex_invoc.log'):
    """ log invocation of actor for easy reexecution with:
        eval $(cat /tmp/reflex_invoc.log)
    """
    import os, sys
    #if 'REFLEX_DEBUG' not in os.environ:
    #    return
    try:
        with open(fn, 'w') as f:
            path = ""
            s = " ".join(['"%s"' % x for x in sys.argv])
            try:
                import reflex
                path = '"%s:"' % os.path.dirname(reflex.__file__)
            except ImportError:
                pass
            f.write('PYTHONPATH=%s$PYTHONPATH python %s\n' % (path, s)) 
    except:
        pass




# import the needed modules
try:
  import matplotlib
  import reflex_plot_widgets
  from pipeline_product import *
  matplotlib.use('WXAgg')
  import_sucess = 'true'

#NOTE for developers: 
# -If you want to modify the current script to cope
#  with different parameters, this is the function to modify:
#  setInteractiveParameters()
# -If you want to modify the current script to read different data from
#  the input FITS, this is the function to modify:
#  readFitsData()                  (from class DataPlotterManager) 
# -If you want to modify the current script to modify the plots (using the same
#  data),  this is the function to modify:
#  plotProductsGraphics()          (from class DataPlotterManager)
# -If you want to modify the text that appears in the "Help" button,
#  this is the function to modify:
#  setWindowHelp()
# -If you want to modify the title of the window, modify this function:
#  setWindowTitle()


  #This class deals with the specific details of data reading and final plotting.
  class DataPlotterManager:
    # This function will read all the columns, images and whatever is needed
    # from the products. The variables , self.plot_x, self.plot_y, etc...
    # are used later in function plotProductsGraphics().
    # Add/delete these variables as you need (only that plotProductsGraphics()
    # has to use the same names).
    # You can also create some additional variables (like statistics) after
    # reading the files.
    # If you use a control variable (self.xxx_found), you can modify 
    # later on the layout of the plotting window based on the presence of 
    # given input files. 
    # sof contains all the set of frames
    def readFitsData(self, fitsFiles):
      #Control variable to ced_heck if the interesting files where at the input
      self.red_sci_point_mos_found = False
      self.sci_info_tab_mos_found = False
      self.cubify = False

      #Read all the products
      from collections import defaultdict
      frames = defaultdict(list)

      list_hdus_sci_redl = list()
      list_hdus_sci_redu = list()
      list_hdus_err_sci_redl = list()
      list_hdus_err_sci_redu = list()
      list_hdus_sci_info_tab_redu = list()
      list_sci_redl = list()
      list_sci_redu = list()
      list_err_sci_redl = list()
      list_err_sci_redu = list()
      list_sci_info_tab_redl = list()
      list_sci_info_tab_redu = list()

      for frame in fitsFiles:
        if frame == '' :
          continue
        pipe_product = PipelineProduct(frame)
        hdulist = pipe_product.hdulist()  
        hdus_item = dict()
        #hdulist = frame.hdulist()
        category = frame.category
        frames[category].append(frame)
        #print category, frame.name

        if category == "MWXB_SCI_REDL" :
           list_sci_redl.append(frame)
           hdus_item[category]=hdulist
           list_hdus_sci_redl.append(hdus_item[category])
           #print category

        elif category == "MWXB_SCI_REDU" :
           list_sci_redu.append(frame)
           hdus_item[category]=hdulist
           list_hdus_sci_redu.append(hdus_item[category])
           #print category

        elif category == "ERR_MWXB_SCI_REDL" :
           list_err_sci_redl.append(frame)
           hdus_item[category]=hdulist
           list_hdus_err_sci_redl.append(hdus_item[category])           
           #print category

        elif category == "ERR_MWXB_SCI_REDU" :
           list_err_sci_redu.append(frame)
           hdus_item[category]=hdulist
           list_hdus_err_sci_redu.append(hdus_item[category])
           #print category

        elif category == "FIB_SCI_INFO_TAB_REDL":
           self.sci_info_tab_mos_found = True  
           list_sci_info_tab_redl.append(frame)
           hdus_item[category]=hdulist
           hdu_sci_info_tab_redl = hdus_item[category]

           #get WLEN and SIMCAL setting
           table = pipeline_product.PipelineProduct(frame)
           self.wlen = table.all_hdu[0].header['ESO INS GRAT2 WLEN']
           self.is_simcal = table.all_hdu[0].header['ESO OCS SIMCAL']

        elif "FIB_SCI_INFO_TAB_REDU" in frames :
           self.sci_info_tab_mos_found = True  
           list_sci_info_tab_redu.append(frame)
           hdus_item[category]=hdulist
           list_hdus_sci_info_tab_redu.append(hdus_item[category])
 
      naxis=list_hdus_sci_redl[0][0].header['NAXIS']
      if  naxis == 2  :
          nfl = list_hdus_sci_redl[0].data.shape[0]
          nfu = list_hdus_sci_redu[0].data.shape[0]
          self.cubify = True
          print("Sci Fibre cube format")
      else :
          nfl = len(list_sci_redl)
          nfu = len(list_sci_redu)
          print("Sci fibre normal format")
      
      #print "nf=", nf , "cubify = ", self.cubify

      # rarely happens that no same fibres are extracted on both
      # REDL/U arms. The following trick allow the python display
      # not to crush
      #nf = min(nfl,nfu)
      self.nfibres =  min(nfl,nfu)
      self.list_spec_raw = []
      self.list_snr_raw = []
      #i=1
  
      # Read fibre info table a first time to know if the fibre is a
      # simultaneous calibration fibre or allocated to object or sky
      self.list_obs_obj_cal = []
      self.list_obs_obj_type = []
      for i in range(9):
         self.list_obs_obj_cal.append(hdu_sci_info_tab_redl[1].data['Retractor_pt1'][i])
         self.list_obs_obj_type.append(hdu_sci_info_tab_redl[1].data['TYPE'][i])

     
      # Read fibre info table a second time to get object id and
      # corresponding RA,DEC only for simultaneous calibration fibre
      # and fibres allocated to object or sky      
      self.list_obs_objects = []
      self.list_obs_obj_ra = []
      self.list_obs_obj_dec = []
      for i in range(9):
         if self.list_obs_obj_cal[i] == "Calibration"  and self.is_simcal == 1 :
           self.list_obs_objects.append(hdu_sci_info_tab_redl[1].data['OBJECT'][i])

           self.list_obs_obj_ra.append(hdu_sci_info_tab_redl[1].data['RA'][i])
           self.list_obs_obj_dec.append(hdu_sci_info_tab_redl[1].data['DEC'][i])
    
         if self.list_obs_obj_type[i] == "U" or self.list_obs_obj_type[i] == "S" :
           self.list_obs_objects.append(hdu_sci_info_tab_redl[1].data['OBJECT'][i])
           self.list_obs_obj_ra.append(hdu_sci_info_tab_redl[1].data['RA'][i])
           self.list_obs_obj_dec.append(hdu_sci_info_tab_redl[1].data['DEC'][i])
          
           
      #print self.list_obs_objects

      self.list_used_fibres = []
      for i in range(self.nfibres):
        #table=list_hdus_sci_info_tab_redl[0].data[i,:]
        self.red_sci_point_mos_found = True
        if self.cubify == True :
           print("cube case")
           #import pdb
           #pdb.set_trace()
           self.list_spec_raw.append(UvesRedSpectrum(list_hdus_sci_redl[0].data[i,:],list_hdus_sci_redu[0].data[i,:],list_hdus_err_sci_redl[0].data[i,:],list_hdus_err_sci_redu[0].data[i,:]))
        else :
           #print "normal case i=",i,self.list_obs_obj_cal[i],self.list_obs_obj_type[i]
           self.list_used_fibres.append(i)

           #print "add to list",  i,self.list_obs_objects[i] 
           self.list_spec_raw.append(FlamesUvesSpectrum(frames["MWXB_SCI_REDL"],
                                                  frames["MWXB_SCI_REDU"],i,
                                                  frames["ERR_MWXB_SCI_REDL"],
                                                  frames["ERR_MWXB_SCI_REDU"]))
           #self.list_snr_raw.append(FlamesUvesSN(frames["MWXB_SCI_REDL"],
           #                                       frames["MWXB_SCI_REDU"],i,
           #                                       frames["ERR_MWXB_SCI_REDL"],
           #                                       frames["ERR_MWXB_SCI_REDU"]))
 
        #print "no of fibres", self.nfibres
        #print "allocated fibres", self.list_used_fibres


    def plotWidgets(self) :
        widgets = list()
        #labels= [0,1,2,3,4,5,6,7]
        labels = self.list_used_fibres
        # Files Selector radiobutton
        self.radiobutton = reflex_plot_widgets.InteractiveRadioButtons(self.fibre_selector, self.setFibreSelectCallBack, labels, 0, 
                title="")
        widgets.append(self.radiobutton)

        return widgets

    def setFibreSelectCallBack(self, fibre_id) :
        fibre_id = int(fibre_id)
        self.list_subplot_spec[0].cla()
        #self.list_subplot_snr[0].cla()
        #print self.list_spec_raw[fibre_id]
        self.plotData(fibre_id)

    # This function creates all the subplots. It is responsible for the plotting 
    # layouts. 
    # There can different layouts, depending on the availability of data
    # Note that subplot(I,J,K) means the Kth plot in a IxJ grid 
    # Note also that the last one is actually a box with text, no graphs.
    def addSubplots(self, figure):
      #gs = gridspec.GridSpec(3, 1)
      #self.fibre_selector = figure.add_subplot(gs[0,0])
      if self.red_sci_point_mos_found == True :
        self.list_subplot_spec     = []
        self.list_subplot_snr       = []

        nc=1
        nf=self.nfibres
        nraws=2
        #print "nf=",nf
        i=self.list_used_fibres[0]
        #for i in range(nf) :
        self.list_subplot_spec.append(figure.add_subplot(nraws,nc,nc*i+1))
        #self.list_subplot_snr.append(figure.add_subplot(nraws,nc,nc*i+2))

      else : 
        self.subtext_nodata     = figure.add_subplot(1,1,1)
        
      #self.fibre_selector = figure.add_axes([0.85, 0.85, 0.14, 0.14])
      self.fibre_selector = figure.add_axes([0.85, 0.25, 0.14, 0.14])
      tip="Frame Id:\n (Mouse \n left\n button)"
      self.fibre_selector.text(0.35, 0.50, tip, color='black', fontsize=12, ha='left', va='center', alpha=1.0)
    def plotData(self,fibre_id):
        self.setInfoObj(fibre_id)
        self.plotSpectrum(fibre_id)
        #self.plotSNRSpectrum(fibre_id)

    def setInfoCommon(self):
        self.title_common   = 'Extracted and Merged Spectrum (REDL,REDU).\n'

        self.title_snr   = 'Signal to Noise Ratio of Extracted and Merged Spectrum (REDL,REDU)' 
        self.tooltip_spec = """Plot of the extracted and merged spectrum of the science object (blue line) as total flux (ADU) versus wavelength (Ang). 
The +-1 sigma uncertainties are plotted as the light blue region encompassing the object spectrum (and bounded by the light grey spectra). 
Note that this spectrum is not flux calibrated."""

        self.tooltip_snr = """Plot of the Signal to Noise ratio associated to the extracted and merged spectrum of the science object (blue line) as total flux (ADU) versus wavelength (Ang). 
The +-1 sigma uncertainties are plotted as the light blue region encompassing the object spectrum (and bounded by the light grey spectra). 
Note that this spectrum is not flux calibrated."""

    def setInfoSetting(self):
        if self.is_simcal == 1 :
           self.info_setting = 'SimCal %#.3d mode. '% (self.wlen) 
        else :
           self.info_setting = 'OzPoz %#.3d mode. '% (self.wlen) 

    def setInfoObj(self,fibre_id):
        if self.list_obs_obj_cal[fibre_id] == "Calibration" and self.is_simcal == 1 :
           self.info_obj='Simultaneous calibration fibre'
        else :
           self.info_obj='Target: '+self.list_obs_objects[fibre_id]+' RA=%#.4g DEC=%#.4g'% (self.list_obs_obj_ra[fibre_id],self.list_obs_obj_dec[fibre_id])

    def plotSpectrum(self,fibre_id):
        self.title_spec=self.title_common+self.info_setting+self.info_obj
        title=self.title_spec
        tooltip=self.tooltip_spec
        self.list_spec_raw[fibre_id].plot(self.list_subplot_spec[0], title, tooltip)

    def plotSNRSpectrum(self,fibre_id):
        title=self.title_snr
        tooltip=self.tooltip_snr
        self.list_snr_raw[fibre_id].plot(self.list_subplot_snr[0], title, tooltip)

    # This is the function that makes the plots.
    # Add new plots or delete them using the given scheme.
    # The data has been already stored in self.plot_x, self.plot_xdif, etc ...
    # It is mandatory to add a tooltip variable to each subplot.
    # One might be tempted to merge addSubplots() and plotProductsGraphics().
    # There is a reason not to do it: addSubplots() is called only once at
    # startup, while plotProductsGraphics() is called always there is a resize.
    def plotProductsGraphics(self):
      if self.red_sci_point_mos_found == True :
        i=self.list_used_fibres[0] 

        #First subpanel: a spectrum
        if self.red_sci_point_mos_found == True :
           self.setInfoCommon()
           self.setInfoSetting()
           self.setInfoObj(i)
           self.plotData(i)


      else :
        #Data not found info
        self.subtext_nodata.set_axis_off()
        self.text_nodata = """Science object spectra not found in the products:
For Red data:  PRO.CATG=MWXB_SCI_REDL, MWXB_SCI_REDU"""
        self.subtext_nodata.text(0.1, 0.6, self.text_nodata, color='#11557c', fontsize=18,
                                 ha='left', va='center', alpha=1.0)
        self.subtext_nodata.tooltip='Science object spectra not found in the products'

  
    # This function specifies which are the parameters that should be presented
    # in the window to be edited.
    # Note that the parameter has to be also in the in_sop port (otherwise it 
    # won't appear in the window) 
    # The descriptions are used to show a tooltip. They should match one to one
    # with the parameter list 
    # Note also that parameters have to be prefixed by the 'recipe name:'
    def setInteractiveParameters(self):
      paramList = list()
#      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='debug','general','Whether or not to save intermediate results to local directory. [FALSE]'))
#      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='plotter','general',
#      """Any plots produced by the recipe are redirected to the
#         command specified by this parameter. The plotting command must
#         contain the substring 'gnuplot' and must be able to parse
#         gnuplot syntax on its standard input. Valid examples of such
#         a command may include 'gnuplot -persist' and 'cat >
#         mygnuplot$$.gp'. A finer control of the plotting options can
#         be obtained by writing an executable script, e.g. my_gnuplot.pl, that 
#         executes gnuplot after setting the desired gnuplot options
#         (e.g. set terminal pslatex color). To turn off plotting, set
#         this parameter to 'no'. [no]"""))

#      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='process_chip',group='general',"""For
#      RED arm data proces the redl, redu, or both chip(s). <both |
#      redl | redu | REDL | REDU> [both]"""))

      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='ext_method',group='extract',description="""Extraction method. <std | opt | fst | fop> [opt]"""))

#      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='cor_max_fnd',group='extract',description="""Find correlation maximum?. <N | Y> [Y]"""))

      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='cor_max_ rng',group='extract',description="""Correlation range size?. [6.0]"""))

      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='cor_max_ pnt',group='extract',description="""Correlation sampling points?. [25]"""))

      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='cor_def_off',group='extract',description="""Correlation center offset?. [0.0]"""))

      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='corvel_iter',group='extract',description="""Velocity correlation iteration number (SimCal). [1]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='bias_method',group='extract',description="""Bias subtraction method. <M | V | N> [M]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='bias_value',group='extract',description="""Bias value (only if bias_method = V). [200]"""))
#      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='cubify_sw',group='general',description="""Cubify switch. <Y | N> [N]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='filt_sw',group='extract',description="""Filter switch. <none | median> [none]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='bkg_max_io_win',group='background',description="""Background window number in each full inter order. [500]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='bkg_xy_win_sz_x',group='background',description="""x maximum size of each background window:. [6]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='bkg_xy_win_sz_y',group='background',description="""y maximum size of each background window:. [2]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='pixel_thresh_max',group='BP',description="""Pixel saturation threshold max. [55000]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='pixel_thresh_min',group='BP',description="""Pixel saturation threshold min. [-20]""")) 
#     paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='input_fmt_cube',group='general',description="""Input data format. [TRUE]"""))
#     paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='output_fmt_cube',group='general',description="""Output data format. [FALSE]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='drs_k_s_thre',group='general',description="""Kappa sigma threshold. [10.0]"""))
#      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='drs_base_name',group='general',description="""Base name for science products. [fxb]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='drs_maxyshift',group='extract',description="""Half width of the interval to scan for correlation, when determining y shift. [3.0]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='drs_ext_w_siz',group='extract',description="""Integration window size good: 10 (if fibre deconvolution works fine). [10.0]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='rebin.wavestep',group='rebin',description="""The bin size (in w.l.u.) in wavelength space. If negative, a step size of 2/3 * ( average pixel size ) is used. [-1.0]"""))
      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='rebin.scale',group='rebin',
      description="""Whether
      or not to multiply by the factor dx/dlambda (pixels per
      wavelength) during the rebinning. This option is disabled as
      default in concordance with the method used in the MIDAS
      pipeline. This option should be set to true to convert the
      observed flux (in pixel-space) to a flux per wavelength (in 
      wavelength-space). [FALSE]"""))

      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='merge',group='merge',
      description="""Order merging method. If 'optimal', the flux in the
      overlapping region is set to the (optimally computed, using the 
      uncertainties) average of single order spectra. If 'sum', the
      flux in the overlapping region is computed as the sum of the
      single order spectra. If flat-fielding is done, method 'optimal'
      is recommended, otherwise 'sum'. <optimal | sum> [optimal]"""))

      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='merge_delt1',group='merge',
      description="""Order merging left hand (short wavelength) cut. To reduce the
      amount of order overlapping regions we allow to cut short and
      long wavelength ranges. This may reduce the ripple possibly
      introduced by the order merging. Suggested values are: 10 
      (W<=390), 12 (390<W<=437, 520<W<=564), 14 (437<W<=520, 564<W). [0.0]"""))

      paramList.append(RecipeParameter(recipe='flames_obs_scired',displayName='merge_delt2',group='merge',
      description="""Order merging right hand (long wavelength) cut. To reduce the
      amount of order overlapping regions we allow to cut short and
      long wavelength ranges. This may reduce the ripple possibly
      introduced by the order merging. Suggested values is 4. [0.0]"""))


      return paramList

    def setWindowHelp(self):
      help_text = """
In this window, the user should check that the science object extracted spectrum is of a good quality by using the pan and zoom buttons at the top-left of this window. 
Attempt to optimise the S/N of the extracted spectrum as a function of spectral order (the upper plot of the bottom-left plots) by choosing different parameter values and re-running the pipeline recipe."""
      return help_text

    def setWindowTitle(self):
      title = 'Uves Interactive Spectrum Extraction'
      return title


except ImportError:
  import_sucess = 'false'
  print("Error importing modules pyfits, wx, matplotlib, numpy")

#This is the 'main' function
if __name__ == '__main__':
  log_invocation()

  # import reflex modules
  from reflex import *
  from reflex_interactive_app import *
  from pipeline_display import *

  # import UVES reflex modules
  from flames_plot_common import *

  # Create interactive application
  interactive_app = PipelineInteractiveApp(enable_init_sop=True)

  #Check if import failed or not
  if import_sucess == 'false' :
    interactive_app.setEnableGUI('false')

  #Open the interactive window if enabled
  if interactive_app.isGUIEnabled() :
    #Get the specific functions for this window
    dataPlotManager = DataPlotterManager()
    interactive_app.setPlotManager(dataPlotManager)
    interactive_app.showGUI()
  else :
    interactive_app.passProductsThrough()

  # print outputs
  interactive_app.print_outputs()

  sys.exit()