Timing tools - Revision history

Grinberg at 11:13, 15 April 2021

2021-04-15T11:13:28Z

Grinberg at 11:11, 15 April 2021

2021-04-15T11:11:34Z

Koenig at 10:26, 21 December 2020

2020-12-21T10:26:54Z

Gokus: /* Lag-energy spectrum for archival XMM-Newton data */

2020-08-11T15:14:17Z

Lag-energy spectrum for archival XMM-Newton data

Gokus: /* Lag-energy spectrum for archival XMM-Newton data */

2020-08-11T15:11:49Z

Lag-energy spectrum for archival XMM-Newton data

Gokus: /* Lag-energy spectrum for archival XMM-Newton data */

2020-08-10T15:19:18Z

Lag-energy spectrum for archival XMM-Newton data

Koenig at 11:50, 21 February 2020

2020-02-21T11:50:17Z

Koenig at 11:32, 21 February 2020

2020-02-21T11:32:19Z

Koenig at 11:21, 21 February 2020

2020-02-21T11:21:40Z

Koenig at 11:12, 21 February 2020

2020-02-21T11:12:47Z

@@ Line 21: / Line 21: @@
 === Calculating the PSD ===
-Given two PCA lightcurves, <pre>50110-01-45-00_437_AND_good_14off_excl_0-10.lc</pre> and <pre>50110-01-45-00_437_AND_good_14off_excl_11-13.lc</pre> (which have been checked for gaps and found to be gapless <ref>A useful function for this is <tt>split_lc_at_gaps</tt></ref>. Alternatively - but please read carefully what the function does - there is the function <tt>segment_lc_for_psd<tt> that does the job for you and prepares your lightcurve for PSD calculation the way foucalc does it.) we want to calculate and bin the PSD, lags and coherence.
+Given two PCA lightcurves, <pre>50110-01-45-00_437_AND_good_14off_excl_0-10.lc</pre> and <pre>50110-01-45-00_437_AND_good_14off_excl_11-13.lc</pre> (which have been checked for gaps and found to be gapless <ref>A useful function for this is <tt>split_lc_at_gaps</tt>. Alternatively - but please read carefully what the function does - there is the function <tt>segment_lc_for_psd</tt> that does the job for you and prepares your lightcurve for PSD calculation the way foucalc does it.</ref>) we want to calculate and bin the PSD, lags and coherence.
 <pre>

@@ Line 21: / Line 21: @@
 === Calculating the PSD ===
-Given two PCA lightcurves, <pre>50110-01-45-00_437_AND_good_14off_excl_0-10.lc</pre> and <pre>50110-01-45-00_437_AND_good_14off_excl_11-13.lc</pre> (which have been checked for gaps and found to be gapless <ref>A useful function for this is <tt>split_lc_at_gaps</tt></ref> ) we want to calculate and bin the PSD, lags and coherence.
+Given two PCA lightcurves, <pre>50110-01-45-00_437_AND_good_14off_excl_0-10.lc</pre> and <pre>50110-01-45-00_437_AND_good_14off_excl_11-13.lc</pre> (which have been checked for gaps and found to be gapless <ref>A useful function for this is <tt>split_lc_at_gaps</tt></ref>. Alternatively - but please read carefully what the function does - there is the function <tt>segment_lc_for_psd<tt> that does the job for you and prepares your lightcurve for PSD calculation the way foucalc does it.) we want to calculate and bin the PSD, lags and coherence.
 <pre>

← Older revision		Revision as of 10:26, 21 December 2020
Line 216:		Line 216:

	The coherence function is 1 at low frequencies showing a perfect (linear) correlation of the two lightcurves (values above 1 are possible due to counting noise, see Sect. 3.4.1 and 3.5.3 in Katja's diss.). As a final step, we calculate the lag-energy spectrum, by extracting lightcurves in logarithmically binned energy intervals, and calculating their timelags.		The coherence function is 1 at low frequencies showing a perfect (linear) correlation of the two lightcurves (values above 1 are possible due to counting noise, see Sect. 3.4.1 and 3.5.3 in Katja's diss.). As a final step, we calculate the lag-energy spectrum, by extracting lightcurves in logarithmically binned energy intervals, and calculating their timelags.
