Luigi Tibaldo / ctools

Commit 97cd3d1c4a43099898a5c6856fd043ab0d0b5339

Authored by Jürgen Knödlseder
1 parent 7f487530

Implement multiple Test sources in comlixmap

Inline Side-by-side
Showing 1 changed file with 53 additions and 36 deletions
modules/comscripts/comlixmap.py
  View file @97cd3d1
... ... @@ -41,6 +41,7 @@ class comlixmap(ctools.cslikelihood):
41 41 # Initialise members
42 42 self._maps = []
43 43 self._map_names = []
  44 + self._srcnames = []
44 45  
45 46 # Return
46 47 return
... ... @@ -58,11 +59,12 @@ class comlixmap(ctools.cslikelihood):
58 59 if self.obs().models().is_empty():
59 60 self.obs().models(self['inmodel'].filename())
60 61  
61   - # Get source name and check if source exists in model
62   - srcname = self['srcname'].string()
63   - if not self.obs().models().contains(srcname):
64   - msg = 'Source "%s" not found in models.' % srcname
65   - raise RuntimeError(msg)
  62 + # Get test source names and check if test sources exists in model
  63 + self._srcnames = self['srcname'].string().split(';')
  64 + for srcname in self._srcnames:
  65 + if not self.obs().models().contains(srcname):
  66 + msg = 'Source "%s" not found in models.' % srcname
  67 + raise RuntimeError(msg)
66 68  
67 69 # Query parameters
68 70 self['like_accuracy'].real()
... ... @@ -160,16 +162,20 @@ class comlixmap(ctools.cslikelihood):
160 162 nxpix = self['nxpix'].integer()
161 163 nypix = self['nypix'].integer()
162 164  
163   - # Compute number of free model parameters, exclusing 'RA' and 'DEC'
164   - nfree = 0
165   - source = self.obs().models()[self['srcname'].string()]
166   - for par in source:
167   - if par.name() != 'RA' and par.name() != 'DEC' and par.is_free():
168   - self._map_names.append(par.name())
169   - nfree += 1
  165 + # Initialise number of free parameters
  166 + nfree = 0
  167 +
  168 + # Compute number of free model parameters for all test sources,
  169 + # exclusing 'RA' and 'DEC'
  170 + for srcname in self._srcnames:
  171 + source = self.obs().models()[srcname]
  172 + for par in source:
  173 + if par.name() != 'RA' and par.name() != 'DEC' and par.is_free():
  174 + self._map_names.append(par.name())
  175 + nfree += 1
170 176  
171 177 # Initialise sky maps
172   - for i in range(nfree+1):
  178 + for i in range(nfree+len(self._srcnames)):
173 179 self._maps.append(gammalib.GSkyMap(proj, coordsys, xref, yref,
174 180 -binsz, binsz, nxpix, nypix))
175 181  
... ... @@ -192,15 +198,18 @@ class comlixmap(ctools.cslikelihood):
192 198 ts : float
193 199 TS of test source
194 200 """
195   - # Set and fix position of test source
196   - source = self.obs().models()[self['srcname'].string()]
197   - source['RA'].value(dir.ra_deg())
198   - source['DEC'].value(dir.dec_deg())
199   - source['RA'].fix()
200   - source['DEC'].fix()
  201 + # Loop over test source names
  202 + for srcname in self._srcnames:
201 203  
202   - # Remove response cache for test source
203   - self.obs().remove_response_cache(self['srcname'].string())
  204 + # Set and fix position of test source
  205 + source = self.obs().models()[srcname]
  206 + source['RA'].value(dir.ra_deg())
  207 + source['DEC'].value(dir.dec_deg())
  208 + source['RA'].fix()
  209 + source['DEC'].fix()
  210 +
  211 + # Remove response cache for test source
  212 + self.obs().remove_response_cache(srcname)
204 213  
205 214 # Get parameters and initialise some variables
206 215 niter = self['max_iter'].integer()
... ... @@ -240,7 +249,9 @@ class comlixmap(ctools.cslikelihood):
240 249 ts = self._final_model_fit()
241 250  
242 251 # Log TS
243   - self._log_value(gammalib.NORMAL, 'TS value', ts)
  252 + for i, srcname in enumerate(self._srcnames):
  253 + key = 'TS %s' % srcname
  254 + self._log_value(gammalib.NORMAL, key, ts[i])
244 255  
245 256 # Return TS
246 257 return ts
... ... @@ -261,8 +272,10 @@ class comlixmap(ctools.cslikelihood):
261 272 # Run ctlike
262 273 like.run()
263 274  
264   - # Recover TS for test source
265   - ts = like.obs().models()[self['srcname'].string()].ts()
  275 + # Recover TS for all test sources
  276 + ts = []
  277 + for srcname in self._srcnames:
  278 + ts.append(like.obs().models()[srcname].ts())
266 279  
267 280 # Return TS
268 281 return ts
... ... @@ -289,9 +302,9 @@ class comlixmap(ctools.cslikelihood):
289 302 # Log header
290 303 self._log_header1(gammalib.NORMAL, 'Initialise models')
291 304  
292   - # Make sure that tscalc flag is only set for test source
  305 + # Make sure that tscalc flag is only set for test sources
293 306 for source in self.obs().models():
294   - if source.name() == self['srcname'].string():
  307 + if source.name() in self._srcnames:
295 308 source.tscalc(True)
296 309 else:
297 310 source.tscalc(False)
... ... @@ -332,16 +345,18 @@ class comlixmap(ctools.cslikelihood):
332 345 # Compute TS
333 346 ts = self._compute_ts(dir)
334 347  
335   - # Store TS value in map
336   - self._maps[0][ipix] = ts
  348 + # Store TS values in sky map
  349 + for i, value in enumerate(ts):
  350 + self._maps[i][ipix] = value
337 351  
338   - # Extract fitted model parameters
339   - ipar = 1
340   - source = self.obs().models()[self['srcname'].string()]
341   - for par in source:
342   - if par.name() != 'RA' and par.name() != 'DEC' and par.is_free():
343   - self._maps[ipar][ipix] = par.value()
344   - ipar += 1
  352 + # Store fitted model parameters in sky maps
  353 + ipar = len(ts)
  354 + for srcname in self._srcnames:
  355 + source = self.obs().models()[srcname]
  356 + for par in source:
  357 + if par.name() != 'RA' and par.name() != 'DEC' and par.is_free():
  358 + self._maps[ipar][ipix] = par.value()
  359 + ipar += 1
345 360  
346 361 # Return
347 362 return
... ... @@ -376,8 +391,10 @@ class comlixmap(ctools.cslikelihood):
376 391 map.write(fits)
377 392  
378 393 # Set extension name for all maps
  394 + for i, srcname in enumerate(self._srcnames):
  395 + fits[i].extname(srcname)
379 396 for i, name in enumerate(self._map_names):
380   - fits[i+1].extname(name)
  397 + fits[i+len(self._srcnames)].extname(name)
381 398  
382 399 # Save FITS file
383 400 fits.saveto(outmap, self['clobber'].boolean())
... ...