Konstancja Satalecka / ctools

Commit ac3d2d26bd2b41c6dc79d95e08a9000566caa441

Authored by Jürgen Knödlseder
1 parent bb3fc5da

Change interface. Now use a specific source from the inout model.

Inline Side-by-side
Showing 2 changed files with 182 additions and 149 deletions
doc/source/reference_manual/cstsdist.rst
  View file @ac3d2d2
... ... @@ -16,30 +16,49 @@ Monte-Carlo simulations.
16 16 General parameters
17 17 ------------------
18 18  
19   -``(outfile = ts.dat) [file]``
  19 +``(inobs = NONE) [file]``
  20 + Input event list, counts cube or observation definition XML file.
  21 +
  22 +``inmodel [file]``
  23 + Input source model XML file.
  24 +
  25 +``srcname [string]``
  26 + Name of the source in the source model XML file which should be used
  27 + for sensitivity computation.
  28 +
  29 +``outfile [file]``
20 30 ASCII file containing the TS distribution values.
21   -
22   -``ntrials [integer]``
23   - Number of Monte Carlo samples.
24   -
  31 +
  32 +``expcube [file]``
  33 + Exposure cube file (only needed for stacked analysis).
  34 +
  35 +``psfcube [file]``
  36 + PSF cube file (only needed for stacked analysis).
  37 +
  38 +``bkgcube [file]``
  39 + Background cube file (only needed for stacked analysis).
  40 +
25 41 ``caldb [string]``
26 42 Calibration database.
27 43  
28 44 ``irf [string]``
29 45 Instrumental response function.
30 46  
31   -``type <point|gauss|shell|disk> [string]``
32   - Source model type.
33   -
34   -``(index = -2.48) [real]``
35   - Spectral index for source model.
36   -
37   -``offset [real]``
38   - Source offset angle (in degrees).
39   -
40   -``bkg [string]``
41   - Background model file function (none=power law for subarray E).
  47 +``(deadc = 0.95) [real]``
  48 + Average deadtime correction factor.
  49 +
  50 +``(edisp = no) [boolean]``
  51 + Apply energy dispersion to response computation.
  52 +
  53 +``ntrials [integer]``
  54 + Number of Monte Carlo samples.
  55 +
  56 +``ra [real]``
  57 + Right Ascension of CTA pointing (J2000, in degrees).
42 58  
  59 +``dec [real]``
  60 + Declination of CTA pointing (J2000, in degrees).
  61 +
43 62 ``emin [real]``
44 63 Lower energy limit (in TeV).
45 64  
... ... @@ -49,18 +68,27 @@ General parameters
49 68 ``enumbins [integer]``
50 69 Number of energy bins (0=unbinned).
51 70  
52   -``duration [real]``
53   - Effective exposure time (in seconds).
  71 +``(tmin = 0.0) [real]``
  72 + Start time (in seconds).
54 73  
55   -``(rad = 5.0) [real]``
56   - Radius of ROI (in degrees).
  74 +``tmax [real]``
  75 + Stop time (in seconds).
57 76  
58   -``(npix = 200) [integer]``
  77 +``npix [integer]``
59 78 Number of pixels for binned analysis.
60 79  
61   -``(binsz = 0.05) [real]``
  80 +``binsz [real]``
62 81 Pixel size for binned analysis.
63 82  
  83 +``(rad = 5.0) [real]``
  84 + Radius of ROI (in degrees).
  85 +
  86 +``(pattern = single) [string]``
  87 + Pattern for pointing simulation (single/four).
  88 +
  89 +``(offset = 1.5) [real]``
  90 + Observation pattern offset (in degrees).
  91 +
64 92  
65 93 Standard parameters
66 94 -------------------
... ... @@ -86,6 +114,9 @@ Standard parameters
86 114 ``(mode = ql) [string]``
87 115 Mode of automatic parameters (default is "ql", i.e. "query and learn").
88 116  
  117 +``(logfile = cstsdist.log) [string]``
  118 + Log filename.
  119 +
89 120  
90 121 Related tools
91 122 -------------
... ...
scripts/cstsdist.py
  View file @ac3d2d2
... ... @@ -95,27 +95,29 @@ class cstsdist(ctools.cscript):
95 95 # Create default parfile
96 96 pars = gammalib.GApplicationPars()
97 97 pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file"))
98   - pars.append(gammalib.GApplicationPar("outfile","f","h","ts.dat","","","Output file name"))
  98 + pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
  99 + pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
  100 + pars.append(gammalib.GApplicationPar("outfile","f","a","ts.dat","","","Output file name"))
