xrs_scans.py 98.6 KB
Newer Older
christoph's avatar
christoph committed
1
2
3
#!/usr/bin/python
# Filename: xrs_scans.py

4
5
6
7
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

christoph's avatar
christoph committed
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
#/*##########################################################################
#
# The XRStools software package for XRS spectroscopy
#
# Copyright (c) 2013-2014 European Synchrotron Radiation Facility
#
# This file is part of the XRStools XRS spectroscopy package developed at
# the ESRF by the DEC and Software group.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
#############################################################################*/
__author__ = "Christoph J. Sahle - ESRF"
__contact__ = "christoph.sahle@esrf.fr"
__license__ = "MIT"
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"

41

christoph's avatar
christoph committed
42
import numpy as np
christoph's avatar
christoph committed
43
from . import xrs_utilities, math_functions, xrs_fileIO, xrs_rois
christoph's avatar
christoph committed
44
import h5py
christoph's avatar
christoph committed
45
import os
christoph's avatar
christoph committed
46

christoph's avatar
christoph committed
47
from itertools import groupby
48
from scipy import optimize
49
from scipy.interpolate import Rbf
50
from matplotlib import pylab as plt
christoph's avatar
christoph committed
51
from scipy import ndimage
52

christoph's avatar
christoph committed
53
54
55
56
57
58
59
60
61
62
63
64
65
# try to import the fast PyMCA parsers
try:
    import PyMca5.PyMcaIO.EdfFile as EdfIO
    import PyMca5.PyMcaIO.specfilewrapper as SpecIO
    use_PyMca = True
except:
    use_PyMca = False

print( ' >>>>>>>>  use_PyMca ' , use_PyMca)

__metaclass__ = type # new style classes

class Scan:
66
67
68
69
70
71
72
73
74
75
76
77
78
79
    """ **Scan**

    Class for manipulating scan data from the Hydra and Fourc spectrometers. 

    All relevant information from the SPEC- and EDF-files are organized
    in instances of this class.

    Attributes:
        edf_mats      (np.array): Array containing all 2D images that belong to the scan.
        number             (int): Scan number as in the SPEC file.
        scan_type       (string): Keyword, used later to group scans (add similar scans, etc.).
        energy        (np.array): Array containing the energy axis (1D).
        monitor       (np.array): Array containing the monitor signal (1D).
        counters    (dictionary): Counters with assiciated data from the SPEC file.
mirone's avatar
mirone committed
80
        motors            (dictionary): Motor positions as found in the SPEC file header.
81
82
83
84
85
86
87
88
89
90
91
92
93
        eloss         (np.array): Array of the energy loss scale.
        signals       (np.array): Array of signals extracted from the ROIs.
        errors        (np.array): Array of Poisson errors.
        cenom             (list): Center of mass for each ROI (used if scan is an elastic line scan).
        signals_pw        (list): Pixelwise (PW) data, one array of PW data per ROI.
        errors_pw         (list): Pixelwise (PW) Poisson errors, one array of PW errors per ROI.
        cenom_pw          (list): Center of mass for each pixel.
        signals_pw_interp (list): Interpolated signals for each pixel.
        Ein              (float): Incident energy, used if energy2 is scanned.
        raw_signals       (dict): Dictionary of raw data (summed, line-by-line, pixel-by-pixel).
        raw_errors        (dict): Dictionary of raw erros (summed, line-by-line, pixel-by-pixel).

    """
94
    normalizationDict={} 
95
96
97
98
99
100
101
    def __init__( self ):
        self.edfmats    = np.array([])
        self.scan_number= None
        self.scan_type  = None
        self.energy     = np.array([])
        self.monitor    = np.array([])
        self.counters   = {}
mirone's avatar
mirone committed
102
        self.motors     = {}
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
        self.eloss      = np.array([])
        self.signals    = np.array([])
        self.errors     = np.array([])
        self.cenom      = []
        self.signals_pw = []
        self.errors_pw  = []
        self.cenom_pw   = []
        self.signals_pw_interp = []
        self.Ein        = None
        self.scan_motor = None

        self.raw_signals = {}
        self.raw_errors  = {}

        self.__signals_normalized__ = False

    def load( self, path, SPECfname, EDFprefix, EDFname, EDFpostfix, scan_number, \
                direct=False, roi_obj=None, scaling=None, scan_type='generic', \
                en_column=None, moni_column='izero', method='sum', comp_factor=None,\
122
                rot_angles=None, clean_edf_stack=False, cenom_dict=None, storeInsets=False ):
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
        """ **load**

        Parse SPEC-file and EDF-files for loading a scan.

        Note:
            If 'direct' is 'True' all EDF-files will be deleted 
            after application of the ROIs.

        Args:
            path        (str): Absolute path to directory in which the SPEC-file is located.
            SPECfname   (str): SPEC-file name.
            EDFprefix   (str): Prefix for the EDF-files.
            EDFpostfix  (str): Postfix for the EDF-files.
            scan_number (int): Scan number of the scan to be loaded.
            direct  (boolean): If 'True', all EDF-files will be deleted after loading the scan.
            method      (str): Keyword specifying the selected choice of data treatment:
                can be 'sum', 'row', 'pixel', or 'column'. Default is 'sum'.

        """
        print( 'Parsing EDF- and SPEC-files of scan No. %s.' % scan_number)

        self.scan_number = scan_number

        # load SPEC-file
        fname = os.path.join(path , SPECfname)
        if use_PyMca == True:
            spec_data, self.motors, self.counters, lables = xrs_fileIO.PyMcaSpecRead_my(fname,scan_number)
        else:
            spec_data, self.motors, self.counters = xrs_fileIO.SpecRead(fname,scan_number)

