Dmitry Malyshev / ctools

Commit 8085434844f98fec0e64fce3729f5e718ce65901

Authored by Jürgen Knödlseder
1 parent 003f01cc

Refactor cscripts and add some functions to obsutils (#1828) 

The obsutils.sim() function does not correctly handle stacked observations. It appends stacked response cubes to a counts cube in case that multiple observations have been simulated and stacked analysis was requested.

The stacked response cube computation is done using the obsutils.get_stacked_response() function.

The function obsutils.set_ts_model() has been added to setup a model for TS determination. This is used by some of the cscripts, and avoids code duplication.

Finally, the cstsdist and cspull scripts have been refactored.
Inline Side-by-side
Showing 6 changed files with 701 additions and 579 deletions
cscripts/cspull.py
  View file @8085434
... ... @@ -36,25 +36,19 @@ class cspull(ctools.cscript):
36 36 # Constructor
37 37 def __init__(self, *argv):
38 38 """
39   - Constructor.
  39 + Constructor
40 40 """
41 41 # Set name
42 42 self._name = 'cspull'
43   - self._version = '1.1.0'
  43 + self._version = '1.2.0'
44 44  
45 45 # Initialise some members
46 46 self._edisp = False
47   - self._exposure = None
48   - self._psfcube = None
49   - self._bckcube = None
50   - self._edispcube = None
51   - self._stackmodels = None
52 47 self._coordsys = 'CEL'
53 48 self._proj = 'TAN'
54 49 self._log_clients = False # Static parameter
55 50 self._offset = 0.0
56 51 self._ntrials = 0
57   - self._outfile = gammalib.GFilename()
58 52 self._npix = 0
59 53 self._binsz = 0.0
60 54 self._pattern = 'single'
... ... @@ -62,11 +56,10 @@ class cspull(ctools.cscript):
62 56 self._seed = 1
63 57 self._chatter = 2
64 58  
65   - # Initialise observation container from constructor arguments.
  59 + # Initialise observation container from constructor arguments
66 60 self._obs, argv = self._set_input_obs(argv)
67 61  
68   - # Initialise application by calling the appropriate class
69   - # constructor.
  62 + # Initialise application by calling the appropriate class constructor
70 63 self._init_cscript(argv)
71 64  
72 65 # Return
... ... @@ -108,7 +101,6 @@ class cspull(ctools.cscript):
108 101 self._obs.models(self['inmodel'].filename())
109 102  
110 103 # Read other parameters
111   - self._outfile = self['outfile'].filename()
112 104 self._ntrials = self['ntrials'].integer()
113 105 self._edisp = self['edisp'].boolean()
114 106 self._offset = self['offset'].real()
... ... @@ -116,6 +108,7 @@ class cspull(ctools.cscript):
116 108 self._chatter = self['chatter'].integer()
117 109  
118 110 # Query some parameters
  111 + self['outfile'].filename()
119 112 self['profile'].boolean()
120 113  
121 114 # Write input parameters into logger
... ... @@ -135,6 +128,9 @@ class cspull(ctools.cscript):
135 128 obs : `~gammalib.GObservations`
136 129 Observation container.
137 130 """
  131 + # Set duration
  132 + duration = self['tmax'].real() - self['tmin'].real()
  133 +
138 134 # Setup observation definition list
139 135 obsdeflist = obsutils.set_obs_patterns(self._pattern,
140 136 ra=self['ra'].real(),
... ... @@ -144,7 +140,7 @@ class cspull(ctools.cscript):
144 140 # Create list of observations
145 141 obs = obsutils.set_obs_list(obsdeflist,
146 142 tstart=self['tmin'].real(),
147   - duration=self['tmax'].real()-self['tmin'].real(),
  143 + duration=duration,
148 144 deadc=self['deadc'].real(),
149 145 emin=self['emin'].real(),
150 146 emax=self['emax'].real(),
... ... @@ -155,128 +151,6 @@ class cspull(ctools.cscript):
155 151 # Return observation container
156 152 return obs
157 153  
158   - def _set_stacked_irf(self):
159   - """
160   - Prepare stacked IRFs for stacked analysis
161   - """
162   - # Write header into logger
163   - if self._logTerse():
164   - self._log('\n')
165   - self._log.header1('Compute stacked response')
166   -
167   - # Set number of energy bins to at least 30 per energy decade
168   - enumbins = int((math.log10(self['emax'].real()) -
169   - math.log10(self['emin'].real())) * 30.0)
170   - if self._enumbins > enumbins:
171   - enumbins = self._enumbins
172   -
173   - # Compute spatial binning for point spread function and
174   - # energy dispersion cubes
175   - binsz = 10.0 * self['binsz'].real()
176   - nxpix = self['npix'].integer() // 10 # Make sure result is int
177   - nypix = self['npix'].integer() // 10 # Make sure result is int
178   - if nxpix < 2:
179   - nxpix = 2
180   - if nypix < 2:
181   - nypix = 2
182   -
183   - # Get stacked exposure
184   - expcube = ctools.ctexpcube(self._obs)
185   - expcube['incube'] = 'NONE'
186   - expcube['usepnt'] = True
187   - expcube['ebinalg'] = 'LOG'
188   - expcube['binsz'] = self._binsz
189   - expcube['nxpix'] = self._npix
190   - expcube['nypix'] = self._npix
191   - expcube['enumbins'] = enumbins
192   - expcube['emin'] = self['emin'].real()
193   - expcube['emax'] = self['emax'].real()
194   - expcube['coordsys'] = self._coordsys
195   - expcube['proj'] = self._proj
196   - expcube['chatter'] = self._chatter
197   - expcube['debug'] = self._logDebug()
198   - expcube.run()
199   -
200   - # Notify exposure computation
201   - if self._logTerse():
202   - self._log('Computed exposure cube\n')
203   -
204   - # Get stacked Psf
205   - psfcube = ctools.ctpsfcube(self._obs)
206   - psfcube['incube'] = 'NONE'
207   - psfcube['usepnt'] = True
208   - psfcube['ebinalg'] = 'LOG'
209   - psfcube['binsz'] = binsz
210   - psfcube['nxpix'] = nxpix
211   - psfcube['nypix'] = nypix
212   - psfcube['enumbins'] = enumbins
213   - psfcube['emin'] = self['emin'].real()
214   - psfcube['emax'] = self['emax'].real()
215   - psfcube['coordsys'] = self._coordsys
216   - psfcube['proj'] = self._proj
217   - psfcube['chatter'] = self._chatter
218   - psfcube['debug'] = self._logDebug()
219   - psfcube.run()
220   -
221   - # Notify Psf computation
222   - if self._logTerse():
223   - self._log('Computed point spread function cube\n')
224   -
225   - # Optionally get stacked Edisp
226   - if self._edisp:
227   - edispcube = ctools.ctedispcube(self._obs)
228   - edispcube['incube'] = 'NONE'
229   - edispcube['usepnt'] = True
230   - edispcube['ebinalg'] = 'LOG'
231   - edispcube['binsz'] = binsz
232   - edispcube['nxpix'] = nxpix
233   - edispcube['nypix'] = nypix
234   - edispcube['enumbins'] = enumbins
235   - edispcube['emin'] = self['emin'].real()
236   - edispcube['emax'] = self['emax'].real()
237   - edispcube['coordsys'] = self._coordsys
238   - edispcube['proj'] = self._proj
239   - edispcube['chatter'] = self._chatter
240   - edispcube['debug'] = self._logDebug()
241   - edispcube.run()
242   -
243   - # Store result
244   - self._edispcube = edispcube.edispcube().copy()
245   -
246   - # Logging
247   - if self._logTerse():
248   - self._log('Computed energy dispersion cube\n')
249   -
250   - # Get stacked background
251   - bkgcube = ctools.ctbkgcube(self._obs)
252   - bkgcube['incube'] = 'NONE'
253   - bkgcube['usepnt'] = True
254   - bkgcube['ebinalg'] = 'LOG'
255   - bkgcube['binsz'] = self._binsz
256   - bkgcube['nxpix'] = self._npix
257   - bkgcube['nypix'] = self._npix
258   - bkgcube['enumbins'] = enumbins
259   - bkgcube['emin'] = self['emin'].real()
260   - bkgcube['emax'] = self['emax'].real()
261   - bkgcube['coordsys'] = self._coordsys
262   - bkgcube['proj'] = self._proj
263   - bkgcube['chatter'] = self._chatter
264   - bkgcube['debug'] = self._logDebug()
265   - bkgcube.run()
266   -
267   - # Notify background cube computation
268   - if self._logTerse():
269   - self._log('Computed background cube\n')
270   -
271   - # Store results
272   - self._exposure = expcube.expcube().copy()
273   - self._psfcube = psfcube.psfcube().copy()
274   - self._bckcube = bkgcube.bkgcube().copy()
275   - self._stackmodels = bkgcube.models().copy()
276   -
277   - # Return
278   - return
279   -
280 154 def _obs_string(self, obs):
281 155 """
282 156 Generate summary string for observation
... ... @@ -310,7 +184,7 @@ class cspull(ctools.cscript):
310 184 # Return summary string
311 185 return text
312 186  
313   - def _trial(self, seed, stacked):
  187 + def _trial(self, seed):
314 188 """
315 189 Compute the pull for a single trial
316 190  
... ... @@ -318,8 +192,6 @@ class cspull(ctools.cscript):
318 192 ----------
319 193 seed : int
320 194 Random number generator seed
321   - stacked : bool
322   - Use stacked analysis
323 195  
324 196 Returns
325 197 -------
... ... @@ -354,16 +226,6 @@ class cspull(ctools.cscript):
354 226 debug=self._logDebug(),
355 227 chatter=self._chatter)
356 228  
357   - # Set response for a stacked observation
358   - if stacked:
359   - if self._edisp:
360   - obs[0].response(self._exposure, self._psfcube,
361   - self._edispcube, self._bckcube)
362   - else:
363   - obs[0].response(self._exposure, self._psfcube,
364   - self._bckcube)
365   - obs.models(self._stackmodels)
366   -
367 229 # Determine number of events in simulation
368 230 nevents = 0.0
369 231 for run in obs:
... ... @@ -427,10 +289,10 @@ class cspull(ctools.cscript):
427 289 # Set parameter name
428 290 name = model_name+'_'+par.name()
429 291  
430   - # Append parameter, Pull_parameter and Unc_parameter
  292 + # Append parameter, Pull_parameter and e_parameter
431 293 colnames.append(name)
432 294 colnames.append('Pull_'+name)
433   - colnames.append('Unc_'+name)
  295 + colnames.append('e_'+name)
434 296  
435 297 # Compute pull
436 298 fitted_value = par.value()
... ... @@ -444,7 +306,7 @@ class cspull(ctools.cscript):
444 306 # Store results
445 307 values[name] = fitted_value
446 308 values['Pull_'+name] = pull
447   - values['Unc_'+name] = error
  309 + values['e_'+name] = error
448 310  
449 311 # Write result
450 312 if self._logNormal():
... ... @@ -481,13 +343,10 @@ class cspull(ctools.cscript):
481 343 self._log.header1(gammalib.number('Observation',len(self._obs)))
482 344 self._log(str(self._obs))
483 345 self._log('\n')
484   -
485   - # Signal if stacked analysis is requested
486   - stacked = self._obs.size() > 1 and self._enumbins > 0
487   -
488   - # If several observations and binned: prepare stacked irfs
489   - if stacked:
490   - self._set_stacked_irf()
  346 + if self._logExplicit():
  347 + for obs in self._obs:
  348 + self._log(str(obs))
  349 + self._log('\n')
491 350  
492 351 # Write header
493 352 if self._logTerse():
... ... @@ -498,18 +357,21 @@ class cspull(ctools.cscript):
498 357 for seed in range(self._ntrials):
499 358  
500 359 # Make a trial and add initial seed
501   - result = self._trial(seed + self._seed, stacked)
  360 + result = self._trial(seed + self._seed)
  361 +
  362 + # Set output filename
  363 + outfile = self['outfile'].filename().url()
502 364  
503 365 # Write out result immediately
504 366 if seed == 0:
505   - f = open(self._outfile.url(), 'w')
  367 + f = open(outfile, 'w')
506 368 writer = csv.DictWriter(f, result['colnames'])
507 369 headers = {}
508 370 for n in result['colnames']:
509 371 headers[n] = n
510 372 writer.writerow(headers)
511 373 else:
512   - f = open(self._outfile.url(), 'a')
  374 + f = open(outfile, 'a')
513 375 writer = csv.DictWriter(f, result['colnames'])
514 376 writer.writerow(result['values'])
515 377 f.close()
... ...
cscripts/cstsdist.par
  View file @8085434
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: cstsdist.par #
8   -# Version: 1.1.0 #
  8 +# Version: 1.2.0 #
9 9 # #
10 10 # Author: Juergen Knoedlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -36,6 +36,8 @@ outfile, f, a, ts.dat,,, &quot;Output Test Statistic distribution file&quot;
36 36 ntrials, i, a, 10,,, "Number of trials"
37 37 ra, r, a, 83.63,0,360, "RA of pointing (deg)"
38 38 dec, r, a, 22.01,-90,90, "Dec of pointing (deg)"
  39 +coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
  40 +proj, s, a, TAN,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
39 41 emin, r, a, 0.1,,, "Lower energy limit (TeV)"
40 42 emax, r, a, 100.0,,, "Upper energy limit (TeV)"
41 43 enumbins, i, a, 0,,, "Number of energy bins (0=unbinned)"
... ...
cscripts/cstsdist.py
  View file @8085434
... ... @@ -30,13 +30,7 @@ from cscripts import obsutils
30 30 # ============== #
31 31 class cstsdist(ctools.cscript):
32 32 """
33   - Generates Test Statistic distribution for a model.
34   -
35   -
36   - This class implements the TS distribution generation script. It derives
37   - from the ctools.cscript class which provides support for parameter
38   - files, command line arguments, and logging. In that way the Python
39   - script behaves just as a regular ctool.
  33 + Generates Test Statistic distribution for a model
40 34 """
41 35  
42 36 # Constructor
... ... @@ -45,24 +39,17 @@ class cstsdist(ctools.cscript):
45 39 Constructor.
46 40 """
47 41 # Set name
48   - self._name = "cstsdist"
49   - self._version = "1.1.0"
  42 + self._name = 'cstsdist'
  43 + self._version = '1.2.0'
50 44  
51 45 # Initialise some members
52   - self._obs = gammalib.GObservations()
53   - self._pattern = "single"
54   - self._srcname = ""
55   - self._enumbins = 0
56   - self._npix = 0
57   - self._binsz = 0.0
58   - self._outfile = gammalib.GFilename("ts.dat")
59   - self._ntrials = 10
60   - self._edisp = False
61   - self._debug = False
62   - self._log_clients = False
63   -
64   - # Initialise application by calling the appropriate class
65   - # constructor.
  46 + self._srcname = ''
  47 + self._log_clients = False
  48 +
  49 + # Initialise observation container from constructor arguments
  50 + self._obs, argv = self._set_input_obs(argv)
  51 +
  52 + # Initialise application by calling the appropriate class constructor
66 53 self._init_cscript(argv)
67 54  
68 55 # Return
... ... @@ -74,165 +61,116 @@ class cstsdist(ctools.cscript):
74 61 """
75 62 Get parameters from parfile and setup the observation.
76 63 """
77   -
78 64 # Set observation if not done before
79 65 if self._obs == None or self._obs.size() == 0:
80 66 self._obs = self._get_observations()
81 67  
82   - # Check for requested pattern and use above
83   - # observation parameters to set wobble pattern
84   - self._pattern = self["pattern"].string()
85   - if self._pattern == "four":
  68 + # Check for requested pattern and use above observation parameters
  69 + # to set wobble pattern
  70 + if self['pattern'].string() == 'four':
86 71 self._obs = self._set_obs()
87 72  
88 73 # Get source name
89   - self._srcname = self["srcname"].string()
90   -
91   - # Get number of energy bins
92   - self._enumbins = self["enumbins"].integer()
  74 + self._srcname = self['srcname'].string()
93 75  
94   - # Read parameters for binned if requested
95   - if not self._enumbins == 0:
96   - self._npix = self["npix"].integer()
97   - self._binsz = self["binsz"].real()
  76 + # Query parameters for binned if requested
  77 + if self['enumbins'].integer() != 0:
  78 + self['npix'].integer()
  79 + self['binsz'].real()
  80 + self['coordsys'].string()
  81 + self['proj'].string()
98 82  
99 83 # Set models if we have none
100 84 if self._obs.models().size() == 0:
101   - self._obs.models(self["inmodel"].filename())
  85 + self._obs.models(self['inmodel'].filename())
102 86  
103   - # Get other parameters
104   - self._outfile = self["outfile"].filename()
105   - self._ntrials = self["ntrials"].integer()
106   - self._edisp = self["edisp"].boolean()
107   - self._debug = self["debug"].boolean()
  87 + # Query parameters
  88 + self['edisp'].boolean()
  89 + self['ntrials'].integer()
  90 + self['outfile'].filename()
  91 + self['debug'].boolean()
108 92  
109 93 # Write input parameters into logger
110 94 if self._logTerse():
111 95 self._log_parameters()
112   - self._log("\n")
  96 + self._log('\n')
113 97  
114 98 # Return
115 99 return
116 100  
117 101 def _set_obs(self):
118 102 """
119   - Set an observation container.
  103 + Set an observation container
120 104  
121   - Returns:
122   - Observation container.
  105 + Returns
  106 + -------
  107 + obs : `~gammalib.GObservations`
  108 + Observation container
123 109 """
  110 + # Set duration
  111 + duration = self['tmax'].real() - self['tmin'].real()
124 112  
125 113 # Setup observation definition list
126   - obsdeflist = obsutils.set_obs_patterns(self._pattern,
127   - ra=self["ra"].real(),
128   - dec=self["dec"].real(),
129   - offset=self["offset"].real())
  114 + obsdeflist = obsutils.set_obs_patterns(self['pattern'].string(),
  115 + ra = self['ra'].real(),
  116 + dec = self['dec'].real(),
  117 + offset = self['offset'].real())
130 118  
131 119 # Create list of observations
132 120 obs = obsutils.set_obs_list(obsdeflist,
133   - tstart=self["tmin"].real(),
134   - duration=self["tmax"].real()-self["tmin"].real(),
135   - deadc=self["deadc"].real(),
136   - emin=self["emin"].real(),
137   - emax=self["emax"].real(),
138   - rad=self["rad"].real(),
139   - irf=self["irf"].string(),
140   - caldb=self["caldb"].string())
  121 + tstart = self['tmin'].real(),
  122 + duration = duration,
  123 + deadc = self['deadc'].real(),
  124 + emin = self['emin'].real(),
  125 + emax = self['emax'].real(),
  126 + rad = self['rad'].real(),
  127 + irf = self['irf'].string(),
  128 + caldb = self['caldb'].string())