99 101 pars.append(gammalib.GApplicationPar("expcube","s","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
100 102 pars.append(gammalib.GApplicationPar("psfcube","s","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
  103 + pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)"))
101 104 pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database"))
102 105 pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function"))
103   - pars.append(gammalib.GApplicationPar("ra","r","h","83.63","0","360","RA of pointing (deg)"))
104   - pars.append(gammalib.GApplicationPar("dec","r","h","22.01","-90","90","Dec of pointing (deg)"))
  106 + pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))
  107 + pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
  108 + pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
  109 + pars.append(gammalib.GApplicationPar("ra","r","a","83.63","0","360","RA of pointing (deg)"))
  110 + pars.append(gammalib.GApplicationPar("dec","r","a","22.01","-90","90","Dec of pointing (deg)"))
105 111 pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit (TeV)"))
106 112 pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit (TeV)"))
  113 + pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
107 114 pars.append(gammalib.GApplicationPar("tmin","r","h","0.0","","","Start time (MET in s)"))
108 115 pars.append(gammalib.GApplicationPar("tmax","r","a","1800.0","","","Observation duration (in s)"))
109   - pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
110   - pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
111   - pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
112   - pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))
  116 + pars.append(gammalib.GApplicationPar("npix","i","a","200","","","Number of pixels for binned"))
  117 + pars.append(gammalib.GApplicationPar("binsz","r","a","0.05","","","Pixel size for binned (deg/pixel)"))
