Veronika Vodeb / ctools

Commit 110dc2797a91a3f183656ea6fe83f3d0a06267e0

Authored by Jürgen Knödlseder
1 parent 7863238e

Integrate code

Inline Side-by-side
Showing 3 changed files with 198 additions and 176 deletions
cscripts/csscs.par
  View file @110dc27
... ... @@ -59,7 +59,6 @@ calc_ulim, b, h, no,,, "Compute upper limit for each map bin"
59 59 fix_srcs, b, h, yes,,, "Fix normalisations of other skymodels"
60 60 fix_bkg, b, h, no,,, "Fix background normalisation"
61 61  
62   -
63 62 #
64 63 # Standard parameters
65 64 #====================
... ...
cscripts/csscs.py
  View file @110dc27
... ... @@ -26,9 +26,9 @@ from cscripts import mputils
26 26 from cscripts import obsutils
27 27  
28 28  
29   -# ================= #
  29 +# =========== #
30 30 # csscs class #
31   -# ================= #
  31 +# =========== #
32 32 class csscs(ctools.csobservation):
33 33 """
34 34 Spectral component separation script
... ... @@ -49,11 +49,11 @@ class csscs(ctools.csobservation):
49 49 self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
50 50  
51 51 # Initialise data members
52   - self._nthreads = 0
53   - self._srcnames = []
54   - self._method = None
55   - self._fits = None
56   - self._excl_reg_map = None # Exclusion region map for on/off analysis
  52 + self._nthreads = 0
  53 + self._srcnames = []
  54 + self._method = None
  55 + self._fits = None
  56 + self._excl_reg_map = None # Exclusion region map for on/off analysis
57 57 self._excl_reg_name = None
58 58  
59 59 # Return
... ... @@ -70,13 +70,12 @@ class csscs(ctools.csobservation):
70 70 Pickled instance
71 71 """
72 72 # Set pickled dictionary
73   - # Set pickled dictionary
74 73 state = {'base': ctools.csobservation.__getstate__(self),
75   - 'nthreads': self._nthreads,
76   - 'method': self._method,
77   - 'srcnames': self._srcnames,
78   - 'fits': self._fits,
79   - 'excl_reg_map': self._excl_reg_map,
  74 + 'nthreads': self._nthreads,
  75 + 'srcnames': self._srcnames,
  76 + 'method': self._method,
  77 + 'fits': self._fits,
  78 + 'excl_reg_map': self._excl_reg_map,
80 79 'excl_reg_name' : self._excl_reg_name}
81 80  
82 81 # Return pickled dictionary
... ... @@ -93,11 +92,11 @@ class csscs(ctools.csobservation):
93 92 """
94 93 # Set state
95 94 ctools.csobservation.__setstate__(self, state['base'])
96   - self._nthreads = state['nthreads']
97   - self._srcnames = state['srcnames']
98   - self._method = state['method']
99   - self._fits = state['fits']
100   - self._excl_reg_map = state['excl_reg_map']
  95 + self._nthreads = state['nthreads']
  96 + self._srcnames = state['srcnames']
  97 + self._method = state['method']
  98 + self._fits = state['fits']
  99 + self._excl_reg_map = state['excl_reg_map']
101 100 self._excl_reg_name = state['excl_reg_name']
102 101  
103 102 # Return
... ... @@ -107,11 +106,12 @@ class csscs(ctools.csobservation):
107 106 def _get_onoff_parameters(self):
108 107 """
109 108 Get On/Off analysis parameters.
110   - This is done here rather then in ctools.is_onoff
111   - because the exclusion map needs to be handled differently
112   - as an automatic parameter but queried only if not already set
113   - and we need to verify is not empty.
  109 +
  110 + This is done here rather then in ctools.is_onoff because the exclusion
  111 + map needs to be handled differently as an automatic parameter but
  112 + queried only if not already set and we need to verify is not empty.
114 113 Also many parameters are already queried for all methods
  114 +
