Commit 02b0b1c2158e5031d00bf2844990e3532b84f1cf

Authored by Michael Mayer
1 parent cd6bebbf

Added simulation part of the documentation

doc/source/user_manual/iact_analysis/analysis.rst
... ... @@ -74,10 +74,10 @@ For this purpose one can simply pass an invalid coordinate to :ref:`csfindobs`.
74 74  
75 75 The selection by region will simply be omitted.
76 76  
77   -Note
78   -^^^^
79   -On default, :ref:`csfindobs` only select data of highest quality (i.e. QUALITY=0). Specifying the hidden parameter e.g. ``min_qual=1``
80   -allows to select all data with looser quality criteria.
  77 +.. note::
  78 +
  79 + On default, :ref:`csfindobs` only select data of highest quality (i.e. QUALITY=0). Specifying the hidden parameter e.g. ``min_qual=1``
  80 + allows to select all data with looser quality criteria.
81 81  
82 82 Create an observation list
83 83 --------------------------
... ...
doc/source/user_manual/iact_analysis/high_level.rst
... ... @@ -219,7 +219,7 @@ Example unbinned observation container
219 219 Size of the X axis in pixels [100]
220 220 Size of the Y axis in pixels [100]
221 221 Pixel size (deg/pixel) [0.02]
222   - Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUB]
  222 + Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUBDIV]
223 223 Output residual map file [resmap.fits]
224 224  
225 225 Example binned/stacked observation container
... ... @@ -235,7 +235,7 @@ The tool behaves differently if the input observation parameter is a count cube.
235 235 Input PSF cube file (only needed for stacked analysis) [psfcube.fits]
236 236 Input background cube file (only needed for stacked analysis) [bkgcube.fits]
237 237 Input model XML file [binned_results.xml]
238   - Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUB]
  238 + Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUBDIV]
239 239 Output residual map file [resmap.fits]
240 240  
241 241 Note that the tool also queries for an input model cube file (which we set to ``NONE`` here). This is very convenient in case the model cube has already been precomputed
... ...
doc/source/user_manual/iact_analysis/simulation.rst
... ... @@ -2,29 +2,299 @@
2 2  
3 3 Simulating IACT data
4 4 ====================
5   -
6   -
  5 +This section describes how to run simulations with current IACT data. In particular, a guide how to simulate data for a given set of IRFs, or a given set of
  6 +observations is provided. In addition, tools to compute sensitivity curves or pull distributions are presented.
  7 +All tools are based on the event list simulator :ref:`ctobssim`.
7 8  
8 9 Simulations for a specific observation container
9 10 ------------------------------------------------
  11 +On default, :ref:`ctobssim` builds a new observation container by querying relevant user parameters. In case
  12 +an observation XML container is already in place, the hidden parameter ``inobs`` can be specified (:ref:`see here how to create an
  13 +observation container <sec_analysis>`). It might be useful to run :ref:`ctselect` beforehand to ensure that the simulation knows about
  14 +energy thresholds, field of view, etc. If an unselected observation container is used, the simulation tool will query for global energy
  15 +boundaries and a global RoI radius. The following example illustrates the usage of an already selected observation container.
  16 +
  17 +.. code-block:: bash
  18 +
  19 + $ ctobssim inobs=selected_obs.xml
  20 + Input model XML file [crab_models.xml]
  21 + Output event data file or observation definition XML file [sim_obs.xml]
  22 +
  23 +Note that :ref:`ctobssim` writes a new event list file per observation into the current directory. If the observation container
  24 +is large, it might be convenient to create a new folder and specify the ``prefix`` keyword (similar to :ref:`ctselect`):
  25 +
  26 +.. code-block:: bash
  27 +
  28 + $ mkdir simulated
  29 + $ ctobssim inobs=selected_obs.xml prefix=simulated/sim_events_
  30 + ...
  31 +
  32 +The newly created observation container can be used in the same way as a normal observation container, i.e. it can be analysed with
  33 +:ref:`ctlike` etc.
  34 +
  35 +.. note::
  36 +
  37 + On default, :ref:`ctobssim` uses the same seed in every execution, i.e. the tool will always return the same simulated event list, regardless
  38 + when, where and how it was executed. To steer this the user can provide the hidden integer parameter ``seed``:
  39 +
  40 + .. code-block:: bash
  41 +
  42 + $ ctobssim inobs=selected_obs.xml seed=2
  43 +
  44 + This might be important if simulations are run in parallel.