141 129  
142 130 # Return observation container
143 131 return obs
144 132  
145   - def _set_models(self, fitspat=False, fitspec=False):
  133 + def _trial(self, seed):
146 134 """
147   - Set full model and background model.
  135 + Create the TS for a single trial
148 136  
149   - Kwargs:
150   - fitspat: Fit spatial parameter (default: False).
151   - fitspec: Fit spectral parameters (default: False).
152   -
153   - Returns:
154   - Tuple containing full model and background model.
155   - """
156   - # Retrieve full model from observation container
157   - full_model = self._obs.models().copy()
158   -
159   - # Get source model
160   - model = full_model[self._srcname]
161   -
162   - # Check that model has a Prefactor
163   - if not model.has_par("Prefactor"):
164   - msg = "Model \""+self._srcname+"\" has no parameter "+\
165   - "\"Prefactor\". Only spectral models with a "+\
166   - "\"Prefactor\" parameter are supported."
167   - raise RuntimeError(msg)
168   -
169   - # Fit or fix spatial parameters
170   - if fitspat:
171   - if model.has_par("RA"):
172   - model["RA"].free()
173   - if model.has_par("DEC"):
174   - model["DEC"].free()
175   - if model.has_par("Sigma"):
176   - model["Sigma"].free()
177   - if model.has_par("Radius"):
178   - model["Radius"].free()
179   - if model.has_par("Width"):
180   - model["Width"].free()
181   - else:
182   - if model.has_par("RA"):
183   - model["RA"].fix()
184   - if model.has_par("DEC"):
185   - model["DEC"].fix()
186   - if model.has_par("Sigma"):
187   - model["Sigma"].fix()
188   - if model.has_par("Radius"):
189   - model["Radius"].fix()
190   - if model.has_par("Width"):
191   - model["Width"].fix()
192   -
193   - # Fit or fix spectral parameters
194   - if fitspec:
195   - if model.has_par("Index"):
196   - model["Index"].free()
197   - if model.has_par("Cutoff"):
198   - model["Cutoff"].free()
199   - else:
200   - if model.has_par("Index"):
201   - model["Index"].fix()
202   - if model.has_par("Cutoff"):
203   - model["Cutoff"].fix()
204   -
205   - # Create background model
206   - bkg_model = full_model.copy()
207   - bkg_model.remove(self._srcname)
208   -
209   - # Return models
210   - return full_model, bkg_model
211   -
212   - def _trial(self, seed, full_model, bkg_model):
213   - """
214   - Create the TS for a single trial.
  137 + Parameters
  138 + ----------
  139 + seed : int
  140 + Random number generator seed
215 141  
216   - Args:
217   - seed: Random number generator seed
218   - full_model: Full model
219   - bkg_model: Background model
220   -
221   - Returns:
  142 + Returns
  143 + -------
  144 + result : dict
222 145 Result dictionary
223 146 """
224 147 # Write header
225 148 if self._logExplicit():
226   - self._log.header2("Trial "+str(seed+1))
  149 + self._log.header2('Trial '+str(seed+1))
  150 +
  151 + # Set default binned parameters
  152 + coordsys = 'CEL'
  153 + proj = 'TAN'
  154 + npix = 0
  155 + binsz = 0.0
  156 +
  157 + # If binned analysis is requested the read the binned parameters
  158 + if self['enumbins'].integer() != 0:
  159 + npix = self['npix'].integer()
  160 + binsz = self['binsz'].real()
  161 + coordsys = self['coordsys'].string()
  162 + proj = self['proj'].string()
