Commit a63b9b9249bf6083b38125855a687a59591e562b
1 parent
b8ffb0f0
Implement Off runs (unit test not working)
Showing
5 changed files
with
1554 additions
and
68 deletions
cscripts/csphagen.par
| @@ -17,6 +17,7 @@ | @@ -17,6 +17,7 @@ | ||
| 17 | # Input and output parameters | 17 | # Input and output parameters |
| 18 | #============================ | 18 | #============================ |
| 19 | inobs, f, a, obs.xml,,, "Input event list or observation definition XML file" | 19 | inobs, f, a, obs.xml,,, "Input event list or observation definition XML file" |
| 20 | +inobsoff, f, a, obs_off.xml,,, "Input Off event list or observation definition XML file" | ||
| 20 | inmodel, f, a, NONE,,, "Input model definition XML file (if NONE, use point source)" | 21 | inmodel, f, a, NONE,,, "Input model definition XML file (if NONE, use point source)" |
| 21 | srcname, s, a, Crab,,, "Source name" | 22 | srcname, s, a, Crab,,, "Source name" |
| 22 | caldb, s, a, prod2,,, "Calibration database" | 23 | caldb, s, a, prod2,,, "Calibration database" |
| @@ -46,10 +47,10 @@ width, r, a, 0.2,0,180, "Width of source region rectangle (deg)" | @@ -46,10 +47,10 @@ width, r, a, 0.2,0,180, "Width of source region rectangle (deg)" | ||
| 46 | height, r, a, 0.2,0,180, "Height of source region rectangle (deg)" | 47 | 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 | posang, r, a, 0.0,-360,360, "Position angle of source region rectangle (deg)" |
| 48 | srcregfile, f, a, on.reg,,, "Source region file (ds9 or FITS WCS)" | 49 | srcregfile, f, a, on.reg,,, "Source region file (ds9 or FITS WCS)" |
| 49 | -bkgmethod, s, a, REFLECTED,REFLECTED|CUSTOM,, "Method for background estimation" | 50 | +bkgmethod, s, a, REFLECTED,REFLECTED|CUSTOM|OFF,, "Method for background estimation" |
| 50 | bkgregfile, f, a, bkg.reg,,, "Background regions file (ds9 or FITS WCS)" | 51 | bkgregfile, f, a, bkg.reg,,, "Background regions file (ds9 or FITS WCS)" |
| 51 | -bkgregmin, i, h, 2,1,1000, "Minimum number of background regions" | ||
| 52 | -bkgregskip, i, h, 1,0,10, "Number of background regions to skip next to On region" | 52 | +bkgregmin, i, h, 2,1,1000, "Minimum number of background regions for REFLECTED method" |
| 53 | +bkgregskip, i, h, 1,0,10, "Number of background regions to skip next to On region for REFLECTED method" | ||
| 53 | use_model_bkg, b, h, yes,,, "Use background model?" | 54 | use_model_bkg, b, h, yes,,, "Use background model?" |
| 54 | maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation" | 55 | maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation" |
| 55 | stack, b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?" | 56 | stack, b, a, no,,, "Stack multiple observations into single PHA, ARF and RMF files?" |
cscripts/csphagen.py
| @@ -43,6 +43,7 @@ class csphagen(ctools.csobservation): | @@ -43,6 +43,7 @@ class csphagen(ctools.csobservation): | ||
| 43 | self._init_csobservation(self.__class__.__name__, ctools.__version__, argv) | 43 | self._init_csobservation(self.__class__.__name__, ctools.__version__, argv) |
| 44 | 44 | ||
| 45 | # Initialise other variables | 45 | # Initialise other variables |
| 46 | + self._obs_off = gammalib.GObservations() | ||
| 46 | self._ebounds = gammalib.GEbounds() | 47 | self._ebounds = gammalib.GEbounds() |
| 47 | self._etruebounds = gammalib.GEbounds() | 48 | self._etruebounds = gammalib.GEbounds() |
| 48 | self._src_dir = gammalib.GSkyDir() | 49 | self._src_dir = gammalib.GSkyDir() |
| @@ -74,6 +75,7 @@ class csphagen(ctools.csobservation): | @@ -74,6 +75,7 @@ class csphagen(ctools.csobservation): | ||
| 74 | """ | 75 | """ |
| 75 | # Set pickled dictionary | 76 | # Set pickled dictionary |
| 76 | state = {'base' : ctools.csobservation.__getstate__(self), | 77 | state = {'base' : ctools.csobservation.__getstate__(self), |
| 78 | + 'obs_off' : self._obs_off, | ||
| 77 | 'ebounds' : self._ebounds, | 79 | 'ebounds' : self._ebounds, |
| 78 | 'etruebounds' : self._etruebounds, | 80 | 'etruebounds' : self._etruebounds, |
| 79 | 'src_dir' : self._src_dir, | 81 | 'src_dir' : self._src_dir, |
| @@ -103,6 +105,7 @@ class csphagen(ctools.csobservation): | @@ -103,6 +105,7 @@ class csphagen(ctools.csobservation): | ||
| 103 | Pickled instance | 105 | Pickled instance |
| 104 | """ | 106 | """ |
| 105 | ctools.csobservation.__setstate__(self, state['base']) | 107 | ctools.csobservation.__setstate__(self, state['base']) |
| 108 | + self._obs_off = state['obs_off'] | ||
| 106 | self._ebounds = state['ebounds'] | 109 | self._ebounds = state['ebounds'] |
| 107 | self._etruebounds = state['etruebounds'] | 110 | self._etruebounds = state['etruebounds'] |
| 108 | self._src_dir = state['src_dir'] | 111 | self._src_dir = state['src_dir'] |
| @@ -205,21 +208,17 @@ class csphagen(ctools.csobservation): | @@ -205,21 +208,17 @@ class csphagen(ctools.csobservation): | ||
| 205 | # Return region container | 208 | # Return region container |
| 206 | return regs | 209 | return regs |
| 207 | 210 | ||
| 208 | - def _get_parameters_bkgmethod_reflected(self): | 211 | + def _get_source_parameters(self): |
| 209 | """ | 212 | """ |
| 210 | - Get parameters for REFLECTED background method | 213 | + Get parameters to define source/On region |
| 211 | """ | 214 | """ |
| 215 | + | ||
| 212 | # Get source shape | 216 | # Get source shape |
| 213 | self._srcshape = self['srcshape'].string() | 217 | self._srcshape = self['srcshape'].string() |
| 214 | 218 | ||
| 215 | # Query source direction | 219 | # Query source direction |
| 216 | self._query_src_direction() | 220 | self._query_src_direction() |
| 217 | 221 | ||
| 218 | - # Query minimum number of background regions and number of | ||
| 219 | - # background regions to skip next to On region | ||
| 220 | - self['bkgregmin'].integer() | ||
| 221 | - self['bkgregskip'].integer() | ||
| 222 | - | ||
| 223 | # If source shape is a circle the append GSkyRegionCircle | 222 | # If source shape is a circle the append GSkyRegionCircle |
| 224 | if self._srcshape == 'CIRCLE': | 223 | if self._srcshape == 'CIRCLE': |
| 225 | 224 | ||
| @@ -240,8 +239,21 @@ class csphagen(ctools.csobservation): | @@ -240,8 +239,21 @@ class csphagen(ctools.csobservation): | ||
| 240 | self._reg_height, | 239 | self._reg_height, |
| 241 | self._reg_posang)) | 240 | self._reg_posang)) |
| 242 | 241 | ||
| 243 | - # Query usage of background model | ||
| 244 | - self['use_model_bkg'].boolean() | 242 | + # Return |
| 243 | + return | ||
| 244 | + | ||
| 245 | + def _get_parameters_bkgmethod_reflected(self): | ||
| 246 | + """ | ||
| 247 | + Get parameters for REFLECTED background method | ||
| 248 | + """ | ||
| 249 | + | ||
| 250 | + # Query parameters for source/On region definition | ||
| 251 | + self._get_source_parameters() | ||
| 252 | + | ||
| 253 | + # Query minimum number of background regions and | ||
| 254 | + # number of background regions to skip next to On region | ||
| 255 | + self['bkgregmin'].integer() | ||
| 256 | + self['bkgregskip'].integer() | ||
| 245 | 257 | ||
| 246 | # Return | 258 | # Return |
| 247 | return | 259 | return |
| @@ -284,6 +296,58 @@ class csphagen(ctools.csobservation): | @@ -284,6 +296,58 @@ class csphagen(ctools.csobservation): | ||
| 284 | # Return | 296 | # Return |
| 285 | return | 297 | return |
| 286 | 298 | ||
| 299 | + def _get_parameters_bkgmethod_off(self): | ||
| 300 | + """ | ||
| 301 | + Get parameters for OFF background method | ||
| 302 | + | ||
| 303 | + Raises | ||
| 304 | + ------ | ||
| 305 | + RuntimeError | ||
| 306 | + On and Off observations must have same size | ||
| 307 | + """ | ||
| 308 | + | ||
| 309 | + # Set up Off observations | ||
| 310 | + # If there are no Off observations in container then load them via user parameters | ||
| 311 | + if self.obs_off().is_empty(): | ||
| 312 | + filename = self['inobsoff'].filename() | ||
| 313 | + # If it is a FITS file load into observation and append to container | ||
| 314 | + if gammalib.GFilename(filename).is_fits(): | ||
| 315 | + self._obs_off.append(gammalib.GCTAObservation(filename)) | ||
| 316 | + # Otherwise load file into container | ||
| 317 | + else: | ||
| 318 | + self._obs_off.load(filename) | ||
| 319 | + | ||
| 320 | + # Check that size of On and Off observations are the same, otherwise throw error | ||
| 321 | + if self.obs().size() != self._obs_off.size(): | ||
| 322 | + raise RuntimeError('On and Off observations must have the same size') | ||
| 323 | + | ||
| 324 | + # Loop through observations | ||
| 325 | + for obs in self.obs_off(): | ||
| 326 | + # Check that observation is event list, otherwise throw error | ||
| 327 | + if obs.eventtype() != "EventList": | ||
| 328 | + raise RuntimeError('Off observations must be event lists') | ||
| 329 | + # Check that they have response, otherwise assign based on user parameter | ||
| 330 | + if obs.has_response() == False: | ||
| 331 | + database = self["caldb"].string() | ||
| 332 | + irf = self["irf"].string() | ||
| 333 | + # xml element for response | ||
| 334 | + response_xml = "parameter name=\"Calibration\"" +\ | ||
| 335 | + " database=\"" + database + "\"" +\ | ||
| 336 | + " response=\"" + irf + "\"" | ||
| 337 | + print(response_xml) | ||
| 338 | + xml = gammalib.GXmlElement() | ||
| 339 | + xml.append(response_xml) | ||
| 340 | + print(xml) | ||
| 341 | + # define response and assign to observation | ||
| 342 | + response = gammalib.GCTAResponseIrf(xml) | ||
| 343 | + obs.response(response) | ||
| 344 | + | ||
| 345 | + # Query parameters for source/On region definition | ||
| 346 | + self._get_source_parameters() | ||
| 347 | + | ||
| 348 | + # Return | ||
| 349 | + return | ||
| 350 | + | ||
| 287 | def _get_parameters_bkgmethod(self): | 351 | def _get_parameters_bkgmethod(self): |
| 288 | """ | 352 | """ |
| 289 | Get background method parameters | 353 | Get background method parameters |
| @@ -296,9 +360,12 @@ class csphagen(ctools.csobservation): | @@ -296,9 +360,12 @@ class csphagen(ctools.csobservation): | ||
| 296 | self._get_parameters_bkgmethod_reflected() | 360 | self._get_parameters_bkgmethod_reflected() |
| 297 | elif bkgmethod == 'CUSTOM': | 361 | elif bkgmethod == 'CUSTOM': |
| 298 | self._get_parameters_bkgmethod_custom() | 362 | self._get_parameters_bkgmethod_custom() |
| 363 | + elif bkgmethod == 'OFF': | ||
| 364 | + self._get_parameters_bkgmethod_off() | ||
| 299 | 365 | ||
| 300 | # Query parameters that are needed for all background methods | 366 | # Query parameters that are needed for all background methods |
| 301 | self['maxoffset'].real() | 367 | self['maxoffset'].real() |
| 368 | + self['use_model_bkg'].boolean() | ||
| 302 | 369 | ||
| 303 | # Return | 370 | # Return |
| 304 | return | 371 | return |
| @@ -351,9 +418,9 @@ class csphagen(ctools.csobservation): | @@ -351,9 +418,9 @@ class csphagen(ctools.csobservation): | ||
| 351 | else: | 418 | else: |
| 352 | self._has_exclusion = False | 419 | self._has_exclusion = False |
| 353 | 420 | ||
| 354 | - # Query remaining parameters | ||
| 355 | - self['use_model_bkg'].boolean() | ||
| 356 | - self['stack'].boolean() | 421 | + # If there are multiple observations query whether to stack them |
| 422 | + if self.obs().size() > 1: | ||
| 423 | + self['stack'].boolean() | ||
| 357 | 424 | ||
| 358 | # Query ahead output parameters | 425 | # Query ahead output parameters |
| 359 | if (self._read_ahead()): | 426 | if (self._read_ahead()): |
| @@ -361,11 +428,6 @@ class csphagen(ctools.csobservation): | @@ -361,11 +428,6 @@ class csphagen(ctools.csobservation): | ||
| 361 | self['outmodel'].filename() | 428 | self['outmodel'].filename() |
| 362 | self['prefix'].string() | 429 | self['prefix'].string() |
| 363 | 430 | ||
| 364 | - # If there are no observations in container then get them from the | ||
| 365 | - # parameter file | ||
| 366 | - if self.obs().is_empty(): | ||
| 367 | - self.obs(self._get_observations(False)) | ||
| 368 | - | ||
| 369 | # Get background method parameters (have to come after setting up of | 431 | # Get background method parameters (have to come after setting up of |
| 370 | # observations) | 432 | # observations) |
| 371 | self._get_parameters_bkgmethod() | 433 | self._get_parameters_bkgmethod() |
| @@ -459,14 +521,6 @@ class csphagen(ctools.csobservation): | @@ -459,14 +521,6 @@ class csphagen(ctools.csobservation): | ||
| 459 | if self._srcshape == 'CIRCLE': | 521 | if self._srcshape == 'CIRCLE': |
| 460 | msg += ' (circle rad=%.3f).' % (self._rad) | 522 | msg += ' (circle rad=%.3f).' % (self._rad) |
| 461 | self._log_string(gammalib.NORMAL, msg) | 523 | self._log_string(gammalib.NORMAL, msg) |
| 462 | - | ||
| 463 | - # Skip observation if it is pointed too far from the source | ||
| 464 | - elif offset >= self['maxoffset'].real(): | ||
| 465 | - msg = ' Skip because observation is pointed at %.3f deg >= '\ | ||
| 466 | - '"maxoffset=%.3f" from source.' \ | ||
| 467 | - % (offset, self['maxoffset'].real()) | ||
| 468 | - self._log_string(gammalib.NORMAL, msg) | ||
| 469 | - | ||
| 470 | # ... otherwise | 524 | # ... otherwise |
| 471 | else: | 525 | else: |
| 472 | posang = pnt_dir.posang_deg(self._src_dir) | 526 | posang = pnt_dir.posang_deg(self._src_dir) |
| @@ -474,7 +528,6 @@ class csphagen(ctools.csobservation): | @@ -474,7 +528,6 @@ class csphagen(ctools.csobservation): | ||
| 474 | 528 | ||
| 475 | # Determine number of background regions to skip | 529 | # Determine number of background regions to skip |
| 476 | N_skip = self['bkgregskip'].integer() | 530 | N_skip = self['bkgregskip'].integer() |
| 477 | - N_start = 1 + N_skip | ||
| 478 | N_lim = 1 + 2 * N_skip | 531 | N_lim = 1 + 2 * N_skip |
| 479 | 532 | ||
| 480 | # Compute the angular separation of reflected regions wrt | 533 | # Compute the angular separation of reflected regions wrt |
| @@ -532,13 +585,78 @@ class csphagen(ctools.csobservation): | @@ -532,13 +585,78 @@ class csphagen(ctools.csobservation): | ||
| 532 | regions.append(region) | 585 | regions.append(region) |
| 533 | dphi += alpha | 586 | dphi += alpha |
| 534 | 587 | ||
| 535 | - # Log number of reflected regions | ||
| 536 | - msg = ' Use %d reflected regions.' % (regions.size()) | ||
| 537 | - self._log_string(gammalib.NORMAL, msg) | 588 | + # Check again that we have enough background regions |
| 589 | + # now that we have checked for overlap with exclusion region | ||
| 590 | + if regions.size() >= self['bkgregmin'].integer(): | ||
| 591 | + # Log number of reflected regions | ||
| 592 | + msg = ' Use %d reflected regions.' % (regions.size()) | ||
| 593 | + self._log_string(gammalib.NORMAL, msg) | ||
| 594 | + # Otherwise log observation skipped and return empty region container | ||
| 595 | + else: | ||
| 596 | + msg = ' Skip because the number %d of regions' \ | ||
| 597 | + 'for background estimation not overlapping ' \ | ||
| 598 | + 'with the exclusion region is smaller than ' \ | ||
| 599 | + '"bkgregmin"=%d.' % (bkg_reg.size(), self['bkgregmin'].integer()) | ||
| 600 | + self._log_string(gammalib.NORMAL, msg) | ||
| 601 | + regions = gammalib.GSkyRegions() | ||
| 538 | 602 | ||
| 539 | # Return reflected regions | 603 | # Return reflected regions |
| 540 | return regions | 604 | return regions |
| 541 | 605 | ||
| 606 | + def _instrument_regions(self, obs, obs_off): | ||
| 607 | + """ | ||
| 608 | + Calculate background region in Off observation that corresponds to the source region | ||
| 609 | + in the On observation in instrument coordinates | ||
| 610 | + | ||
| 611 | + Parameters | ||
| 612 | + ---------- | ||
| 613 | + obs : `~gammalib.GCTAObservation()` | ||
| 614 | + On CTA observation | ||
| 615 | + obs_off : `~gammalib.GCTAObservation()` | ||
| 616 | + Off CTA observation | ||
| 617 | + | ||
| 618 | + Returns | ||
| 619 | + ------- | ||
| 620 | + regions : `~gammalib.GSkyRegions` | ||
| 621 | + Container with background region | ||
| 622 | + """ | ||
| 623 | + # Initialise region container | ||
| 624 | + regions = gammalib.GSkyRegions() | ||
| 625 | + | ||
| 626 | + # Convert source position in On observation to instrument coordinates | ||
| 627 | + instdir = obs.pointing().instdir(self._src_dir) | ||
| 628 | + | ||
| 629 | + # Convert instrument position to sky direction for Off observation | ||
| 630 | + off_dir = obs_off.pointing().skydir(instdir) | ||
| 631 | + | ||
| 632 | + # Build region according to shape specified by user | ||
| 633 | + # If circle | ||
| 634 | + if self._srcshape == 'CIRCLE': | ||
| 635 | + region = gammalib.GSkyRegionCircle(off_dir, self._rad) | ||
| 636 | + # ... otherwise if rectangle | ||
| 637 | + elif self._srcshape == 'RECT': | ||
| 638 | + # Instrument coordinates take sky direction as reference | ||
| 639 | + # so no need to change the position angle | ||
| 640 | + region = gammalib.GSkyRegionRectangle(off_dir, | ||
| 641 | + self._reg_width, | ||
| 642 | + self._reg_height, | ||
| 643 | + self._reg_posang) | ||
| 644 | + | ||
| 645 | + # Check if background region overlaps with exclusion region | ||
| 646 | + is_valid = True | ||
| 647 | + if self._has_exclusion: | ||
| 648 | + if self._excl_reg.overlaps(region): | ||
| 649 | + # Signal region overlap | ||
| 650 | + msg = ' Background region overlaps with ' \ | ||
| 651 | + 'exclusion region.' | ||
| 652 | + self._log_string(gammalib.EXPLICIT, msg) | ||
| 653 | + is_valid = False | ||
| 654 | + # ... otherwise store region in container | ||
| 655 | + regions.append(region) | ||
| 656 | + | ||
| 657 | + # Return | ||
| 658 | + return regions | ||
| 659 | + | ||
| 542 | def _set_models(self, results): | 660 | def _set_models(self, results): |
| 543 | """ | 661 | """ |
| 544 | Set models for On/Off fitting | 662 | Set models for On/Off fitting |
| @@ -587,7 +705,7 @@ class csphagen(ctools.csobservation): | @@ -587,7 +705,7 @@ class csphagen(ctools.csobservation): | ||
| 587 | # observations | 705 | # observations |
| 588 | for result in results: | 706 | for result in results: |
| 589 | if model.is_valid(result['instrument'], result['id']): | 707 | if model.is_valid(result['instrument'], result['id']): |
| 590 | - if result['bkg_reg'].size() >= self['bkgregmin'].integer(): | 708 | + if result['bkg_reg'].size() > 0: |
| 591 | use_model = True | 709 | use_model = True |
| 592 | break | 710 | break |
| 593 | 711 | ||
| @@ -699,7 +817,7 @@ class csphagen(ctools.csobservation): | @@ -699,7 +817,7 @@ class csphagen(ctools.csobservation): | ||
| 699 | # Return energy boundaries | 817 | # Return energy boundaries |
| 700 | return ebounds | 818 | return ebounds |
| 701 | 819 | ||
| 702 | - def _set_background_regions(self, obs): | 820 | + def _set_background_regions(self, obs, obs_off=None): |
| 703 | """ | 821 | """ |
| 704 | Set background regions for an observation | 822 | Set background regions for an observation |
| 705 | 823 | ||
| @@ -716,7 +834,7 @@ class csphagen(ctools.csobservation): | @@ -716,7 +834,7 @@ class csphagen(ctools.csobservation): | ||
| 716 | # Initialise empty background regions for this observation | 834 | # Initialise empty background regions for this observation |
| 717 | bkg_reg = gammalib.GSkyRegions() | 835 | bkg_reg = gammalib.GSkyRegions() |
| 718 | 836 | ||
| 719 | - # If reflected background is requested then created reflected | 837 | + # If reflected background is requested then create reflected |
| 720 | # background regions | 838 | # background regions |
| 721 | if self['bkgmethod'].string() == 'REFLECTED': | 839 | if self['bkgmethod'].string() == 'REFLECTED': |
| 722 | bkg_reg = self._reflected_regions(obs) | 840 | bkg_reg = self._reflected_regions(obs) |
| @@ -728,6 +846,11 @@ class csphagen(ctools.csobservation): | @@ -728,6 +846,11 @@ class csphagen(ctools.csobservation): | ||
| 728 | elif self['bkgmethod'].string() == 'CUSTOM': | 846 | elif self['bkgmethod'].string() == 'CUSTOM': |
| 729 | bkg_reg = obs.off_regions().copy() | 847 | bkg_reg = obs.off_regions().copy() |
| 730 | 848 | ||
| 849 | + # ... otherwise if dedicated Off runs are use then | ||
| 850 | + # use background region that correspond to the same instrument coordinates | ||
| 851 | + if self['bkgmethod'].string() == 'OFF': | ||
| 852 | + bkg_reg = self._instrument_regions(obs,obs_off) | ||
| 853 | + | ||
| 731 | # Return background regions | 854 | # Return background regions |
| 732 | return bkg_reg | 855 | return bkg_reg |
| 733 | 856 | ||
| @@ -750,6 +873,13 @@ class csphagen(ctools.csobservation): | @@ -750,6 +873,13 @@ class csphagen(ctools.csobservation): | ||
| 750 | bkg_reg = None | 873 | bkg_reg = None |
| 751 | obs = self.obs()[i] | 874 | obs = self.obs()[i] |
| 752 | 875 | ||
| 876 | + # Retrieve dedicated Off observation if it exists | ||
| 877 | + if not self.obs_off().is_empty(): | ||
| 878 | + obs_off = self.obs_off()[i] | ||
| 879 | + # Otherwise use the same as On | ||
| 880 | + else: | ||
| 881 | + obs_off = self.obs()[i] | ||
| 882 | + | ||
| 753 | # Log header | 883 | # Log header |
| 754 | self._log_header3(gammalib.NORMAL,'%s observation "%s"' % \ | 884 | self._log_header3(gammalib.NORMAL,'%s observation "%s"' % \ |
| 755 | (obs.instrument(), obs.id())) | 885 | (obs.instrument(), obs.id())) |
| @@ -768,32 +898,43 @@ class csphagen(ctools.csobservation): | @@ -768,32 +898,43 @@ class csphagen(ctools.csobservation): | ||
| 768 | msg = ' Background model not used, assume constant backround rate.' | 898 | msg = ' Background model not used, assume constant backround rate.' |
| 769 | self._log_string(gammalib.NORMAL, msg) | 899 | self._log_string(gammalib.NORMAL, msg) |
| 770 | 900 | ||
| 771 | - # Set background regions for this observation | ||
| 772 | - bkg_reg = self._set_background_regions(obs) | 901 | + # Get offset angle of source |
| 902 | + pnt_dir = obs.pointing().dir() | ||
| 903 | + offset = pnt_dir.dist_deg(self._src_dir) | ||
| 773 | 904 | ||
| 774 | - # If there are background regions then create On/Off observation | ||
| 775 | - # and append it to the output container | ||
| 776 | - if bkg_reg.size() >= self['bkgregmin'].integer(): | 905 | + # Skip observation if it is pointed too far from the source |
| 906 | + if offset >= self['maxoffset'].real(): | ||
| 907 | + msg = ' Skip because observation is pointed at %.3f deg >= ' \ | ||
| 908 | + '"maxoffset=%.3f" from source.' \ | ||
| 909 | + % (offset, self['maxoffset'].real()) | ||
| 910 | + self._log_string(gammalib.NORMAL, msg) | ||
| 911 | + # ... otherwise continue to process | ||
| 912 | + else: | ||
| 777 | 913 | ||
| 778 | - # Create On/Off observation | ||
| 779 | - onoff = gammalib.GCTAOnOffObservation(obs, | ||
| 780 | - self._models, | ||
| 781 | - self._srcname, | ||
| 782 | - self._etruebounds, | ||
| 783 | - self._ebounds, | ||
| 784 | - self._src_reg, | ||
| 785 | - bkg_reg, | ||
| 786 | - use_model_bkg) | 914 | + # Set background regions for this observation |
| 915 | + bkg_reg = self._set_background_regions(obs,obs_off) | ||
| 787 | 916 | ||
| 788 | - # Set On/Off observation ID | ||
| 789 | - onoff.id(obs.id()) | 917 | + # If there are any background regions then create On/Off observation |
| 918 | + # and append it to the output container | ||
| 919 | + if bkg_reg.size() >= 0: | ||
| 790 | 920 | ||
| 791 | - # Otherwise log observation skipped | ||
| 792 | - else: | ||
| 793 | - msg = ' Skip because the number %d of regions ' \ | ||
| 794 | - 'for background estimation is smaller than ' \ | ||
| 795 | - '"bkgregmin"=%d.' % (bkg_reg.size(), self['bkgregmin'].integer()) | ||
| 796 | - self._log_string(gammalib.NORMAL, msg) | 921 | + # Create On/Off observation |
| 922 | + onoff = gammalib.GCTAOnOffObservation(obs, obs_off, | ||
| 923 | + self._models, | ||
| 924 | + self._srcname, | ||
| 925 | + self._etruebounds, | ||
| 926 | + self._ebounds, | ||
| 927 | + self._src_reg, | ||
| 928 | + bkg_reg, | ||
| 929 | + use_model_bkg) | ||
| 930 | + | ||
| 931 | + # Set On/Off observation ID | ||
| 932 | + onoff.id(obs.id()) | ||
| 933 | + | ||
| 934 | + # Otherwise log observation skipped | ||
| 935 | + else: | ||
| 936 | + msg = ' Skip because no valid Off regions could be determined' | ||
| 937 | + self._log_string(gammalib.NORMAL, msg) | ||
| 797 | 938 | ||
| 798 | # Construct dictionary with results | 939 | # Construct dictionary with results |
| 799 | result = {'onoff' : onoff, | 940 | result = {'onoff' : onoff, |
| @@ -854,6 +995,8 @@ class csphagen(ctools.csobservation): | @@ -854,6 +995,8 @@ class csphagen(ctools.csobservation): | ||
| 854 | 995 | ||
| 855 | # Write observation into logger | 996 | # Write observation into logger |
| 856 | self._log_observations(gammalib.NORMAL, self.obs(), 'Observation') | 997 | self._log_observations(gammalib.NORMAL, self.obs(), 'Observation') |
| 998 | + if not self.obs_off().is_empty(): | ||
| 999 | + self._log_observations(gammalib.NORMAL, self._obs_off, 'Off Observation') | ||
| 857 | 1000 | ||
| 858 | # Set true energy bins | 1001 | # Set true energy bins |
| 859 | self._etruebounds = self._etrue_ebounds() | 1002 | self._etruebounds = self._etrue_ebounds() |
| @@ -1032,6 +1175,27 @@ class csphagen(ctools.csobservation): | @@ -1032,6 +1175,27 @@ class csphagen(ctools.csobservation): | ||
| 1032 | # Return | 1175 | # Return |
| 1033 | return self._excl_reg | 1176 | return self._excl_reg |
| 1034 | 1177 | ||
| 1178 | + def obs_off(self, obs=None): | ||
| 1179 | + """ | ||
| 1180 | + Return and optionally set the Off observations | ||
| 1181 | + | ||
| 1182 | + Parameters | ||
| 1183 | + ---------- | ||
| 1184 | + obs : `~gammalib.GCTAObservations` | ||
| 1185 | + Off observations container | ||
| 1186 | + | ||
| 1187 | + Returns | ||
| 1188 | + ------- | ||
| 1189 | + observation container : `~gammalib.GCTAObservations` | ||
| 1190 | + Off observations container | ||
| 1191 | + """ | ||
| 1192 | + # If an observation container is provided then set the Off observations ... | ||
| 1193 | + if obs is not None: | ||
| 1194 | + self._obs_off = obs | ||
| 1195 | + | ||
| 1196 | + # Return | ||
| 1197 | + return self._obs_off | ||
| 1198 | + | ||
| 1035 | 1199 | ||
| 1036 | # ======================== # | 1200 | # ======================== # |
| 1037 | # Main routine entry point # | 1201 | # Main routine entry point # |
doc/source/users/reference_manual/csphagen.rst
| @@ -17,17 +17,16 @@ the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF), | @@ -17,17 +17,16 @@ the OGIP format normally used in X-ray astronomy (PHA, ARF, RMF), | ||
| 17 | If the ``bkgmethod`` parameter is ``REFLECTED``, the script will use an On region | 17 | If the ``bkgmethod`` parameter is ``REFLECTED``, the script will use an On region |
| 18 | shape defined by the ``srcshape`` parameter and generate reflected regions by | 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 | 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 | 20 | +pointing direction within the same observation. So far circular and |
| 21 | +rectangular regions are supported. Circular regions are defined by a | ||
| 22 | +center that is specified by the ``ra`` and ``dec`` parameters if | ||
| 23 | +``coordsys=CEL`` or the ``glon`` and ``glat`` parameters if | ||
| 23 | ``coordsys=GAL``, and a radius that is specified by the ``rad`` parameter. | 24 | ``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 | 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 | ``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 | a ``height`` parameter and a ``posang`` parameter that specified the position |
| 27 | angle measured counter clockwise from celestial North. Note that the ``height`` | 28 | 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 | 29 | +axis is pointing in the direction of the position angle. If less than |
| 31 | ``bkgregmin`` reflected regions can be found because the offset angle between | 30 | ``bkgregmin`` reflected regions can be found because the offset angle between |
| 32 | source and pointing direction is too small, the observation will be skipped. | 31 | 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 | 32 | Regions to be excluded for background determination can be specified using a |
| @@ -44,7 +43,7 @@ respectively. Both parameters accept either a ds9 region file or a sky map where | @@ -44,7 +43,7 @@ respectively. Both parameters accept either a ds9 region file or a sky map where | ||
| 44 | all non-zero pixel values define a region. In case that an observation definition | 43 | all non-zero pixel values define a region. In case that an observation definition |
| 45 | XML file is specified as input using the ``inobs``, background region files can | 44 | XML file is specified as input using the ``inobs``, background region files can |
| 46 | be specified separately for each observation using the ``OffRegions`` parameter | 45 | be specified separately for each observation using the ``OffRegions`` parameter |
| 47 | -in the format | 46 | +in the format. |
| 48 | 47 | ||
| 49 | .. code-block:: xml | 48 | .. code-block:: xml |
| 50 | 49 | ||
| @@ -53,6 +52,16 @@ in the format | @@ -53,6 +52,16 @@ in the format | ||
| 53 | ... | 52 | ... |
| 54 | </observation> | 53 | </observation> |
| 55 | 54 | ||
| 55 | +Finally, if ``bkgmethod=OFF`` is specified, the background regions | ||
| 56 | +will be taken from independent Off observations given by the | ||
| 57 | +``inobsoff`` parameter. The On region can be defined in the same way | ||
| 58 | +as for ``bkgmethod=REFLECTED``, and the Off regions are the | ||
| 59 | +corresponding regions in instrument coordinates from the Off observations. | ||
| 60 | + | ||
| 61 | +For all background estimation methods only observations with a | ||
| 62 | +pointing direction at an offset smaller than ``maxoffset`` from the | ||
| 63 | +source are used. | ||
| 64 | + | ||
| 56 | :ref:`csphagen` supports multiprocessing. By default each event list will be treated | 65 | :ref:`csphagen` supports multiprocessing. By default each event list will be treated |
| 57 | in parallel over as many processes as the number of CPUs available on your | 66 | in parallel over as many processes as the number of CPUs available on your |
| 58 | machine. The maximum number of parallel processes can be set by the user through the | 67 | machine. The maximum number of parallel processes can be set by the user through the |
| @@ -78,6 +87,11 @@ General parameters | @@ -78,6 +87,11 @@ General parameters | ||
| 78 | ``inobs [file]`` | 87 | ``inobs [file]`` |
| 79 | Input event list or observation definition XML file. | 88 | Input event list or observation definition XML file. |
| 80 | 89 | ||
| 90 | +``inobsoff [file]`` | ||
| 91 | + Input event list or observation definition XML file for dedicated | ||
| 92 | + Off observations. The number of observations must be the same as | ||
| 93 | + for ``inobs``. | ||
| 94 | + | ||
| 81 | ``inmodel [file]`` | 95 | ``inmodel [file]`` |
| 82 | Input model XML file (if ``NONE`` a point source at the centre of the | 96 | Input model XML file (if ``NONE`` a point source at the centre of the |
| 83 | source region is used). | 97 | source region is used). |
| @@ -172,7 +186,7 @@ General parameters | @@ -172,7 +186,7 @@ General parameters | ||
| 172 | ``srcregfile [file]`` | 186 | ``srcregfile [file]`` |
| 173 | Source region file (ds9 or FITS WCS map). | 187 | Source region file (ds9 or FITS WCS map). |
| 174 | 188 | ||
| 175 | -``bkgmethod <REFLECTED|CUSTOM> [string]`` | 189 | +``bkgmethod <REFLECTED|CUSTOM|OFF> [string]`` |
| 176 | Method for background estimation: | 190 | Method for background estimation: |
| 177 | 191 | ||
| 178 | - ``REFLECTED``: background is evaluated in regions with the same shape as | 192 | - ``REFLECTED``: background is evaluated in regions with the same shape as |
| @@ -184,18 +198,25 @@ General parameters | @@ -184,18 +198,25 @@ General parameters | ||
| 184 | observations specified in the observation definition XML file the name of | 198 | observations specified in the observation definition XML file the name of |
| 185 | the region file will be extracted from the ``OffRegions`` parameter that | 199 | the region file will be extracted from the ``OffRegions`` parameter that |
| 186 | needs to be specified for each observation in the observation definition | 200 | needs to be specified for each observation in the observation definition |
| 187 | - XML file. Off region files can be either ds9 region files or FITS WCS maps. | 201 | + XML file. Off region files can be either ds9 region files or |
| 202 | + FITS WCS maps. | ||
| 203 | + | ||
| 204 | + - ``OFF``: background is evaluated from dedicated Off | ||
| 205 | + observations in regions that correspond to the source region | ||
| 206 | + in the On observations in instrument coordinates. | ||
| 188 | 207 | ||
| 189 | ``bkgregfile [file]`` | 208 | ``bkgregfile [file]`` |
| 190 | Background regions file (ds9 or FITS WCS map). | 209 | Background regions file (ds9 or FITS WCS map). |
| 191 | 210 | ||
| 192 | ``(bkgregmin = 2) [integer]`` | 211 | ``(bkgregmin = 2) [integer]`` |
| 193 | - Minimum number of background regions that are required for an observation. | 212 | + Minimum number of background regions that are required for an |
| 213 | + observation when the ``REFLECTED`` background estimation is used. | ||
| 194 | If this number of background regions is not available the observation is | 214 | If this number of background regions is not available the observation is |
| 195 | skipped. | 215 | skipped. |
| 196 | 216 | ||
| 197 | ``(bkgregskip = 1) [integer]`` | 217 | ``(bkgregskip = 1) [integer]`` |
| 198 | - Number of background regions that should be skipped next to the On regions. | 218 | + Number of background regions that should be skipped next to the On |
| 219 | + region when the ``REFLECTED`` background estimation is used. | ||
| 199 | Typically, one region is skipped so that the Off regions are taken sufficiently | 220 | Typically, one region is skipped so that the Off regions are taken sufficiently |
| 200 | distant from the On region, but in some cases it may be useful to keep the | 221 | distant from the On region, but in some cases it may be useful to keep the |
| 201 | background regions next to the On region. | 222 | background regions next to the On region. |
test/data/background_events.fits
0 โ 100644
No preview for this file type
test/test_csphagen.py
| @@ -46,6 +46,7 @@ class Test(test): | @@ -46,6 +46,7 @@ class Test(test): | ||
| 46 | self._myevents1 = self._datadir + '/crab_offaxis1.fits' | 46 | self._myevents1 = self._datadir + '/crab_offaxis1.fits' |
| 47 | self._myevents2 = self._datadir + '/crab_offaxis2.fits' | 47 | self._myevents2 = self._datadir + '/crab_offaxis2.fits' |
| 48 | self._exclusion = self._datadir + '/crab_exclusion.fits' | 48 | self._exclusion = self._datadir + '/crab_exclusion.fits' |
| 49 | + self._offevents = self._datadir + '/background_events.fits' | ||
| 49 | self._nreg_with_excl = 5 | 50 | self._nreg_with_excl = 5 |
| 50 | self._nreg_wo_excl = 8 | 51 | self._nreg_wo_excl = 8 |
| 51 | self._nreg_bkg_reg = 5 | 52 | self._nreg_bkg_reg = 5 |
| @@ -414,6 +415,101 @@ class Test(test): | @@ -414,6 +415,101 @@ class Test(test): | ||
| 414 | self._check_output('csphagen_py8', self._nbins, 9) | 415 | self._check_output('csphagen_py8', self._nbins, 9) |
| 415 | self._check_outobs('csphagen_py8_obs.xml', 1) | 416 | self._check_outobs('csphagen_py8_obs.xml', 1) |
| 416 | 417 | ||
| 418 | + # Setup csphagen for test with dedicated Off observations | ||
| 419 | + phagen = cscripts.csphagen() | ||
| 420 | + phagen['inobs'] = self._myevents1 | ||
| 421 | + phagen['inobsoff'] = self._offevents | ||
| 422 | + phagen['inmodel'] = 'NONE' | ||
| 423 | + phagen['caldb'] = self._caldb | ||
| 424 | + phagen['irf'] = self._irf | ||
| 425 | + phagen['ebinalg'] = 'LOG' | ||
| 426 | + phagen['emin'] = 0.1 | ||
| 427 | + phagen['emax'] = 100.0 | ||
| 428 | + phagen['enumbins'] = self._nbins | ||
| 429 | + phagen['bkgmethod'] = 'OFF' | ||
| 430 | + phagen['etruemin'] = 0.05 | ||
| 431 | + phagen['etruemax'] = 150.0 | ||
| 432 | + phagen['etruebins'] = 5 | ||
| 433 | + phagen['stack'] = False | ||
| 434 | + phagen['outobs'] = 'csphagen_py9_obs.xml' | ||
| 435 | + phagen['outmodel'] = 'csphagen_py9_model.xml' | ||
| 436 | + phagen['prefix'] = 'csphagen_py9' | ||
| 437 | + phagen['logfile'] = 'csphagen_py9.log' | ||
| 438 | + phagen['chatter'] = 2 | ||
| 439 | + | ||
| 440 | + # Execute script | ||
| 441 | + phagen.execute() | ||
| 442 | + | ||
| 443 | + # Check output | ||
| 444 | + self._check_output('csphagen_py9', self._nbins, 1) | ||
| 445 | + self._check_outobs('csphagen_py9_obs.xml', 1) | ||
| 446 | + | ||
| 447 | + # Setup csphagen for test with dedicated Off observations | ||
| 448 | + # and exclusion region | ||
| 449 | + phagen = cscripts.csphagen() | ||
| 450 | + phagen['inobs'] = self._myevents1 | ||
| 451 | + phagen['inobsoff'] = self._offevents | ||
| 452 | + phagen['inmodel'] = 'NONE' | ||
| 453 | + phagen['caldb'] = self._caldb | ||
| 454 | + phagen['irf'] = self._irf | ||
| 455 | + phagen['ebinalg'] = 'LOG' | ||
| 456 | + phagen['emin'] = 0.1 | ||
| 457 | + phagen['emax'] = 100.0 | ||
| 458 | + phagen['enumbins'] = self._nbins | ||
| 459 | + phagen['inexclusion'] = self._exclusion | ||
| 460 | + phagen['bkgmethod'] = 'OFF' | ||
| 461 | + phagen['etruemin'] = 0.05 | ||
| 462 | + phagen['etruemax'] = 150.0 | ||
| 463 | + phagen['etruebins'] = 5 | ||
| 464 | + phagen['stack'] = False | ||
| 465 | + phagen['outobs'] = 'csphagen_py10_obs.xml' | ||
| 466 | + phagen['outmodel'] = 'csphagen_py10_model.xml' | ||
| 467 | + phagen['prefix'] = 'csphagen_py10' | ||
| 468 | + phagen['logfile'] = 'csphagen_py10.log' | ||
| 469 | + phagen['chatter'] = 2 | ||
| 470 | + | ||
| 471 | + # Execute script | ||
| 472 | + phagen.execute() | ||
| 473 | + | ||
| 474 | + # Check output | ||
| 475 | + self._check_output('csphagen_py10', self._nbins, 1) | ||
| 476 | + self._check_outobs('csphagen_py10_obs.xml', 1) | ||
| 477 | + | ||
| 478 | + # Test with multiple input observations and dedicated Off runs | ||
| 479 | + # Create Off observation container | ||
| 480 | + obs_off = gammalib.GObservations() | ||
| 481 | + for s, events in enumerate([self._offevents, self._offevents]): | ||
| 482 | + run = gammalib.GCTAObservation(events) | ||
| 483 | + run.id(str(s + 1)) | ||
| 484 | + run.response(self._irf, gammalib.GCaldb('cta', self._caldb)) | ||
| 485 | + obs_off.append(run) | ||
| 486 | + | ||
| 487 | + # Setup csphagen | ||
| 488 | + phagen = cscripts.csphagen(obs) | ||
| 489 | + phagen.obs_off(obs_off) | ||
| 490 | + phagen['inmodel'] = 'NONE' | ||
| 491 | + phagen['ebinalg'] = 'LOG' | ||
| 492 | + phagen['emin'] = 0.1 | ||
| 493 | + phagen['emax'] = 100.0 | ||
| 494 | + phagen['enumbins'] = self._nbins | ||
| 495 | + phagen['bkgmethod'] = 'OFF' | ||
| 496 | + phagen['etruemin'] = 0.05 | ||
| 497 | + phagen['etruemax'] = 150.0 | ||
| 498 | + phagen['etruebins'] = 5 | ||
| 499 | + phagen['stack'] = False | ||
| 500 | + phagen['outobs'] = 'csphagen_py11_obs.xml' | ||
| 501 | + phagen['outmodel'] = 'csphagen_py11_model.xml' | ||
| 502 | + phagen['prefix'] = 'csphagen_py11' | ||
| 503 | + phagen['logfile'] = 'csphagen_py11.log' | ||
| 504 | + phagen['chatter'] = 2 | ||
| 505 | + | ||
| 506 | + # Execute script | ||
| 507 | + phagen.execute() | ||
| 508 | + | ||
| 509 | + # Check output | ||
| 510 | + self._check_output('csphagen_py11', self._nbins, 1) | ||
| 511 | + self._check_outobs('csphagen_py11_obs.xml', 1) | ||
| 512 | + | ||
| 417 | # Return | 513 | # Return |
| 418 | return | 514 | return |
| 419 | 515 |