10 45  
11 46 Simulating observation time for a given set of IRFs
12 47 ---------------------------------------------------
  48 +There are other use cases to run simulations, for instance when there is no observation container available yet. This might be the case e.g
  49 +for an observation proposal: one would want to use a certain set of IRFs to simulate e.g. 50 hours of data to check whether a source
  50 +of interest can be detected. For this purpose it is necessary to first find a set of suitable IRFs to be used for
  51 +simulations. This can be easily done using :ref:`csfindobs` and using the ``expression`` parameter. The following example extracts
  52 +observations with at a certain zenith angle range (20-25 degrees) and azimuth angle range (90-100 degrees):
  53 +
  54 +.. code-block:: bash
  55 +
  56 + $ csfindobs expression="AZ_PNT>90.0&&AZ_PNT<100.0&&ZEN_PNT>20.0&&ZEN_PNT<25.0"
  57 + Name of FITS production (Run csiactdata to view your options) [fits-prod-name]
  58 + Right ascension [83.6331] NONE
  59 + Runlist outfile [zen20_az90.lis]
  60 +
  61 +Of course, further details could be specified, e.g. one might want to check for a certain range of muon efficieny, or any other parameter which impacts IRFs.
  62 +The output runlist subsequently has to be converted to an observation XML file using :ref:`csiactobs` (:ref:`see here how to execute this tool <sec_analysis>`).
  63 +In the following it is assumed that the observation XML file of the selected representative observations is produced and called ``myobs.xml``.
13 64  
14 65 Create caldb entry
15 66 ^^^^^^^^^^^^^^^^^^
  67 +Now, one of the observations in this container have to be added to the local calibration database. This database is used
  68 +by ctools to extract IRFs for simulations.
  69 +
  70 +.. code-block:: bash
  71 +
  72 + $ csobs2caldb
  73 + Input observation definition file [myobs.xml]
  74 + Response output name (e.g. Zenith50) [zen20_az90]
  75 +
  76 +Note that the response name has been specified to "zen20_az90". Any name can be chosen which represents the observation which was inserted into the database.
  77 +On default, :ref:`csobs2caldb` uses the first observation in the provided XML file. This can be modifyied by specifying the hidden integer parameter index, e.g.:
  78 +
  79 +.. code-block:: bash
  80 +
  81 + $ csobs2caldb index=1
  82 +
  83 +This would use the second entry instead of the first from the observation container to extract the IRFs.
  84 +
  85 +.. note::
  86 +
  87 + The ``index`` parameter starts counting from zero.
  88 +
  89 +Check caldb entries
  90 +^^^^^^^^^^^^^^^^^^^
  91 +In order to check if the calibration database has successfully been added, the tool :ref:`cscaldb` can be used to inspect the calibration database:
  92 +
  93 +.. code-block:: bash
  94 +
  95 + $ cscaldb debug=yes
  96 + ...
  97 + 2016-02-25T15:57:11: +==============+
  98 + 2016-02-25T15:57:11: | Mission: cta |
  99 + 2016-02-25T15:57:11: +==============+
  100 + 2016-02-25T15:57:11: === Response functions in database "hess" ===
  101 + 2016-02-25T15:57:11: zen20_az90
  102 + 2016-02-25T15:57:11:
  103 + 2016-02-25T15:57:11: === Response functions in database "prod2" ===
  104 + 2016-02-25T15:57:11: North_0.5h
  105 + 2016-02-25T15:57:11: North_50h
  106 + 2016-02-25T15:57:11: North_5h
  107 + 2016-02-25T15:57:11: South_0.5h
  108 + 2016-02-25T15:57:11: South_50h
  109 + 2016-02-25T15:57:11: South_5h
  110 + ...
  111 +
  112 +Note that the observation of the actual instrument ("hess" in this example) was added to the CTA mission. This is mainly for structural and simplicity reasons
  113 +and might be improved eventually. Nevertheless, this caldb entry can now be used to simulate any observation duration.