mirone's avatar
mirone committed
153
154
155
156
        # if isinstance(self.motors,dict):
        #     print(" CONVERSIONE 1 ")
        #     self.motors = list( self.motors.items()    )
            
157
158
159
160
161
162
163
164
165
166
        # assign values, energy only if en_column is specified, first counter in SPECfile otherwise
        if en_column:
            self.energy     = np.array(self.counters[en_column.lower()])
            self.scan_motor = en_column.lower()
        else:
            self.energy     = np.array(self.counters[lables[0].lower()])
            self.scan_motor = lables[0].lower()

        # normalization
        the_moni        = np.array(self.counters[moni_column.lower()])
operator for beamline's avatar
operator for beamline committed
167
168
169
170
        
        if moni_column.lower() == 'izero': 
            the_moni  *= np.mean(self.counters['seconds'])
        
171
        if moni_column not in self.normalizationDict:
christoph's avatar
christoph committed
172
            self.normalizationDict[moni_column.lower()] = np.mean(the_moni)/  np.mean(self.counters['seconds'])
operator for beamline's avatar
operator for beamline committed
173
    
174
        self.monitor    = the_moni/self.normalizationDict[moni_column.lower()]
operator for beamline's avatar
operator for beamline committed
175

176
177
178
179
180
181
182
183
184
185
186
        # assign the scan type
        self.scan_type  = scan_type

        # load EDF-files
        if use_PyMca == True:
            self.edfmats = xrs_fileIO.ReadEdfImages_PyMca( self.counters['ccdno'], path, EDFprefix, EDFname, EDFpostfix)
        else:
            self.edfmats = xrs_fileIO.ReadEdfImages_my( self.counters['ccdno'], path, EDFprefix, EDFname, EDFpostfix )

        # remove totally saturated images
        if clean_edf_stack:
187
            self.edfmats = edf_cleaner(self.edfmats, 1.0e8)
188
189
190

        # apply ROIs (if applicable)
        if direct and isinstance( roi_obj, xrs_rois.roi_object ):
myron's avatar
myron committed
191
            self.get_raw_signals( roi_obj, method=method, scaling=scaling, rot_angles=rot_angles , storeInsets = storeInsets)
192

193
194
195
196
197
198
199
200
201
202
203
            if cenom_dict:
                if method == 'row':
                    self.get_signals( method='row', comp_factor=comp_factor, scaling=scaling )
                elif method == 'sum':
                    self.get_signals( method='sum', scaling=scaling )
                elif method == 'pixel':
                    self.get_signals( method='pixel', cenom_dict=cenom_dict )
                elif method == 'pixel2':
                    self.get_signals( method='pixel2', cenom_dict=cenom_dict )
                #elif method == 'column':
                #   self.get_signals( method='column', cenom_dict=cenom_dict )
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227

    def assign( self, edf_arrays, scan_number, energy_scale, monitor_signal, counters, \
                motor_positions, specfile_data, scan_type='generic' ):
        """ **assign**

        Method to group together existing data from a scan (for backward compatibility).

        Args:
            edf_arrays     (np.array): Array of all 2D images that belong to the scan.
            scan_number         (int): Number under which this scan can be found in the SPEC file. 
            energy_scale   (np.array): Array of the energy axis.
            monitor_signal (np.array): Array of the monitor signal.
            counters     (dictionary): Counters with associated data from the SPEC file.
            motor_positions    (list): Motor positions as found in the SPEC file header.
            specfile_data  (np.array): Matrix with all data as found in the SPEC file.
            scantype            (str): Keyword, used later to group scans (add similar scans, etc.).

        """
        self.edfmats    = np.array(edf_arrays)
        self.scan_number= scan_number
        self.energy     = np.array(energy_scale)
        self.monitor    = np.array(monitor_signal)
        self.counters   = counters
        self.motors     = motor_positions
mirone's avatar
mirone committed
228
229
230
231
232
233

        # if isinstance(self.motors,dict):
        #     print(" CONVERSIONE 1 in assign")
        #     self.motors = list( self.motors.items()    )

        
234
235
236
237
238
239
240
241
242
243
244
245
        self.spec_data  = specfile_data
        self.scan_type  = scan_type

    def save_hdf5( self, fname ):
        """ **save_hdf5**
        Save a scan in an HDF5 file.
        Note:
            HDF5 files are strange for overwriting files.

        Args:
            fname (str): Path and filename for the HDF5 file.
        """
myron's avatar
myron committed
246
247
248
249
        if isinstance(fname, h5py.Group):
            f=fname
        else:
            f = h5py.File(fname, "w")
250

251
252
253
254
255
256
257
258
        h5_md = f.require_group("motorDict")
        for mn, mv in self.motors.items():
            h5_md[mn] = mv
            
        h5_md = f.require_group("counters")
        for mn, mv in self.counters.items():
            h5_md[mn] = mv
            
myron's avatar
myron committed
259
        for attr in ['edfmats', 'scan_number', 'energy', 'monitor', 'scan_type','signals','errors']:
260
            f[attr] = eval( 'self.' + attr )
261
            
myron's avatar
myron committed
262
263
264
265
266
267
268
269
270
        for key in self.used_masks.keys():
            hgroup = f.require_group(key)
            pos, mask  = self.used_masks[key]
            hgroup["mask"] = mask
            hgroup["mask_pos"] = pos
            if hasattr( self, "insets") :
                hgroup["insets"]  = self.insets  [key]
        if not isinstance(fname, h5py.Group):
            f.close()
271
272
273
274
275
276
277

    def load_hdf5( self, fname ):
        """ **load_hdf5**
        Load a scan from an HDF5 file.
        Args:
            fname (str): Filename of the HDF5 file.
        """
278