115 115 TODO: verify if this can be made more uniform with other scripts
116 116 """
117 117 # Exclusion map
... ... @@ -126,22 +126,22 @@ class csscs(ctools.csobservation):
126 126 raise RuntimeError(msg)
127 127  
128 128 # Other csphagen parameters
129   - self["enumbins"].integer()
130   - if self["bkgmethod"].string() == "REFLECTED":
131   - self["bkgregmin"].integer()
  129 + self['enumbins'].integer()
  130 + if self['bkgmethod'].string() == 'REFLECTED':
  131 + self['bkgregmin'].integer()
132 132 self['bkgregskip'].integer()
133   - self["maxoffset"].real()
134   - self["etruemin"].real()
135   - self["etruemax"].real()
136   - self["etruebins"].integer()
  133 + self['maxoffset'].real()
  134 + self['etruemin'].real()
  135 + self['etruemax'].real()
  136 + self['etruebins'].integer()
137 137  
  138 + # Return
138 139 return
139 140  
140 141 def _get_parameters(self):
141 142 """
142 143 Get parameters from parfile
143 144 """
144   -
145 145 # Set observation if not done before
146 146 if self.obs().is_empty():
147 147 self.obs(self._get_observations())
... ... @@ -152,8 +152,8 @@ class csscs(ctools.csobservation):
152 152 # Collect number of unbinned, binned and On/Off observations in
153 153 # observation container
154 154 n_unbinned = 0
155   - n_binned = 0
156   - n_onoff = 0
  155 + n_binned = 0
  156 + n_onoff = 0
157 157 for obs in self.obs():
158 158 if obs.classname() == 'GCTAObservation':
159 159 if obs.eventtype() == 'CountsCube':
... ... @@ -162,7 +162,7 @@ class csscs(ctools.csobservation):
162 162 n_unbinned += 1
163 163 elif obs.classname() == 'GCTAOnOffObservation':
164 164 n_onoff += 1
165   - n_cta = n_unbinned + n_binned + n_onoff
  165 + n_cta = n_unbinned + n_binned + n_onoff
166 166 n_other = self.obs().size() - n_cta
167 167  
168 168 # Check if there are On/Off or non-IACT observations
... ... @@ -170,7 +170,8 @@ class csscs(ctools.csobservation):
170 170 msg = 'On/Off or non-CTA observations found. ' \
171 171 'csscs does not support this type of observations.'
172 172 raise RuntimeError(msg)
173   - # Otherwise check if there is a mix of binned and unbinned
  173 +
  174 + # ... otherwise check if there is a mix of binned and unbinned
174 175 elif n_unbinned > 0 and n_binned > 0:
175 176 msg = 'Mix of unbinned and binned CTA observations ' \
176 177 'found in observation container. csscs does not ' \
... ... @@ -190,12 +191,12 @@ class csscs(ctools.csobservation):
190 191 self['nypix'].integer()
191 192 self['binsz'].real()
192 193  
193   - # Radius of ROI for component separation
  194 + # Query radius of ROI for component separation
194 195 self['rad'].real()
195 196  
196   - # Energy boundaries
197   - self["emin"].real()
198   - self["emax"].real()
  197 + # Query energy boundaries
  198 + self['emin'].real()
  199 + self['emax'].real()
199 200  
200 201 # If we have unbinned observations query analysis method
201 202 if n_unbinned > 0:
... ... @@ -210,9 +211,9 @@ class csscs(ctools.csobservation):
210 211 # Query target source names
211 212 srcnames = self['srcnames'].string()
212 213  
213   - # Fashion target sources into Python list
  214 + # Fashion target sources into Python list. Strip leading and trailing
  215 + # spaces from names.
214 216 self._srcnames = srcnames.split(';')
215   - # Strip leading and trailing spaces from names
216 217 for s, name in enumerate(self._srcnames):
217 218 self._srcnames[s] = name.strip()
218 219  
... ... @@ -226,18 +227,23 @@ class csscs(ctools.csobservation):
226 227  
227 228 # For On/Off analysis check the number of gamma-ray sources in the model
228 229 if self._method == 'ONOFF':
  230 +
  231 + # Initialise number of sources
229 232 nsources = 0
  233 +
230 234 # Loop over models and verify if model is of type GModelSky
231 235 for model in self.obs().models():
232 236 if model.classname() == 'GModelSky':
233 237 nsources += 1
  238 +
234 239 # If there are background gamma-ray sources in the model
235 240 # throw runtime error
236 241 if nsources > len(self._srcnames):
237 242 msg = 'Background gamma-ray sources found in the model. ' \
238 243 'On/Off analysis does not support this feature.'
239 244 raise RuntimeError(msg)
240   - # Otherwise continue
  245 +
  246 + # ... otherwise continue
241 247 else:
242 248 pass
243 249  
... ... @@ -268,14 +274,16 @@ class csscs(ctools.csobservation):
268 274 Returns
269 275 -------
270 276 skymap : `~gammalib.GSkyMap`
  277 + Sky map
271 278 """
272   -
  279 + # Create sky map
