<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Reading µs-ALEX data from Photon-HDF5 with h5py
In this notebook we show how to read a µs-ALEX smFRET measurement stored in
Photon-HDF5 format
using python and a few common scientific libraries (numpy, h5py, matplotlib).
Specifically, we show how to load timestamps, build an alternation histogram
and select photons in the donor and acceptor excitation periods.
See also a similar notebook using 
pytables instead of h5py.
Step2:  1. Utility functions
Here we define an utility function to print HDF5 file contents
Step3:  2. Open the data file
Let assume we have a Photon-HDF5 file at the following location
Step4:  We can open the file, as a normal HDF5 file
Step5:  The object h5file is a pytables file reference. The root group is accessed with h5file.root.
3. Print the content
Let's start by taking a look at the file content
Step6:  We see the typical Photon-HDF5 structure. In particular the field description provides a short description of the measurement and acquisition_duration tells that the acquisition lasted 600 seconds.
As an example let's take a look at the content of the sample group
Step7:  Finally, we define a shortcut to the photon_data group to save some typing later
Step8:  4. Reading the data
First, we make sure the file contains the right type of measurement
Step9:  Ok, tha's what we espect. 
Now we can load all the timestamps (including timestamps unit) and detectors arrays
Step10:  We may want to check the excitation wavelengths used in the measurement. This information is found in the setup group
Step11:  Now, let's load the definitions of donor/acceptor channel and excitation periods
Step12:  These numbers define the donor and acceptor alternation periods as shown below
Step13:  Now that the data has been loaded we can plot an alternation histogram using matplotlib
Step14:  6. Timestamps in different excitation periods
We conclude by showing, as an example, how to create arrays of timestamps containing only donor or acceptor exitation photons.
<ASSISTANT_TASK:>
Python Code:
from __future__ import division, print_function  # only needed on py2
%matplotlib inline
import numpy as np
import h5py
import matplotlib.pyplot as plt
Explanation: Reading µs-ALEX data from Photon-HDF5 with h5py
In this notebook we show how to read a µs-ALEX smFRET measurement stored in
Photon-HDF5 format
using python and a few common scientific libraries (numpy, h5py, matplotlib).
Specifically, we show how to load timestamps, build an alternation histogram
and select photons in the donor and acceptor excitation periods.
See also a similar notebook using 
pytables instead of h5py.
End of explanation
def print_children(group):
    Print all the sub-groups in `group` and leaf-nodes children of `group`.
    Parameters:
        data_file (h5py HDF5 file object): the data file to print
    
    for name, value in group.items():
        if isinstance(value, h5py.Group):
            content = '(Group)'
        else:
            content = value[()]
        print(name)
        print('    Content:     %s' % content)
        print('    Description: %s\n' % value.attrs['TITLE'].decode())