mirone's avatar
mirone committed
279
        doclose = False
280
281
282
        if isinstance(fname, h5py.Group):
            f=fname
        else:
mirone's avatar
mirone committed
283
            doclose = True
284
            f = h5py.File(fname, "r")
285

mirone's avatar
mirone committed
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
        for attr in ['edfmats', 'scan_number', 'energy',  'scan_type', 'signals','errors']:
            setattr( self , attr ,f[attr][()])
    
        self.motors = {}
        if "motorDict"  in f:
            subGroup = f["motorDict"]
            for mname in subGroup:
                self.motors[ mname]  =  subGroup[mname][()]    
                    
        self.counters = {}
        if  "counters" in f:
            subGroup = f["counters"]
            for mname in subGroup:
                    self.counters[mname] =      subGroup[mname][()]  
                    
301
302
303
304
305
        self.used_masks = {}
        for key in f:
            if str(key)[:3] == "ROI" :
                mygroup = f[key]
                self.used_masks[key] = [mygroup["mask_pos"][:].tolist(),  np.array(mygroup["mask"][:])] 
mirone's avatar
mirone committed
306
307
        if doclose :
            f.close()
308
309

        
310

myron's avatar
myron committed
311
    def get_raw_signals( self, roi_obj, method='sum', scaling=None, rot_angles=None, storeInsets=False ):
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
        """ **get_raw_signals**

        Applies given ROIs to EDF-images. 

        Applies the provided ROIs to the EDF-images in the specified mannar:
        summing, line-by-line, or pixel-by-pixel. Depending on the choice, the
        resulting data array is 2D (sum), 3D (line-by-line), or 4D (pixel-by-pixel).
        The scanned direction is always the first dimension of the resulting data 
        matrix.

        Args:
            roi_obj (instance): Instance of the 'XRStools.xrs_rois.roi_object' class defining the ROIs.
            method    (string): Keyword specifying the selected choice of data treatment:
                can be 'sum', 'row', 'pixel', or 'column'. Default is 'sum'.
            scaling (np.array): Array of float-type scaling factors (factor for each ROI).
        Returns:
            None if there are not EDF-files to apply the ROIs to.

        """

myron's avatar
myron committed
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350

        self.used_masks={}

        if storeInsets:
            self.insets={}
            for key, (pos, M) in roi_obj.red_rois.items():
                S     = M.shape
                self.insets[key] =  np.zeros(   [  len(self.edfmats), S[0] , S[1] ],  self.edfmats.dtype   ) 


        for ii in range(len(self.edfmats)):
            for key, (pos, M) in roi_obj.red_rois.items():
                S     = M.shape
                inset = (slice( pos[0], pos[0]+(S[0]) ), slice( pos[1], pos[1]+(S[1]) ))
                self.used_masks[key] = (   pos, M )
                if storeInsets:
                    self.insets[key][ii] = self.edfmats[ii, inset[0], inset[1]] * (M/M.max())

                
351
352
353
354
        # sum
        if method == 'sum':
            print('selected method is \'sum\': summing up pixels from each ROI.')
            signals = {} # dict (one entry per ROI, with vector (one entry per energy point))
christoph's avatar
christoph committed
355
            errors  = {} # sqrt of the sum of counts
356
357
            for key, (pos, M) in roi_obj.red_rois.items():
                signals[key] = np.zeros((len(self.energy)))
christoph's avatar
christoph committed
358
                errors[key]  = np.zeros((len(self.energy)))
359
360
361
362
363
364

            for ii in range(len(self.edfmats)):
                ind = 0
                for key, (pos, M) in roi_obj.red_rois.items():
                    S     = M.shape
                    inset = (slice( pos[0], pos[0]+(S[0]) ), slice( pos[1], pos[1]+(S[1]) ))
myron's avatar
myron committed
365
                    
366
                    signals[key][ii] = np.sum( self.edfmats[ii, inset[0], inset[1]] * (M/M.max()))
christoph's avatar
christoph committed
367
                    errors[key][ii]  = np.sqrt(signals[key][ii])
368
                    signals[key][ii] /= self.monitor[ii]
christoph's avatar
christoph committed
369
                    errors[key][ii]  /= self.monitor[ii]
370
371
372
373
374
375
                    ind += 1

        # row
        elif method == 'row':
            print('selected method is \'row\': summing over non-dispersive direction for each ROI.')
            signals = {} # dict (one entry per ROI, with 2D matrix (energy vs row))
christoph's avatar
christoph committed
376
            errors  = {} # sqrt of the sum of counts
377
378
            rot_angles_dict = {} # put possible rotation angles into dict
            counter = 0
379
            for key, (pos, M) in sorted(roi_obj.red_rois.items()   ):
380
                signals[key] = np.zeros((len(self.energy), M.shape[0]))
christoph's avatar
christoph committed
381
                errors[key]  = np.zeros((len(self.energy), M.shape[0]))
382
383
384
385
386
387
388
389
390
                if rot_angles is not None:
                    rot_angles_dict[key] = rot_angles[counter]
                    counter += 1

            for ii in range(len(self.edfmats)):
                ind = 0
                for key, (pos, M) in roi_obj.red_rois.items():
                    S     = M.shape
                    inset = (slice( pos[0], pos[0]+(S[0]) ), slice( pos[1], pos[1]+(S[1]) ))