16 114  
17 115 Simulate
18 116 ^^^^^^^^
19 117  
  118 +.. code-block:: bash
  119 +
  120 + $ ctobssim
  121 + RA of pointing (degrees) (0-360) [83.63]
  122 + Dec of pointing (degrees) (-90-90) [22.01]
  123 + Radius of FOV (degrees) (0-180) [2.5]
  124 + Start time (MET in s) [0.0]
  125 + End time (MET in s) [18000.0]
  126 + Lower energy limit (TeV) [0.5]
  127 + Upper energy limit (TeV) [50]
  128 + Calibration database [hess]
  129 + Instrument response function [zen20_az90]
  130 + Input model XML file [input_model.xml]
  131 + Output event data file or observation definition XML file [simulated_events.fits]
  132 +
  133 +Here, a total duration of 5 hours at the crab position was simulated using the newly created caldb entry. Note that the choice of the energy range is up to the user.
  134 +One should bear in mind that simulations only make sense in energy ranges where the IRFs are well-defined. The simulation
  135 +assumes one observation of defined duration here. Accordingly, only one FITS event list is written as output. The input XML model
  136 +must now contain only one background model and (of course) all sky models that are needed. Here is an example of an input model
  137 +e.g. for the Crab Nebula:
  138 +
  139 +.. code-block:: xml
  140 +
  141 + <?xml version="1.0" standalone="no"?>
  142 + <source_library title="source library">
  143 + <source name="Crab" type="PointSource">
  144 + <spectrum type="PowerLaw">
  145 + <parameter name="Prefactor" scale="1e-17" value="3.5" min="1e-07" max="1000.0" free="1"/>
  146 + <parameter name="Index" scale="-1" value="2.5" min="0.0" max="+5.0" free="1"/>
  147 + <parameter name="Scale" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
  148 + </spectrum>
  149 + <spatialModel type="SkyDirFunction">
  150 + <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
  151 + <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/>
  152 + </spatialModel>
  153 + </source>
  154 + <source name="CTABackgroundModel" type="CTAIrfBackground" instrument="CTA">
  155 + <spectrum type="PowerLaw">
  156 + <parameter name="Prefactor" scale="1.0" value="1.0" min="1e-3" max="1e+3" free="1"/>
  157 + <parameter name="Index" scale="1.0" value="0.0" min="-5.0" max="+5.0" free="1"/>
  158 + <parameter name="Scale" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
  159 + </spectrum>
  160 + </source>
  161 + </source_library>
  162 +
  163 +Note that since the caldb entry (:ref:`csobs2caldb`) has been produced under the mission "CTA", the event list is now considered a CTA observation (regardless of the input instrument).
  164 +Therefore, it is important to have "CTA" specified in the instrument attribute of the background model XML component. Of course the addition of any sky model can be
  165 +done without a problem.
  166 +
  167 +The output of this kind of simulation is a single FITS event list file that can be used for further analysis.
  168 +For instance one could run :ref:`ctlike` in the following way:
  169 +
  170 +.. code-block:: bash
  171 +
  172 + $ ctlike
  173 + Input event list, counts cube or observation definition XML file [simulated_events.fits]
  174 + Calibration database [hess]
  175 + Instrument response function [zen20_az90]
  176 + Input model XML file [input_model.xml]
  177 + Output model XML file [simulated_fit_results.xml]
  178 +
  179 +Note that :ref:`ctlike` automatically queries for the IRF database in case no full observation container is provided. Any other tool taking
  180 +unbinned observations as input will behave the same way. Now all means to run simulations with current IACTs are at hand.
  181 +The following sections will describe tools for running higher level simulation analysis.
  182 +
