ctools / ctools

Commit aeaeaf4178a8af0b066518481b66c20453a07510

Authored by Jürgen Knödlseder
1 parent 414b1c9b

Support DRW in comobsadd

Inline Side-by-side
Showing 8 changed files with 132 additions and 19 deletions
ChangeLog
  View file @aeaeaf4
1   -2023-01-18
  1 +2023-02-15
2 2  
3 3 * Version 2.1.0 released
4 4 ========================
5 5  
  6 + Add support for DRW weighting cubes (#4209)
6 7 Prefit models without test source in comlixmap (#4201)
7 8 Add inmap parameter to comlixmap (#4187)
8 9 Add comobsconv script (#4159)
... ...
1 1 New Features and Important Changes in ctools 2.1.0
2 2  
3   -18 January 2023
  3 +15 February 2023
4 4  
5 5  
6 6 Introduction
... ... @@ -342,17 +342,22 @@ parameter uncertainties in the resulting TS map FITS file (#4187).
342 342  
343 343 comobsadd - Combine COMPTEL observations
344 344 ----------------------------------------
345   -None
  345 +Added DRW to the cubes that will be combined by the script. The combination
  346 +will be done using a simple addition, alike the combination of the DRB cubes
  347 +(#4209).
346 348  
347 349  
348 350 comobsback - Generate background model for COMPTEL observations
349 351 ---------------------------------------------------------------
350   -None
  352 +Added BGDLIXF background generation method that makes use of the DRW instead
  353 +a DRG for background modelling. Set "outfolder" parameter to automatic (#4209).
351 354  
352 355  
353 356 comobsbin - Bin COMPTEL observations
354 357 ------------------------------------
355   -None
  358 +Implement computation of weighting cubes (DRW) for improved background modelling.
  359 +Added "timebin" parameter to provide control over the time binning used for
  360 +DRW computation (#4209).
356 361  
357 362  
358 363 comobsconv - Convolve models with COMPTEL response
... ...
README.md
  View file @aeaeaf4
1 1 ctools information
2 2 ==================
3   -* Version: 2.1.0.dev (18 January 2023)
  3 +* Version: 2.1.0.dev (15 February 2023)
4 4 * GammaLib dependency: 2.1.0.dev
5 5  
6 6 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/)
... ...
doc/source/admin/release_history/2.1.rst
  View file @aeaeaf4
... ... @@ -10,7 +10,7 @@ ctools 2.1 is a major release that adds significant functionality.
10 10  
11 11 In particular, this release provides:
12 12  
13   -*
  13 +* support for COMPTEL DRW weighting cubes for improved background modelling
14 14  
15 15  
16 16 Bug fixes
... ... @@ -28,6 +28,8 @@ Bug fixes
28 28 Improvements
29 29 ------------
30 30  
  31 +* [`4209 <https://cta-redmine.irap.omp.eu/issues/4209>`_] -
  32 + Add support for DRW weighting cubes
31 33 * [`4201 <https://cta-redmine.irap.omp.eu/issues/4201>`_] -
32 34 Prefit models without test source in :ref:`comlixmap`
33 35  
... ...
modules/comscripts/comobsadd.py
  View file @aeaeaf4
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # Add 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
... ... @@ -43,6 +43,7 @@ class comobsadd(ctools.csobservation):
43 43 self._xml = gammalib.GXml()
44 44 self._dres = []
45 45 self._drbs = []
  46 + self._drws = []
46 47 self._drg = []
47 48 self._drx = []
48 49 self._nspuse = 0
... ... @@ -144,6 +145,7 @@ class comobsadd(ctools.csobservation):
144 145 if ebin['emin'] == emin and ebin['emax'] == emax:
145 146 ebin['dres'].append(o.drename())
146 147 ebin['drbs'].append(o.drbname())
  148 + ebin['drws'].append(o.drwname())
147 149 ebin['drgs'].append(o.drgname())
148 150 ebin['drxs'].append(o.drxname())
149 151 ebin['gti'].extend(o.events().dre().gti())
... ... @@ -157,6 +159,7 @@ class comobsadd(ctools.csobservation):
157 159 'emax' : emax,
158 160 'dres' : [o.drename()],
159 161 'drbs' : [o.drbname()],
  162 + 'drws' : [o.drwname()],
160 163 'drgs' : [o.drgname()],
161 164 'drxs' : [o.drxname()],
162 165 'iaq' : o.response().rspname(),
... ... @@ -185,22 +188,24 @@ class comobsadd(ctools.csobservation):
185 188 self._log_value(gammalib.NORMAL, ' Phase correction', ebin['phasecor'])
186 189 self._log_value(gammalib.NORMAL, ' Response', ebin['iaq'])
187 190  
188   - # Create sky map for DRE, DRG and DRB
  191 + # Create sky map for DRE, DRB, DRW and DRG
189 192 cube = gammalib.GSkyMap(self['proj'].string(), coordsys,
190 193 chi0, psi0, -dchi, dpsi, nchi, npsi)
191 194  
192 195 # Create sky map for DRX
193 196 expo = gammalib.GSkyMap('CAR', 'GAL', 0.0, 0.0, -1.0, 1.0, 360, 180)
194 197  
195   - # Allocate DRE, DRG, DRB and DRX
  198 + # Allocate DRE, DRB, DRW, DRG and DRX
196 199 dre = gammalib.GCOMDri(cube, 0.0, dphibar, nphibar)
197 200 drb = gammalib.GCOMDri(cube, 0.0, dphibar, nphibar)
  201 + drw = gammalib.GCOMDri(cube, 0.0, dphibar, nphibar)
198 202 self._drg = gammalib.GCOMDri(cube, 0.0, dphibar, nphibar)
199 203 self._drx = gammalib.GCOMDri(expo)
200 204  
201   - # Set GTIs of DRE, DRB, DRG and DRX
  205 + # Set GTIs of DRE, DRB, DRW, DRG and DRX
202 206 dre.gti(ebin['gti'])
203 207 drb.gti(ebin['gti'])
  208 + drw.gti(ebin['gti'])
204 209 self._drg.gti(ebin['gti'])
205 210 self._drx.gti(ebin['gti'])
206 211  
... ... @@ -218,7 +223,7 @@ class comobsadd(ctools.csobservation):
218 223 self._drx.num_used_superpackets(ebins[0]['nspuse'])
219 224 self._drx.num_skipped_superpackets(ebins[0]['nspskp'])
220 225  
221   - # Create DRE and DRB for each energy bin and append an XML entry
  226 + # Create DRE, DRB and DRW for each energy bin and append an XML entry
222 227 for ebin in ebins:
223 228  
224 229 # Set observation name and id
... ... @@ -230,6 +235,7 @@ class comobsadd(ctools.csobservation):
230 235 ebounds.append(ebin['emin'], ebin['emax'])
231 236 dre.ebounds(ebounds)
232 237 drb.ebounds(ebounds)
  238 + drw.ebounds(ebounds)
233 239  
234 240 # Set ToF correction
235 241 dre.tof_correction(ebin['tofcor'])
... ... @@ -244,22 +250,28 @@ class comobsadd(ctools.csobservation):
244 250 drb.num_superpackets(ebin['nspinp'])
245 251 drb.num_used_superpackets(ebin['nspuse'])
246 252 drb.num_skipped_superpackets(ebin['nspskp'])
  253 + drw.num_superpackets(ebin['nspinp'])
  254 + drw.num_used_superpackets(ebin['nspuse'])
  255 + drw.num_skipped_superpackets(ebin['nspskp'])
247 256  
248   - # Set DRE and DRB filenames
  257 + # Set DRE, DRB and DRW filenames
249 258 dre.name('%s/%s_%s_dre.fits' % (self['outfolder'].string(), prefix, id))
250 259 drb.name('%s/%s_%s_drb.fits' % (self['outfolder'].string(), prefix, id))
  260 + drw.name('%s/%s_%s_drw.fits' % (self['outfolder'].string(), prefix, id))
251 261  
252   - # Put DRE and DRB in lists
  262 + # Put DRE, DRB and DRW in lists
253 263 self._dres.append(dre.copy())
254 264 self._drbs.append(drb.copy())
  265 + self._drws.append(drw.copy())
255 266  
256 267 # Append observation
257 268 xml_obs = xml_list.append('observation name="%s" id="%s" instrument="COM"' % \
258 269 (name, id))
259 270  
260   - # Append DRE, DRB, DRG and DRX files
  271 + # Append DRE, DRB, DRW, DRG and DRX files
261 272 xml_obs.append('parameter name="DRE" file="%s"' % (dre.name()))
262 273 xml_obs.append('parameter name="DRB" file="%s"' % (drb.name()))
  274 + xml_obs.append('parameter name="DRW" file="%s"' % (drw.name()))
263 275 xml_obs.append('parameter name="DRG" file="%s"' % (self._drg.name()))
264 276 xml_obs.append('parameter name="DRX" file="%s"' % (self._drx.name()))
265 277  
... ... @@ -287,6 +299,9 @@ class comobsadd(ctools.csobservation):
287 299 # Combine DRBs
288 300 self._combine_drbs(ebins)
289 301  
  302 + # Combine DRWs
  303 + self._combine_drws(ebins)
  304 +
290 305 # Combine DRGs
291 306 self._combine_drgs(ebins)
292 307  
... ... @@ -460,6 +475,82 @@ class comobsadd(ctools.csobservation):
460 475 # Return
461 476 return
462 477  
  478 + def _combine_drws(self, ebins):
  479 + """
  480 + Combine DRWs
  481 +
  482 + The DRWs are combined by dividing the original pixels by their solid
  483 + angles and multipliying with the solid angle of the combined DRW. In
  484 + that way the solid angle changes due to different projections are
  485 + correctly taken into account.
  486 +
  487 + Parameters
  488 + ----------
  489 + ebins : list of dict
  490 + List with input DRIs
  491 + """
  492 + # Write header
  493 + self._log_header2(gammalib.NORMAL, 'Combine DRWs')
  494 +
  495 + # Loop over all energy bins
  496 + for iebin, ebin in enumerate(ebins):
  497 +
  498 + # Write header
  499 + self._log_header3(gammalib.NORMAL, 'Energy bin %s - %s' %
  500 + (str(ebin['emin']), str(ebin['emax'])))
  501 +
  502 + # Get dimensions of combined data space
  503 + npix = self._drws[iebin].nchi() * self._drws[iebin].npsi()
  504 + nphibar = self._drws[iebin].nphibar()
  505 +
  506 + # Set combined data space map
  507 + map = self._drws[iebin].map()
  508 +
  509 + # Loop over DRWs
  510 + for drwname in ebin['drws']:
  511 +
  512 + # Log DRW file name
  513 + self._log_value(gammalib.NORMAL, 'Filename', drwname.url())
  514 +
  515 + # Load DRW
  516 + drw = gammalib.GCOMDri(drwname)
  517 +
  518 + # Divide pixels by solid angle so that the DRW returns the
  519 + # number of background counts per steradian
  520 + nmaps = drw.map().nmaps()
  521 + npixels = drw.map().npix()
  522 + for ipixel in range(npixels):
  523 + omega = drw.map().solidangle(ipixel)
  524 + for imap in range(nmaps):
  525 + drw.map()[ipixel,imap] /= omega
  526 +
  527 + # Loop over combined data space Chi/Psi pixels
  528 + for ipix in range(npix):
  529 +
  530 + # Get Chi/Psi direction of pixel
  531 + dir = map.inx2dir(ipix)
  532 +
  533 + # If Chi/Psi is contained in DRW then extract value.
  534 + # Put the check in a try-except block since there
  535 + # may be coordinate transformation issues at this
  536 + # point
  537 + try:
  538 + if drw.map().contains(dir):
  539 +
  540 + # Get solid angle of target map to convert the
  541 + # DRB values into absolute counts
  542 + omega = map.solidangle(ipix)
  543 +
  544 + # Loop over combined data space Phibar layers
  545 + for iphibar in range(nphibar):
  546 + map[ipix, iphibar] += drw.map()(dir, iphibar) * omega
  547 +
  548 + except RuntimeError:
  549 + pass
  550 +
  551 + # Return
  552 + return
  553 +
463 554 def _combine_drgs(self, ebins):
464 555 """
465 556 Combine DRGs
... ... @@ -630,6 +721,11 @@ class comobsadd(ctools.csobservation):
630 721 drb.save(drb.name(), self['clobber'].boolean())
631 722 self._log_value(gammalib.NORMAL, 'DRB file', drb.name())
632 723  
  724 + # Save DRWs
  725 + for drw in self._drws:
  726 + drw.save(drw.name(), self['clobber'].boolean())
  727 + self._log_value(gammalib.NORMAL, 'DRW file', drw.name())
  728 +
633 729 # Save DRG
634 730 self._drg.save(self._drg.name(), self['clobber'].boolean())
635 731 self._log_value(gammalib.NORMAL, 'DRG file', self._drg.name())
... ...
modules/comscripts/comobsback.par
  View file @aeaeaf4
... ... @@ -19,7 +19,7 @@
19 19 inobs, f, a, obs.xml,,, "Input observation definition file"
20 20 inmodel, f, a, models.xml,,, "Input model definition file"
21 21 suffix, s, h, ,,, "Suffix for DRB files"
22   -outfolder, s, h, dri,,, "Output folder for DRB files"
  22 +outfolder, s, a, dri,,, "Output folder for DRB files"
23 23 outobs, f, a, obs_back.xml,,, "Output observation definition file"
24 24  
25 25 #
... ...
modules/comscripts/comobsbin.par
  View file @aeaeaf4
... ... @@ -49,6 +49,7 @@ psdmin, i, h, 0,0,110, &quot;Minimum PSD value&quot;
49 49 psdmax, i, h, 110,0,110, "Maximum PSD value"
50 50 zetamin, r, h, 5.0,0.0,10.0, "Minimum Earth horizon - Phibar (zeta) angle (deg)"
51 51 fpmtflag, i, h, 0,0,2, "Handling of D2 modules with failed PMT flag (0: exclude, 1: include, 2: exclude PMT)"
  52 +timebin, r, h, 300.0,,, "Internal time bin for DRW computation (sec)"
52 53 d1use, s, h, 1111111,,, "D1 module usage (1: use, 0: don't use)"
53 54 d2use, s, h, 11111111111111,,, "D2 module usage (1: use, 0: don't use)"
54 55  
... ...
modules/comscripts/comobsbin.py
  View file @aeaeaf4
... ... @@ -86,6 +86,7 @@ class comobsbin(ctools.csobservation):
86 86 self['psdmax'].integer()
87 87 self['zetamin'].real()
88 88 self['fpmtflag'].integer()
  89 + self['timebin'].real()
89 90  
90 91 # Get D1 and D2 module usage strings
91 92 d1use = self['d1use'].string()
... ... @@ -510,7 +511,7 @@ class comobsbin(ctools.csobservation):
510 511 drwfile_global = gammalib.GFilename(drwname_global)
511 512  
512 513 # Write header
513   - self._log_header3(gammalib.NORMAL, 'Compute DRW for %.3f - %.3f MeV' % \
  514 + self._log_header3(gammalib.NORMAL, 'Check DRW for %.3f - %.3f MeV' % \
514 515 (ebounds.emin(i).MeV(), ebounds.emax(i).MeV()))
515 516  
516 517 # If DRW file exists in global datastore then do nothing
... ... @@ -534,14 +535,21 @@ class comobsbin(ctools.csobservation):
534 535 # Append filename to list
535 536 drwfiles.append(drwfile)
536 537  
  538 + # Log appending
  539 + self._log_value(gammalib.NORMAL, 'Append DRW for computation', drwfile.url())
  540 +
537 541 # Append DRW filename in output folder
538 542 drwnames.append(drwname)
539 543  
540 544 # If there are DRW files to compute then compute and save them now
541 545 if len(drwfiles) > 0:
542 546  
  547 + # Write header
  548 + self._log_header3(gammalib.NORMAL, 'Compute DRWs')
  549 +
543 550 # Compute DRWs
544   - drws.compute_drws(obs, self._select, self['zetamin'].real())
  551 + drws.compute_drws(obs, self._select, self['zetamin'].real(),
  552 + self['timebin'].real())
545 553  
546 554 # Phibar normalise DRWs to DREs
547 555 for i in range(ebounds.size()):
... ... @@ -572,7 +580,7 @@ class comobsbin(ctools.csobservation):
572 580 for i, drwfile in enumerate(drwfiles):
573 581  
574 582 # Write header
575   - self._log_header3(gammalib.NORMAL, 'Created DRW for %.3f - %.3f MeV' % \
  583 + self._log_header3(gammalib.NORMAL, 'Save DRW for %.3f - %.3f MeV' % \
576 584 (ebounds.emin(i).MeV(), ebounds.emax(i).MeV()))
577 585  
578 586 # Save DRW
... ...