ctools

The actual code in GCTAOnOffObservation::compute_arf is:

            // Initialize totals
            double totsolid = 0.0;
            double totpsf   = 0.0;

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

                // Get pointer on sky region map
                const GSkyRegionMap* on_map = static_cast<const GSkyRegionMap*>(on[k]);

                // Loop over pixels in On region map and integrate effective
                // area
                for (int j = 0; j < on_map->nonzero_indices().size(); ++j) {

                    // Get pixel index
                    int pixidx = on_map->nonzero_indices()[j];

                    // Get direction to pixel center
                    GSkyDir pixdir = on_map->map().inx2dir(pixidx);

                    // Get solid angle subtended by this pixel
                    double pixsolid = on_map->map().solidangle(pixidx);

                    // Compute position of pixel centre in instrument coordinates
                    double theta = obsdir.dist(pixdir);
                    double phi   = obsdir.posang(pixdir);

                    // Add up effective area
                    m_arf[i] += response->aeff(theta,
                                               phi,
                                               zenith,
                                               azimuth,
                                               logEtrue) * pixsolid;

                    // Add pixel solid angle to total for averaging later
                    totsolid += pixsolid;

                    // Compute offset angle to source
                    double delta = srcdir.dist(pixdir);

                    // Integrate PSF
                    totpsf += response->psf(delta,
                                            theta,
                                            phi,
                                            zenith,
                                            azimuth,
                                            logEtrue) * pixsolid;

                } // endfor: looped over all pixels in region map

            } // endfor: looped over all regions

            // Average effective area over solid angle
            if (totsolid > 0.0) {
                m_arf[i] /= totsolid;
            }

            // Correct effective area by containment fraction
            if (totpsf >= 0.0 && totpsf <= 1.0) {
                m_arf[i] *= totpsf;
            }

This is not accurate as the effective area should in fact be weighted by the PSF, meaning that the regions where the PSF is large count more. Here the respective code (which is in fact simpler):

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

                // Get pointer on sky region map
                const GSkyRegionMap* on_map = static_cast<const GSkyRegionMap*>(on[k]);

                // Loop over pixels in On region map and integrate effective
                // area
                for (int j = 0; j < on_map->nonzero_indices().size(); ++j) {

                    // Get pixel index
                    int pixidx = on_map->nonzero_indices()[j];

                    // Get direction to pixel center
                    GSkyDir pixdir = on_map->map().inx2dir(pixidx);

                    // Get solid angle subtended by this pixel
                    double pixsolid = on_map->map().solidangle(pixidx);

                    // Compute position of pixel centre in instrument coordinates
                    double theta = obsdir.dist(pixdir);
                    double phi   = obsdir.posang(pixdir);

                    // Compute offset angle to source
                    double delta = srcdir.dist(pixdir);

                    // Get PSF value times the solid angle
                    double psf = response->psf(delta,
                                               theta,
                                               phi,
                                               zenith,
                                               azimuth,
                                               logEtrue) * pixsolid;

                    // Add up effective area
                    m_arf[i] += response->aeff(theta,
                                               phi,
                                               zenith,
                                               azimuth,
                                               logEtrue) * psf;

                } // endfor: looped over all pixels in region map

            } // endfor: looped over all regions

Below the comparison between the old computation (top) and the new computation (bottom). The pull distributions were generated for 10 reconstructed energy bins between 0.1 and 100 TeV, a on region radius of 0.2 deg and a deadtime correction factor of 1.0.

The mean of the Prefactor is now closer to zero.

Below the corresponding histograms after 100 pulls:

I checked that the value of deadc has no impact on the pull distribution by running 100 pulls for deadc=0.5. The pull histograms and pull evolutions are shown below. Everything looks ok.

I also checked the impact of the On region radius. Before I used 0.2 deg, below results for 0.1 deg (top) and 0.3 deg (bottom). Also here everything is okay.

I also tested the larger energy range from 50 GeV - 100 TeV with 40 energy bins instead of 10. Also looks good here.

I redid the same energy binning but now using WSTAT as fit statistic. The command was

$ cspull.py inobs=obs.xml statistic=WSTAT logfile=cspull_newarf_wstat_0.05-100_b40.log

The result is shown below. Obviously, the spread of the pull distributions is too large. I’m wondering whether the issue arises tue to many bins with Non and/or Noff = 0.

Here the WSTAT result with 10 energy bins between 100 GeV and 100 TeV. Results are still not as expected. Bias is small, but width is too large.

I consider this issue as fixed, but opened issue #2266 to follow-up on the WSTAT problem.