Karim Sabri / ctools

Commit e0f8913415367049fbebbb8113377267990aba7c

Authored by Luigi Tibaldo
1 parent d0da27ce

minimally working for binned and unbinned analysis

Inline Side-by-side
Showing 3 changed files with 610 additions and 246 deletions
cscripts/Makefile.am
  View file @e0f8913
... ... @@ -43,78 +43,42 @@
43 43 # Scripts to be distributed in $(pyexecdir)/cscripts
44 44  
45 45 pythonscripts = $(top_builddir)/cscripts/__init__.py \
46   -
47 46 $(srcdir)/csbkgmodel.py \
48   -
49 47 $(srcdir)/cscaldb.py \
50   -
51 48 $(srcdir)/csebins.py \
52   -
53 49 $(srcdir)/csinfo.py \
54   -
55 50 $(srcdir)/cslightcrv.py \
56   -
57 51 $(srcdir)/csmodelinfo.py \
58   -
59 52 $(srcdir)/csmodelmerge.py \
60   -
61 53 $(srcdir)/csmodelselect.py \
62   -
63 54 $(srcdir)/csmodelsois.py \
64   -
65 55 $(srcdir)/csobs2caldb.py \
66   -
67 56 $(srcdir)/csobsdef.py \
68   -
69 57 $(srcdir)/csobsinfo.py \
70   -
71 58 $(srcdir)/csobsselect.py \
72   -
73 59 $(srcdir)/csphagen.py \
74   -
75 60 $(srcdir)/csphasecrv.py \
76   -
77 61 $(srcdir)/cspull.py \
78   -
79 62 $(srcdir)/csresmap.py \
80   -
81 63 $(srcdir)/csresspec.py \
82   -
83 64 $(srcdir)/csroot2caldb.py \
84   -
85 65 $(srcdir)/cssens.py \
86   -
87 66 $(srcdir)/csspec.py \
88   -
89 67 $(srcdir)/cssrcdetect.py \
90   -
91 68 $(srcdir)/cstsdist.py \
92   -
93 69 $(srcdir)/cstsmapmerge.py \
94   -
95 70 $(srcdir)/cstsmapsplit.py \
96   -
97 71 $(srcdir)/csviscube.py \
98   -
99 72 $(srcdir)/csscs.py \
100 73 $(srcdir)/csworkflow.py \
101   -
102 74 $(srcdir)/csfindobs.py \
103   -
104 75 $(srcdir)/csiactcopy.py \
105   -
106 76 $(srcdir)/csiactdata.py \
107   -
108 77 $(srcdir)/csiactobs.py \
109   -
110 78 $(srcdir)/obsutils.py \
111   -
112 79 $(srcdir)/modutils.py \
113   -
114 80 $(srcdir)/calutils.py \
115   -
116 81 $(srcdir)/ioutils.py \
117   -
118 82 $(srcdir)/mputils.py
119 83  
120 84  
... ... @@ -122,64 +86,35 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
122 86 # Parfiles to be distributed in $(prefix)/syspfiles
123 87  
124 88 parfiles = $(srcdir)/csbkgmodel.par \
125   -
126 89 $(srcdir)/cscaldb.par \
127   -
128 90 $(srcdir)/csebins.par \
129   -
130 91 $(srcdir)/cslightcrv.par \
131   -
132 92 $(srcdir)/csmodelinfo.par \
133   -
134 93 $(srcdir)/csmodelmerge.par \
135   -
136 94 $(srcdir)/csmodelselect.par \
137   -
138 95 $(srcdir)/csmodelsois.par \
139   -
140 96 $(srcdir)/csobs2caldb.par \
141   -
142 97 $(srcdir)/csobsdef.par \
143   -
144 98 $(srcdir)/csobsinfo.par \
145   -
146 99 $(srcdir)/csobsselect.par \
147   -
148 100 $(srcdir)/csphagen.par \
149   -
150 101 $(srcdir)/csphasecrv.par \
151   -
152 102 $(srcdir)/cspull.par \
153   -
154 103 $(srcdir)/csresmap.par \
155   -
156 104 $(srcdir)/csresspec.par \
157   -
158 105 $(srcdir)/csroot2caldb.par \
159   -
160 106 $(srcdir)/cssens.par \
161   -
162 107 $(srcdir)/csspec.par \
163   -
164 108 $(srcdir)/cssrcdetect.par \
165   -
166 109 $(srcdir)/cstsdist.par \
167   -
168 110 $(srcdir)/cstsmapmerge.par \
169   -
170 111 $(srcdir)/cstsmapsplit.par \
171   -
172 112 $(srcdir)/csviscube.par \
173   -
174 113 $(srcdir)/csscs.par \
175 114 $(srcdir)/csworkflow.par \
176   -
177 115 $(srcdir)/csfindobs.par \
178   -
179 116 $(srcdir)/csiactcopy.par \
180   -
181 117 $(srcdir)/csiactdata.par \
182   -
183 118 $(srcdir)/csiactobs.par
184 119  
185 120  
... ... @@ -187,66 +122,36 @@ parfiles = $(srcdir)/csbkgmodel.par \
187 122 # Scripts for which symbolic links will be created in bindir
188 123  
189 124 linkscripts = csbkgmodel \
190   -
191 125 cscaldb \
192   -
193 126 csebins \
194   -
195 127 csinfo \
196   -
197 128 cslightcrv \
198   -
199 129 csmodelinfo \
200   -
201 130 csmodelmerge \
202   -
203 131 csmodelselect \
204   -
205 132 csmodelsois \
206   -
207 133 csobs2caldb \
208   -
209 134 csobsdef \
210   -
211 135 csobsinfo \
212   -
213 136 csobsselect \
214   -
215 137 csphagen \
216   -
217 138 csphasecrv \
218   -
219 139 cspull \
220   -
221 140 csresmap \
222   -
223 141 csresspec \
224   -
225 142 csroot2caldb \
226   -
227 143 cssens \
228   -
229 144 csspec \
230   -
231 145 cssrcdetect \
232   -
233 146 cstsdist \
234   -
235 147 cstsmapmerge \
236   -
237 148 cstsmapsplit \
238   -
239 149 csviscube \
240   -
241 150 csscs \
242 151 csworkflow \
243   -
244 152 csfindobs \
245   -
246 153 csiactcopy \
247   -
248 154 csiactdata \
249   -
250 155 csiactobs
251 156  
252 157  
... ... @@ -254,78 +159,42 @@ linkscripts = csbkgmodel \
254 159 # Test scripts to be distributed in $(pyexecdir)/cscripts/tests
255 160  
256 161 testscripts = $(srcdir)/tests/__init__.py \
257   -
258 162 $(top_srcdir)/test/testing.py \
259   -
260 163 $(top_srcdir)/test/test_python_cscripts.py \
261   -
262 164 $(top_srcdir)/test/test_cscript.py \
263   -
264 165 $(top_srcdir)/test/test_cscaldb.py \
265   -
266 166 $(top_srcdir)/test/test_csbkgmodel.py \
267   -
268 167 $(top_srcdir)/test/test_csebins.py \
269   -
270 168 $(top_srcdir)/test/test_csfindobs.py \
271   -
272 169 $(top_srcdir)/test/test_csinfo.py \
273   -
274 170 $(top_srcdir)/test/test_cslightcrv.py \
275   -
276 171 $(top_srcdir)/test/test_csmodelinfo.py \
277   -
278 172 $(top_srcdir)/test/test_csmodelmerge.py \
279   -
280 173 $(top_srcdir)/test/test_csmodelselect.py \
281   -
282 174 $(top_srcdir)/test/test_csmodelsois.py \
283   -
284 175 $(top_srcdir)/test/test_csobs2caldb.py \
285   -
286 176 $(top_srcdir)/test/test_csobsdef.py \
287   -
288 177 $(top_srcdir)/test/test_csobsinfo.py \
289   -
290 178 $(top_srcdir)/test/test_csobsselect.py \
291   -
292 179 $(top_srcdir)/test/test_csphagen.py \
293   -
294 180 $(top_srcdir)/test/test_csphasecrv.py \
295   -
296 181 $(top_srcdir)/test/test_cspull.py \
297   -
298 182 $(top_srcdir)/test/test_csresmap.py \
299   -
300 183 $(top_srcdir)/test/test_csresspec.py \
301   -
302 184 $(top_srcdir)/test/test_csroot2caldb.py \
303   -
304 185 $(top_srcdir)/test/test_cssens.py \
305   -
306 186 $(top_srcdir)/test/test_csspec.py \
307   -
308 187 $(top_srcdir)/test/test_cssrcdetect.py \
309   -
310 188 $(top_srcdir)/test/test_cstsdist.py \
311   -
312 189 $(top_srcdir)/test/test_cstsmapmerge.py \
313   -
314 190 $(top_srcdir)/test/test_cstsmapsplit.py \
315   -
316 191 $(top_srcdir)/test/test_csviscube.py \
317   -
318 192 $(top_srcdir)/test/test_csscs.py \
319 193 $(top_srcdir)/test/test_csworkflow.py \
320   -
321 194 $(top_srcdir)/test/test_obsutils.py \
322   -
323 195 $(top_srcdir)/test/test_mputils.py \
324   -
325 196 $(top_srcdir)/test/test_csiactcopy.py \
326   -
327 197 $(top_srcdir)/test/test_csiactdata.py \
328   -
329 198 $(top_srcdir)/test/test_csiactobs.py
330 199  
331 200  
... ... @@ -333,153 +202,79 @@ testscripts = $(srcdir)/tests/__init__.py \
333 202 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data
334 203  
335 204 testdata = $(top_srcdir)/test/data/crab.xml \
336   -
337 205 $(top_srcdir)/test/data/obs_unbinned.xml \
338   -
339 206 $(top_srcdir)/test/data/obs_unbinned_two.xml \
340   -
341 207 $(top_srcdir)/test/data/obs_stacked.xml \
342   -
343 208 $(top_srcdir)/test/data/obs_stacked_two.xml \
344   -
345 209 $(top_srcdir)/test/data/crab_events.fits \
346   -
347 210 $(top_srcdir)/test/data/crab_cntmap.fits \
348   -
349 211 $(top_srcdir)/test/data/crab_expcube.fits.gz \
350   -
351 212 $(top_srcdir)/test/data/crab_psfcube.fits.gz \
352   -
353 213 $(top_srcdir)/test/data/crab_edispcube.fits.gz \
354   -
355 214 $(top_srcdir)/test/data/crab_bkgcube.fits.gz \
356   -
357 215 $(top_srcdir)/test/data/crab_onoff.xml \
358   -
359 216 $(top_srcdir)/test/data/crab_offaxis1.fits.gz \
360   -
361 217 $(top_srcdir)/test/data/crab_offaxis2.fits.gz \
362   -
363 218 $(top_srcdir)/test/data/crab_exclusion.fits.gz \
364   -
365 219 $(top_srcdir)/test/data/crab_src_reg.reg \
366   -
367 220 $(top_srcdir)/test/data/crab_bkg_reg.reg \
368   -
369 221 $(top_srcdir)/test/data/tsmap_0.fits \
370   -
371 222 $(top_srcdir)/test/data/tsmap_1.fits \
372   -
373 223 $(top_srcdir)/test/data/exclusion.reg \
374   -
375 224 $(top_srcdir)/test/data/pntdef_minimal.dat \
376   -
377 225 $(top_srcdir)/test/data/pntdef_maximal.dat \
378   -
379 226 $(top_srcdir)/test/data/models.txt \
380   -
381 227 $(top_srcdir)/test/data/tsmaps.txt \
382   -
383 228 $(top_srcdir)/test/data/model_cube_background1.xml \
384   -
385 229 $(top_srcdir)/test/data/model_cube_background2.xml \
386   -
387 230 $(top_srcdir)/test/data/crab_bkgcube.xml \
388   -
389 231 $(top_srcdir)/test/data/irf.root \
390   -
391 232 $(top_srcdir)/test/data/irf_file.fits.gz \
392   -
393 233 $(top_srcdir)/test/data/workflow.xml \
394   -
395 234 $(top_srcdir)/test/data/butterfly.txt \
396   -
397 235 $(top_srcdir)/test/data/pha.fits \
398   -
399 236 $(top_srcdir)/test/data/pha_channels.fits \
400   -
401 237 $(top_srcdir)/test/data/pull.dat \
402   -
403 238 $(top_srcdir)/test/data/sensitivity.dat \
404   -
405 239 $(top_srcdir)/test/data/sensitivity_integral.dat \
406   -
407 240 $(top_srcdir)/test/data/skymap.fits \
408   -
409 241 $(top_srcdir)/test/data/spectrum.fits \
410   -
411 242 $(top_srcdir)/test/data/residual.fits \
412   -
413 243 $(top_srcdir)/test/data/residual_onoff.fits \
414   -
415 244 $(top_srcdir)/test/data/lightcurve.fits \
416   -
417 245 $(top_srcdir)/test/data/ts.dat \
418   -
419 246 $(top_srcdir)/test/data/models_spatial.xml \
420   -
421 247 $(top_srcdir)/test/data/models_spectral.xml \
422   -
423 248 $(top_srcdir)/test/data/model_filefunction.txt \
424   -
425 249 $(top_srcdir)/test/data/model_map.fits.gz \
426   -
427 250 $(top_srcdir)/test/data/model_mapcube.fits.gz \
428   -
429 251 $(top_srcdir)/test/data/iact_runlist.dat \
430   -
431 252 $(top_srcdir)/test/data/iact_runlist_1.dat \
432   -
433 253 $(top_srcdir)/test/data/iact_runlist_2.dat \
434   -
435 254 $(top_srcdir)/test/data/model_temporal_phasecurve.xml \
436   -
437 255 $(top_srcdir)/test/data/model_temporal_phasecurve.fits \
438   -
439 256 $(top_srcdir)/test/data/phased_events.fits.gz \
440   -
441 257 $(top_srcdir)/test/data/onoff_arf.fits.gz \
442   -
443 258 $(top_srcdir)/test/data/onoff_arf1.fits.gz \
444   -
445 259 $(top_srcdir)/test/data/onoff_arf2.fits.gz \
446   -
447 260 $(top_srcdir)/test/data/onoff_model.xml \
448   -
449 261 $(top_srcdir)/test/data/onoff_obs.xml \
450   -
451 262 $(top_srcdir)/test/data/onoff_obs1.xml \
452   -
453 263 $(top_srcdir)/test/data/onoff_obs2.xml \
454   -
455 264 $(top_srcdir)/test/data/onoff_obs_stacked.xml \
456   -
457 265 $(top_srcdir)/test/data/onoff_pha_off.fits.gz \
458   -
459 266 $(top_srcdir)/test/data/onoff_pha_off1.fits.gz \
460   -
461 267 $(top_srcdir)/test/data/onoff_pha_off2.fits.gz \
462   -
463 268 $(top_srcdir)/test/data/onoff_pha_on.fits.gz \
464   -
465 269 $(top_srcdir)/test/data/onoff_pha_on1.fits.gz \
466   -
467 270 $(top_srcdir)/test/data/onoff_pha_on2.fits.gz \
468   -
469 271 $(top_srcdir)/test/data/onoff_rmf.fits.gz \
470   -
471 272 $(top_srcdir)/test/data/onoff_rmf1.fits.gz \
472   -
473 273 $(top_srcdir)/test/data/onoff_rmf2.fits.gz \
474   -
475 274 $(top_srcdir)/test/data/model_crab_radialacceptance.xml \
476   -
477 275 $(top_srcdir)/test/data/sigdist_map.fits.gz \
478   -
479 276 $(top_srcdir)/test/data/variability.fits.gz \
480   -
481 277 $(top_srcdir)/test/data/viscube.fits.gz \
482   -
483 278 $(top_srcdir)/test/data/hess_bkg_lookup.fits
484 279  
485 280  
... ... @@ -487,9 +282,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \
487 282 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata
488 283  
489 284 testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \
490   -
491 285 $(top_srcdir)/test/data/iactdata/master.json \
492   -
493 286 $(top_srcdir)/test/data/iactdata/obs-index.fits.gz
494 287  
495 288  
... ... @@ -497,23 +290,14 @@ testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \
497 290 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata/events
498 291  
499 292 testiactevents = $(top_srcdir)/test/data/iactdata/events/irf_file.fits.gz \
500   -
501 293 $(top_srcdir)/test/data/iactdata/events/events_0.fits.gz \
502   -
503 294 $(top_srcdir)/test/data/iactdata/events/events_1.fits.gz \
504   -
505 295 $(top_srcdir)/test/data/iactdata/events/events_2.fits.gz \
506   -
507 296 $(top_srcdir)/test/data/iactdata/events/events_3.fits.gz \
508   -
509 297 $(top_srcdir)/test/data/iactdata/events/events_4.fits.gz \
510   -
511 298 $(top_srcdir)/test/data/iactdata/events/events_5.fits.gz \
512   -
513 299 $(top_srcdir)/test/data/iactdata/events/events_6.fits.gz \
514   -
515 300 $(top_srcdir)/test/data/iactdata/events/events_7.fits.gz \
516   -
517 301 $(top_srcdir)/test/data/iactdata/events/events_8.fits.gz
518 302  
519 303  
... ... @@ -557,9 +341,7 @@ testiactevents_DATA = $(testiactevents)
557 341 # Files to be distributed in the tar file
558 342  
559 343 EXTRA_DIST = $(binscripts) $(pythonscripts) $(parfiles) \
560   -
561 344 $(testscripts) $(testdata) $(testiactdata) \
562   -
563 345 $(testiactevents)
564 346  
565 347  
... ... @@ -569,23 +351,15 @@ EXTRA_DIST = $(binscripts) $(pythonscripts) $(parfiles) \
569 351 install-data-hook:
570 352  
571 353 @binlist='$(binscripts)';\
572   -
573 354 for p in $$binlist; do \
574   -
575 355 chmod +x "$(DESTDIR)$(bindir)/$$p"; \
576   -
577 356 done
578 357  
579 358 @linklist='$(linkscripts)';\
580   -
581 359 for p in $$linklist; do \
582   -
583 360 rm -rf "$(DESTDIR)$(bindir)/$$p"; \
584   -
585 361 cp "$(DESTDIR)$(pyexecdir)/cscripts/$$p.py" "$(DESTDIR)$(bindir)/$$p"; \
586   -
587 362 chmod +x "$(DESTDIR)$(bindir)/$$p"; \
588   -
589 363 done
590 364  
591 365  
... ... @@ -595,11 +369,8 @@ install-data-hook:
595 369 uninstall-hook:
596 370  
597 371 @linklist='$(linkscripts)';\
598   -
599 372 for p in $$linklist; do \
600   -
601 373 rm -rf "$(DESTDIR)$(bindir)/$$p"; \
602   -
603 374 done
604 375  
605 376  
... ... @@ -609,21 +380,13 @@ uninstall-hook:
609 380 clean-local:
610 381  
611 382 @if ! [ '@top_srcdir@' -ef '@top_builddir@' ]; then \
612   -
613 383 rm -rf $(top_builddir)/cscripts/cs*.py; \
614   -
615 384 rm -rf $(top_builddir)/cscripts/obsutils.py; \
616   -
617 385 rm -rf $(top_builddir)/cscripts/modutils.py; \
618   -
619 386 rm -rf $(top_builddir)/cscripts/calutils.py; \
620   -
621 387 rm -rf $(top_builddir)/cscripts/ioutils.py; \
622   -
623 388 rm -rf $(top_builddir)/cscripts/mputils.py; \
624   -
625 389 rm -rf $(top_builddir)/cscripts/*.par; \
626   -
627 390 fi
628 391  
629 392 rm -rf $(top_builddir)/cscripts/*.pyc
... ...
cscripts/csscs.par
  View file @e0f8913
... ... @@ -14,13 +14,56 @@
14 14 # #
15 15 #############################################################################
16 16 #
17   -# General parameters
18   -#===================
19   -# Your parameters go here
  17 +# Input and output parameters
  18 +#============================
  19 +inobs, f, a, events.fits,,, "Input event list, counts cube, or observation definition XML file"
  20 +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
  21 +srcnames, s, a, Crab,,, "Semicolon separated list of sources"
  22 +expcube, f, a, NONE,,, "Input exposure cube file"
  23 +psfcube, f, a, NONE,,, "Input PSF cube file"
  24 +edispcube, f, a, NONE,,, "Input energy dispersion cube file"
  25 +bkgcube, f, a, NONE,,, "Input background cube file"
  26 +caldb, s, a, prod2,,, "Calibration database"
  27 +irf, s, a, South_0.5h,,, "Instrument response function"
  28 +inexclusion, f, h, NONE,,, "Exclusion regions"
  29 +edisp, b, h, no,,, "Apply energy dispersion?"
  30 +outfile, f, a, map.fits,,, "Output map"
  31 +
  32 +#
  33 +# Script parameters
  34 +#==================
  35 +coordsys, s, a, CEL,CEL|GAL,, "Coordinate System"
  36 +proj, s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
  37 +xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
  38 +yref, r, a, 22.51,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
  39 +nxpix, i, a, 20,,, "Size of the X axis in pixels"
  40 +nypix, i, a, 20,,, "Size of the Y axis in pixels"
  41 +binsz, r, a, 0.1,,, "Pixel size (deg/pixel)"
  42 +roisz, r, a, 2.,1.42,1000, "Size of the region of interest (numbers of binsz)"
  43 +statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
  44 +method, s, a, 3D, 3D|ONOFF,, "Analysis method"
  45 +ebinalg, s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins"
  46 +emin, r, a, 0.1,,, "Lower energy limit (TeV)"
  47 +emax, r, a, 100.0,,, "Upper energy limit (TeV)"
  48 +enumbins, i, a, 20,,, "Number of energy bins"
  49 +bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation"
  50 +srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region"
  51 +bkgregmin, i, h, 2,1,1000, "Minimum number of background regions"
  52 +use_model_bkg, b, h, yes,,, "Use background model?"
  53 +maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation"
  54 +etruemin, r, h, 0.01,,, "Minimum true energy (TeV)"
  55 +etruemax, r, h, 300.0,,, "Maximum true energy (TeV)"
  56 +etruebins, i, h, 30,,, "Number of bins per decade for true energy bins"
  57 +calc_ts, b, h, yes,,, "Compute TS value in each bin"
  58 +calc_ulim, b, h, yes,,, "Compute upper limit in each bin"
  59 +fix_srcs, b, h, yes,,, "Fix other skymodel parameters"
  60 +fix_bkg, b, h, no,,, "Fix background parameters"
20 61  
21 62 #
22 63 # Standard parameters
23 64 #====================
  65 +nthreads, i, h, 0,0,1000, "Number of parallel processes (0=use all available CPUs)"
  66 +publish, b, h, no,,, "Publish output map on VO Hub?"
24 67 chatter, i, h, 2,0,4, "Chattiness of output"
25 68 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
26 69 debug, b, h, no,,, "Activate debugging mode?"
... ...
cscripts/csscs.py
  View file @e0f8913
... ... @@ -21,6 +21,9 @@
21 21 import sys
22 22 import gammalib
23 23 import ctools
  24 +from cscripts import mputils
  25 +from cscripts import obsutils
  26 +import math
24 27  
25 28  
26 29 # ================= #
... ... @@ -45,6 +48,13 @@ class csscs(ctools.csobservation):
45 48 # Initialise application by calling the base class constructor
46 49 self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
47 50  
  51 + # Initialise data members
  52 + self._nthreads = 0
  53 + self._roisz = 0.
  54 + self._srcnames = []
  55 + self._fits = None
  56 + self._excl_reg_map = None # Exclusion region map for on/off analysis
  57 +
48 58 # Return
49 59 return
50 60  
... ... @@ -59,7 +69,13 @@ class csscs(ctools.csobservation):
59 69 Pickled instance
60 70 """
61 71 # Set pickled dictionary
62   - state = {'base' : ctools.csobservation.__getstate__(self)}
  72 + # Set pickled dictionary
  73 + state = {'base' : ctools.csobservation.__getstate__(self),
  74 + 'nthreads' : self._nthreads,
  75 + 'roisz' : self._roisz,
  76 + 'srcnames' : self._srcnames,
  77 + 'fits' : self._fits,
  78 + 'excl_reg_map' : self._excl_reg_map}
63 79  
64 80 # Return pickled dictionary
65 81 return state
... ... @@ -75,6 +91,11 @@ class csscs(ctools.csobservation):
75 91 """
76 92 # Set state
77 93 ctools.csobservation.__setstate__(self, state['base'])
  94 + self._nthreads = state['nthreads']
  95 + self._roisz = state['roisz']
  96 + self._srcnames = state['srcnames']
  97 + self._fits = state['fits']
  98 + self._excl_reg_map = state['excl_reg_map']
78 99  
79 100 # Return
80 101 return
... ... @@ -85,14 +106,497 @@ class csscs(ctools.csobservation):
85 106 """
86 107 Get parameters from parfile
87 108 """
88   - # TODO: Add code to query all relevant parameters
  109 +
  110 + # Set observation if not done before
  111 + if self.obs().is_empty():
  112 + self.obs(self._get_observations())
  113 +
  114 + # Collect number of unbinned, binned and On/Off observations in
  115 + # observation container
  116 + n_unbinned = 0
  117 + n_binned = 0
  118 + n_onoff = 0
  119 + for obs in self.obs():
  120 + if obs.classname() == 'GCTAObservation':
  121 + if obs.eventtype() == 'CountsCube':
  122 + n_binned += 1
  123 + else:
  124 + n_unbinned += 1
  125 + elif obs.classname() == 'GCTAOnOffObservation':
  126 + n_onoff += 1
  127 +
  128 + # Check that we have at least one unbinned or binned CTA observation
  129 + if n_unbinned == 0 and n_binned == 0:
  130 + msg = 'No unbinned or binned CTA observations found. ' \
  131 + 'Please provide at least one unbinned or binned ' \
  132 + 'CTA observation.'
  133 + raise RuntimeError(msg)
  134 +
  135 + # Set observation statistic
  136 + self._set_obs_statistic(gammalib.toupper(self['statistic'].string()))
  137 +
  138 + # Set models if there are none in the container
  139 + if self.obs().models().size() == 0:
  140 + self.obs().models(self['inmodel'].filename())
  141 +
  142 + # Query exclusion region
  143 + self["inexclusion"].query()
  144 +
  145 + # Query the map definition parameters
  146 + self['xref'].real()
  147 + self['yref'].real()
  148 + self['coordsys'].string()
  149 + self['proj'].string()
  150 + self['nxpix'].integer()
  151 + self['nypix'].integer()
  152 +
  153 + # Compute the size of the ROIs for the analysis in deg
  154 + self._roisz = self['binsz'].real() * self['roisz'].real()
  155 +
  156 + # Query energy boundaries
  157 + self["emin"].real()
  158 + self["emax"].real()
  159 +
  160 + # If we have unbinned observations query the analysis method
  161 + if n_unbinned > 0:
  162 + # If method is ONOFF query some csphagen parameters
  163 + if self['method'].string() == 'ONOFF':
  164 + pass
  165 + # TODO: implement Onoff analysis
  166 + # self["enumbins"].integer()
  167 + # self['use_model_bkg'].boolean()
  168 + # if self["bkgmethod"].string() == "REFLECTED":
  169 + # self["srcshape"].string()
  170 + # self["bkgregmin"].integer()
  171 + # self["maxoffset"].real()
  172 + # self["etruemin"].real()
  173 + # self["etruemax"].real()
  174 + # self["etruebins"].integer()
  175 +
  176 + # Query target source names
  177 + srcnames = self['srcnames'].string()
  178 +
  179 + # Fashion target sources into Python list
  180 + self._srcnames = srcnames.split(';')
  181 + # Strip leading and trailing spaces from names
  182 + for s, name in enumerate(self._srcnames):
  183 + self._srcnames[s] = name.strip()
  184 +
  185 + # Verify that the target sources are in the model
  186 + src_in_model = True
  187 + for name in self._srcnames:
  188 + src_in_model = src_in_model and self.obs().models().contains(name)
  189 + if src_in_model == False:
  190 + msg = 'Not all target sources are present in input model.'
  191 + raise RuntimeError(msg)
  192 +
  193 + # Query the hidden parameters, just in case
  194 + self['edisp'].boolean()
  195 + self['calc_ulim'].boolean()
  196 + self['calc_ts'].boolean()
  197 + self['fix_bkg'].boolean()
  198 + self['fix_srcs'].boolean()
  199 +
  200 + # Read ahead output parameters
  201 + if self._read_ahead():
  202 + self['outfile'].query()
89 203  
90 204 # Write input parameters into logger
91 205 self._log_parameters(gammalib.TERSE)
92 206  
  207 + # Set number of processes for multiprocessing
  208 + self._nthreads = mputils.nthreads(self)
  209 +
  210 + # Return
  211 + return
  212 +
  213 + def _create_map(self):
  214 + """
  215 + Create a gammalib.GSkyMap object to store the output
  216 +
  217 + Returns
  218 + -------
  219 + skymap : `~gammalib.GSkyMap`
  220 + """
  221 +
  222 + skymap = gammalib.GSkyMap(self['proj'].string(), self['coordsys'].string(),
  223 + self['xref'].real(), self['yref'].real(),
  224 + -self['binsz'].real(), self['binsz'].real(),
  225 + self['nxpix'].integer(), self['nypix'].integer(),
  226 + 1)
  227 +
  228 + return skymap
  229 +
  230 + def _adjust_model_pars(self):
  231 + """
  232 + Adjust model parameters
  233 + """
  234 + # Write header
  235 + self._log_header1(gammalib.TERSE, 'Adjust model parameters')
  236 +
  237 + # Adjust model parameters dependent on input user parameters
  238 + for model in self.obs().models():
  239 +
  240 + # Initialise TS flag for all models to false
  241 + model.tscalc(False)
  242 +
  243 + # Log model name
  244 + self._log_header3(gammalib.EXPLICIT, model.name())
  245 +
  246 + # Freeze all parameters except the normalization
  247 + # which is assumed to be the first spectral parameter
  248 + normpar = model.spectral()[0]
  249 + for par in model:
  250 + if par.name() == normpar.name():
  251 + pass
  252 + else:
  253 + if par.is_free():
  254 + self._log_string(gammalib.EXPLICIT, ' Fixing "' + par.name() + '"')
  255 + par.fix()
  256 +
  257 + # Deal with the target sources
  258 + if model.name() in self._srcnames:
  259 +
  260 + # Free the normalisation parameter which is assumed to be
  261 + # the first spectral parameter
  262 + if normpar.is_fixed():
  263 + self._log_string(gammalib.EXPLICIT, ' Freeing "'+normpar.name()+'"')
  264 + normpar.free()
  265 +
  266 + # Optionally compute Test Statistic value
  267 + if self['calc_ts'].boolean():
  268 + model.tscalc(True)
  269 +
  270 + # Deal with background models
  271 + elif self['fix_bkg'].boolean() and not model.classname() == 'GModelSky':
  272 +
  273 + if normpar.is_free():
  274 + self._log_string(gammalib.EXPLICIT, ' Fixing "'+normpar.name()+'"')
  275 + normpar.fix()
  276 +
  277 +
  278 + # Deal with background sources
  279 + elif self['fix_srcs'].boolean() and model.classname() == 'GModelSky':
  280 +
  281 + if normpar.is_free():
  282 + self._log_string(gammalib.EXPLICIT, ' Fixing "'+normpar.name()+'"')
  283 + normpar.fix()
  284 +
93 285 # Return
94 286 return
95 287  
  288 + def _mask_cube(self,obs,ra,dec,rad):
  289 + """
  290 + Mask cube observation
  291 +
  292 + Parameters
  293 + ----------
  294 + obs : `gammalib.GCTAObservation` Input observation of type cube
  295 + ra : `float` R.A. (deg)
  296 + dec : `float` Dec (deg)
  297 + rad : `float` radius (deg)
  298 +
  299 + Returns
  300 + -------
  301 + new_obs : `gammalib.GCTAObservation`
  302 + Output observation with masked cube
  303 + """
  304 +
  305 + # Put observation in container
  306 + container = gammalib.GObservations()
  307 + container.append(obs)
  308 +
  309 + # Filter cube according to ROI and user energy range
  310 + cubemask = ctools.ctcubemask(container)
  311 + cubemask['ra'] = ra
  312 + cubemask['dec'] = dec
  313 + cubemask['rad'] = rad
  314 + cubemask['regfile'] = 'NONE'
  315 + cubemask['emin'] = self['emin'].real()
  316 + cubemask['emax'] = self['emax'].real()
  317 + cubemask.run()
  318 +
  319 + # Get deep copy of filtered observation
  320 + new_obs = cubemask.obs()[0].copy()
  321 +
  322 + # Return
  323 + return new_obs
  324 +
  325 + def _mask_evlist(self,obs,ra,dec,rad):
  326 + """
  327 + Mask event list
  328 +
  329 + Parameters
  330 + ----------
  331 + obs : `gammalib.GCTAObservation` Input observation of type cube
  332 + ra : `float` R.A. (deg)
  333 + dec : `float` Dec (deg)
  334 + rad : `float` radius (deg)
  335 +
  336 + Returns
  337 + -------
  338 + new_obs : `gammalib.GCTAObservation`
  339 + Output observation with masked event list
  340 + """
  341 +
  342 + # Put observation in container
  343 + container = gammalib.GObservations()
  344 + container.append(obs)
  345 +
  346 + # Filter event list according to ROI and user energy range
  347 + select = ctools.ctselect(container)
  348 + select['ra'] = ra
  349 + select['dec'] = dec
  350 + select['rad'] = rad
  351 + select['emin'] = self['emin'].real()
  352 + select['emax'] = self['emax'].real()
  353 + select['tmin'] = 'INDEF'
  354 + select['tmax'] = 'INDEF'
  355 + select.run()
  356 +
  357 + # Get deep copy of filtered observation
  358 + new_obs = select.obs()[0].copy()
  359 +
  360 + # Return
  361 + return new_obs
  362 +
  363 + def _mask_onoff(self,obs,ra,dec,rad):
  364 + """
  365 + Create On/Off observation
  366 + with On region matching the required mask
  367 +
  368 + Parameters
  369 + ----------
  370 + obs : `gammalib.GCTAObservation` Input observation of type cube
  371 + ra : `float` R.A. (deg)
  372 + dec : `float` Dec (deg)
  373 + rad : `float` radius (deg)
  374 +
  375 + Returns
  376 + -------
  377 + new_obs : `gammalib.GCTAOnoffObservation`
  378 + Output On/Off observation
  379 + """
  380 +
  381 + # Put observation in container
  382 + container = gammalib.GObservations()
  383 + container.append(obs)
  384 +
  385 + onoff_obs = obsutils.get_onoff_obs(self, container, ra=ra, dec=dec, rad=rad)
  386 +
  387 + return onoff_obs[0]
  388 +
  389 + def _mask_observations(self,ra,dec,rad):
  390 + """
  391 + Create observations restricted to circular ROI
  392 +
  393 + Parameters
  394 + ----------
  395 + ra : `float` R.A. (deg)
  396 + dec : `float` Dec (deg)
  397 + rad : `float` radius (deg)
  398 +
  399 + Returns
  400 + -------
  401 + new_obs : `~gammalib.GObservations`
  402 + Observations in circular ROI
  403 + """
  404 +
  405 + # Create output observation container
  406 + new_obs = gammalib.GObservations()
  407 +
  408 + # Loop over input observations
  409 + for obs in self.obs():
  410 + if obs.classname() == 'GCTAObservation':
  411 + if obs.eventtype() == 'CountsCube':
  412 + masked_obs = self._mask_cube(obs,ra,dec,rad)
  413 + else:
  414 + if self['method'].string() == '3D':
  415 + masked_obs = self._mask_evlist(obs,ra,dec,rad)
  416 + elif self['method'].string() == 'ONOFF':
  417 + pass
  418 + # TODO: implement Onoff analysis
  419 + # masked_obs = self._mask_onoff(obs, ra, dec, rad)
  420 + new_obs.append(masked_obs)
  421 + # Skip On-Off and non-CTA observations
  422 + else:
  423 + pass
  424 +
  425 + # Append models to output observations
  426 + new_obs.models(self.obs().models())
  427 +
  428 + # Return
  429 + return new_obs
  430 +
  431 + def _pixel_analysis(self,ra,dec):
  432 + """
  433 + Performs analysis over the region of interest
  434 + corresponding to a single pixel of output map
  435 +
  436 + Parameters
  437 + ----------
  438 + inx : `int` pixel index
  439 + ra : `float` R.A. (deg)
  440 + dec : `float` Dec (deg)
  441 +
  442 + Returns
  443 + -------
  444 + result : `dict` Results
  445 + """
  446 +
  447 + # Write header for spatial bin
  448 + msg = 'Spatial bin centred on (R.A.,Dec) = (%f,%f) deg' % (ra,dec)
  449 + self._log_header2(gammalib.EXPLICIT, msg)
  450 +
  451 + # Initialize results
  452 + result = {}
  453 + for name in self._srcnames:
  454 + result[name] = {'flux' : 0.0,
  455 + 'flux_err': 0.0,
  456 + 'TS' : 0.0,
  457 + 'ulimit' : 0.0}
  458 +
  459 + # Mask observations
  460 + self._log_header3(gammalib.EXPLICIT, 'Masking observations')
  461 + masked_obs = self._mask_observations(ra,dec,self._roisz)
  462 +
  463 + # Set up likelihood analysis
  464 + self._log_header3(gammalib.EXPLICIT, 'Performing fit in energy bin')
  465 + like = ctools.ctlike(masked_obs)
  466 + like['edisp'] = self['edisp'].boolean()
  467 + like['nthreads'] = 1 # Avoids OpenMP conflict
  468 +
  469 + # If chatter level is verbose and debugging is requested then
  470 + # switch also on the debug model in ctlike
  471 + if self._logVerbose() and self._logDebug():
  472 + like['debug'] = True
  473 +
  474 + # Perform maximum likelihood fit
  475 + like.run()
  476 +
  477 + # Write model results for explicit chatter level
  478 + self._log_string(gammalib.EXPLICIT, str(like.obs().models()))
  479 +
  480 + # Prepare objects for flux extraction
  481 + # ROI
  482 + centre = gammalib.GSkyDir()
  483 + centre.radec_deg(ra,dec)
  484 + roi = gammalib.GSkyRegionCircle(centre,self._roisz)
  485 + # Energy boundaries
  486 + emin = gammalib.GEnergy(self['emin'].real(),'TeV')
  487 + emax = gammalib.GEnergy(self['emax'].real(), 'TeV')
  488 +
  489 + # Continue only if log-likelihood is non-zero
  490 + if like.obs().logL() != 0.0:
  491 +
  492 + # Loop over target sources
  493 + for name in self._srcnames:
  494 +
  495 + # Get source
  496 + source = like.obs().models()[name]
  497 +
  498 + # Get flux
  499 + # Integrate spectral model between emin and emax
  500 + flux = source.spectral().flux(emin,emax)
  501 + # Multiply by flux from spatial component in ROI
  502 + flux *= source.spatial().flux(roi)
  503 + # Divide by solid angle of ROI
  504 + flux /= 1 - math.cos(math.radians(self._roisz))
  505 + result[name]['flux'] = flux
  506 +
  507 + # Get flux error
  508 + # Normalization parameter
  509 + normpar = source.spectral()[0]
  510 + # Only normalization parameter free
  511 + # Relative error on flux is same as on normalization
  512 + # Avoid zero division error
  513 + if normpar.value() > 0.:
  514 + flux_error = flux * normpar.error()/normpar.value()
  515 + else:
  516 + flux_error = 0.
  517 + result[name]['flux_error'] = flux_error
  518 +
  519 + # If requested get TS
  520 + if self['calc_ts'].boolean():
  521 + result[name]['TS'] = source.ts()
  522 +
  523 + # Upper limit
  524 + # TODO
  525 +
  526 + else:
  527 + value = 'Likelihood is zero. Bin is skipped.'
  528 + self._log_value(gammalib.TERSE, '(R.A.,Dec) = (%f,%f) deg' % (ra,dec), value)
  529 +
  530 + # Return result
  531 + return result
  532 +
  533 + def _fill_fits(self,results):
  534 + """
  535 + Fill FITS object to store the results
  536 +
  537 + Parameters
  538 + ----------
  539 + result : `dict` Results
  540 +
  541 + Returns
  542 + -------
  543 + fits : `~gammalib.GFits` Fits object
  544 + """
  545 +
  546 + # Create new Fits object
  547 + fits = gammalib.GFits()
  548 +
  549 + # If more than one source of interest
  550 + # append empty primary HDU
  551 + if len(self._srcnames) > 1:
  552 + fits.append(gammalib.GFitsImageDouble())
  553 +
  554 + # Prepare some header cards
  555 + # Flux units
  556 + fcard = gammalib.GFitsHeaderCard('BUNIT', 'ph/cm2/s/sr' , 'Photon flux')
  557 + txt = 'Correlation radius %f deg' %(self._roisz)
  558 + roiszcard = gammalib.GFitsHeaderCard('COMMENT', '' , txt)
  559 +
  560 + # Loop over target sources
  561 + for s, name in enumerate(self._srcnames):
  562 +
  563 + # Create minimal set of skymaps to store fluxes and flux errors
  564 + fmap = self._create_map()
  565 + errmap = self._create_map()
  566 +
  567 + #Create map for TS
  568 + if self['calc_ts'].boolean():
  569 + tsmap = self._create_map()
  570 +
  571 + # Loop over pixels and fill maps
  572 + for inx in range(fmap.npix()):
  573 +
  574 + # Fill flux and flux error
  575 + fmap[inx,0] = results[inx][name]['flux']
  576 + errmap[inx, 0] = results[inx][name]['flux_error']
  577 +
  578 + # If requested fill TS
  579 + if self['calc_ts'].boolean():
  580 + tsmap[inx, 0] = results[inx][name]['TS']
  581 +
  582 + # Write maps to Fits
  583 + fhdu = fmap.write(fits, name + ' FLUX')
  584 + errhdu = errmap.write(fits, name + ' FLUX ERROR')
  585 +
  586 + # Add units and comments to headers
  587 + fhdu.header().append(fcard)
  588 + errhdu.header().append(fcard)
  589 + fhdu.header().append(roiszcard)
  590 + errhdu.header().append(roiszcard)
  591 +
  592 + # If requested write TS map to fits
  593 + if self['calc_ts'].boolean():
  594 + tshdu = tsmap.write(fits, name + ' TS')
  595 + tshdu.header().append(roiszcard)
  596 +
  597 + # Return
  598 + return fits
  599 +
96 600  
97 601 # Public methods
98 602 def run(self):
... ... @@ -106,23 +610,77 @@ class csscs(ctools.csobservation):
106 610 # Get parameters
107 611 self._get_parameters()
108 612  
109   - # TODO: Your code goes here
  613 + # Write input observation container into logger
  614 + self._log_observations(gammalib.NORMAL, self.obs(), 'Input observation')
  615 +
  616 + # Create template for output map, to use to extract pixel values
  617 + outmap = self._create_map()
  618 +
  619 + # Adjust model parameters dependent on input user parameters
  620 + self._adjust_model_pars()
  621 +
  622 + # Loop over pixels and extract results
  623 + results = {}
  624 + for inx in range(outmap.npix()):
  625 +
  626 + # Determine pixel centre
  627 + pixdir = outmap.inx2dir(inx)
  628 + ra = pixdir.ra_deg()
  629 + dec = pixdir.dec_deg()
  630 +
  631 + # Analyse ROI
  632 + result = self._pixel_analysis(ra,dec)
  633 + results[inx] = result
  634 +
  635 + # Fill output Fits
  636 + self._fits = self._fill_fits(results)
110 637  
111 638 # Return
112 639 return
113 640  
114 641 def save(self):
115 642 """
116   - Save something
  643 + Save maps to Fits flie
117 644 """
118 645 # Write header
119   - self._log_header1(gammalib.TERSE, 'Save something')
  646 + self._log_header1(gammalib.TERSE, 'Save maps')
  647 +
  648 + # Continue only if FITS file is valid
  649 + if self._fits != None:
  650 +
  651 + # Get outmap parameter
  652 + outfile = self['outfile'].filename()
120 653  
121   - # TODO: Your code goes here
  654 + # Log file name
  655 + self._log_value(gammalib.NORMAL, 'Output file', outfile.url())
  656 +
  657 + # Save spectrum
  658 + self._fits.saveto(outfile, self['clobber'].boolean())
122 659  
123 660 # Return
124 661 return
125 662  
  663 + def exclusion_map(self, object=None):
  664 + """
  665 + Return and optionally set the exclusion regions map
  666 +
  667 + Parameters
  668 + ----------
  669 + object : `~gammalib.GSkyRegion` or `~gammalib.GSkyMap` or `~gammalib.GFilename`
  670 + Exclusion regions object
  671 +
  672 + Returns
  673 + -------
  674 + region : `~gammalib.GSkyRegionMap`
  675 + Exclusion regions map
  676 + """
  677 + # If a regions object is provided then set the regions ...
  678 + if object is not None:
  679 + self._excl_reg_map = gammalib.GSkyRegionMap(object)
  680 +
  681 + # Return
  682 + return self._excl_reg_map
  683 +
126 684  
127 685 # ======================== #
128 686 # Main routine entry point #
... ...