ctools

h2. 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:

<pre>
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
</pre>

h2. 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.

h2. 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).

h2. Place A: ctselect::get_parameters():

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

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

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

<pre>
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);

</pre>

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.

h3. Fix:

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

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

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.)

h2. 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:

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

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:

<pre>
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);

...

</pre>

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).

h3. Fix:

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

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

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

h2. 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.

h2. 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?