227 163  
228 164 # Simulate events
229 165 sim = obsutils.sim(self._obs,
230   - nbins=self._enumbins,
231   - seed=seed,
232   - binsz=self._binsz,
233   - npix=self._npix,
234   - log=self._log_clients,
235   - debug=self._debug)
  166 + nbins = self['enumbins'].integer(),
  167 + seed = seed,
  168 + proj = proj,
  169 + coord = coordsys,
  170 + binsz = binsz,
  171 + npix = npix,
  172 + log = self._log_clients,
  173 + debug = self['debug'].boolean())
236 174  
237 175 # Determine number of events in simulation
238 176 nevents = 0.0
... ... @@ -241,110 +179,70 @@ class cstsdist(ctools.cscript):
241 179  
242 180 # Write simulation results
243 181 if self._logExplicit():
244   - self._log.header3("Simulation")
245   - self._log.parformat("Number of simulated events")
  182 + self._log.header3('Simulation')
  183 + self._log.parformat('Number of simulated events')
246 184 self._log(nevents)
247   - self._log("\n")
248   -
249   - # Fit background only
250   - sim.models(bkg_model)
251   - like_bgm = obsutils.fit(sim,
252   - log=self._log_clients,
253   - debug=self._debug)
254   - result_bgm = like_bgm.obs().models().copy()
255   - LogL_bgm = like_bgm.opt().value()
256   - npred_bgm = like_bgm.obs().npred()
257   -
258   - # Write background fit results
259   - if self._logExplicit():
260   - self._log.header3("Background model fit")
261   - self._log.parformat("log likelihood")
262   - self._log(LogL_bgm)
263   - self._log("\n")
264   - self._log.parformat("Number of predicted events")
265   - self._log(npred_bgm)
266   - self._log("\n")
267   - for model in result_bgm:
268   - self._log.parformat("Model")
269   - self._log(model.name())
270   - self._log("\n")
271   - for par in model:
272   - self._log(str(par)+"\n")
273   -
274   - # Fit background and test source
275   - sim.models(full_model)
276   - like_all = obsutils.fit(sim,
277   - log=self._log_clients,
278   - debug=self._debug)
279   - result_all = like_all.obs().models().copy()
280   - LogL_all = like_all.opt().value()
281   - npred_all = like_all.obs().npred()
282   - ts = 2.0*(LogL_bgm-LogL_all)
283   -
284   - # Write background and test source fit results
  185 + self._log('\n')
  186 +
  187 + # Fit model
  188 + fit = obsutils.fit(sim, log = self._log_clients,
  189 + debug = self['debug'].boolean())
  190 +
  191 + # Get model fitting results
  192 + logL = fit.opt().value()
  193 + npred = fit.obs().npred()
  194 + models = fit.obs().models()
  195 + model = models[self._srcname]
  196 + ts = model.ts()
  197 +
  198 + # Write fit results, either explicit or normal
