Commit 4768b0eae6668a8f03853805a52fa07f746ada2c
1 parent
7c9d9804
Add compulbin script (#3973)
Showing
9 changed files
with
518 additions
and
7 deletions
ChangeLog
NEWS
| 1 | 1 | New Features and Important Changes in ctools 2.0.0 |
| 2 | 2 | |
| 3 | -1 December 2021 | |
| 3 | +9 January 2022 | |
| 4 | 4 | |
| 5 | 5 | |
| 6 | 6 | Introduction |
| ... | ... | @@ -442,6 +442,11 @@ comobssim - Simulate COMPTEL observations |
| 442 | 442 | Added script (#3738). |
| 443 | 443 | |
| 444 | 444 | |
| 445 | +compulbin - Generate pulsar phase curve | |
| 446 | +--------------------------------------- | |
| 447 | +Added script (#3973). | |
| 448 | + | |
| 449 | + | |
| 445 | 450 | comsrcdetect - Detect source in Test Statistic map |
| 446 | 451 | -------------------------------------------------- |
| 447 | 452 | Added script (#3738). | ... | ... |
README.md
| 1 | 1 | ctools information |
| 2 | 2 | ================== |
| 3 | -* Version: 2.0.0.dev (1 December 2021) | |
| 3 | +* Version: 2.0.0.dev (9 January 2022) | |
| 4 | 4 | * GammaLib dependency: 2.0.0.dev |
| 5 | 5 | |
| 6 | 6 | [](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/) | ... | ... |
cscripts/cspull.py
| ... | ... | @@ -386,7 +386,7 @@ class cspull(ctools.csobservation): |
| 386 | 386 | table = gammalib.GFitsBinTable(nrows) |
| 387 | 387 | table.extname('PULL_DISTRIBUTION') |
| 388 | 388 | |
| 389 | - # Add keywors for compatibility with gammalib.GMWLSpectrum | |
| 389 | + # Add keywords for compatibility with gammalib.GMWLSpectrum | |
| 390 | 390 | table.card('INSTRUME', 'CTA', 'Name of Instrument') |
| 391 | 391 | table.card('TELESCOP', 'CTA', 'Name of Telescope') |
| 392 | 392 | ... | ... |
modules/comscripts/Makefile.am
| 1 | 1 | ############################################################################# |
| 2 | 2 | # Makefile for comscripts Python scripts # |
| 3 | 3 | # ------------------------------------------------------------------------- # |
| 4 | -# Copyright (C) 2021 Juergen Knoedlseder # | |
| 4 | +# Copyright (C) 2021-2022 Juergen Knoedlseder # | |
| 5 | 5 | # ------------------------------------------------------------------------- # |
| 6 | 6 | # # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify # |
| ... | ... | @@ -30,6 +30,7 @@ pythonscripts = $(top_builddir)/modules/comscripts/__init__.py \ |
| 30 | 30 | $(srcdir)/comobsres.py \ |
| 31 | 31 | $(srcdir)/comobsselect.py \ |
| 32 | 32 | $(srcdir)/comobssim.py \ |
| 33 | + $(srcdir)/compulbin.py \ | |
| 33 | 34 | $(srcdir)/comsrcdetect.py |
| 34 | 35 | |
| 35 | 36 | # Parfiles to be distributed in $(prefix)/syspfiles |
| ... | ... | @@ -42,6 +43,7 @@ parfiles = $(srcdir)/comlixfit.par \ |
| 42 | 43 | $(srcdir)/comobsres.par \ |
| 43 | 44 | $(srcdir)/comobsselect.par \ |
| 44 | 45 | $(srcdir)/comobssim.par \ |
| 46 | + $(srcdir)/compulbin.par \ | |
| 45 | 47 | $(srcdir)/comsrcdetect.par |
| 46 | 48 | |
| 47 | 49 | # Scripts for which symbolic links will be created in bindir |
| ... | ... | @@ -54,6 +56,7 @@ linkscripts = comlixfit \ |
| 54 | 56 | comobsres \ |
| 55 | 57 | comobsselect \ |
| 56 | 58 | comobssim \ |
| 59 | + compulbin \ | |
| 57 | 60 | comsrcdetect |
| 58 | 61 | |
| 59 | 62 | # Test scripts to be distributed in $(pyexecdir)/comscripts/tests | ... | ... |
modules/comscripts/__init__.py.in
| 1 | 1 | # ========================================================================== |
| 2 | 2 | # comscripts Python module |
| 3 | 3 | # |
| 4 | -# Copyright (C) 2021 Juergen Knoedlseder | |
| 4 | +# Copyright (C) 2021-2022 Juergen Knoedlseder | |
| 5 | 5 | # |
| 6 | 6 | # This program is free software: you can redistribute it and/or modify |
| 7 | 7 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -26,6 +26,7 @@ __all__ = ['comlixfit', |
| 26 | 26 | 'comobsmodel', |
| 27 | 27 | 'comobsres', |
| 28 | 28 | 'comobssim', |
| 29 | + 'compulbin', | |
| 29 | 30 | 'comsrcdetect'] |
| 30 | 31 | |
| 31 | 32 | # Import |
| ... | ... | @@ -38,6 +39,7 @@ from comscripts.comobsbin import comobsbin |
| 38 | 39 | from comscripts.comobsmodel import comobsmodel |
| 39 | 40 | from comscripts.comobsres import comobsres |
| 40 | 41 | from comscripts.comobssim import comobssim |
| 42 | +from comscripts.compulbin import compulbin | |
| 41 | 43 | from comscripts.comsrcdetect import comsrcdetect |
| 42 | 44 | |
| 43 | 45 | # Add test function | ... | ... |
modules/comscripts/comobsbin.py
| ... | ... | @@ -304,7 +304,7 @@ class comobsbin(ctools.csobservation): |
| 304 | 304 | self._dre_suffix += '_%s' % (d2use) |
| 305 | 305 | self._drg_suffix += '_%s' % (d2use) |
| 306 | 306 | |
| 307 | - # Log input observations | |
| 307 | + # Log selection set | |
| 308 | 308 | self._log_string(gammalib.NORMAL, str(self._select)) |
| 309 | 309 | |
| 310 | 310 | # Return | ... | ... |
modules/comscripts/compulbin.par
0 → 100644
| 1 | +############################################################################# | |
| 2 | +# # | |
| 3 | +# COMPTEL Science Analysis Tools # | |
| 4 | +# # | |
| 5 | +# ------------------------------------------------------------------------- # | |
| 6 | +# # | |
| 7 | +# File: compulbin.par # | |
| 8 | +# # | |
| 9 | +# Author: Juergen Knoedlseder # | |
| 10 | +# Jurgen.Knodlseder@irap.omp.eu # | |
| 11 | +# IRAP # | |
| 12 | +# # | |
| 13 | +# Purpose: Parameter file for the compulbin script # | |
| 14 | +# # | |
| 15 | +############################################################################# | |
| 16 | +# | |
| 17 | +# Input and output parameters | |
| 18 | +#============================ | |
| 19 | +inobs, f, a, obs.xml,,, "Input observation definition file" | |
| 20 | +outfile, f, a, pulsar.fits,,, "Output pulsar phase bins file" | |
| 21 | + | |
| 22 | +# | |
| 23 | +# Script parameters | |
| 24 | +#================== | |
| 25 | +psrname, s, a, PSR B0531+21,,, "Pulsar name" | |
| 26 | +ephemerides, f, a, ephem_crab.psrtime,,, "Pulsar ephemerides file" | |
| 27 | +emin, r, a, 0.75,0.05,50.0, "Minimum energy (MeV)" | |
| 28 | +emax, r, a, 30.0,0.05,50.0, "Maximum energy (MeV)" | |
| 29 | +pnumbins, i, a, 33,1,200, "Number of phase bins" | |
| 30 | +armmin, r, h, -3.0,,, "Minimum ARM (deg)" | |
| 31 | +armmax, r, h, 3.0,,, "Maximum ARM (deg)" | |
| 32 | +phimin, r, h, 0.0,,, "Minimum Phibar angle (deg)" | |
| 33 | +phimax, r, h, 50.0,,, "Maximum Phibar angle (deg)" | |
| 34 | +psdmin, i, h, 0,0,110, "Minimum PSD value" | |
| 35 | +psdmax, i, h, 110,0,110, "Maximum PSD value" | |
| 36 | +zetamin, r, h, 5.0,0.0,10.0, "Minimum Earth horizon - Phibar (zeta) angle (deg)" | |
| 37 | +fpmtflag, i, h, 0,0,2, "Handling of D2 modules with failed PMT flag (0: exclude, 1: include, 2: exclude PMT)" | |
| 38 | +d1use, s, h, 1111111,,, "D1 module usage (1: use, 0: don't use)" | |
| 39 | +d2use, s, h, 11111111111111,,, "D2 module usage (1: use, 0: don't use)" | |
| 40 | + | |
| 41 | +# | |
| 42 | +# Standard parameters | |
| 43 | +#==================== | |
| 44 | +chatter, i, h, 2,0,4, "Chattiness of output" | |
| 45 | +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" | |
| 46 | +debug, b, h, no,,, "Activate debugging mode?" | |
| 47 | +mode, s, h, ql,,, "Mode of automatic parameters" | |
| 48 | +logfile, f, h, compulbin.log,,, "Log filename" | ... | ... |
modules/comscripts/compulbin.py
0 → 100755
| 1 | +#! /usr/bin/env python | |
| 2 | +# ========================================================================== | |
| 3 | +# Generate pulsar phase bins | |
| 4 | +# | |
| 5 | +# Copyright (C) 2022 Juergen Knoedlseder | |
| 6 | +# | |
| 7 | +# This program is free software: you can redistribute it and/or modify | |
| 8 | +# it under the terms of the GNU General Public License as published by | |
| 9 | +# the Free Software Foundation, either version 3 of the License, or | |
| 10 | +# (at your option) any later version. | |
| 11 | +# | |
| 12 | +# This program is distributed in the hope that it will be useful, | |
| 13 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 14 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| 15 | +# GNU General Public License for more details. | |
| 16 | +# | |
| 17 | +# You should have received a copy of the GNU General Public License | |
| 18 | +# along with this program. If not, see <http://www.gnu.org/licenses/>. | |
| 19 | +# | |
| 20 | +# ========================================================================== | |
| 21 | +import sys | |
| 22 | +import glob | |
| 23 | +import os | |
| 24 | +import gammalib | |
| 25 | +import ctools | |
| 26 | + | |
| 27 | + | |
| 28 | +# =============== # | |
| 29 | +# compulbin class # | |
| 30 | +# =============== # | |
| 31 | +class compulbin(ctools.csobservation): | |
| 32 | + """ | |
| 33 | + Generate pulsar phase bins | |
| 34 | + """ | |
| 35 | + # Constructor | |
| 36 | + def __init__(self, *argv): | |
| 37 | + """ | |
| 38 | + Constructor | |
| 39 | + """ | |
| 40 | + # Initialise application by calling the base class constructor | |
| 41 | + self._init_csobservation(self.__class__.__name__, ctools.__version__, argv) | |
| 42 | + | |
| 43 | + # Initialise members | |
| 44 | + self._select = gammalib.GCOMSelection() | |
| 45 | + self._pnumbins = 0 | |
| 46 | + self._phasebin = 1.0 | |
| 47 | + self._bins = [] | |
| 48 | + self._fits = None | |
| 49 | + | |
| 50 | + # Return | |
| 51 | + return | |
| 52 | + | |
| 53 | + # Private methods | |
| 54 | + def _get_parameters(self): | |
| 55 | + """ | |
| 56 | + Get parameters from parfile | |
| 57 | + """ | |
| 58 | + # Set observation if not done before | |
| 59 | + if self.obs().is_empty(): | |
| 60 | + self.obs().load(self['inobs'].filename()) | |
| 61 | + | |
| 62 | + # Query parameters | |
| 63 | + self['psrname'].string() | |
| 64 | + self['ephemerides'].filename() | |
| 65 | + self['emin'].real() | |
| 66 | + self['emax'].real() | |
| 67 | + self._pnumbins = self['pnumbins'].integer() | |
| 68 | + self['armmin'].real() | |
| 69 | + self['armmax'].real() | |
| 70 | + self['phimin'].real() | |
| 71 | + self['phimax'].real() | |
| 72 | + self['zetamin'].real() | |
| 73 | + self['fpmtflag'].integer() | |
| 74 | + self['psdmin'].integer() | |
| 75 | + self['psdmax'].integer() | |
| 76 | + | |
| 77 | + # Get D1 and D2 module usage strings | |
| 78 | + d1use = self['d1use'].string() | |
| 79 | + d2use = self['d2use'].string() | |
| 80 | + | |
| 81 | + # Check D1 and D2 module usage strings | |
| 82 | + if len(d1use) != 7: | |
| 83 | + msg = 'Incorrect length %d of D1 usage string. String needs to have 7 digits.' % \ | |
| 84 | + len(d1use) | |
| 85 | + raise RuntimeError(msg) | |
| 86 | + if len(d2use) != 14: | |
| 87 | + msg = 'Incorrect length %d of D2 usage string. String needs to have 14 digits.' % \ | |
| 88 | + len(d2use) | |
| 89 | + raise RuntimeError(msg) | |
| 90 | + | |
| 91 | + # Query ahead output model filename | |
| 92 | + if self._read_ahead(): | |
| 93 | + self['outfile'].filename() | |
| 94 | + | |
| 95 | + # Write input parameters into logger | |
| 96 | + self._log_parameters(gammalib.TERSE) | |
| 97 | + | |
| 98 | + # Return | |
| 99 | + return | |
| 100 | + | |
| 101 | + def _set_selection_set(self): | |
| 102 | + """ | |
| 103 | + Set selection set | |
| 104 | + """ | |
| 105 | + # Write header | |
| 106 | + self._log_header1(gammalib.NORMAL, 'Set selection set') | |
| 107 | + | |
| 108 | + # Clear selection set | |
| 109 | + self._select.clear() | |
| 110 | + | |
| 111 | + # Get pulsar ephemerides | |
| 112 | + pulsar = gammalib.GPulsar(self['ephemerides'].filename(), self['psrname'].string()) | |
| 113 | + | |
| 114 | + # Set pulsar ephemerides | |
| 115 | + self._select.pulsar(pulsar) | |
| 116 | + | |
| 117 | + # Set PSD interval | |
| 118 | + self._select.psd_min(self['psdmin'].integer()) | |
| 119 | + self._select.psd_max(self['psdmax'].integer()) | |
| 120 | + | |
| 121 | + # Set handling of D2 modules with failed PMT flag | |
| 122 | + self._select.fpmtflag(self['fpmtflag'].integer()) | |
| 123 | + | |
| 124 | + # Set D1 module usage flags | |
| 125 | + d1use = self['d1use'].string() | |
| 126 | + for i in range(7): | |
| 127 | + if d1use[i] == '1': | |
| 128 | + self._select.use_d1(i,True) | |
| 129 | + else: | |
| 130 | + self._select.use_d1(i,False) | |
| 131 | + | |
| 132 | + # Set D2 module usage flags | |
| 133 | + d2use = self['d2use'].string() | |
| 134 | + for i in range(14): | |
| 135 | + if d2use[i] == '1': | |
| 136 | + self._select.use_d2(i,True) | |
| 137 | + else: | |
| 138 | + self._select.use_d2(i,False) | |
| 139 | + | |
| 140 | + # Log selection set | |
| 141 | + self._log_string(gammalib.NORMAL, str(self._select)) | |
| 142 | + | |
| 143 | + # Return | |
| 144 | + return | |
| 145 | + | |
| 146 | + def _bin_observation(self, obs, zetamin=5.0): | |
| 147 | + """ | |
| 148 | + Bin events in phase bins | |
| 149 | + | |
| 150 | + Parameters | |
| 151 | + ---------- | |
| 152 | + obs : `~gammalib.GCOMObservation` | |
| 153 | + Unbinned COMPTEL observation | |
| 154 | + zetamin : float, optional | |
| 155 | + Minimum zeta angle | |
| 156 | + """ | |
| 157 | + # Initialise event index | |
| 158 | + ievent = 0 | |
| 159 | + | |
| 160 | + # Initialise statistics | |
| 161 | + num_used_superpackets = 0 | |
| 162 | + num_skipped_superpackets = 0 | |
| 163 | + num_processed = 0 | |
| 164 | + num_event_outside_sp = 0 | |
| 165 | + num_energy_too_low = 0 | |
| 166 | + num_energy_too_high = 0 | |
| 167 | + num_phibar_too_small = 0 | |
| 168 | + num_phibar_too_large = 0 | |
| 169 | + num_eha_too_small = 0 | |
| 170 | + num_arm_too_small = 0 | |
| 171 | + num_arm_too_large = 0 | |
| 172 | + num_bin_too_small = 0 | |
| 173 | + num_bin_too_large = 0 | |
| 174 | + num_used_events = 0 | |
| 175 | + | |
| 176 | + # Get selection parameters | |
| 177 | + emin = gammalib.GEnergy(self['emin'].real(), 'MeV') | |
| 178 | + emax = gammalib.GEnergy(self['emax'].real(), 'MeV') | |
| 179 | + armmin = self['armmin'].real() | |
| 180 | + armmax = self['armmax'].real() | |
| 181 | + phimin = self['phimin'].real() | |
| 182 | + phimax = self['phimax'].real() | |
| 183 | + | |
| 184 | + # Get Good Time Intervals and reduce them to the validity range | |
| 185 | + # of the ephemerides | |
| 186 | + tim = obs.tim() | |
| 187 | + tim.reduce(self._select.pulsar().validity()) | |
| 188 | + | |
| 189 | + # Loop over Orbit Aspect Data | |
| 190 | + for oad in obs.oads(): | |
| 191 | + | |
| 192 | + # Skip superpacket if it is not fully enclosed within the COMPTEL | |
| 193 | + # Good Time Interval | |
| 194 | + if not tim.contains(oad.tstart()) or not tim.contains(oad.tstop()): | |
| 195 | + num_skipped_superpackets += 1 | |
| 196 | + continue | |
| 197 | + | |
| 198 | + # Increment superpacket counter | |
| 199 | + num_used_superpackets += 1 | |
| 200 | + | |
| 201 | + # Collect all events in the superpacket | |
| 202 | + while ievent < obs.events().size(): | |
| 203 | + | |
| 204 | + # Get event | |
| 205 | + event = obs.events()[ievent] | |
| 206 | + | |
| 207 | + # Break event loop if end of superpacket was reached | |
| 208 | + if event.time() > oad.tstop(): | |
| 209 | + break | |
| 210 | + | |
| 211 | + # Increment event counter | |
| 212 | + ievent += 1 | |
| 213 | + | |
| 214 | + # Skip event if it lies before the superpacket start | |
| 215 | + if event.time() < oad.tstart(): | |
| 216 | + num_event_outside_sp += 1 | |
| 217 | + continue | |
| 218 | + | |
| 219 | + # Increment number of processed events | |
| 220 | + num_processed += 1 | |
| 221 | + | |
| 222 | + # Skip event if it lies outside energy range | |
| 223 | + if event.energy() < emin: | |
| 224 | + num_energy_too_low += 1 | |
| 225 | + continue | |
| 226 | + elif event.energy() > emax: | |
| 227 | + num_energy_too_high += 1 | |
| 228 | + continue | |
| 229 | + | |
| 230 | + # Apply event selection | |
| 231 | + if not self._select.use_event(event): | |
| 232 | + continue | |
| 233 | + | |
| 234 | + # Skip event if it lies outside Phibar range | |
| 235 | + if event.phibar() < phimin: | |
| 236 | + num_phibar_too_small += 1 | |
| 237 | + continue | |
| 238 | + elif event.phibar() > phimax: | |
| 239 | + num_phibar_too_large += 1 | |
| 240 | + continue | |
| 241 | + | |
| 242 | + # Skip event if it comes from the Earth horizon | |
| 243 | + ehamin = event.phibar() + zetamin | |
| 244 | + if event.eha() < ehamin: | |
| 245 | + num_eha_too_small += 1 | |
| 246 | + continue | |
| 247 | + | |
| 248 | + # Get pulsar ephemerides for event time | |
| 249 | + ephemerides = self._select.pulsar().ephemerides(event.time()) | |
| 250 | + | |
| 251 | + # Compute Phigeo | |
| 252 | + phigeo = event.dir().dir().dist_deg(ephemerides.dir()) | |
| 253 | + | |
| 254 | + # Make ARM selection | |
| 255 | + arm = event.phibar() - phigeo | |
| 256 | + if arm < armmin: | |
| 257 | + num_arm_too_small += 1 | |
| 258 | + continue | |
| 259 | + elif arm > armmax: | |
| 260 | + num_arm_too_large += 1 | |
| 261 | + continue | |
| 262 | + | |
| 263 | + # Convert event time to Solar System Barycentre time | |
| 264 | + tdelta = obs.bvcs().tdelta(ephemerides.dir(), event.time()) | |
| 265 | + time = event.time() + tdelta | |
| 266 | + | |
| 267 | + # Compute pulsar phase | |
| 268 | + phase = ephemerides.phase(time) | |
| 269 | + | |
| 270 | + # Fill event in phase bin | |
| 271 | + ibin = int(phase / self._phasebin) | |
| 272 | + if ibin < 0: | |
| 273 | + num_bin_too_small += 1 | |
| 274 | + elif ibin >= self._pnumbins: | |
| 275 | + num_bin_too_large += 1 | |
| 276 | + else: | |
| 277 | + self._bins[ibin] += 1.0 | |
| 278 | + num_used_events += 1 | |
| 279 | + | |
| 280 | + # Break superpacket loop if there are no more events | |
| 281 | + if ievent >= obs.events().size(): | |
| 282 | + break | |
| 283 | + | |
| 284 | + # If not all events are exhausted then add them to the events outside | |
| 285 | + # superpackets | |
| 286 | + if ievent < obs.events().size(): | |
| 287 | + num_event_outside_sp += (obs.events().size() - ievent) | |
| 288 | + | |
| 289 | + # Log statistics | |
| 290 | + self._log_header3(gammalib.NORMAL, 'Event binning statistics') | |
| 291 | + self._log_value(gammalib.NORMAL, 'Number of superpackets', obs.oads().size()) | |
| 292 | + self._log_value(gammalib.NORMAL, 'Used superpackets', num_used_superpackets) | |
| 293 | + self._log_value(gammalib.NORMAL, 'Skipped superpackets', num_skipped_superpackets) | |
| 294 | + self._log_value(gammalib.NORMAL, 'Total number of events', obs.events().size()) | |
| 295 | + self._log_value(gammalib.NORMAL, 'Processed events', num_processed) | |
| 296 | + self._log_value(gammalib.NORMAL, 'Used events', num_used_events) | |
| 297 | + self._log_value(gammalib.NORMAL, 'Events outside superpacket', num_event_outside_sp) | |
| 298 | + self._log_value(gammalib.NORMAL, 'Energy too low', num_energy_too_low) | |
| 299 | + self._log_value(gammalib.NORMAL, 'Energy too high', num_energy_too_high) | |
| 300 | + self._log_value(gammalib.NORMAL, 'Phibar too small', num_phibar_too_small) | |
| 301 | + self._log_value(gammalib.NORMAL, 'Phibar too large', num_phibar_too_large) | |
| 302 | + self._log_value(gammalib.NORMAL, 'Earth hor. angle too small', num_eha_too_small) | |
| 303 | + self._log_value(gammalib.NORMAL, 'ARM too small', num_arm_too_small) | |
| 304 | + self._log_value(gammalib.NORMAL, 'ARM too large', num_arm_too_large) | |
| 305 | + self._log_value(gammalib.NORMAL, 'Phase bin index too small', num_bin_too_small) | |
| 306 | + self._log_value(gammalib.NORMAL, 'Phase bin index too large', num_bin_too_large) | |
| 307 | + | |
| 308 | + # Return | |
| 309 | + return | |
| 310 | + | |
| 311 | + def _create_fits(self): | |
| 312 | + """ | |
| 313 | + Creates FITS file with phase bins | |
| 314 | + """ | |
| 315 | + # Create FITS table columns | |
| 316 | + phase_min = gammalib.GFitsTableFloatCol('PHASE_MIN', self._pnumbins) | |
| 317 | + phase_max = gammalib.GFitsTableFloatCol('PHASE_MAX', self._pnumbins) | |
| 318 | + nevents = gammalib.GFitsTableDoubleCol('NEVENTS', self._pnumbins) | |
| 319 | + nevents.unit('counts') | |
| 320 | + | |
| 321 | + # Fill FITS table columns | |
| 322 | + for i in range(self._pnumbins): | |
| 323 | + phase_min[i] = i * self._phasebin | |
| 324 | + phase_max[i] = (i+1) * self._phasebin | |
| 325 | + nevents[i] = self._bins[i] | |
| 326 | + | |
| 327 | + # Create FITS Table with extension "PHASECURVE" | |
| 328 | + table = gammalib.GFitsBinTable(self._pnumbins) | |
| 329 | + table.extname('PHASECURVE') | |
| 330 | + | |
| 331 | + # Add keywords | |
| 332 | + table.card('INSTRUME', 'CGRO', 'Name of Instrument') | |
| 333 | + table.card('TELESCOP', 'COMPTEL', 'Name of Telescope') | |
| 334 | + | |
| 335 | + # Stamp header | |
| 336 | + self._stamp(table) | |
| 337 | + | |
| 338 | + # Add script keywords | |
| 339 | + table.card('PSRNAME', self['psrname'].string(), 'Pulsar name') | |
| 340 | + table.card('EPHEM', self['ephemerides'].filename().url(), 'Ephemerides file') | |
| 341 | + table.card('EMIN', self['emin'].real(), '[MeV] Minimum energy') | |
| 342 | + table.card('EMAX', self['emax'].real(), '[MeV] Maximum energy') | |
| 343 | + table.card('ARMMIN', self['armmin'].real(), '[deg] Minimum angular resolution measure') | |
| 344 | + table.card('ARMMAX', self['armmax'].real(), '[deg] Maximum angular resolution measure') | |
| 345 | + table.card('PHIMIN', self['phimin'].real(), '[deg] Minimum Phibar') | |
| 346 | + table.card('PHIMAX', self['phimax'].real(), '[deg] Maximum Phibar') | |
| 347 | + table.card('ZETAMIN', self['zetamin'].real(), '[deg] Minimum zeta angle') | |
| 348 | + table.card('FPMTFLAG', self['fpmtflag'].integer(), 'Handling of failed PMTs') | |
| 349 | + table.card('PSDMIN', self['psdmin'].integer(), 'Minimum PSD channel') | |
| 350 | + table.card('PSDMAX', self['psdmax'].integer(), 'Maximum PSD channel') | |
| 351 | + table.card('D1USE', self['d1use'].string(), 'D1 modules usage') | |
| 352 | + table.card('D2USE', self['d2use'].string(), 'D2 modules usage') | |
| 353 | + | |
| 354 | + # Append filled columns to fits table | |
| 355 | + table.append(phase_min) | |
| 356 | + table.append(phase_max) | |
| 357 | + table.append(nevents) | |
| 358 | + | |
| 359 | + # Create the FITS file now | |
| 360 | + self._fits = gammalib.GFits() | |
| 361 | + self._fits.append(table) | |
| 362 | + | |
| 363 | + # Return table | |
| 364 | + return | |
| 365 | + | |
| 366 | + | |
| 367 | + # Public methods | |
| 368 | + def run(self): | |
| 369 | + """ | |
| 370 | + Run the script | |
| 371 | + """ | |
| 372 | + # Switch screen logging on in debug mode | |
| 373 | + if self._logDebug(): | |
| 374 | + self._log.cout(True) | |
| 375 | + | |
| 376 | + # Get parameters | |
| 377 | + self._get_parameters() | |
| 378 | + | |
| 379 | + # Initialise phase bins | |
| 380 | + self._phasebin = 1.0 / self._pnumbins | |
| 381 | + self._bins = [0.0 for i in range(self._pnumbins)] | |
| 382 | + | |
| 383 | + # Log header | |
| 384 | + self._log_header1(gammalib.NORMAL, 'Input observations') | |
| 385 | + | |
| 386 | + # Log input observations | |
| 387 | + self._log_string(gammalib.NORMAL, str(self.obs())) | |
| 388 | + | |
| 389 | + # Set selection set | |
| 390 | + self._set_selection_set() | |
| 391 | + | |
| 392 | + # Write header | |
| 393 | + self._log_header1(gammalib.NORMAL, 'Generate pulsar phase bins') | |
| 394 | + | |
| 395 | + # Loop over all input observations | |
| 396 | + for obs in self.obs(): | |
| 397 | + | |
| 398 | + # Write header | |
| 399 | + self._log_header2(gammalib.NORMAL, self._get_obs_header(obs)) | |
| 400 | + | |
| 401 | + # Log observation | |
| 402 | + self._log_string(gammalib.NORMAL, str(obs)) | |
| 403 | + | |
| 404 | + # Generate phase bins | |
| 405 | + self._bin_observation(obs) | |
| 406 | + | |
| 407 | + # Create FITS file | |
| 408 | + self._create_fits() | |
| 409 | + | |
| 410 | + # Return | |
| 411 | + return | |
| 412 | + | |
| 413 | + def save(self): | |
| 414 | + """ | |
| 415 | + Save phase bin file | |
| 416 | + """ | |
| 417 | + # Write header | |
| 418 | + self._log_header1(gammalib.TERSE, 'Save phase bin file') | |
| 419 | + | |
| 420 | + # Continue only if FITS file is valid | |
| 421 | + if self._fits != None: | |
| 422 | + | |
| 423 | + # Get output filename | |
| 424 | + outfile = self['outfile'].filename() | |
| 425 | + | |
| 426 | + # If file exists and clobber flag is false then raise an exception | |
| 427 | + if outfile.exists() and not self['clobber'].boolean(): | |
| 428 | + msg = ('Cannot save "'+outfile.url()+'": File already exists. ' | |
| 429 | + 'Use parameter clobber=yes to allow overwriting of files.') | |
| 430 | + raise RuntimeError(msg) | |
| 431 | + | |
| 432 | + # ... otherwise log filename and save file | |
| 433 | + else: | |
| 434 | + # Log filename | |
| 435 | + self._log_value(gammalib.NORMAL, 'Phase bin file', outfile.url()) | |
| 436 | + | |
| 437 | + # Save phase bins | |
| 438 | + self._fits.saveto(outfile, self['clobber'].boolean()) | |
| 439 | + | |
| 440 | + # Return | |
| 441 | + return | |
| 442 | + | |
| 443 | + | |
| 444 | +# ======================== # | |
| 445 | +# Main routine entry point # | |
| 446 | +# ======================== # | |
| 447 | +if __name__ == '__main__': | |
| 448 | + | |
| 449 | + # Create instance of application | |
| 450 | + app = compulbin(sys.argv) | |
| 451 | + | |
| 452 | + # Execute application | |
| 453 | + app.execute() | ... | ... |