Dmitry Malyshev / ctools

Commit 1536c6d691fc3dd09c438fdb7b31ddee2dcbc596

Authored by Jürgen Knödlseder
1 parent 48427fd1

Update Python documentation.

Inline Side-by-side
Showing 2 changed files with 78 additions and 45 deletions
doc/source/contents.rst
  View file @1536c6d
... ... @@ -10,9 +10,9 @@ Documentation
10 10 :maxdepth: 2
11 11  
12 12 about
  13 + download
13 14 user_manual/index
14 15 reference_manual/index
15 16 faq
16 17 issues
17   - download
18 18 develop
19 19 \ No newline at end of file
... ...
doc/source/user_manual/getting_started/python.rst
  View file @1536c6d
... ... @@ -28,20 +28,20 @@ way as if it were executed from the command line.
28 28  
29 29 .. code-block:: python
30 30  
31   - import ctools
32   - sim = ctools.ctobssim()
33   - sim["inmodel"] = "${CTOOLS}/share/models/crab.xml"
34   - sim["outevents"] = "events.fits"
35   - sim["caldb"] = "prod2"
36   - sim["irf"] = "South_50h"
37   - sim["ra"] = 83.63
38   - sim["dec"] = 22.01
39   - sim["rad"] = 5.0
40   - sim["tmin"] = 0
41   - sim["tmax"] = 1800
42   - sim["emin"] = 0.1
43   - sim["emax"] = 100
44   - sim.execute()
  31 + >>> import ctools
  32 + >>> sim = ctools.ctobssim()
  33 + >>> sim["inmodel"] = "${CTOOLS}/share/models/crab.xml"
  34 + >>> sim["outevents"] = "events.fits"
  35 + >>> sim["caldb"] = "prod2"
  36 + >>> sim["irf"] = "South_0.5h"
  37 + >>> sim["ra"] = 83.63
  38 + >>> sim["dec"] = 22.01
  39 + >>> sim["rad"] = 5.0
  40 + >>> sim["tmin"] = 0.0
  41 + >>> sim["tmax"] = 1800.0
  42 + >>> sim["emin"] = 0.1
  43 + >>> sim["emax"] = 100.0
  44 + >>> sim.execute()
45 45  
46 46 Alternatively, you may "run" the :ref:`ctobssim` tool using
47 47  
... ... @@ -76,7 +76,14 @@ To visualise this container, type:
76 76 Livetime ..................: 1710 s
77 77 Deadtime correction .......: 0.95
78 78 User energy range .........: undefined
79   - ...
  79 + === GCTAPointing ===
  80 + Pointing direction ........: (RA,Dec)=(83.63,22.01)
  81 + === GCTAResponseIrf ===
  82 + Caldb mission .............: cta
  83 + Caldb instrument ..........: prod2
  84 + Response name .............: South_0.5h
  85 + Energy dispersion .........: Not used
  86 + Save energy range .........: undefined
80 87  
81 88 There is one CTA observation in the container and to visualise the events
82 89 in that observation you may type:
... ... @@ -85,7 +92,7 @@ in that observation you may type:
85 92  
86 93 >>> print(sim.obs()[0].events())
87 94 === GCTAEventList ===
88   - Number of events ..........: 31928
  95 + Number of events ..........: 45370
89 96 Time interval .............: 51544.5 - 51544.5 days
90 97 === GEbounds ===
91 98 Number of intervals .......: 1
... ... @@ -118,9 +125,9 @@ terrestrial times (TT) of all events.
118 125 >>> for event in events:
119 126 ... print(event)
120 127 ...
121   - Dir=RA=83.5759, DEC=21.9757 [-0.000597839,-0.000875483] Energy=328.153 GeV Time=-3.15576e+08 s (TT)
122   - Dir=RA=83.5775, DEC=22.0128 [4.88943e-05,-0.000849504] Energy=353.264 GeV Time=-3.15576e+08 s (TT)
123   - Dir=RA=83.4302, DEC=21.931 [-0.00137616,-0.0032343] Energy=117.999 GeV Time=-3.15576e+08 s (TT)
  128 + Dir=RA=83.5743, DEC=21.9746 [-0.000617357,-0.000902194] Energy=328.153 GeV Time=-3.15576e+08 s (TT)
  129 + Dir=RA=83.5756, DEC=22.0127 [4.78931e-05,-0.000879886] Energy=353.264 GeV Time=-3.15576e+08 s (TT)
  130 + Dir=RA=83.6811, DEC=22.0034 [-0.000114612,0.000826825] Energy=550.35 GeV Time=-3.15576e+08 s (TT)
