xrs_read.py

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from six.moves import range
from six.moves import zip
from six.moves import input
#!/usr/bin/python
# Filename: xrs_read.py

#/*##########################################################################
#
# 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"

from . import xrs_rois, xrs_scans, xrs_utilities, math_functions, xrs_fileIO, roifinder_and_gui

import h5py
import scipy.io
import traceback
import sys
import os
import numpy as np
import array as arr
import pickle
import matplotlib.pyplot as plt

from numpy import array
from itertools import groupby
from scipy.integrate import trapz
from scipy.interpolate import interp1d, Rbf
from scipy import signal, optimize
from scipy.ndimage import measurements

# 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

def print_citation_message():
	"""Prints plea for citing the XRStools article when using this software.

	"""
	print ('                                                                                ')
	print (' ############################# Welcome to XRStools #############################')
	print (' # If you are using this software, please cite the following work:             #')
	print (' # Ch.J. Sahle, A. Mirone, J. Niskanen, J. Inkinen, M. Krisch, and S. Huotari: #') 
	print (' # "Planning, performing, and analyzing X-ray Raman scattering experiments."   #')
	print (' # Journal of Synchrotron Radiation 22, No. 2 (2015): 400-409.                 #')
	print (' ###############################################################################')
	print ('                                                                                ')

print_citation_message()

