Commit 679001d775501df8c13b97fbe4bc9d38c89ad697
1 parent
9d09558a
ctbkgcube now requires an input counts cube (#2763)
Showing
10 changed files
with
118 additions
and
310 deletions
ChangeLog
| 1 | -2018-11-16 | |
| 1 | +2018-11-22 | |
| 2 | 2 | |
| 3 | 3 | * Version 1.6.0 released |
| 4 | 4 | ======================== |
| 5 | 5 | |
| 6 | + ctbkgcube now requires an input counts cube (#2763) | |
| 6 | 7 | CTA background rates are now per ontime (#2758) |
| 7 | 8 | Apply event bin weighting in ctmodel for stacked computation (#2721) |
| 8 | 9 | Enable ctbin to produce counts cubes for a joint binned analysis (#2714) | ... | ... |
NEWS
| 1 | 1 | New Features and Important Changes in ctools 1.6.0 |
| 2 | 2 | |
| 3 | -16 November 2018 | |
| 3 | +22 November 2018 | |
| 4 | 4 | |
| 5 | 5 | |
| 6 | 6 | Introduction |
| ... | ... | @@ -83,6 +83,11 @@ Copy the time reference from the first input observation (#2054). |
| 83 | 83 | |
| 84 | 84 | ctbkgcube - CTA background cube generation |
| 85 | 85 | ------------------------------------------ |
| 86 | +The tool now required a counts cube on input, either specified via the | |
| 87 | +ctbkgcube::cntcube() method or the "incube" parameter. This assures that the | |
| 88 | +background cube has the same dimension as the counts cube. All spatial and | |
| 89 | +spectral binning parameters were removed from the parameter file (#2763). | |
| 90 | + | |
| 86 | 91 | Add ctbkgcube::cntcube() method so that a counts cube be specified as |
| 87 | 92 | instance before running the tool. This avoids storing the counts cube on |
| 88 | 93 | disk if the counts cube definition should be used to generate a background | ... | ... |
cscripts/obsutils.py
| ... | ... | @@ -185,14 +185,10 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, |
| 185 | 185 | # proper background cube. |
| 186 | 186 | cntcube = binning.cube() |
| 187 | 187 | |
| 188 | - # Get stacked response (use pointing for map centre) | |
| 189 | - response = get_stacked_response(obssim.obs(), None, None, | |
| 190 | - binsz=binsz, nxpix=npix, nypix=npix, | |
| 191 | - emin=emin, emax=emax, enumbins=nbins, | |
| 188 | + # Get stacked response | |
| 189 | + response = get_stacked_response(obssim.obs(), cntcube, | |
| 192 | 190 | edisp=edisp, |
| 193 | - coordsys=coord, proj=proj, | |
| 194 | 191 | addbounds=addbounds, |
| 195 | - cntcube=cntcube, | |
| 196 | 192 | log=log, debug=debug, |
| 197 | 193 | chatter=chatter) |
| 198 | 194 | |
| ... | ... | @@ -519,50 +515,24 @@ def set_observations(ra, dec, rad, tstart, duration, emin, emax, irf, caldb, |
| 519 | 515 | # ==================== # |
| 520 | 516 | # Get stacked response # |
| 521 | 517 | # ==================== # |
| 522 | -def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, | |
| 523 | - emin=0.1, emax=100.0, enumbins=20, edisp=False, | |
| 524 | - coordsys='GAL', proj='TAN', addbounds=False, | |
| 525 | - cntcube=None, log=False, debug=False, chatter=2): | |
| 518 | +def get_stacked_response(obs, cntcube, edisp=False, addbounds=False, | |
| 519 | + log=False, debug=False, chatter=2): | |
| 526 | 520 | """ |
| 527 | 521 | Get stacked response cubes |
| 528 | 522 | |
| 529 | - The number of energies bins are set to at least 30 bins per decade, and | |
| 530 | - the "enumbins" parameter is only used if the number of bins is larger | |
| 531 | - than 30 bins per decade. | |
| 532 | - | |
| 533 | - If the xref or yref arguments are "None" the response cube centre will be | |
| 534 | - determined from the pointing information in the observation container. | |
| 523 | + The number of energies bins are set to at least 30 bins per decade, unless | |
| 524 | + the counts cube has more energy bins per decade. | |
| 535 | 525 | |
| 536 | 526 | Parameters |
| 537 | 527 | ---------- |
| 538 | 528 | obs : `~gammalib.GObservations` |
| 539 | 529 | Observation container |
| 540 | - xref : float | |
| 541 | - Right Ascension or Galactic longitude of response centre (deg) | |
| 542 | - yref : float | |
| 543 | - Declination or Galactic latitude of response centre (deg) | |
| 544 | - binsz : float, optional | |
| 545 | - Pixel size (deg/pixel) | |
| 546 | - nxpix : int, optional | |
| 547 | - Number of pixels in X direction | |
| 548 | - nypix : int, optional | |
| 549 | - Number of pixels in Y direction | |
| 550 | - emin : float, optional | |
| 551 | - Minimum energy (TeV) | |
| 552 | - emax : float, optional | |
| 553 | - Maximum energy (TeV) | |
| 554 | - enumbins : int, optional | |
| 555 | - Number of energy bins | |
| 530 | + cntcube : `~gammalib.GCTAEventCube` | |
| 531 | + Counts cube | |
| 556 | 532 | edisp : bool, optional |
| 557 | 533 | Apply energy dispersion? |
| 558 | - coordsys : str, optional | |
| 559 | - Coordinate system | |
| 560 | - proj : str, optional | |
| 561 | - Projection | |
| 562 | 534 | addbounds : bool, optional |
| 563 | 535 | Add boundaries at observation energies |
| 564 | - cntcube : `~gammalib.GCTAEventCube`, optional | |
| 565 | - Counts cube | |
| 566 | 536 | log : bool, optional |
| 567 | 537 | Create log file(s) |
| 568 | 538 | debug : bool, optional |
| ... | ... | @@ -575,26 +545,32 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, |
| 575 | 545 | result : dict |
| 576 | 546 | Dictionary of response cubes |
| 577 | 547 | """ |
| 578 | - # If no xref and yref arguments have been specified then use the pointing | |
| 579 | - # information | |
| 580 | - if xref == None or yref == None: | |
| 581 | - usepnt = True | |
| 582 | - else: | |
| 583 | - usepnt = False | |
| 548 | + # Derive binning parameters from counts cube | |
| 549 | + xref = cntcube.counts().projection().crval(0) | |
| 550 | + yref = cntcube.counts().projection().crval(1) | |
| 551 | + binsz = cntcube.counts().projection().cdelt(1) | |
| 552 | + coordsys = cntcube.counts().projection().coordsys() | |
| 553 | + proj = cntcube.counts().projection().code() | |
| 554 | + nxpix = cntcube.nx() | |
| 555 | + nypix = cntcube.ny() | |
| 556 | + emin = cntcube.emin().TeV() | |
| 557 | + emax = cntcube.emax().TeV() | |
| 558 | + ebins = cntcube.ebins() | |
| 559 | + | |
| 560 | + # Translate coordinate system | |
| 561 | + if coordsys == 'EQU': | |
| 562 | + coordsys = 'CEL' | |
| 584 | 563 | |
| 585 | 564 | # Set energy limits, with larger etrue limits for the case that energy |
| 586 | 565 | # dispersion is requested |
| 587 | 566 | if edisp: |
| 588 | - _emin = 0.5 * emin | |
| 589 | - _emax = 1.5 * emax | |
| 590 | - else: | |
| 591 | - _emin = emin | |
| 592 | - _emax = emax | |
| 567 | + emin *= 0.5 | |
| 568 | + emax *= 1.5 | |
| 593 | 569 | |
| 594 | 570 | # Set number of energy bins to at least 30 per energy decade |
| 595 | - _enumbins = int((math.log10(_emax) - math.log10(_emin)) * 30.0) | |
| 596 | - if enumbins > _enumbins: | |
| 597 | - _enumbins = enumbins | |
| 571 | + enumbins = int((math.log10(emax) - math.log10(emin)) * 30.0) | |
| 572 | + if ebins > enumbins: | |
| 573 | + enumbins = ebins | |
| 598 | 574 | |
| 599 | 575 | # Compute spatial binning for point spread function and energy dispersion |
| 600 | 576 | # cubes. The spatial binning is 10 times coarser than the spatial binning |
| ... | ... | @@ -606,22 +582,20 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, |
| 606 | 582 | # Create exposure cube |
| 607 | 583 | expcube = ctools.ctexpcube(obs) |
| 608 | 584 | expcube['incube'] = 'NONE' |
| 609 | - expcube['usepnt'] = usepnt | |
| 610 | 585 | expcube['ebinalg'] = 'LOG' |
| 586 | + expcube['xref'] = xref | |
| 587 | + expcube['yref'] = yref | |
| 611 | 588 | expcube['binsz'] = binsz |
| 612 | 589 | expcube['nxpix'] = nxpix |
| 613 | 590 | expcube['nypix'] = nypix |
| 614 | - expcube['enumbins'] = _enumbins | |
| 615 | - expcube['emin'] = _emin | |
| 616 | - expcube['emax'] = _emax | |
| 591 | + expcube['enumbins'] = enumbins | |
| 592 | + expcube['emin'] = emin | |
| 593 | + expcube['emax'] = emax | |
| 617 | 594 | expcube['coordsys'] = coordsys |
| 618 | 595 | expcube['proj'] = proj |
| 619 | 596 | expcube['addbounds'] = addbounds |
| 620 | 597 | expcube['debug'] = debug |
| 621 | 598 | expcube['chatter'] = chatter |
| 622 | - if not usepnt: | |
| 623 | - expcube['xref'] = xref | |
| 624 | - expcube['yref'] = yref | |
| 625 | 599 | if log: |
| 626 | 600 | expcube.logFileOpen() |
| 627 | 601 | expcube.run() |
| ... | ... | @@ -629,22 +603,20 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, |
| 629 | 603 | # Create point spread function cube |
| 630 | 604 | psfcube = ctools.ctpsfcube(obs) |
| 631 | 605 | psfcube['incube'] = 'NONE' |
| 632 | - psfcube['usepnt'] = usepnt | |
| 633 | 606 | psfcube['ebinalg'] = 'LOG' |
| 607 | + psfcube['xref'] = xref | |
| 608 | + psfcube['yref'] = yref | |
| 634 | 609 | psfcube['binsz'] = psf_binsz |
| 635 | 610 | psfcube['nxpix'] = psf_nxpix |
| 636 | 611 | psfcube['nypix'] = psf_nypix |
| 637 | - psfcube['enumbins'] = _enumbins | |
| 638 | - psfcube['emin'] = _emin | |
| 639 | - psfcube['emax'] = _emax | |
| 612 | + psfcube['enumbins'] = enumbins | |
| 613 | + psfcube['emin'] = emin | |
| 614 | + psfcube['emax'] = emax | |
| 640 | 615 | psfcube['coordsys'] = coordsys |
| 641 | 616 | psfcube['proj'] = proj |
| 642 | 617 | psfcube['addbounds'] = addbounds |
| 643 | 618 | psfcube['debug'] = debug |
| 644 | 619 | psfcube['chatter'] = chatter |
| 645 | - if not usepnt: | |
| 646 | - psfcube['xref'] = xref | |
| 647 | - psfcube['yref'] = yref | |
| 648 | 620 | if log: |
| 649 | 621 | psfcube.logFileOpen() |
| 650 | 622 | psfcube.run() |
| ... | ... | @@ -653,26 +625,9 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, |
| 653 | 625 | # for the counts cube since no interpolation should be actually done |
| 654 | 626 | # for the background cube. |
| 655 | 627 | bkgcube = ctools.ctbkgcube(obs) |
| 656 | - if cntcube != None: | |
| 657 | - bkgcube.cntcube(cntcube) | |
| 658 | - else: | |
| 659 | - bkgcube['incube'] = 'NONE' | |
| 660 | - bkgcube['ebinalg'] = 'LOG' | |
| 661 | - bkgcube['binsz'] = binsz | |
| 662 | - bkgcube['nxpix'] = nxpix | |
| 663 | - bkgcube['nypix'] = nypix | |
| 664 | - bkgcube['enumbins'] = enumbins | |
| 665 | - bkgcube['emin'] = emin | |
| 666 | - bkgcube['emax'] = emax | |
| 667 | - bkgcube['coordsys'] = coordsys | |
| 668 | - bkgcube['proj'] = proj | |
| 669 | - if not usepnt: | |
| 670 | - bkgcube['xref'] = xref | |
| 671 | - bkgcube['yref'] = yref | |
| 672 | - bkgcube['usepnt'] = usepnt | |
| 673 | - bkgcube['addbounds'] = addbounds | |
| 674 | - bkgcube['debug'] = debug | |
| 675 | - bkgcube['chatter'] = chatter | |
| 628 | + bkgcube.cntcube(cntcube) | |
| 629 | + bkgcube['debug'] = debug | |
| 630 | + bkgcube['chatter'] = chatter | |
| 676 | 631 | if log: |
| 677 | 632 | bkgcube.logFileOpen() |
| 678 | 633 | bkgcube.run() |
| ... | ... | @@ -681,22 +636,20 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, |
| 681 | 636 | if edisp: |
| 682 | 637 | edispcube = ctools.ctedispcube(obs) |
| 683 | 638 | edispcube['incube'] = 'NONE' |
| 684 | - edispcube['usepnt'] = usepnt | |
| 685 | 639 | edispcube['ebinalg'] = 'LOG' |
| 640 | + edispcube['xref'] = xref | |
| 641 | + edispcube['yref'] = yref | |
| 686 | 642 | edispcube['binsz'] = psf_binsz |
| 687 | 643 | edispcube['nxpix'] = psf_nxpix |
| 688 | 644 | edispcube['nypix'] = psf_nypix |
| 689 | - edispcube['enumbins'] = _enumbins | |
| 690 | - edispcube['emin'] = _emin | |
| 691 | - edispcube['emax'] = _emax | |
| 645 | + edispcube['enumbins'] = enumbins | |
| 646 | + edispcube['emin'] = emin | |
| 647 | + edispcube['emax'] = emax | |
| 692 | 648 | edispcube['coordsys'] = coordsys |
| 693 | 649 | edispcube['proj'] = proj |
| 694 | 650 | edispcube['addbounds'] = addbounds |
| 695 | 651 | edispcube['debug'] = debug |
| 696 | 652 | edispcube['chatter'] = chatter |
| 697 | - if not usepnt: | |
| 698 | - edispcube['xref'] = xref | |
| 699 | - edispcube['yref'] = yref | |
| 700 | 653 | if log: |
| 701 | 654 | edispcube.logFileOpen() |
| 702 | 655 | edispcube.run() |
| ... | ... | @@ -756,26 +709,14 @@ def get_stacked_obs(cls, obs): |
| 756 | 709 | |
| 757 | 710 | # Store counts cube so that we can use it to build the background |
| 758 | 711 | # cube |
| 759 | - cube = cntcube.cube().copy() | |
| 712 | + cube = cntcube.cube() | |
| 760 | 713 | |
| 761 | 714 | # Write header |
| 762 | 715 | if cls._logExplicit(): |
| 763 | 716 | cls._log.header3('Creating stacked response') |
| 764 | 717 | |
| 765 | 718 | # Get stacked response |
| 766 | - response = get_stacked_response(obs, | |
| 767 | - cls['xref'].real(), | |
| 768 | - cls['yref'].real(), | |
| 769 | - binsz=cls['binsz'].real(), | |
| 770 | - nxpix=cls['nxpix'].integer(), | |
| 771 | - nypix=cls['nypix'].integer(), | |
| 772 | - emin=cls['emin'].real(), | |
| 773 | - emax=cls['emax'].real(), | |
| 774 | - enumbins=cls['enumbins'].integer(), | |
| 775 | - edisp=cls['edisp'].boolean(), | |
| 776 | - coordsys=cls['coordsys'].string(), | |
| 777 | - proj=cls['proj'].string(), | |
| 778 | - cntcube=cube) | |
| 719 | + response = get_stacked_response(obs, cube, edisp=cls['edisp'].boolean()) | |
| 779 | 720 | |
| 780 | 721 | # Retrieve a new oberservation container |
| 781 | 722 | new_obs = cntcube.obs().copy() | ... | ... |
doc/source/admin/release_history/1.6.rst
| ... | ... | @@ -43,6 +43,8 @@ Bug fixes |
| 43 | 43 | Improvements |
| 44 | 44 | ------------ |
| 45 | 45 | |
| 46 | +* [`2763 <https://cta-redmine.irap.omp.eu/issues/2763>`_] - | |
| 47 | + :ref:`ctbkgcube` now requires an input counts cube | |
| 46 | 48 | * [`2714 <https://cta-redmine.irap.omp.eu/issues/2714>`_] - |
| 47 | 49 | Enable :ref:`ctbin` to produce counts cubes for a joint binned analysis. |
| 48 | 50 | Rename ``outcube`` to ``outobs`` parameter. | ... | ... |
doc/source/users/reference_manual/ctbkgcube.rst
| ... | ... | @@ -18,21 +18,14 @@ background counts in each background cube bin. |
| 18 | 18 | |
| 19 | 19 | :ref:`ctbkgcube` requires on input the event list or observation definition file |
| 20 | 20 | that has been used in the generation of the counts cube using :doc:`ctbin`. |
| 21 | - | |
| 22 | -To assure consistency between an existing counts cube and the | |
| 23 | -corresponding background cube, it is recommended to specify the counts | |
| 24 | -cube as ``incube`` parameter. This will instruct ctbkgcube to extract the | |
| 25 | -background cube definition (such as sky coordinates and projection, number | |
| 26 | -of pixels, pixel scale, energy binning) from the counts cube. | |
| 27 | - | |
| 28 | -Note that there should at least be 25 bins per energy decade in the background | |
| 29 | -cube to assure that the energy dependence of the background rate is sufficiently | |
| 30 | -well sampled (in particular at low energies). | |
| 21 | +To assure consistency between an existing counts cube and the | |
| 22 | +corresponding background cube, the counts cube needs also to be specified via | |
| 23 | +the ``incube`` parameter. | |
| 31 | 24 | |
| 32 | 25 | :ref:`ctbkgcube` generates a background cube FITS file comprising two extensions. |
| 33 | 26 | The primary extension contains a 3-dimensional image that contains the |
| 34 | -background cube values. The next extension named ``ENERGIES`` contains a | |
| 35 | -binary table that defines the energies of the background cube. | |
| 27 | +background cube values. The next extension named ``EBOUNDS`` contains a | |
| 28 | +binary table that defines the energy boundaries of the background cube. | |
| 36 | 29 | |
| 37 | 30 | :ref:`ctbkgcube` generates also an output model XML file that can serve as input |
| 38 | 31 | for a maximum likelihood analysis. The output model XML file is a copy of |
| ... | ... | @@ -53,12 +46,12 @@ General parameters |
| 53 | 46 | ``inobs [file]`` |
| 54 | 47 | Input event list or observation definition XML file. |
| 55 | 48 | |
| 56 | -``inmodel [file]`` | |
| 57 | - Input model XML file. | |
| 58 | - | |
| 59 | 49 | ``incube [file]`` |
| 60 | 50 | Input counts cube file to extract background cube definition. |
| 61 | 51 | |
| 52 | +``inmodel [file]`` | |
| 53 | + Input model XML file. | |
| 54 | + | |
| 62 | 55 | ``caldb [string]`` |
| 63 | 56 | Calibration database. |
| 64 | 57 | |
| ... | ... | @@ -71,55 +64,6 @@ General parameters |
| 71 | 64 | ``outmodel [file]`` |
| 72 | 65 | Output model XML file. |
| 73 | 66 | |
| 74 | -``ebinalg <FILE|LIN|LOG> [string]`` | |
| 75 | - Algorithm for defining energy bins. For ``FILE``, the energy bins are defined | |
| 76 | - in a FITS file that is specified by the ``ebinfile`` parameter, for ``LIN`` | |
| 77 | - and ``LOG`` there will be ``enumbins`` energy bins spaced linearly or | |
| 78 | - logarithmically between ``emin`` and ``emax``, respectively. | |
| 79 | - | |
| 80 | -``emin [real]`` | |
| 81 | - Lower energy value for first energy bin (in TeV) if ``LIN`` or ``LOG`` | |
| 82 | - energy binning algorithms are used. | |
| 83 | - | |
| 84 | -``emax [real]`` | |
| 85 | - Upper energy value for last energy bin (in TeV) if ``LIN`` or ``LOG`` | |
| 86 | - energy binning algorithms are used. | |
| 87 | - | |
| 88 | -``enumbins [integer]`` | |
| 89 | - Number of energy bins if ``LIN`` or ``LOG`` energy binning algorithms are | |
| 90 | - used. | |
| 91 | - | |
| 92 | -``ebinfile [file]`` | |
| 93 | - Name of the file containing the energy binning definition if ``ebinalg=FILE``. | |
| 94 | - You may use :ref:`csebins` to generate a file with appropriate energy binning. | |
| 95 | - | |
| 96 | -``(addbounds = no) [boolean]`` | |
| 97 | - Add energies to the background cube at the observation energy boundaries? | |
| 98 | - | |
| 99 | -``(usepnt = no) [boolean]`` | |
| 100 | - Use CTA pointing direction for cube centre instead of xref/yref parameters? | |
| 101 | - | |
| 102 | -``nxpix [integer]`` | |
| 103 | - Number of cube bins in Right Ascension or Galactic longitude. | |
| 104 | - | |
| 105 | -``nypix [integer]`` | |
| 106 | - Number of cube bins in Declination or Galactic latitude. | |
| 107 | - | |
| 108 | -``binsz [real]`` | |
| 109 | - Cube bin size (in degrees/pixel). | |
| 110 | - | |
| 111 | -``coordsysL <CEL|GAL> [string]`` | |
| 112 | - Coordinate system (CEL - celestial, GAL - galactic). | |
| 113 | - | |
| 114 | -``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]`` | |
| 115 | - Projection method. | |
| 116 | - | |
| 117 | -``xref [real]`` | |
| 118 | - Right Ascension / Galactic longitude of cube centre (J2000, in degrees). | |
| 119 | - | |
| 120 | -``yref [real]`` | |
| 121 | - Declination / Galactic latitude of cube centre (J2000, in degrees). | |
| 122 | - | |
| 123 | 67 | |
| 124 | 68 | Standard parameters |
| 125 | 69 | ------------------- | ... | ... |
src/ctbkgcube/ctbkgcube.cpp
| ... | ... | @@ -438,7 +438,6 @@ void ctbkgcube::init_members(void) |
| 438 | 438 | // Initialise user parameters |
| 439 | 439 | m_outcube.clear(); |
| 440 | 440 | m_outmodel.clear(); |
| 441 | - m_addbounds = true; | |
| 442 | 441 | m_publish = false; |
| 443 | 442 | m_chatter = static_cast<GChatter>(2); |
| 444 | 443 | |
| ... | ... | @@ -463,7 +462,6 @@ void ctbkgcube::copy_members(const ctbkgcube& app) |
| 463 | 462 | // Copy user parameters |
| 464 | 463 | m_outmodel = app.m_outmodel; |
| 465 | 464 | m_outcube = app.m_outcube; |
| 466 | - m_addbounds = app.m_addbounds; | |
| 467 | 465 | m_publish = app.m_publish; |
| 468 | 466 | m_chatter = app.m_chatter; |
| 469 | 467 | |
| ... | ... | @@ -500,17 +498,9 @@ void ctbkgcube::get_parameters(void) |
| 500 | 498 | // event lists, but do not accept counts cubes. |
| 501 | 499 | setup_observations(m_obs, true, true, false); |
| 502 | 500 | |
| 503 | - // If there is a valid event cube than use it for construction. | |
| 504 | - // Otherwise check if "incube" file name is valid then use the specified | |
| 505 | - // counts cube. Otherwise create a counts cube from the user | |
| 506 | - // parameters | |
| 501 | + // If no counts cube exists then load one from the "incube" filename | |
| 507 | 502 | if (m_cube.size() == 0) { |
| 508 | - if ((*this)["incube"].is_valid()) { | |
| 509 | - m_cube = GCTAEventCube((*this)["incube"].filename()); | |
| 510 | - } | |
| 511 | - else { | |
| 512 | - m_cube = create_cube(m_obs); | |
| 513 | - } | |
| 503 | + m_cube = GCTAEventCube((*this)["incube"].filename()); | |
| 514 | 504 | } |
| 515 | 505 | |
| 516 | 506 | // Define background cube |
| ... | ... | @@ -520,9 +510,8 @@ void ctbkgcube::get_parameters(void) |
| 520 | 510 | setup_models(m_obs); |
| 521 | 511 | |
| 522 | 512 | // Get remaining parameters |
| 523 | - m_addbounds = (*this)["addbounds"].boolean(); | |
| 524 | - m_publish = (*this)["publish"].boolean(); | |
| 525 | - m_chatter = static_cast<GChatter>((*this)["chatter"].integer()); | |
| 513 | + m_publish = (*this)["publish"].boolean(); | |
| 514 | + m_chatter = static_cast<GChatter>((*this)["chatter"].integer()); | |
| 526 | 515 | |
| 527 | 516 | // If needed later, query output filenames now |
| 528 | 517 | if (read_ahead()) { | ... | ... |
src/ctbkgcube/ctbkgcube.hpp
| ... | ... | @@ -73,7 +73,6 @@ protected: |
| 73 | 73 | // User parameters |
| 74 | 74 | GFilename m_outcube; //!< Filename of output cube |
| 75 | 75 | GFilename m_outmodel; //!< Filename of output XML model |
| 76 | - bool m_addbounds; //!< Add energies at boundaries? | |
| 77 | 76 | bool m_publish; //!< Publish background cube? |
| 78 | 77 | GChatter m_chatter; //!< Chattiness |
| 79 | 78 | ... | ... |
src/ctbkgcube/ctbkgcube.par
| ... | ... | @@ -17,36 +17,14 @@ |
| 17 | 17 | # General parameters |
| 18 | 18 | #=================== |
| 19 | 19 | inobs, f, a, NONE,,, "Input event list or observation definition XML file" |
| 20 | +incube, f, a, cntcube.fits,,, "Input counts cube file to extract background cube definition" | |
| 20 | 21 | inmodel, f, a, NONE,,, "Input model definition XML file" |
| 21 | -incube, f, a, NONE,,, "Input counts cube file to extract background cube definition" | |
| 22 | 22 | caldb, s, a, prod2,,, "Calibration database" |
| 23 | 23 | irf, s, a, South_0.5h,,, "Instrument response function" |
| 24 | 24 | outcube, f, a, bkgcube.fits,,, "Output background cube file" |
| 25 | 25 | outmodel, f, a, NONE,,, "Output model definition XML file" |
| 26 | 26 | |
| 27 | 27 | # |
| 28 | -# Energy binning parameters | |
| 29 | -#========================== | |
| 30 | -ebinalg, s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy bins" | |
| 31 | -emin, r, a, 0.1,,, "Lower energy limit (TeV)" | |
| 32 | -emax, r, a, 100.0,,, "Upper energy limit (TeV)" | |
| 33 | -enumbins, i, a, 20,1,200, "Number of energy bins" | |
| 34 | -ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" | |
| 35 | -addbounds, b, h, no,,, "Add energies at observation energy boundaries?" | |
| 36 | - | |
| 37 | -# | |
| 38 | -# Spatial binning parameters | |
| 39 | -#=========================== | |
| 40 | -usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" | |
| 41 | -nxpix, i, a, 200,,, "Size of the X axis in pixels" | |
| 42 | -nypix, i, a, 200,,, "Size of the Y axis in pixels" | |
| 43 | -binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" | |
| 44 | -coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" | |
| 45 | -proj, s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method" | |
| 46 | -xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" | |
| 47 | -yref, r, a, 22.51,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" | |
| 48 | - | |
| 49 | -# | |
| 50 | 28 | # Standard parameters |
| 51 | 29 | #==================== |
| 52 | 30 | publish, b, h, no,,, "Publish sky map on VO Hub?" | ... | ... |
test/test_ctbkgcube.py
| ... | ... | @@ -67,14 +67,11 @@ class Test(test): |
| 67 | 67 | |
| 68 | 68 | # Setup ctbkgcube command |
| 69 | 69 | cmd = ctbkgcube+' inobs="'+self._events+'"'+ \ |
| 70 | + ' incube="'+self._cntcube+'"'+ \ | |
| 70 | 71 | ' inmodel="'+self._model+'"'+ \ |
| 71 | - ' incube="NONE"'+ \ | |
| 72 | 72 | ' outcube="ctbkgcube_cmd1.fits"'+ \ |
| 73 | 73 | ' outmodel="ctbkgcube_cmd1.xml"'+ \ |
| 74 | 74 | ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ |
| 75 | - ' ebinalg="LOG" emin=1.0 emax=100.0 enumbins=5'+ \ | |
| 76 | - ' nxpix=10 nypix=10 binsz=0.4'+ \ | |
| 77 | - ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ | |
| 78 | 75 | ' logfile="ctbkgcube_cmd1.log" chatter=1' |
| 79 | 76 | |
| 80 | 77 | # Check if execution was successful |
| ... | ... | @@ -86,14 +83,11 @@ class Test(test): |
| 86 | 83 | |
| 87 | 84 | # Setup ctbkgcube command |
| 88 | 85 | cmd = ctbkgcube+' inobs="event_file_that_does_not_exist.fits"'+ \ |
| 86 | + ' incube="'+self._cntcube+'"'+ \ | |
| 89 | 87 | ' inmodel="'+self._model+'"'+ \ |
| 90 | - ' incube="NONE"'+ \ | |
| 91 | 88 | ' outcube="ctbkgcube_cmd2.fits"'+ \ |
| 92 | 89 | ' outmodel="ctbkgcube_cmd2.xml"'+ \ |
| 93 | 90 | ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ |
| 94 | - ' ebinalg="LOG" emin=1.0 emax=100.0 enumbins=5'+ \ | |
| 95 | - ' nxpix=10 nypix=10 binsz=0.4'+ \ | |
| 96 | - ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ | |
| 97 | 91 | ' logfile="ctbkgcube_cmd2.log" debug=yes chatter=1' |
| 98 | 92 | |
| 99 | 93 | # Check if execution failed |
| ... | ... | @@ -153,21 +147,10 @@ class Test(test): |
| 153 | 147 | |
| 154 | 148 | # Now set ctbkgcube parameters |
| 155 | 149 | bkgcube['inobs'] = self._events |
| 150 | + bkgcube['incube'] = self._cntcube | |
| 156 | 151 | bkgcube['inmodel'] = self._model |
| 157 | - bkgcube['incube'] = 'NONE' | |
| 158 | 152 | bkgcube['caldb'] = self._caldb |
| 159 | 153 | bkgcube['irf'] = self._irf |
| 160 | - bkgcube['ebinalg'] = 'LOG' | |
| 161 | - bkgcube['emin'] = 1.0 | |
| 162 | - bkgcube['emax'] = 100.0 | |
| 163 | - bkgcube['enumbins'] = 5 | |
| 164 | - bkgcube['nxpix'] = 10 | |
| 165 | - bkgcube['nypix'] = 10 | |
| 166 | - bkgcube['binsz'] = 0.4 | |
| 167 | - bkgcube['coordsys'] = 'CEL' | |
| 168 | - bkgcube['proj'] = 'CAR' | |
| 169 | - bkgcube['xref'] = 83.63 | |
| 170 | - bkgcube['yref'] = 22.01 | |
| 171 | 154 | bkgcube['outcube'] = 'ctbkgcube_py3.fits' |
| 172 | 155 | bkgcube['outmodel'] = 'ctbkgcube_py3.xml' |
| 173 | 156 | bkgcube['logfile'] = 'ctbkgcube_py3.log' |
| ... | ... | @@ -189,15 +172,12 @@ class Test(test): |
| 189 | 172 | self._check_models(cpy_bkgcube.models()) |
| 190 | 173 | |
| 191 | 174 | # Execute copy of ctbkgcube tool again, now with a higher chatter |
| 192 | - # level than before. In addition, use counts cube to define the | |
| 193 | - # background cube | |
| 194 | - cpy_bkgcube['incube'] = self._cntcube | |
| 195 | - cpy_bkgcube['outcube'] = 'ctbkgcube_py4.fits' | |
| 196 | - cpy_bkgcube['outmodel'] = 'ctbkgcube_py4.xml' | |
| 197 | - cpy_bkgcube['logfile'] = 'ctbkgcube_py4.log' | |
| 198 | - cpy_bkgcube['chatter'] = 3 | |
| 199 | - cpy_bkgcube['publish'] = True | |
| 200 | - cpy_bkgcube['addbounds'] = True | |
| 175 | + # level than before. In addition, publish the background cube. | |
| 176 | + cpy_bkgcube['outcube'] = 'ctbkgcube_py4.fits' | |
| 177 | + cpy_bkgcube['outmodel'] = 'ctbkgcube_py4.xml' | |
| 178 | + cpy_bkgcube['logfile'] = 'ctbkgcube_py4.log' | |
| 179 | + cpy_bkgcube['chatter'] = 3 | |
| 180 | + cpy_bkgcube['publish'] = True | |
| 201 | 181 | cpy_bkgcube.logFileOpen() # Make sure we get a log file |
| 202 | 182 | cpy_bkgcube.execute() |
| 203 | 183 | |
| ... | ... | @@ -218,25 +198,13 @@ class Test(test): |
| 218 | 198 | |
| 219 | 199 | # Set-up ctbkgcube from observation container |
| 220 | 200 | bkgcube = ctools.ctbkgcube(obs) |
| 221 | - bkgcube['incube'] = '' | |
| 222 | - bkgcube['caldb'] = self._caldb | |
| 223 | - bkgcube['irf'] = self._irf | |
| 224 | - bkgcube['ebinalg'] = 'LOG' | |
| 225 | - bkgcube['emin'] = 0.2 | |
| 226 | - bkgcube['emax'] = 150.0 | |
| 227 | - bkgcube['enumbins'] = 5 | |
| 228 | - bkgcube['nxpix'] = 10 | |
| 229 | - bkgcube['nypix'] = 10 | |
| 230 | - bkgcube['binsz'] = 0.4 | |
| 231 | - bkgcube['coordsys'] = 'CEL' | |
| 232 | - bkgcube['proj'] = 'CAR' | |
| 233 | - bkgcube['xref'] = 83.63 | |
| 234 | - bkgcube['yref'] = 22.01 | |
| 235 | - bkgcube['outcube'] = 'ctbkgcube_py5.fits' | |
| 236 | - bkgcube['outmodel'] = 'ctbkgcube_py5.xml' | |
| 237 | - bkgcube['logfile'] = 'ctbkgcube_py5.log' | |
| 238 | - bkgcube['addbounds'] = True | |
| 239 | - bkgcube['chatter'] = 4 | |
| 201 | + bkgcube['incube'] = self._cntcube | |
| 202 | + bkgcube['caldb'] = self._caldb | |
| 203 | + bkgcube['irf'] = self._irf | |
| 204 | + bkgcube['outcube'] = 'ctbkgcube_py5.fits' | |
| 205 | + bkgcube['outmodel'] = 'ctbkgcube_py5.xml' | |
| 206 | + bkgcube['logfile'] = 'ctbkgcube_py5.log' | |
| 207 | + bkgcube['chatter'] = 4 | |
| 240 | 208 | |
| 241 | 209 | # Execute ctbkgcube tool |
| 242 | 210 | bkgcube.logFileOpen() # Make sure we get a log file |
| ... | ... | @@ -306,6 +274,6 @@ class Test(test): |
| 306 | 274 | """ |
| 307 | 275 | # Check number of models |
| 308 | 276 | self.test_value(models.size(), nmodels, 'Check number of models') |
| 309 | - | |
| 277 | + | |
| 310 | 278 | # Return |
| 311 | 279 | return | ... | ... |
test/test_obsutils.py
| ... | ... | @@ -379,7 +379,7 @@ class Test(test): |
| 379 | 379 | |
| 380 | 380 | # Setup one CTA observation |
| 381 | 381 | res = obsutils.set_obs(pnt, emin=1.0, emax=10.0) |
| 382 | - | |
| 382 | + | |
| 383 | 383 | # Check result |
| 384 | 384 | self.test_value(res.eventtype(), 'EventList', 'Check event type') |
| 385 | 385 | self.test_value(res.events().ebounds().emin().TeV(), 1.0, |
| ... | ... | @@ -498,94 +498,75 @@ class Test(test): |
| 498 | 498 | # Set-up observation container |
| 499 | 499 | obs = self._setup_sim() |
| 500 | 500 | |
| 501 | - # Get stacked response with xref/yref set | |
| 502 | - res = obsutils.get_stacked_response(obs, 83.6331, 22.0145, | |
| 503 | - coordsys='CEL', | |
| 504 | - binsz=0.1, nxpix=10, nypix=10, | |
| 505 | - emin=1.0, emax=100.0, enumbins=5, | |
| 506 | - edisp=False) | |
| 501 | + # Set-up counts cube | |
| 502 | + map = gammalib.GSkyMap('CAR','CEL',83.6331,22.0145,0.1,0.1,10,10,5) | |
| 503 | + emin = gammalib.GEnergy(0.1, 'TeV') | |
| 504 | + emax = gammalib.GEnergy(100.0, 'TeV') | |
| 505 | + ebds = gammalib.GEbounds(5, emin, emax) | |
| 506 | + tmin = gammalib.GTime(0.0) | |
| 507 | + tmax = gammalib.GTime(1000.0) | |
| 508 | + gti = gammalib.GGti(tmin, tmax) | |
| 509 | + cntcube = gammalib.GCTAEventCube(map, ebds, gti) | |
| 510 | + | |
| 511 | + # Get stacked response | |
| 512 | + res = obsutils.get_stacked_response(obs, cntcube, edisp=False) | |
| 507 | 513 | |
| 508 | 514 | # Check result |
| 509 | 515 | self.test_value(res['expcube'].cube().npix(), 100, |
| 510 | 516 | 'Check number of exposure cube pixels') |
| 511 | - self.test_value(res['expcube'].cube().nmaps(), 61, | |
| 517 | + self.test_value(res['expcube'].cube().nmaps(), 91, | |
| 512 | 518 | 'Check number of exposure cube maps') |
| 513 | 519 | self.test_value(res['psfcube'].cube().npix(), 4, |
| 514 | 520 | 'Check number of PSF cube pixels') |
| 515 | - self.test_value(res['psfcube'].cube().nmaps(), 12200, | |
| 521 | + self.test_value(res['psfcube'].cube().nmaps(), 18200, | |
| 516 | 522 | 'Check number of PSF cube maps') |
| 517 | 523 | self.test_value(res['bkgcube'].cube().npix(), 100, |
| 518 | 524 | 'Check number of background cube pixels') |
| 519 | 525 | self.test_value(res['bkgcube'].cube().nmaps(), 5, |
| 520 | 526 | 'Check number of background cube maps') |
| 521 | 527 | |
| 522 | - # Get stacked response with None xref/yref | |
| 523 | - res = obsutils.get_stacked_response(obs, None, None, | |
| 524 | - coordsys='CEL', | |
| 525 | - binsz=0.1, nxpix=10, nypix=10, | |
| 526 | - emin=1.0, emax=100.0, enumbins=5, | |
| 527 | - edisp=False) | |
| 528 | + # Get stacked response with energy dispersion | |
| 529 | + res = obsutils.get_stacked_response(obs, cntcube, edisp=True) | |
| 528 | 530 | |
| 529 | 531 | # Check result |
| 530 | 532 | self.test_value(res['expcube'].cube().npix(), 100, |
| 531 | 533 | 'Check number of exposure cube pixels') |
| 532 | - self.test_value(res['expcube'].cube().nmaps(), 61, | |
| 534 | + self.test_value(res['expcube'].cube().nmaps(), 105, | |
| 533 | 535 | 'Check number of exposure cube maps') |
| 534 | 536 | self.test_value(res['psfcube'].cube().npix(), 4, |
| 535 | 537 | 'Check number of PSF cube pixels') |
| 536 | - self.test_value(res['psfcube'].cube().nmaps(), 12200, | |
| 538 | + self.test_value(res['psfcube'].cube().nmaps(), 21000, | |
| 537 | 539 | 'Check number of PSF cube maps') |
| 538 | 540 | self.test_value(res['bkgcube'].cube().npix(), 100, |
| 539 | 541 | 'Check number of background cube pixels') |
| 540 | 542 | self.test_value(res['bkgcube'].cube().nmaps(), 5, |
| 541 | 543 | 'Check number of background cube maps') |
| 544 | + self.test_value(res['edispcube'].cube().npix(), 4, | |
| 545 | + 'Check number of energy dispersion cube pixels') | |
| 546 | + self.test_value(res['edispcube'].cube().nmaps(), 10500, | |
| 547 | + 'Check number of energy dispersion cube maps') | |
| 548 | + | |
| 549 | + # Set-up counts cube with large number of energy bins | |
| 550 | + ebds = gammalib.GEbounds(100, emin, emax) | |
| 551 | + cntcube = gammalib.GCTAEventCube(map, ebds, gti) | |
| 542 | 552 | |
| 543 | 553 | # Get stacked response with xref/yref set and large number of energy bins |
| 544 | - res = obsutils.get_stacked_response(obs, 83.6331, 22.0145, | |
| 545 | - coordsys='CEL', | |
| 546 | - binsz=0.1, nxpix=10, nypix=10, | |
| 547 | - emin=1.0, emax=100.0, enumbins=100, | |
| 548 | - edisp=False) | |
| 554 | + res = obsutils.get_stacked_response(obs, cntcube, edisp=False) | |
| 549 | 555 | |
| 550 | 556 | # Check result |
| 551 | 557 | self.test_value(res['expcube'].cube().npix(), 100, |
| 552 | 558 | 'Check number of exposure cube pixels') |
| 553 | - self.test_value(res['expcube'].cube().nmaps(), 101, | |
| 559 | + self.test_value(res['expcube'].cube().nmaps(), 91, | |
| 554 | 560 | 'Check number of exposure cube maps') |
| 555 | 561 | self.test_value(res['psfcube'].cube().npix(), 4, |
| 556 | 562 | 'Check number of PSF cube pixels') |
| 557 | - self.test_value(res['psfcube'].cube().nmaps(), 20200, | |
| 563 | + self.test_value(res['psfcube'].cube().nmaps(), 18200, | |
| 558 | 564 | 'Check number of PSF cube maps') |
| 559 | 565 | self.test_value(res['bkgcube'].cube().npix(), 100, |
| 560 | 566 | 'Check number of background cube pixels') |
| 561 | 567 | self.test_value(res['bkgcube'].cube().nmaps(), 100, |
| 562 | 568 | 'Check number of background cube maps') |
| 563 | 569 | |
| 564 | - # Get stacked response with energy dispersion | |
| 565 | - res = obsutils.get_stacked_response(obs, 83.6331, 22.0145, | |
| 566 | - coordsys='CEL', | |
| 567 | - binsz=0.1, nxpix=10, nypix=10, | |
| 568 | - emin=1.0, emax=100.0, enumbins=5, | |
| 569 | - edisp=True) | |
| 570 | - | |
| 571 | - # Check result | |
| 572 | - self.test_value(res['expcube'].cube().npix(), 100, | |
| 573 | - 'Check number of exposure cube pixels') | |
| 574 | - self.test_value(res['expcube'].cube().nmaps(), 75, | |
| 575 | - 'Check number of exposure cube maps') | |
| 576 | - self.test_value(res['psfcube'].cube().npix(), 4, | |
| 577 | - 'Check number of PSF cube pixels') | |
| 578 | - self.test_value(res['psfcube'].cube().nmaps(), 15000, | |
| 579 | - 'Check number of PSF cube maps') | |
| 580 | - self.test_value(res['bkgcube'].cube().npix(), 100, | |
| 581 | - 'Check number of background cube pixels') | |
| 582 | - self.test_value(res['bkgcube'].cube().nmaps(), 5, | |
| 583 | - 'Check number of background cube maps') | |
| 584 | - self.test_value(res['edispcube'].cube().npix(), 4, | |
| 585 | - 'Check number of energy dispersion cube pixels') | |
| 586 | - self.test_value(res['edispcube'].cube().nmaps(), 7500, | |
| 587 | - 'Check number of energy dispersion cube maps') | |
| 588 | - | |
| 589 | 570 | # Return |
| 590 | 571 | return |
| 591 | 572 | ... | ... |