ctools

The issue came from GCTAOnOffObservation(const GObservations& obs) which did not correctly handle the case where the IRF background template is absent. I modified the code so that the handling is now done correctly.

Here is the ctlike result after having done the change:

2018-10-04T08:33:36: +=================================+
2018-10-04T08:33:36: | Maximum likelihood optimisation |
2018-10-04T08:33:36: +=================================+
2018-10-04T08:33:36:  >Iteration   0: -logL=209.552, Lambda=1.0e-03
2018-10-04T08:33:36:  >Iteration   1: -logL=97.156, Lambda=1.0e-03, delta=112.396, step=1.0e+00, max(|grad|)=266.954813 [Index:3]
2018-10-04T08:33:36:  >Iteration   2: -logL=14.417, Lambda=1.0e-04, delta=82.740, step=1.0e+00, max(|grad|)=69.552529 [Index:3]
2018-10-04T08:33:36:  >Iteration   3: -logL=11.287, Lambda=1.0e-05, delta=3.130, step=1.0e+00, max(|grad|)=41.225603 [Index:3]
2018-10-04T08:33:36:  >Iteration   4: -logL=10.609, Lambda=1.0e-06, delta=0.678, step=1.0e+00, max(|grad|)=-0.660643 [Index:3]
2018-10-04T08:33:36:  >Iteration   5: -logL=10.604, Lambda=1.0e-07, delta=0.005, step=1.0e+00, max(|grad|)=0.563913 [Index:3]
2018-10-04T08:33:36:  
2018-10-04T08:33:36:  +==================+
2018-10-04T08:33:36:  | Curvature matrix |
2018-10-04T08:33:36:  +==================+
2018-10-04T08:33:36:  === GMatrixSparse ===
2018-10-04T08:33:36:   Number of rows ............: 6
2018-10-04T08:33:36:   Number of columns .........: 6
2018-10-04T08:33:36:   Number of nonzero elements : 4
2018-10-04T08:33:36:   Number of allocated cells .: 516
2018-10-04T08:33:36:   Memory block size .........: 512
2018-10-04T08:33:36:   Sparse matrix fill ........: 0.111111111111111
2018-10-04T08:33:36:   Pending element ...........: 0
2018-10-04T08:33:36:   Fill stack size ...........: 0 (none)
2018-10-04T08:33:36:   0, 0, 0, 0, 0, 0
2018-10-04T08:33:36:   0, 0, 0, 0, 0, 0
2018-10-04T08:33:36:   0, 0, 0.152622527932884, -13.5947132096962, 0, 0
2018-10-04T08:33:36:   0, 0, -13.5947132096962, 1313.28746416818, 0, 0
2018-10-04T08:33:36:   0, 0, 0, 0, 0, 0
2018-10-04T08:33:36:   0, 0, 0, 0, 0, 0
2018-10-04T08:33:36: 
2018-10-04T08:33:36: +=========================================+
2018-10-04T08:33:36: | Maximum likelihood optimisation results |
2018-10-04T08:33:36: +=========================================+
2018-10-04T08:33:36: === GOptimizerLM ===
2018-10-04T08:33:36:  Optimized function value ..: 10.604
2018-10-04T08:33:36:  Absolute precision ........: 0.005
2018-10-04T08:33:36:  Acceptable value decrease .: 2
2018-10-04T08:33:36:  Optimization status .......: converged
2018-10-04T08:33:36:  Number of parameters ......: 6
2018-10-04T08:33:36:  Number of free parameters .: 2
2018-10-04T08:33:36:  Number of iterations ......: 5
2018-10-04T08:33:36:  Lambda ....................: 1e-08
2018-10-04T08:33:36:  Maximum log likelihood ....: -10.604
2018-10-04T08:33:36:  Observed events  (Nobs) ...: 534.000
2018-10-04T08:33:36:  Predicted events (Npred) ..: 534.505 (Nobs - Npred = -0.504815279792638)
2018-10-04T08:33:36: === GModels ===
2018-10-04T08:33:36:  Number of models ..........: 1
2018-10-04T08:33:36:  Number of parameters ......: 6
2018-10-04T08:33:36: === GModelSky ===
2018-10-04T08:33:36:  Name ......................: Crab
2018-10-04T08:33:36:  Instruments ...............: all
2018-10-04T08:33:36:  Instrument scale factors ..: unity
2018-10-04T08:33:36:  Observation identifiers ...: all
2018-10-04T08:33:36:  Model type ................: PointSource
2018-10-04T08:33:36:  Model components ..........: "PointSource" * "PowerLaw" * "Constant" 
2018-10-04T08:33:36:  Number of parameters ......: 6
2018-10-04T08:33:36:  Number of spatial par's ...: 2
2018-10-04T08:33:36:   RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
2018-10-04T08:33:36:   DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
2018-10-04T08:33:36:  Number of spectral par's ..: 3
2018-10-04T08:33:36:   Prefactor ................: 4.6694286848533e-15 +/- 9.16906665122574e-16 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2018-10-04T08:33:36:   Index ....................: -2.66112738324749 +/- 0.0988449649199173 [-0,-5]  (free,scale=-1,gradient)
2018-10-04T08:33:36:   PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
2018-10-04T08:33:36:  Number of temporal par's ..: 1
2018-10-04T08:33:36:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2018-10-04T08:33:36: 
2018-10-04T08:33:36: +==============+
2018-10-04T08:33:36: | Save results |
2018-10-04T08:33:36: +==============+
2018-10-04T08:33:36:  Model definition file .....: result_crab_plaw_stacked.xml
2018-10-04T08:33:36:  Covariance matrix file ....: NONE
2018-10-04T08:33:36: 
2018-10-04T08:33:36: Application "ctlike" terminated after 1 wall clock seconds, consuming 0.012137 seconds of CPU time.