Explanation: 1. Utility functions
Here we define an utility function to print HDF5 file contents:
End of explanation
filename = '../data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'
Explanation: 2. Open the data file
Let assume we have a Photon-HDF5 file at the following location:
End of explanation
h5file = h5py.File(filename)
Explanation: We can open the file, as a normal HDF5 file
End of explanation
print_children(h5file)
Explanation: The object h5file is a pytables file reference. The root group is accessed with h5file.root.
3. Print the content
Let's start by taking a look at the file content:
End of explanation
print_children(h5file['sample'])
Explanation: We see the typical Photon-HDF5 structure. In particular the field description provides a short description of the measurement and acquisition_duration tells that the acquisition lasted 600 seconds.
As an example let's take a look at the content of the sample group:
End of explanation
photon_data = h5file['photon_data']
Explanation: Finally, we define a shortcut to the photon_data group to save some typing later:
End of explanation
photon_data['measurement_specs']['measurement_type'][()].decode()
Explanation: 4. Reading the data
First, we make sure the file contains the right type of measurement:
End of explanation
timestamps = photon_data['timestamps'][:]
timestamps_unit = photon_data['timestamps_specs']['timestamps_unit'][()]
detectors = photon_data['detectors'][:]
print('Number of photons: %d' % timestamps.size)
print('Timestamps unit:   %.2e seconds' % timestamps_unit)
print('Detectors:         %s' % np.unique(detectors))
Explanation: Ok, tha's what we espect. 
Now we can load all the timestamps (including timestamps unit) and detectors arrays:
End of explanation
h5file['setup']['excitation_wavelengths'][:]
Explanation: We may want to check the excitation wavelengths used in the measurement. This information is found in the setup group:
End of explanation
donor_ch = photon_data['measurement_specs']['detectors_specs']['spectral_ch1'][()]
acceptor_ch = photon_data['measurement_specs']['detectors_specs']['spectral_ch2'][()]
print('Donor CH: %d     Acceptor CH: %d' % (donor_ch, acceptor_ch))
alex_period = photon_data['measurement_specs']['alex_period'][()]
donor_period = photon_data['measurement_specs']['alex_excitation_period1'][()]
offset = photon_data['measurement_specs']['alex_offset'][()]
acceptor_period = photon_data['measurement_specs']['alex_excitation_period2'][()]
print('ALEX period:     %d  \nOffset:         %4d      \nDonor period:    %s      \nAcceptor period: %s' % \
      (alex_period, offset, donor_period, acceptor_period))
Explanation: Now, let's load the definitions of donor/acceptor channel and excitation periods:
End of explanation
timestamps_donor = timestamps[detectors == donor_ch]
timestamps_acceptor = timestamps[detectors == acceptor_ch]
Explanation: These numbers define the donor and acceptor alternation periods as shown below:
$$2180 < \widetilde{t} < 3900 \qquad \textrm{donor period}$$
$$200 < \widetilde{t} < 1800 \qquad \textrm{acceptor period}$$
where $\widetilde{t}$ represent the (timestamps - offset) MODULO alex_period.
For more information
please refer to the measurements_specs section
of the Reference Documentation.
5. Plotting the alternation histogram
Let start by separating timestamps from donor and acceptor channels:
End of explanation
fig, ax = plt.subplots()
ax.hist((timestamps_acceptor - offset) % alex_period, bins=100, alpha=0.8, color='red', label='donor')
ax.hist((timestamps_donor - offset) % alex_period, bins=100, alpha=0.8, color='green', label='acceptor')
ax.axvspan(donor_period[0], donor_period[1], alpha=0.3, color='green')
ax.axvspan(acceptor_period[0], acceptor_period[1], alpha=0.3, color='red')
ax.set_xlabel('(timestamps - offset) MOD alex_period')
ax.set_title('ALEX histogram')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False);
Explanation: Now that the data has been loaded we can plot an alternation histogram using matplotlib:
End of explanation
timestamps_mod = (timestamps - offset) % alex_period
donor_excitation = (timestamps_mod < donor_period[1])*(timestamps_mod > donor_period[0])
acceptor_excitation = (timestamps_mod < acceptor_period[1])*(timestamps_mod > acceptor_period[0])
timestamps_Dex = timestamps[donor_excitation]
timestamps_Aex = timestamps[acceptor_excitation]
fig, ax = plt.subplots()
ax.hist((timestamps_Dex - offset) % alex_period, bins=np.arange(0, alex_period, 40), alpha=0.8, color='green', label='D_ex')
ax.hist((timestamps_Aex - offset) % alex_period, bins=np.arange(0, alex_period, 40), alpha=0.8, color='red', label='A_ex')
ax.set_xlabel('(timestamps - offset) MOD alex_period')
ax.set_title('ALEX histogram (selected periods only)')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False);
#plt.close('all')
Explanation: 6. Timestamps in different excitation periods
We conclude by showing, as an example, how to create arrays of timestamps containing only donor or acceptor exitation photons.
End of explanation
<END_TASK>