{{lastupdated_at}} by {{lastupdated_by}} h1. GCTAResponse {{>toc}} h2. Diffuse model response The diffuse model response is implemented by the methods @GCTAResponse::irf_diffuse@ and @GCTAResponse::npred_diffuse@. h3. @GCTAResponse::irf_diffuse@ @GCTAResponse::irf_diffuse@ performs an integration over the true photon arrival direction of the sky map times the instrument response function. The integration is performed in the reference of the observed photon direction. The offset angle (or theta) integration is performed using the helper class @cta_irf_diffuse_kern_theta@ while the azimuth angle (or phi) integration is performed using the helper class @cta_irf_diffuse_kern_phi@. Some analysis has been done to estimate the required precision of the integration. For this purpose, a 180000 s observation of Cen A has been simulated using @ctobssim@, with a deadtime correction of 0.889. The ROI radius was set to 7.5 deg, the energy range spanned 0.1-100 TeV. The events were then subsequently binned using @ctbin@ (100 bins of 0.05 times 0.05 deg). The number of events binned by gtbin was 1117428, hence the expected number of events for a integration time of 1800 s should be 11174.28. Below are the model computation results for a 1800 s integration period and the same observation parameters as used for the simulations. The computation was performed using the @test_response_irf@ function in the @test_CTA.cpp@ file. As the test uses event binning, we also studied the impact of the number of energy bins on the result: | *Theta* | *Phi* | *Ebins* | *Events* | *Uncertainty* | *Time* | | 1e-4 | 1e-4 | 5 | 12387.8 | +10.859% | 1m7.4s | | 1e-4 | 1e-4 | 10 | 11325.4 | +1.352% | 2m49.1s | | 1e-4 | 1e-4 | 20 | 11219.2 | +0.402% | 5m24.1s | | 1e-4 | 1e-4 | 30 | 11207.3 | +0.296% | 8m46.6s | | 1e-3 | 1e-3 | 30 | 11207.4 | +0.296% | 3m5.0s | | 1e-2 | 1e-2 | 30 | 11207.2 | +0.295% | 1m43.1s | | 1e-1 | 1e-1 | 30 | 11207.2 | +0.295% | 1m41.9s | Obviously, there is a clear impact of the number of energy bins used on the precision of the model computation, so we stick to 30 energy bins for the evaluation of the integration precision. It is interesting to notice that the uncertainty in the model computation is basically independent of the integration precision. This can be understood by the fact that the integration is done over the size of the PSF, and if the model and effective do not vary drastically over this size, the integration converges quickly. To illustrate that even with a low precision of 1e-2 and 1e-3 satisfactory images are obtained, model maps of the Cen A map are shown below. The FITS files of the model maps with all ten energy layers are attached to the page for reference. | Integration precision = 1e-2 | Integration precision = 1e-3 | | !cena_model_eps2.png! | !cena_model_eps3.png! | For the moment (18/07/2012) we set the precisions to 1e-2 for the Phi and Theta integrations. Reducing the precision to an even smaller value does not seem to improve the performances. h3. @GCTAResponse::npred_diffuse@ @GCTAResponse::npred_diffuse@ performs an integration of the sky map times the point source Npred method over the region of interest. The offset angle integration is performed using the helper class @cta_npred_diffuse_kern_theta@ while the azimuth angle integration is performed using the helper class @cta_npred_diffuse_kern_phi@. Some analysis has been done to estimate the required precision of the integration. For this purpose a 1800 s observation of Cen A has been simulated, with a deadtime correction of 0.889. The ROI radius was set to 5.0 deg, the energy range spanned 0.1-100 TeV. In doing 10 @ctobssim@ runs, the mean/median number of simulated events was 11232/11239. Doing a single @ctobssim@ run of 180000 s, the number of simulated events was 1124003. The expected number of events should thus be 11240. Here now the results obtained from the Npred computation for the same observation parameters. The computation was performed using the @test_response_npred@ function in the @test_CTA.cpp@ file: | *Theta* | *Phi* | *Events* | *Uncertainty* | *Time* | | 1e-6 | 1e-6 | 11240.6 | +0.005% | 6m28s | | 1e-6 | 1e-4 | 11237.9 | -0.019% | 1m4.6s | | *1e-4* | *1e-4* | *11212.3* | *-0.246%* | *0m6.2s* | | 1e-5 | 1e-5 | 11238.8 | -0.011% | 0m49.2s | | 5e-5 | 5e-5 | 11225.2 | -0.132% | 0m11.7s | | 1e-4 | 1e-5 | 11246.2 | +0.055% | 0m16.5s | | 1e-5 | 1e-4 | 11236.2 | -0.034% | 0m15.4s | | 1e-3 | 1e-3 | 11135.2 | -0.932% | 0m1.3s | For the moment (18/07/2012) we set the precisions to 1e-4 for the Phi and Theta integrations. Eventually, we may reduce them further to 1e-3, but we should first see how this globally affects the fit results. Note that the integration over energy which is performed using the @GObservation::npred_spec@ method uses a precision of 1e-6. Decreasing the precision down to 1e-5 results in an estimated number of 11239.4 events while reducing the computing time to 0m35s. We should think about reducing globally the precision of the Romberg integrator (to maybe 1e-5). 5/4/2013: Here some fit results depending on precision. What is show is the fitted pre factor and index. The true pre factor was 1e-7, the true index was 2.100. The number of simulated Cen A events were 2414, the number of background events were 366, so the total number of events is 2780. | *Theta* | *Phi* | *Prefactor* | *Index* | *Npred* | | 1e-5 | 1e-5 | 0.442e-7 +/- 0.168e-7 | 2.067 +/- 0.037 | 2780 | | *1e-4* | *1e-4* | 0.426e-7 +/- 0.162e-7 | 2.065 +/- 0.037 | 2780 | | 1e-3 | 1e-3 | 0.448e-7 +/- 0.171e-7 | 2.070 +/- 0.037 | 2780 | I also check the precision in @GObservation::npred_spec@: | *npred_spec* | *Evaluations* | *Events* | | 1e-6 | 128 | 11212.3 | | 1e-5 | 64 | 11212.4 | | 1e-4 | 32 | 11217.2 | | 1e-3 | 16 | 11388.3 | This indicates that (at least for a power law) a value of 1e-5 is okay. h3. Diffuse model testing h4. Diffuse model only The implementation of the diffuse model response has been tested using a Cen A simulation. The integration time has been set to 1800 s using the default deadc of @ctobssim@. Events have been simulated for an acceptance cone of 7.5 deg around Cen A. The events spanned an energy range of 0.1-100 TeV. A total of 11956 events has been simulated. Using @ctselect@, events have then be selected within a radius of 5 deg. This resulted in 11956 events, i.e. the same number as simulated (apparently, the Cen A sky map is smaller than the 5 deg ROI). The model used for the simulation was the following:
Using @ctlike@, the following results have been obtained: 
2012-07-18T12:55:14: === GOptimizerLM ===
2012-07-18T12:55:14:  Optimized function value ..: 93564.5
2012-07-18T12:55:14:  Absolute precision ........: 1e-06
2012-07-18T12:55:14:  Optimization status .......: converged
2012-07-18T12:55:14:  Number of parameters ......: 5
2012-07-18T12:55:14:  Number of free parameters .: 2
2012-07-18T12:55:14:  Number of iterations ......: 3
2012-07-18T12:55:14:  Lambda ....................: 1e-06
2012-07-18T12:55:14:  Maximum log likelihood ....: -93564.5
2012-07-18T12:55:14:  Npred .....................: 11956
2012-07-18T12:55:14: === GModels ===
2012-07-18T12:55:14:  Number of models ..........: 1
2012-07-18T12:55:14:  Number of parameters ......: 5
2012-07-18T12:55:14: === GModelDiffuseSource ===
2012-07-18T12:55:14:  Name ......................: Cen A lobes
2012-07-18T12:55:14:  Instruments ...............: all
2012-07-18T12:55:14:  Model type ................: "SpatialMap" * "PowerLaw" * "Constant"
2012-07-18T12:55:14:  Number of parameters ......: 5
2012-07-18T12:55:14:  Number of spatial par's ...: 1
2012-07-18T12:55:14:   Prefactor ................: 1 [0.001,1000]  (fixed,scale=1,gradient)
2012-07-18T12:55:14:  Number of spectral par's ..: 3
2012-07-18T12:55:14:   Prefactor ................: 1.8199e-07 +/- 9.56868e-09 [1e-14,0.0001] ph/cm2/s/MeV (free,scale=1e-07,gradient)
2012-07-18T12:55:14:   Index ....................: -2.10617 +/- 0.00565118 [-5,5]  (free,scale=1,gradient)
2012-07-18T12:55:14:   PivotEnergy ..............: 100 [10,1e+06] MeV (fixed,scale=1,gradient)
2012-07-18T12:55:14:  Number of temporal par's ..: 1
2012-07-18T12:55:14:   Constant .................: 1 (relative value) (fixed,scale=1,gradient)
As expected for a good fit, Npred is equal to the observed number of events. The Prefactor is at +0.94 sigma above the true value, the index is at -1.09 sigma below the true value. This offset is reasonable. As next step the events have been binned using @ctbin@ (30 energy slices from 0.1-100 TeV, 100 x 100 pixels with binsize of 0.05 pixels). This resulted in 11884 events binned into the counts map. @ctlike@ was then employed to fit the counts map. The following results have been obtained: 
2012-07-18T13:14:15: === GOptimizerLM ===
2012-07-18T13:14:15:  Optimized function value ..: 24976.5
2012-07-18T13:14:15:  Absolute precision ........: 1e-06
2012-07-18T13:14:15:  Optimization status .......: converged
2012-07-18T13:14:15:  Number of parameters ......: 5
2012-07-18T13:14:15:  Number of free parameters .: 2
2012-07-18T13:14:15:  Number of iterations ......: 3
2012-07-18T13:14:15:  Lambda ....................: 1e-06
2012-07-18T13:14:15:  Maximum log likelihood ....: -24976.5
2012-07-18T13:14:15:  Npred .....................: 11884
2012-07-18T13:14:15: === GModels ===
2012-07-18T13:14:15:  Number of models ..........: 1
2012-07-18T13:14:15:  Number of parameters ......: 5
2012-07-18T13:14:15: === GModelDiffuseSource ===
2012-07-18T13:14:15:  Name ......................: Cen A lobes
2012-07-18T13:14:15:  Instruments ...............: all
2012-07-18T13:14:15:  Model type ................: "SpatialMap" * "PowerLaw" * "Constant"
2012-07-18T13:14:15:  Number of parameters ......: 5
2012-07-18T13:14:15:  Number of spatial par's ...: 1
2012-07-18T13:14:15:   Prefactor ................: 1 [0.001,1000]  (fixed,scale=1,gradient)
2012-07-18T13:14:15:  Number of spectral par's ..: 3
2012-07-18T13:14:15:   Prefactor ................: 1.82128e-07 +/- 9.61141e-09 [1e-14,0.0001] ph/cm2/s/MeV (free,scale=1e-07,gradient)
2012-07-18T13:14:15:   Index ....................: -2.10646 +/- 0.0056721 [-5,5]  (free,scale=1,gradient)
2012-07-18T13:14:15:   PivotEnergy ..............: 100 [10,1e+06] MeV (fixed,scale=1,gradient)
2012-07-18T13:14:15:  Number of temporal par's ..: 1
2012-07-18T13:14:15:   Constant .................: 1 (relative value) (fixed,scale=1,gradient)
Again, Npred is equal to the observed number of events. The Prefactor and index are very similar (but not identical) to the results obtained using an unbinned analysis. h4. Diffuse and background model Now I fitted a diffuse model on top of the background model. The XML file used for the model is shown below:
There were 581 counts for Cen A and 3250 counts in the background, giving a total of 3831 counts. The @ctselect@ step left again the total number of counts unchanged. The (unbinned) @ctlike@ run gave the following result: 
2012-07-18T13:31:51: === GOptimizerLM ===
2012-07-18T13:31:51:  Optimized function value ..: 34312.9
2012-07-18T13:31:51:  Absolute precision ........: 1e-06
2012-07-18T13:31:51:  Optimization status .......: converged
2012-07-18T13:31:51:  Number of parameters ......: 10
2012-07-18T13:31:51:  Number of free parameters .: 5
2012-07-18T13:31:51:  Number of iterations ......: 4
2012-07-18T13:31:51:  Lambda ....................: 1e-07
2012-07-18T13:31:51:  Maximum log likelihood ....: -34312.9
2012-07-18T13:31:51:  Npred .....................: 3831
2012-07-18T13:31:51: === GModels ===
2012-07-18T13:31:51:  Number of models ..........: 2
2012-07-18T13:31:51:  Number of parameters ......: 10
2012-07-18T13:31:51: === GModelDiffuseSource ===
2012-07-18T13:31:51:  Name ......................: Cen A lobes
2012-07-18T13:31:51:  Instruments ...............: all
2012-07-18T13:31:51:  Model type ................: "SpatialMap" * "PowerLaw" * "Constant"
2012-07-18T13:31:51:  Number of parameters ......: 5
2012-07-18T13:31:51:  Number of spatial par's ...: 1
2012-07-18T13:31:51:   Prefactor ................: 1 [0.001,1000]  (fixed,scale=1,gradient)
2012-07-18T13:31:51:  Number of spectral par's ..: 3
2012-07-18T13:31:51:   Prefactor ................: 5.49351e-16 +/- 4.03622e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2012-07-18T13:31:51:   Index ....................: -2.52938 +/- 0.0527701 [-0,-5]  (free,scale=-1,gradient)
2012-07-18T13:31:51:   PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2012-07-18T13:31:51:  Number of temporal par's ..: 1
2012-07-18T13:31:51:   Constant .................: 1 (relative value) (fixed,scale=1,gradient)
2012-07-18T13:31:51: === GCTAModelRadialAcceptance ===
2012-07-18T13:31:51:  Name ......................: Background
2012-07-18T13:31:51:  Instruments ...............: CTA
2012-07-18T13:31:51:  Model type ................: "Gaussian" * "PowerLaw" * "Constant"
2012-07-18T13:31:51:  Number of parameters ......: 5
2012-07-18T13:31:51:  Number of radial par's ....: 1
2012-07-18T13:31:51:   Sigma ....................: 3.00231 +/- 0.0376892 [0.01,10] deg2 (free,scale=1,gradient)
2012-07-18T13:31:51:  Number of spectral par's ..: 3
2012-07-18T13:31:51:   Prefactor ................: 6.56254e-05 +/- 2.06902e-06 [0,0.001] ph/cm2/s/MeV (free,scale=1e-06,gradient)
2012-07-18T13:31:51:   Index ....................: -1.82475 +/- 0.017447 [-0,-5]  (free,scale=-1,gradient)
2012-07-18T13:31:51:   PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2012-07-18T13:31:51:  Number of temporal par's ..: 1
2012-07-18T13:31:51:   Constant .................: 1 (relative value) (fixed,scale=1,gradient)
Also here, Npred corresponds to the observed number of counts. The prefactor is at -0.51 sigma and the index is at -0.94 sigma from the expected model values, which is pretty reasonable.