ctools  2.0.0
 All Classes Namespaces Files Functions Variables Macros Pages
csebins.py
Go to the documentation of this file.
1 #! /usr/bin/env python
2 # ==========================================================================
3 # Generates energy boundaries for stacked analysis
4 #
5 # Copyright (C) 2017-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 math
23 import gammalib
24 import ctools
25 
26 
27 # ============= #
28 # csebins class #
29 # ============= #
30 class csebins(ctools.csobservation):
31  """
32  Generates energy boundaries for stacked analysis
33  """
34 
35  # Constructor
36  def __init__(self, *argv):
37  """
38  Constructor
39  """
40  # Initialise application by calling the appropriate class constructor
41  self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
42 
43  # Set members
44  self._ebounds = gammalib.GEbounds()
45 
46  # Return
47  return
48 
49 
50  # Private methods
51  def _get_parameters(self):
52  """
53  Get parameters from parfile
54  """
55  # Initialise observation container if it is currently empty
56  if self.obs().size() == 0:
57 
58  # If an observation definition file was provided then load it,
59  # otherwise build a single observation with response information
60  if self['inobs'].is_valid():
61  self.obs(gammalib.GObservations(self['inobs'].filename()))
62  else:
63  cta = gammalib.GCTAObservation()
64  mission = gammalib.tolower(self['instrument'].string())
65  caldb = gammalib.GCaldb(mission, self['caldb'].string())
66  rsp = gammalib.GCTAResponseIrf(self['irf'].string(), caldb)
67  cta.instrument(mission)
68  cta.response(rsp)
69  self.obs().append(cta)
70 
71  # Query input parameters
72  self['emin'].real()
73  self['emax'].real()
74  self['aeffthres'].real()
75  self['bkgthres'].real()
76 
77  # Query ahead output model filename
78  if self._read_ahead():
79  self['outfile'].filename()
80 
81  # Write input parameters into logger
82  self._log_parameters(gammalib.TERSE)
83 
84  # Return
85  return
86 
87  def _insert_energy(self, responses, energy, comment):
88  """
89  Set energy boundaries for one CTA observation
90 
91  Parameters
92  ----------
93  responses : list of dict
94  List of response dictionaries
95  energy : float
96  Energy to insert (TeV)
97  comment : str
98  Reason for energy insertion
99  """
100  # Set thresholds of all response components
101  logE = math.log10(energy)
102  for rsp in responses:
103  rsp['thres'] = rsp['irf'](logE, 0.0, 0.0)
104 
105  # Convert energy into gammalib energy
106  eng = gammalib.GEnergy(energy, 'TeV')
107 
108  # Set insertion message
109  msg = str(energy)+' TeV'
110  if len(comment) > 0:
111  msg += ' ('+comment+')'
112 
113  # If energy boundaries are empty then append an energy interval
114  # with zero size ...
115  if self._ebounds.size() == 0:
116 
117  # Append energy interval with zero size
118  self._ebounds.append(eng, eng)
119 
120  # Log first energy boundary
121  self._log_value(gammalib.NORMAL, 'First boundary', msg)
122 
123  # ... otherwise append an interval
124  elif self._ebounds.size() > 0:
125 
126  # If current boundaries have zero size then recover boundary energy,
127  # remove boundary, and append a boundary of non-zero size
128  if self._ebounds.emin() == self._ebounds.emax():
129 
130  # Recover old energy, remove boundary and append new boundary
131  emin = self._ebounds.emin().copy() # Because emin becomes zero after remove
132  emax = eng
133  self._ebounds.remove(0)
134  self._ebounds.append(emin, emax)
135 
136  # Log other energy boundary
137  self._log_value(gammalib.NORMAL, 'Boundary', msg)
138 
139  else:
140 
141  # Use old emax as new emin. Only append if interval is non zero
142  emin = self._ebounds.emax()
143  if emin < eng:
144  self._ebounds.append(emin, eng)
145 
146  # Log other energy boundary
147  self._log_value(gammalib.NORMAL, 'Boundary', msg)
148 
149  # Return
150  return
151 
152  def _set_ebounds(self, obs):
153  """
154  Set energy boundaries for one CTA observation
155 
156  Parameters
157  ----------
158  obs : `~gammalib.GObservations`
159  CTA observation
160  """
161  # Get parameters
162  emin = self['emin'].real()
163  emax = self['emax'].real()
164  aeffthres = self['aeffthres'].real()
165  bkgthres = self['bkgthres'].real()
166 
167  # Build list of response dictionaries
168  responses = []
169 
170  # Loop over all observations
171  for run in obs:
172 
173  # Get response name
174  mission = run.response().caldb().mission()
175  instrument = run.response().caldb().instrument()
176  rspname = run.response().rspname()
177  name = '%s::%s::%s' % (mission, instrument, rspname)
178 
179  # If response name exists already then examine next observation
180  exists = False
181  for rsp in responses:
182  if rsp['name'] == name:
183  exists = True
184  break
185  if exists:
186  continue
187 
188  # Append effective area and background template to list
189  aeff = {'name': name, 'type': 'Aeff',
190  'irf': run.response().aeff(), 'thres': 0.0}
191  bkg = {'name': name, 'type': 'Background',
192  'irf': run.response().background(), 'thres': 0.0}
193  responses.append(aeff)
194  responses.append(bkg)
195 
196  # Log response name
197  self._log_value(gammalib.NORMAL, 'Append response', '%s (%s) [%s]' %
198  (rspname, instrument, mission))
199 
200  # Setup energy vector in log10(TeV)
201  logEs = [math.log10(1.0e-3*float(i)) for i in range(int(emin*1000),int(emax*1000))]
202 
203  # Loop over all energies
204  for logE in logEs:
205 
206  # Loop over all response components
207  for rsp in responses:
208 
209  # Get IRF value
210  irf = rsp['irf'](logE, 0.0, 0.0)
211 
212  # Get threshold
213  if rsp['type'] == 'Aeff':
214  threshold = aeffthres
215  else:
216  threshold = bkgthres
217 
218  # If threshold is zero and response is positive then insert
219  # an energy ...
220  if rsp['thres'] == 0.0 and irf > 0.0:
221  self._insert_energy(responses, math.pow(10.0, logE), rsp['type'])
222  break
223 
224  # ... otherwise if threshold is not zero then check whether
225  # the fractional change exceeds the threshold
226  elif rsp['thres'] != 0.0:
227  f = irf/rsp['thres']
228  if 1.0-f > threshold or f-1.0 > threshold:
229  self._insert_energy(responses, math.pow(10.0, logE), rsp['type'])
230  break
231 
232  # Insert energy at end of range
233  self._insert_energy(responses, emax, 'End of range')
234 
235  # Log number of energy boundaries
236  self._log_value(gammalib.NORMAL, 'Number of boundaries', self._ebounds.size())
237 
238  # Return
239  return
240 
241 
242  # Public methods
243  def process(self):
244  """
245  Process the script
246  """
247  # Get parameters
248  self._get_parameters()
249 
250  # Write observation into logger
251  self._log_observations(gammalib.NORMAL, self.obs(), 'Input observation')
252 
253  # Write header
254  self._log_header1(gammalib.TERSE, 'Define energy boundaries')
255 
256  # Define energy boundaries for observation container
257  self._set_ebounds(self.obs())
258 
259  # Return
260  return
261 
262  def save(self):
263  """
264  Save energy boundaries
265  """
266  # Write header
267  self._log_header1(gammalib.TERSE, 'Save energy boundaries')
268 
269  # Get energy boundaries filename
270  outfile = self['outfile'].filename()
271 
272  # Log file name
273  self._log_value(gammalib.NORMAL, 'Energy boundaries file', outfile.url())
274 
275  # Save energy boundaries
276  self._ebounds.save(outfile, self._clobber())
277 
278  # Stamp energy boundaries
279  self._stamp(outfile)
280 
281  # Return
282  return
283 
284  def ebounds(self):
285  """
286  Return energy boundaries
287 
288  Returns
289  -------
290  ebounds : `~gammalib.GEbounds()`
291  Energy boundaries
292  """
293  # Return
294  return self._ebounds
295 
296 
297 # ======================== #
298 # Main routine entry point #
299 # ======================== #
300 if __name__ == '__main__':
301 
302  # Create instance of application
303  app = csebins(sys.argv)
304 
305  # Execute application
306  app.execute()