285 199 if self._logExplicit():
286   - self._log.header3("Background and test source model fit")
287   - self._log.parformat("Test statistics")
  200 + self._log.header3('Test source model fit')
  201 + self._log.parformat('Test statistics')
288 202 self._log(ts)
289   - self._log("\n")
290   - self._log.parformat("log likelihood")
291   - self._log(LogL_all)
292   - self._log("\n")
293   - self._log.parformat("Number of predicted events")
294   - self._log(npred_all)
295   - self._log("\n")
296   - for model in result_all:
297   - self._log.parformat("Model")
  203 + self._log('\n')
  204 + self._log.parformat('log likelihood')
  205 + self._log(logL)
  206 + self._log('\n')
  207 + self._log.parformat('Number of predicted events')
  208 + self._log(npred)
  209 + self._log('\n')
  210 + for model in models:
  211 + self._log.parformat('Model')
298 212 self._log(model.name())
299   - self._log("\n")
  213 + self._log('\n')
300 214 for par in model:
301   - self._log(str(par)+"\n")
302   -
303   - # Write result
304   - elif self._logTerse():
305   - self._log.parformat("Trial "+str(seed))
306   - self._log("TS=")
  215 + self._log(str(par)+'\n')
  216 + elif self._logNormal():
  217 + self._log.parformat('Trial '+str(seed))
  218 + self._log('TS=')
307 219 self._log(ts)
308   - self._log(" Prefactor=")
309   - self._log(result_all[self._srcname]["Prefactor"].value())
310   - self._log("+/-")
311   - self._log(result_all[self._srcname]["Prefactor"].error())
312   - self._log("\n")
  220 + self._log(' Prefactor=')
  221 + self._log(model['Prefactor'].value())
  222 + self._log('+/-')
  223 + self._log(model['Prefactor'].error())
  224 + self._log('\n')
313 225  
314 226 # Initialise results
315 227 colnames = []
316 228 values = {}
317 229  
318 230 # Set TS value
319   - colnames.append("TS")
320   - values["TS"] = ts
321   -
322   - # Set logL for background fit
323   - colnames.append("LogL_bgm")
324   - values["LogL_bgm"] = LogL_bgm
325   -
326   - # Set logL for full fit
327   - colnames.append("LogL_all")
328   - values["LogL_all"] = LogL_all
  231 + colnames.append('TS')
  232 + values['TS'] = ts
329 233  
330 234 # Set Nevents
331   - colnames.append("Nevents")
332   - values["Nevents"] = nevents
333   -
334   - # Set Npred for background fit
335   - colnames.append("Npred_bkg")
336   - values["Npred_bkg"] = npred_bgm
  235 + colnames.append('Nevents')
  236 + values['Nevents'] = nevents
337 237  
338   - # Set Npred for full fit
339   - colnames.append("Npred_all")
340   - values["Npred_all"] = npred_all
  238 + # Set Npred
  239 + colnames.append('Npred')
  240 + values['Npred'] = npred
341 241  
342 242 # Gather free full fit parameters
343   - for i in range(result_all.size()):
344   - model = result_all[i]
  243 + for model in models:
345 244 model_name = model.name()
346   - for k in range(model.size()):
347   - par = model[k]
  245 + for par in model:
348 246 if par.is_free():
349 247  
350 248 # Set parameter name
... ... @@ -355,7 +253,7 @@ class cstsdist(ctools.cscript):
355 253 values[name] = par.value()
356 254  
357 255 # Append error
358   - name = "Unc_"+name
  256 + name = 'e_'+name
359 257 colnames.append(name)
360 258 values[name] = par.error()
361 259  
... ... @@ -369,7 +267,7 @@ class cstsdist(ctools.cscript):
369 267 # Public methods
370 268 def run(self):
371 269 """
372   - Run the script.
  270 + Run the script
373 271 """
374 272 # Switch screen logging on in debug mode
375 273 if self._logDebug():
... ... @@ -378,55 +276,51 @@ class cstsdist(ctools.cscript):
378 276 # Get parameters
379 277 self._get_parameters()
380 278  
381   - # Initialise models
382   - full_model, bkg_model = self._set_models()
  279 + # Set models for this observation
  280 + self.models(obsutils.set_ts_model(self._obs.models(), self._srcname))
383 281  
384 282 # Write models into logger
385 283 if self._logTerse():
386   - self._log("\n")
387   - self._log.header1("Models")
388   - self._log.header2("Background model")
389   - self._log(str(bkg_model))
390   - self._log("\n\n")
391   - self._log.header2("Full model")
392   - self._log(str(full_model))
393   - self._log("\n")
  284 + self._log('\n')
  285 + self._log.header1('Models')
  286 + self._log(str(self._obs.models()))
  287 + self._log('\n')
394 288  
395 289 # Write observation into logger
396 290 if self._logTerse():
397   - self._log("\n")
398   - self._log.header1("Observation")
  291 + self._log('\n')
  292 + self._log.header1(gammalib.number('Observation',len(self._obs)))
399 293 self._log(str(self._obs))
400   - self._log("\n")
401   -
402   - # Write models into logger
403   - if self._logTerse():
404   - self._log("\n")
405   - self._log.header1("Test model")
406   - self._log(str(full_model))
407   - self._log("\n")
  294 + self._log('\n')
  295 + if self._logExplicit():
  296 + for obs in self._obs:
  297 + self._log(str(obs))
  298 + self._log('\n')
408 299  
409 300 # Write header
410 301 if self._logTerse():
411   - self._log("\n")
412   - self._log.header1("Generate TS distribution")
  302 + self._log('\n')
  303 + self._log.header1('Generate TS distribution')
413 304  
414 305 # Loop over trials
415   - for seed in range(self._ntrials):
  306 + for seed in range(self['ntrials'].integer()):
416 307  
417 308 # Make a trial
418   - result = self._trial(seed, full_model, bkg_model)
  309 + result = self._trial(seed)
  310 +
  311 + # Set output filename
  312 + outfile = self['outfile'].filename().url()
419 313  
420 314 # Write out result immediately
421 315 if seed == 0:
422   - f = open(self._outfile.url(), 'w')
  316 + f = open(outfile, 'w')
423 317 writer = csv.DictWriter(f, result['colnames'])
424 318 headers = {}
425 319 for n in result['colnames']:
426 320 headers[n] = n
427 321 writer.writerow(headers)
428 322 else:
429   - f = open(self._outfile.url(), 'a')
  323 + f = open(outfile, 'a')
430 324 writer = csv.DictWriter(f, result['colnames'])
431 325 writer.writerow(result['values'])
432 326 f.close()
... ... @@ -436,17 +330,22 @@ class cstsdist(ctools.cscript):
436 330  
437 331 def models(self, models):
438 332 """
439   - Set model.
  333 + Set model
  334 +
  335 + Parameters
  336 + ----------
  337 + models : `~gammalib.GModels`
  338 + Model container
440 339 """
441 340 # Copy models
442   - self._obs.models(models.copy())
  341 + self._obs.models(models)
443 342  
444 343 # Return
445 344 return
446 345  
447 346 def execute(self):
448 347 """
449   - Execute the script.
  348 + Execute the script
450 349 """
451 350 # Open logfile
452 351 self.logFileOpen()
... ...
cscripts/obsutils.py
  View file @8085434
1 1 # ==========================================================================
2   -# CTA observation handling support functions.
  2 +# CTA observation handling support functions