myron's avatar
myron committed
391
392
                    
                    
393
394
395
396
397
398
399
400
401
402
403
404
405
                    # rotate raw_signals and raw_errors if method is 'line' and angles are provided
                    if rot_angles:
                        if len(roi_obj.red_rois) is not len(rot_angles):
                            print('Only %d rotation angles provided for %d ROIs. Will end here.'%(len(rot_angles), len(self.raw_signals)))
                            return

                        # rotate images before summation
                        orig_slice    = self.edfmats[ii, inset[0], inset[1]] * (M/M.max())
                        slice_for_sum = ndimage.interpolation.rotate( orig_slice, rot_angles_dict[key],\
                                        reshape=False, order=0, mode='constant' )
                    else:
                        slice_for_sum = self.edfmats[ii, inset[0], inset[1]] * (M/M.max())
                    signals[key][ii,:] = np.sum( slice_for_sum , axis=1)
christoph's avatar
christoph committed
406
                    errors[key][ii,:]  = np.sqrt(signals[key][ii,:])
407
                    signals[key][ii,:] /= self.monitor[ii]
christoph's avatar
christoph committed
408
                    errors[key][ii,:]  /= self.monitor[ii]
409
410
411
412
413
414
                    ind += 1

        # column
        elif method == 'column':
            print('selected method is \'column\': summing over dispersive direction for each ROI.')
            signals = {} # dict (one entry per ROI, with 2D matrix (energy vs row))
christoph's avatar
christoph committed
415
            errors  = {} # sqrt of the sum of counts
416
417
            for key, (pos, M) in roi_obj.red_rois.items():
                signals[key] = np.zeros((len(self.energy), M.shape[1]))
christoph's avatar
christoph committed
418
                errors[key]  = np.zeros((len(self.energy), M.shape[1]))
419
420
421
422
423
424

            for ii in range(len(self.edfmats)):
                ind = 0
                for key, (pos, M) in roi_obj.red_rois.items():
                    S     = M.shape
                    inset = (slice( pos[0], pos[0]+(S[0]) ), slice( pos[1], pos[1]+(S[1]) ))
myron's avatar
myron committed
425
426

                    
427
                    signals[key][ii,:] = np.sum( self.edfmats[ii, inset[0], inset[1]] * (M/M.max()), axis=0)
christoph's avatar
christoph committed
428
                    errors[key][ii,:]  = np.sqrt( signals[key][ii,:] )
429
                    signals[key][ii,:] /= self.monitor[ii]
christoph's avatar
christoph committed
430
                    errors[key][ii,:]  /= self.monitor[ii]
431
432
433
434
435
436
                    ind += 1

        # pixel
        elif method == 'pixel' or method == 'pixel2':
            print('selected method is \'pixel\': returning ROI pixel-by-pixel.')
            signals = {} # dict (one entry per ROI, with 3D matrix (energy vs pixel_0 vs pixel_1))
christoph's avatar
christoph committed
437
            errors  = {} # sqrt of the sum of counts
438
439
            for key, (pos, M) in roi_obj.red_rois.items():
                signals[key] = np.zeros((len(self.energy), M.shape[0], M.shape[1]))
christoph's avatar
christoph committed
440
                errors[key]  = np.zeros((len(self.energy), M.shape[0], M.shape[1]))
441
442
443
444
445
446

            for ii in range(len(self.edfmats)):
                ind = 0
                for key, (pos, M) in roi_obj.red_rois.items():
                    S     = M.shape
                    inset = (slice( pos[0], pos[0]+(S[0]) ), slice( pos[1], pos[1]+(S[1]) ))
myron's avatar
myron committed
447
448

                    
449
                    signals[key][ii,:,:] = self.edfmats[ii, inset[0], inset[1]] * (M/M.max())
christoph's avatar
christoph committed
450
                    errors[key][ii,:,:]  = np.sqrt( signals[key][ii,:,:]  )
451
                    signals[key][ii,:,:] /= self.monitor[ii]
christoph's avatar
christoph committed
452
                    errors[key][ii,:,:]  /= self.monitor[ii]
453
454
455
456
457
458
459
460
461
462
463
                    ind += 1

        # unknown method
        else:
            print( 'Unknown integration method. Use either \'sum\', \'row\', or \'pixel\'.' )
            return

        # set normalization 
        self.__signals_normalized__ = True

        # assign
464
        for key,ii in zip(sorted(roi_obj.red_rois, key = lambda x:  int(''.join(filter(str.isdigit, str(x) ))) ), range(len(roi_obj.red_rois))):
465
466
            if np.any(scaling):
                self.raw_signals[key] = signals[key] * scaling[ii]
christoph's avatar
christoph committed
467
                self.raw_errors[key]  = errors[key]  * scaling[ii]
468
469
            else:
                self.raw_signals[key] = signals[key]
christoph's avatar
christoph committed
470
                self.raw_errors[key]  = errors[key]
471

christoph's avatar
christoph committed
472
        # delete EDF-files (these are obsolete after pixel-by-pixel integration)
mirone's avatar
mirone committed
473
        print( 'Deleting ++ EDF-files of scan No. %s.' % self.scan_number )
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
        self.edfmats = np.array([])

    def get_signals( self, method='sum', cenom_dict=None, comp_factor=None, scaling=None, PIXEL_SIZE=0.055 ):
        """ **get_signals**

        Turns pixel-, column- or sum-wise raw-data into data.

        Takes the raw-data after application of the ROIs and applies the chosen compensation
        scheme.

        Args:
            method        (str): Keyword specifying the selected choice of data treatment:
                can be 'sum', 'row', or 'pixel'. Default is 'sum'.
            cenom_dict   (dict): Dictionary with one entry per ROI holding information about
                the center of mass of the according elastic line.
            comp_factor (float): Factor used in the RIXS-style line-by-line compensation.
            scaling  (np.array): Array of float-type scaling factors (factor for each ROI).

        """

494
        if method == 'sum':
495
496
            self.signals = np.zeros(( len(self.energy), len(self.raw_signals) ))
            self.errors  = np.zeros(( len(self.energy), len(self.raw_signals) ))