I also added a NaN check to GCTAOnOffObservation::wstat_value so that problems can in the future be tracked down more easily. I will no merge in the code.

There is still a difference in the normalisation. Here the results for the joint analysis

2018-10-08T11:45:04:   Prefactor ................: 5.04442310978829e-17 +/- 3.89446356042083e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2018-10-08T11:45:04:   Index ....................: -2.63791542087657 +/- 0.0845343261931618 [-0,-5]  (free,scale=-1,gradient)
2018-10-08T11:45:04:   PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)

2018-10-08T13:01:40:   Prefactor ................: 2.05784483532498e-16 +/- 1.63719727869171e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2018-10-08T13:01:40:   Index ....................: -2.63551070265515 +/- 0.0862634071633871 [-0,-5]  (free,scale=-1,gradient)
2018-10-08T13:01:40:   PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)

It turns out that the 4 H.E.S.S. runs combine observations with a different offset angle of the Crab and hence the background rate of the Off regions vary between the first two and the second two runs. We always need an assumption about the spatial background rate distribution, and in absence of a template the probably best thing to do is to assume that the background rate is constant. That’s what the actual code does. This assumption is however violated in the present case, and I see no obvious solution to this problem. Stacking seems simply not possible in the absence of a background template!

It turned out that the exposure time was not correctly written by the GPha write method (it was in fact overwritten by the exposure time of the first stacked run). I corrected this bug and now the results look more reasonable:

2018-10-08T13:50:18:   Prefactor ................: 5.07679333256999e-17 +/- 4.09478853758408e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2018-10-08T13:50:18:   Index ....................: -2.64578609753126 +/- 0.0880065744992247 [-0,-5]  (free,scale=-1,gradient)
2018-10-08T13:50:18:   PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)

Concerning csresspec, the issue is that you need to run the script using the statistic=wstat option:

$ csresspec statistic=wstat components=yes
Input event list, counts cube, or observation definition XML file [obs_crab_onoff_stacked.xml] 
Input model definition XML file [results.xml] 
Residuals computation algorithm (SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE) [SIGNIFICANCE] 
Output residual spectrum file [resspec.fits] 
$CTOOLS/share/examples/python/show_residuals.py resspec.fits

Merged into devel. Luigi, Andreas, can you tell me how this fix works on your side?