3 3 #
4 4 # Copyright (C) 2011-2016 Juergen Knoedlseder
5 5 #
... ... @@ -30,6 +30,14 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
30 30 """
31 31 Simulate events for all observations in the container
32 32  
  33 + Simulate events for all observations using ctobssim. If the number of
  34 + energy bins is positive, the events are then binned in a counts cube
  35 + using ctbin. If multiple observations are simulated, the counts cube is
  36 + a stacked cube and the corresponding response cubes are computed using
  37 + ctexpcube, ctpsfcube, ctbkgcube and optionally ctedispcube. The response
  38 + cubes are attached to the first observation in the container, which
  39 + normally is the observation with the counts cube.
  40 +
33 41 Parameters
34 42 ----------
35 43 obs : `~gammalib.GObservations`
... ... @@ -42,7 +50,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
42 50 Chatter level
43 51 edisp : bool, optional
44 52 Apply energy dispersion?
45   - seed : int, integer
  53 + seed : int, optional
46 54 Seed value for simulations
47 55 emin : float, optional
48 56 Minimum energy of counts cube for binned (TeV)
... ... @@ -64,7 +72,6 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
64 72 obs : `~gammalib.GObservations`
65 73 Observation container filled with simulated events
66 74 """
67   -
68 75 # Allocate ctobssim application and set parameters
69 76 sim = ctools.ctobssim(obs)
70 77 sim['seed'] = seed
... ... @@ -125,6 +132,48 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
125 132 # the container and bin the events in counts maps
126 133 bin.run()
127 134  
  135 + # If we have multiple input observations then create stacked response
  136 + # cubes and append them to the observation
  137 + if len(sim.obs()) > 1:
  138 +
  139 + # First get the cube. We need this because ctbin builds a cube on
  140 + # the fly. We should avoid this and make the cube part of ctbin
  141 + cube = bin.cube()
  142 +
  143 + # Get xref, yref, coordsys and proj from counts cube
  144 + _xref = cube.counts().projection().crval(0)
  145 + _yref = cube.counts().projection().crval(1)
  146 + _coordsys = cube.counts().projection().coordsys()
  147 + _proj = cube.counts().projection().code()
  148 +
  149 + # Kluge: map "EQU" to "CEL" coordinate system (not sure why the
  150 + # GSkyProjection::coordsys returns "EQU" instead of "CEL"
  151 + if _coordsys == 'EQU':
  152 + _coordsys = 'CEL'
  153 +
  154 + # Get stacked response
  155 + response = get_stacked_response(sim.obs(), _xref, _yref,
  156 + binsz=binsz, nxpix=npix, nypix=npix,
  157 + emin=emin, emax=emax, enumbins=nbins,
  158 + edisp=edisp,
  159 + coordsys=_coordsys, proj=_proj,
  160 + log=log, debug=debug,
  161 + chatter=chatter)
  162 +
  163 + # Set stacked response
  164 + if edisp:
  165 + bin.obs()[0].response(response['expcube'],
  166 + response['psfcube'],
  167 + response['edispcube'],
  168 + response['bkgcube'])
  169 + else:
  170 + bin.obs()[0].response(response['expcube'],
  171 + response['psfcube'],
  172 + response['bkgcube'])
  173 +
  174 + # Set new models
  175 + bin.obs().models(response['models'])
  176 +
128 177 # Make a deep copy of the observation that will be returned
129 178 # (the ctbin object will go out of scope one the function is
130 179 # left)
... ... @@ -149,33 +198,36 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
149 198 # ================ #
150 199 def fit(obs, log=False, debug=False, chatter=2, edisp=False):
151 200 """
152   - Perform maximum likelihood fitting of observations in the container.
153   -
154   - Parameters:
155   - obs - Observation container
156   - Keywords:
157   - log - Create log file(s)
158   - debug - Create screen dump
159   - chatter - Chatter level
160   - edisp - Apply energy dispersion?
  201 + Perform maximum likelihood fitting of observations in the container
  202 +
  203 + Parameters
  204 + ----------
  205 + obs : `~gammalib.GObservations`
  206 + Observation container
  207 + log : bool, optional
  208 + Create log file(s)
  209 + debug : bool, optional
  210 + Create console dump?
  211 + chatter : int, optional
  212 + Chatter level
  213 + edisp : bool, optional
  214 + Apply energy dispersion?
  215 +
  216 + Returns
  217 + -------
  218 + like : `~ctools.ctlike`
  219 + ctlike application
161 220 """
162   - # Allocate ctlike application
  221 + # Allocate ctlike application and set parameters
163 222 like = ctools.ctlike(obs)
  223 + like['debug'] = debug
  224 + like['chatter'] = chatter
  225 + like['edisp'] = edisp
164 226  
165 227 # Optionally open the log file
166 228 if log:
167 229 like.logFileOpen()
168 230  
169   - # Optionally switch-on debugging model
170   - if debug:
171   - like["debug"] = True
172   -
173   - # Set chatter level
174   - like["chatter"] = chatter
175   -
176   - # Optionally apply energy dispersion
177   - like["edisp"] = edisp
178   -
179 231 # Run ctlike application.
180 232 like.run()
181 233  
... ... @@ -188,34 +240,36 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False):
188 240 # ============================================================== #
189 241 def cterror(obs, srcname, log=False, debug=False, chatter=2):
190 242 """
191   - Perform maximum likelihood fitting of observations in the container.
192   -
193   - Parameters:
194   - obs - Observation container
195   - srcname - Source name
196   - Keywords:
197   - log - Create log file(s)
198   - debug - Create screen dump
199   - chatter - Chatter level
200   - edisp - Apply energy dispersion?
  243 + Fit observations and determine errors using likelihood profile
  244 +
  245 + Parameters
  246 + ----------
  247 + obs : `~gammalib.GObservations`
  248 + Observation container
  249 + srcname : str
  250 + Source name
  251 + log : bool, optional
  252 + Create log file(s)
  253 + debug : bool, optional
  254 + Create console dump?
  255 + chatter : int, optional
  256 + Chatter level
  257 +
  258 + Returns
  259 + -------
  260 + error : `~ctools.cterror`
  261 + cterror application
201 262 """
202   - # Allocate cterror application
  263 + # Allocate cterror application and set parameters
203 264 error = ctools.cterror(obs)
204   -
205   - # Set cterror parameters
206   - error["srcname"] = srcname
  265 + error['srcname'] = srcname
  266 + error['debug'] = debug
  267 + error['chatter'] = chatter
207 268  
208 269 # Optionally open the log file
209 270 if log:
210 271 error.logFileOpen()
211 272  
212   - # Optionally switch-on debugging model
213   - if debug:
214   - error["debug"] = True
215   -
216   - # Set chatter level
217   - error["chatter"] = chatter
218   -
219 273 # Run cterror application.
220 274 error.run()
221 275  
... ... @@ -226,24 +280,39 @@ def cterror(obs, srcname, log=False, debug=False, chatter=2):
226 280 # ================= #
227 281 # Create counts map #
228 282 # ================= #
229   -def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
  283 +def cntmap(obs, proj='TAN', coord='GAL', xval=0.0, yval=0.0, \
230 284 binsz=0.05, nxpix=200, nypix=200, \
231   - outname="cntmap.fits"):
  285 + outname='cntmap.fits'):
232 286 """
233   - Creates a counts map by combining the events of all observations.
234   - The counts map will be a summed map over all energies.
235   -
236   - Parameters:
237   - obs - Observation container
238   - Keywords:
239   - proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN)
240   - coord - Coordinate type (GAL, CEL) (default: GAL)
241   - xval - Reference longitude value [deg] (default: 0.0)
242   - yval - Reference latitude value [deg] (default: 0.0)
243   - binsz - Pixel size [deg/pixel] (default: 0.05)
244   - nxpix - Number of pixels in X direction (default: 200)
245   - nypix - Number of pixels in Y direction (default: 200)
246   - outname - Counts map FITS filename (default: cntmap.fits)
  287 + Creates a counts map by combining the events of all observations
  288 +
  289 + The counts map will be a summed map over all energies
  290 +
  291 + Parameters
  292 + ----------
  293 + obs : `~gammalib.GObservations`
  294 + Observation container without events
  295 + proj : str, optional
  296 + Projection for binned simulation
  297 + coord : str, optional
  298 + Coordinate system for binned simulation
  299 + xval : float, optional
  300 + Counts map centre of X direction (deg)
  301 + yval : float, optional
  302 + Counts map centre of Y direction (deg)
  303 + binsz : float, optional
  304 + Pixel size for binned simulation (deg/pixel)
  305 + nxpix : int, optional
  306 + Number of pixels in X direction
  307 + nypix : int, optional
  308 + Number of pixels in Y direction
  309 + outname : str, optional
  310 + Counts map file name
  311 +
  312 + Returns
  313 + -------
  314 + map : `~gammalib.GSkyMap`
  315 + Counts map
247 316 """
248 317 # Allocate counts map
249 318 map = gammalib.GSkyMap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
... ... @@ -276,26 +345,42 @@ def cntmap(obs, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
276 345 # ================ #
277 346 # Create model map #
278 347 # ================ #
279   -def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
  348 +def modmap(obs, eref=0.1, proj='TAN', coord='GAL', xval=0.0, yval=0.0, \
280 349 binsz=0.05, nxpix=200, nypix=200, \
281   - outname="modmap.fits"):
  350 + outname='modmap.fits'):