113 118 pars.append(gammalib.GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
114   - pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
115   - pars.append(gammalib.GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
116   - pars.append(gammalib.GApplicationPar("index","r","h","-2.48","","","Spectral index"))
117   - pars.append(gammalib.GApplicationPar("offset","r","a","0.0","","","Source offset angle (deg)"))
118   - pars.append(gammalib.GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function"))
  119 + pars.append(gammalib.GApplicationPar("pattern","s","h","single","","","Observation pattern (single/four)"))
  120 + pars.append(gammalib.GApplicationPar("offset","r","h","1.5","","","Observation pattern offset (deg)"))
119 121 pars.append_standard()
120 122 pars.append(gammalib.GApplicationPar("logfile","f","h","cstsdist.log","","","Log filename"))
121 123 pars.save(parfile)
... ... @@ -131,6 +133,15 @@ class cstsdist(ctools.cscript):
131 133 # Set observation if not done before
132 134 if self.obs == None or self.obs.size() == 0:
133 135 self.obs = self.get_observations()
  136 +
  137 + # Check for requested pattern and use above
  138 + # observation parameters to set wobble pattern
  139 + self.m_pattern = self["pattern"].string()
  140 + if self.m_pattern == "four":
  141 + self.obs = self.set_obs()
  142 +
  143 + # Get source name
  144 + self.m_srcname = self["srcname"].string()
134 145  
135 146 # Get number of energy bins
136 147 self.m_enumbins = self["enumbins"].integer()
... ... @@ -143,36 +154,21 @@ class cstsdist(ctools.cscript):
143 154 # Set dummy values (required by obsutils)
144 155 self.m_npix = 0
145 156 self.m_binsz = 0.0
  157 +
  158 + # Set models if we have none
  159 + if self.obs.models().size() == 0:
  160 + self.obs.models(self["inmodel"].filename())
146 161  
147 162 # Get other parameters
148   - self.m_outfile = self["outfile"].filename()
149   - self.m_ntrials = self["ntrials"].integer()
150   - self.m_type = self["type"].string()
151   - self.m_index = self["index"].real()
152   - self.m_offset = self["offset"].real()
153   - self.m_bkg = self["bkg"].string()
  163 + self.m_outfile = self["outfile"].filename()
  164 + self.m_ntrials = self["ntrials"].integer()
  165 + self.m_edisp = self["edisp"].boolean()
  166 + self.m_offset = self["offset"].real()
154 167  
155 168 # Set some fixed parameters
156   - self.m_log = False # Logging in client tools
157   - self.m_debug = False # Debugging in client tools
158   -
159   - pnt = gammalib.GSkyDir()
160   - pnt.radec_deg(self["ra"].real(),self["dec"].real())
161   -
162   - # Initialise models. Note that we centre the point source at the
163   - # center of our observation, so we're onaxis.
164   - self.bkg_model = gammalib.GModels()
165   - self.full_model = gammalib.GModels()
166   - self.bkg_model.append(self.set_bkg_model())
167   - self.full_model.append(self.set_bkg_model())
168   - self.full_model.append(self.set_src_model(pnt.l_deg(), pnt.b_deg()+self.m_offset, \
169   - flux=0.010, \
170   - type=self.m_type, \
171   - index=self.m_index))
172   -
173   - # Attach background model to observation container
174   - self.obs.models(self.bkg_model)
175   -
  169 + self.m_log = False # Logging in client tools
  170 + self.m_debug = self["debug"].boolean() # Debugging in client tools
  171 +
176 172 # Return
177 173 return
178 174  
... ... @@ -211,6 +207,20 @@ class cstsdist(ctools.cscript):
211 207 if self.logTerse():
212 208 self.log_parameters()
213 209 self.log("\n")
  210 +
  211 + # Initialise models
  212 + full_model, bkg_model = self.set_models()
  213 +
  214 + # Write models into logger
  215 + if self.logTerse():
  216 + self.log("\n")
  217 + self.log.header1("Models")
  218 + self.log.header2("Background model")
  219 + self.log(str(bkg_model))
  220 + self.log("\n\n")
  221 + self.log.header2("Full model")
  222 + self.log(str(full_model))
  223 + self.log("\n")
214 224  
215 225 # Write observation into logger
216 226 if self.logTerse():
... ... @@ -235,7 +245,7 @@ class cstsdist(ctools.cscript):
235 245 for seed in range(self.m_ntrials):
236 246  
237 247 # Make a trial
238   - result = self.trial(seed)
  248 + result = self.trial(seed, full_model, bkg_model)
239 249  
240 250 # Write out result immediately
241 251 if seed == 0:
... ... @@ -251,97 +261,89 @@ class cstsdist(ctools.cscript):
251 261 # Return
252 262 return
253 263  
254   - def set_bkg_model(self, fitsigma=False):
  264 + def set_models(self, fitpos=False, fitspec=False):
255 265 """
256   - Setup CTA background model.
  266 + Set full and background model.
257 267 """
258   - # Define radial component
259   - radial = gammalib.GCTAModelRadialGauss(3.0)
260   - if fitsigma:
261   - radial["Sigma"].free()
262   - else:
263   - radial["Sigma"].fix()
264   -
265   - # Define spectral component
266   - spectrum = gammalib.GModelSpectralFunc(self.m_bkg, 1.0)
  268 + # Retrieve full model from observation container
  269 + full_model = self.obs.models().copy()
  270 +
  271 + # Get source model
  272 + model = full_model[self.m_srcname]
267 273  
268   - # Create background model
269   - model = gammalib.GCTAModelRadialAcceptance(radial, spectrum)
270   - model.name("Background")
271   - model.instruments("CTA")
272   -
273   - # Return background model
274   - return model
275   -
276   - def set_src_model(self, l, b, flux=1.0, index=-2.48, \
277   - type="point", sigma=1.0, radius=1.0, width=0.1, \
278   - fitpos=False, fitidx=False):
279   - """
280   - Returns a single source with Crab-like spectrum. The source flux
281   - can be scaled in Crab units. The Crab spectrum is based on MAGIC
282   - observations (Albert et al. 2008, ApJ, 674, 1037).
  274 + # Check that model has a Prefactor
  275 + if not model.has_par("Prefactor"):
  276 + msg = "Model \""+self.m_srcname+"\" has no parameter \"Prefactor\"."+ \
  277 + " Only spectral models with a \"Prefactor\" parameter are supported."
  278 + raise gammalib.GException.invalid_value("cssens", msg)
283 279  
284   - Parameters:
285   - l - Galactic longitude of source location [deg]
286   - b - Galactic latitude of source location [deg]
287   - Keywords:
288   - flux - Source flux [Crabs]
289   - index - Spectral index
290   - type - Source type ("point", "gauss", "disk", "shell")
291   - sigma - Gaussian sigma (for type="gauss")
292   - radius - Disk or shell inner radius [deg] (for type="disk" and type="shell")
293   - width - Shell width [deg] (for type="shell")
294   - fitpos - Fit position? (default: True)
295   - fitidx - Fit index? (default: True)
296   - """
297   - # Set source location
298   - location = gammalib.GSkyDir()
299   - location.lb_deg(l, b)
300   -
301   - # Set source spectrum
302   - spectrum = gammalib.GModelSpectralPlaw(flux*5.7e-16, index, gammalib.GEnergy(0.3, "TeV"))
303   - if fitidx:
304   - spectrum["Index"].free()
  280 + # Fit or fix spatial parameters
  281 + if fitpos:
  282 + if model.has_par("RA"):
  283 + model["RA"].free()
  284 + if model.has_par("DEC"):
  285 + model["DEC"].free()
  286 + if model.has_par("Sigma"):
  287 + model["Sigma"].free()
  288 + if model.has_par("Radius"):
  289 + model["Radius"].free()
  290 + if model.has_par("Width"):
  291 + model["Width"].free()
305 292 else:
306   - spectrum["Index"].fix()
307   -
308   - # Set source
309   - if type == "point":
310   - spatial = gammalib.GModelSpatialPointSource(location)
311   - if fitpos:
312   - spatial[0].free()
313   - spatial[1].free()
314   - elif type == "gauss":
315   - spatial = gammalib.GModelSpatialRadialGauss(location, sigma)
316   - if fitpos:
317   - spatial[0].free()
318   - spatial[1].free()
319   - elif type == "disk":
320   - spatial = gammalib.GModelSpatialRadialDisk(location, radius)
321   - if fitpos:
322   - spatial[0].free()
323   - spatial[1].free()
324   - elif type == "shell":
325   - spatial = gammalib.GModelSpatialRadialShell(location, radius, width)
326   - if fitpos:
327   - spatial[0].free()
328   - spatial[1].free()
  293 + if model.has_par("RA"):
  294 + model["RA"].fix()
  295 + if model.has_par("DEC"):
  296 + model["DEC"].fix()
  297 + if model.has_par("Sigma"):
  298 + model["Sigma"].fix()
  299 + if model.has_par("Radius"):
  300 + model["Radius"].fix()
  301 + if model.has_par("Width"):
  302 + model["Width"].fix()
  303 +
  304 + # Fit or fix spectral parameters
  305 + if fitspec:
  306 + if model.has_par("Index"):
  307 + model["Index"].free()
  308 + if model.has_par("Cutoff"):
  309 + model["Cutoff"].free()
329 310 else:
330   - self.log("ERROR: Unknown source type '"+type+"'.\n")
331   - return None
332   - source = gammalib.GModelSky(spatial, spectrum)
  311 + if model.has_par("Index"):
  312 + model["Index"].fix()
  313 + if model.has_par("Cutoff"):
  314 + model["Cutoff"].fix()
333 315  
334   - # Set source name
335   - source.name("Test")
  316 + # Create background model
  317 + bkg_model = full_model.copy()
  318 + bkg_model.remove(self.m_srcname)
  319 +
  320 + # Return models
  321 + return full_model, bkg_model
  322 +
  323 + def set_obs(self):
  324 + """
  325 + Returns an observation container with a set of CTA observations.
336 326  
337   - # Turn TS value computation off
338   - # since we compute it by hand later on
339   - source.tscalc(False)
340   -
341   - # Return source
342   - return source
  327 + Keywords:
  328 + """
  329 +
  330 + # Setup observation definition list
  331 + obsdeflist = obsutils.set_obs_patterns(self.m_pattern, \
  332 + ra=self["ra"].real(), dec=self["dec"].real(), \
  333 + offset=self["offset"].real())
  334 +
  335 + # Create list of observations
  336 + obs = obsutils.set_obs_list(obsdeflist, \
  337 + tstart=self["tmin"].real(), duration=self["tmax"].real()-self["tmin"].real(), \
  338 + deadc=self["deadc"].real(), \
  339 + emin=self["emin"].real(), emax=self["emax"].real(), \
  340 + rad=self["rad"].real(), \
  341 + irf=self["irf"].string(), caldb=self["caldb"].string())
343 342  
344   - def trial(self, seed):
  343 + # Return observation container
  344 + return obs
  345 +
  346 + def trial(self, seed, full_model, bkg_model):
345 347 """
346 348 Create the TS for a single trial.
347 349  
... ... @@ -373,9 +375,9 @@ class cstsdist(ctools.cscript):
373 375 self.log("\n")
374 376  
375 377 # Fit background only
376   - sim.models(self.bkg_model)
  378 + sim.models(bkg_model)
377 379 like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
378   - result_bgm = like_bgm.obs().models()
  380 + result_bgm = like_bgm.obs().models().copy()
379 381 LogL_bgm = like_bgm.opt().value()
380 382 npred_bgm = like_bgm.obs().npred()
381 383  
... ... @@ -396,9 +398,9 @@ class cstsdist(ctools.cscript):
396 398 self.log(str(par)+"\n")
397 399  
398 400 # Fit background and test source
399   - sim.models(self.full_model)
  401 + sim.models(full_model)
400 402 like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
401   - result_all = like_all.obs().models()
  403 + result_all = like_all.obs().models().copy()
402 404 LogL_all = like_all.opt().value()
403 405 npred_all = like_all.obs().npred()
404 406 ts = 2.0*(LogL_bgm-LogL_all)
... ... @@ -428,9 +430,9 @@ class cstsdist(ctools.cscript):
428 430 self.log("TS=")
429 431 self.log(ts)
430 432 self.log(" Prefactor=")
431   - self.log(result_all["Test"]["Prefactor"].value())
  433 + self.log(result_all[self.m_srcname]["Prefactor"].value())
432 434 self.log("+/-")
433   - self.log(result_all["Test"]["Prefactor"].error())
  435 + self.log(result_all[self.m_srcname]["Prefactor"].error())
434 436 self.log("\n")
435 437  
436 438 # Initialise results
... ...