1) With the current devel version I obtain comparable results for joint and stacked analysis of the Crab. I still recommend to implement the weighting of BACKSCAL on N_OFF for stacking without background template. This will allow for more general cases to be accommodated. Josh told me that the same strategy is used in VERITAS as well.
2) I still have problems with calculating stacked residuals from multiple On/Off observations that were analyzed jointly. I’m copying here my code and the results

observations = 'obs_crab_selected.xml'
models = 'crab_models.xml'
onoff_obs = 'obs_crab_onoff_joint.xml'

phagen = cscripts.csphagen()

phagen['inobs'] = observations
phagen['inmodel'] = models
phagen['srcname'] = 'Crab'
phagen['outobs'] = onoff_obs
phagen['ebinalg'] = 'LOG'
phagen['emin'] = 0.66
phagen['emax'] = 30.
phagen['enumbins'] = 20
phagen['coordsys'] = 'CEL'
phagen['ra'] = 83.63
phagen['dec'] = 22.01
phagen['rad'] = 0.2
phagen['bkgmethod'] = 'REFLECTED'
phagen['use_irf_bkg'] = False
phagen['maxoffset'] = 2.0
phagen['stack'] = False

phagen.execute()

fit_models = 'crab_fit_models_joint.xml'

like = ctools.ctlike()

like['statistic'] = 'WSTAT'
like['inobs'] = onoff_obs
like['inmodel'] = models
like['outmodel'] = fit_models

like.execute()

residuals = 'residuals_joint_stacked.fits'
res = cscripts.csresspec()

res['statistic'] = 'WSTAT'
res['inobs'] = onoff_obs
res['inmodel'] = fit_models
res['components'] = True
res['algorithm'] = 'SIGNIFICANCE'
res['outfile'] = residuals
res['stack'] = True

res.execute()

The resulting residuals are

I’m using WSTAT in csresspec. The background model is identically 0. Please, let me know if there are any issues, but otherwise could it be something in csresspec that does not work properly without a background template?

Just to recap where we are for the 9th coding sprint:

The GCTAOnOffObservation::compute_bgd computes the BACKRESP column in the Off spectrum, which without a background model is just the solid angle of the Off region region for a given observations. Note that in the commit before, the values of the BACKRESP column were simply set to 0 in case that no background model exists.

The GCTAOnOffObservation::compute_alpha computes the alpha factor, which without a background model is just the energy independent ratio between the solid angle of the On and Off regions.

The wstat method is not using the BACKRESP column information, it is however used when stacking in the GCTAOnOffObservation::GCTAOnOffObservation(const GObservations& obs) constructor:

// Compute background scaling factor contribution
for (int i = 0; i < m_on_spec.size(); ++i) {
    double  background = onoff->off_spec()["BACKRESP"][i];
    double  alpha      = onoff->on_spec().backscal(i);
    double  scale      = alpha * background * exposure;
    m_on_spec.backscal(i, m_on_spec.backscal(i) + scale);
}

// Compute background scaling factor contribution
for (int i = 0; i < m_on_spec.size(); ++i) {
    double  background = onoff->off_spec()["BACKRESP"][i];
    double  alpha      = onoff->on_spec().backscal(i);
    double  scale      = 0.0;
    if (background > 0.0) {
        scale = alpha * background * exposure;
    }
    else {
        scale = alpha * onoff->off_spec()[i];
    }
    m_on_spec.backscal(i, m_on_spec.backscal(i) + scale);
}

// Compute stacked background scaling factor
for (int i = 0; i < m_on_spec.size(); ++i) {
    double norm = m_off_spec["BACKRESP"][i];
    if (norm > 0.0) {
        double scale = m_on_spec.backscal(i) / norm;
        m_on_spec.backscal(i,scale);
    }
    else {
        double scale = m_on_spec.backscal(i) / m_off_spec[i];
        m_on_spec.backscal(i,scale);
    }
}

        // ... otherwise
        else {

            // Loop over regions
            for (int k = 0; k < off.size(); ++k) {