273 280 skymap = gammalib.GSkyMap(self['proj'].string(), self['coordsys'].string(),
274 281 self['xref'].real(), self['yref'].real(),
275 282 -self['binsz'].real(), self['binsz'].real(),
276 283 self['nxpix'].integer(), self['nypix'].integer(),
277 284 1)
278 285  
  286 + # Return sky map
279 287 return skymap
280 288  
281 289 def _adjust_model_pars(self):
... ... @@ -298,13 +306,15 @@ class csscs(ctools.csobservation):
298 306 # Thus, if we have more than one target we will set them all
299 307 # to isotropic
300 308 if self._method == 'ONOFF' and len(self._srcnames) > 1:
  309 +
301 310 # Check if it is a source
302 311 try:
303 312 gammalib.GModelSky(model)
304 313 # In this case check if the spatial model is already isotropic
305 314 if model.spatial() == 'DiffuseIsotropic':
306 315 pass
307   - # Otherwise change it to isotropic
  316 +
  317 + # ... otherwise change it to isotropic
308 318 else:
309 319 msg = ' Setting spatial model to diffuse isotropic'
310 320 self._log_string(gammalib.EXPLICIT, msg)
... ... @@ -360,17 +370,20 @@ class csscs(ctools.csobservation):
360 370  
361 371 Parameters
362 372 ----------
363   - obs : `gammalib.GCTAObservation` Input observation of type cube
364   - ra : `float` R.A. (deg)
365   - dec : `float` Dec (deg)
366   - rad : `float` radius (deg)
  373 + obs : `gammalib.GCTAObservation`
  374 + Input observation of type cube
  375 + ra : float
  376 + Right Ascension (deg)
  377 + dec : float
  378 + Declination (deg)
  379 + rad : float
  380 + Radius (deg)
367 381  
368 382 Returns
369 383 -------
370 384 new_obs : `gammalib.GCTAObservations`
371 385 Output observations with masked cubes
372 386 """
373   -
374 387 # Filter cube according to ROI and user energy range
375 388 cubemask = ctools.ctcubemask(obs)
376 389 cubemask['ra'] = ra
... ... @@ -400,27 +413,30 @@ class csscs(ctools.csobservation):
400 413  
401 414 Parameters
402 415 ----------
403   - obs : `gammalib.GCTAObservations` Input observations of type event list
404   - ra : `float` R.A. (deg)
405   - dec : `float` Dec (deg)
406   - rad : `float` radius (deg)
  416 + obs : `gammalib.GCTAObservations`
  417 + Input observations of type event list
  418 + ra : float
  419 + Right Ascension (deg)
  420 + dec : float
  421 + Declination (deg)
  422 + rad : float
  423 + Radius (deg)
407 424  
408 425 Returns
409 426 -------
410 427 new_obs : `gammalib.GCTAObservations`
411 428 Output observations with masked event lists
412 429 """
413   -
414 430 # Filter event list according to ROI and user energy range
415 431 select = ctools.ctselect(obs)
416 432 select['usepnt'] = False
417   - select['ra'] = ra
418   - select['dec'] = dec
419   - select['rad'] = rad
420   - select['emin'] = self['emin'].real()
421   - select['emax'] = self['emax'].real()
422   - select['tmin'] = 'INDEF'
423   - select['tmax'] = 'INDEF'
  433 + select['ra'] = ra
  434 + select['dec'] = dec
  435 + select['rad'] = rad
  436 + select['emin'] = self['emin'].real()
  437 + select['emax'] = self['emax'].real()
  438 + select['tmin'] = 'INDEF'
  439 + select['tmax'] = 'INDEF'
424 440  
425 441 # If chatter level is verbose and debugging is requested then
426 442 # switch also on the debug model in ctcubemask
... ... @@ -442,16 +458,18 @@ class csscs(ctools.csobservation):
442 458  
443 459 Parameters
444 460 ----------
445   - ra : `float` R.A. (deg)
446   - dec : `float` Dec (deg)
447   - rad : `float` radius (deg)
  461 + ra : float
  462 + Right Ascension (deg)
  463 + dec : float
  464 + Declination (deg)
  465 + rad : float
  466 + Radius (deg)