497
            for key,ii in zip(sorted(self.raw_signals, key = lambda x:  int(''.join(filter(str.isdigit, str(x) ))) ), range(len(self.raw_signals))):
498
499
500
501
502
503
504
                if not self.__signals_normalized__:
                    self.signals[:,ii] = self.raw_signals[key]/self.monitor
                    self.errors[:,ii]  = self.raw_errors[key]/self.monitor
                else:
                    self.signals[:,ii] = self.raw_signals[key]
                    self.errors[:,ii]  = self.raw_errors[key]

505
        elif method == 'pixel':
506
507
            # assign sign of compensation
            if self.scan_motor == 'energy':
508
                comp_direction = 1.0
509
510
511
512
513
514
515
            elif self.scan_motor == 'anal energy':
                comp_direction = -1.0
            else:
                print('Unknown energy motor, will break here.')
                return
            self.signals = np.zeros(( len(self.energy), len(self.raw_signals) ))
            self.errors  = np.zeros(( len(self.energy), len(self.raw_signals) ))
516
            for key,ii in zip(sorted(self.raw_signals, key = lambda x:  int(''.join(filter(str.isdigit, str(x) ))) ), range(len(self.raw_signals))):
517
518
519
520
                S = cenom_dict[key].shape
                master_cenom =  cenom_dict[key][int(S[0]/2.),int(S[1]/2.)]
                for dim1 in range(self.raw_signals[key].shape[1]):
                    for dim2 in range(self.raw_signals[key].shape[2]):
christoph's avatar
christoph committed
521
                        x    = self.energy + ( comp_direction*cenom_dict[key][dim1,dim2] - comp_direction*master_cenom )
522
523
524
525
526
527
528
529
530
531
532
533
                        y    = self.raw_signals[key][:,dim1,dim2]
                        rbfi = Rbf( x, y, function='linear' )
                        self.signals[:,ii] += rbfi(self.energy)
                        y    = self.raw_errors[key][:,dim1,dim2]
                        rbfi = Rbf( x, y, function='linear' )
                        self.errors[:,ii] += rbfi(self.energy)**2
                if not self.__signals_normalized__:
                    self.errors[:,ii] = np.sqrt( self.errors[:,ii] )/self.monitor
                    self.signals[:,ii] /= self.monitor
                else:
                    self.errors[:,ii] = np.sqrt( self.errors[:,ii] )

534
        elif method == 'pixel2':
535
536
            # assign sign of compensation
            if self.scan_motor == 'energy':
537
                comp_direction = 1.0
538

539
540
541
542
543
544
545
546
            elif self.scan_motor == 'anal energy':
                comp_direction = -1.0
            else:
                print('Unknown energy motor, will break here.')
                return
            self.signals = np.zeros(( len(self.energy), len(self.raw_signals) ))
            self.errors  = np.zeros(( len(self.energy), len(self.raw_signals) ))
            energy = self.energy #* 1e3 # energy in eV
operator for beamline's avatar
operator for beamline committed
547
            ## meanmon = np.mean(self.monitor)
548
            for key,ii in zip(sorted(self.raw_signals, key = lambda x:  int(''.join(filter(str.isdigit, str(x) ))) ), range(len(self.raw_signals))):
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
                S = cenom_dict[key].shape
                #the_hist = np.histogram(cenom_dict[key][cenom_dict[key]>0.0 ], bins=10)
                #master_cenom =  np.average(the_hist[1][1:], weights= the_hist[0])#cenom_dict[key][int(S[0]/2.),int(S[1]/2.)]
                master_cenom = np.mean(cenom_dict[key][cenom_dict[key]>0.0 ])
                y  = self.raw_signals[key]
                yn = y# (y.T/self.monitor).T * meanmon
                yc = np.zeros(( len(energy),S[0]*S[1] ))
                for dim1 in range(S[0]):
                    for dim2 in range(S[1]):
                        sort  = np.argsort(energy)
                        if cenom_dict[key][dim1,dim2] > 0.0:
                            yc[:,dim1*dim2] = np.interp( energy[sort] -(comp_direction*(cenom_dict[key][dim1,dim2] - \
                                                    master_cenom)), energy[sort], yn[sort,dim1,dim2], \
                                                    left=float('nan'), right=float('nan') )
                        else:
                            yc[:,dim1*dim2] = yn[sort,dim1,dim2]
                self.signals[:,ii] = xrs_utilities.nonzeroavg(yc)
                self.errors[:,ii]  = np.sqrt(self.signals[:,ii])
                # and normalize to I0
                if not self.__signals_normalized__:
                    self.signals[:,ii] /= self.monitor
                    self.errors[:,ii] /= self.monitor

572
        elif method == 'row':
573

574
575
576
577
578
579
580
581
582
583
584
            # assign sign of compensation !!! test this !!!
            if self.scan_motor == 'energy':
                comp_direction = 1.0
            elif self.scan_motor == 'anal energy':
                comp_direction = -1.0
            else:
                print('Unknown energy motor, will break here.')
                return
            self.signals = np.zeros(( len(self.energy), len(self.raw_signals) ))
            self.errors  = np.zeros(( len(self.energy), len(self.raw_signals) ))
            energy = self.energy * 1e3 # energy in eV
operator for beamline's avatar
operator for beamline committed
585
            ### meanmon = np.mean(self.monitor)
586
            for key,ii in zip(sorted(self.raw_signals, key = lambda x:  int(''.join(filter(str.isdigit, str(x) ))) ), range(len(self.raw_signals))):
587
                y  = self.raw_signals[key]
