Commit 854ebb52715e1d9497171404fbf826f567bd9cbf

Authored by Luigi Tibaldo
1 parent ad59645c

Improved algorithm to find reflected regions in presence of exclusion regions

Showing 1 changed file with 20 additions and 7 deletions
cscripts/csphagen.py
... ... @@ -374,7 +374,7 @@ class csphagen(ctools.csobservation):
374 374 # Determine number of background regions to skip
375 375 N_skip = self['bkgregskip'].integer()
376 376 N_start = 1 + N_skip
377   - N_lim = 1 + 2*N_skip
  377 + N_lim = 1 + 2 * N_skip
378 378  
379 379 # Compute the angular separation of reflected regions wrt
380 380 # camera center. The factor 1.05 ensures background regions
... ... @@ -397,14 +397,11 @@ class csphagen(ctools.csobservation):
397 397 # regions
398 398 else:
399 399  
400   - # Log appending of reflected regions
401   - msg = ' Use %d reflected regions.' % (N-N_lim)
402   - self._log_string(gammalib.NORMAL, msg)
403   -
404 400 # Append reflected regions
405 401 alpha = 360.0 / N
406   - for s in range(N_start, N - N_skip):
407   - dphi = s * alpha
  402 + dphi_max = 360. - alpha * (1 + N_skip)
  403 + dphi = alpha * (1 + N_skip)
  404 + while dphi <= dphi_max:
408 405 ctr_dir = pnt_dir.clone()
409 406 ctr_dir.rotate_deg(posang + dphi, offset)
410 407 region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
... ... @@ -413,10 +410,19 @@ class csphagen(ctools.csobservation):
413 410 msg = ' Reflected region overlaps with '\
414 411 'exclusion region.'
415 412 self._log_string(gammalib.EXPLICIT, msg)
  413 + # if region overlaps with exclusion region
  414 + # try to increment by 10% of angular step
  415 + dphi += 0.1 * alpha
416 416 else:
417 417 regions.append(region)
  418 + dphi += alpha
418 419 else:
419 420 regions.append(region)
  421 + dphi += alpha
  422 +
  423 + # Log number of reflected regions
  424 + msg = ' Use %d reflected regions.' % (regions.size())
  425 + self._log_string(gammalib.NORMAL, msg)
420 426  
421 427 # Return reflected regions
422 428 return regions
... ... @@ -670,6 +676,13 @@ class csphagen(ctools.csobservation):
670 676 # Set On/Off observation ID
671 677 onoff.id(obs.id())
672 678  
  679 + # Otherwise log observation skipped
  680 + else:
  681 + msg = ' Skip because the number %d of regions ' \
  682 + 'for background estimation is smaller than ' \
  683 + '"bkgregmin"=%d.' % (bkg_reg.size(), self['bkgregmin'].integer())
  684 + self._log_string(gammalib.NORMAL, msg)
  685 +
673 686 # Construct dictionary with results
674 687 result = {'onoff' : onoff,
675 688 'bkg_reg' : bkg_reg,
... ...