448 467  
449 468 Returns
450 469 -------
451 470 new_obs : `~gammalib.GObservations`
452 471 Observations in circular ROI
453 472 """
454   -
455 473 # Determine type of masking according
456 474 # to observation type and analysis method
457 475 if self._method == 'UNBINNED':
... ... @@ -468,16 +486,18 @@ class csscs(ctools.csobservation):
468 486  
469 487 def _pixel_analysis(self, inx):
470 488 """
471   - Performs analysis over the region of interest
472   - corresponding to a single pixel of output map
  489 + Performs analysis over the region of interest corresponding to a
  490 + single pixel of output map
473 491  
474 492 Parameters
475 493 ----------
476   - inx : `int` pixel index
  494 + inx : int
  495 + Pixel index
477 496  
478 497 Returns
479 498 -------
480   - result : `dict` Results
  499 + result : dict
  500 + Results
481 501 """
482 502  
483 503 # Create template for output map to extract pixel coordinates
... ... @@ -485,8 +505,8 @@ class csscs(ctools.csobservation):
485 505  
486 506 # Determine pixel centre coordinates
487 507 pixdir = outmap.inx2dir(inx)
488   - ra = pixdir.ra_deg()
489   - dec = pixdir.dec_deg()
  508 + ra = pixdir.ra_deg()
  509 + dec = pixdir.dec_deg()
490 510  
491 511 # Write header for spatial bin
492 512 msg = 'Spatial bin (%d) centred on (R.A.,Dec) = (%f,%f) deg' % (inx, ra, dec)
... ... @@ -510,7 +530,7 @@ class csscs(ctools.csobservation):
510 530 # Set up likelihood analysis
511 531 self._log_header3(gammalib.NORMAL, 'Performing fit in region')
512 532 like = ctools.ctlike(masked_obs)
513   - like['edisp'] = self['edisp'].boolean()
  533 + like['edisp'] = self['edisp'].boolean()
514 534 like['nthreads'] = 1 # Avoids OpenMP conflict
515 535  
516 536 # If chatter level is verbose and debugging is requested then
... ... @@ -532,6 +552,7 @@ class csscs(ctools.csobservation):
532 552 centre = gammalib.GSkyDir()
533 553 centre.radec_deg(ra, dec)
534 554 roi = gammalib.GSkyRegionCircle(centre, self['rad'].real())
  555 +
535 556 # Energy boundaries
536 557 emin = gammalib.GEnergy(self['emin'].real(), 'TeV')
537 558 emax = gammalib.GEnergy(self['emax'].real(), 'TeV')
... ... @@ -545,18 +566,20 @@ class csscs(ctools.csobservation):
545 566 # Get flux
546 567 # Integrate spectral model between emin and emax
547 568 flux = source.spectral().flux(emin, emax)
  569 +
548 570 # Calculate correction factor
549 571 # Spatial model flux over ROI divided by ROI solid angle
550   - corr_factor = source.spatial().flux(roi)
551   - corr_factor /= gammalib.twopi * (
552   - 1 - math.cos(math.radians(self['rad'].real())))
  572 + corr_factor = source.spatial().flux(roi)
  573 + corr_factor /= gammalib.twopi * (1 - math.cos(math.radians(self['rad'].real())))
  574 +
553 575 # Multiply flux by correction factor
554   - flux *= corr_factor
  576 + flux *= corr_factor
555 577 result[name]['flux'] = flux
556 578  
557 579 # Get flux error
558 580 # Normalization parameter
559 581 normpar = source.spectral()[0]
  582 +
560 583 # Only normalization parameter free
561 584 # Relative error on flux is same as on normalization
562 585 # Avoid zero division error
... ... @@ -618,17 +641,17 @@ class csscs(ctools.csobservation):
618 641  
619 642 Parameters
620 643 ----------
621   - hdu : `~gammalib.GFitsHDU` HDU
  644 + hdu : `~gammalib.GFitsHDU`
  645 + HDU
622 646 """
623   -
624 647 # Flux units
625 648 if is_flux:
626 649 hdu.card('BUNIT', 'ph/cm2/s/sr', 'Photon flux')
627 650  
628 651 # Minimum and maximum energy
629   - hdu.card("E_MIN", self['emin'].real(), "[TeV] Lower energy boundary")
630   - hdu.card("E_MAX", self['emax'].real(), "[TeV] Upper energy boundary")
631   - hdu.card("EUNIT", "TeV", "Units for E_MIN and E_MAX")
  652 + hdu.card('E_MIN', self['emin'].real(), '[TeV] Lower energy boundary')
  653 + hdu.card('E_MAX', self['emax'].real(), '[TeV] Upper energy boundary')
  654 + hdu.card('EUNIT', 'TeV', 'Units for E_MIN and E_MAX')
