lag_energy_spectrum: Replace lc_ref with elo/hi_ref

src/data/timing/lag_energy_spectrum.sl

@@ -27,8 +27,9 @@ private define verifyLightcurveTime(lc1, lc2){
 
 private define getLightcurve(lc, dt, dimseg){ 
   %{{{
-  variable data = fits_read_table(lc);
-  
+
+  variable data = typeof(lc) == String_Type ? fits_read_table(lc) : lc;
+
   if (struct_field_exists(data, "time") * struct_field_exists(data, "rate") != 1) {
     throw UsageError, sprintf("ERROR (%s): %s does not have a time and rate column! Exiting...", _function_name(), lc);
   }
@@ -45,7 +46,7 @@ private define calculateSummedReferenceLightcurve(lc_list, index, dt, dimseg){
 
   %% todo: clean up this horrible if-else mess!
   
-  variable channel_of_interest = fits_read_table(lc_list[index]);
+  variable channel_of_interest = typeof(lc_list[index]) == String_Type ? fits_read_table(lc_list[index]) : lc_list[index];
 
   variable jj, lc;
   variable summed_rate = Double_Type[length(channel_of_interest.rate)];
@@ -53,9 +54,6 @@ private define calculateSummedReferenceLightcurve(lc_list, index, dt, dimseg){
   if (qualifier_exists("elo") * qualifier_exists("ehi") * qualifier_exists("elo_ref") * qualifier_exists("ehi_ref") == 1) {
     %% If we have energy information about the ref LC we want to
     %% only subtract the COI in case it's within the ref LC borders.
-    %% 
-    %% todo: this reproduces the code in the main function, need to
-    %% be exported into a getReferenceLightcurve function
     variable elo_ref = qualifier("elo_ref");
     variable ehi_ref = qualifier("ehi_ref");
     variable elo = qualifier("elo");
@@ -70,23 +68,23 @@ private define calculateSummedReferenceLightcurve(lc_list, index, dt, dimseg){
 						  _function_name(), index, elo[index], ehi[index], elo_ref, ehi_ref);
 	continue;
       } else {
-	lc = fits_read_table(lc_list[jj]);
+	lc = typeof(lc_list[jj]) == String_Type ? fits_read_table(lc_list[jj]) : lc_list[jj];
 	verifyLightcurveTime(channel_of_interest, lc);
 	if (qualifier("verbose", 0) > 1) vmessage("  ***(%s): Summing LC (index %d, %g-%g keV) to Ref. LC (%g-%g keV).",
 						  _function_name(), jj, elo[jj], ehi[jj], elo_ref, ehi_ref);
 	summed_rate += lc.rate;
       }
     }	     
-  } else {    
+  } else {
     %% If no energy information about the reference lightcurve is
     %% provided, we assume that it spans the full energy range.
     _for jj (0, length(lc_list)-1, 1){
       if (jj == index){
 	if (qualifier("verbose", 0) > 1) vmessage("  ***(%s): Not adding COI LC (index %d) to Ref. LC, assuming that the Ref. LC goes over all energies.",
-						  index, _function_name());
+						  _function_name(), index);
 	continue;
       } else {
-	lc = fits_read_table(lc_list[jj]);
+	lc = typeof(lc_list[jj]) == String_Type ? fits_read_table(lc_list[jj]) : lc_list[jj];
 	verifyLightcurveTime(channel_of_interest, lc);
 	summed_rate += lc.rate;
       }
@@ -169,10 +167,6 @@ private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){
 		 "Analyzed lightcurves:", lc_list]
     };
     
-    if (qualifier_exists("lc_ref")) {
-      hist.comment = [hist.comment, sprintf("Reference LC: %g-%g keV", qualifier("elo_ref", "?"), qualifier("ehi_ref", "?"))];
-    }
-    
     fits_write_binary_table(fp, extname, sdata, keys, hist);
   }
   
@@ -232,10 +226,6 @@ define lag_energy_spectrum(lc_list, dimseg) {
 %    
 %    returns: lag-energy-spectrum as a struct with fields lag, errlag
 %\qualifiers{
-%    \qualifier{lc_ref}{Filename of a (broad) reference lightcurve. If not
-%                set, the reference lightcurve will consist of the
-%                summed up energy lightcurves except for the current 
-%                channel-of-interest.}
 %    \qualifier{dt}{Time resolution of the lightcurve in seconds
 %                [default: TIMEDEL keyword]}
 %    \qualifier{deadtime}{Detector deadtime in seconds [default: 0.0s]}
@@ -246,8 +236,8 @@ define lag_energy_spectrum(lc_list, dimseg) {
 %    \qualifier{verbose}{Increase verbosity [default=0]}
 %    \qualifier{outfile}{String_Type. If set, write various timing products
 %                into a FITS file.}
-%    \qualifier{elo, ehi}{Double_Type: Energy grid, only used for writing outfile}
-%    
+%    \qualifier{elo, ehi}{Double_Type: Energy grid}
+%    \qualifier{elo_ref, ehi_ref}{Boundaries of reference lightcurve}
 %}
 %\notes
 %    * The default error computation is Nowak et al., 1999, ApJ, 510,
@@ -316,9 +306,13 @@ define lag_energy_spectrum(lc_list, dimseg) {
 %    }
 %\seealso{segment_lc_for_psd, foucalc, colacal}
 %!%-
-  variable dt = qualifier("dt", fits_read_key(lc_list[0], "TIMEDEL"));
-  if (dt == NULL) return sprintf("ERROR (%s): Couldn't receive the time resolution from TIMEDEL, need dt qualifier!", _function_name());
-
+  try {
+    variable dt = qualifier_exists("dt") ? qualifier("dt") : fits_read_key(lc_list[0], "TIMEDEL");
+  } catch TypeMismatchError:
+  {
+    return vmessage("ERROR (%s): Couldn't receive the time resolution from TIMEDEL, need dt qualifier!", _function_name());
+  }
+  
   variable f_lo = qualifier("f_lo", [1/(dt*dimseg)]);
   variable f_hi = qualifier("f_hi", [1/(2*dt)]);  % Nyquist frequency
   variable verbose = qualifier("verbose", 0);
@@ -356,12 +350,6 @@ define lag_energy_spectrum(lc_list, dimseg) {
     freqres_errlags   = Array_Type[N_energies]
   };
 
-  %% Load reference band here, such that it is not loaded again for each energy channel
-  variable lc_ref;
-  if (qualifier_exists("lc_ref")){
-    lc_ref = getLightcurve(qualifier("lc_ref"), dt, dimseg;; __qualifiers);
-  }
-  
   %% Loop through energy-resolved lightcurves (channel-of-interests)
   variable ii, jj;
   _for ii (0, N_energies-1, 1){
@@ -369,38 +357,14 @@ define lag_energy_spectrum(lc_list, dimseg) {
     if (verbose > 0) vmessage("  ***(%s): Calculate time lags for %s", _function_name(), lc_list[ii]);
     
     variable channel_of_interest = getLightcurve(lc_list[ii], dt, dimseg;; __qualifiers);
-    
-    variable foucalc_quals = struct{normtype=normtype, deadtime=deadtime};
-    
-    if (qualifier_exists("lc_ref")) {
-      %% If a reference lightcurve is given we need to determine
-      %% whether the current channel-of-interest is within the
-      %% reference band. In this case we have to correct the Poisson
-      %% noise in the cross spectrum. For this we need the energy
-      %% information of the channels. 
-      if (qualifier_exists("elo") * qualifier_exists("ehi") * qualifier_exists("elo_ref") * qualifier_exists("ehi_ref") != 1) {
-	throw ReadError, "ReadError: If a reference lightcurve is given one must also pass the elo, ehi,\n" +
-	                 "           elo_ref, ehi_ref qualifiers for the Poisson noise correction of the CPD.";
-      }
-      
-      variable elo_ref = qualifier("elo_ref");
-      variable ehi_ref = qualifier("ehi_ref");
-      variable elo = qualifier("elo")[ii];
-      variable ehi = qualifier("ehi")[ii];
-      
-      if ((elo >= elo_ref) and (ehi <= ehi_ref)) {
-	if (verbose > 0) vmessage("  ***(%s): Subject LC (%g-%g keV) is within reference LC (%g-%g keV). Correcting for Poisson noise.",
-				  _function_name(), elo, ehi, elo_ref, ehi_ref);
-	foucalc_quals = struct_combine(foucalc_quals, "correct_correlated_Poisson_noise_in_CPD");
-      }
-    } else {
-      lc_ref = calculateSummedReferenceLightcurve(lc_list, ii, dt, dimseg;; __qualifiers);     
-    }
+
+    variable lc_ref = calculateSummedReferenceLightcurve(lc_list, ii, dt, dimseg;; __qualifiers);     
 
     variable fouquant = foucalc(struct{time  = channel_of_interest.time,
                                        rate1 = lc_ref.rate,
                                        rate2 = channel_of_interest.rate},
-				dimseg;; foucalc_quals);
+				dimseg; normtype=normtype, deadtime=deadtime);
+
     variable numseg = fouquant.numavgall[0];
     
     dat.freq[ii] = fouquant.freq;