Veronika Vodeb / ctools

Commit fba90429514122919d7f376fa68aad2b308c2727

Authored by Jürgen Knödlseder
1 parent 20133af4

Integrate code (#3248)

Inline Side-by-side
Showing 3 changed files with 91 additions and 120 deletions
cscripts/csphagen.par
  View file @fba9042
... ... @@ -43,8 +43,8 @@ glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (
43 43 glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
44 44 rad, r, a, 0.2,0,180, "Radius of source region circle (deg)"
45 45 width, r, a, 0.2,0,180, "Width of source region rectangle (deg)"
46   -height, r, a, 0.2,0,180, "Width of source region rectangle (deg)"
47   -posang, r, a, 0,-360,360, "Position angle of source region rectangle ccw from North (deg)"
  46 +height, r, a, 0.2,0,180, "Height of source region rectangle (deg)"
  47 +posang, r, a, 0.0,-360,360, "Position angle of source region rectangle (deg)"
48 48 srcregfile, f, a, on.reg,,, "Source region file (ds9 or FITS WCS)"
49 49 bkgmethod, s, a, REFLECTED,REFLECTED|CUSTOM,, "Method for background estimation"
50 50 bkgregfile, f, a, bkg.reg,,, "Background regions file (ds9 or FITS WCS)"
... ...
cscripts/csphagen.py
  View file @fba9042
... ... @@ -149,6 +149,34 @@ class csphagen(ctools.csobservation):
149 149 # Return
150 150 return
151 151  
  152 + def _compute_posang(self, pnt_dir, a, b):
  153 + """
  154 + Compute the difference in position angle wrt the pointing in degrees
  155 +
  156 + Parameters
  157 + ----------
  158 + pnt_dir : `~gammalib.GSkyDir`
  159 + Pointing direction
  160 + a : `~gammalib.GSkyDir`
  161 + First sky direction
  162 + a : `~gammalib.GSkyDir`
  163 + Second sky direction
  164 +
  165 + Returns
  166 + -------
  167 + posang : float
  168 + Position angle (degrees)
  169 + """
  170 + # Compute position angles
  171 + posang_a = pnt_dir.posang_deg(a) % 360
  172 + posang_b = pnt_dir.posang_deg(b) % 360
  173 +
  174 + # Compute difference
  175 + posang = abs(posang_a - posang_b)
  176 +
  177 + # Return position angle
  178 + return posang
  179 +
152 180 def _get_regions(self, filename):
153 181 """
154 182 Get regions from DS9 file or FITS file
... ... @@ -192,26 +220,24 @@ class csphagen(ctools.csobservation):
192 220 self['bkgregmin'].integer()
193 221 self['bkgregskip'].integer()
194 222  
195   - # Set source region (so far only CIRCLE is supported)
  223 + # If source shape is a circle the append GSkyRegionCircle
196 224 if self._srcshape == 'CIRCLE':
197 225  
198 226 # Set circular source region
199 227 self._rad = self['rad'].real()
200 228 self._src_reg.append(gammalib.GSkyRegionCircle(self._src_dir, self._rad))
201 229  
  230 + # ... otherwise if source shape is a rectangle then append GSkyRegionRect
202 231 elif self._srcshape == 'RECT':
203 232  
204 233 # Set rectangular source region
205 234 self._reg_width = self['width'].real()
206 235 self._reg_height = self['height'].real()
207 236 self._reg_posang = self['posang'].real()
208   - self._src_reg.append(gammalib.GSkyRegionRect(self._src_dir
209   - ,self._reg_width
210   - ,self._reg_height
211   - ,self._reg_posang
212   - )
213   - )
214   -
  237 + self._src_reg.append(gammalib.GSkyRegionRect(self._src_dir,
  238 + self._reg_width,
  239 + self._reg_height,
  240 + self._reg_posang))
215 241  
216 242 # Query usage of background model
217 243 self['use_model_bkg'].boolean()
... ... @@ -357,118 +383,49 @@ class csphagen(ctools.csobservation):
357 383 model = gammalib.GModelSky(spatial, spectral)
358 384 model.name('Dummy')
359 385 self._models.append(model)
360   - self._srcname = 'Dummy'
  386 + self._srcname = 'Dummy'
361 387 self['use_model_bkg'].boolean(False)
362 388  
363 389 # Return
364 390 return
365 391  
366 392 def _compute_region_separation(self, pnt_dir):
367   - """ Compute the separation angle for reflected off regions in radians.
  393 + """
  394 + Compute the separation angle for reflected off regions in radians
368 395  
369 396 Returns
370 397 -------
371   - angle : `float`
372   - Separation angle of two off regions (rad)
  398 + angle : float
  399 + Separation angle of two off regions (radians)
373 400 """
374 401 # Initialise the result
375   - separation = -1
  402 + separation = -1.0
376 403  
377 404 # Compute offset of reflected regions to pointing position
378 405 offset = pnt_dir.dist_deg(self._src_dir)
379 406  
  407 + # If shape is a circle then compute apparent diameter of the circle
  408 + # as separation
380 409 if self._srcshape == 'CIRCLE':
381   -
382   - # Compute apparent diameter of the circle
383 410 separation = 2.0 * self._rad / offset
384 411  
  412 + # ... otherwise if shape is a rectangle then compute the opening angle
  413 + # towards combinations of rectangle corners. This method overestimates
  414 + # the real need of space between the ectangles, so the method may be
  415 + # optimised to gain more off regions! Anyway, it is assured that the
  416 + # off regions will never overlap.
385 417 elif self._srcshape == 'RECT':
386 418  
387   - # Compute the opening angle towards combinations of rectangle corners.
388   - # This method overestimates the real need of space between the
389   - # rectangles, so the method may be optimised to gain more off regions!
390   - # Anyway, it is assured that the off regions will never overlap.
391   -
392   - # Get the corners
  419 + # Get the sky directions of the corners of the rectangle
393 420 cs = [self._src_reg[0].get_corner(icorner) for icorner in range(4)]
394 421  
395   - # Auxillary, local method for opening angle computation
396   - def __compute_posang(pnt_dir,a,b):
397   - """ Compute the difference in position angle to the pointing in
398   - degrees.
399   - """
400   - posang_a = pnt_dir.posang_deg(a) %360
401   - posang_b = pnt_dir.posang_deg(b) %360
402   - posang = abs(posang_a - posang_b)
403   - return posang
404   -
405 422 # Compute the 6 opening angles
406 423 combinations = [[0,1], [0,2], [0,3], [1,2], [1,3], [2,3]]
407   - angles = [__compute_posang(pnt_dir, cs[i], cs[j]) for i,j in combinations]
  424 + angles = [self._compute_posang(pnt_dir, cs[i], cs[j]) \
  425 + for i,j in combinations]
408 426  
409 427 # The desired separation is the maximum opening angle
410   - separation = max(angles) *gammalib.deg2rad
411   -
412   -
413   - ## Different approach:
414   - ## Find the separation angle for which two rectangles to not overlap
415   - ## Issue: Scanning not optimal.
416   - # posang = pnt_dir.posang_deg(self._src_dir)
417   - # dphi_incr = 0.1
418   - # dphi = dphi_incr
419   - # print("Src ctr direction: ra={} dec={}".format(self._src_dir.ra_deg(), self._src_dir.dec_deg()))
420   - # print(self._src_reg)
421   - # while dphi<360:
422   -
423   - # ctr_dir = pnt_dir.clone()
424   - # ctr_dir.rotate_deg(posang + dphi, offset)
425   - # print("Query ctr direction: ra={} dec={}".format(ctr_dir.ra_deg(), ctr_dir.dec_deg()))
426   -
427   - # # Create test rectangle
428   - # reg = gammalib.GSkyRegionRect(ctr_dir
429   - # ,self._reg_width
430   - # ,self._reg_height
431   - # ,self._reg_posang + dphi
432   - # )
433   - # print(reg)
434   -
435   - # # Check test rectangle corner containment in reference rectangle
436   - # is_corner_contained = self._src_reg.overlaps(reg)
437   - # # is_corner_contained = False
438   - # # for icorner in range(4):
439   - # # corner_dir = reg.get_corner(icorner)
440   - # # is_corner_contained |= self._src_reg.contains(corner_dir)
441   -
442   - # # Check if any corner was contained
443   - # if not is_corner_contained:
444   - # diameter = dphi
445   - # break
446   -
447   - # # Increment dphi
448   - # dphi += dphi_incr
449   -
450   -
451   - ## Very different approach: Compute the width along the circular
452   - ## region placement path.
453   - ## Issue: Long regions will overlap in the inner part.
454   - # def l(w,h,phi):
455   - # """ Super complicated method to compute the height along dec-axis
456   - # of a rotated rectangle.
457   - # """
458   - # phi = np.deg2rad(phi) % (2*np.pi)
459   - # phi_c = np.arctan2(w,h)
460   - # if phi>=(0.5*np.pi) and (phi<np.pi):
461   - # phi = np.pi-phi
462   - # elif (phi>=np.pi) and (phi<(1.5*np.pi)):
463   - # phi -=np.pi
464   - # elif (phi>=(1.5*np.pi)) and (phi<(2*np.pi)):
465   - # phi = 2*np.pi - phi
466   -
467   - # if (0<=phi) and (phi<phi_c):
468   - # result = h/2/np.cos(phi)
469   - # else:
470   - # result = w/2/np.sin(phi)
471   - # return 2*result
  428 + separation = max(angles) * gammalib.deg2rad
472 429  
473 430 # Return
474 431 return separation
... ... @@ -551,11 +508,10 @@ class csphagen(ctools.csobservation):
551 508 region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
552 509 elif self._srcshape == 'RECT':
553 510 # Adjust the posang of the rectangle correspondingly
554   - region = gammalib.GSkyRegionRect(ctr_dir
555   - ,self._reg_width
556   - ,self._reg_height
557   - ,self._reg_posang + dphi
558   - )
  511 + region = gammalib.GSkyRegionRect(ctr_dir,
  512 + self._reg_width,
  513 + self._reg_height,
  514 + self._reg_posang + dphi)
559 515 if self._has_exclusion:
560 516 if self._excl_reg.overlaps(region):
561 517  
... ...
doc/source/users/reference_manual/csphagen.rst
  View file @fba9042
... ... @@ -16,22 +16,27 @@ the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF),
16 16  
17 17 If the ``bkgmethod`` parameter is ``REFLECTED``, the script will use an On region
18 18 shape defined by the ``srcshape`` parameter and generate reflected regions by
19   -placing regions of the same shape that are reflected with respect to the pointing
20   -direction. So far only circular regions are supported, with a center that is
21   -specified by the ``ra`` and ``dec`` parameters if ``coordsys=CEL`` or the ``glon``
22   -and ``glat`` parameters if ``coordsys=GAL``, and a radius that is specified by the
23   -``rad`` parameter. Reflected regions are only generated for offset angles with
24   -respect to the pointing direction that are smaller than the value specified by
25   -the ``maxoffset`` parameter. If less than ``bkgregmin`` reflected regions can be
26   -found because the offset angle between source and pointing direction is too small,
27   -the observation will be skipped. Regions to be excluded for background
28   -determination can be specified using a FITS mask file that is specified using the
29   -``inexclusion`` parameter. A mask is a sky map where all pixel values that
30   -differ from zero will be used as exclusion region. If the FITS file
31   -contains multiple extensions and the user does not specify which one
32   -to use, the extention ``EXCLUSION`` will be used if available (e.g.,
33   -if you use a skymap generated by :ref:`ctskymap`).
34   -Otherwise the primary extension will be used.
  19 +placing regions of the same shape that are reflected with respect to the
  20 +pointing direction. So far circular and rectangular regions are supported.
  21 +Circular regions are defined by a center that is specified by the ``ra`` and
  22 +``dec`` parameters if ``coordsys=CEL`` or the ``glon`` and ``glat`` parameters if
  23 +``coordsys=GAL``, and a radius that is specified by the ``rad`` parameter.
  24 +Rectangular regions are defined by a center, again specified by either the
  25 +``ra`` and ``dec`` or the ``glon`` and ``glat`` parameters, a ``width`` parameter,
  26 +a ``height`` parameter and a ``posang`` parameter that specified the position
  27 +angle measured counter clockwise from celestial North. Note that the ``height``
  28 +axis is pointing in the direction of the position angle. Reflected regions are
  29 +only generated for offset angles with respect to the pointing direction that are
  30 +smaller than the value specified by the ``maxoffset`` parameter. If less than
  31 +``bkgregmin`` reflected regions can be found because the offset angle between
  32 +source and pointing direction is too small, the observation will be skipped.
  33 +Regions to be excluded for background determination can be specified using a
  34 +FITS mask file that is specified using the ``inexclusion`` parameter. A mask is
  35 +a sky map where all pixel values that differ from zero will be used as exclusion
  36 +region. If the FITS file contains multiple extensions and the user does not
  37 +specify which one to use, the extention ``EXCLUSION`` will be used if available
  38 +(e.g., if you use a skymap generated by :ref:`ctskymap`). Otherwise the primary
  39 +extension will be used.
35 40  
36 41 Alternatively, if ``bkgmethod=CUSTOM`` is specified, the source and background
37 42 regions will be specified using the ``srcregfile`` and ``bkgregfile`` parameters,
... ... @@ -132,9 +137,9 @@ General parameters
132 137 Exponent of the power law for ``POW`` energy binning. An exponent of 1.0
133 138 corresponds to a logarithmic energy binning.
134 139  
135   -``(srcshape = CIRCLE) <CIRCLE> [string]``
136   - Shape of the source region. So far only ``CIRCLE`` exists which defines a
137   - circular region around given position.
  140 +``srcshape <CIRCLE|RECT> [string]``
  141 + Shape of the source region. ``CIRCLE`` defines a circular region around given
  142 + position, ``RECT`` defines a rectangle centred on a given position.
138 143  
139 144 ``coordsys <CEL|GAL> [string]``
140 145 Coordinate system (CEL - celestial, GAL - galactic).
... ... @@ -154,6 +159,16 @@ General parameters
154 159 ``rad [real]``
155 160 Radius of source region circle (deg).
156 161  
  162 +``width [real]``
  163 + Width of source region rectangle (deg).
  164 +
  165 +``height [real]``
  166 + Height of source region rectangle (deg).
  167 +
  168 +``posang [real]``
  169 + Position angle of source region rectangle, counted counterclockwise from
  170 + celestial North (deg).
  171 +
157 172 ``srcregfile [file]``
158 173 Source region file (ds9 or FITS WCS map).
159 174  
... ...