632 655  
633 656 # ROI size
634 657 hdu.card('ROISZ', self['rad'].real(),
... ... @@ -642,14 +665,15 @@ class csscs(ctools.csobservation):
642 665  
643 666 # Exclusion map
644 667 if self._excl_reg_name is not None:
645   - hdu.card("EXCLMAP", self._excl_reg_name.url(), "Exclusion map name")
  668 + hdu.card('EXCLMAP', self._excl_reg_name.url(), 'Exclusion map name')
646 669  
647 670 # Model background?
648   - hdu.card("BKGMODEL", self['use_model_bkg'].boolean(), 'Use background model?')
  671 + hdu.card('BKGMODEL', self['use_model_bkg'].boolean(), 'Use background model?')
649 672  
650 673 # Fit statistic
651   - hdu.card("STAT", self['statistic'].string(), 'Fit statistic')
  674 + hdu.card('STAT', self['statistic'].string(), 'Fit statistic')
652 675  
  676 + # Return
653 677 return
654 678  
655 679 def _fill_fits(self, results):
... ... @@ -658,13 +682,14 @@ class csscs(ctools.csobservation):
658 682  
659 683 Parameters
660 684 ----------
661   - result : `dict` Results
  685 + result : `dict`
  686 + Results
662 687  
663 688 Returns
664 689 -------
665   - fits : `~gammalib.GFits` Fits object
  690 + fits : `~gammalib.GFits`
  691 + FITS object
666 692 """
667   -
668 693 # Create new Fits object
669 694 fits = gammalib.GFits()
670 695  
... ... @@ -677,18 +702,17 @@ class csscs(ctools.csobservation):
677 702 # Create skymaps to store results
678 703 # Maps for TS and upper limits created and filled to avoid if statements
679 704 # Will not be saved in the end if not requested
680   - fmap = self._create_map()
  705 + fmap = self._create_map()
681 706 errmap = self._create_map()
682   - tsmap = self._create_map()
683   - ulmap = self._create_map()
  707 + tsmap = self._create_map()
  708 + ulmap = self._create_map()
684 709  
685 710 # Loop over pixels and fill maps
686 711 for inx in range(fmap.npix()):
687   - # Fill maps
688   - fmap[inx, 0] = results[inx][name]['flux']
  712 + fmap[inx, 0] = results[inx][name]['flux']
689 713 errmap[inx, 0] = results[inx][name]['flux_error']
690   - tsmap[inx, 0] = results[inx][name]['TS']
691   - ulmap[inx, 0] = results[inx][name]['ulimit']
  714 + tsmap[inx, 0] = results[inx][name]['TS']
  715 + ulmap[inx, 0] = results[inx][name]['ulimit']
692 716  
693 717 # Uppercase name to adhere to gammalib convention
694 718 Name = gammalib.toupper(name)
... ... @@ -738,12 +762,16 @@ class csscs(ctools.csobservation):
738 762  
739 763 # Proceed only if fits has been filled
740 764 if self._fits != None:
  765 +
741 766 # Check that the name is in the list of target sources
742 767 if name in self._srcnames:
  768 +
743 769 # Uppercase name
744 770 Name = gammalib.toupper(name)
  771 +
745 772 # Assemble HDU name
746 773 hduname = Name + ' ' + quantity
  774 +
747 775 # Try to fetch HDU
748 776 try:
749 777 hdu = self._fits[hduname]
... ... @@ -753,7 +781,8 @@ class csscs(ctools.csobservation):
753 781 except:
754 782 msg = 'HDU "' + hduname + '" not found in FITS container.'
755 783 raise RuntimeError(msg)
756   - # Otherwise throw name error
  784 +
  785 + # ... otherwise throw name error
757 786 else:
758 787 print('ERROR: Source "' + name + '" not found in list of target sources.')
759 788 raise NameError(name)
... ... @@ -843,7 +872,7 @@ class csscs(ctools.csobservation):
843 872 region : `~gammalib.GSkyRegionMap`
844 873 Exclusion regions map
845 874 """
846   - # If a regions object is provided then set the regions ...
  875 + # If a regions object is provided then set the regions
847 876 if object is not None:
848 877 self._excl_reg_map = gammalib.GSkyRegionMap(object)
849 878  
... ... @@ -856,7 +885,7 @@ class csscs(ctools.csobservation):
856 885  
857 886 Returns
858 887 -------
859   - fits : `~gammalib.GFits()`
  888 + fits : `~gammalib.GFits`
860 889 FITS file containing the maps
861 890 """
862 891 # Return
... ... @@ -876,9 +905,10 @@ class csscs(ctools.csobservation):
876 905 skymap : `~gammalib.GSkyMap`
877 906 Flux skymap
878 907 """
879   -
  908 + # Get flux skymap
880 909 skymap = self._get_skymap(name, 'FLUX')
881 910  
  911 + # Return flux sky map
882 912 return skymap
883 913  
884 914 def flux_error(self, name):
... ... @@ -895,9 +925,10 @@ class csscs(ctools.csobservation):
895 925 skymap : `~gammalib.GSkyMap`
896 926 Flux error skymap
897 927 """
898   -
  928 + # Get flux error sky map
899 929 skymap = self._get_skymap(name, 'FLUX ERROR')
900 930  
  931 + # Return flux error sky map
901 932 return skymap
902 933  
903 934 def ts(self, name):
... ... @@ -914,17 +945,17 @@ class csscs(ctools.csobservation):
914 945 skymap : `~gammalib.GSkyMap`
915 946 TS skymap
916 947 """
917   -
918 948 # Check that TS computation is requested
919 949 if self['calc_ts'].boolean():
920 950 skymap = self._get_skymap(name, 'TS')
921   - # Otherwise throw error
  951 +
  952 + # ... otherwise raise error
922 953 else:
923   - msg = 'TS computation not requested. ' \
924   - 'Change user parameter calc_ts to True ' \
925   - 'to obtain TS.'
  954 + msg = 'TS computation not requested. Change user parameter ' \
  955 + '"calc_ts" to True to obtain TS.'
926 956 raise RuntimeError(msg)
927 957  
  958 + # Return TS map
928 959 return skymap
929 960  
930 961 def ulimit(self, name):
... ... @@ -941,17 +972,17 @@ class csscs(ctools.csobservation):
941 972 skymap : `~gammalib.GSkyMap`
942 973 Flux upper limit skymap
943 974 """
944   -
945 975 # Check that upper limit computation is requested
946 976 if self['calc_ulim'].boolean():
947 977 skymap = self._get_skymap(name, 'FLUX UPPER LIMIT')
948   - # Otherwise throw error
  978 +
  979 + # ... otherwise raise error
949 980 else:
950   - msg = 'Upper limit computation not requested. ' \
951   - 'Change user parameter calc_ulimit to True ' \
952   - 'to obtain upper limit.'
  981 + msg = 'Upper limit computation not requested. Change user parameter ' \
  982 + '"calc_ulimit" to True to obtain upper limit.'
953 983 raise RuntimeError(msg)
954 984  
  985 + # Return upper limit map
955 986 return skymap
956 987  
957 988  
... ... @@ -959,6 +990,7 @@ class csscs(ctools.csobservation):
959 990 # Main routine entry point #
960 991 # ======================== #
961 992 if __name__ == '__main__':
  993 +
962 994 # Create instance of application
963 995 app = csscs(sys.argv)
964 996  
... ...
doc/source/users/reference_manual/csscs.rst
  View file @110dc27
1 1 .. _csscs:
2 2  
3 3 csscs
4   -==============================
  4 +=====
5 5  
6 6 Performs spectral component separation.
7 7  
... ... @@ -9,50 +9,43 @@ Performs spectral component separation.
9 9 Synopsis
10 10 --------
11 11  
12   -This script performs a spectral component separation, that is,
13   -produces a flux skymap of one or multiple sources by running a likelihood
14   -analysis over multiple adjacent circular regions of interest (ROI)
15   -based on prior knowledge of the sources' spectra. The centre of the
16   -ROIs is displaced to coincide with the centre of each pixel in the
17   -output flux map, while the radius is set through the user parameter ``rad``.
18   -
19   -The model fit per ROI is performed using :doc:`ctlike`, and only the
20   -normalisations of the model components are treated as free parameters.
21   -The script provides the possibility to fix sources other than the
22   -source(s) of interest (hidden parameter ``fix_srcs``) or to fix the
23   -normalisation of the background
24   -model component(s) (hidden parameter ``fix_bkg``).
25   -
26   -The script computes the source flux and its uncertainty in each
27   -ROI/map pixel, as well as the significance of the source detection
28   -(for faster execution the significance computation can be
29   -deactivated using the hidden parameter ``calc_ts``). Optionally, the
30   -script also computes an upper flux limit that is particularly useful
31   -in case that the source is not significantly detected within an
32   -ROI (hidden parameter ``calc_ulim``).
33   -
34   -The script accepts in input event lists, binned data, or observation
35   -definition XML files. Only data
36   -within an energy interval spanned by ``emin`` and ``emax`` are
37   -considered. For event lists two analysis methods are proposed:
38   -unbinned or On/Off.
39   -
40   -For the On/Off method :doc:`csphagen` is
41   -used to prepare On/Off observations for every ROI
42   -analysed. In On/Off mode the user is required to provide an exclusion
43   -map to avoid using regions of significant gamma-ray emission in the
44   -evaluation of the background (exclusion maps can be generate using
45   -:doc:`ctskymap`). Furthemore, On/Off analysis including
46   -sources other than the source(s) of interest is not supported,
47   -and if more than one source of interest is
48   -present the response will be calculated based on the
49   -hypothesis that emission from the sources of interest is uniform over
50   -each ROI.
51   -
52   -:ref:`csscs` generates a FITS file with empty primary extension.
53   -Additional extensions provide for every source of interest a skymap of
54   -flux, flux error, and, if requested, detection signficance and flux
55   -upper limits.
  12 +This script performs a spectral component separation, that is, produces a flux
  13 +skymap of one or multiple sources by running a likelihood analysis over multiple
  14 +adjacent circular regions of interest (ROI) based on prior knowledge of the
  15 +sources' spectra. The centre of the ROIs is displaced to coincide with the
  16 +centre of each pixel in the output flux map, while the radius is set through
  17 +the user parameter ``rad``.
  18 +
  19 +The model fit per ROI is performed using :doc:`ctlike`, and only the normalisations
  20 +of the model components are treated as free parameters. The script provides the
  21 +possibility to fix sources other than the source(s) of interest (hidden
  22 +parameter ``fix_srcs``) or to fix the normalisation of the background model
  23 +component(s) (hidden parameter ``fix_bkg``).
  24 +
  25 +The script computes the source flux and its uncertainty in each ROI/map pixel,
  26 +as well as the significance of the source detection (for faster execution the
  27 +significance computation can be deactivated using the hidden parameter
  28 +``calc_ts``). Optionally, the script also computes an upper flux limit that is
  29 +particularly useful in case that the source is not significantly detected within
  30 +an ROI (hidden parameter ``calc_ulim``).
  31 +
  32 +The script accepts in input event lists, binned data, or observation definition
  33 +XML files. Only data within an energy interval spanned by ``emin`` and ``emax`` are
  34 +considered. For event lists two analysis methods are proposed: unbinned or
  35 +On/Off.
  36 +
  37 +For the On/Off method :doc:`csphagen` is used to prepare On/Off observations for
  38 +every ROI analysed. In On/Off mode the user is required to provide an exclusion
  39 +map to avoid using regions of significant gamma-ray emission in the evaluation
  40 +of the background (exclusion maps can be generate using :doc:`ctskymap`). Furthemore,
  41 +On/Off analysis including sources other than the source(s) of interest is not
  42 +supported, and if more than one source of interest is present the response will
  43 +be calculated based on the hypothesis that emission from the sources of interest
  44 +is uniform over each ROI.
  45 +
  46 +:ref:`csscs` generates a FITS file with empty primary extension. Additional
  47 +extensions provide for every source of interest a skymap of flux, flux error,
  48 +and, if requested, detection signficance and flux upper limits.
56 49  
57 50 General parameters
58 51 ------------------
... ... @@ -64,9 +57,8 @@ General parameters
64 57 Input model XML file.
65 58  
66 59 ``srcnames [string]``
67   - Semicolon-separated list of names of the sources in the source
68   - model XML file which should be considered for the spectral
69   - component separation.
  60 + Semicolon-separated list of names of the sources in the source model XML
  61 + file which should be considered for the spectral component separation.
70 62  
71 63 ``expcube [file]``
72 64 Input exposure cube file.
... ... @@ -88,10 +80,10 @@ General parameters
88 80  
89 81 ``inexclusion [file]``
90 82 FITS file containing a WCS map that defines sky regions not to be used for
91   - background estimation (where map value != 0). If the file
92   - contains multiple extensions the user may specify which one
93   - to use. Otherwise, the extention ``EXCLUSION`` will be used if
94   - available, or else the primary extension will be used.
  83 + background estimation (where map value != 0). If the file contains multiple
  84 + extensions the user may specify which one to use. Otherwise, the extention
  85 + ``EXCLUSION`` will be used if available, or else the primary extension will
  86 + be used.
95 87  
96 88 ``(edisp = no) [boolean]``
97 89 Apply energy dispersion to response computation?
... ... @@ -122,9 +114,9 @@ General parameters
122 114 degrees).
123 115  
124 116 ``rad [real]``
125   - Radius of region of interest for component separation (deg). Sets
126   - the correlation scale between neighbour pixels in the output maps. Must
127   - be at least sqrt(2) times ``binsz`` for full coverage of input data.
  117 + Radius of region of interest for component separation (deg). Sets the
  118 + correlation scale between neighbour pixels in the output maps. Must be at
  119 + least sqrt(2) times ``binsz`` for full coverage of input data.
128 120  
129 121 ``emin [real]``
130 122 Minimum energy (in TeV).
... ... @@ -133,17 +125,16 @@ General parameters
133 125 Maximum energy (in TeV).
134 126  
135 127 ``method <UNBINNED|ONOFF> [string]``
136   - For input event lists selects between ``UNBINNED`` analysis
137   - (3D spatial/energy likelihood without binning) and ``ONOFF``
138   - analysis (1D likelihood with background from Off regions).
  128 + For input event lists selects between ``UNBINNED`` analysis (3D spatial/energy
  129 + likelihood without binning) and ``ONOFF`` analysis (1D likelihood with
  130 + background from Off regions).
139 131  
140 132 ``(statistic = DEFAULT) <DEFAULT|CSTAT|WSTAT|CHI2> [string]``
141 133 Optimization statistic. ``DEFAULT`` uses the default statistic for all
142 134 observations, which is ``CSTAT`` or the statistic specified in the
143   - observation definition XML file. ``CSTAT`` uses the C statistic for
144   - all observations, ``WSTAT`` uses the W statistic for On/Off
145   - observations, and ``CHI2`` uses the Chi squared statistic for
146   - binned or stacked observations.
  135 + observation definition XML file. ``CSTAT`` uses the C statistic for all
  136 + observations, ``WSTAT`` uses the W statistic for On/Off observations, and
  137 + ``CHI2`` uses the Chi squared statistic for binned or stacked observations.
147 138  
148 139 ``(calc_ts = yes) [boolean]``
149 140 Compute TS value for each map bin?
... ... @@ -164,10 +155,9 @@ On/Off parameters
164 155 Number of energy bins.
165 156  
166 157 ``(bkgmethod = REFLECTED) [string]``
167   - Method for background estimation.
168   - ``REFLECTED:`` background evaluated in regions with the same shape as
169   - source region reflected w.r.t. pointing direction for each
170   - observation.
  158 + Method for background estimation. ``REFLECTED:`` background evaluated in regions
  159 + with the same shape as source region reflected w.r.t. pointing direction for
  160 + each observation.
171 161  
172 162 ``(srcshape = CIRCLE) [string]``
173 163 Shape of the source region.
... ... @@ -235,7 +225,8 @@ Standard parameters
235 225 Specifies whether an existing output file should be overwritten.
236 226  
237 227 ``(debug = no) [boolean]``
238   - Enables debug mode. In debug mode the executable will dump any log file output to the console.
  228 + Enables debug mode. In debug mode the executable will dump any log file
  229 + output to the console.
239 230  
240 231 ``(mode = ql) [string]``
241 232 Mode of automatic parameters (default is ``ql``, i.e. "query and learn").
... ...