282 351 """
283   - Make model map for a given reference energy by combining all observations.
  352 + Make model map for a given reference energy by combining all observations
  353 +
284 354 The model map will be evaluated for a given reference energy 'eref' and will
285 355 be given in units of [counts/(sr MeV s)].
286 356  
287   - Parameters:
288   - obs - Observation container
289   - Keywords:
290   - eref - Reference energy for which model is created [TeV] (default: 0.1)
291   - proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN)
292   - coord - Coordinate type (GAL, CEL) (default: GAL)
293   - xval - Reference longitude value [deg] (default: 0.0)
294   - yval - Reference latitude value [deg] (default: 0.0)
295   - binsz - Pixel size [deg/pixel] (default: 0.05)
296   - nxpix - Number of pixels in X direction (default: 200)
297   - nypix - Number of pixels in Y direction (default: 200)
298   - outname - Model map FITS filename (default: modmap.fits)
  357 + Parameters
  358 + ----------
  359 + obs : `~gammalib.GObservations`
  360 + Observation container without events
  361 + eref : float, optional
  362 + Reference energy for which model is created (TeV)
  363 + proj : str, optional
  364 + Projection for binned simulation
  365 + coord : str, optional
  366 + Coordinate system for binned simulation
  367 + xval : float, optional
  368 + Counts map centre of X direction (deg)
  369 + yval : float, optional
  370 + Counts map centre of Y direction (deg)
  371 + binsz : float, optional
  372 + Pixel size for binned simulation (deg/pixel)
  373 + nxpix : int, optional
  374 + Number of pixels in X direction
  375 + nypix : int, optional
  376 + Number of pixels in Y direction
  377 + outname : str, optional
  378 + Model map file name
  379 +
  380 + Returns
  381 + -------
  382 + map : `~gammalib.GSkyMap`
  383 + Model map
299 384 """
300 385 # Allocate model map
301 386 map = gammalib.GSkyMap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
... ... @@ -338,53 +423,46 @@ def modmap(obs, eref=0.1, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
338 423 return map
339 424  
340 425  
341   -# ===================================== #
342   -# Set one CTA observation (old version) #
343   -# ===================================== #
344   -def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
345   - emin=0.1, emax=100.0, rad=5.0, \
346   - irf="South_50h", caldb="prod2"):
347   - """
348   - Obsolete function, use set_obs instead.
349   - """
350   - # Print warning
351   - print("Warning: obsutils.set is obsolete, use obsutils.set_obs instead.")
352   -
353   - # Call new function
354   - obs_cta = set_obs(pntdir, tstart=tstart, duration=duration, deadc=deadc, \
355   - emin=emin, emax=emax, rad=rad, \
356   - irf=irf, caldb=caldb)
357   -
358   - # Return CTA observation
359   - return obs_cta
360   -
361   -
362 426 # ======================= #
363 427 # Set one CTA observation #
364 428 # ======================= #
365 429 def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
366 430 emin=0.1, emax=100.0, rad=5.0, \
367   - irf="South_50h", caldb="prod2", id="000000"):
  431 + irf='South_50h', caldb='prod2', id='000000'):
368 432 """
369   - Set a single CTA observation.
  433 + Set a single CTA observation
370 434  
371 435 The function sets a single CTA observation containing an empty CTA
372 436 event list. By looping over this function you can add CTA observations
373 437 to the observation container.
374 438  
375   - Args:
376   - pntdir: Pointing direction [GSkyDir]
377   -
378   - Kwargs:
379   - tstart: Start time (seconds) (default: 0.0)
380   - duration: Duration of observation (seconds) (default: 1800.0)
381   - deadc: Deadtime correction factor (default: 0.95)
382   - emin: Minimum event energy (TeV) (default: 0.1)
383   - emax: Maximum event energy (TeV) (default: 100.0)
384   - rad: ROI radius used for analysis (deg) (default: 5.0)
385   - irf: Instrument response function (default: "South_50h")
386   - caldb: Calibration database path (default: "prod2")
387   - id: Run identifier (default: "000000")
  439 + Parameters
  440 + ----------
  441 + pntdir : `~gammalib.GSkyDir`
  442 + Pointing direction
  443 + tstart : float, optional
  444 + Start time (seconds)
  445 + duration : float, optional
  446 + Duration of observation (seconds)
  447 + deadc : float, optional
  448 + Deadtime correction factor
  449 + emin : float, optional
  450 + Minimum event energy (TeV)
  451 + emax : float, optional
  452 + Maximum event energy (TeV)
  453 + rad : float, optional
  454 + ROI radius used for analysis (deg)
  455 + irf : str, optional
  456 + Instrument response function
  457 + caldb : str, optional
  458 + Calibration database path
  459 + id : str, optional
  460 + Observation identifier
  461 +
  462 + Returns
  463 + -------
  464 + obs : `~gammalib.GCTAObservation`
  465 + CTA observation
388 466 """
389 467 # Allocate CTA observation
390 468 obs_cta = gammalib.GCTAObservation()
... ... @@ -394,7 +472,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
394 472 if (gammalib.dir_exists(caldb)):
395 473 db.rootdir(caldb)
396 474 else:
397   - db.open("cta", caldb)
  475 + db.open('cta', caldb)
398 476  
399 477 # Set pointing direction
400 478 pnt = gammalib.GCTAPointing()
... ... @@ -413,8 +491,8 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
413 491 gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration))
414 492  
415 493 # Set energy boundaries
416   - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"),
417   - gammalib.GEnergy(emax, "TeV"))
  494 + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, 'TeV'),
  495 + gammalib.GEnergy(emax, 'TeV'))
418 496  
419 497 # Allocate event list
420 498 events = gammalib.GCTAEventList()
... ... @@ -440,10 +518,11 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
440 518 # Set list of CTA observations #
441 519 # ============================ #
442 520 def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
443   - emin=0.1, emax=100.0, rad=5.0, \
444   - irf="South_50h", caldb="prod2"):
  521 + emin=0.1, emax=100.0, rad=5.0, \
  522 + irf='South_50h', caldb='prod2'):
445 523 """
446   - Returns an observation container filled with a list of CTA observations.
  524 + Returns an observation container filled with a list of CTA observations
  525 +
447 526 The list is defined by the obsdeflist parameter which is a dictionnary
448 527 containing the mandatory keywords 'ra' and 'dec' that specify the
449 528 pointing direction for a given observation. Optional keyword give control
... ... @@ -451,27 +530,31 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
451 530 energy range, etc. If an optional keyword is not specified, the function
452 531 keyword is used instead.
453 532  
454   - Parameters:
455   - obsdeflist - Observation definition list [{'ra': x.xx, 'dec': x.xx}]
456   - The directory can take the following optional keywords:
457   - - duration: Duration of CTA observation [seconds]
458   - - deadc: Deadtime correction factor
459   - - emin: Minimum event energy [TeV]
460   - - emax: Maximum event energy [TeV]
461   - - rad: ROI radius used for analysis [deg]
462   - - irf: Instrument response function
463   - - caldb: Calibration database path
464   - Optional keywords overwrite keywords specified in the
465   - function call.
466   - Keywords:
467   - tstart - Start time [seconds] (default: 0.0)
468   - duration - Duration of one CTA observation [seconds] (default: 1800.0)
469   - deadc - Deadtime correction factor (default: 0.95)
470   - emin - Minimum event energy [TeV] (default: 0.1)
471   - emax - Maximum event energy [TeV] (default: 100.0)
472   - rad - ROI radius used for analysis [deg] (default: 5.0)
473   - irf - Instrument response function (default: cta_dummy_irf)
474   - caldb - Calibration database path (default: "dummy")
  533 + Parameters
  534 + ----------
  535 + obsdeflist : list of dict
  536 + Observation definition list
  537 + tstart : float, optional
  538 + Start time (seconds)
  539 + duration : float, optional
  540 + Duration of observation (seconds)
  541 + deadc : float, optional
  542 + Deadtime correction factor
  543 + emin : float, optional
  544 + Minimum event energy (TeV)
  545 + emax : float, optional
  546 + Maximum event energy (TeV)
  547 + rad : float, optional
  548 + ROI radius used for analysis (deg)
  549 + irf : str, optional
  550 + Instrument response function
  551 + caldb : str, optional
  552 + Calibration database path
  553 +
  554 + Returns
  555 + -------
  556 + obs : `~gammalib.GObservations`
  557 + Observation container filled with CTA observation
475 558 """
476 559 # Initialise empty observation container
477 560 obs = gammalib.GObservations()
... ... @@ -530,7 +613,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
530 613 obs_irf = irf
531 614  
532 615 # Generate identifier string
533   - id = "%6.6d" % obs_id
  616 + id = '%6.6d' % obs_id