class Hydra:
    """Main class for handling XRS data from ID20's multi-analyzer spectrometer 'Hydra'.

    This class is intended to read SPEC- and according EDF-files and generate spectra from 
    multiple individual energy loss scans.

    Note:
        Hydra is the name of the multi-analyzer x-ray Raman scattering spectrometer at ESRF's
        ID20 beamline. This class has been adopted specifically for this spectrometer.

        If you are using this program, please cite the following work:

        Sahle, Ch J., A. Mirone, J. Niskanen, J. Inkinen, M. Krisch, and S. Huotari. 
        "Planning, performing and analyzing X-ray Raman scattering experiments." 
        Journal of Synchrotron Radiation 22, No. 2 (2015): 400-409.

    Args:
          * path        (str): Absolute path to directory holding the data.
          * SPECfname   (str): Name of the SPEC-file ('hydra' is the default).
          * EDFprefix   (str): Prefix for the EDF-files ('/edf/' is the default).
          * EDFname     (str): Filename of the EDF-files ('hydra_' is the default).
          * EDFpostfix  (str): Postfix for the EDF-files ('.edf' is the default).
          * en_column   (str): Counter mnemonic for the energy motor ('energy' is the default).
          * moni_column (str): Mnemonic for the monitor counter ('izero' is the default).

    Attributes:
          * path        (str): Absolute path to directory holding the data.
          * SPECfname   (str): Name of the SPEC-file ('hydra' is the default).
          * EDFprefix   (str): Prefix for the EDF-files ('/edf/' is the default).
          * EDFname     (str): Filename of the EDF-files ('hydra_' is the default).
          * EDFpostfix  (str): Postfix for the EDF-files ('.edf' is the default).
          * en_column   (str): Counter mnemonic for the energy motor ('energy' is the default).
          * moni_column (str): Mnemonic for the monitor counter ('izero' is the default).
          * scans        (dict): Dictionary holding all loaded scans.
          * scan_numbers (list): List with scan number of all loaded scans.
          * eloss    (np.array): Array with the common energy loss scale for all analyzers.
          * energy   (np.array): Array with the common energy scale for all analyzers.
          * signals  (np.array): Array with the signals for all analyzers (one column per anayzer).
          * errors   (np.array): Array with the poisson errors for all analyzers (one column per anayzer).
          * qvalues      (list): List of momentum transfer values for all analyzers.
          * groups       (dict): Dictionary of groups of scans (instances of the 'scangroup' class,
            such as 2 'elastic', or 5 'edge1', etc.).
          * cenom        (list): List of center of masses of the elastic lines.
          * E0            (int): Elastic line energy value, mean value of all center of masses.
          * tth          (list): List of all scattering angles (one value for each ROI).
          * resolution   (list): List of FWHM of the elastic lines (one for each analyzer).
          * comp_factor (float): Compensation factor for line-by-line energy dispersion compensation.
          * cenom_dict   (dict): Dictionary holding center-of-masses for of the elastic line.
          * raw_signals  (dict): Dictionary holding pixel- or line-wise signals.
          * raw_errors   (dict): Dictionary holding pixel- or line-wise Poisson errors.
          * TTH_OFFSETS1 (np.array): Two-Theta offsets between individual analyzers inside each analyzer 
            module in one direction (horizontal for V-boxes, vertical for H-boxes).
          * TTH_OFFSETS2 (np.array): Two-Theta offsets between individual analyzers inside each analyzer 
            module in one direction (horizontal for H-boxes, vertical for V-boxes).
          * roi_obj (instance): Instance of the roi_object class from the xrs_rois module defining all
            ROIs for the current dataset (default is 'None').

    """

    def __init__( self, path, SPECfname='hydra', EDFprefix='/edf/', EDFname='hydra_', \
                        EDFpostfix='.edf', en_column='energy', moni_column='izero' ):

        self.path        = path
        self.SPECfname   = SPECfname

        if not os.path.isfile(os.path.join(path, SPECfname)):
            raise Exception('IOError! No such file or directory.')

        self.EDFprefix   = EDFprefix
        self.EDFname     = EDFname
        self.EDFpostfix  = EDFpostfix
        self.en_column   = en_column.lower()
        self.moni_column = moni_column.lower()

        self.scans        = {}
        self.scan_numbers = []

        self.eloss    = np.array([])
        self.energy   = np.array([])
        self.signals  = np.array([])
        self.errors   = np.array([])
        self.q_values = []
        self.groups   = {}
        self.cenom    = []
        self.E0       = []
        self.tth      = []
        self.resolution  = []
        self.comp_factor = None
        self.cenom_dict  = {}
        self.raw_signals = {}
        self.raw_errors  = {} 

        self.TTH_OFFSETS1 = np.array([5.0, 0.0, -5.0, 5.0, 0.0, -5.0, 5.0, 0.0, -5.0, 5.0, 0.0, -5.0])
        self.TTH_OFFSETS2 = np.array([-9.71, -9.75, -9.71, -3.24, -3.25, -3.24, 3.24, 3.25, 3.24, 9.71, 9.75, 9.71]) 
        self.PIXEL_SIZE   = 0.055 # pixel size in mm

        self.roi_obj = None

        print_citation_message()

    def save_state_hdf5( self, file_name, group_name, comment="" ):
        """ **save_state_hdf5**

        Save the status of the current instance in an HDF5 file.

        Args:
            file_name  (str): Path and file name for the HDF5-file to be created.
            group_name (str): Group name under which to store status in the HDF5-file.
            comment    (str): Optional comment (no comment is default).

        """
        h5 = h5py.File(file_name,"a")

        h5.require_group(group_name)
        h5group =  h5[group_name]

        for key in self.__dict__.keys():
            if key in h5group.keys():
                raise Exception( 'Data \'' + key + '\' already present in  ' + file_name + ':' + group_name )
            else:
                h5group[key] = getattr( self, key )

        h5group["comment"]  = comment
        h5.flush()
        h5.close()

    def load_state_hdf5( self, file_name, group_name ):
        """ **load_state_hdf5**

        Load the status of an instance from an HDF5 file.

        Args:
            file_name  (str): Path and filename for the HDF5-file to be created.
            group_name (str): Group name under which to store status in the HDF5-file.

        """
        h5 = h5py.File( file_name,"r" )

        h5group =  h5[group_name]
        keys    = {"eloss":array, "energy":array, "signals":array, "errors":array, "q_values":array, 
                    "cenom":array, "E0":float, "tth":array, "resolution":array }

        for key in keys:
            setattr(self, key, keys[key](array(h5group[key])))

        h5.flush()
        h5.close()

    def set_roiObj( self,roiobj ):
        """ **set_roiObj**

        Assign an instance of the 'roi_object' class to the current data set.

        Args:
            roiobj (instance): Instance of the 'roi_object' class holding all
                information about the definition of the ROIs.

        """
        self.roi_obj = roiobj

    def load_scan( self, scan_numbers, scan_type='generic', direct=True, scaling=None, method='sum' ):
        """**load_scan**

        Load a single or multiple scans.
        Note:
            When composing a spectrum later, scans of the same scan_type will  
            be averaged over. Scans with scan type 'elastic' or long in their
            names are recognized and will be treated specially.

        Args:
              * scan_numbers (int or list): Integer or iterable of scan numbers to be loaded.
              * scan_type            (str): String describing the scan to be loaded (e.g. 'edge1' or 'K-edge').
              * direct           (boolean): Flag, 'True' if the EDF-files should be deleted after loading/integrating the scan.


        """

        # make sure scannumbers are iterable (list)
        numbers = []
        if not isinstance(scan_numbers,list):
            numbers.append(scan_numbers)
        else:
            numbers = scan_numbers

        # make sure there is a cenom_dict available if direct=True AND method='sum'
        if direct and method=='pixel' and not self.cenom_dict:
            print('Please run the get_compensation_factor method first for pixel-wise compensation.')
            return

        # go throught list of scan_numbers and load scans		
        for number in numbers:

            # create a name for each scan
            scan_name = 'Scan%03d' % number

            # create an instance of the Scan class
            scan = xrs_scans.Scan()

            # load the scan
            scan.load( self.path, self.SPECfname, self.EDFprefix, self.EDFname, self.EDFpostfix, number, \
                direct=direct, roi_obj=self.roi_obj, scaling=scaling, scan_type=scan_type, \
                en_column=self.en_column, moni_column=self.moni_column, method=method, \
                cenom_dict=self.cenom_dict, comp_factor=self.comp_factor )

            # add it to the scans dict
            self.scans[scan_name] = scan

            # add the number to the list of scan numbers
            if not number in self.scan_numbers:
                self.scan_numbers.extend([number])

    def load_loop( self, beg_nums, num_of_regions, direct=True, method='sum' ):
        """ **load_loop**

        Loads a whole loop of scans based on their starting numbers and 
        the number of single scans in the loop.

        Args:
        beg_nums      (list): List of scan numbers of the first scans in each loop.
        num_of_regions (int): Number of scans in each loop.

        """

        type_names = []
        for n in range(num_of_regions):
            type_names.append('edge'+str(n+1))

        numbers = []
        for n in range(len(beg_nums)):
            for m in range(num_of_regions):
                numbers.append((beg_nums[n]+m))

        type_names = type_names*len(beg_nums)

        for n in range(len(type_names)):
            number = []
            number.append(numbers[n])
            self.load_scan( number, type_names[n], direct=True, method=method )

    def delete_scan( self, scan_numbers ):
        """ **delete_scan**

        Deletes scans from the dictionary of scans.

        Args:
            scan_numbers (int or list): Integer or list of integers (SPEC scan numbers) to be deleted.

        """

        # make sure scannumbers are iterable (list)
        numbers = []
        if not isinstance(scan_numbers,list):
            numbers.append(scan_numbers)
        else:
            numbers = scan_numbers

        # delete the scan
        for number in numbers:
            scan_name = 'Scan%03d' % number
            del(self.scans[scan_name])
            self.scan_numbers.remove(number)

    def get_raw_data( self, method='sum', scaling=None ):
        """ **get_raw_data**

        Applies the ROIs to the EDF-files. 

        This extracts the raw-data from the EDF-files subject to three 
        different methods: 'sum' will simply sum up all pixels inside a
        ROI, 'row' will sum up over the dispersive direction, and 'pixel'
        will not sum at all.

        Args:
            method (str): Keyword specifying the selected choice of data treatment:
                can be 'sum', 'row', or 'pixel'. Default is 'sum'.

        """

        if not self.roi_obj:
            print ( 'Did not find a ROI object, please set one first.' )
            return

        for scan in self.scans:
            if len(self.scans[scan].edfmats):
                print ( "Integrating " + scan + " using method \'" + method + '\'.')
                self.scans[scan].get_raw_signals( self.roi_obj , method=method, scaling=scaling )


    def get_data_new(self, method='sum', scaling=None):
        """ **get_data_new**

        Applies the ROIs to the EDF-files. 

        This extracts the raw-data from the EDF-files subject to three 
        different methods: 'sum' will simply sum up all pixels inside a
        ROI, 'row' will sum up over the dispersive direction, and 'pixel'
        will not sum at all.

        Args:
            method (str): Keyword specifying the selected choice of data treatment:
                can be 'sum', 'row', or 'pixel'. Default is 'sum'.

        """

        if not self.cenom_dict:
            print ( 'No compensation factor/center-of-mass info found, please provide first.' )
            return

        for scan in self.scans:
            self.scans[scan].get_signals( method=method, cenom_dict=self.cenom_dict, comp_factor=self.comp_factor, scaling=scaling )

    def get_data(self):
        """ **get_data**

        Applies the ROIs and sums up intensities.

        Returns:
            'None', if no ROI object is available.

        """

        if not self.roi_obj:
            print ( 'Did not find a ROI object, please set one first.' )
            return

        for scan in self.scans:
            if len(self.scans[scan].edfmats):
                print ( "Integrating " + scan )
                self.scans[scan].apply_rois( self.roi_obj )

    def get_data_pw(self):
        """ **get_data_pw**

        Extracts intensities for each pixel in each ROI.

        Returns:
            'None', if no ROI object is available.

        """

        if not np.any(self.roi_obj.indices):
            print ( 'Did not find a ROI object, please set one first.' )
            return

        for scan in self.scans:
            if len(self.scans[scan].edfmats):
                print ( "Integrating pixelwise " + scan )
                self.scans[scan].apply_rois_pw(self.roi_obj)

    def SumDirect( self, scan_numbers, index=None ):
        """ **SumDirect**

        Creates a summed 2D image of a given scan or list of scans.

        Args:
            scan_numbers (int or list): Scan number or list of scan numbers to be added up.

        Returns:
            A 2D np.array of the same size as the detector with the summed image.

        """

        # make sure scannumbers are iterable (list)
        numbers = []
        if not isinstance(scan_numbers,list):
            numbers.append(scan_numbers)
        else:
            numbers = scan_numbers

        im_sum    = None
        en_column = None # uses first column in SPEC file as scanned motor
        for number in numbers:
            scan = xrs_scans.Scan()
            scan.load(self.path, self.SPECfname, self.EDFprefix, self.EDFname, self.EDFpostfix, number, \
                    direct=False, roi_obj=None, scaling=None, scan_type='generic', \
                    en_column=en_column, moni_column=self.moni_column)

            if im_sum is None:
                im_sum1 = np.zeros(scan.edfmats[0].shape ,"f")
                im_sum2 = np.zeros(scan.edfmats[0].shape ,"f")
            if not index:
                im_sum1[:] += scan.edfmats.sum(axis=0)
            else:
                im_sum1[:] += scan.edfmats[0:index,:,:].sum(axis=0)
                im_sum2[:] += scan.edfmats[index:,:,:].sum(axis=0)

        if not index:
            return im_sum1
        else:
            return im_sum2-im_sum1

    def get_eloss_new(self, method='sum'):
        """ **get_eloss_new**

        Defines the energy loss scale for all ROIs and applies dispersion 
        compensation if applicable.

        Args:
            method (str): Keyword describing which dispersion compensation 
                method to use. Possible choices are 'sum' (no compensation), 
                'pixel' (pixel-by-pixel compensation), or 'row' (line-by-line)
                compensation.

        """

        # make sure an elastic line was loaded
        if not 'elastic' in self.groups:
            print( 'Please load/integrate at least one elastic scan first!' )
            return

        # make sure there is data
        elif not np.any(self.groups['elastic'].raw_signals):
            self.get_raw_data( method=method )

        # make sure there is center of masses available
        elif not self.cenom_dict and  method != 'row':
            for scan in self.scans:
                if self.scans[scan].get_type() == 'elastic':
                    print('GETIING COMP FACTOR!!!')
                    self.get_compensation_factor( self.scans[scan].scan_number, method=method )

        first_key = list(self.raw_signals.keys())[0]
        # 'sum'
        if method == 'sum':
            # master eloss scale in eV is the one of the first ROI
            self.signals = np.zeros((len(self.energy),len(self.cenom_dict)))
            self.errors  = np.zeros((len(self.energy),len(self.cenom_dict)))
            master_eloss = (self.energy - np.median([self.cenom_dict[key] for key in  self.cenom_dict]))*1.0e3
            self.E0      = np.median([self.cenom_dict[key] for key in  self.cenom_dict])
            for key,ii in zip(sorted(self.cenom_dict), range(len(self.cenom_dict))):
                # signals
                x = ( self.energy - self.cenom_dict[key] )*1.0e3
                y = self.raw_signals[key][:]
                #try:
                #	rbfi = Rbf( x, y, function='linear' )
                #	self.signals[:,ii] = rbfi( master_eloss )
                #except:
                self.signals[:,ii] = np.interp( master_eloss, x, y )
                # errors
                y = self.raw_errors[key][:]
                #try:
                #	rbfi = Rbf( x, y, function='linear' )
                #	self.errors[:,ii] = rbfi( master_eloss )
                #except:
                self.errors[:,ii] = np.interp( master_eloss, x, y )
            self.eloss = master_eloss

        # 'pixel'
        elif method ==  'pixel':
            # master eloss scale in eV is the one of central pixel in first ROI
            self.signals = np.zeros((len(self.energy),len(self.cenom_dict)))
            self.errors  = np.zeros((len(self.energy),len(self.cenom_dict)))
            master_eloss = ( self.energy - np.median(self.cenom_dict[first_key][self.cenom_dict[first_key] > 0.0]) )*1.0e3
            self.E0      = np.median(self.cenom_dict[first_key][self.cenom_dict[first_key] > 0.0])
            for key,ii in zip(sorted(self.cenom_dict), range(len(self.cenom_dict))):
                print ('Pixel-by-pixel compensation for ' + key +'.')
                signal = np.zeros(len(master_eloss))
                error  = np.zeros(len(master_eloss))
                for dim1 in range(self.cenom_dict[key].shape[0]):
                    for dim2 in range(self.cenom_dict[key].shape[1]):
                        x = ( self.energy - self.cenom_dict[key][dim1, dim2] )*1.0e3
                        # signals
                        y = self.raw_signals[key][:, dim1, dim2]
                        if np.any(y)>0.0:
                            #print "Y AMAX", signal.max()
                            #rbfi    = Rbf( x, y, function='linear' )
                            rbfi = interp1d(x, y,bounds_error=False, fill_value=0.0)
                            #print "rbf AMAX", Rbf( x, y, function='linear' )
                            #signal += rbfi( master_eloss )
                            signal += rbfi( master_eloss )
                            #print "SIGNAL AMAX", signal.max()
                            # errors
                            y = self.raw_errors[key][:, dim1, dim2]
                            rbfi    = Rbf( x, y, function='linear' )
                            error  += rbfi( master_eloss )**2
                self.signals[:,ii] = signal
                self.errors[:,ii]  = np.sqrt(error)
            self.eloss = master_eloss

        # 'row'
        elif method == 'row':
            self.signals = np.zeros((len(self.energy),len(self.roi_obj.red_rois)))
            self.errors  = np.zeros((len(self.energy),len(self.roi_obj.red_rois)))
            energy = self.energy * 1e3 # energy in eV

            for key,ii in zip(sorted(self.raw_signals), range(len(self.raw_signals))):
                y  = self.raw_signals[key]
                meanii = len(range(y.shape[1]))/2
                yc = np.zeros_like(y)
                for jj in range(y.shape[1]):
                    sort  = np.argsort(energy)
                    yc[:,jj] = np.interp( energy[sort] + (jj-meanii)*self.comp_factor*self.PIXEL_SIZE, energy[sort], y[sort,jj], 
                                        left=float('nan'), right=float('nan') )
                self.signals[:,ii] = xrs_utilities.nonzeroavg(yc)
                self.errors[:,ii]  = np.sqrt(self.signals[:,ii])

        else:
            print('Method \''+method+'\' not supported, use either \'sum\', \'pixel\', or \'row\'.')
            return

    def get_eloss( self ):
        """ **get_eloss**

        Finds the energy loss scale for all ROIs by calculating the center 
        of mass (COM) for each ROI's elastic line. Calculates the resolution 
        function (FWHM) of the elastic lines.

        """

        # make sure an elastic line was loaded
        if not 'elastic' in self.groups:
            print( 'Please load/integrate at least one elastic scan first!' )
            return

        # make sure there is data
        elif not self.groups['elastic'].signals.shape[1] == len(self.roi_obj.indices):
            self.get_data()

        else:
            # reset values, in case this is run several times
            self.cenom      = []
            self.resolution = []
            Origin          = None
            valid_cenoms    = []

            for n in range(self.groups['elastic'].signals.shape[1]):
                # find the center of mass for each ROI
                cofm = xrs_utilities.find_center_of_mass(self.groups['elastic'].energy,self.groups['elastic'].signals[:,n])
                self.cenom.append(cofm)
                if self.there_is_a_valid_roi_at( n ):
                    valid_cenoms.append(cofm)
                    if Origin is None:
                        Origin = cofm

                # find the FWHM/resolution for each ROI
                en_scale  = (self.groups['elastic'].energy - self.cenom[n])*1e3
                intensity = self.groups['elastic'].signals_orig[:,n]
                FWHM, x0 = xrs_utilities.fwhm(en_scale,intensity)

            # find master E0
            self.E0 = np.mean(valid_cenoms)

            # define last eloss scale as the 'master' scale for all ROIs
            self.eloss = (self.energy - cofm)*1e3 # energy loss in eV

            # define eloss-scale for each ROI and interpolate onto the 'master' eloss-scale 
            for n in range(self.signals.shape[1]):
                # inserting zeros at beginning and end of the vectors to avoid interpolation errors
                x = (self.energy-self.cenom[n])*1e3
                x = np.insert(x,0,-1e10)
                x = np.insert(x,-1,1e10)
                y = self.signals[:,n]
                y = np.insert(y,0,0)
                y = np.insert(y,-1,0)
                f = interp1d(x,y, bounds_error=False,fill_value=0.0)
                self.signals[:,n] = f(self.eloss)

            # do the same for the errors 
            for n in range(self.signals.shape[1]):
                # inserting zeros at beginning and end of the vectors to avoid interpolation errors
                x = (self.energy-self.cenom[n])*1e3
                x = np.insert(x,0,-1e10)
                x = np.insert(x,-1,1e10)
                y = self.errors[:,n]
                y = np.insert(y,0,0)
                y = np.insert(y,-1,0)
                f = interp1d(x,y, bounds_error=False,fill_value=0.0)
                self.errors[:,n] = f(self.eloss)

    def get_spectrum( self, include_elastic=False, abs_counts=False ):
        """ **get_spectrum**

        Constructs a spectrum based on the scans loaded so far. Defines the energy 
        loss scale based on the elastic lines.

        Args:
            include_elastic (boolean): Boolean flag, does not include the elastic line if
                set to 'False' (this is the default).
            abs_counts      (boolean): Boolean flag, constructs the spectrum in absolute
                counts if set to 'True' (default is 'False')

        """

        # find the groups 
        all_groups = xrs_scans.findgroups(self.scans)

        # append scans
        for group in all_groups:
            one_group = xrs_scans.makegroup_nointerp(group, absCounts=abs_counts)
            self.groups[one_group.get_type()] = one_group

        self.energy, self.signals, self.errors = xrs_scans.appendScans(self.groups,include_elastic)

        # define the energy loss scale
        self.get_eloss()

    def get_spectrum_new( self, method='sum', include_elastic=False, abs_counts=False, interpolation=False ):
        """ **get_spectrum_new**

        Constructs a spectrum from the scans loaded so far based on the chosen method: 
        detector pixels can either be summed up (method='sum'), pixel-by-pixel 
        compensation can be applied (method='pixel'), or a line-by-line compensation
        scheme can be applied (method='row').

        The energy loss scale will be defined in the process and the data will be
        interpolated onto this scale using norm-conserving wavelet interpolation.

        Args:
            method              (str): Keyword describing the kind of integration scheme
                to be used (possible values are 'sum', 'pixel', or 'row'), default is 'sum'.
            include_elastic (boolean): Boolean flag, does not include the elastic line if
                set to 'False' (this is the default).
            abs_counts      (boolean): Boolean flag, constructs the spectrum in absolute
                counts if set to 'True' (default is 'False')
            interpolation   (boolean): Boolean flag, if True, signals are interpolated 
                onto energy grid of the first scan in each group of scans.

        """

        if not method in ['pixel', 'sum', 'row']:
            print('Unknown integration method. Use either \'pixel\', \'sum\', or \'row\'')
            return

        # make sure there is an elastic line available
        elastic_number = None
        for key in self.scans:
            if self.scans[key].scan_type == 'elastic':
                elastic_number = self.scans[key].scan_number
                print ('Using scan No. %d for CENOMs.'%elastic_number)
                break
            else:
                pass
            
        if not elastic_number:
            print( 'Please load/integrate at least one elastic scan first!' )
            return

        # get compensation factor/CENOM for an elastic line
        if method in [ 'sum' , 'pixel']  and not self.cenom_dict:
            self.get_compensation_factor( elastic_number, method=method )

        if method == 'row' and not self.comp_factor:

            self.get_compensation_factor( elastic_number, method=method )

        # get raw data
        for key in self.scans:
            if not self.scans[key].raw_signals:
                self.scans[key].get_raw_signals( self.roi_obj, method=method )

        # sum up similar scans
        # find all groups of scans
        all_groups = xrs_scans.findgroups(self.scans)

        # initiate groups
        self.groups = {}

        # sum up similar scans
        for group in all_groups:
            self.groups[group[0].get_type()] = xrs_scans.sum_scans_to_group( group, method=method, interp=interpolation )

        # stitch groups together into a spectrum
        spectrum = xrs_scans.stitch_groups_to_spectrum( self.groups, method=method, include_elastic=include_elastic )
        self.energy      = spectrum.energy
        self.signals     = spectrum.signals
        self.errors      = spectrum.errors
        self.raw_signals = spectrum.raw_signals
        self.raw_errors  = spectrum.raw_errors

        # define energy loss scale and apply compensation if applicable
        self.get_eloss_new( method=method )

    def get_q_values( self, inv_angstr=False, energy_loss=None ):
        """ **get_q_values**

        Calculates the momentum transfer for each analyzer.

        Args:
            inv_angstr (boolean): Boolean flag, if 'True' momentum transfers are calculated in 
                inverse Angstroms.
            energy_loss  (float): Energy loss value at which the momentum transfer is to be 
                calculated. If 'None' is given, the momentum transfer is calculated for every 
                energy loss point of the spectrum.

        Returns:
            If an energy loss value is passed, the function returns the momentum transfers
            at this energy loss value for each analyzer crystal.

        """

        # one q-value per analyzer and energy loss point
        q_vals = np.zeros_like(self.signals)
        if inv_angstr:
            for n in range( self.signals.shape[1] ):
                q_vals[:,n] = xrs_utilities.momtrans_inva(self.E0+self.eloss/1e3,self.E0,self.tth[n]) 
        else:
            for n in range( self.signals.shape[1] ):
                q_vals[:,n] = xrs_utilities.momtrans_au(self.E0+self.eloss/1e3,self.E0,self.tth[n])
        self.q_values = q_vals

        # return all q-values if a specific energy loss is given
        if energy_loss:
            ind = np.abs(self.eloss - energy_loss).argmin()
            return self.q_values[ind,:]

    def copy_edf_files( self, scan_numbers, dest_dir ):
        """ **copy_edf_files**

        Copies all EDF-files from given scan_numbers into given directory.

        Args:
          * scan_numbers (int or list) = Integer or list of integers defining the scan 
            numbers of scans to be copied.
          * dest_dir             (str) = String with absolute path for the destination.

        """

        import shutil

        # make scan_numbers iterable
        numbers = []
        if not isinstance(scan_numbers,list):
            numbers.append(scan_numbers)
        else:
            numbers = scan_numbers
        fname = self.path + self.SPECfname

        # go through the scans, find the EDF-files and copy them
        for nscan in numbers:
            if use_PyMca:
                data, motors, counters = xrs_fileIO.PyMcaSpecRead(fname,nscan)
            else:
                data, motors, counters = xrs_fileIO.SpecRead(fname,nscan)

            for m in range(len(counters['ccdno'])):
                ccdnumber = counters['ccdno'][m]
                edfname   = self.path + self.EDFprefix + self.EDFname + "%04d" % ccdnumber + self.EDFpostfix
                shutil.copy2(edfname, dest_dir)

    def dump_spectrum_ascii( self, file_name, header='' ):
        """ **dump_spectrum_ascii**

        Stores the energy loss and signals in a txt-file.

        Args:
            filename (str): Path and filename to the file to be written.

        """

        data = np.zeros((len(self.eloss),self.signals.shape[1]*2+1))
        data[:,0] = self.eloss
        col = 1
        for ii in range(self.signals.shape[1]):
            data[:,col] = self.signals[:,ii]
            col += 1
            data[:,col] = self.errors[:,ii]
            col += 1

        np.savetxt( file_name, data, header=header )

    def dump_spectrum_hdf5( self, file_name, group_name, comment='' ):
        """ **dump_spectrum_hdf5**

        Writes the summed spectrum into an HDF5 file.

        Args:
            file_name  (str): Path and file name for the HDF5-file to be created.
            group_name (str): Group name under which to store status in the HDF5-file.
            comment    (str): Optional comment (no comment is default).

        """

        h5 = h5py.File(file_name,"a")
        h5.require_group(group_name)
        h5group =  h5[group_name]

        keys = ['energy', 'eloss', 'signals', 'errors']
        for key in keys:
            h5group[key] = getattr( self, key )

        h5group["comment"]  = comment
        h5.flush()
        h5.close()

    def dump_scans_ascii( self, scan_numbers, pre_fix, f_name, post_fix='.dat', header='' ):
        """ **dump_scans_ascii**

        Produce ASCII-type files with columns of energy, signal, and Poisson error.

        Args:
            scan_numbers (int or list): SPEC scan numbers of scans to be safed in ASCII format.
            pre_fix  (str): Path to directory where files should be written into.
            f_name   (str): Base name for the files to be written.
            post_fix (str): Extention for the files to be written (default is '.dat').

        """

        # make sure scan_numbers are iterable
        numbers = []
        if not isinstance(scan_numbers,list):
            numbers.append(scan_numbers)
        else:
            numbers = scan_numbers

        for number in numbers:
            scan_name = 'Scan%03d' % number
            if not scan_name in self.scans.keys():
                print ('Scan No. %03d is currently not loaded.'% number)
                return
            x = self.scans[scan_name].energy
            y = self.scans[scan_name].signals
            z = self.scans[scan_name].errors
            data = np.zeros((x.shape[0], y.shape[1]+z.shape[1]+1 ))
            data[:,0] = x
            data[:,1:1+y.shape[1]] = y
            data[:,1+y.shape[1]:1+y.shape[1]+z.shape[1]] = z
            file_name = pre_fix + f_name + '_' + scan_name + post_fix
            np.savetxt(file_name, data, header=header)

    def print_scan_length( self,scan_numbers ):
        """ **print_scan_length**

        Print out the numper of points in given scans.

        Args:
            scan_numbers (int or list): Scan number or list of scan numbers.

        """

        # make scan_numbers iterable
        numbers = []
        if not isinstance(scan_numbers,list):
            numbers.append(scan_numbers)
        else:
            numbers = scan_numbers

        # print out the number of points for all scan_numbers
        for ii in numbers:
            name = 'Scan%03d' % ii
            print( 'Length of scan %03d ' %ii + ' is ' + str(len(self.scans[name].energy)) + '.')

    def get_tths( self, rvd=None, rvu=None, rvb=None, rhl=None, rhr=None, rhb=None, order=[0,1,2,3,4,5] ):
        """ **get_tths**

        Calculates the scattering angles for all analyzer crystals based on 
        the mean angle of the analyzer modules.

        Args:
              * rhl  (float): Mean scattering angle of the HL module (default is 0.0).
              * rhr  (float): Mean scattering angle of the HR module (default is 0.0).
              * rhb  (float): Mean scattering angle of the HB module (default is 0.0).
              * rvd  (float): Mean scattering angle of the VD module (default is 0.0).
              * rvu  (float): Mean scattering angle of the VU module (default is 0.0).
              * rvb  (float): Mean scattering angle of the VB module (default is 0.0).
              * order (list): List of integers (0-5) that describe the order of modules in which the 
                ROIs were defined (default is VD, VU, VB, HR, HL, HB; i.e. [0,1,2,3,4,5]).

        """
        # reset all tth values
        self.tth   = []

        # try to grab from positions from SPEC-file (first scan in list)
        scan_name = list(self.scans.keys())[0]
        if not rvd:
            rvd = self.scans[scan_name].motors['RVD']
        if not rvu:
            rvu = self.scans[scan_name].motors['RVU']
        if not rvb:
            rvb = self.scans[scan_name].motors['RVB']
        if not rhr:
            rhr = self.scans[scan_name].motors['RHR']
        if not rhl:
            rhl = self.scans[scan_name].motors['RHL']
        if not rhb:
            rhb = self.scans[scan_name].motors['RHB']

        # horizontal modules
        # HL (motor name rhl)
        v_angles = self.TTH_OFFSETS1 
        h_angles = self.TTH_OFFSETS2 + rhl
        HLtth    = []
        for n in range(len(h_angles)):
            HLtth.append( np.arccos(np.cos(np.radians(h_angles[n]))*np.cos(np.radians(v_angles[n])))*180.0/np.pi)

        # HR (motor name rhr)
        v_angles = self.TTH_OFFSETS1 
        h_angles = self.TTH_OFFSETS2 + rhr
        HRtth    = []
        for n in range(len(h_angles)):
            HRtth.append( np.arccos(np.cos(np.radians(h_angles[n]))*np.cos(np.radians(v_angles[n])))*180.0/np.pi)

        # HB (motor name rhb)
        v_angles = self.TTH_OFFSETS1 
        h_angles = self.TTH_OFFSETS2 + rhb
        HBtth = []
        for n in range(len(h_angles)):
            HBtth.append( np.arccos(np.cos(np.radians(h_angles[n]))*np.cos(np.radians(v_angles[n])))*180.0/np.pi)

        # vertical modules
        # VD
        v_angles = self.TTH_OFFSETS2 + rvd
        h_angles = self.TTH_OFFSETS1 
        VDtth    = []
        for n in range(len(h_angles)):
            VDtth.append( np.arccos(np.cos(np.radians(v_angles[n]))*np.cos(np.radians(h_angles[n])))*180.0/np.pi)

        # VU
        v_angles = self.TTH_OFFSETS2 + rvu
        h_angles = self.TTH_OFFSETS1 
        VUtth    = []
        for n in range(len(h_angles)):
            VUtth.append( np.arccos(np.cos(np.radians(v_angles[n]))*np.cos(np.radians(h_angles[n])))*180.0/np.pi)