20 183 Create pull distributions
21 184 -------------------------
  185 +Pull distributions inspect the capability of the science tools to find source parameters in the data that were injected into the simulations.
  186 +The tool :ref:`cspull` will run several times the sequence of simulating and fitting the data. After each sequence, the pull value of the
  187 +fitted parameters is written into an ASCII file. The pull value ``p`` of parameter ``X`` is defined as follows:
  188 +
  189 +.. math::
  190 + p = \frac{X_{\rm fitted} - X_{\rm true}}{\Delta X_{\rm fitted}}
  191 +
  192 +Where :math:`\Delta X_{\rm fitted}` is the derived error of the fitted parameter :math:`X_{\rm fitted}`. After several trials, the pull value should follow
  193 +a normal distribution. Knowing this, the tool may help to distinguish if a source is detectable and can be normally reconstructed.
  194 +The tool can either work on an input observation container or on the caldb entry that was created above. The number of trials is queried by the tool, too.
  195 +
  196 +Example: observation container
  197 +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  198 +
  199 +.. code-block:: bash
  200 +
  201 + $ cspull inobs=selected_obs.xml
  202 + Number of energy bins (0=unbinned) [0]
  203 + Input model XML file [crab_results.xml]
  204 + Output pull distribution file [pull.dat]
  205 + Number of trials [100]
  206 +
  207 +Note that the tool asks for the energy binning. It will use an unbinned analysis if zero bins are provided. Otherwise
  208 +it will query for binning parameters and will convert the simulated observations into a binned/stacked observation by running
  209 +:ref:`ctbin`, :ref:`ctexpcube`, :ref:`ctpsfcube` and :ref:`ctbkgcube` in order before fitting the simulated data with :ref:`ctlike`.
  210 +
  211 +Example: caldb entry
  212 +^^^^^^^^^^^^^^^^^^^^
  213 +
  214 +.. code-block:: bash
  215 +
  216 + $ cspull
  217 + RA of pointing (deg) (0-360) [83.6331]
  218 + Dec of pointing (deg) (-90-90) [22.0145]
  219 + Duration (in s) [18000.0]
  220 + Lower energy limit (TeV) [0.5]
  221 + Upper energy limit (TeV) [50]
  222 + Calibration database [hess]
  223 + Instrument response function [zen20_az90]
  224 + Number of energy bins (0=unbinned) [0]
  225 + Input model XML file [input_models.xml]
  226 + Output pull distribution file [pull.dat]
  227 + Number of trials [100]
  228 +
  229 +Again, the caldb entry ``zen20_az90`` which was created above is used to simulate 5 hours of data. The data is simulated and analysed 100
  230 +times in a row and the output is written into the file ``pull.dat``.
22 231  
23 232 Visualise pull distributions
24 233 ----------------------------
  234 +The output ASCII file of :ref:`cspull` can be inspected by the user using own scripts. As a starting point, the ``example`` folder
  235 +contains a script that can be used to visualise the pull distribution histogram:
  236 +
  237 +.. code-block:: bash
  238 +
  239 + $ python $CTOOLS/examples/show_pull_histogram.py pull.dat Pull_Crab_Prefactor 50