534 617  
535 618 # Set CTA observation
536 619 obs_cta = set_obs(pntdir, tstart=obs_start, duration=obs_duration, \
... ... @@ -553,26 +636,38 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
553 636 # ============================ #
554 637 def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
555 638 """
556   - Sets a number of standard patterns.
557   -
558   - Parameters:
559   - pattern - Observation pattern. Possible options are:
560   - - "single": single pointing
561   - - "four": four pointings 'offset' around pattern centre
  639 + Sets a number of standard patterns
  640 +
  641 + Parameters
  642 + ----------
  643 + pattern : str
  644 + Observation pattern ("single", "four")
  645 + ra : float, optional
  646 + Right Ascension of pattern centre (deg)
  647 + dec : float, optional
  648 + Declination of pattern centre (deg)
  649 + offset : float, optional
  650 + Offset from pattern centre (deg)
562 651  
563   - Keywords:
564   - ra - Right Ascension of pattern centre [deg] (default: 83.6331)
565   - dec - Declination of pattern centre [deg] (default: 22.0145)
566   - offset - Offset from pattern centre [deg] (default: 1.5)
  652 + Returns
  653 + -------
  654 + obsdeflist : list
  655 + Observation definition list
567 656 """
568 657 # Initialise observation definition list
569 658 obsdeflist = []
570 659  
571   - # Add patterns
572   - if pattern == "single":
  660 + # If the pattern is a single observation then append the Right Ascension
  661 + # and Declination to the observation definition list
  662 + if pattern == 'single':
573 663 obsdef = {'ra': ra, 'dec': dec}
574 664 obsdeflist.append(obsdef)
575   - elif pattern == "four":
  665 +
  666 + # ... otherwise, if the pattern is four observations then append four
  667 + # observations offset by a certain amount from the pattern centre to the
  668 + # observation definition list
  669 + elif pattern == 'four':
  670 +
576 671 # Set pattern centre
577 672 centre = gammalib.GSkyDir()
578 673 centre.radec_deg(ra, dec)
... ... @@ -582,8 +677,232 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
582 677 pntdir = centre.copy()
583 678 pntdir.rotate_deg(phi, offset)
584 679 obsdeflist.append({'ra': pntdir.ra_deg(), 'dec': pntdir.dec_deg()})
  680 +
  681 + # ... otherwise we have an unknown pattern
585 682 else:
586   - print("Warning: Observation pattern '"+str(pattern)+"' not recognized.")
  683 + print('Warning: Observation pattern "'+str(pattern)+'" not recognized.')
587 684  
588 685 # Return observation definition list
589 686 return obsdeflist
  687 +
  688 +
  689 +# ======================== #
  690 +# Set model for TS fitting #
  691 +# ======================== #
  692 +def set_ts_model(models, srcname, ra=None, dec=None, fitspat=False, fitspec=False):
  693 + """
  694 + Set model for TS fitting
  695 +
  696 + Parameters
  697 + ----------
  698 + models : `~gammalib.GModels`
  699 + Input model container
  700 + srcname : str
  701 + Test source name
  702 + ra : float, optional
  703 + Right Ascension of test source (deg)
  704 + dec : float, optional
  705 + Declination of test source (deg)
  706 + fitspat : bool, optional
  707 + Fit spatial parameter?
  708 + fitspec : bool, optional
  709 + Fit spectral parameters?
  710 +
  711 + Returns
  712 + -------
  713 + model : `~gammalib.GModels`
  714 + Model container for TS fitting
  715 + """
  716 + # Create a clone of the input model
  717 + outmodels = models.copy()
  718 +
  719 + # Disable TS computation for all model components (will enable the
  720 + # test source later)
  721 + for model in outmodels:
  722 + model.tscalc(False)
  723 +
  724 + # Get source model and enable TS computation
  725 + model = outmodels[srcname]
  726 + model.tscalc(True)
  727 +
  728 + # If source model has no "Prefactor" parameter then raise an exception
  729 + if not model.has_par('Prefactor'):
  730 + msg = ('Model "%s" has no parameter "Prefactor". Only spectral '
  731 + 'models with a "Prefactor" parameter are supported.' %
  732 + srcname)
  733 + raise RuntimeError(msg)
  734 +
  735 + # Set position of test source
  736 + if ra != None and dec != None:
  737 + if model.has_par('RA') and model.has_par('DEC'):
  738 + model['RA'].value(ra)
  739 + model['DEC'].value(dec)
  740 +
  741 + # Set possible spatial and spectral parameters
  742 + spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA',
  743 + 'MinorRadius', 'MajorRadius']
  744 + spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy',
  745 + 'InverseCutoffEnergy']
  746 +
  747 + # Fit or fix spatial parameters
  748 + for par in spatial:
  749 + if model.has_par(par):
  750 + if fitspat:
  751 + model[par].free()
  752 + else:
  753 + model[par].fix()
  754 +
  755 + # Fit or fix spectral parameters
  756 + for par in spectral:
  757 + if model.has_par(par):
  758 + if fitspec:
  759 + model[par].free()
  760 + else:
  761 + model[par].fix()
  762 +
  763 + # Return model container
  764 + return outmodels
  765 +
  766 +
  767 +# ==================== #
  768 +# Get stacked response #
  769 +# ==================== #
  770 +def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
  771 + emin=0.1, emax=100.0, enumbins=20, edisp=False,
  772 + coordsys='GAL', proj='TAN',
  773 + log=False, debug=False, chatter=2):
  774 + """
  775 + Get stacked response cubes
  776 +
  777 + Parameters
  778 + ----------
  779 + obs : `~gammalib.GObservations`
  780 + Observation container
  781 + xref : float
  782 + Right Ascension or Galactic longitude of response centre (deg)
  783 + yref : float
  784 + Declination or Galactic latitude of response centre (deg)
  785 + binsz : float, optional
  786 + Pixel size (deg/pixel)
  787 + nxpix : int, optional
  788 + Number of pixels in X direction
  789 + nypix : int, optional
  790 + Number of pixels in Y direction
  791 + emin : float, optional
  792 + Minimum energy (TeV)
  793 + emax : float, optional
  794 + Maximum energy (TeV)
  795 + enumbins : int, optional
  796 + Number of energy bins
  797 + edisp : bool, optional
  798 + Apply energy dispersion?
  799 + coordsys : str, optional
  800 + Coordinate system
  801 + proj : str, optional
  802 + Projection
  803 + log : bool, optional
  804 + Create log file(s)
  805 + debug : bool, optional
  806 + Create console dump?
  807 + chatter : int, optional
  808 + Chatter level
  809 +
  810 + Returns
  811 + -------
  812 + result : dict
  813 + Dictionary of response cubes
  814 + """
  815 + # Compute spatial binning for point spread function and
  816 + # energy dispersion cubes
  817 + psf_binsz = 10.0 * binsz
  818 + psf_nxpix = min(nxpix // 10, 2) # Make sure result is int
  819 + psf_nypix = min(nypix // 10, 2) # Make sure result is int
  820 +
  821 + # Create exposure cube
  822 + expcube = ctools.ctexpcube(obs)
  823 + expcube['incube'] = 'NONE'
  824 + expcube['usepnt'] = False
  825 + expcube['ebinalg'] = 'LOG'
  826 + expcube['xref'] = xref
  827 + expcube['yref'] = yref
  828 + expcube['binsz'] = binsz
  829 + expcube['nxpix'] = nxpix
  830 + expcube['nypix'] = nypix
  831 + expcube['enumbins'] = enumbins
  832 + expcube['emin'] = emin
  833 + expcube['emax'] = emax
  834 + expcube['coordsys'] = coordsys
  835 + expcube['proj'] = proj
  836 + if log:
  837 + expcube.logFileOpen()
  838 + expcube.run()
  839 +
  840 + # Create point spread function cube
  841 + psfcube = ctools.ctpsfcube(obs)
  842 + psfcube['incube'] = 'NONE'
  843 + psfcube['usepnt'] = False
  844 + psfcube['ebinalg'] = 'LOG'
  845 + psfcube['xref'] = xref
  846 + psfcube['yref'] = yref
  847 + psfcube['binsz'] = psf_binsz
  848 + psfcube['nxpix'] = psf_nxpix
  849 + psfcube['nypix'] = psf_nypix
  850 + psfcube['enumbins'] = enumbins
  851 + psfcube['emin'] = emin
  852 + psfcube['emax'] = emax
  853 + psfcube['coordsys'] = coordsys
  854 + psfcube['proj'] = proj
  855 + if log:
  856 + psfcube.logFileOpen()
  857 + psfcube.run()
  858 +
  859 + # Create background cube
  860 + bkgcube = ctools.ctbkgcube(obs)
  861 + bkgcube['incube'] = 'NONE'
  862 + bkgcube['usepnt'] = False
  863 + bkgcube['ebinalg'] = 'LOG'
  864 + bkgcube['xref'] = xref
  865 + bkgcube['yref'] = yref
  866 + bkgcube['binsz'] = binsz
  867 + bkgcube['nxpix'] = nxpix
  868 + bkgcube['nypix'] = nypix
  869 + bkgcube['enumbins'] = enumbins
  870 + bkgcube['emin'] = emin
  871 + bkgcube['emax'] = emax
  872 + bkgcube['coordsys'] = coordsys
  873 + bkgcube['proj'] = proj
  874 + if log:
  875 + bkgcube.logFileOpen()
  876 + bkgcube.run()
  877 +
  878 + # If energy dispersion is requested then create energy dispersion cube
  879 + if edisp:
  880 + edispcube = ctools.ctedispcube(obs)
  881 + edispcube['incube'] = 'NONE'
  882 + edispcube['usepnt'] = False
  883 + edispcube['ebinalg'] = 'LOG'
  884 + edispcube['xref'] = xref
  885 + edispcube['yref'] = yref
  886 + edispcube['binsz'] = psf_binsz
  887 + edispcube['nxpix'] = psf_nxpix
  888 + edispcube['nypix'] = psf_nypix
  889 + edispcube['enumbins'] = enumbins
  890 + edispcube['emin'] = emin
  891 + edispcube['emax'] = emax
  892 + edispcube['coordsys'] = coordsys
  893 + edispcube['proj'] = proj
  894 + if log:
  895 + edispcube.logFileOpen()
  896 + edispcube.run()
  897 +
  898 + # Build response dictionary
  899 + response = {}
  900 + response['expcube'] = expcube.expcube().copy()
  901 + response['psfcube'] = psfcube.psfcube().copy()
  902 + response['bkgcube'] = bkgcube.bkgcube().copy()
  903 + response['models'] = bkgcube.models().copy()
  904 + if edisp:
  905 + response['edispcube'] = edispcube.expcube().copy()
  906 +
  907 + # Return response cubes
  908 + return response
... ...
doc/source/users/reference_manual/cstsdist.rst
  View file @8085434
... ... @@ -75,6 +75,12 @@ General parameters
75 75 ``dec [real]``
76 76 Declination of CTA pointing (J2000, in degrees).
77 77  
  78 +``coordsys <CEL|GAL> [string]``
  79 + Coordinate system (CEL - celestial, GAL - galactic).
  80 +
  81 +``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
  82 + Projection method.
  83 +
78 84 ``emin [real]``
79 85 Lower energy limit (in TeV).
80 86  
... ...
test/test_cstsdist.py
  View file @8085434
... ... @@ -76,7 +76,7 @@ class Test(test):
76 76 ' enumbins=0 tmax=1800.0 rad=5.0 npix=200 npix=200'+ \
77 77 ' binsz=0.05'+ \
78 78 ' outfile="cstsdist_cmd1.dat"'+ \
79   - ' logfile="cstsdist_cmd1.log" chatter=1'
  79 + ' logfile="cstsdist_cmd1.log" chatter=2'
80 80  
81 81 # Check if execution of wrong command fails
82 82 self.test_assert(self._execute('command_that_does_not_exist') != 0,
... ... @@ -117,7 +117,6 @@ class Test(test):
117 117 tsdist['srcname'] = 'Crab'
118 118 tsdist['caldb'] = self._caldb
119 119 tsdist['irf'] = self._irf
120   - tsdist['outfile'] = 'cstsdist_py1.dat'
121 120 tsdist['ntrials'] = 1
122 121 tsdist['ra'] = 83.63
123 122 tsdist['dec'] = 22.01
... ... @@ -126,10 +125,9 @@ class Test(test):
126 125 tsdist['enumbins'] = 0
127 126 tsdist['tmax'] = 1800.0
128 127 tsdist['rad'] = 5.0
129   - tsdist['npix'] = 200
130   - tsdist['binsz'] = 0.05
  128 + tsdist['outfile'] = 'cstsdist_py1.dat'
131 129 tsdist['logfile'] = 'cstsdist_py1.log'
132   - tsdist['chatter'] = 2
  130 + tsdist['chatter'] = 3
133 131  
134 132 # Run cstsdist script
135 133 tsdist.logFileOpen() # Make sure we get a log file
... ... @@ -139,20 +137,56 @@ class Test(test):
139 137 # Check pull distribution file
140 138 self._check_result_file('cstsdist_py1.dat')
141 139  
  140 + # Set-up binned cstsdist in four observation wobble mode
  141 + tsdist = cscripts.cstsdist()
  142 + tsdist['inmodel'] = self._model
  143 + tsdist['srcname'] = 'Crab'
  144 + tsdist['caldb'] = self._caldb
  145 + tsdist['irf'] = self._irf
  146 + tsdist['ntrials'] = 1
  147 + tsdist['ra'] = 83.63
  148 + tsdist['dec'] = 22.01
  149 + tsdist['coordsys'] = 'CEL'
  150 + tsdist['proj'] = 'TAN'
  151 + tsdist['emin'] = 0.1
  152 + tsdist['emax'] = 100.0
  153 + tsdist['enumbins'] = 10
  154 + tsdist['tmax'] = 1800.0
  155 + tsdist['rad'] = 5.0
  156 + tsdist['npix'] = 50
  157 + tsdist['binsz'] = 0.1
  158 + tsdist['pattern'] = 'four'
  159 + tsdist['outfile'] = 'cstsdist_py2.dat'
  160 + tsdist['logfile'] = 'cstsdist_py2.log'
  161 + tsdist['chatter'] = 4
  162 +
  163 + # Execute cstsdist script
  164 + tsdist.execute()
  165 +
  166 + # Check pull distribution file
  167 + self._check_result_file('cstsdist_py2.dat', ncols=9)
  168 +
142 169 # Return
143 170 return
144 171  
145 172 # Check result file
146   - def _check_result_file(self, filename):
  173 + def _check_result_file(self, filename, ncols=11):
147 174 """
148 175 Check result file
  176 +
  177 + Parameters
  178 + ----------
  179 + filename : str
  180 + Name of result file
  181 + ncols : int, optional
  182 + Required number of columns
149 183 """
150 184 # Open result file as CSV file
151 185 results = gammalib.GCsv(filename, ',')
152 186  
153 187 # Check dimensions
154   - self.test_value(results.nrows(), 2, 'Check for 2 rows in TS file')
155   - self.test_value(results.ncols(), 14, 'Check for 14 columns in TS file')
  188 + self.test_value(results.nrows(), 2, 'Check rows in TS file')
  189 + self.test_value(results.ncols(), ncols, 'Check columns in TS file')
156 190  
157 191 # Return
158 192 return
... ...