<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 In this notebook, we will analyze a NuSTAR data of the black hole X-ray binary H1743-322 with Stingray. Here we assume that the user has already reduced the data with the official pipeline and ran barycorr or other tools to refer the event times to the solar system barycenter.
Step1:  Quicklook NuSTAR data with Stingray
Let us load the data from two event lists corresponding to the two detectors onboard NuSTAR. fmt='hea' indicates event data produced by HEAsoft tools or compatible (e.g. XMM-Newton).
Step2:  For the sake of a quicklook, let us join the two event lists
Step3:  Let us calculate the light curve and plot it. 
In red, we show the bad time intervals, when the satellite was not acquiring valid data due to Earth occultation, SAA passages, or other issues.
Step4:  The light curve shows some long-term variability. Let us look at the colors. First of all, let us check that the events contain the energy of each photon. This should be the case, because NuSTAR data, together with XMM and NICER, are very well understood by Stingray and the calibration is done straight away.
Step5:  Other missions might have all_ev.energy set to None. In which case, one needs to use all_ev.pi and express the energy through the PI channels (See the HENDRICS documentation for more advanced calibration using the rmf files).
Also, we notice that some GTIs do not catch all bad intervals (see how the light curve drops close to GTI borders). We make a more aggressive GTI filtering now
Step6:  Despite some light curve variability, the hardness ratio seems pretty stable during the observation.
Let us now look at the power density spectrum. Notice that we are using a sampling time of 0.001 s, meaning that we will investigate the power spectrum up to 500 Hz
Step7:  Nice Quasi-periodic oscillations there! Note that at high frequencies the white noise level increases. This is not real variability, but an effect of dead time. The easiest way to get a flat periodogram at high frequencies is using the cospectrum instead of the power density spectrum. For this, we use separately the events from the two detectors. The cospectrum calculation is slightly slower than the power spectrum.
For an accurate way to correct the power density spectrum from dead time, see the documentation of stingray.deadtime and the Frequency Amplitude Difference (FAD) correction.
Step8:  To improve the plot, we can rebin the data logarithmically
<ASSISTANT_TASK:>
Python Code:
%load_ext autoreload
%autoreload 2
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from stingray.powerspectrum import AveragedPowerspectrum, DynamicalPowerspectrum
from stingray.crossspectrum import AveragedCrossspectrum
from stingray.events import EventList
from stingray.lightcurve import Lightcurve
from stingray.gti import create_gti_from_condition
Explanation: In this notebook, we will analyze a NuSTAR data of the black hole X-ray binary H1743-322 with Stingray. Here we assume that the user has already reduced the data with the official pipeline and ran barycorr or other tools to refer the event times to the solar system barycenter.
End of explanation
evA = EventList.read('nustar_A_src.evt', 'hea')
evB = EventList.read('nustar_B_src.evt', 'hea')
Explanation: Quicklook NuSTAR data with Stingray
Let us load the data from two event lists corresponding to the two detectors onboard NuSTAR. fmt='hea' indicates event data produced by HEAsoft tools or compatible (e.g. XMM-Newton).
End of explanation
all_ev = evA.join(evB)
Explanation: For the sake of a quicklook, let us join the two event lists
End of explanation
lc = all_ev.to_lc(100)
plt.figure(figsize=(12, 7))
plt.plot(lc.time, lc.counts)
bad_time_intervals = list(zip(lc.gti[:-1, 1], lc.gti[1:, 0]))
for b in bad_time_intervals:
    plt.axvspan(b[0], b[1], color='r', alpha=0.5, zorder=10)
plt.ylim([5000, 6500])
Explanation: Let us calculate the light curve and plot it. 
In red, we show the bad time intervals, when the satellite was not acquiring valid data due to Earth occultation, SAA passages, or other issues.
End of explanation
all_ev.energy
Explanation: The light curve shows some long-term variability. Let us look at the colors. First of all, let us check that the events contain the energy of each photon. This should be the case, because NuSTAR data, together with XMM and NICER, are very well understood by Stingray and the calibration is done straight away.
End of explanation
new_gti = create_gti_from_condition(lc.time, lc.counts > 5200)
all_ev.gti = new_gti
evA.gti = new_gti
evB.gti = new_gti
lc.gti = new_gti
hard = (all_ev.energy > 10) & (all_ev.energy < 79)
soft = (all_ev.energy > 3) & (all_ev.energy < 5)
hard_ev = all_ev.apply_mask(hard)
soft_ev = all_ev.apply_mask(soft)
hard_lc = hard_ev.to_lc(200)
soft_lc = soft_ev.to_lc(200)
hard_lc.apply_gtis()
soft_lc.apply_gtis()
hardness_ratio = hard_lc.counts / soft_lc.counts
intensity = hard_lc.counts + soft_lc.counts
plt.figure()
plt.scatter(hardness_ratio, intensity)
plt.xlabel("Hardness")
plt.ylabel("Counts")
Explanation: Other missions might have all_ev.energy set to None. In which case, one needs to use all_ev.pi and express the energy through the PI channels (See the HENDRICS documentation for more advanced calibration using the rmf files).
Also, we notice that some GTIs do not catch all bad intervals (see how the light curve drops close to GTI borders). We make a more aggressive GTI filtering now:
End of explanation
pds = AveragedPowerspectrum.from_events(all_ev, segment_size=256, dt=0.001, norm='leahy')
plt.figure(figsize=(10,7))
plt.loglog(pds.freq, pds.power)
plt.xlabel("Frequency")
plt.ylabel("Power (Leahy)")
Explanation: Despite some light curve variability, the hardness ratio seems pretty stable during the observation.
Let us now look at the power density spectrum. Notice that we are using a sampling time of 0.001 s, meaning that we will investigate the power spectrum up to 500 Hz
End of explanation
cs = AveragedCrossspectrum.from_events(evA, evB, segment_size=256, dt=0.001, norm='leahy')
plt.figure(figsize=(10,7))
plt.semilogx(cs.freq, cs.power.real)
Explanation: Nice Quasi-periodic oscillations there! Note that at high frequencies the white noise level increases. This is not real variability, but an effect of dead time. The easiest way to get a flat periodogram at high frequencies is using the cospectrum instead of the power density spectrum. For this, we use separately the events from the two detectors. The cospectrum calculation is slightly slower than the power spectrum.
For an accurate way to correct the power density spectrum from dead time, see the documentation of stingray.deadtime and the Frequency Amplitude Difference (FAD) correction.
End of explanation
cs_reb = cs.rebin_log(0.02)
plt.figure(figsize=(10,7))
plt.loglog(cs_reb.freq, cs_reb.power.real)
plt.ylim([1e-3, None])
plt.xlabel("Frequency")
plt.ylabel("Cospectrum Power")
Explanation: To improve the plot, we can rebin the data logarithmically
End of explanation
<END_TASK>