Bug #1329

ctselect improperly applies time cuts when event lists do not use the CTA reference (MET) time

Added by Kelley-Hoskins Nathan over 9 years ago. Updated over 9 years ago.

Status:ClosedStart date:10/02/2014
Priority:NormalDue date:
Assigned To:Kelley-Hoskins Nathan% Done:

100%

Category:-
Target version:00-08-00
Duration:

Description

What was I trying to do?

I’ve been loading veritas event lists (and a dummy hess event list) and feeding them to ctselect in a python script, with the goal of breaking up observations into n-second-long chunks. Example code:

import gammalib
import ctools

ob = gammalib.GCTAObservation('VERITAS')
ob.load( '/path/to/my/eventlist.fits')
ob.id( '12345' )
obs = gammalib.GObservations()
obs.append(ob)

sel = ctools.ctselect( obs )
sel['tmin'].real( 100.0)
sel['tmax'].real( 230.0)
sel['ra'  ].real( -1.0 )
sel['dec' ].real( -1.0 )
sel['rad' ].real( -1.0 )
sel['emin'].real(  0.1 )
sel['emax'].real( 100.0)
sel.run()
obs_cut = sel.obs()

# obs has 1400 events
# obs_cut has 0 events

However: ctselect has been cutting all events,

irregardless of the time cuts you give it, in the veritas and hess event lists. But, if you give ctselect a simulated cta events list, the cuts work fine. After poking and prodding it for a few weeks, I’ve finally pinned down why it does this.

Cause?

Its due to two separate issues in ctselect.cpp, where ctselect assumes any input GObservations uses cta’s reference time (MET) when saving/using time cut boundaries. My veritas reference time (MET) is something else, and HESS as well is has a third mjd reference time (MET).

Place A: ctselect::get_parameters():

When you first save your tmin/tmax parameters to a ctselect:

import ctools
sel=ctools.ctselect()
sel['tmin'].real(100.0)
sel['tmax'].real(200.0)

they call these lines in ctselect::get_parameters() :

from $CTOOLS/src/ctselect/ctselect.cpp : 
ctselect::get_parameters(void) :

    m_tmin = (*this)["tmin"].real();
    m_tmax = (*this)["tmax"].real();

    ...

    // Set time interval with input times given in CTA reference
    // time (in seconds)
    m_timemin.set(m_tmin, m_cta_ref);
    m_timemax.set(m_tmax, m_cta_ref);

This takes your tmin/tmax parameters and saves them to m_tmin/m_tmax, before converting them to cta time and saving them to m_timemin/m_timemax .

To properly apply this conversion, the user needs to convert their cut time parameters tmin/tmax into the the reference time cta is using. This means the above python code will only allow events between 100.0 and 200.0 seconds since the cta reference time . Since VERITAS and HESS have picked their own reference times, and cta’s reference time may change in the future, VERITAS and HESS don’t have an automatic way to convert their time cuts times to whatever time reference ctselect is using.

Fix:

One way to combat this would be to add a public function to ctselect that can return m_cta_ref, aka:

in $CTOOLS/src/ctselect/ctselect.cpp:
GTimeReference ctselect::reference() {
    return m_cta_ref ;
}

as well as adding more documentation to the parameters tmin and tmax, indicating these are times in seconds-since-cta-reference-time (I couldn’t find much documentation on them when you run ctselect in interactive mode, it tells you CTA MET time, but for a VERITAS or HESS event-list-loading-script, it needs to know what time reference cta is using, so the problem still stands). At first I (wrongly) assumed they were times in seconds since the beginning of the observation :( ).