25 240  
26   -Create sensitiviy curve
27   ------------------------
  241 +Run the script without arguments to see the usage. In case a wrong parameter name was provided, the tool will print available parameters on the screen.
  242 +The last integer parameter defines the binning of the histogram. In addition, the tool shows the standard normal distribution which should
  243 +match the pull distribution for a high number of trials if everything went well.
  244 +
  245 +Create sensitivity curve
  246 +------------------------
  247 +A very important means to illustrate if a source can be detected with a certain significance is a sensitivity curve.
  248 +The tool :ref:`cssens` can be used to compute either an integral, or a differential sensitivity curve for a given observation container.
  249 +Analogous to :ref:`cspull`, the tool works either on an unbinned observation container or using a caldb entry.
  250 +The script :ref:`cssens` computes the differential or integrated CTA sensitivity using maximum likelihood fitting of a test source.
  251 +The differential sensitivity is determined for a number of energy bins, the integral sensitivity is determined for a number
  252 +of energy thresholds. The test source is fitted to simulated data using ctlike to determine itโ€™s detection significance as a
  253 +function of source flux. The source flux is then varied until the source significance achieves a given level, specified by the
  254 +(hidden) significance parameter ``sigma``. To damp variations between individual Monte Carlo simulations, a sliding average is
  255 +applied in the significance computation (controlled by hidden ``num_avg parameter``). Note that this procedure might take a while to compute.
  256 +To switch from diffferential to integral sensitivity, the hidden parameter ``type="Integral"`` has to be specified. It should be noted that
  257 +the results are written to a file ``sensitivity.dat``. To change the output file name, the hidden parameter ``outfile=mysensitivity.dat``
  258 +can be provided.
  259 +
  260 +
  261 +Example: observation container
  262 +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  263 +
  264 +.. code-block:: bash
  265 +
  266 + $ cssens inobs=selected_obs.xml
  267 + Input model XML file [crab_results.xml]
  268 + Source name [Crab]
  269 + Lower energy limit (TeV) [0.5]
  270 + Upper energy limit (TeV) [50]
  271 + Number of energy bins for differential sensitivity computation [21]
  272 +
  273 +Example: caldb entry
  274 +^^^^^^^^^^^^^^^^^^^^
  275 +
  276 +.. code-block:: bash
  277 +
  278 + $ cssens
  279 + Calibration database [hess]
  280 + Instrument response function [zen20_az90]
  281 + Effective exposure time (s) [18000.0]
  282 + Radius of ROI (deg) [2.5]
  283 + Input model XML file [input_models.xml]
  284 + Source name [Crab]
  285 + Lower energy limit (TeV) [0.5]
  286 + Upper energy limit (TeV) [50]
  287 + Number of energy bins for differential sensitivity computation [21]
  288 +
  289 +Note that the input model in the latter case only requires to contain one background model (see section above). Using this setup,
  290 +the user can specify the observation time for the sensitivity curve. In this example, again, 5 hours are used.
28 291  
29 292 Visualise sensitivity curves
30   -----------------------------
31 293 \ No newline at end of file
  294 +----------------------------
  295 +In order to show a sensitivity curve, there is also a script in the example folder:
  296 +
  297 +.. code-block:: bash
  298 +
  299 + $ python $CTOOLS/examples/show_sensitivity.py sensitivity.dat
  300 +
  301 +This will display a sensitivity curve for the applied settings.
... ...
doc/source/user_manual/iact_analysis/tips.rst
... ... @@ -2,39 +2,304 @@
2 2  
3 3 Tips, Tricks and FAQs
4 4 =====================
  5 +This section shows some of the wrinkles that could fascilitate the users' life when running analyses.
5 6  
6 7 Visualise observations
7 8 ----------------------
  9 +A very common task is to inspect the observation container and visualise several quantities.
  10 +
  11 +Create a skymap of counts
  12 +^^^^^^^^^^^^^^^^^^^^^^^^^
  13 +The tool :ref:`ctskymap` allows to create a single map of counts contained inside an observation container.
  14 +In contrast to :ref:`ctbin`, ctskymap takes into account all events that fall into the skymap energy range and spatial range.
  15 +Since :ref:`ctbin` is indended for binned/stacked analysis, it respects energy thresholds, RoI selections etc.
  16 +Accordingly, :ref:`ctskymap` should be used for a simple visualisation of the counts.
  17 +
  18 +
  19 +.. code-block:: bash
  20 +
  21 + $ ctskymap debug=yes
  22 + Input event list or observation definition XML file [events.fits] obs.xml
  23 + Lower energy limit (TeV) [0.1]
  24 + Upper energy limit (TeV) [100.0]
  25 + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63]
  26 + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01]
  27 + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR]
  28 + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
  29 + Image scale (in degrees/pixel) [0.02]
  30 + Size of the X axis in pixels [200]
  31 + Size of the Y axis in pixels [200]
  32 + Output skymap file [skymap.fits]
  33 +
  34 +Create region file of pointings
  35 +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  36 +In particular, one might want to know how the pointings of the observation container distribute across the sky.
  37 +For this purpose, :ref:`csobsinfo` provides the hidden parameter ``ds9file`` which signals to produce a region file that
  38 +contains these information:
  39 +
  40 +.. code-block:: bash
  41 +
  42 + $ csobsinfo ds9file=pointings.reg
  43 + Parfile "csobsinfo.par" not found. Create default parfile.
  44 + Event list, counts cube, or observation definition file [obs.xml]
  45 +
  46 +Together with the previously produced skymap, one can plot the pointings using e.g. ds9:
  47 +
  48 +.. code-block:: bash
  49 +
  50 + $ ds9 skymap.fits -regions pointings.reg
