24 from cscripts
import obsutils
25 from cscripts
import modutils
26 from cscripts
import ioutils
27 from cscripts
import mputils
35 Generates Test Statistic distribution for a model
44 self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
59 Extend ctools.csobservation getstate method to include some members
67 state = {
'base' : ctools.csobservation.__getstate__(self),
79 Extend ctools.csobservation setstate method to include some members
86 ctools.csobservation.__setstate__(self, state[
'base'])
88 self.
_fits = state[
'fits']
90 self.
_model = state[
'model']
99 Get parameters from parfile and setup the observation
102 if self.obs().size() == 0:
103 self.obs(self._get_observations())
106 self._set_obs_statistic(gammalib.toupper(self[
'statistic'].string()))
109 self.
_srcname = self[
'srcname'].string()
112 if self.obs().
models().size() == 0:
113 self.obs().
models(self[
'inmodel'].filename())
116 self[
'edisp'].boolean()
117 self[
'ntrials'].integer()
118 self[
'debug'].boolean()
121 if self._read_ahead():
122 self[
'outfile'].filename()
125 self._log_parameters(gammalib.TERSE)
135 Return a simulated observation container
140 Random number generator seed
144 sim : `~gammalib.GObservations`
145 Simulated observation container
149 if self.obs().size() == 1
and self.obs()[0].eventtype() ==
'CountsCube':
153 model = ctools.ctmodel(self.obs())
154 model[
'debug'] = self[
'debug'].boolean()
155 model[
'chatter'] = self[
'chatter'].integer()
157 self.
_model = model.cube().copy()
160 ran = gammalib.GRan()
163 counts = self._model.counts().copy()
166 for i
in range(counts.npix()):
167 counts[i] = ran.poisson(counts[i])
170 sim = self.obs().copy()
173 sim[0].events().counts(counts)
178 sim = obsutils.sim(self.obs(),
181 debug = self[
'debug'].boolean(),
189 Create the TS for a single trial
194 Random number generator seed
202 self._log_header2(gammalib.EXPLICIT,
'Trial %d' % (seed+1))
205 sim = self.
_sim(seed)
210 nevents += run.events().number()
213 self._log_header3(gammalib.EXPLICIT,
'Simulation')
214 self._log_value(gammalib.EXPLICIT,
'Number of simulated events', nevents)
217 fit = ctools.ctlike(sim)
219 fit[
'debug'] = self[
'debug'].boolean()
220 fit[
'chatter'] = self[
'chatter'].integer()
224 logL = fit.opt().value()
225 npred = fit.obs().npred()
226 models = fit.obs().
models()
231 if self._logExplicit():
232 self._log_header3(gammalib.EXPLICIT,
'Test source model fit')
233 self._log_value(gammalib.EXPLICIT,
'Test statistics', ts)
234 self._log_value(gammalib.EXPLICIT,
'log likelihood', logL)
235 self._log_value(gammalib.EXPLICIT,
'Number of predicted events', npred)
237 self._log_value(gammalib.EXPLICIT,
'Model', model.name())
239 self._log_string(gammalib.EXPLICIT, str(par))
240 elif self._logNormal():
241 prefactor = modutils.normalisation_parameter(model)
242 name =
'Trial %d' % seed
243 value =
'TS=%.3f %s=%e +/- %e' % \
244 (ts, prefactor.name(), prefactor.value(), prefactor.error())
245 self._log_value(gammalib.TERSE, name, value)
252 colnames.append(
'TS')
256 colnames.append(
'Nevents')
257 values[
'Nevents'] = nevents
260 colnames.append(
'Npred')
261 values[
'Npred'] = npred
265 model_name = model.name()
270 name = model_name +
'_' + par.name()
273 colnames.append(name)
274 values[name] = par.value()
278 colnames.append(name)
279 values[name] = par.error()
282 result = {
'colnames': colnames,
'values': values}
289 Create FITS file from results
293 results : list of dict
294 List of result dictionaries
298 for colname
in results[0][
'colnames']:
299 if colname !=
'TS' and colname !=
'Nevents' and \
301 headers.append(colname)
305 ts = gammalib.GFitsTableDoubleCol(
'TS', nrows)
306 nevents = gammalib.GFitsTableDoubleCol(
'NEVENTS', nrows)
307 npred = gammalib.GFitsTableDoubleCol(
'NPRED', nrows)
309 nevents.unit(
'counts')
312 for header
in headers:
313 name = gammalib.toupper(header)
314 column = gammalib.GFitsTableDoubleCol(name, nrows)
316 columns.append(column)
319 for i, result
in enumerate(results):
320 ts[i] = result[
'values'][
'TS']
321 nevents[i] = result[
'values'][
'Nevents']
322 npred[i] = result[
'values'][
'Npred']
323 for k, column
in enumerate(columns):
324 column[i] = result[
'values'][headers[k]]
327 table = gammalib.GFitsBinTable(nrows)
328 table.extname(
'TS_DISTRIBUTION')
331 table.card(
'INSTRUME',
'CTA',
'Name of Instrument')
332 table.card(
'TELESCOP',
'CTA',
'Name of Telescope')
338 table.card(
'NTRIALS', self[
'ntrials'].integer(),
'Number of trials')
339 table.card(
'STAT', self[
'statistic'].string(),
'Optimization statistic')
340 table.card(
'EDISP', self[
'edisp'].boolean(),
'Use energy dispersion?')
344 table.append(nevents)
346 for column
in columns:
350 self.
_fits = gammalib.GFits()
351 self._fits.append(table)
369 self._log_observations(gammalib.NORMAL, self.obs(),
'Input observation')
372 self._log_models(gammalib.NORMAL, self.obs().
models(),
'Input model')
375 self._log_header1(gammalib.TERSE,
'Generate TS distribution')
378 ntrials = self[
'ntrials'].integer()
385 args = [(self,
'_trial', i)
for i
in range(ntrials)]
386 poolresults = mputils.process(self.
_nthreads, mputils.mpfunc, args)
389 for i
in range(ntrials):
394 results.append(poolresults[i][0])
395 self._log_string(gammalib.TERSE, poolresults[i][1][
'log'],
False)
404 results.append(result)
418 models : `~gammalib.GModels`
429 Save TS distribution FITS file
432 self._log_header1(gammalib.TERSE,
'Save TS distribution')
435 if self.
_fits !=
None:
438 outfile = self[
'outfile'].filename()
441 self._log_value(gammalib.NORMAL,
'TS distribution file', outfile.url())
444 self._fits.saveto(outfile, self[
'clobber'].boolean())
451 Return TS distribution FITS file
454 FITS file containing TS distribution
463 if __name__ ==
'__main__':