Power spectrum of real data with fftpack on log axis

I already read many discussion about this topic (comparison between lomb-scargle and fft , Plotting power spectrum in python, Scipy/Numpy FFT Frequency Analysis, and many others), but still can't manage it, so I need some tips. I have a list of photon events (detections vs time), the data are available here. The columns are time, counts , errors, and counts in different energy bands (you can ignore them). I know the source has a periodicity around 8.9 days = 1.3*10^-6 Hz. I would like to plot the Power spectrum density showing a peak at this frequency (on a log x-axis, possibly). It would also be nice if I can avoid the half part of the plot (symmetric). This is my code till now, not so far but still something:

import numpy as np
from scipy.fftpack import fft, rfft, fftfreq
import pylab as plt

x,y = np.loadtxt('datafile.txt', usecols = (0,1), unpack=True)
y = y - y.mean() # Removes the large value at the 0 frequency that we don't care about

f_range = np.linspace(10**(-7), 10**(-5), 1000)
W = fftfreq(y.size, d=x[1]-x[0])

plt.subplot(2,1,1)
plt.plot(x,y)
plt.xlabel('Time (days)')

f_signal = fft(y)
plt.subplot(2,1,2)
plt.plot(W, abs(f_signal))
plt.xlabel('Frequency (Hz)')

Here the (useless) plot produced: DFT

Solution

Here is an improved version of the code above:

import pyfits
import numpy as np
from scipy.fftpack import fft, rfft, fftfreq
import pylab as plt

x,y = np.loadtxt('data.txt', usecols = (0,1), unpack=True)
y = y - y.mean()

W = fftfreq(y.size, d=(x[1]-x[0])*86400)
plt.subplot(2,1,1)
plt.plot(x,y)
plt.xlabel('Time (days)')

f_signal = fft(y)
plt.subplot(2,1,2)
plt.plot(W, abs(f_signal)**2)
plt.xlabel('Frequency (Hz)')

plt.xscale('log')
plt.xlim(10**(-6), 10**(-5))
plt.show()

And here the plot produced (correctly): fft The highest peak is the peak I was trying to reproduce. The second peak is also expected, but with less power (as it is, indeed). If rfft is used instead of fft (and rfftfreq instead of fftfreq) the same plot is reproduced (in that case, the frequencies values, instead of the module, can be used numpy.fft.rfft)

I don't want to block the topic, so I will ask here: And how can I retrieve the frequencies of the peaks? Would be great to plot the frequencies by side the peaks.