−
−	~~<pre>~~
−	~~refeband=[540,10080]; % (eV), lower and upper bound of energy reference band~~
−	~~reffreqband=[0.034,0.12]; %% (Hz), lower and upper bound of frequency reference band, as in Fig. 2(a) of Uttley+11~~
−
−	~~fnameref=sprintf("%sfull_sd_0.1s_%04d-%d.lc",datadir,refeband[0],refeband[1]);~~
−	~~lcref = fits_read_table(fnameref);~~
−	~~lcrefgti = fits_read_table(fnameref+"[GTI00006]"); % GTI correction~~
−	~~filtref = filter_gti(lcref.time, lcrefgti);~~
−	~~timelcref = lcref.time[filtref]-tstart;~~
−	~~cratelcref = lcref.counts[filtref]/dt;~~
−
−	~~(elo,ehi)=log_grid(500,10000,14); % define logarithmic energy bins~~
−
−	~~tlags=Double_Type[length(elo)];~~
−	~~tlagerrs=Double_Type[length(elo)];~~
−
−	~~_for ii (0,length(elo)-1,1) {~~
−	~~lcname=sprintf("full_sd_0.1s_%04.0f-%04.0f.lc",elo[ii],ehi[ii]);~~
−	~~if (ehi[ii]>9999) { lcname=sprintf("full_sd_0.1s_%04.0f-%05.0f.lc",elo[ii],ehi[ii]); };~~
−
−	~~lc = fits_read_table(lcname);~~
−	~~lcgti = fits_read_table(lcname+"[GTI00006]");~~
−	~~filtgti = filter_gti(lc.time, lcgti);~~
−
−	~~timelc=lc.time[filtgti]-tstart; %% Do GTI correction~~
−	~~cratelc=lc.counts[filtgti]/dt;~~
−
−	~~%% Calculate PSD, cross-spectrum (CPD) and lags with foucalc~~
−	~~%% Note that we subtract off the small band LC from the reference~~
−	~~%% LC. Uttley+11 call this a Channel-of-interest (CI) correction.~~
−	~~%% This is needed to account for the Poisson noise which is~~
−	~~%% correlated for both LCs. See also Uttley+11 Sect 3.2.~~
−	~~r=foucalc(struct{time=timelcref, rate1=cratelcref-cratelc, rate2=cratelc}, dimseg ; normtype="Miyamoto");~~
−
−	~~%% Filter the CPD / timelag array for the given frequency range~~
−	~~filt=where(reffreqband[0]<=r.freq<=reffreqband[1]);~~
−
−	~~%% Compute average time lag as arithmetic weighted mean~~
−	~~tlags[ii]=weighted_mean(r.lag12[filt] ; err=r.errlag12[filt]);~~
−	~~tlagerrs[ii]=1./sqrt(sum(1/r.errlag12[filt]^2));~~
−	}
−
−	~~ylin;hplot_with_err(elo/1000., ehi/1000., tlags, tlagerrs);~~
−	~~</pre>~~

@@ Line 126: / Line 126: @@
 </pre>
-The time information in the light curve is not given as a grid, however, a time grid with connecting high to low bins is necessary for the filter_gti function to filter out all time bins with an unsufficient fractional exposure. Compared to the default binning width of 100s, telemetry drop-outs can last for only a few seconds, hence only datapoints with times values that fall directly in the gaps between the GTIs, are sorted out. This can be avoided by inserting a grid, which we have to create from the information in the light curve FITS-file. Something that might be important for some analyses is the information about at which position in the grid the given time data points have to be set (either beginning, middle or the end). This is instrument dependent, and should be usually found either in the header or the manual.
+The time information in the light curve is not given as a grid, however, a time grid with connecting high to low bins is necessary for the filter_gti function to filter out all time bins with an unsufficient fractional exposure. Compared to the default binning width of 100s, telemetry drop-outs can last for only a few seconds, hence only datapoints with times values that fall directly in the gaps between the GTIs, are sorted out. This can be avoided by inserting a grid, which we have to create from the information in the light curve FITS-file. Something that might be important for some analyses is the information about at which position in the grid the given time data points have to be set (either beginning, middle or the end). This is instrument dependent, and should be usually found either in the header or the manual. For XMM data, the data points mark the beginning of each bin (trel = 0).
 <pre>
 %% The binning width is encoded in the header with the keyword TIMEDEL
 variable tdiff = fits_read_key("full_sd_0.1s_0500-0900.lc[1]", "TIMEDEL");