124 131 ...
125 132  
126 133 .. note::
... ... @@ -141,7 +148,7 @@ terrestrial times (TT) of all events.
141 148 >>> sim.run()
142 149 >>> print(sim.obs()[0].events())
143 150 === GCTAEventList ===
144   - Number of events ..........: 15964
  151 + Number of events ..........: 22685
145 152 Time interval .............: 51544.5 - 51544.5 days
146 153 === GEbounds ===
147 154 Number of intervals .......: 1
... ... @@ -180,7 +187,7 @@ To check how the fit went you can inspect the optimiser used by
180 187  
181 188 >>> print(like.opt())
182 189 === GOptimizerLM ===
183   - Optimized function value ..: 101048.382
  190 + Optimized function value ..: 152299.205
184 191 Absolute precision ........: 0.005
185 192 Acceptable value decrease .: 2
186 193 Optimization status .......: converged
... ... @@ -213,8 +220,8 @@ member of the observation container:
213 220 RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
214 221 DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
215 222 Number of spectral par's ..: 3
216   - Prefactor ................: 5.80565e-16 +/- 1.11147e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
217   - Index ....................: -2.49266 +/- 0.0163476 [-0,-5] (free,scale=-1,gradient)
  223 + Prefactor ................: 5.81591e-16 +/- 1.01725e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  224 + Index ....................: -2.50618 +/- 0.0153538 [-0,-5] (free,scale=-1,gradient)
218 225 PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
219 226 Number of temporal par's ..: 1
220 227 Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
... ... @@ -226,8 +233,8 @@ member of the observation container:
226 233 Model type ................: "PowerLaw" * "Constant"
227 234 Number of parameters ......: 4
228 235 Number of spectral par's ..: 3
229   - Prefactor ................: 1.01701 +/- 0.0216563 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
230   - Index ....................: 0.0230864 +/- 0.0123006 [-5,5] (free,scale=1,gradient)
  236 + Prefactor ................: 0.998788 +/- 0.0118567 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
  237 + Index ....................: 0.00990034 +/- 0.00735768 [-5,5] (free,scale=1,gradient)
231 238 PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
232 239 Number of temporal par's ..: 1
233 240 Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
... ... @@ -243,8 +250,8 @@ This is easy from Python:
243 250 >>> like.run()
244 251 >>> print(like.obs().models())
245 252 ...
246   - RA .......................: 83.6341 +/- 0.000958947 [-360,360] deg (free,scale=1)
247   - DEC ......................: 22.0139 +/- 0.000870912 [-90,90] deg (free,scale=1)
  253 + RA .......................: 83.6334 +/- 0.000917669 [-360,360] deg (free,scale=1)
  254 + DEC ......................: 22.0144 +/- 0.000836554 [-90,90] deg (free,scale=1)
248 255  
249 256 The ``like.obs().models()`` method provides the model container, using the
250 257 ``["Crab"]`` operator we access the Crab model in that container and using
... ... @@ -286,7 +293,7 @@ the ``fix()`` method.
286 293 >>> print(like.obs())
287 294 === GObservations ===
288 295 Number of observations ....: 1
289   - Number of predicted events : 15964
  296 + Number of predicted events : 22685
290 297  
291 298 is okay as the ``ctools.ctlike(sim.obs())`` constructor will create
292 299 a copy of the observation container that lives within the :ref:`ctlike`
... ... @@ -312,15 +319,22 @@ ctools provides the Python module ``obsutils`` that may further simplify
312 319 your analysis efforts.
313 320 ``obsutils`` is a Python script that makes use of the GammaLib and ctools
314 321 modules to create standard analysis steps.
  322 +As all Python scripts, ``obsutils`` is part of the cscripts module that
  323 +is imported using
  324 +
  325 +.. code-block:: python
  326 +
  327 + >>> import cscripts
  328 +
315 329 Here an example of how to use ``obsutils``:
316 330  
317 331 .. code-block:: python
318 332  
319 333 >>> import gammalib
320 334 >>> import ctools
321   - >>> from ctools import obsutils
  335 + >>> from cscripts import obsutils