8 51  
9 52 Plot zenith angle distribution of observations
10 53 ----------------------------------------------
  54 +The tool :ref:`csobsinfo` dumps more information into the log file that might be of interest to catch in a script.
  55 +Therefore, one could simply use this tool via python and plot some distributions:
  56 +
  57 +.. code-block:: python
  58 +
  59 + # Run script
  60 + import cscripts
  61 + info = cscripts.csobsinfo()
  62 + info["inobs"] = "selected_obs.xml"
  63 + info.run()
  64 +
  65 + # Plot data
  66 + import matplotlib.pyplot as plt
  67 + plt.hist(info.zeniths(), bins=10, range=[0,90])
  68 + plt.show()
  69 +
  70 +This will plot a histogram showing the zenith angle distribution of the observation container. In addition there are more methods
  71 +to access information in python:
  72 +
  73 +.. code-block:: python
  74 +
  75 + # csobsinfo can also return:
  76 + info.azimuths() # list of azimuth values
  77 + info.offsets() # list of offset values (only computed if hidden parameter offset=yes was specified)
  78 + info.ebounds() # gammalib.GEbounds object of energy ranges
  79 + info.gti() # gammalib.GGti object containing good time intervals
  80 +
11 81  
12 82 Visualise models
13 83 ----------------
  84 +Similar to :ref:`csobsinfo`, there is :ref:`csmodelinfo` to find out more about the content of a model XML file. In order to show the
  85 +position and sizes of the model on top of a skymap, this tool has the hidden parameter ``ds9file``, too.
  86 +
  87 +.. code-block:: bash
  88 +
  89 + csmodelinfo ds9file=models.reg
  90 + Input model XML file [$CTOOLS/share/models/crab.xml]
  91 +
  92 +There are several options regarding the color, text and other attributes of the region file. To see a full list,
  93 +:ref:`visit the reference page <csmodelinfo>`.
  94 +Analogous to the pointings of the observation container, the models can be visualised using e.g. ds9:
  95 +
  96 +.. code-block:: bash
  97 +
  98 + $ ds9 skymap.fits -regions models.reg
  99 +
  100 +Manipulating models in python
  101 +-----------------------------
  102 +This step will give some guidance how to work with model XML files in python and how to manipulate their content
  103 +
  104 +.. code-block:: python
  105 +
  106 + import gammalib
  107 +
  108 + # Open model file
  109 + models = gammalib.GModels("$CTOOLS/share/models/crab.xml")
  110 +
  111 + # Access a model component
  112 + src = models["Crab"]
  113 + print(src)
  114 +
  115 + # Retrieve the spectral or spatial component
  116 + print(src.spectral())
  117 + print(src.spatial())
  118 +
  119 + # Print the spectral parameter "Prefactor"
  120 + prefactor = src.spectral()["Prefactor"]
  121 + print("Prefactor:",prefactor.value(),"+-",prefactor.error())
  122 +
  123 + # Set the prefactor value
  124 + prefactor.value(3.5e-17)
  125 +
  126 + # Loop over models and fix all parameters
  127 + for model in models:
  128 +
  129 + # Loop over model parameters
  130 + for par in model:
  131 +
  132 + # Fix parameter
  133 + par.fix()
  134 +
  135 + # Release one specific parameter
  136 + models["Crab"]["Prefactor"].free()
  137 +
  138 + # Set the parameter range
  139 + models["Crab"]["Prefactor"].min(1e-18)
  140 + models["Crab"]["Prefactor"].max(1e-16)
  141 +
  142 + # ... or in one step:
  143 + models["Crab"]["Prefactor].range(1e-18, 1e-16)
  144 +
  145 + # Remove model from container
  146 + models.remove("Crab")
  147 +
  148 + # Save model to another XML file
  149 + models.save("mymodels.xml")
  150 +
  151 +
  152 +Retrieve likelihood values from :ref:`ctlike`
  153 +---------------------------------------------
  154 +For some purposes, it might be useful to retrieve the fitted likelihood value and other results of the fit with :ref:`ctlike`.
  155 +
  156 +.. code-block:: python
  157 +
  158 + import gammalib
  159 + import ctools
  160 +
  161 + # Create ctlike
  162 + like = ctools.ctlike()
  163 + like["inobs"] = "selected_obs.xml"
  164 + like["inmodel"] = "mymodels.xml"
  165 + like.run()
  166 +
  167 + # Get optimizer from ctlike
  168 + opt = like.opt()
  169 + print(opt)
  170 +
  171 + # Get likelihood value from optimizer
  172 + print(opt.value())
  173 +
  174 + # Get fit status:
  175 + print(opt.status())
  176 + # 0: converged
  177 + # 1: stalled
  178 + # 2: singular curvature matrix encountered
  179 + # 3: curvature matrix not positive definite
  180 + # 4: errors are inaccurate
14 181  
15 182 Speed up analysis
16 183 -----------------
  184 +In some cases it may occur that the analysis takes very long. There are several reasons that can slow down the analysis quite severly.
  185 +Not all issues can be tackled. Nevertheless, here is a list of actions that have proven to speed up the fit.
17 186  
18   -Create count map with one energy bin
19   -------------------------------------
  187 +* Consider switching to binned analysis if the observation time exceeds 50-100 hours.
  188 +* Diffuse and extended sources in the RoI are computational expensive. However, in unbinned analysis they impact the speed much more.
  189 +* Usage of energy dispersion in unbinned analysis could also cause a long fit. Try switching it off if speed is important.
  190 +* Switch to a machine that supports OpenMP. This allows to parallelise the fit onto several cores.
20 191  
21   -Compute excess map
22   -------------------
  192 +Compute excess maps
  193 +-------------------
  194 +A very important means to visualise images of sources are excess maps. The tool :ref:`csresmap` can also be used to create such a map.
  195 +The user, however, has to manually remove the source(s) of interest from the input model XML file.
  196 +
  197 +.. code-block:: python
  198 +
  199 + import gammalib
  200 + models = gammalib.GModels("crab_models.xml")
  201 + models.remove("Crab")
  202 + models.save("models_without_crab.xml")
  203 +
  204 +Of course, this can also be done by editing the XML file with an editor.
  205 +
  206 +Subsequently, :ref:`csresmap` can be executed using the ``algorithm=SUB`` parameter:
  207 +
  208 +.. code-block:: bash
  209 +
  210 + $ csresmap
  211 + Input event list, counts cube, or observation definition XML file [selected_obs.xml]
  212 + Input model XML file [models_without_crab.xml]
  213 + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63]
  214 + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01]
  215 + Coordinate System (CEL|GAL) [CEL]
  216 + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR]
  217 + Size of the X axis in pixels [100]
  218 + Size of the Y axis in pixels [100]
  219 + Pixel size (deg/pixel) [0.02]
  220 + Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUB]
  221 + Output residual map file [excessmap.fits]
