ctools  2.0.0
 All Classes Namespaces Files Functions Variables Macros Pages
modutils.py
Go to the documentation of this file.
1 # ==========================================================================
2 # Utility functions for model handling
3 #
4 # Copyright (C) 2016-2019 Juergen Knoedlseder
5 #
6 # This program is free software: you can redistribute it and/or modify
7 # it under the terms of the GNU General Public License as published by
8 # the Free Software Foundation, either version 3 of the License, or
9 # (at your option) any later version.
10 #
11 # This program is distributed in the hope that it will be useful,
12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 # GNU General Public License for more details.
15 #
16 # You should have received a copy of the GNU General Public License
17 # along with this program. If not, see <http://www.gnu.org/licenses/>.
18 #
19 # ==========================================================================
20 #import gammalib
21 #import ctools
22 
23 
24 # ============================================ #
25 # Return model for TS fitting of a test source #
26 # ============================================ #
27 def test_source(models, srcname, ra=None, dec=None, fitspat=False,
28  fitspec=False):
29  """
30  Return model for TS fitting of a test source
31 
32  Parameters
33  ----------
34  models : `~gammalib.GModels`
35  Input model container
36  srcname : str
37  Test source name
38  ra : float, optional
39  Right Ascension of test source (deg)
40  dec : float, optional
41  Declination of test source (deg)
42  fitspat : bool, optional
43  Fit spatial parameters of all models?
44  fitspec : bool, optional
45  Fit spectral parameters of all models?
46 
47  Returns
48  -------
49  model : `~gammalib.GModels`
50  Model container for TS fitting
51  """
52  # Create a clone of the input model
53  outmodels = models.copy()
54 
55  # Disable TS computation for all model components (will enable the
56  # test source later)
57  for model in outmodels:
58  model.tscalc(False)
59 
60  # Get source model and enable TS computation
61  model = outmodels[srcname]
62  model.tscalc(True)
63 
64  # Get spectral model component
65  spectral = model.spectral()
66 
67  # If source model has no normalisation parameter then raise an exception
68  if not (spectral.has_par('Normalization') or
69  spectral.has_par('Value') or
70  spectral.has_par('Prefactor') or
71  spectral.has_par('Integral') or
72  spectral.has_par('PhotonFlux') or
73  spectral.has_par('EnergyFlux')):
74  msg = ('Model "%s" has no flux normalisation parameter. Only spectral '
75  'models with a flux normalisation parameter (one of '
76  '"Normalization", "Value", "Prefactor", "Integral", '
77  '"PhotonFlux", and "EnergyFlux") are supported.' %
78  srcname)
79  raise RuntimeError(msg)
80 
81  # Set position of test source
82  if ra != None and dec != None:
83  if model.has_par('RA') and model.has_par('DEC'):
84  model['RA'].value(ra)
85  model['DEC'].value(dec)
86 
87  # Set possible spatial and spectral parameters
88  spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
89  'MinorRadius', 'MajorRadius']
90  spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
91  'InverseCutoffEnergy']
92 
93  # Fit or fix spatial parameters
94  for par in spatial:
95  if model.has_par(par):
96  if fitspat:
97  model[par].free()
98  else:
99  model[par].fix()
100 
101  # Fit or fix spectral parameters
102  for par in spectral:
103  if model.has_par(par):
104  if fitspec:
105  model[par].free()
106  else:
107  model[par].fix()
108 
109  # Return model container
110  return outmodels
111 
112 
113 # ============================================= #
114 # Return spectral model normalisation parameter #
115 # ============================================= #
117  """
118  Return spectral model normalisation parameter
119 
120  Parameters
121  ----------
122  model : `~gammalib.GModel`
123  Input model
124 
125  Returns
126  -------
127  par : `~gammalib.GModelPar`
128  Spectral model normalisation parameter
129  """
130  # Initialise normalisation parameter as None
131  par = None
132 
133  # Set list of possible normalisation parameters
134  norm_pars = ['Normalization', 'Value', 'Prefactor', 'Integral',
135  'PhotonFlux', 'EnergyFlux']
136 
137  # Get spectral model component
138  spectral = model.spectral()
139 
140  # Find normalisation parameter
141  for norm_par in norm_pars:
142  if spectral.has_par(norm_par):
143  par = spectral[norm_par]
144 
145  # If no normalisation parameter then raise an exception
146  if par == None:
147  msg = ('Model "%s" has no flux normalisation parameter. Only spectral '
148  'models with a flux normalisation parameter (one of '
149  '"Normalization", "Value", "Prefactor", "Integral", '
150  '"PhotonFlux", and "EnergyFlux") are supported.' %
151  model.name())
152  raise RuntimeError(msg)
153 
154  # Return normalisation parameter
155  return par
156 
157 
158 # ================================= #
159 # Returns attributes for ds9 string #
160 # ================================= #
161 def ds9attributes(model, free_color='green', fixed_color='magenta',
162  width=2, fontfamily='helvetica', fontsize=12,
163  fontweight='normal', fontslant='roman'):
164  """
165  Returns attributes for DS9 string
166 
167  Parameters
168  ----------
169  model : `~gammalib.GModel`
170  Model
171  free_color : str, optional
172  Color for sources with free parameters (any DS9 color or hex code)
173  fixed_color : str, optional
174  Color for source without free parameters (any DS9 color or hex code)
175  width : int, optional
176  Line width for regions
177  fontfamily : str, optional
178  Font family for source labels (helvetica,times,courier)
179  fontsize : int, optional
180  Font size for source labels
181  fontweight : str, optional
182  Font weight for source labels (normal,bold)
183  fontslant : str, optional
184  Font slant for source labels (roman,italic)
185 
186  Returns
187  -------
188  attributes : str
189  DS9 attribute string
190  """
191  # Determine region color. A model where all parameters are fixed will
192  # get the color defined by "fixed_color", a model with at least one
193  # free parameter will get the color defined by "free_color".
194  color = fixed_color
195  for par in model:
196  if par.is_free():
197  color = free_color
198  break
199 
200  # Set DS9 attributes
201  attributes = (' color=%s width=%d font="%s %d %s %s"' %
202  (color, width, fontfamily, fontsize, fontweight, fontslant))
203 
204  # Return attributes
205  return attributes
206 
207 
208 # =========================== #
209 # Convert model to ds9 string #
210 # =========================== #
211 def model2ds9string(model, pnt_type='cross', pnt_mark_size=12,
212  show_labels=True, show_ext_type=True,
213  free_color='green', fixed_color='magenta',
214  width=2, fontfamily='helvetica', fontsize=12,
215  fontweight='normal', fontslant='roman'):
216  """
217  Converts model into a DS9 region string
218 
219  Parameters
220  ----------
221  model : `~gammalib.GModel`
222  Model
223  pnt_type : str, optional
224  Marker type for point sources (circle,box,diamond,cross,x,arrow,boxcircle)
225  pnt_mark_size : integer, optional
226  Marker size for point sources
227  show_labels : boolean, optional
228  Add source labels?
229  show_ext_type : boolean, optional
230  Show type of extended model in source name?
231  free_color : str, optional
232  Color for sources with free parameters (any DS9 color or hex code)
233  fixed_color : str, optional
234  Color for source without free parameters (any DS9 color or hex code)
235  width : int, optional
236  Line width for regions
237  fontfamily : str, optional
238  Font family for source labels (helvetica,times,courier)
239  fontsize : int, optional
240  Font size for source labels
241  fontweight : str, optional
242  Font weight for source labels (normal,bold)
243  fontslant : str, optional
244  Font slant for source labels (roman,italic)
245 
246  Returns
247  -------
248  ds9string : str
249  DS9 region string
250  msg : str
251  Error message
252  """
253  # Initialise DS9 region string and message string
254  ds9string = ''
255  msg = ''
256 
257  # Retrieve model sky direction. The model is skipped in case that
258  # it does not provide a sky direction.
259  is_valid = True
260  try:
261  modelpos = model.spatial().dir()
262  except AttributeError:
263  msg = ('Skip model "%s" since it has no sky direction.\n' %
264  model.name())
265  is_valid = False
266 
267  # Continue only if sky direction was found
268  if is_valid:
269 
270  # Retrieve model type and name
271  modeltype = model.type()
272  modelname = model.name()
273 
274  # Handle point source
275  if modeltype == 'PointSource':
276 
277  # Append point with Right Ascension and Declination to the DS9
278  # string. The type of the point is specified by the "pnt_type"
279  # parameter, the size of the point by the "pnt_mark_size"
280  # parameter
281  ds9string += ('point(%.6f,%.6f) # point=%s %d' %
282  (modelpos.ra_deg(), modelpos.dec_deg(),
283  pnt_type, pnt_mark_size))
284 
285  # Handle extended sources
286  elif modeltype == "ExtendedSource":
287 
288  # Retrieve spatial model
289  spatial = model.spatial()
290  classname = spatial.classname()
291 
292  # Handle radial sources
293  if 'Radial' in classname:
294 
295  # Retrieve short name of model class (e.g. "Gauss"
296  # or "Disk)
297  shorttype = classname.split('Radial')[-1]
298 
299  # Handle Disk and Shell model
300  if (classname == 'GModelSpatialRadialDisk' or
301  classname == 'GModelSpatialRadialShell'):
302  size = spatial.radius()
303 
304  # Handle Gauss Model
305  elif classname == 'GModelSpatialRadialGauss':
306  size = spatial.sigma()
307 
308  # Skip if source is unknown
309  else:
310  msg = ('Skip model "%s" since the radial model "%s" '
311  'is unknown.\n' % (model.name(), classname))
312  is_valid = False
313 
314  # Append circle to DS9 string
315  if is_valid:
316  ds9string += ('circle(%.6f,%.6f,%.6f") #' %
317  (modelpos.ra_deg(), modelpos.dec_deg(),
318  size*3600.0))
319 
320  # Handle elliptical sources
321  elif 'Elliptical' in classname:
322 
323  # Retrieve short name and source size
324  shorttype = classname.split('Elliptical')[-1]
325  size1 = spatial.semimajor()
326  size2 = spatial.semiminor()
327  angle = spatial.posangle()
328 
329  # Append ellipse to DS9 string
330  ds9string += ('ellipse(%.6f,%.6f,%.6f",%.6f",%.6f) #' %
331  (modelpos.ra_deg(), modelpos.dec_deg(),
332  size1*3600.0, size2*3600.0, angle+90.0))
333 
334  # Skip if source is neither radial nor elliptical
335  else:
336  msg = ('Skip model "%s" since the model "%s" is neither '
337  'a point source, a radial source, nor an elliptical '
338  'source.\n' % (model.name(), classname))
339  is_valid = False
340 
341  # Add short model type to modelname
342  if show_ext_type:
343  modelname +=' ('+shorttype+')'
344 
345  # Add DS9 attributes
346  if is_valid:
347  ds9string += ds9attributes(model, free_color=free_color,
348  fixed_color=fixed_color,
349  width=width,
350  fontfamily=fontfamily,
351  fontsize=fontsize,
352  fontweight=fontweight,
353  fontslant=fontslant)
354  if show_labels:
355  ds9string += ' text={'+modelname+'}'
356  else:
357  ds9string = ''
358 
359  # Return ds9 and message strings
360  return ds9string, msg
361 
362 
363 # ========================= #
364 # Save models into ds9 file #
365 # ========================= #
366 def models2ds9file(models, filename,
367  pnt_type='cross', pnt_mark_size=12,
368  show_labels=True, show_ext_type=True,
369  free_color='green', fixed_color='magenta',
370  width=2, fontfamily='helvetica', fontsize=12,
371  fontweight='normal', fontslant='roman'):
372  """
373  Save models in a DS9 region file
374 
375  Parameters
376  ----------
377  model : `~gammalib.GModel`
378  Model
379  filename : str
380  DS9 region file name
381  pnt_type : str, optional
382  Marker type for point sources (circle,box,diamond,cross,x,arrow,boxcircle)
383  pnt_mark_size : integer, optional
384  Marker size for point sources
385  show_labels : boolean, optional
386  Add source labels?
387  show_ext_type : boolean, optional
388  Show type of extended model in source name?
389  free_color : str, optional
390  Color for sources with free parameters (any DS9 color or hex code)
391  fixed_color : str, optional
392  Color for source without free parameters (any DS9 color or hex code)
393  width : int, optional
394  Line width for regions
395  fontfamily : str, optional
396  Font family for source labels (helvetica,times,courier)
397  fontsize : int, optional
398  Font size for source labels
399  fontweight : str, optional
400  Font weight for source labels (normal,bold)
401  fontslant : str, optional
402  Font slant for source labels (roman,italic)
403 
404  Returns
405  -------
406  errors : str
407  Error message
408  """
409  # Initialise error string
410  errors = ''
411 
412  # Open file
413  f = open(filename, 'w')
414 
415  # Write coordinate system
416  f.write('fk5\n')
417 
418  # Loop over models
419  for model in models:
420 
421  # Continue only if point source or extended source model
422  if (model.type() == 'PointSource' or
423  model.type() == 'ExtendedSource'):
424  line, msg = model2ds9string(model,
425  pnt_type=pnt_type,
426  pnt_mark_size=pnt_mark_size,
427  show_labels=show_labels,
428  show_ext_type=show_ext_type,
429  free_color=free_color,
430  fixed_color=fixed_color,
431  width=width,
432  fontfamily=fontfamily,
433  fontsize=fontsize,
434  fontweight=fontweight,
435  fontslant=fontslant)
436  if len(line):
437  f.write(line+'\n')
438  elif len(msg):
439  errors += msg+'\n'
440 
441  # Logging for diffuse components
442  elif model.type() == 'DiffuseSource':
443  errors += 'Skipping diffuse model "'+model.name()+'".\n'
444 
445  # Logging for background components
446  else:
447  errors += 'Skipping background model "'+model.name()+'".\n'
448 
449  # Close file
450  f.close()
451 
452  # Return error message string
453  return errors
def normalisation_parameter
Definition: modutils.py:116