24 from cscripts
import obsutils
25 from cscripts
import ioutils
26 from cscripts
import mputils
34 Generates spectra in phase bins
36 This script computes spectra by performing a maximum likelihood fit
37 using :doc:`ctlike` in a series of phase bins for pulsars. The phase bins
38 can be either specified in an ASCII file or as an interval divided into
39 equally sized phase bins.
42 >>> phcrv = csphasecrv()
44 >>> ... (querying for parameters) ...
45 >>> phcrv = phrv.phasecrv()
46 Generates phase fits and retrieves dictionary with best fit models.
48 >>> phcrv = csphasecrv(obs)
50 >>> ... (querying for parameters) ...
51 Generates phase fits from the observations
52 and saves results in XML files.
61 self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
80 Extend ctools.csobservation getstate method to include some members
88 state = {
'base' : ctools.csobservation.__getstate__(self),
103 Extend ctools.csobservation setstate method to include some members
110 ctools.csobservation.__setstate__(self, state[
'base'])
113 self.
_onoff = state[
'onoff']
115 self.
_fits = state[
'fits']
126 Get parameters from parfile
130 self._setup_observations(self.obs(),
True,
True,
False)
133 self._set_obs_statistic(gammalib.toupper(self[
'statistic'].string()))
136 if self.obs().models().size() == 0:
137 self.obs().models(self[
'inmodel'].filename())
143 self.
_onoff = self._is_onoff()
146 self.
_srcname = self[
'srcname'].string()
154 self[
'edisp'].boolean()
157 if self._read_ahead():
158 self[
'outfile'].filename()
161 self._log_parameters(gammalib.TERSE)
184 algorithm = self[
'phbinalg'].string()
188 if algorithm ==
'FILE':
191 filename = self[
'phbinfile'].filename()
197 csv = gammalib.GCsv(filename)
198 for i
in range(csv.nrows()):
199 phmin = csv.real(i,0)
200 phmax = csv.real(i,1)
201 phbins.append([phmin,phmax])
205 elif algorithm ==
'LIN':
208 nbins = self[
'phbins'].integer()
211 ph_step = 1.0 / float(nbins)
212 for i
in range(nbins):
214 phmax = (i+1) * ph_step
215 phbins.append([phmin,phmax])
222 Return list of free parameter names
227 List of free parameter names.
233 for par
in self.obs().models()[self.
_srcname]:
235 names.append(par.name())
242 Creates FITS binary table containing phase curve results
246 results : list of dict
247 List of result dictionaries
251 table : `~gammalib.GFitsBinTable`
252 FITS binary table containing phase curve
258 table = gammalib.GFitsBinTable(nrows)
259 table.extname(
'PHASECURVE')
262 ioutils.append_result_column(table, results,
'PHASE_MIN',
'',
'phmin')
263 ioutils.append_result_column(table, results,
'PHASE_MAX',
'',
'phmax')
266 ioutils.append_model_par_column(table, self.obs().models()[self.
_srcname],
274 Create FITS file object from fit results
283 for i
in range(len(self.
_phbins)):
290 result = {
'phmin': phmin,
296 phname = str(phmin)+
'-'+str(phmax)
302 result[
'values'][par] = source[par].value()
303 result[
'values'][
'e_'+par] = source[par].error()
308 result[
'values'][par] = 0.
309 result[
'values'][
'e_'+par] = 0.
312 results.append(result)
318 self.
_fits = gammalib.GFits()
319 self._fits.append(table)
322 self._stamp(self.
_fits)
329 Saved phase fit results into a fits file
332 self._log_header1(gammalib.TERSE,
'Save phase FITS file')
338 outfile = self[
'outfile'].filename()
341 self._log_value(gammalib.NORMAL,
'Phase curve file', outfile.url())
344 self._fits.saveto(outfile, self._clobber())
351 Save phase fit results into one XML model each bin
354 self._log_header1(gammalib.TERSE,
'Save phase XML files')
357 outname = self[
'outfile'].filename()
363 phname = str(phbin[0])+
'-'+str(phbin[1])
364 replace_name =
'_phase_'+phname+
'.xml'
365 outname_phase = (outname.file()).replace(
'.fits', replace_name)
366 outfile = outname.path() + outname_phase
372 name =
'Phase [%s] file' % phname
373 self._log_value(gammalib.NORMAL, name, outfile)
380 Analyse one phase bin
384 phbin : list of float
389 self._log_header2(gammalib.TERSE,
'PHASE %f - %f' %
390 (phbin[0], phbin[1]))
393 select = ctools.ctselect(self.obs())
394 select[
'emin'] = self[
'emin'].real()
395 select[
'emax'] = self[
'emax'].real()
396 select[
'tmin'] =
'UNDEFINED'
397 select[
'tmax'] =
'UNDEFINED'
398 select[
'rad'] =
'UNDEFINED'
399 select[
'ra'] =
'UNDEFINED'
400 select[
'dec'] =
'UNDEFINED'
401 select[
'expr'] =
'PHASE>' + str(phbin[0]) +
' && PHASE<' + str(phbin[1])
405 phstr = str(phbin[0]) +
'-' + str(phbin[1])
408 for i
in range(0, select.obs().size()):
409 oldid = select.obs()[i].id()
410 select.obs()[i].id(oldid +
'_' + phstr)
415 obs = obsutils.get_onoff_obs(self, select.obs(), nthreads=1)
421 obs = obsutils.get_stacked_obs(self, select.obs())
424 self._log_header3(gammalib.EXPLICIT,
'Fitting the data')
433 like = ctools.ctlike(obs)
434 like[
'edisp'] = self[
'edisp'].boolean()
440 for model
in like.obs().models():
441 scaled_norm = model[
'Normalization'].value() / (phbin[1] - phbin[0])
442 model[
'Normalization'].value(scaled_norm)
445 fitmodels = like.obs().models().copy()
449 self._log_string(gammalib.TERSE,
450 'PHASE %f - %f: no observations available'
451 ' for fitting' % (phbin[0], phbin[1]))
454 fitmodels = gammalib.GModels()
457 result = {
'phstr': phstr,
'fitmodels': fitmodels}
472 self._log_observations(gammalib.NORMAL, self.obs(),
'Observation')
478 self._log_header1(gammalib.TERSE,
'Generate phase curve')
484 args = [(self,
'_phase_bin', phbin)
for phbin
in self.
_phbins]
485 poolresults = mputils.process(self.
_nthreads, mputils.mpfunc, args)
488 for i
in range(len(self.
_phbins)):
489 self.
_fitmodels[poolresults[i][0][
'phstr']] = poolresults[i][0][
'fitmodels']
490 self._log_string(gammalib.TERSE, poolresults[i][1][
'log'],
False)
496 self.
_fitmodels[result[
'phstr']] = result[
'fitmodels']
502 if self[
'publish'].boolean():
510 Saves results to fits and XML
531 self._log_header1(gammalib.TERSE,
'Publish phase curve')
534 if self._fits.contains(
'PHASECURVE'):
538 user_name = self._name()
543 self._log_value(gammalib.NORMAL,
'Phase curve name', user_name)
546 self._fits.publish(
'PHASECURVE', user_name)
553 Return dictionary with best fit models
560 Return and optionally set the exclusion regions map
564 object : `~gammalib.GSkyRegion` or `~gammalib.GSkyMap` or `~gammalib.GFilename`
565 Exclusion regions object
569 region : `~gammalib.GSkyRegionMap`
570 Exclusion regions map
573 if object
is not None:
583 if __name__ ==
'__main__':