322 336 >>> pattern = obsutils.set_obs_patterns("four", ra=83.63, dec=22.01, offset=1.0)
323   - >>> obs = obsutils.set_obs_list(pattern, duration=1800, emin=0.1, emax=100.0, rad=5.0, caldb="prod2", irf="South_50h")
  337 + >>> obs = obsutils.set_obs_list(pattern, duration=1800, emin=0.1, emax=100.0, rad=5.0, caldb="prod2", irf="South_0.5h")
324 338 >>> print(obs)
325 339 === GObservations ===
326 340 Number of observations ....: 4
... ... @@ -335,7 +349,7 @@ Here an example of how to use ``obsutils``:
335 349 Number of parameters ......: 10
336 350 ...
337 351  
338   -The module is imported using the ``from ctools import obsutils`` directive.
  352 +The module is imported using the ``from cscripts import obsutils`` directive.
339 353 The ``obsutils.set_obs_patterns()`` function will create a pointing
340 354 pattern of four observations located at offset angles of 1 degree from the
341 355 nominal location of the Crab nebula.
... ... @@ -354,10 +368,13 @@ The ``obsutils.sim()`` function then simulates the event data, the
354 368 functions will be added in the future.
355 369  
356 370  
357   -Access common analysis results
358   -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
359   -Here are some examples how to access your analysis results in python. In the following it is
360   -assumed that you have a ctlike object called 'like' which is setup and runs:
  371 +Access analysis results
  372 +~~~~~~~~~~~~~~~~~~~~~~~
  373 +
  374 +Here are some examples that show how to access your analysis results in
  375 +python.
  376 +In the following it is assumed that you have a ctlike object called
  377 +``like`` which is setup and runs:
361 378  
362 379 .. code-block:: python
363 380  
... ... @@ -372,9 +389,9 @@ assumed that you have a ctlike object called 'like' which is setup and runs:
372 389  
373 390 .. code-block:: python
374 391  
375   - obs = like.obs() # Get observations object
376   - obs.models()["Crab"]["Prefactor"].value() # This returns the actual fit value
377   - obs.models()["Crab"]["Prefactor"].error() # This returns the actual fit error
  392 + obs = like.obs() # Get observations object
  393 + obs.models()["Crab"]["Prefactor"].value() # This returns the actual fit value
  394 + obs.models()["Crab"]["Prefactor"].error() # This returns the actual fit error
378 395  
379 396 * Open XML model file from python and print models on screen:
380 397  
... ... @@ -385,16 +402,25 @@ assumed that you have a ctlike object called 'like' which is setup and runs:
385 402  
386 403 * Likelihood value
387 404  
388   - .. code-block:: python
  405 + .. code-block:: python
389 406  
390 407 like.opt().value() # Returns likelihood value
391 408  
392   -* Covariance Matrix
  409 +* Curvature Matrix (aka Hessian)
393 410  
394   - .. code-block:: python
  411 + .. code-block:: python
395 412  
396   - covariance = like.opt().curvature().invert() # Return GMatrix object
397   - print(covariance)
  413 + curvature = like.opt().curvature().invert() # Return GMatrix object
  414 + print(curvature)
  415 +
  416 + To get the covariance matrix, the curvature matrix needs to be
  417 + inverted:
  418 +
  419 + .. code-block:: python
  420 +
  421 + covariance = curvature.invert()
  422 + print(covariance)
  423 +
398 424  
399 425 Modify and work on XML models
400 426 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
... ... @@ -415,7 +441,8 @@ Modify and work on XML models
415 441  
416 442 .. code-block:: python
417 443  
418   - models["Crab"]["Prefactor"].is_free()
  444 + if models["Crab"]["Prefactor"].is_free():
  445 + ...
419 446  
420 447 * Set parameter value
421 448  
... ... @@ -430,6 +457,12 @@ Modify and work on XML models
430 457 models["Crab"]["Index"].min(-5.0)
431 458 models["Crab"]["Index"].max(-1.0)
432 459  
  460 + or quicker
  461 +
  462 + .. code-block:: python
  463 +
  464 + models["Crab"]["Index"].range(-5.0,-1.0)
  465 +
433 466 * Loop over models and parameters (e.g. fix all parameters)
434 467  
435 468 .. code-block:: python
... ...