23 222  
24   -Split TS map computation
25   -------------------------
26 223  
27 224 Compute flux maps
28 225 -----------------
  226 +The creation of a flux map is much more expensive in terms of computing than an excess map. In the section about
  227 +:ref:`high level analysis tools <sec_high_level>`, the usage of :ref:`cttsmap` was demonstrated. This tool can also
  228 +be used to generate flux maps. For this, the user has to incoorporate a test source in the XML model file. This source
  229 +should consist of a spectral component using a PowerLaw2 model:
  230 +
  231 +.. code-block:: xml
  232 +
  233 + <source name="TestSource" type="PointSource">
  234 + <spectrum type="PowerLaw2">
  235 + <parameter scale="1e-07" name="Integral" min="1e-07" max="1000.0" value="1.0" free="1"/>
  236 + <parameter scale="1.0" name="Index" min="-5.0" max="+5.0" value="-2.0" free="1"/>
  237 + <parameter scale="1.0" name="LowerLimit" min="10.0" max="1000000.0" value="100.0" free="0"/>
  238 + <parameter scale="1.0" name="UpperLimit" min="10.0" max="1000000.0" value="500000.0" free="0"/>
  239 + </spectrum>
  240 + <spatialModel type="SkyDirFunction">
  241 + <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
  242 + <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
  243 + </spatialModel>
  244 + </source>
