ctools / ctools

Commit 414b1c9b3c7ab1eed0d823d76cb33998715449c9

Authored by Jürgen Knödlseder
1 parent 1cf2fb6d

Add BGDLIXF algorithm

Inline Side-by-side
Showing 2 changed files with 198 additions and 4 deletions
modules/comscripts/comobsback.par
  View file @414b1c9
... ... @@ -25,7 +25,7 @@ outobs, f, a, obs_back.xml,,, "Output observation definition file"
25 25 #
26 26 # Script parameters
27 27 #==================
28   -bkgmethod, s, a, BGDLIXE,PHINOR|BGDLIXA|BGDLIXE,, "Method for background computation"
  28 +bkgmethod, s, a, BGDLIXE,PHINOR|BGDLIXA|BGDLIXE|BGDLIXF,, "Method for background computation"
29 29 nrunav, i, h, 3,,, "Number of Chi/Psi bins used for running average"
30 30 navgr, i, h, 5,,, "Number of Chi/Psi bins used for averaging"
31 31 nincl, i, h, 15,,, "Number of Phibar layers to include"
... ...
modules/comscripts/comobsback.py
  View file @414b1c9
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # Generate background model for COMPTEL observations
4 4 #
5   -# Copyright (C) 2021-2022 Juergen Knoedlseder
  5 +# Copyright (C) 2021-2023 Juergen Knoedlseder
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -716,6 +716,195 @@ class comobsback(ctools.csobservation):
716 716 # Return
717 717 return drbname
718 718  
  719 + def _generate_drb_bgdlixF(self, drename, drwname, drm, debug=False):
  720 + """
  721 + Generate DRB from DRW and DRE file with BGDLIX method F
  722 +
  723 + The method generates a background model using the BGDLIX method F,
  724 + taking into account a source model DRM so that the method can be
  725 + used for an iterative SRCLIX approach.
  726 +
  727 + Parameters
  728 + ----------
  729 + drename : str
  730 + DRE filename
  731 + drwname : str
  732 + DRW filename
  733 + drm : `~gammalib.GCOMDri`
  734 + Source DRM
  735 + debug : bool
  736 + Write some files for debugging (default: False)
  737 +
  738 + Returns
  739 + -------
  740 + drbname : str
  741 + DRB filename
  742 + """
  743 + # Get suffix and outfolder
  744 + suffix = self['suffix'].string()
  745 + outfolder = self['outfolder'].string()
  746 +
  747 + # Get task parameters
  748 + navgr = self['navgr'].integer()
  749 + nincl = self['nincl'].integer()
  750 + nexcl = self['nexcl'].integer()
  751 +
  752 + # Set DRB filename
  753 + drbname = '%s/%s' % (outfolder, os.path.basename(drename.replace('dre',
  754 + 'drb-bgdlixF-na%d-ni%d-ne%d%s' %
  755 + (navgr, nincl, nexcl, suffix))))
  756 + drbfile = gammalib.GFilename(drbname)
  757 +
  758 + # Write header
  759 + self._log_header3(gammalib.NORMAL, 'Compute DRB using BGDLIX method F')
  760 + self._log_value(gammalib.NORMAL, 'Bins used for averaging', navgr)
  761 + self._log_value(gammalib.NORMAL, 'Phibar layers to include', nincl)
  762 + self._log_value(gammalib.NORMAL, 'Phibar layers to exclude', nexcl)
  763 +
  764 + # If DRB file exists then do nothing
  765 + if drbfile.exists():
  766 + self._log_value(gammalib.NORMAL, 'DRB file exists', drbfile.url())
  767 +
  768 + # ... otherwise compute and save it
  769 + else:
  770 +
  771 + # Set DRE and DRW filenames
  772 + drefile = gammalib.GFilename(drename)
  773 + drwfile = gammalib.GFilename(drwname)
  774 +
  775 + # Load DRE and DRW
  776 + dre = gammalib.GCOMDri(drefile)
  777 + drw = gammalib.GCOMDri(drwfile)
  778 +
  779 + # Get dataspace dimensions
  780 + nchi = dre.nchi()
  781 + npsi = dre.npsi()
  782 + nphibar = dre.nphibar()
  783 + npix = nchi * npsi
  784 +
  785 + # Precompute half running average lengths
  786 + navgr2 = int(navgr/2)
  787 + nexcl2 = int(nexcl/2)
  788 + nincl2 = int(nincl/2)
  789 +
  790 + # Initialise SCRAT and DRB
  791 + scrat = drw.copy()
  792 + drb = drw.copy()
  793 +
  794 + # Do Phibar normalization of DRW and store result in SCRAT.
  795 + for i3 in range(nphibar):
  796 + sum_dre = 0.0
  797 + sum_drw = 0.0
  798 + sum_drm = 0.0
  799 + for i2 in range(npsi):
  800 + for i1 in range(nchi):
  801 + ipix = i3*npix + i2*nchi + i1
  802 + sum_dre += dre[ipix]
  803 + sum_drw += drw[ipix]
  804 + sum_drm += drm[ipix]
  805 + if sum_drw > 0.0:
  806 + scale = (sum_dre-sum_drm) / sum_drw
  807 + for i2 in range(npsi):
  808 + for i1 in range(nchi):
  809 + ipix = i3*npix + i2*nchi + i1
  810 + scrat[ipix] = drw[ipix] * scale
  811 +
  812 + # Debug: save B^0_L
  813 + if debug:
  814 + srcatname = drbname.replace('drb', 'scrat1')
  815 + srcatfile = gammalib.GFilename(srcatname)
  816 + scrat.save(srcatfile, self['clobber'].boolean())
  817 +
  818 + # Make background model
  819 + for i1 in range(nchi):
  820 + for i2 in range(npsi):
  821 +
  822 + # Loop over Phibar
  823 + for i3 in range(nphibar):
  824 +
  825 + # Determine index range for sum. There are two sums that
  826 + # go over [isel1, iex1[ and [ixe2, isel2[
  827 + iex1 = i3 - nexcl2
  828 + iex2 = i3 + nexcl2 + 1
  829 + if nexcl == 0:
  830 + iex2 -= 1
  831 + isel1 = max((i3 - nincl2), 0)
  832 + isel2 = min((i3 + nincl2 + 1), nphibar)
  833 + if isel1 == 0:
  834 + isel2 = min(nincl, nphibar)
  835 + if isel2 == nphibar:
  836 + isel1 = max(nphibar - nincl, 0)
  837 +
  838 + # Initialise sums
  839 + sum_dre = 0.0
  840 + sum_drm = 0.0
  841 + sum_scrat = 0.0
  842 +
  843 + # Compute average over small patch in Chi/Psi
  844 + for j1 in range(max(0, i1 - navgr2), min(nchi, i1 + navgr2 + 1)):
  845 + for j2 in range(max(0, i2 - navgr2), min(npsi, i2 + navgr2 + 1)):
  846 +
  847 + # Take sum over [isel1, iex1[
  848 + if iex1 >= 1:
  849 + for j3 in range(isel1, iex1):
  850 + jpix = j1 + j2*nchi + j3*npix
  851 + sum_dre += dre[jpix]
  852 + sum_drm += drm[jpix]
  853 + sum_scrat += scrat[jpix]
  854 +
  855 + # Take sum over [iex2, isel2[
  856 + if iex2 < nphibar:
  857 + for j3 in range(iex2, isel2):
  858 + jpix = j1 + j2*nchi + j3*npix
  859 + sum_dre += dre[jpix]
  860 + sum_drm += drm[jpix]
  861 + sum_scrat += scrat[jpix]
  862 +
  863 + # Re-normalise SCRAT to derive DRB
  864 + ipix = i1 + i2*nchi + i3*npix
  865 + if sum_scrat > 0.0:
  866 + drb[ipix] = scrat[ipix] * (sum_dre-sum_drm)/sum_scrat
  867 + else:
  868 + drb[ipix] = 0.0
  869 +
  870 + # Optionally Phibar normalise DRB
  871 + if self['phinor'].boolean():
  872 + for i3 in range(nphibar):
  873 + sum_dre = 0.0
  874 + sum_drb = 0.0
  875 + sum_drm = 0.0
  876 + for i in range(npix):
  877 + ipix = i3*npix + i
  878 + sum_dre += dre[ipix]
  879 + sum_drb += drb[ipix]
  880 + sum_drm += drm[ipix]
  881 + if sum_drb > 0.0:
  882 + norm = (sum_dre-sum_drm) / sum_drb
  883 + for i in range(npix):
  884 + ipix = i3*npix + i
  885 + drb[ipix] *= norm
  886 +
  887 + # Get statistics
  888 + sum_dre = 0.0
  889 + sum_drm = 0.0
  890 + sum_drb = 0.0
  891 + for i in range(dre.size()):
  892 + sum_dre += dre[i]
  893 + sum_drm += drm[i]
  894 + sum_drb += drb[i]
  895 + self._log_value(gammalib.NORMAL, 'Events in DRE', sum_dre)
  896 + self._log_value(gammalib.NORMAL, 'Events in DRM', sum_drm)
  897 + self._log_value(gammalib.NORMAL, 'Events in DRB', sum_drb)
  898 +
  899 + # Save DRB
  900 + drb.save(drbfile, self['clobber'].boolean())
  901 +
  902 + # Log creation
  903 + self._log_value(gammalib.NORMAL, 'DRB file created', drbfile.url())
  904 +
  905 + # Return
  906 + return drbname
  907 +