-variable trel = 0.5;
+variable trel = 0.;
 variable tlo_a = a.time - trel*tdiff;

@@ Line 134: / Line 134: @@
 variable tlo_a = a.time - trel*tdiff;
-variable thi_a = a.time + trel*tdiff;
+variable thi_a = a.time + (1.-trel)*tdiff;
 variable tlo_b = b.time - trel*tdiff;
-variable thi_b = b.time + trel*tdiff;
+variable thi_b = b.time + (1.-trel)*tdiff;

@@ Line 90: / Line 90: @@
 == Lag-energy spectrum for archival XMM-Newton data ==
-[--- under construction --- written by OK, please ask if something is unclear!]
+[--- under construction --- written by OK, please ask if something is unclear!; last edit by AG (10 Aug 2020) regarding the correction of telemetry drop-outs]
 The following is an example how to compute PSDs, a lag-frequency spectrum, and a lag-energy spectrum on archival XMM-Newton data of GX 339-4 <ref>The corresponding publication is [https://ui.adsabs.harvard.edu/abs/2011MNRAS.414L..60U/abstract Uttley, P. e al, MNRAS, Vol. 414, I. 1, pp. L60-L64 (2011)]</ref>.
@@ Line 110: / Line 110: @@
 </pre>
-Now we can start the analysis. Let's load the data and correct them for the telemetry drop-outs (encoded in the Good Time Intervals):
+Now we can start the analysis. Let's load the data and the Good Time Intervals (GTI) information to be able to correct the light curve for the telemetry drop-outs:
 <pre>
@@ Line 124: / Line 124: @@
 tstart=fits_read_key("full_sd_0.1s_0500-0900.lc[1]", "TSTART");
 %% Subtract TSTART and do GTI correction
-(filta,filtb)     = filter_gti(a.time, gtia), filter_gti(b.time,gtib);
+variable minimum_frac = 0.95;
 (timea,timeb)     = a.time[filta]-tstart,b.time[filtb]-tstart;
 (countsa,countsb) = a.counts[filta],b.counts[filtb];

@@ Line 207: / Line 207: @@
 filtref     = filter_gti(lcref.time, lcrefgti);
 timelcref   = lcref.time[filtref]-tstart;
-countslcref = lcref.counts[filtref]/dt;
+cratelcref = lcref.counts[filtref]/dt;
-(lo,hi)=log_grid(500,10000,14); % define logarithmic energy bins
+(elo,ehi)=log_grid(500,10000,14); % define logarithmic energy bins
 tlags=Double_Type[length(elo)];
@@ Line 215: / Line 215: @@
 _for ii (0,length(elo)-1,1) {
-   lcname=datadir+sprintf("full_sd_0.1s_%04.0f-%04.0f.lc",elo[ii],ehi[ii]);
+   lcname=sprintf("full_sd_0.1s_%04.0f-%04.0f.lc",elo[ii],ehi[ii]);
-   if (ehi[ii]>9999) { lcname=datadir+sprintf("full_sd_0.1s_%04.0f-%05.0f.lc",elo[ii],ehi[ii]); };
+   if (ehi[ii]>9999) { lcname=sprintf("full_sd_0.1s_%04.0f-%05.0f.lc",elo[ii],ehi[ii]); };
    lc = fits_read_table(lcname);
@@ Line 223: / Line 223: @@
    timelc=lc.time[filtgti]-tstart; %% Do GTI correction
-   countslc=lc.counts[filtgti]/dt;
+   cratelc=lc.counts[filtgti]/dt;
    %% Calculate PSD, cross-spectrum (CPD) and lags with foucalc
@@ Line 238: / Line 238: @@
    tlags[ii]=weighted_mean(r.lag12[filt] ; err=r.errlag12[filt]);
    tlagerrs[ii]=1./sqrt(sum(1/r.errlag12[filt]^2));
-hplot(elo/1000., ehi/1000., tlags, tlagerrs);
+ylin;hplot_with_err(elo/1000., ehi/1000., tlags, tlagerrs);
 </pre>

@@ Line 140: / Line 140: @@
 <pre>
 %% Choose segment length
-dimseg=8192; %% (bins), should be a power of 2 for FFT
+dimseg=8192; % (bins), should be a power of 2 for FFT
 %% Compute PSDs, CPD, and time lags with foucalc
 r=foucalc( struct{time=timea, rate1=countsa/dt, rate2=countsb/dt}, dimseg ; normtype="Miyamoto");
-xlog;ylog;plot_with_err(r.freq,r.signormpsd1*r.freq,r.errnormpsd1*r.freq);
+xlog;plot_with_err(r.freq,r.signormpsd1*r.freq,r.errnormpsd1*r.freq);
 %% There is a warning that there are gaps in the LC (due to the telemetry dropouts) if verbose is set in foucalc. We ignore that for this simple analysis.
@@ Line 157: / Line 158: @@
 <pre>
 %% Choose Delta f/f, the logarithmic frequency resolution
-dff=0.15
+dff=0.15;
 reb_psd1 = rebin_fouquan(r.freq, r.signormpsd1, r.numavgall ; logfreq=dff);
 reb_err1 = rebin_fouquan(r.freq, r.errnormpsd1, r.numavgall ; logfreq=dff);
@@ Line 164: / Line 165: @@
 %% Need to weight the PSD error with the square root of the number of values the averaging was performed over
-hplot_with_err( reb_psd1.freq_lo, reb_psd1.freq_hi, reb_psd1.value*reb_psd1.freq_lo, reb_err1.value/sqrt(reb_err1.n)*reb_psd1.freq_lo );
+ylog;hplot_with_err( reb_psd1.freq_lo, reb_psd1.freq_hi, reb_psd1.value*reb_psd1.freq_lo, reb_err1.value/sqrt(reb_err1.n)*reb_psd1.freq_lo );
 </pre>
@@ Line 180: / Line 181: @@
 <pre>
 %% Rebin time lag results
-variable reb_lag12 = rebin_fouquan(r.freq, r.lag12, r.numavgall ; logfreq=dff);
+reb_lag12 = rebin_fouquan(r.freq, r.lag12, r.numavgall ; logfreq=dff);
-variable reb_errlag12 = rebin_fouquan(r.freq, r.errlag12, r.numavgall ; logfreq=dff);
+reb_errlag12 = rebin_fouquan(r.freq, r.errlag12, r.numavgall ; logfreq=dff);
 hplot_with_err(reb_lag12.freq_lo, reb_lag12.freq_hi, reb_lag12.value, reb_errlag12.value/sqrt(reb_errlag12.n) );
 %% Rebin coherence function
-variable reb_cof12 = rebin_fouquan(r.freq, r.cof12, r.numavgall ; logfreq=dff);
+reb_cof12 = rebin_fouquan(r.freq, r.cof12, r.numavgall ; logfreq=dff);
-variable reb_errcof12 = rebin_fouquan(r.freq, r.errcof12, r.numavgall ; logfreq=dff);
+reb_errcof12 = rebin_fouquan(r.freq, r.errcof12, r.numavgall ; logfreq=dff);
 hplot_with_err( reb_cof12.freq_lo, reb_cof12.freq_hi, reb_cof12.value, reb_errcof12.value/sqrt(reb_errcof12.n) );
 </pre>

@@ Line 113: / Line 113: @@
 <pre>
 %% Load the lightcurve and GTI file
 a=fits_read_table("full_sd_0.1s_0500-0900.lc");
 b=fits_read_table("full_sd_0.1s_2000-3000.lc");
-gtia=fits_read_table("full_sd_0.1s_0500-0900.lc[GTI0006]");
+gtia=fits_read_table("full_sd_0.1s_0500-0900.lc[GTI00006]");
-gtib=fits_read_table("full_sd_0.1s_2000-3000.lc[GTI0006]");
+gtib=fits_read_table("full_sd_0.1s_2000-3000.lc[GTI00006]");
 %% Get the start time such that LC starts at 0s