Change request #2659
Do not require a background model template for On/Off fitting with wstat
| Status: | Closed | Start date: | 08/02/2018 | |
|---|---|---|---|---|
| Priority: | Normal | Due date: | ||
| Assigned To: |  Knödlseder Jürgen | % Done: | 100% | |
| Category: | - | |||
| Target version: | 1.6.0 | |||
| Duration: |
Description
On/Off fitting with wstat does not need to assess the absolute background rate using the GCTAOnOffObservation::N_bgd() method, but only requires access to the background scaling factor alpha, computed by GCTAOnOffObservation::compute_alpha() and stored in the On spectrum. Although the later currently requires a background template to cope with spatial variations in the background rate, one could in the absence of a background template simply assume the same background rate in the On and Off regions, which is the standard analysis method when taking for example the reflected region method.
The usage of the IRF background template should therefore be optional, and the code should also be able to cope with an IRF that has no background template.
For the first, a flag should be added to the GCTAOnOffObservation constructor and GCTAOnOffObservation::set() method to
GCTAOnOffObservation::GCTAOnOffObservation(const GCTAObservation& obs,
const GModelSpatial& spatial,
const GEbounds& etrue,
const GEbounds& ereco,
const GSkyRegions& on,
const GSkyRegions& off,
const bool& use_irf_bkg = true)
void GCTAOnOffObservation::set(const GCTAObservation& obs,
const GModelSpatial& spatial,
const GSkyRegions& on,
const GSkyRegions& off,
const bool& use_irf_bkg = true)
csphagen can take control over the IRF background template usage.
The flag needs then also to be added to
void GCTAOnOffObservation::compute_bgd(const GCTAObservation& obs,
const GSkyRegions& off,
const bool& use_irf_bkg = true)
BACKRESP column of the Off spectrum with zeros in case that the flag is set to false. In that way the client knows that no background rate information is available, and GCTAOnOffObservation::N_bgd() will simply return zeros.
The flag needs also to be added to
void GCTAOnOffObservation::compute_alpha(const GCTAObservation& obs,
const GSkyRegions& on,
const GSkyRegions& off)
const bool& use_irf_bkg = true)
alpha factor simply reflects the difference in the solid angles between On and Off region. Note that in this case, alpha is energy independent, hence no loop over reconstructed energy is actually needed.
The dealing with IRFs without template should then be straight forward, the code should enter the same branch in absence of the background template as with use_irf_bkg = false.
Recurrence
No recurrence.
History
#1
Updated by Knödlseder Jürgen over 5 years ago
- Subject changed from Do not require a background model template for On/Off fitting with @wstat@ to Do not require a background model template for On/Off fitting with wstat
- Status changed from New to In Progress
- % Done changed from 0 to 50
I implemented the use_irf_bkg parameter in GCTAOnOffObservations, yet the code does not autodetect the absence of the IRF background template. The code works in a way that no background template is needed if use_irf_bkg=false. Hence if no background template exists, the code will throw an exception, except if use_irf_bkg=false.
#2
Updated by Knödlseder Jürgen over 5 years ago
I added a hidden use_irf_bkg parameter to csphagen to make use of the new functionality. Running now csphagen with use_irf_bkg=no will ignore the IRF background template and assume that the background rates in On and Off regions are identical.
Here a usage example:
$ csphagen prefix=onoff_const use_irf_bkg=no Input event list or observation definition XML file [events.fits] Calibration database [prod2] Instrument response function [South_0.5h] Input model definition XML file (if NONE, use point source) [$CTOOLS/share/models/crab.xml] Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] Start value for first energy bin in TeV [0.1] Stop value for last energy bin in TeV [100.0] Number of energy bins [20] Stack multiple observations into single PHA, ARF and RMF files? [no] Output observation definition XML file [onoff_obs_ref.xml] onoff_obs_const.xml Method for background estimation (REFLECTED|CUSTOM) [REFLECTED] Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] Right Ascension of source region centre (deg) (0-360) [83.63] Declination of source region centre (deg) (-90-90) [22.01] Radius of source region circle (deg) (0-180) [0.2]
Fitting the resulting data using the cstat statistic will result in a singular curvature matrix:
2018-08-03T08:46:58: +=========================================+ 2018-08-03T08:46:58: | Maximum likelihood optimisation results | 2018-08-03T08:46:58: +=========================================+ 2018-08-03T08:46:58: === GOptimizerLM === 2018-08-03T08:46:58: Optimized function value ..: 0.000 2018-08-03T08:46:58: Absolute precision ........: 0.005 2018-08-03T08:46:58: Acceptable value decrease .: 2 2018-08-03T08:46:58: Optimization status .......: singular curvature matrix encountered 2018-08-03T08:46:58: Number of parameters ......: 10 2018-08-03T08:46:58: Number of free parameters .: 4 2018-08-03T08:46:58: Number of iterations ......: 101 2018-08-03T08:46:58: Lambda ....................: 0.001 2018-08-03T08:46:58: Maximum log likelihood ....: -0.000 2018-08-03T08:46:58: Observed events (Nobs) ...: 3086.000 2018-08-03T08:46:58: Predicted events (Npred) ..: 0.000 (Nobs - Npred = 3086)
Fitting with
wstat however works:2018-08-03T08:47:34: +=========================================+ 2018-08-03T08:47:34: | Maximum likelihood optimisation results | 2018-08-03T08:47:34: +=========================================+ 2018-08-03T08:47:34: === GOptimizerLM === 2018-08-03T08:47:34: Optimized function value ..: 9.713 2018-08-03T08:47:34: Absolute precision ........: 0.005 2018-08-03T08:47:34: Acceptable value decrease .: 2 2018-08-03T08:47:34: Optimization status .......: converged 2018-08-03T08:47:34: Number of parameters ......: 10 2018-08-03T08:47:34: Number of free parameters .: 4 2018-08-03T08:47:34: Number of iterations ......: 2 2018-08-03T08:47:34: Lambda ....................: 1e-05 2018-08-03T08:47:34: Maximum log likelihood ....: -9.713 2018-08-03T08:47:34: Observed events (Nobs) ...: 3086.000 2018-08-03T08:47:34: Predicted events (Npred) ..: 3086.872 (Nobs - Npred = -0.871702810702573) 2018-08-03T08:47:34: === GModels === 2018-08-03T08:47:34: Number of models ..........: 2 2018-08-03T08:47:34: Number of parameters ......: 10 2018-08-03T08:47:34: === GModelSky === 2018-08-03T08:47:34: Name ......................: Crab 2018-08-03T08:47:34: Instruments ...............: all 2018-08-03T08:47:34: Instrument scale factors ..: unity 2018-08-03T08:47:34: Observation identifiers ...: all 2018-08-03T08:47:34: Model type ................: PointSource 2018-08-03T08:47:34: Model components ..........: "PointSource" * "PowerLaw" * "Constant" 2018-08-03T08:47:34: Number of parameters ......: 6 2018-08-03T08:47:34: Number of spatial par's ...: 2 2018-08-03T08:47:34: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) 2018-08-03T08:47:34: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) 2018-08-03T08:47:34: Number of spectral par's ..: 3 2018-08-03T08:47:34: Prefactor ................: 5.78794235100673e-16 +/- 1.16894725927285e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) 2018-08-03T08:47:34: Index ....................: -2.44692021286776 +/- 0.0166311257112595 [-0,-5] (free,scale=-1,gradient) 2018-08-03T08:47:34: PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient) 2018-08-03T08:47:34: Number of temporal par's ..: 1 2018-08-03T08:47:34: Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
For comparison, here the Crab model results for
wstat without ignoring the IRF background template. The results are very similar:2018-08-03T08:45:43: === GModelSky === 2018-08-03T08:45:43: Name ......................: Crab 2018-08-03T08:45:43: Instruments ...............: all 2018-08-03T08:45:43: Instrument scale factors ..: unity 2018-08-03T08:45:43: Observation identifiers ...: all 2018-08-03T08:45:43: Model type ................: PointSource 2018-08-03T08:45:43: Model components ..........: "PointSource" * "PowerLaw" * "Constant" 2018-08-03T08:45:43: Number of parameters ......: 6 2018-08-03T08:45:43: Number of spatial par's ...: 2 2018-08-03T08:45:43: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) 2018-08-03T08:45:43: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) 2018-08-03T08:45:43: Number of spectral par's ..: 3 2018-08-03T08:45:43: Prefactor ................: 5.78797906772499e-16 +/- 1.16894538223252e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) 2018-08-03T08:45:43: Index ....................: -2.44692374255593 +/- 0.0166310559574436 [-0,-5] (free,scale=-1,gradient) 2018-08-03T08:45:43: PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient) 2018-08-03T08:45:43: Number of temporal par's ..: 1 2018-08-03T08:45:43: Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
#4
Updated by Knödlseder Jürgen over 5 years ago
- Status changed from In Progress to Closed