Implement ON-OFF fitting
After the work done during the first coding sprint, classes now exist to hold most data required for an ON-OFF classical analysis (regions, ON and OFF spectra, ARF and RMF arrays,...). What remains to be implemented now is the optimization part. This should allow both the fitting of a given spectral shape and the determination of flux points.
As a starting point, it was decided to add to the gammalib only the mimimum that is required, and to leave higher-level functionalities to specific ctools. These functionalities can then be moved to gammalib if it seems more appropriate.
Anticipating on a possible evolution of the GObservations/GObservation, the optimization function will be attached to the observation and not to the container of observations (contrary to the initial design of the gammalib).
#1 Updated by Martin Pierrick about 10 years ago
- % Done changed from 0 to 50
Created a dedicated branch and pushed to it a version of the code that compiles (on my machine). Still remains to be tested through python scripts. Also, the branch forks a bit behind the current state of the code. Will do the merging later.
The likelihood function is at the observation level, and the call of the optimizer is done at the container level (now the standard also for the rest of the gammalib). The likelihood function slightly differs in its arguments from the other likelihood functions. Can homogeneize that later if deemed necessary.
No unit test written so far. I will first try to play with this new functionality.
#2 Updated by Martin Pierrick about 10 years ago
- File test_OnOff.py added
Here is a script to test the ON-OFF fitting (and associated functionalities by the way). The script simulates two wobble observations with one ON region and any number of OFF regions distributed around the centre of each pointing (rough method). It allows to compare ON-OFF analysis and unbinned analysis over the entire field of view. So far, it runs and tells that no fit is performed for the ON-OFF analysis. Problem of singular matrix. Inspection under way...
#4 Updated by Martin Pierrick about 10 years ago
- File test_OnOff.py added
Attached a new version of the test script for the ON-OFF analysis. Corrected a few bugs and added many prints for checking. The gammalib code was also modified but currently yields a segmentation fault. Hopefully this will be solved during the second coding sprint...
#6 Updated by Martin Pierrick about 10 years ago
- File test_OnOff.py added
Latest version of the script.
A few notes:
- OFF region creation is using a rough method that works only close to equator
- Python interface at the moment does not allow using derived classes specific methods for GSkyRegion objects
- Hard-coded path to data base
- ON-OFF fit not performed properly (gammalib problem, not the script)
#7 Updated by Martin Pierrick about 10 years ago
Summary of the second coding sprint at DESY:
- ON-OFF fitting debugged, now works and give reasonable results (on just one test case)
- Code layout was made more similar to that of the global fit case (the standard gammalib/ctool one)
- GSkyRegion derived classes were debugged (constructors, DS9 file saving)
- GSkyRegionRing class created and GSkyRegionMap almost done
- Some test cases for the GSkyRegionRing class
- Draft python scripts for ON-OFF analysis ctools: ctonoffreg, ctonoffirf, ctonofffit
- Modified GCTAPointing:
- Implemented a GHorizDir class
- Reads in a pointing fits file, describing the pointing position in alt/az as a function of time
- Added dir_horiz(time) that interpolates the pointing table and returns alt az
Remains to be done for the ON-OFF analysis:
- Add handling of exclusion map (a GSkyRegion* member in GCTAOnOffObservation)
- Add energy dependence (ON region size varying with energy; could be computed from a table to be loaded)
- Make effective area computation more general (to allow ring OFF regions...; the integration can be done in a ctool so we would just need storage)
- Update alpha parameter calculation (very rough at the moment)
- Precomputation of OFF region solid angle (done too many times in likelihood computation at the moment)
- Write test cases for the ON-OFF analysis
- Have csscripts or ctools modules for a complete analysis
#12 Updated by Martin Pierrick about 10 years ago
Not many tests performed so far, just the comparison of fitting a power-law spectrum by unbinned and ON-OFF analyses, for a Crab-like source (so a strong one). What remains to be done is testing the code for a weaker source, and determining spectral points using a NodeFunction.
#14 Updated by Martin Pierrick almost 8 years ago
Summary of progress after 5th coding sprint at DESY:
- Caught up with the development of the rest code over the past two years
- Replaced the old background model by GCTAModelIrfBackground (currently no other choice)
- Implemented clean calculation of alpha (ratio of ON to OFF effective areas)
- Implemented energy-dependent alpha
- Implemented dependencies on theta, phi, zenith, azimuth (instrument coordinates)
- Added a precalculation of the background rate in OFF regions (method compute_bgd() in GCTAOnOffObservation)
- Revised/corrected calculation of expected ON and OFF counts
- No modification to the likelihood calculation part
- The code runs and returns meaningful results
- ... but pull distributions show that true values are not exactly recovered on average
In the only tested case so far, a simple point source with power-law spectrum and two wobble observations, the average prefactor over 100 runs is about 5% below the true value. The photon index is also slightly below true value but the deviation is not as large. As a comparison, an unbinned analysis over the same setup leads to perfect recovery of the true values (on average).
- Investigate the systematic offset in prefactor mentioned above
- Use new class Gfilename where needed
- Move includes from hpp to cpp (new practice)
- There is currently no check on region overlap; it is left to the user to input correct parameters
- The effective area and background rates are currently taken at the center of circular regions and not integrated over their full extent; this is probably alright for small regions but it needs to be improved upon for the general case
- The code currently handles only collections of circular regions; at least a ring region should be added.
That’s all folks !
#15 Updated by Martin Pierrick almost 8 years ago
Attached two scripts:
1- test_OnOff.py for testing the ON-OFF analysis on a single run: observations are simulated for two wobble pointings on either side of a point source with power-law spectrum; ON-OFF analysis is setup and executed; unbinned analysis is performed for comparison
2- pull_OnOff.py to perform a series of such runs and plot pull distributions
Both scripts are self-contained.
#18 Updated by Knödlseder Jürgen almost 8 years ago
- Target version set to 1.1.0
- % Done changed from 70 to 80
I started with integrating the code.
First I went over the code and adapt it to the coding conventions. I also needed to adjust the so far existing unit test at it was broken due to the additional of the models parameter in the
I also removed the
GCTAOnOffObservations class as it is no longer needed. Now
GCTAOnOffObservation objects can be simply appended to the
GObservations container (see #1669). This means that
GCTAOnOffObservation now derives from
- there is no complete unit test that checks the On/Off fitting functionality. Such a test should be added to the
test_CTA.pyscript (there is already a unit test, but it does so far not perform On/Off fitting)
- the computation of the Hessian in
GCTAOnOffObservation::likelihood()looks strange. Could you please write down the intended formula in the documentation (can use LaTeX style), and then check whether the formula is indeed implemented?
GCTAOnOffObservation::compute_alphamethod actually computes the ratio between effective areas, not exposures. Am I correct that it should compute the ratio of exposures? What about the solid angle?
- in the
GCTAOnOffObservation::compute_rmfmethod, if a weighting is done, should it not be done over effective area times livetime? Where is the time coming in? Could you write down the formula in the documentation what the method should be doing? Also, should the
aeffmethod not use true instead of reconstructed energy?
- shoudn’t the
model_off()methods be private (it’s always simpler to make a private method public when really needed then vice versa)
- are the
Furthermore I see that the PHA data format requires an EXPOSURE keyword (the integration time (in seconds) for the PHA data, assumed to be corrected for deadtime, data drop-outs etc.). So this is a kind of livetime. I think we should add this attribute to the
GPha class, and then the storage of ontime and offtime is not longer needed. By the way I see
m_ontime = obs.livetime(); m_offtime = obs.livetime();which means that you store in fact exposure as ontime, which is not very correct.
#19 Updated by Knödlseder Jürgen almost 8 years ago
I pushed a branch with the current code that has been rebased on devel into the gammalib repo: 1044-implement-on-off-fitting
Before now merging in the code, could you add a unit test so that we are sure that we test the likelihood fitting when running make check?