29 245  
30   -Accessing fit parameters in python
31   -----------------------------------
  246 +The further content of the rest of the XML model file depends on the user requirements:
32 247  
33   -Accessing upper limits in python
34   ---------------------------------
  248 +* For an absolute flux map, all other sky models should be removed
  249 +* For a residual flux map, all other sky models should be kept
35 250  
36   -Manipulate models in python
  251 +The tool :ref:`cttsmap` will create one skymap per free spectral parameter in the model. In the result FITS file,
  252 +there will be an extension called "Integral" that contains the flux map. For the above XML example, the source name
  253 +"TestSource" should be specified to :ref:`cttsmap`.
  254 +
  255 +Speed up TS map computation
37 256 ---------------------------
  257 +Since the computation of a TS map can be extremely time consuming, the option of splitting the computation into several jobs
  258 +is supported. This might be of particular interest if the user has access to a batch farm.
  259 +
  260 +Split TS map computation
  261 +^^^^^^^^^^^^^^^^^^^^^^^^
  262 +For the purpose of job splitting, the hidden parameter ``binmin`` and ``binmax`` were included in the tool. These are integer parameter
  263 +that specify which bins should be computed. For instance, if the map should consist of 30x30 =90 pixels, the use could e.g. execute
  264 +
  265 +.. code-block:: bash
  266 +
  267 + $ cttsmap binmin=0 binmax=299
  268 + ...
  269 + $ cttsmap binmin=300 binmax=599
  270 + ...
  271 + $ cttsmap binmin=600 binmax=899
  272 + ...
  273 +
  274 +on different machines.
  275 +
  276 +.. note::
  277 +
  278 + The output file name of individual jobs can be the same. If the parameters are specified, the output file name will be appended by the suffix
  279 + "_<binmin>_<binmax>" to distinguish between the sliced TS maps.
  280 +
  281 +Merge sliced TS maps
  282 +^^^^^^^^^^^^^^^^^^^^
  283 +The tool :ref:`cstsmapmerge` is intended to take care of merging the files that were produced while splitting the TS map computation.
  284 +There are several options to pass the files to be merged as arguments:
  285 +
  286 +* a space-separated list of file names (e.g. tsmap1.fits tsmap2.fits)
  287 +* a comma-separated list of file names (e.g. tsmap1.fits,tsmap2.fits)
  288 +* a wildcard string (e.g. tsmap*.xml)
  289 +* an ASCII file containing the file names, one file per line (e.g. @mymaps.txt)
  290 +
  291 +In this example the ASCII file method is presented:
  292 +
  293 +.. code-block:: bash
  294 +
  295 + # Put slice files into an ascii file
  296 + $ ls tsmap_*.fits > tsmapfiles.txt
  297 +
  298 + # Run cstsmapmerge
  299 + $ cstsmapmerge
  300 + Input TS map FITS files [@tsmapfiles.txt]
  301 + Output TS map FITS file [mytsmap.fits]
  302 +
38 303  
39 304 Creating a python analysis pipeline
40 305 -----------------------------------
41 306 \ No newline at end of file
... ...