719 908  
720 909 # Public methods
721 910 def process(self):
... ... @@ -766,20 +955,25 @@ class comobsback(ctools.csobservation):
766 955 # Write header
767 956 self._log_header2(gammalib.NORMAL, self._get_obs_header(obs))
768 957  
769   - # Get DRE and DRG names
  958 + # Get DRE name
770 959 drename = obs.drename().url()
771   - drgname = obs.drgname().url()
772 960  
773 961 # Compute DRM
774 962 drm = obs.drm(models)
775 963  
776 964 # Generate background
777 965 if method == 'PHINOR':
  966 + drgname = obs.drgname().url()
778 967 drbname = self._generate_drb_phinor(drename, drgname, drm)
779 968 elif method == 'BGDLIXA':
  969 + drgname = obs.drgname().url()
780 970 drbname = self._generate_drb_bgdlixA(drename, drgname, drm)
781 971 elif method == 'BGDLIXE':
  972 + drgname = obs.drgname().url()
782 973 drbname = self._generate_drb_bgdlixE(drename, drgname, drm)
  974 + elif method == 'BGDLIXF':
  975 + drwname = obs.drwname().url()
  976 + drbname = self._generate_drb_bgdlixF(drename, drwname, drm)
783 977  
784 978 # Set DRB name
785 979 obs.drbname(drbname)
... ...