operator for beamline's avatar
operator for beamline committed
588
                yn = (y.T/self.monitor).T #### * meanmon CHECK THIS REMOVAL
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
                meanii = len(range(yn.shape[1]))//2
                yc = np.zeros_like(y)
                for jj in range(yn.shape[1]):
                    sort  = np.argsort(energy)
                    yc[:,jj] = np.interp( energy[sort] + (jj-meanii)*comp_direction*comp_factor*PIXEL_SIZE, energy[sort], yn[sort,jj], 
                                        left=float('nan'), right=float('nan') )
                self.signals[:,ii] = xrs_utilities.nonzeroavg(yc)
                self.errors[:,ii]  = np.sqrt(self.signals[:,ii])
                # and normalize to I0
                if not self.__signals_normalized__:
                    self.signals[:,ii] /= self.monitor
                    self.errors[:,ii] /= self.monitor

        else:
            print( 'Unknown integration method. Use either \'sum\', \'row\', or \'pixel\'.' )
            return

        # set normalization
        self.__signals_normalized__ = True

    def apply_rois( self, roi_obj, scaling=None ):
        """ **apply_rois**

        Sums up intensities in each ROI.

        Note:
            Old function, keeping for backward compatibility.

        Args:
            roi_obj (instance): Instance of the 'XRStools.xrs_rois.roi_object' class defining the ROIs.
            scaling (np.array): Array of float-type scaling factors (factor for each ROI).

        Returns:
            None if there are not EDF-files to apply the ROIs to.

        """

        # make sure there is data loaded
        if self.edfmats.size == 0:
            print( 'Please load some EDF-files first.' )
            return

        # apply ROIs
        signals = np.zeros((len(self.energy), len(roi_obj.red_rois)))
        for ii in range(len(self.edfmats)):
            ind = 0
            for key, (pos, M) in roi_obj.red_rois.items():
                S     = M.shape
                inset = (slice( pos[0], pos[0]+(S[0]) ), slice( pos[1], pos[1]+(S[1]) ))
                signals[ii,ind] = np.sum( self.edfmats[ii, inset[0], inset[1]] * (M/M.max()))
                ind += 1

        # assign
        self.signals = signals
        self.errors  = np.sqrt(signals)

        # apply scaling (if applicable)
        if np.any(scaling):
            # make sure, there is one scaling factor for each ROI
            assert len(scaling) == signals.shape[1]
            for ii in range(signals.shape[1]):
                self.signals[:,ii] *= scaling[ii]
                self.errors[:,ii]  *= scaling[ii]

    def apply_rois_pw( self,roi_obj, scaling=None ):
        """ **apply_rois_pw**

        Pixel-wise reading of the ROIs' pixels into a list of arrays.

        I.e. each n-pixel ROI will have n Spectra, saved in a 2D array.

        Args:
        roi_obj (instance): Instance of the 'XRStools.xrs_rois.roi_object' class defining the ROIs.
        scaling (list) or (np.array): Array or list of float-type scaling factors (one factor for each ROI).

        """
        data   = [] # list of 2D arrays (energy vs. intensity for each pixel inside a single ROI) 
        errors = []
        for n in range(len(roi_obj.indices)): # each ROI
            roidata = np.zeros((len(self.edfmats),len(roi_obj.indices[n]))) # 2D np array energy vs pixels in current roi
            for m in range(len(self.edfmats)): # each energy point along the scan
                for l in range(len(roi_obj.indices[n])): # each pixel on the detector
                    roidata[m,l] = self.edfmats[m,roi_obj.indices[n][l][0],roi_obj.indices[n][l][1]]
            data.append(roidata) # list which contains one array (energy point, pixel) per ROI
            errors.append(np.sqrt(roidata))

        self.signals_pw = data
        self.errors_pw  = errors

        # apply scaling (if applicable)
        if np.any(scaling):
            # make sure, there is one scaling factor for each ROI
            assert len(scaling) == len(roi_obj.indices)
            for ii in range(len(roi_obj.indices)):
                self.signals_pw[ii] *= scaling[ii]
                self.errors_pw[ii]  *= scaling[ii]

    def append_scan( self, scan, method='sum', where='right' ):
        """ **append_scan**

        Appends scan to the current scan, either at higher energies 
        (where='right') or at lower energies (where='left').

        Args:
            scan       (obj): Object of the Scan class.
            method     (str): Keyword specifying the selected choice of data treatment:
                can be 'sum', 'row', or 'pixel'. Default is 'sum'.
            where      (str): Keyword specifying if the scan should be appended
697
698
                at lower (where='left') or highger energies (where='righ'). Default is 
                'right'.        
699
 
700
701
        """

702
        if where == 'right':
703

704
705
706
707
708
            self.energy  = np.append( self.energy , scan.energy )
            self.monitor = np.append( self.monitor, scan.monitor )
            for key in self.raw_signals:
                self.raw_signals[key] = np.append( self.raw_signals[key], scan.raw_signals[key], axis=0 )
                self.raw_errors[key]  = np.append( self.raw_errors[key], scan.raw_errors[key], axis=0 )
709
        if where == 'left':
710

711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
            self.energy  = np.append( scan.energy,  self.energy )
            self.monitor = np.append( scan.monitor, self.monitor )
            for key in self.raw_signals:
                self.raw_signals[key] = np.append( scan.raw_signals[key], self.raw_signals[key], axis=0 )
                self.raw_errors[key]  = np.append( scan.raw_errors[key], self.raw_errors[key], axis=0 )

    def insert_scan( self, scan, method='sum', where=None ):
        """ **insert_scan**

        Inserts another scan into the current instance.

        Args:
            scan       (obj): Object of the Scan class.
            method     (str): Keyword specifying the selected choice of data treatment:
                can be 'sum', 'row', or 'pixel'. Default is 'sum'.
            where     (list): Optional tuple of energy values (high and low (in keV))
                for where to insert the scan. By default (None), the lowest and highest
                energy values of the given scan will be used. 

        """

        # find where given scan should be inserted
        if not where:
            low_inds  = np.where( self.energy < scan.get_E_start() )[0]
            high_inds = np.where( self.energy > scan.get_E_end() )[0]
        else:
            low_inds  = np.where( self.energy < where[0] )[0]