Alternatives:
  • changing tmin and tmax to accept GTimes(), or
  • change tmin and tmax params so that they are full MJD times ( sel['tmin’].real(54232.3384919), etc), instead of seconds-since-reference-mjd times. ( < I kind of like this solution the most, but I suspect this change will mess up some existing examples, and the interactive mode for apps will get more awkward.)

Place B: ctselect::select_events(...):

When time cuts are applied, they are assembled (m_timemin/m_timemax from above) into a string, which is added onto the fits filename, which is then opened like a normal fits filename:

double cmin = 224.0 ;
double cmax = 337.0 ;
GFits file("/path/to/myfile.fits[EVENTS][TIME >= "+tostring(cmin)+" && TIME <= "+tostring(cmax)]") ;

ctselect assembles the “filename[expression], and cfitsio handles the expression when opening the file. For example, the above string would look in table “EVENTS” and cut on column “TIME”. That TIME column is in whatever time reference is specified in its EVENTS header MJDREF However, when ctselect assembles 'cmin’ and 'cmax’ in the lines below:

from $CTOOLS/src/ctselect/ctselect.cpp : 
ctselect::select_events(GCTAObservation* obs, const std::string& filename) :

        ...

        // this reference for filtering.
        double tmin = gti.tstart().convert(m_cta_ref);
        double tmax = gti.tstop().convert(m_cta_ref);

        ...

        // Format time with sufficient accuracy and add to selection string
        sprintf(cmin, "%.8f", tmin);
        sprintf(cmax, "%.8f", tmax);

        ...

ctselect converts these times to seconds-since-cta-reference-time, (rather than seconds-since-fits-file-reference-time) before feeding them to the fits file expression. Since the cut times cmin/cmax are in one reference, and the event times in the fits file are in another reference, all events get cut, unless you are very very lucky with your event times (I never was).

Fix:

To fix this, I think the 'convert’ lines could be replaced with something like:

        // this reference for filtering.
        double tmin = gti.tstart().convert( obs->events().reference() );
        double tmax = gti.tstop().convert( obs->events().reference() );

but I’m not sure where in the observations the event-list-wide time reference is stored.

Impact?

Both of these problems need to be addressed before the time cuts work properly in ctselect. I believe this problem will also crop up when ctselect is applied to any GObservations object that has a non-cta time reference, not just those loaded from fits files, which makes me think this issue is a good one to fix in the name of multi-observatory studies.

Why did you make this report so long?

I’ve suggested some changes above, but since ctselect performs such an important task in ctools, I thought I’d ask if anyone has suggestions or alternative solutions for dealing with these?


Recurrence

No recurrence.

History

#1 Updated by Kelley-Hoskins Nathan over 9 years ago

  • Description updated (diff)
  • Assigned To set to Kelley-Hoskins Nathan

#2 Updated by Knödlseder Jürgen over 9 years ago

Thanks for bringing this important issue up. This is indeed an annoying issue, and I’m sorry that it prevented you for going ahead. In fact, my primary concern was to implement a time reference correctly, but I left somehow the precise definition of the tmin and tmax parameters in ctselect aside (and I probably forgot about the fact that every of the existing IACT will have a different time reference system).

I think I would opt for “Place B”, as one could image that ctselect gets on input an observation definition XML file that mixes data sets with different time references (e.g. VERITAS & MAGIC data). ctselect should then take the times in a standard (and human readable) format, and it may be good to change tmin and tmax to a string to allow different formats, such as UTC strings alike 2014-09-10T10:24:41 or MJD in seconds. I would suggest to add methods GTime::utc() to digest the conversion of UTC strings by GTime. Some code would be needed to analyze the time string and to call the appropriate GTime method, we can put such code in the GTools file.

The observation wide time reference is stored in the GTI of the event list or event cube, see for example GCTAEventList::read(). Hence the code

double tmin = gti.tstart().convert(gti.reference());
double tmax = gti.tstop().convert(gti.reference());

should work, I guess.

#3 Updated by Kelley-Hoskins Nathan over 9 years ago

  • % Done changed from 0 to 100

For Place A, I’ve added a function to expose the time reference ctselect is using, as branch
origin/1329-ctselect-expose-timref .

For Place B, that change you suggested works fine, and handles different fits files properly. I pushed it as branch
origin/1329-ctselect-fits-time-conversion

if it hasn’t already been done.

#4 Updated by Kelley-Hoskins Nathan over 9 years ago

  • Status changed from New to Pull request

I think 'Pull request” is how I request the branches origin/1329-ctselect* get merged into devel?

#5 Updated by Knödlseder Jürgen over 9 years ago

Pull request is indeed the right thing to do.

While merging I recognized that there is still

    // Set time interval with input times given in CTA reference
    // time (in seconds)
    m_timemin.set(m_tmin, m_cta_ref);
    m_timemax.set(m_tmax, m_cta_ref);

in ctselect::get_parameters, hence the tmin and tmax time values are still in CTA reference time, right?

#6 Updated by Kelley-Hoskins Nathan over 9 years ago

Yes, tmin and tmax are still in CTA reference time. I pushed another branch, origin/1329-ctselect-expose-timeref , which adds a function ctselect::reference(), that should let a user get around the problem mentioned at Place A.

This exposed time reference would still allow an external app that uses ctselect to be able to properly convert their times for tmin and tmax, even if the CTA reference time is changed in the future to something else.

import gammalib
import ctools
tcutmin = gammalib.GTime()
tcutmin.mjd(54763.212452) # starting mjd
tcutmax = gammalib.GTime()
tcutmax.mjd(56372.382392) # ending mjd
sel = ctools.ctselect()
sel['tmin'].real( tcutmin.convert( sel.reference() ) )
sel['tmax'].real( tcutmax.convert( sel.reference() ) )

So, I think its ok that 'm_timemin/m_timemax’ are still set to reference 'm_cta_ref’ in ctselect::set_parameters() .

#7 Updated by Knödlseder Jürgen over 9 years ago

  • Status changed from Pull request to Closed
  • Target version set to 00-08-00

Okay, a added the method to ctselect. I called it by the way time_reference() as reference is too general (who knows what references we may get in the future). I still think that in the long run all times should be input in ctools in a specific reference format, such as MJD and/or UTC. This would make the need for the method probably obsolete in the future.

Also available in: Atom PDF