ctools / ctools

Commit 1cf2fb6d972999179533fb110cfa7d286f2c980a

Authored by Jürgen Knödlseder
1 parent d3460e53

Implement computation to weighting cubes (DRW)

Inline Side-by-side
Showing 1 changed file with 111 additions and 11 deletions
modules/comscripts/comobsbin.py
  View file @1cf2fb6
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # Bin COMPTEL observations
4 4 #
5   -# Copyright (C) 2019-2022 Juergen Knoedlseder
  5 +# Copyright (C) 2019-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
... ... @@ -44,6 +44,7 @@ class comobsbin(ctools.csobservation):
44 44 self._select = gammalib.GCOMSelection()
45 45 self._drx_suffix = ''
46 46 self._drg_suffix = ''
  47 + self._drw_suffix = ''
47 48 self._dre_suffix = ''
48 49 self._global_datastore = '$COMDATA/datastore'
49 50  
... ... @@ -335,8 +336,8 @@ class comobsbin(ctools.csobservation):
335 336  
336 337 Returns
337 338 -------
338   - drxname, drgname, drenames : tuple of str
339   - Names of DRX, DRG and DREs
  339 + drxname, drgname, drwnames, drenames : tuple of str
  340 + Names of DRX, DRG, DRWs and DREs
340 341 """
341 342 # Get pointing direction
342 343 pnt = self._get_pointing_direction(obs)
... ... @@ -390,12 +391,17 @@ class comobsbin(ctools.csobservation):
390 391 # Create sky map for DRX
391 392 expo = gammalib.GSkyMap('CAR', 'GAL', 0.0, 0.0, 1.0, 1.0, 360, 180)
392 393  
393   - # Allocate DRE, DRG and DRX
  394 + # Allocate DRE, DRW, DRG and DRX
394 395 dre = gammalib.GCOMDri(cube, 0.0, dphibar, nphibar)
  396 + drw = gammalib.GCOMDri(cube, 0.0, dphibar, nphibar)
395 397 drg = gammalib.GCOMDri(cube, 0.0, dphibar, nphibar)
396 398 drx = gammalib.GCOMDri(expo)
397 399  
398   - # Set DRG and DRX filenames
  400 + # Allocate container for DRWs and initialise filename list
  401 + drws = gammalib.GCOMDris()
  402 + drwfiles = []
  403 +
  404 + # Set DRG and DRX filenames in output folder
399 405 drxname = '%s/%s_drx%s.fits' % (self['outfolder'].string(), obs.id(), self._drx_suffix)
400 406 drgname = '%s/%s%s_drg%s.fits' % (self['outfolder'].string(), obs.id(), dri_prefix, self._drg_suffix)
401 407 drxfile = gammalib.GFilename(drxname)
... ... @@ -439,7 +445,7 @@ class comobsbin(ctools.csobservation):
439 445 # Generate one DRE for each energy boundary
440 446 for i in range(ebounds.size()):
441 447  
442   - # Set DRE filename
  448 + # Set DRE filename in output folder
443 449 drename = '%s/%s%s_dre%s_%6.6d-%6.6dkeV.fits' % \
444 450 (self['outfolder'].string(), obs.id(), dri_prefix, self._dre_suffix,
445 451 ebounds.emin(i).keV(), ebounds.emax(i).keV())
... ... @@ -460,7 +466,7 @@ class comobsbin(ctools.csobservation):
460 466 drename = drename_global
461 467 self._log_value(gammalib.NORMAL, 'Global DRE file exists', drefile_global.url())
462 468  
463   - # If DRE file exists then do nothing
  469 + # If DRE file exists in output folder then do nothing
464 470 elif drefile.exists():
465 471 self._log_value(gammalib.NORMAL, 'Local DRE file exists', drefile.url())
466 472  
... ... @@ -482,11 +488,104 @@ class comobsbin(ctools.csobservation):
482 488 # Log DRE
483 489 self._log_string(gammalib.NORMAL, str(dre))
484 490  
485   - # Append DRE filename
  491 + # Append DRE filename in output folder
486 492 drenames.append(drename)
487 493  
  494 + # Initialse list of DRW names
  495 + drwnames = []
  496 +
  497 + # Generate one DRW for each energy boundary
  498 + for i in range(ebounds.size()):
  499 +
  500 + # Set DRW filename in output folder
  501 + drwname = '%s/%s%s_drw%s_%6.6d-%6.6dkeV.fits' % \
  502 + (self['outfolder'].string(), obs.id(), dri_prefix, self._drw_suffix,
  503 + ebounds.emin(i).keV(), ebounds.emax(i).keV())
  504 + drwfile = gammalib.GFilename(drwname)
  505 +
  506 + # Set DRW filename in global data store
  507 + drwname_global = '%s/%s%s_drw%s_%6.6d-%6.6dkeV.fits' % \
  508 + (self._global_datastore, obs.id(), dri_prefix, self._drw_suffix,
  509 + ebounds.emin(i).keV(), ebounds.emax(i).keV())
  510 + drwfile_global = gammalib.GFilename(drwname_global)
  511 +
  512 + # Write header
  513 + self._log_header3(gammalib.NORMAL, 'Compute DRW for %.3f - %.3f MeV' % \
  514 + (ebounds.emin(i).MeV(), ebounds.emax(i).MeV()))
  515 +
  516 + # If DRW file exists in global datastore then do nothing
  517 + if drwfile_global.exists():
  518 + drwname = drwname_global
  519 + self._log_value(gammalib.NORMAL, 'Global DRW file exists', drwfile_global.url())
  520 +
  521 + # If DRW file exists in output folder then do nothing
  522 + elif drwfile.exists():
  523 + self._log_value(gammalib.NORMAL, 'Local DRW file exists', drwfile.url())
  524 +
  525 + # ... otherwise append DRW to container
  526 + else:
  527 +
  528 + # Set DRW energy range
  529 + drw.ebounds(gammalib.GEbounds(ebounds.emin(i), ebounds.emax(i)))
  530 +
  531 + # Append DRW to container
  532 + drws.append(drw)
  533 +
  534 + # Append filename to list
  535 + drwfiles.append(drwfile)
  536 +
  537 + # Append DRW filename in output folder
  538 + drwnames.append(drwname)
  539 +
  540 + # If there are DRW files to compute then compute and save them now
  541 + if len(drwfiles) > 0:
  542 +
  543 + # Compute DRWs
  544 + drws.compute_drws(obs, self._select, self['zetamin'].real())
  545 +
  546 + # Phibar normalise DRWs to DREs
  547 + for i in range(ebounds.size()):
  548 +
  549 + # Load DRE
  550 + dre = gammalib.GCOMDri(drenames[i])
  551 +
  552 + # Get dataspace dimensions
  553 + nchi = dre.nchi()
  554 + npsi = dre.npsi()
  555 + nphibar = dre.nphibar()
  556 + npix = nchi * npsi
  557 +
  558 + # Phibar normalise DRW
  559 + for iphibar in range(nphibar):
  560 + sum_dre = 0.0
  561 + sum_drw = 0.0
  562 + for ipix in range(npix):
  563 + index = ipix + iphibar*npix
  564 + sum_dre += dre[index]
  565 + sum_drw += drws[i][index]
  566 + if sum_drw > 0:
  567 + for ipix in range(npix):
  568 + index = ipix + iphibar*npix
  569 + drws[i][index] *= sum_dre / sum_drw
  570 +
  571 + # Save DRWs
  572 + for i, drwfile in enumerate(drwfiles):
  573 +
  574 + # Write header
  575 + self._log_header3(gammalib.NORMAL, 'Created DRW for %.3f - %.3f MeV' % \
  576 + (ebounds.emin(i).MeV(), ebounds.emax(i).MeV()))
  577 +
  578 + # Save DRW
  579 + drws[i].save(drwfile, True)
  580 +
  581 + # Log creation
  582 + self._log_value(gammalib.NORMAL, 'Local DRW file created', drwfile.url())
  583 +
  584 + # Log DRW
  585 + self._log_string(gammalib.NORMAL, str(drw))
  586 +
488 587 # Return
489   - return drxname, drgname, drenames
  588 + return drxname, drgname, drwnames, drenames
490 589  
491 590 def _generate_drb(self, drename, drgname):
492 591 """
... ... @@ -743,7 +842,7 @@ class comobsbin(ctools.csobservation):
743 842 self._log_string(gammalib.NORMAL, str(obs))
744 843  
745 844 # Generate DRIs
746   - drxname, drgname, drenames = self._generate_dri(obs, ebounds)
  845 + drxname, drgname, drwnames, drenames = self._generate_dri(obs, ebounds)
747 846  
748 847 # Skip observation if no superpackets were used
749 848 drx = gammalib.GCOMDri(drxname)
... ... @@ -773,9 +872,10 @@ class comobsbin(ctools.csobservation):
773 872 xml_obs = xml_list.append('observation name="%s" id="%s" instrument="COM"' % \
774 873 (name, id))
775 874  
776   - # Append DRE, DRB, DRG and DRX files
  875 + # Append DRE, DRB, DRW, DRG and DRX files
777 876 xml_obs.append('parameter name="DRE" file="%s"' % (drename))
778 877 xml_obs.append('parameter name="DRB" file="%s"' % (drbname))
  878 + xml_obs.append('parameter name="DRW" file="%s"' % (drwnames[i]))
779 879 xml_obs.append('parameter name="DRG" file="%s"' % (drgname))
780 880 xml_obs.append('parameter name="DRX" file="%s"' % (drxname))
781 881  
... ...