738
            high_inds = np.where( self.energy > where[1] )[0]       
739
740

        # 'sum'
741
        if method == 'sum':
742

743
744
745
746
747
748
749
750
751
            self.energy  = np.append( self.energy[low_inds], np.append( scan.energy, self.energy[high_inds] ) )
            self.monitor = np.append( self.monitor[low_inds], np.append( scan.monitor, self.monitor[high_inds] ) )
            for key in self.raw_signals:
                self.raw_signals[key] = np.append( self.raw_signals[key][low_inds], 
                                        np.append( scan.raw_signals[key], self.raw_signals[key][high_inds], axis=0 ), axis=0 )
                self.raw_errors[key] = np.append( self.raw_errors[key][low_inds], 
                                        np.append( scan.raw_errors[key], self.raw_errors[key][high_inds], axis=0 ), axis=0 )

        # 'pixel'
752
        elif method in [ 'pixel' , 'pixel2']:
753

754
755
756
757
758
759
760
761
762
            self.energy  = np.append( self.energy[low_inds], np.append( scan.energy, self.energy[high_inds] ) )
            self.monitor  = np.append( self.monitor[low_inds], np.append( scan.monitor, self.monitor[high_inds] ) )
            for key in self.raw_signals:
                self.raw_signals[key] = np.append( self.raw_signals[key][low_inds,:,:], 
                                        np.append( scan.raw_signals[key], self.raw_signals[key][high_inds,:,:], axis=0 ), axis=0 )
                self.raw_errors[key] = np.append( self.raw_errors[key][low_inds,:,:], 
                                        np.append( scan.raw_errors[key], self.raw_errors[key][high_inds,:,:], axis=0 ), axis=0 )

        # 'row'
763
        elif method == 'row':
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
            self.energy  = np.append( self.energy[low_inds], np.append( scan.energy, self.energy[high_inds] ) )
            self.monitor  = np.append( self.monitor[low_inds], np.append( scan.monitor, self.monitor[high_inds] ) )
            for key in self.raw_signals:
                self.raw_signals[key] = np.append( self.raw_signals[key][low_inds,:], 
                                        np.append( scan.raw_signals[key], self.raw_signals[high_inds,:], axis=0 ), axis=0 )
                self.raw_errors[key] = np.append( self.raw_errors[key][low_inds,:], 
                                        np.append( scan.raw_errors[key], self.raw_errors[high_inds,:], axis=0 ), axis=0 )

        else:
            print( 'Unknown method. Use either \'sum\', \'row\', or \'pixel\'.' )
            return

    def add_scan(self, scan, method='sum', interp=False):
        """ **add_scan**

        Adds signals from a different scan.

        Args:
            scan       (obj): Object of the Scan class.
            method     (str): Keyword specifying the selected choice of data treatment:
                can be 'sum', 'row', or 'pixel'. Default is 'sum'.
            interp (boolean): Boolean specifying if norm-conserving linear wavelet
                interpolation should be used, False by default. 

        """
789
790
791
        # @@@@@@@@@@@@@@@@@@   DA VERIFICARE TUTTO L'ACCOUNTING
        # @@@@@@@@@@@@@@@@@@   changed by christoph: 13/07/2018
        # @@@@@@@@@@@@@@@@@@        please double check
792
793

        if not interp:
794
795
796
797
            
            # sum monitors
            av_monitor = np.sum((self.monitor, scan.monitor), axis=0)

798
            for key in self.raw_signals:
christoph's avatar
christoph committed
799
                if method in [ 'sum', 'row', 'column']:
800
801
802
803
804
805
806
807
808
809
810
811

                    # recover raw counts
                    signals1 = self.raw_signals[key]*self.monitor
                    signals2 = scan.raw_signals[key]*scan.monitor

                    # sum signals
                    av_signals = np.sum((signals1, signals2), axis=0)

                    # errors of summed signals
                    av_errors  = np.sqrt(av_signals)

                    # assign and renormalize
operator for beamline's avatar
operator for beamline committed
812
                    self.raw_signals[key] = av_signals/av_monitor
christoph's avatar
christoph committed
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
                    self.raw_errors[key]  = av_errors/av_monitor

                if method in [ 'pixel', 'pixel2']:

                    # recover raw counts
                    signals1 = self.raw_signals[key]*self.monitor[:,None,None]
                    signals2 = scan.raw_signals[key]*scan.monitor[:,None,None]

                    # sum signals
                    av_signals = np.sum((signals1, signals2), axis=0)
              
                    # errors of summed signals
                    av_errors  = np.sqrt(av_signals)

                    # assign and renormalize
                    self.raw_signals[key] = av_signals/av_monitor[:,None,None]
                    self.raw_errors[key]  = av_errors/av_monitor[:,None,None]
830

831
        if interp:
832
833
834
835
836
837
838
839
840
            # find longest scan dimension
            dim0 = np.amax([len(self.energy), len(scan.energy)])

            # sum the monitor signal (always 1D)
            mm = np.zeros((dim0, 2))*np.nan
            mm[0:len(self.monitor),0] = self.monitor
            mm[0:len(scan.monitor),1] = scan.monitor
            av_monitor = np.nansum( mm , axis=1)
            
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
            # add also unfinished scans
            for key in self.raw_signals:

                # construct a matrix for nanmean depending on the shape of the raw_signals
                if method in ['sum']:
                    # recover raw counts
                    signals1 = self.raw_signals[key]*self.monitor
                    signals2 = scan.raw_signals[key]*scan.monitor
                    # signals
                    yy = np.zeros((dim0, 2))*np.nan
                    yy[0:len(signals1),0] = signals1
                    yy[0:len(signals2),1] = signals2
                    av_signals = np.nansum( yy , axis=1)
                    # errors of summed signals
                    av_errors  = np.sqrt(av_signals)
                    #
                    self.raw_signals[key] = av_signals/av_monitor
                    self.raw_errors[key]  = av_errors/av_monitor

                if method in ['pixel', 'pixel2']:
                    # recover raw counts
                    signals1 = self.raw_signals[key]*self.monitor[:,None,None]
                    signals2 = scan.raw_signals[key]*scan.monitor[:,None,None]
                    # signals
                    y1 = np.zeros((dim0, signals1.shape[1], signals1.shape[2]))*np.nan
                    y2 = np.zeros((dim0, signals1.shape[1], signals1.shape[2]))*np.nan
                    y1[0:len(signals1),:,:] = signals1
                    y2[0:len(signals2),:,:] = signals2
                    av_signals = np.nansum( (y1,y2) , axis=0 )
                    # errors of summed signals
                    av_errors  = np.sqrt( av_signals )
                    #
                    self.raw_signals[key] = av_signals/av_monitor[:,None,None]
                    self.raw_errors[key]  = av_errors/av_monitor[:,None,None]

                if method in ['row']:
                    # recover raw counts
                    signals1 = self.raw_signals[key]*self.monitor[:,None]
                    signals2 = scan.raw_signals[key]*scan.monitor[:,None]
                    # signals
                    y1 = np.zeros((dim0, signals1.shape[1]))*np.nan
                    y2 = np.zeros((dim0, signals1.shape[1]))*np.nan
                    y1[0:len(signals1),:,:] = signals1
                    y2[0:len(signals2),:,:] = signals2
                    av_signals = np.nansum( (y1,y2) , axis=0 )
                    # errors of summed signals
                    av_errors  = np.sqrt( av_signals )
                    #
                    self.raw_signals[key] = av_signals/av_monitor[:,None]
                    self.raw_errors[key]  = av_errors/av_monitor[:,None]

        if interp=='Rbf':
893
894
            #rbfi = Rbf( scan.energy, scan.monitor, function='linear' )
            #self.monitor += rbgi( self.energy )
895
            for key in self.raw_signals:
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
                # remove zeros in errors
                self.raw_errors[key][self.raw_errors[key] == 0 ] = 1.0
                scan.raw_errors[key][scan.raw_errors[key] == 0 ] = 1.0

                # recover raw counts
                signals1 = self.raw_signals[key]*self.monitor
                signals2 = scan.raw_signals[key]*scan.monitor

                # sum interpolated signals
                rbfi = Rbf( scan.energy, signals2, function='linear' )
                av_signals = signals1 + rbfi( self.energy )

                # sum interpolated monitors
                rbfi = Rbf( scan.energy, scan.monitor, function='linear' )
                av_monitor = self.monitor + rbfi( self.energy )

                # errors of summed signals
                av_errors = np.sqrt(av_signals)
914

915
                # assign and renormalize
operator for beamline's avatar
operator for beamline committed
916
917
                self.raw_signals[key] = av_signals/av_monitor 
                self.raw_errors[key]  = av_errors/av_monitor  
918
                
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
    def get_type(self):
        """ **get_type**

        Returns the type of the scan.

        """
        return self.scan_type

    def get_scan_number(self):
        """ **get_scan_number**

        Returns the number of the scan.

        """
        return self.scan_number

    def get_shape(self):
        """ **get_shape**

        Returns the shape of the matrix holding the signals.

        """
        if not np.any(self.signals):
            print( 'please apply the ROIs first.' )
            return
        else:
            return self.signals.shape

    def get_num_of_rois(self):
        """ **get_num_of_rois**

        Returns the number of ROIs applied to the scan.

        """
        if not self.signals.any():
            print( 'please apply the ROIs first.' )
            return
        else:
            return self.signals.shape[1]

    def get_E_start(self):
        """ **get_E_start**

        Returs the first energy value.

        """
        return self.energy[0]

    def get_E_end(self):
        """ **get_E_end**

        Returs the last energy value.

        """
        return self.energy[-1]

    def get_resolution( self, keV2eV=True ):
        """ **get_resolution**

        Returns the ROI-wise resolution based on the 
        xrs_utilities fwhm method.
        """
        resolutions = []
        for ii in range(self.signals.shape[1]):
            x = self.energy
            y = self.signals[:,ii]
            if keV2eV:
                resolutions.append( xrs_utilities.fwhm(x,y)[0]*1.0e3 )
            else:
                resolutions.append( xrs_utilities.fwhm(x,y)[0] )
        return np.array(resolutions), np.mean(resolutions), np.std(resolutions)
christoph's avatar
christoph committed
990

christoph's avatar
christoph committed
991

christoph's avatar
christoph committed
992
class scan:
993
994
995
996
997
998
999
1000
    """
    Container class, holding information of single scans performed with 2D detectors. 
    """
    def __init__(self,edf_arrays,scannumber,energy_scale,monitor_signal,counters,motor_positions,specfile_data,scantype='generic'):
        # rawdata
        self.edfmats  = np.array(edf_arrays)          # numpy array of all 2D images that belong to the scan
        self.number   = scannumber                    # number under which this scan can be found in the SPEC file 
        self.scantype = scantype                      # keyword, later used to group scans (add similar scans, etc.)