Dmitry Malyshev / ctools

Commit d9a9c73fc2157eb85560603620092b56416433c1

Authored by Jürgen Knödlseder
1 parent a448964c

Convert tabs to spaces.

Inline Side-by-side
Showing 5 changed files with 1262 additions and 1262 deletions
scripts/cssens.py
  View file @d9a9c73
... ... @@ -17,7 +17,7 @@
17 17 #
18 18 # The source is modelled as a point source
19 19 #
20   -# Copyright (C) 2011 Jurgen Knodlseder
  20 +# Copyright (C) 2011-2014 Juergen Knoedlseder
21 21 #
22 22 # This program is free software: you can redistribute it and/or modify
23 23 # it under the terms of the GNU General Public License as published by
... ... @@ -45,682 +45,682 @@ import math
45 45 # cssens class #
46 46 # ============ #
47 47 class cssens(GApplication):
48   - """
49   - This class implements the sensitivity computation script. It derives
50   - from the GammaLib::GApplication class which provides support for parameter
51   - files, command line arguments, and logging. In that way the Python
52   - script behaves just as a regular ctool.
53   - """
54   - def __init__(self, *argv):
55   - """
56   - Constructor.
57   - """
58   - # Set name
59   - self.name = "cssens"
60   - self.version = "0.3.0"
61   -
62   - # Make sure that parfile exists
63   - file = self.parfile()
64   -
65   - # Initialise application
66   - if len(argv) == 0:
67   - GApplication.__init__(self, self.name, self.version)
68   - elif len(argv) ==1:
69   - GApplication.__init__(self, self.name, self.version, *argv)
70   - else:
71   - raise TypeError("Invalid number of arguments given.")
72   -
73   - # Set logger properties
74   - self.logFileOpen()
75   - self.log_header()
76   - self.log.date(True)
77   -
78   - # Return
79   - return
80   -
81   - def __del__(self):
82   - """
83   - Destructor.
84   - """
85   - # Write separator into logger
86   - if self.logTerse():
87   - self.log("\n")
88   -
89   - # Return
90   - return
91   -
92   - def parfile(self):
93   - """
94   - Check if parfile exists. If parfile does not exist then create a
95   - default parfile. This kluge avoids shipping the cscript with a parfile.
96   - """
97   - # Set parfile name
98   - parfile = self.name+".par"
99   -
100   - try:
101   - pars = GApplicationPars(parfile)
102   - except:
103   - # Signal if parfile was not found
104   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
105   -
106   - # Create default parfile
107   - pars = GApplicationPars()
108   - pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
109   - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
110   - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
111   - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
112   - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
113   - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
114   - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
115   - pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
116   - pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)"))
117   - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
118   - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
119   - pars.append(GApplicationPar("sigma","r","h","5.0","","","Significance threshold"))
120   - pars.append(GApplicationPar("ts_use","b","h","yes","","","Use TS?"))
121   - pars.append(GApplicationPar("index","r","h","-2.48","","","Assumed spectral index"))
122   - pars.append(GApplicationPar("radius","r","h","0.1","","","Extended source model radius"))
123   - pars.append(GApplicationPar("width","r","h","0.05","","","Extended source model width"))
124   - pars.append(GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations"))
125   - pars.append(GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
126   - pars.append_standard()
127   - pars.save(parfile)
128   -
129   - # Return
130   - return
131   -
132   - def get_parameters(self):
133   - """
134   - Get parameters from parfile.
135   - """
136   - # Get parameters
137   - self.m_outfile = self["outfile"].filename()
138   - self.m_caldb = self["caldb"].string()
139   - self.m_irf = self["irf"].string()
140   - self.m_type = self["type"].string()
141   - self.m_offset = self["offset"].real()
142   - self.m_bkg = self["bkg"].string()
143   - self.m_duration = self["duration"].real()
144   - self.m_roi = self["rad"].real()
145   - self.m_nbins = self["enumbins"].integer()
146   - self.m_npix = self["npix"].integer()
147   - self.m_binsz = self["binsz"].real()
148   - self.m_ts_thres = self["sigma"].real()*self["sigma"].real()
149   - self.m_use_ts = self["ts_use"].boolean()
150   - self.m_index = self["index"].real()
151   - self.m_radius = self["radius"].real()
152   - self.m_width = self["width"].real()
153   - self.m_max_iter = self["max_iter"].integer()
154   - self.m_num_avg = self["num_avg"].integer()
155   -
156   - # Set some fixed parameters
157   - self.m_log = False # Logging in client tools
158   - self.m_debug = False # Debugging in client tools
159   -
160   - # Return
161   - return
162   -
163   - def execute(self):
164   - """
165   - Execute the script.
166   - """
167   - # Run the script
168   - self.run()
169   -
170   - # Return
171   - return
172   -
173   - def run(self):
174   - """
175   - Run the script.
176   - """
177   - # Switch screen logging on in debug mode
178   - if self.logDebug():
179   - self.log.cout(True)
180   -
181   - # Get parameters
182   - self.get_parameters()
183   -
184   - # Write input parameters into logger
185   - if self.logTerse():
186   - self.log_parameters()
187   - self.log("\n")
188   -
189   - # Initialise script
190   - colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens']
191   - results = []
192   -
193   - # Initialise models. Note that we centre the point source at the Galactic
194   - # center as our observation is also centred at the Galactic centre, so
195   - # we're onaxis.
196   - bkg_model = GModels()
197   - full_model = GModels()
198   - bkg_model.append(self.set_bkg_model())
199   - full_model.append(self.set_bkg_model())
200   - full_model.append(self.set_src_model(0.0, \
201   - self.m_offset, \
202   - type=self.m_type, \
203   - index=self.m_index))
204   -
205   - # Write models into logger
206   - if self.logTerse():
207   - self.log("\n")
208   - self.log.header1("Models")
209   - self.log.header2("Background model")
210   - self.log(str(bkg_model))
211   - self.log("\n\n")
212   - self.log.header2("Full model")
213   - self.log(str(full_model))
214   - self.log("\n")
215   -
216   - # Write heder
217   - if self.logTerse():
218   - self.log("\n")
219   - self.log.header1("Sensitivity determination")
220   -
221   - # Loop over energy bands. The number of energy band is still
222   - # fixed ...
223   - for ieng in range(21):
224   -
225   - # Set energies
226   - loge = -1.7 + ieng * 0.2
227   - emean = pow(10.0, loge)
228   - emin = pow(10.0, loge-0.1)
229   - emax = pow(10.0, loge+0.1)
230   -
231   - # Setup observation(s)
232   - obs = self.set_obs(emin=emin, emax=emax)
233   -
234   - # Determine sensitivity
235   - result = self.get_sensitivity(obs, bkg_model, full_model)
236   -
237   - # Write results
238   - if ieng == 0:
239   - f = open(self.m_outfile, 'w')
240   - writer = csv.DictWriter(f, colnames)
241   - writer.writerow(dict((_,_) for _ in colnames))
242   - writer.writerow(result)
243   - f.close()
244   - else:
245   - f = open(self.m_outfile, 'a')
246   - writer = csv.DictWriter(f, colnames)
247   - writer.writerow(result)
248   - f.close()
249   -
250   - # Append results
251   - results.append(result)
252   -
253   - # Return
254   - return
255   -
256   - def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
257   - """
258   - Returns an observation container with a single CTA observation.
259   -
260   - Keywords:
261   - lpnt - Galactic longitude of pointing [deg] (default: 0.0)
262   - bpnt - Galactic latitude of pointing [deg] (default: 0.0)
263   - emin - Minimum energy [TeV] (default: 0.1)
264   - emax - Maximum energy [TeV] (default: 100.0)
265   - """
266   - # Allocate observation container
267   - obs = GObservations()
268   -
269   - # Set single pointing
270   - pntdir = GSkyDir()
271   - pntdir.lb_deg(lpnt, bpnt)
272   -
273   - # Create CTA observation
274   - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
275   - duration=self.m_duration, \
  48 + """
  49 + This class implements the sensitivity computation script. It derives
  50 + from the GammaLib::GApplication class which provides support for parameter
  51 + files, command line arguments, and logging. In that way the Python
  52 + script behaves just as a regular ctool.
  53 + """
  54 + def __init__(self, *argv):
  55 + """
  56 + Constructor.
  57 + """
  58 + # Set name
  59 + self.name = "cssens"
  60 + self.version = "0.3.0"
  61 +
  62 + # Make sure that parfile exists
  63 + file = self.parfile()
  64 +
  65 + # Initialise application
  66 + if len(argv) == 0:
  67 + GApplication.__init__(self, self.name, self.version)
  68 + elif len(argv) ==1:
  69 + GApplication.__init__(self, self.name, self.version, *argv)
  70 + else:
  71 + raise TypeError("Invalid number of arguments given.")
  72 +
  73 + # Set logger properties
  74 + self.logFileOpen()
  75 + self.log_header()
  76 + self.log.date(True)
  77 +
  78 + # Return
  79 + return
  80 +
  81 + def __del__(self):
  82 + """
  83 + Destructor.
  84 + """
  85 + # Write separator into logger
  86 + if self.logTerse():
  87 + self.log("\n")
  88 +
  89 + # Return
  90 + return
  91 +
  92 + def parfile(self):
  93 + """
  94 + Check if parfile exists. If parfile does not exist then create a
  95 + default parfile. This kluge avoids shipping the cscript with a parfile.
  96 + """
  97 + # Set parfile name
  98 + parfile = self.name+".par"
  99 +
  100 + try:
  101 + pars = GApplicationPars(parfile)
  102 + except:
  103 + # Signal if parfile was not found
  104 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  105 +
  106 + # Create default parfile
  107 + pars = GApplicationPars()
  108 + pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
  109 + pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  110 + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  111 + pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
  112 + pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
  113 + pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
  114 + pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  115 + pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
  116 + pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)"))
  117 + pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  118 + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  119 + pars.append(GApplicationPar("sigma","r","h","5.0","","","Significance threshold"))
  120 + pars.append(GApplicationPar("ts_use","b","h","yes","","","Use TS?"))
  121 + pars.append(GApplicationPar("index","r","h","-2.48","","","Assumed spectral index"))
  122 + pars.append(GApplicationPar("radius","r","h","0.1","","","Extended source model radius"))
  123 + pars.append(GApplicationPar("width","r","h","0.05","","","Extended source model width"))
  124 + pars.append(GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations"))
  125 + pars.append(GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
  126 + pars.append_standard()
  127 + pars.save(parfile)
  128 +
  129 + # Return
  130 + return
  131 +
  132 + def get_parameters(self):
  133 + """
  134 + Get parameters from parfile.
  135 + """
  136 + # Get parameters
  137 + self.m_outfile = self["outfile"].filename()
  138 + self.m_caldb = self["caldb"].string()
  139 + self.m_irf = self["irf"].string()
  140 + self.m_type = self["type"].string()
  141 + self.m_offset = self["offset"].real()
  142 + self.m_bkg = self["bkg"].string()
  143 + self.m_duration = self["duration"].real()
  144 + self.m_roi = self["rad"].real()
  145 + self.m_nbins = self["enumbins"].integer()
  146 + self.m_npix = self["npix"].integer()
  147 + self.m_binsz = self["binsz"].real()
  148 + self.m_ts_thres = self["sigma"].real()*self["sigma"].real()
  149 + self.m_use_ts = self["ts_use"].boolean()
  150 + self.m_index = self["index"].real()
  151 + self.m_radius = self["radius"].real()
  152 + self.m_width = self["width"].real()
  153 + self.m_max_iter = self["max_iter"].integer()
  154 + self.m_num_avg = self["num_avg"].integer()
  155 +
  156 + # Set some fixed parameters
  157 + self.m_log = False # Logging in client tools
  158 + self.m_debug = False # Debugging in client tools
  159 +
  160 + # Return
  161 + return
  162 +
  163 + def execute(self):
  164 + """
  165 + Execute the script.
  166 + """
  167 + # Run the script
  168 + self.run()
  169 +
  170 + # Return
  171 + return
  172 +
  173 + def run(self):
  174 + """
  175 + Run the script.
  176 + """
  177 + # Switch screen logging on in debug mode
  178 + if self.logDebug():
  179 + self.log.cout(True)
  180 +
  181 + # Get parameters
  182 + self.get_parameters()
  183 +
  184 + # Write input parameters into logger
  185 + if self.logTerse():
  186 + self.log_parameters()
  187 + self.log("\n")
  188 +
  189 + # Initialise script
  190 + colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens']
  191 + results = []
  192 +
  193 + # Initialise models. Note that we centre the point source at the Galactic
  194 + # center as our observation is also centred at the Galactic centre, so
  195 + # we're onaxis.
  196 + bkg_model = GModels()
  197 + full_model = GModels()
  198 + bkg_model.append(self.set_bkg_model())
  199 + full_model.append(self.set_bkg_model())
  200 + full_model.append(self.set_src_model(0.0, \
  201 + self.m_offset, \
  202 + type=self.m_type, \
  203 + index=self.m_index))
  204 +
  205 + # Write models into logger
  206 + if self.logTerse():
  207 + self.log("\n")
  208 + self.log.header1("Models")
  209 + self.log.header2("Background model")
  210 + self.log(str(bkg_model))
  211 + self.log("\n\n")
  212 + self.log.header2("Full model")
  213 + self.log(str(full_model))
  214 + self.log("\n")
  215 +
  216 + # Write heder
  217 + if self.logTerse():
  218 + self.log("\n")
  219 + self.log.header1("Sensitivity determination")
  220 +
  221 + # Loop over energy bands. The number of energy band is still
  222 + # fixed ...
  223 + for ieng in range(21):
  224 +
  225 + # Set energies
  226 + loge = -1.7 + ieng * 0.2
  227 + emean = pow(10.0, loge)
  228 + emin = pow(10.0, loge-0.1)
  229 + emax = pow(10.0, loge+0.1)
  230 +
  231 + # Setup observation(s)
  232 + obs = self.set_obs(emin=emin, emax=emax)
  233 +
  234 + # Determine sensitivity
  235 + result = self.get_sensitivity(obs, bkg_model, full_model)
  236 +
  237 + # Write results
  238 + if ieng == 0:
  239 + f = open(self.m_outfile, 'w')
  240 + writer = csv.DictWriter(f, colnames)
  241 + writer.writerow(dict((_,_) for _ in colnames))
  242 + writer.writerow(result)
  243 + f.close()
  244 + else:
  245 + f = open(self.m_outfile, 'a')
  246 + writer = csv.DictWriter(f, colnames)
  247 + writer.writerow(result)
  248 + f.close()
  249 +
  250 + # Append results
  251 + results.append(result)
  252 +
  253 + # Return
  254 + return
  255 +
  256 + def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
  257 + """
  258 + Returns an observation container with a single CTA observation.
  259 +
  260 + Keywords:
  261 + lpnt - Galactic longitude of pointing [deg] (default: 0.0)
  262 + bpnt - Galactic latitude of pointing [deg] (default: 0.0)
  263 + emin - Minimum energy [TeV] (default: 0.1)
  264 + emax - Maximum energy [TeV] (default: 100.0)
  265 + """
  266 + # Allocate observation container
  267 + obs = GObservations()
  268 +
  269 + # Set single pointing
  270 + pntdir = GSkyDir()
  271 + pntdir.lb_deg(lpnt, bpnt)
  272 +
  273 + # Create CTA observation
  274 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  275 + duration=self.m_duration, \
276 276 emin=emin, emax=emax, rad=self.m_roi)
277   -
278   - # Append observation to container
279   - obs.append(run)
280   -
281   - # Return observation container
282   - return obs
283   -
284   - def set_bkg_model(self, fitidx=False):
285   - """
286   - Setup standard CTA background model. We use a simple power law here,
287   - scaled to Konrad's E configuration performance table.
288   -
289   - Keywords:
290   - fitidx - Fit spectral index (default: False)
291   - """
292   - # Define radial component
293   - radial = GCTAModelRadialGauss(3.0)
294   -
295   - # Define spectral component. If a background model is given then
296   - # use a file function. Otherwise use a static power law.
297   - if len(self.m_bkg) > 0:
298   - spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
299   - else:
300   - spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV"))
301   - if fitidx:
302   - spectrum["Index"].free()
303   - else:
304   - spectrum["Index"].fix()
305   -
306   - # Create background model
307   - model = GCTAModelRadialAcceptance(radial, spectrum)
308   - model.name("Background")
309   - model.instruments("CTA")
310   -
311   - # Return background model
312   - return model
313   -
314   - def set_src_model(self, l, b, flux=1.0, index=-2.48, \
  277 +
  278 + # Append observation to container
  279 + obs.append(run)
  280 +
  281 + # Return observation container
  282 + return obs
  283 +
  284 + def set_bkg_model(self, fitidx=False):
  285 + """
  286 + Setup standard CTA background model. We use a simple power law here,
  287 + scaled to Konrad's E configuration performance table.
  288 +
  289 + Keywords:
  290 + fitidx - Fit spectral index (default: False)
  291 + """
  292 + # Define radial component
  293 + radial = GCTAModelRadialGauss(3.0)
  294 +
  295 + # Define spectral component. If a background model is given then
  296 + # use a file function. Otherwise use a static power law.
  297 + if len(self.m_bkg) > 0:
  298 + spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
  299 + else:
  300 + spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV"))
  301 + if fitidx:
  302 + spectrum["Index"].free()
  303 + else:
  304 + spectrum["Index"].fix()
  305 +
  306 + # Create background model
  307 + model = GCTAModelRadialAcceptance(radial, spectrum)
  308 + model.name("Background")
  309 + model.instruments("CTA")
  310 +
  311 + # Return background model
  312 + return model
  313 +
  314 + def set_src_model(self, l, b, flux=1.0, index=-2.48, \
315 315 type="point", fitpos=False, fitidx=False):
316   - """
317   - Returns a single source with Crab-like spectrum. The source flux
318   - can be scaled in Crab units. The Crab spectrum is based on MAGIC
319   - observations (Albert et al. 2008, ApJ, 674, 1037).
320   -
321   - Parameters:
322   - l - Galactic longitude of source location [deg]
323   - b - Galactic latitude of source location [deg]
324   - Keywords:
325   - flux - Source flux [Crabs]
326   - index - Spectral index
327   - type - Source type ("point", "gauss", "disk", "shell")
328   - fitpos - Fit position and size? (default: False)
329   - fitidx - Fit index? (default: False)
330   - """
331   - # Set source location
332   - location = GSkyDir()
333   - location.lb_deg(l, b)
334   -
335   - # Set source spectrum
336   - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")
337   - if fitidx:
338   - spectrum["Index"].free()
339   - else:
340   - spectrum["Index"].fix()
341   -
342   - # Set source
343   - if type == "point":
344   - spatial = GModelSpatialPointSource(location)
345   - if fitpos:
346   - spatial["RA"].free()
347   - spatial["DEC"].free()
348   - else:
349   - spatial["RA"].fix()
350   - spatial["DEC"].fix()
351   - elif type == "gauss":
352   - spatial = GModelSpatialRadialGauss(location, self.m_radius)
353   - if fitpos:
354   - spatial["RA"].free()
355   - spatial["DEC"].free()
356   - spatial["Sigma"].free()
357   - else:
358   - spatial["RA"].fix()
359   - spatial["DEC"].fix()
360   - spatial["Sigma"].fix()
361   - elif type == "disk":
362   - spatial = GModelSpatialRadialDisk(location, self.m_radius)
363   - if fitpos:
364   - spatial["RA"].free()
365   - spatial["DEC"].free()
366   - spatial["Radius"].free()
367   - else:
368   - spatial["RA"].fix()
369   - spatial["DEC"].fix()
370   - spatial["Radius"].fix()
371   - elif type == "shell":
372   - spatial = GModelSpatialRadialShell(location, self.m_radius, self.m_width)
373   - if fitpos:
374   - spatial["RA"].free()
375   - spatial["DEC"].free()
376   - spatial["Radius"].free()
377   - spatial["Width"].free()
378   - else:
379   - spatial["RA"].fix()
380   - spatial["DEC"].fix()
381   - spatial["Radius"].fix()
382   - spatial["Width"].fix()
383   - else:
384   - self.log("ERROR: Unknown source type '"+type+"'.\n")
385   - return None
  316 + """
  317 + Returns a single source with Crab-like spectrum. The source flux
  318 + can be scaled in Crab units. The Crab spectrum is based on MAGIC
  319 + observations (Albert et al. 2008, ApJ, 674, 1037).
  320 +
  321 + Parameters:
  322 + l - Galactic longitude of source location [deg]
  323 + b - Galactic latitude of source location [deg]
  324 + Keywords:
  325 + flux - Source flux [Crabs]
  326 + index - Spectral index
  327 + type - Source type ("point", "gauss", "disk", "shell")
  328 + fitpos - Fit position and size? (default: False)
  329 + fitidx - Fit index? (default: False)
  330 + """
  331 + # Set source location
  332 + location = GSkyDir()
  333 + location.lb_deg(l, b)
  334 +
  335 + # Set source spectrum
  336 + spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")
  337 + if fitidx:
  338 + spectrum["Index"].free()
  339 + else:
  340 + spectrum["Index"].fix()
  341 +
  342 + # Set source
  343 + if type == "point":
  344 + spatial = GModelSpatialPointSource(location)
  345 + if fitpos:
  346 + spatial["RA"].free()
  347 + spatial["DEC"].free()
  348 + else:
  349 + spatial["RA"].fix()
  350 + spatial["DEC"].fix()
  351 + elif type == "gauss":
  352 + spatial = GModelSpatialRadialGauss(location, self.m_radius)
  353 + if fitpos:
  354 + spatial["RA"].free()
  355 + spatial["DEC"].free()
  356 + spatial["Sigma"].free()
  357 + else:
  358 + spatial["RA"].fix()
  359 + spatial["DEC"].fix()
  360 + spatial["Sigma"].fix()
  361 + elif type == "disk":
  362 + spatial = GModelSpatialRadialDisk(location, self.m_radius)
  363 + if fitpos:
  364 + spatial["RA"].free()
  365 + spatial["DEC"].free()
  366 + spatial["Radius"].free()
  367 + else:
  368 + spatial["RA"].fix()
  369 + spatial["DEC"].fix()
  370 + spatial["Radius"].fix()
  371 + elif type == "shell":
  372 + spatial = GModelSpatialRadialShell(location, self.m_radius, self.m_width)
  373 + if fitpos:
  374 + spatial["RA"].free()
  375 + spatial["DEC"].free()
  376 + spatial["Radius"].free()
  377 + spatial["Width"].free()
  378 + else:
  379 + spatial["RA"].fix()
  380 + spatial["DEC"].fix()
  381 + spatial["Radius"].fix()
  382 + spatial["Width"].fix()
  383 + else:
  384 + self.log("ERROR: Unknown source type '"+type+"'.\n")
  385 + return None
386 386 source = GModelSky(spatial, spectrum)
387   -
388   - # Set source name
389   - source.name("Test")
390   -
391   - # Return source
392   - return source
393   -
394   - def get_sensitivity(self, obs, bkg_model, full_model):
395   - """
396   - Determine sensitivity for a given observations.
397   -
398   - Parameters:
399   - obs - Observation container
400   - bkg_model - Background model
401   - full_model - Source model
402   - """
403   - # Set TeV->erg conversion factor
404   - tev2erg = 1.6021764
405   -
406   - # Determine energy boundaries from first observation in
407   - # the container
408   - for run in obs:
409   - emin = run.events().ebounds().emin()
410   - emax = run.events().ebounds().emax()
411   - loge = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
412   - erg_mean = math.pow(10.0, loge) * tev2erg
413   - erg_width = (emax.TeV()-emin.TeV()) * tev2erg
414   - break
415   -
416   - # Write header
417   - if self.logTerse():
418   - self.log("\n")
419   - self.log.header2("Energies: "+str(emin)+" - "+str(emax))
420   -
421   - # Initialise loop
422   - flux_value = []
423   - pflux_value = []
424   - eflux_value = []
425   - diffSens_value = []
426   - iter = 0
427   - test_flux = 0.1 # This is the initial test flux in Crab units
428   -
429   - # Loop until we break
430   - while True:
431   -
432   - # Update iteration counter
433   - iter += 1
434   -
435   - # Write header
436   - if self.logExplicit():
437   - self.log.header2("Iteration "+str(iter))
438   -
439   - # Set test flux
440   - full_model["Test"]['Prefactor'].value(test_flux)
441   - obs.models(full_model)
442   -
443   - # Simulate events
444   - sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \
445   - binsz=self.m_binsz, npix=self.m_npix, \
446   - log=self.m_log, debug=self.m_debug)
447   -
448   - # Determine number of events in simulation
449   - nevents = 0.0
450   - for run in sim:
451   - nevents += run.events().number()
452   -
453   - # Write simulation results
454   - if self.logExplicit():
455   - self.log.header3("Simulation")
456   - self.log.parformat("Number of simulated events")
457   - self.log(nevents)
458   - self.log("\n")
459   -
460   - # Fit background only
461   - sim.models(bkg_model)
462   - like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
463   - result_bgm = like.obs().models().copy()
464   - LogL_bgm = like.opt().value()
465   - npred_bgm = like.obs().npred()
466   -
467   - # Assess quality based on a comparison between Npred and Nevents
468   - quality_bgm = npred_bgm-nevents
469   -
470   - # Write background fit results
471   - if self.logExplicit():
472   - self.log.header3("Background model fit")
473   - self.log.parformat("log likelihood")
474   - self.log(LogL_bgm)
475   - self.log("\n")
476   - self.log.parformat("Number of predicted events")
477   - self.log(npred_bgm)
478   - self.log("\n")
479   - self.log.parformat("Fit quality")
480   - self.log(quality_bgm)
481   - self.log("\n")
482   -
483   - # Start over if the fit quality was bad
484   - if abs(quality_bgm) > 3.0:
485   - if self.logExplicit():
486   - self.log("Fit quality outside required range. Start over.\n")
487   - continue
488   -
489   - # Write model fit results
490   - if self.logExplicit():
491   - for model in result_bgm:
492   - self.log.parformat("Model")
493   - self.log(model.name())
494   - self.log("\n")
495   - for par in model:
496   - self.log(str(par)+"\n")
497   -
498   - # Fit background and test source
499   - sim.models(full_model)
500   - like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
501   - result_all = like.obs().models().copy()
502   - LogL_all = like.opt().value()
503   - npred_all = like.obs().npred()
504   - ts = 2.0*(LogL_bgm-LogL_all)
505   -
506   - # Assess quality based on a comparison between Npred and Nevents
507   - quality_all = npred_all-nevents
508   -
509   - # Write background and test source fit results
510   - if self.logExplicit():
511   - self.log.header3("Background and test source model fit")
512   - self.log.parformat("Test statistics")
513   - self.log(ts)
514   - self.log("\n")
515   - self.log.parformat("log likelihood")
516   - self.log(LogL_all)
517   - self.log("\n")
518   - self.log.parformat("Number of predicted events")
519   - self.log(npred_all)
520   - self.log("\n")
521   - self.log.parformat("Fit quality")
522   - self.log(quality_all)
523   - self.log("\n")
524   - #
525   - for model in result_all:
526   - self.log.parformat("Model")
527   - self.log(model.name())
528   - self.log("\n")
529   - for par in model:
530   - self.log(str(par)+"\n")
531   -
532   - # Start over if the fit quality was bad
533   - if abs(quality_all) > 3.0:
534   - if self.logExplicit():
535   - self.log("Fit quality outside required range. Start over.\n")
536   - continue
537   -
538   - # Start over if TS was non-positive
539   - if ts <= 0.0:
540   - if self.logExplicit():
541   - self.log("Non positive TS. Start over.\n")
542   - continue
543   -
544   - # Get fitted Crab, photon and energy fluxes
545   - cflux = result_all["Test"]['Prefactor'].value()
546   - cflux_err = result_all["Test"]['Prefactor'].error()
547   - pflux = result_all["Test"].spectral().flux(emin, emax)
548   - eflux = result_all["Test"].spectral().eflux(emin, emax)
549   -
550   - # Compute differential sensitivity
551   - diffSens = pflux / erg_width * erg_mean*erg_mean
552   -
553   - # Compute flux correction factor based on average TS
554   - correct = 1.0
555   - if self.m_use_ts:
556   - if ts > 0:
557   - correct = math.sqrt(self.m_ts_thres/ts)
558   - else:
559   - if cflux > 0 and cflux_err > 0:
560   - correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err)
561   -
562   - # Compute extrapolated fluxes
563   - flux = correct * cflux
564   - pflux = correct * pflux
565   - eflux = correct * eflux
566   - diffSens = correct * diffSens
567   - flux_value.append(flux)
568   - pflux_value.append(pflux)
569   - eflux_value.append(eflux)
570   - diffSens_value.append(diffSens)
571   -
572   - # Write background and test source fit results
573   - if self.logExplicit():
574   - self.log.parformat("Photon flux")
575   - self.log(pflux)
576   - self.log(" ph/cm2/s\n")
577   - self.log.parformat("Energy flux")
578   - self.log(eflux)
579   - self.log(" erg/cm2/s\n")
580   - self.log.parformat("Crab flux")
581   - self.log(cflux*1000.0)
582   - self.log(" mCrab\n")
583   - self.log.parformat("Differential sensitivity")
584   - self.log(diffSens)
585   - self.log(" erg/cm2/s\n")
586   - for model in result_all:
587   - self.log.parformat("Model")
588   - self.log(model.name())
589   - self.log("\n")
590   - for par in model:
591   - self.log(str(par)+"\n")
592   - elif self.logTerse():
593   - self.log.parformat("Iteration "+str(iter))
594   - self.log("TS=")
595   - self.log(ts)
596   - self.log(" ")
597   - self.log("corr=")
598   - self.log(correct)
599   - self.log(" ")
600   - self.log(pflux)
601   - self.log(" ph/cm2/s = ")
602   - self.log(eflux)
603   - self.log(" erg/cm2/s = ")
604   - self.log(cflux*1000.0)
605   - self.log(" mCrab = ")
606   - self.log(diffSens)
607   - self.log(" erg/cm2/s\n")
608   -
609   - # Compute sliding average of extrapolated fitted prefactor,
610   - # photon and energy flux. This damps out fluctuations and
611   - # improves convergence
612   - flux = 0.0
613   - pflux = 0.0
614   - eflux = 0.0
615   - diffSens = 0.0
616   - num = 0.0
617   - for k in range(self.m_num_avg):
618   - inx = len(flux_value) - k - 1
619   - if inx >= 0:
620   - flux += flux_value[inx]
621   - pflux += pflux_value[inx]
622   - eflux += eflux_value[inx]
623   - diffSens += diffSens_value[inx]
624   - num += 1.0
625   - flux /= num
626   - pflux /= num
627   - eflux /= num
628   - diffSens /= num
629   -
630   - # Compare average flux to last average
631   - if iter > self.m_num_avg:
632   - if test_flux > 0:
633   - ratio = flux/test_flux
634   -
635   - # We have 2 convergence criteria:
636   - # 1. The average flux does not change
637   - # 2. The flux correction factor is small
638   - if ratio >= 0.99 and ratio <= 1.01 and \
639   - correct >= 0.9 and correct <= 1.1:
640   - if self.logTerse():
641   - self.log(" Converged ("+str(ratio)+")\n")
642   - break
643   - else:
644   - if self.logTerse():
645   - self.log(" Flux is zero.\n")
646   - break
647   -
648   - # Use average for next iteration
649   - test_flux = flux
650   -
651   - # Exit loop if number of trials exhausted
652   - if (iter >= self.m_max_iter):
653   - if self.logTerse():
654   - self.log(" Test ended after "+str(self.m_max_iter)+" iterations.\n")
655   - break
656   -
657   - # Write fit results
658   - if self.logTerse():
659   - self.log.header3("Fit results")
660   - self.log.parformat("Test statistics")
661   - self.log(ts)
662   - self.log("\n")
663   - self.log.parformat("Photon flux")
664   - self.log(pflux)
665   - self.log(" ph/cm2/s\n")
666   - self.log.parformat("Energy flux")
667   - self.log(eflux)
668   - self.log(" erg/cm2/s\n")
669   - self.log.parformat("Crab flux")
670   - self.log(cflux*1000.0)
671   - self.log(" mCrab\n")
672   - self.log.parformat("Differential sensitivity")
673   - self.log(diffSens)
674   - self.log(" erg/cm2/s\n")
675   - self.log.parformat("Number of simulated events")
676   - self.log(nevents)
677   - self.log("\n")
678   - self.log.header3("Background and test source model fitting")
679   - self.log.parformat("log likelihood")
680   - self.log(LogL_all)
681   - self.log("\n")
682   - self.log.parformat("Number of predicted events")
683   - self.log(npred_all)
684   - self.log("\n")
685   - for model in result_all:
686   - self.log.parformat("Model")
687   - self.log(model.name())
688   - self.log("\n")
689   - for par in model:
690   - self.log(str(par)+"\n")
691   - self.log.header3("Background model fit")
692   - self.log.parformat("log likelihood")
693   - self.log(LogL_bgm)
694   - self.log("\n")
695   - self.log.parformat("Number of predicted events")
696   - self.log(npred_bgm)
697   - self.log("\n")
698   - for model in result_bgm:
699   - self.log.parformat("Model")
700   - self.log(model.name())
701   - self.log("\n")
702   - for par in model:
703   - self.log(str(par)+"\n")
704   -
705   - # Store result
706   - result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
707   - 'crab': flux, 'pflux': pflux, 'eflux': eflux, \
708   - 'diffSens': diffSens}
709   -
710   - # Return result
711   - return result
  387 +
  388 + # Set source name
  389 + source.name("Test")
  390 +
  391 + # Return source
  392 + return source
  393 +
  394 + def get_sensitivity(self, obs, bkg_model, full_model):
  395 + """
  396 + Determine sensitivity for a given observations.
  397 +
  398 + Parameters:
  399 + obs - Observation container
  400 + bkg_model - Background model
  401 + full_model - Source model
  402 + """
  403 + # Set TeV->erg conversion factor
  404 + tev2erg = 1.6021764
  405 +
  406 + # Determine energy boundaries from first observation in
  407 + # the container
  408 + for run in obs:
  409 + emin = run.events().ebounds().emin()
  410 + emax = run.events().ebounds().emax()
  411 + loge = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
  412 + erg_mean = math.pow(10.0, loge) * tev2erg
  413 + erg_width = (emax.TeV()-emin.TeV()) * tev2erg
  414 + break
  415 +
  416 + # Write header
  417 + if self.logTerse():
  418 + self.log("\n")
  419 + self.log.header2("Energies: "+str(emin)+" - "+str(emax))
  420 +
  421 + # Initialise loop
  422 + flux_value = []
  423 + pflux_value = []
  424 + eflux_value = []
  425 + diffSens_value = []
  426 + iter = 0
  427 + test_flux = 0.1 # This is the initial test flux in Crab units
  428 +
  429 + # Loop until we break
  430 + while True:
  431 +
  432 + # Update iteration counter
  433 + iter += 1
  434 +
  435 + # Write header
  436 + if self.logExplicit():
  437 + self.log.header2("Iteration "+str(iter))
  438 +
  439 + # Set test flux
  440 + full_model["Test"]['Prefactor'].value(test_flux)
  441 + obs.models(full_model)
  442 +
  443 + # Simulate events
  444 + sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \
  445 + binsz=self.m_binsz, npix=self.m_npix, \
  446 + log=self.m_log, debug=self.m_debug)
  447 +
  448 + # Determine number of events in simulation
  449 + nevents = 0.0
  450 + for run in sim:
  451 + nevents += run.events().number()
  452 +
  453 + # Write simulation results
  454 + if self.logExplicit():
  455 + self.log.header3("Simulation")
  456 + self.log.parformat("Number of simulated events")
  457 + self.log(nevents)
  458 + self.log("\n")
  459 +
  460 + # Fit background only
  461 + sim.models(bkg_model)
  462 + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  463 + result_bgm = like.obs().models().copy()
  464 + LogL_bgm = like.opt().value()
  465 + npred_bgm = like.obs().npred()
  466 +
  467 + # Assess quality based on a comparison between Npred and Nevents
  468 + quality_bgm = npred_bgm-nevents
  469 +
  470 + # Write background fit results
  471 + if self.logExplicit():
  472 + self.log.header3("Background model fit")
  473 + self.log.parformat("log likelihood")
  474 + self.log(LogL_bgm)
  475 + self.log("\n")
  476 + self.log.parformat("Number of predicted events")
  477 + self.log(npred_bgm)
  478 + self.log("\n")
  479 + self.log.parformat("Fit quality")
  480 + self.log(quality_bgm)
  481 + self.log("\n")
  482 +
  483 + # Start over if the fit quality was bad
  484 + if abs(quality_bgm) > 3.0:
  485 + if self.logExplicit():
  486 + self.log("Fit quality outside required range. Start over.\n")
  487 + continue
  488 +
  489 + # Write model fit results
  490 + if self.logExplicit():
  491 + for model in result_bgm:
  492 + self.log.parformat("Model")
  493 + self.log(model.name())
  494 + self.log("\n")
  495 + for par in model:
  496 + self.log(str(par)+"\n")
  497 +
  498 + # Fit background and test source
  499 + sim.models(full_model)
  500 + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  501 + result_all = like.obs().models().copy()
  502 + LogL_all = like.opt().value()
  503 + npred_all = like.obs().npred()
  504 + ts = 2.0*(LogL_bgm-LogL_all)
  505 +
  506 + # Assess quality based on a comparison between Npred and Nevents
  507 + quality_all = npred_all-nevents
  508 +
  509 + # Write background and test source fit results
  510 + if self.logExplicit():
  511 + self.log.header3("Background and test source model fit")
  512 + self.log.parformat("Test statistics")
  513 + self.log(ts)
  514 + self.log("\n")
  515 + self.log.parformat("log likelihood")
  516 + self.log(LogL_all)
  517 + self.log("\n")
  518 + self.log.parformat("Number of predicted events")
  519 + self.log(npred_all)
  520 + self.log("\n")
  521 + self.log.parformat("Fit quality")
  522 + self.log(quality_all)
  523 + self.log("\n")
  524 + #
  525 + for model in result_all:
  526 + self.log.parformat("Model")
  527 + self.log(model.name())
  528 + self.log("\n")
  529 + for par in model:
  530 + self.log(str(par)+"\n")
  531 +
  532 + # Start over if the fit quality was bad
  533 + if abs(quality_all) > 3.0:
  534 + if self.logExplicit():
  535 + self.log("Fit quality outside required range. Start over.\n")
  536 + continue
  537 +
  538 + # Start over if TS was non-positive
  539 + if ts <= 0.0:
  540 + if self.logExplicit():
  541 + self.log("Non positive TS. Start over.\n")
  542 + continue
  543 +
  544 + # Get fitted Crab, photon and energy fluxes
  545 + cflux = result_all["Test"]['Prefactor'].value()
  546 + cflux_err = result_all["Test"]['Prefactor'].error()
  547 + pflux = result_all["Test"].spectral().flux(emin, emax)
  548 + eflux = result_all["Test"].spectral().eflux(emin, emax)
  549 +
  550 + # Compute differential sensitivity
  551 + diffSens = pflux / erg_width * erg_mean*erg_mean
  552 +
  553 + # Compute flux correction factor based on average TS
  554 + correct = 1.0
  555 + if self.m_use_ts:
  556 + if ts > 0:
  557 + correct = math.sqrt(self.m_ts_thres/ts)
  558 + else:
  559 + if cflux > 0 and cflux_err > 0:
  560 + correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err)
  561 +
  562 + # Compute extrapolated fluxes
  563 + flux = correct * cflux
  564 + pflux = correct * pflux
  565 + eflux = correct * eflux
  566 + diffSens = correct * diffSens
  567 + flux_value.append(flux)
  568 + pflux_value.append(pflux)
  569 + eflux_value.append(eflux)
  570 + diffSens_value.append(diffSens)
  571 +
  572 + # Write background and test source fit results
  573 + if self.logExplicit():
  574 + self.log.parformat("Photon flux")
  575 + self.log(pflux)
  576 + self.log(" ph/cm2/s\n")
  577 + self.log.parformat("Energy flux")
  578 + self.log(eflux)
  579 + self.log(" erg/cm2/s\n")
  580 + self.log.parformat("Crab flux")
  581 + self.log(cflux*1000.0)
  582 + self.log(" mCrab\n")
  583 + self.log.parformat("Differential sensitivity")
  584 + self.log(diffSens)
  585 + self.log(" erg/cm2/s\n")
  586 + for model in result_all:
  587 + self.log.parformat("Model")
  588 + self.log(model.name())
  589 + self.log("\n")
  590 + for par in model:
  591 + self.log(str(par)+"\n")
  592 + elif self.logTerse():
  593 + self.log.parformat("Iteration "+str(iter))
  594 + self.log("TS=")
  595 + self.log(ts)
  596 + self.log(" ")
  597 + self.log("corr=")
  598 + self.log(correct)
  599 + self.log(" ")
  600 + self.log(pflux)
  601 + self.log(" ph/cm2/s = ")
  602 + self.log(eflux)
  603 + self.log(" erg/cm2/s = ")
  604 + self.log(cflux*1000.0)
  605 + self.log(" mCrab = ")
  606 + self.log(diffSens)
  607 + self.log(" erg/cm2/s\n")
  608 +
  609 + # Compute sliding average of extrapolated fitted prefactor,
  610 + # photon and energy flux. This damps out fluctuations and
  611 + # improves convergence
  612 + flux = 0.0
  613 + pflux = 0.0
  614 + eflux = 0.0
  615 + diffSens = 0.0
  616 + num = 0.0
  617 + for k in range(self.m_num_avg):
  618 + inx = len(flux_value) - k - 1
  619 + if inx >= 0:
  620 + flux += flux_value[inx]
  621 + pflux += pflux_value[inx]
  622 + eflux += eflux_value[inx]
  623 + diffSens += diffSens_value[inx]
  624 + num += 1.0
  625 + flux /= num
  626 + pflux /= num
  627 + eflux /= num
  628 + diffSens /= num
  629 +
  630 + # Compare average flux to last average
  631 + if iter > self.m_num_avg:
  632 + if test_flux > 0:
  633 + ratio = flux/test_flux
  634 +
  635 + # We have 2 convergence criteria:
  636 + # 1. The average flux does not change
  637 + # 2. The flux correction factor is small
  638 + if ratio >= 0.99 and ratio <= 1.01 and \
  639 + correct >= 0.9 and correct <= 1.1:
  640 + if self.logTerse():
  641 + self.log(" Converged ("+str(ratio)+")\n")
  642 + break
  643 + else:
  644 + if self.logTerse():
  645 + self.log(" Flux is zero.\n")
  646 + break
  647 +
  648 + # Use average for next iteration
  649 + test_flux = flux
  650 +
  651 + # Exit loop if number of trials exhausted
  652 + if (iter >= self.m_max_iter):
  653 + if self.logTerse():
  654 + self.log(" Test ended after "+str(self.m_max_iter)+" iterations.\n")
  655 + break
  656 +
  657 + # Write fit results
  658 + if self.logTerse():
  659 + self.log.header3("Fit results")
  660 + self.log.parformat("Test statistics")
  661 + self.log(ts)
  662 + self.log("\n")
  663 + self.log.parformat("Photon flux")
  664 + self.log(pflux)
  665 + self.log(" ph/cm2/s\n")
  666 + self.log.parformat("Energy flux")
  667 + self.log(eflux)
  668 + self.log(" erg/cm2/s\n")
  669 + self.log.parformat("Crab flux")
  670 + self.log(cflux*1000.0)
  671 + self.log(" mCrab\n")
  672 + self.log.parformat("Differential sensitivity")
  673 + self.log(diffSens)
  674 + self.log(" erg/cm2/s\n")
  675 + self.log.parformat("Number of simulated events")
  676 + self.log(nevents)
  677 + self.log("\n")
  678 + self.log.header3("Background and test source model fitting")
  679 + self.log.parformat("log likelihood")
  680 + self.log(LogL_all)
  681 + self.log("\n")
  682 + self.log.parformat("Number of predicted events")
  683 + self.log(npred_all)
  684 + self.log("\n")
  685 + for model in result_all:
  686 + self.log.parformat("Model")
  687 + self.log(model.name())
  688 + self.log("\n")
  689 + for par in model:
  690 + self.log(str(par)+"\n")
  691 + self.log.header3("Background model fit")
  692 + self.log.parformat("log likelihood")
  693 + self.log(LogL_bgm)
  694 + self.log("\n")
  695 + self.log.parformat("Number of predicted events")
  696 + self.log(npred_bgm)
  697 + self.log("\n")
  698 + for model in result_bgm:
  699 + self.log.parformat("Model")
  700 + self.log(model.name())
  701 + self.log("\n")
  702 + for par in model:
  703 + self.log(str(par)+"\n")
  704 +
  705 + # Store result
  706 + result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
  707 + 'crab': flux, 'pflux': pflux, 'eflux': eflux, \
  708 + 'diffSens': diffSens}
  709 +
  710 + # Return result
  711 + return result
712 712  
713 713  
714 714 # ======================== #
715 715 # Main routine entry point #
716 716 # ======================== #
717 717 if __name__ == '__main__':
718   - """
719   - CTA ctlike-based sensitivity estimator.
720   - """
721   - # Create instance of application
722   - app = cssens(sys.argv)
723   -
724   - # Run application
725   - app.execute()
726   -
727 718 \ No newline at end of file
  719 + """
  720 + CTA ctlike-based sensitivity estimator.
  721 + """
  722 + # Create instance of application
  723 + app = cssens(sys.argv)
  724 +
  725 + # Run application
  726 + app.execute()
  727 +
728 728 \ No newline at end of file
... ...
scripts/cstsdist.py
  View file @d9a9c73
... ... @@ -3,7 +3,7 @@
3 3 # This script generates the TS distribution for a particular model based
4 4 # on Monte-Carlo simulations.
5 5 #
6   -# Copyright (C) 2011-2013 Jurgen Knodlseder
  6 +# Copyright (C) 2011-2014 Juergen Knoedlseder
7 7 #
8 8 # This program is free software: you can redistribute it and/or modify
9 9 # it under the terms of the GNU General Public License as published by
... ... @@ -31,486 +31,486 @@ import math
31 31 # cstsdist class #
32 32 # ============== #
33 33 class cstsdist(GApplication):
34   - """
35   - This class implements the TS distribution generation script. It derives
36   - from the GammaLib::GApplication class which provides support for parameter
37   - files, command line arguments, and logging. In that way the Python
38   - script behaves just as a regular ctool.
39   - """
40   - def __init__(self, *argv):
41   - """
42   - Constructor.
43   - """
44   - # Set name
45   - self.name = "cstsdist"
46   - self.version = "0.1.0"
47   -
48   - # Initialise some members
49   - self.obs = None
50   - self.bkg_model = None
51   - self.full_model = None
52   -
53   - # Make sure that parfile exists
54   - file = self.parfile()
  34 + """
  35 + This class implements the TS distribution generation script. It derives
  36 + from the GammaLib::GApplication class which provides support for parameter
  37 + files, command line arguments, and logging. In that way the Python
  38 + script behaves just as a regular ctool.
  39 + """
  40 + def __init__(self, *argv):
  41 + """
  42 + Constructor.
  43 + """
  44 + # Set name
  45 + self.name = "cstsdist"
  46 + self.version = "0.1.0"
  47 +
  48 + # Initialise some members
  49 + self.obs = None
  50 + self.bkg_model = None
  51 + self.full_model = None
  52 +
  53 + # Make sure that parfile exists
  54 + file = self.parfile()
55 55  
56   - # Initialise application
57   - if len(argv) == 0:
58   - GApplication.__init__(self, self.name, self.version)
59   - elif len(argv) ==1:
60   - GApplication.__init__(self, self.name, self.version, *argv)
61   - else:
62   - raise TypeError("Invalid number of arguments given.")
  56 + # Initialise application
  57 + if len(argv) == 0:
  58 + GApplication.__init__(self, self.name, self.version)
  59 + elif len(argv) ==1:
  60 + GApplication.__init__(self, self.name, self.version, *argv)
  61 + else:
  62 + raise TypeError("Invalid number of arguments given.")
63 63  
64   - # Set logger properties
65   - self.log_header()
66   - self.log.date(True)
  64 + # Set logger properties
  65 + self.log_header()
  66 + self.log.date(True)
67 67  
68   - # Return
69   - return
70   -
71   - def __del__(self):
72   - """
73   - Destructor.
74   - """
75   - # Write separator into logger
76   - if self.logTerse():
77   - self.log("\n")
78   -
79   - # Return
80   - return
  68 + # Return
  69 + return
  70 +
  71 + def __del__(self):
  72 + """
  73 + Destructor.
  74 + """
  75 + # Write separator into logger
  76 + if self.logTerse():
  77 + self.log("\n")
  78 +
  79 + # Return
  80 + return
81 81  
82   - def parfile(self):
83   - """
84   - Check if parfile exists. If parfile does not exist then create a
85   - default parfile. This kluge avoids shipping the cscript with a parfile.
86   - """
87   - # Set parfile name
88   - parfile = self.name+".par"
89   -
90   - try:
91   - pars = GApplicationPars(parfile)
92   - except:
93   - # Signal if parfile was not found
94   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
95   -
96   - # Create default parfile
97   - pars = GApplicationPars()
98   - pars.append(GApplicationPar("outfile","f","h","ts.dat","","","Output file name"))
99   - pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))
100   - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
101   - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
102   - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
103   - pars.append(GApplicationPar("index","r","h","-2.48","","","Spectral index"))
104   - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
105   - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
106   - pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
107   - pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
108   - pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
109   - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
110   - pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
111   - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
112   - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
113   - pars.append_standard()
114   - pars.save(parfile)
115   -
116   - # Return
117   - return
118   -
119   - def get_parameters(self):
120   - """
121   - Get parameters from parfile and setup the observation.
122   - """
123   - # Get parameters
124   - self.m_outfile = self["outfile"].filename()
125   - self.m_ntrials = self["ntrials"].integer()
126   - self.m_caldb = self["caldb"].string()
127   - self.m_irf = self["irf"].string()
128   - self.m_type = self["type"].string()
129   - self.m_index = self["index"].real()
130   - self.m_offset = self["offset"].real()
131   - self.m_bkg = self["bkg"].string()
132   - self.m_emin = self["emin"].real()
133   - self.m_emax = self["emax"].real()
134   - self.m_enumbins = self["enumbins"].integer()
135   - self.m_duration = self["duration"].real()
136   - self.m_rad = self["rad"].real()
137   - self.m_npix = self["npix"].integer()
138   - self.m_binsz = self["binsz"].real()
  82 + def parfile(self):
  83 + """
  84 + Check if parfile exists. If parfile does not exist then create a
  85 + default parfile. This kluge avoids shipping the cscript with a parfile.
  86 + """
  87 + # Set parfile name
  88 + parfile = self.name+".par"
  89 +
  90 + try:
  91 + pars = GApplicationPars(parfile)
  92 + except:
  93 + # Signal if parfile was not found
  94 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  95 +
  96 + # Create default parfile
  97 + pars = GApplicationPars()
  98 + pars.append(GApplicationPar("outfile","f","h","ts.dat","","","Output file name"))
  99 + pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))
  100 + pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  101 + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  102 + pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
  103 + pars.append(GApplicationPar("index","r","h","-2.48","","","Spectral index"))
  104 + pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
  105 + pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
  106 + pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
  107 + pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  108 + pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
  109 + pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  110 + pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  111 + pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  112 + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  113 + pars.append_standard()
  114 + pars.save(parfile)
  115 +
  116 + # Return
  117 + return
  118 +
  119 + def get_parameters(self):
  120 + """
  121 + Get parameters from parfile and setup the observation.
  122 + """
  123 + # Get parameters
  124 + self.m_outfile = self["outfile"].filename()
  125 + self.m_ntrials = self["ntrials"].integer()
  126 + self.m_caldb = self["caldb"].string()
  127 + self.m_irf = self["irf"].string()
  128 + self.m_type = self["type"].string()
  129 + self.m_index = self["index"].real()
  130 + self.m_offset = self["offset"].real()
  131 + self.m_bkg = self["bkg"].string()
  132 + self.m_emin = self["emin"].real()
  133 + self.m_emax = self["emax"].real()
  134 + self.m_enumbins = self["enumbins"].integer()
  135 + self.m_duration = self["duration"].real()
  136 + self.m_rad = self["rad"].real()
  137 + self.m_npix = self["npix"].integer()
  138 + self.m_binsz = self["binsz"].real()
139 139  
140   - # Set some fixed parameters
141   - self.m_log = False # Logging in client tools
142   - self.m_debug = False # Debugging in client tools
143   -
144   - # Setup observation
145   - self.obs = self.set_obs(emin=self.m_emin, emax=self.m_emax)
146   -
147   - # Initialise models. Note that we centre the point source at the Galactic
148   - # center as our observation is also centred at the Galactic centre, so
149   - # we're onaxis.
150   - self.bkg_model = GModels()
151   - self.full_model = GModels()
152   - self.bkg_model.append(self.set_bkg_model())
153   - self.full_model.append(self.set_bkg_model())
154   - self.full_model.append(self.set_src_model(0.0, self.m_offset, \
155   - flux=0.010, \
156   - type=self.m_type, \
157   - index=self.m_index))
  140 + # Set some fixed parameters
  141 + self.m_log = False # Logging in client tools
  142 + self.m_debug = False # Debugging in client tools
  143 +
  144 + # Setup observation
  145 + self.obs = self.set_obs(emin=self.m_emin, emax=self.m_emax)
  146 +
  147 + # Initialise models. Note that we centre the point source at the Galactic
  148 + # center as our observation is also centred at the Galactic centre, so
  149 + # we're onaxis.
  150 + self.bkg_model = GModels()
  151 + self.full_model = GModels()
  152 + self.bkg_model.append(self.set_bkg_model())
  153 + self.full_model.append(self.set_bkg_model())
  154 + self.full_model.append(self.set_src_model(0.0, self.m_offset, \
  155 + flux=0.010, \
  156 + type=self.m_type, \
  157 + index=self.m_index))
158 158  
159   - # Attach background model to observation container
160   - self.obs.models(self.bkg_model)
161   -
162   - # Return
163   - return
164   -
165   - def models(self, models):
166   - """
167   - Set model.
168   - """
169   - # Copy models
170   - self.model = models.copy()
171   -
172   - # Return
173   - return
174   -
175   - def execute(self):
176   - """
177   - Execute the script.
178   - """
179   - # Run the script
180   - self.run()
181   -
182   - # Return
183   - return
  159 + # Attach background model to observation container
  160 + self.obs.models(self.bkg_model)
  161 +
  162 + # Return
  163 + return
  164 +
  165 + def models(self, models):
  166 + """
  167 + Set model.
  168 + """
  169 + # Copy models
  170 + self.model = models.copy()
  171 +
  172 + # Return
  173 + return
  174 +
  175 + def execute(self):
  176 + """
  177 + Execute the script.
  178 + """
  179 + # Run the script
  180 + self.run()
  181 +
  182 + # Return
  183 + return
184 184  
185   - def run(self):
186   - """
187   - Run the script.
188   - """
189   - # Switch screen logging on in debug mode
190   - if self.logDebug():
191   - self.log.cout(True)
  185 + def run(self):
  186 + """
  187 + Run the script.
  188 + """
  189 + # Switch screen logging on in debug mode
  190 + if self.logDebug():
  191 + self.log.cout(True)
192 192  
193   - # Get parameters
194   - self.get_parameters()
195   -
196   - # Write input parameters into logger
197   - if self.logTerse():
198   - self.log_parameters()
199   - self.log("\n")
200   -
201   - # Write observation into logger
202   - if self.logTerse():
203   - self.log("\n")
204   - self.log.header1("Observation")
205   - self.log(str(self.obs))
206   - self.log("\n")
  193 + # Get parameters
  194 + self.get_parameters()
  195 +
  196 + # Write input parameters into logger
  197 + if self.logTerse():
  198 + self.log_parameters()
  199 + self.log("\n")
  200 +
  201 + # Write observation into logger
  202 + if self.logTerse():
  203 + self.log("\n")
  204 + self.log.header1("Observation")
  205 + self.log(str(self.obs))
  206 + self.log("\n")
207 207  
208   - # Write models into logger
209   - if self.logTerse():
210   - self.log("\n")
211   - self.log.header1("Test model")
212   - self.log(str(self.full_model))
213   - self.log("\n")
  208 + # Write models into logger
  209 + if self.logTerse():
  210 + self.log("\n")
  211 + self.log.header1("Test model")
  212 + self.log(str(self.full_model))
  213 + self.log("\n")
214 214  
215   - # Write header
216   - if self.logTerse():
217   - self.log("\n")
218   - self.log.header1("Generate TS distribution")
  215 + # Write header
  216 + if self.logTerse():
  217 + self.log("\n")
  218 + self.log.header1("Generate TS distribution")
219 219  
220   - # Loop over trials
221   - for seed in range(self.m_ntrials):
222   -
223   - # Make a trial
224   - result = self.trial(seed)
225   -
226   - # Write out result immediately
227   - if seed == 0:
228   - file = open(self.m_outfile, 'w')
229   - writer = csv.DictWriter(file, result['colnames'])
230   - writer.writerow(dict((_,_) for _ in result['colnames']))
231   - else:
232   - file = open(self.m_outfile, 'a')
233   - writer = csv.DictWriter(file, result['colnames'])
234   - writer.writerow(result['values'])
235   - file.close()
236   -
237   - # Return
238   - return
239   -
240   - def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
241   - """
242   - Returns an observation container with a single CTA observation.
243   -
244   - Keywords:
245   - lpnt - Galactic longitude of pointing [deg] (default: 0.0)
246   - bpnt - Galactic latitude of pointing [deg] (default: 0.0)
247   - emin - Minimum energy [TeV] (default: 0.1)
248   - emax - Maximum energy [TeV] (default: 100.0)
249   - """
250   - # Allocate observation container
251   - obs = GObservations()
252   -
253   - # Set single pointing
254   - pntdir = GSkyDir()
255   - pntdir.lb_deg(lpnt, bpnt)
256   -
257   - # Create CTA observation
258   - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
259   - duration=self.m_duration, \
  220 + # Loop over trials
  221 + for seed in range(self.m_ntrials):
  222 +
  223 + # Make a trial
  224 + result = self.trial(seed)
  225 +
  226 + # Write out result immediately
  227 + if seed == 0:
  228 + file = open(self.m_outfile, 'w')
  229 + writer = csv.DictWriter(file, result['colnames'])
  230 + writer.writerow(dict((_,_) for _ in result['colnames']))
  231 + else:
  232 + file = open(self.m_outfile, 'a')
  233 + writer = csv.DictWriter(file, result['colnames'])
  234 + writer.writerow(result['values'])
  235 + file.close()
  236 +
  237 + # Return
  238 + return
  239 +
  240 + def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
  241 + """
  242 + Returns an observation container with a single CTA observation.
  243 +
  244 + Keywords:
  245 + lpnt - Galactic longitude of pointing [deg] (default: 0.0)
  246 + bpnt - Galactic latitude of pointing [deg] (default: 0.0)
  247 + emin - Minimum energy [TeV] (default: 0.1)
  248 + emax - Maximum energy [TeV] (default: 100.0)
  249 + """
  250 + # Allocate observation container
  251 + obs = GObservations()
  252 +
  253 + # Set single pointing
  254 + pntdir = GSkyDir()
  255 + pntdir.lb_deg(lpnt, bpnt)
  256 +
  257 + # Create CTA observation
  258 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  259 + duration=self.m_duration, \
260 260 emin=emin, emax=emax, rad=self.m_rad)
261   -
262   - # Append observation to container
263   - obs.append(run)
264   -
265   - # Return observation container
266   - return obs
267   -
268   - def set_bkg_model(self, fitsigma=False):
269   - """
270   - Setup CTA background model.
271   - """
272   - # Define radial component
273   - radial = GCTAModelRadialGauss(3.0)
274   - if fitsigma:
275   - radial["Sigma"].free()
276   - else:
277   - radial["Sigma"].fix()
278   -
279   - # Define spectral component
280   - spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
281   -
282   - # Create background model
283   - model = GCTAModelRadialAcceptance(radial, spectrum)
284   - model.name("Background")
285   - model.instruments("CTA")
286   -
287   - # Return background model
288   - return model
289   -
290   - def set_src_model(self, l, b, flux=1.0, index=-2.48, \
  261 +
  262 + # Append observation to container
  263 + obs.append(run)
  264 +
  265 + # Return observation container
  266 + return obs
  267 +
  268 + def set_bkg_model(self, fitsigma=False):
  269 + """
  270 + Setup CTA background model.
  271 + """
  272 + # Define radial component
  273 + radial = GCTAModelRadialGauss(3.0)
  274 + if fitsigma:
  275 + radial["Sigma"].free()
  276 + else:
  277 + radial["Sigma"].fix()
  278 +
  279 + # Define spectral component
  280 + spectrum = GModelSpectralFunc(self.m_bkg, 1.0)
  281 +
  282 + # Create background model
  283 + model = GCTAModelRadialAcceptance(radial, spectrum)
  284 + model.name("Background")
  285 + model.instruments("CTA")
  286 +
  287 + # Return background model
  288 + return model
  289 +
  290 + def set_src_model(self, l, b, flux=1.0, index=-2.48, \
291 291 type="point", sigma=1.0, radius=1.0, width=0.1, \
292 292 fitpos=False, fitidx=False):
293   - """
294   - Returns a single source with Crab-like spectrum. The source flux
295   - can be scaled in Crab units. The Crab spectrum is based on MAGIC
296   - observations (Albert et al. 2008, ApJ, 674, 1037).
  293 + """
  294 + Returns a single source with Crab-like spectrum. The source flux
  295 + can be scaled in Crab units. The Crab spectrum is based on MAGIC
  296 + observations (Albert et al. 2008, ApJ, 674, 1037).
297 297  
298   - Parameters:
299   - l - Galactic longitude of source location [deg]
300   - b - Galactic latitude of source location [deg]
301   - Keywords:
302   - flux - Source flux [Crabs]
303   - index - Spectral index
304   - type - Source type ("point", "gauss", "disk", "shell")
305   - sigma - Gaussian sigma (for type="gauss")
306   - radius - Disk or shell inner radius [deg] (for type="disk" and type="shell")
307   - width - Shell width [deg] (for type="shell")
308   - fitpos - Fit position? (default: True)
309   - fitidx - Fit index? (default: True)
310   - """
311   - # Set source location
312   - location = GSkyDir()
313   - location.lb_deg(l, b)
314   -
315   - # Set source spectrum
316   - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV"))
317   - if fitidx:
318   - spectrum["Index"].free()
319   - else:
320   - spectrum["Index"].fix()
  298 + Parameters:
  299 + l - Galactic longitude of source location [deg]
  300 + b - Galactic latitude of source location [deg]
  301 + Keywords:
  302 + flux - Source flux [Crabs]
  303 + index - Spectral index
  304 + type - Source type ("point", "gauss", "disk", "shell")
  305 + sigma - Gaussian sigma (for type="gauss")
  306 + radius - Disk or shell inner radius [deg] (for type="disk" and type="shell")
  307 + width - Shell width [deg] (for type="shell")
  308 + fitpos - Fit position? (default: True)
  309 + fitidx - Fit index? (default: True)
  310 + """
  311 + # Set source location
  312 + location = GSkyDir()
  313 + location.lb_deg(l, b)
  314 +
  315 + # Set source spectrum
  316 + spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV"))
  317 + if fitidx:
  318 + spectrum["Index"].free()
  319 + else:
  320 + spectrum["Index"].fix()
321 321  
322   - # Set source
323   - if type == "point":
324   - spatial = GModelSpatialPointSource(location)
325   - if fitpos:
326   - spatial[0].free()
327   - spatial[1].free()
328   - elif type == "gauss":
329   - spatial = GModelSpatialRadialGauss(location, sigma)
330   - if fitpos:
331   - spatial[0].free()
332   - spatial[1].free()
333   - elif type == "disk":
334   - spatial = GModelSpatialRadialDisk(location, radius)
335   - if fitpos:
336   - spatial[0].free()
337   - spatial[1].free()
338   - elif type == "shell":
339   - spatial = GModelSpatialRadialShell(location, radius, width)
340   - if fitpos:
341   - spatial[0].free()
342   - spatial[1].free()
343   - else:
344   - self.log("ERROR: Unknown source type '"+type+"'.\n")
345   - return None
346   - source = GModelSky(spatial, spectrum)
347   -
348   - # Set source name
349   - source.name("Test")
350   -
351   - # Return source
352   - return source
353   -
354   - def trial(self, seed):
355   - """
356   - Create the TS for a single trial.
357   -
358   - Parameters:
359   - seed - Random number generator seed
360   - """
361   - # Write header
362   - if self.logExplicit():
363   - self.log.header2("Trial "+str(seed+1))
  322 + # Set source
  323 + if type == "point":
  324 + spatial = GModelSpatialPointSource(location)
  325 + if fitpos:
  326 + spatial[0].free()
  327 + spatial[1].free()
  328 + elif type == "gauss":
  329 + spatial = GModelSpatialRadialGauss(location, sigma)
  330 + if fitpos:
  331 + spatial[0].free()
  332 + spatial[1].free()
  333 + elif type == "disk":
  334 + spatial = GModelSpatialRadialDisk(location, radius)
  335 + if fitpos:
  336 + spatial[0].free()
  337 + spatial[1].free()
  338 + elif type == "shell":
  339 + spatial = GModelSpatialRadialShell(location, radius, width)
  340 + if fitpos:
  341 + spatial[0].free()
  342 + spatial[1].free()
  343 + else:
  344 + self.log("ERROR: Unknown source type '"+type+"'.\n")
  345 + return None
  346 + source = GModelSky(spatial, spectrum)
  347 +
  348 + # Set source name
  349 + source.name("Test")
  350 +
  351 + # Return source
  352 + return source
  353 +
  354 + def trial(self, seed):
  355 + """
  356 + Create the TS for a single trial.
  357 +
  358 + Parameters:
  359 + seed - Random number generator seed
  360 + """
  361 + # Write header
  362 + if self.logExplicit():
  363 + self.log.header2("Trial "+str(seed+1))
364 364  
365   - # Simulate events
366   - sim = obsutils.sim(self.obs, \
367   - nbins=self.m_enumbins, \
368   - seed=seed, \
369   - binsz=self.m_binsz, \
370   - npix=self.m_npix, \
371   - log=self.m_log, debug=self.m_debug)
  365 + # Simulate events
  366 + sim = obsutils.sim(self.obs, \
  367 + nbins=self.m_enumbins, \
  368 + seed=seed, \
  369 + binsz=self.m_binsz, \
  370 + npix=self.m_npix, \
  371 + log=self.m_log, debug=self.m_debug)
372 372  
373   - # Determine number of events in simulation
374   - nevents = 0.0
375   - for run in sim:
376   - nevents += run.events().number()
  373 + # Determine number of events in simulation
  374 + nevents = 0.0
  375 + for run in sim:
  376 + nevents += run.events().number()
377 377  
378   - # Write simulation results
379   - if self.logExplicit():
380   - self.log.header3("Simulation")
381   - self.log.parformat("Number of simulated events")
382   - self.log(nevents)
383   - self.log("\n")
  378 + # Write simulation results
  379 + if self.logExplicit():
  380 + self.log.header3("Simulation")
  381 + self.log.parformat("Number of simulated events")
  382 + self.log(nevents)
  383 + self.log("\n")
384 384  
385   - # Fit background only
386   - sim.models(self.bkg_model)
387   - like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
388   - result_bgm = like_bgm.obs().models()
389   - LogL_bgm = like_bgm.opt().value()
390   - npred_bgm = like_bgm.obs().npred()
  385 + # Fit background only
  386 + sim.models(self.bkg_model)
  387 + like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  388 + result_bgm = like_bgm.obs().models()
  389 + LogL_bgm = like_bgm.opt().value()
  390 + npred_bgm = like_bgm.obs().npred()
391 391  
392   - # Write background fit results
393   - if self.logExplicit():
394   - self.log.header3("Background model fit")
395   - self.log.parformat("log likelihood")
396   - self.log(LogL_bgm)
397   - self.log("\n")
398   - self.log.parformat("Number of predicted events")
399   - self.log(npred_bgm)
400   - self.log("\n")
401   - for model in result_bgm:
402   - self.log.parformat("Model")
403   - self.log(model.name())
404   - self.log("\n")
405   - for par in model:
406   - self.log(str(par)+"\n")
  392 + # Write background fit results
  393 + if self.logExplicit():
  394 + self.log.header3("Background model fit")
  395 + self.log.parformat("log likelihood")
  396 + self.log(LogL_bgm)
  397 + self.log("\n")
  398 + self.log.parformat("Number of predicted events")
  399 + self.log(npred_bgm)
  400 + self.log("\n")
  401 + for model in result_bgm:
  402 + self.log.parformat("Model")
  403 + self.log(model.name())
  404 + self.log("\n")
  405 + for par in model:
  406 + self.log(str(par)+"\n")
407 407  
408   - # Fit background and test source
409   - sim.models(self.full_model)
410   - like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
411   - result_all = like_all.obs().models()
412   - LogL_all = like_all.opt().value()
413   - npred_all = like_all.obs().npred()
414   - ts = 2.0*(LogL_bgm-LogL_all)
415   -
416   - # Write background and test source fit results
417   - if self.logExplicit():
418   - self.log.header3("Background and test source model fit")
419   - self.log.parformat("Test statistics")
420   - self.log(ts)
421   - self.log("\n")
422   - self.log.parformat("log likelihood")
423   - self.log(LogL_all)
424   - self.log("\n")
425   - self.log.parformat("Number of predicted events")
426   - self.log(npred_all)
427   - self.log("\n")
428   - for model in result_all:
429   - self.log.parformat("Model")
430   - self.log(model.name())
431   - self.log("\n")
432   - for par in model:
433   - self.log(str(par)+"\n")
434   -
435   - # Write result
436   - elif self.logTerse():
437   - self.log.parformat("Trial "+str(seed))
438   - self.log("TS=")
439   - self.log(ts)
440   - self.log(" Prefactor=")
441   - self.log(result_all["Test"]["Prefactor"].value())
442   - self.log("+/-")
443   - self.log(result_all["Test"]["Prefactor"].error())
444   - self.log("\n")
445   -
446   - # Initialise results
447   - colnames = []
448   - values = {}
449   -
450   - # Set TS value
451   - colnames.append("TS")
452   - values["TS"] = ts
  408 + # Fit background and test source
  409 + sim.models(self.full_model)
  410 + like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  411 + result_all = like_all.obs().models()
  412 + LogL_all = like_all.opt().value()
  413 + npred_all = like_all.obs().npred()
  414 + ts = 2.0*(LogL_bgm-LogL_all)
  415 +
  416 + # Write background and test source fit results
  417 + if self.logExplicit():
  418 + self.log.header3("Background and test source model fit")
  419 + self.log.parformat("Test statistics")
  420 + self.log(ts)
  421 + self.log("\n")
  422 + self.log.parformat("log likelihood")
  423 + self.log(LogL_all)
  424 + self.log("\n")
  425 + self.log.parformat("Number of predicted events")
  426 + self.log(npred_all)
  427 + self.log("\n")
  428 + for model in result_all:
  429 + self.log.parformat("Model")
  430 + self.log(model.name())
  431 + self.log("\n")
  432 + for par in model:
  433 + self.log(str(par)+"\n")
  434 +
  435 + # Write result
  436 + elif self.logTerse():
  437 + self.log.parformat("Trial "+str(seed))
  438 + self.log("TS=")
  439 + self.log(ts)
  440 + self.log(" Prefactor=")
  441 + self.log(result_all["Test"]["Prefactor"].value())
  442 + self.log("+/-")
  443 + self.log(result_all["Test"]["Prefactor"].error())
  444 + self.log("\n")
  445 +
  446 + # Initialise results
  447 + colnames = []
  448 + values = {}
  449 +
  450 + # Set TS value
  451 + colnames.append("TS")
  452 + values["TS"] = ts
453 453  
454   - # Set logL for background fit
455   - colnames.append("LogL_bgm")
456   - values["LogL_bgm"] = LogL_bgm
  454 + # Set logL for background fit
  455 + colnames.append("LogL_bgm")
  456 + values["LogL_bgm"] = LogL_bgm
457 457  
458   - # Set logL for full fit
459   - colnames.append("LogL_all")
460   - values["LogL_all"] = LogL_all
  458 + # Set logL for full fit
  459 + colnames.append("LogL_all")
  460 + values["LogL_all"] = LogL_all
461 461  
462   - # Set Nevents
463   - colnames.append("Nevents")
464   - values["Nevents"] = nevents
  462 + # Set Nevents
  463 + colnames.append("Nevents")
  464 + values["Nevents"] = nevents
465 465  
466   - # Set Npred for background fit
467   - colnames.append("Npred_bkg")
468   - values["Npred_bkg"] = npred_bgm
  466 + # Set Npred for background fit
  467 + colnames.append("Npred_bkg")
  468 + values["Npred_bkg"] = npred_bgm
469 469  
470   - # Set Npred for full fit
471   - colnames.append("Npred_all")
472   - values["Npred_all"] = npred_all
473   -
474   - # Gather free full fit parameters
475   - for i in range(result_all.size()):
476   - model = result_all[i]
477   - model_name = model.name()
478   - for k in range(model.size()):
479   - par = model[k]
480   - if par.is_free():
481   -
482   - # Set parameter name
483   - name = model_name+"_"+par.name()
484   -
485   - # Append value
486   - colnames.append(name)
487   - values[name] = par.value()
488   -
489   - # Append error
490   - name = "Unc_"+name
491   - colnames.append(name)
492   - values[name] = par.error()
493   -
494   - # Bundle together results
495   - result = {'colnames': colnames, 'values': values}
496   -
497   - # Return
498   - return result
  470 + # Set Npred for full fit
  471 + colnames.append("Npred_all")
  472 + values["Npred_all"] = npred_all
  473 +
  474 + # Gather free full fit parameters
  475 + for i in range(result_all.size()):
  476 + model = result_all[i]
  477 + model_name = model.name()
  478 + for k in range(model.size()):
  479 + par = model[k]
  480 + if par.is_free():
  481 +
  482 + # Set parameter name
  483 + name = model_name+"_"+par.name()
  484 +
  485 + # Append value
  486 + colnames.append(name)
  487 + values[name] = par.value()
  488 +
  489 + # Append error
  490 + name = "Unc_"+name
  491 + colnames.append(name)
  492 + values[name] = par.error()
  493 +
  494 + # Bundle together results
  495 + result = {'colnames': colnames, 'values': values}
  496 +
  497 + # Return
  498 + return result
499 499  
500 500  
501 501 # ======================== #
502 502 # Main routine entry point #
503 503 # ======================== #
504 504 if __name__ == '__main__':
505   - """
506   - Generates TS distribution.
507   - """
508   - # Create instance of application
509   - app = cstsdist(sys.argv)
510   -
511   - # Open logfile
512   - app.logFileOpen()
513   -
514   - # Execute application
515   - app.execute()
516   -
517 505 \ No newline at end of file
  506 + """
  507 + Generates TS distribution.
  508 + """
  509 + # Create instance of application
  510 + app = cstsdist(sys.argv)
  511 +
  512 + # Open logfile
  513 + app.logFileOpen()
  514 +
  515 + # Execute application
  516 + app.execute()
  517 +
518 518 \ No newline at end of file
... ...
scripts/cta_kb_bkgrates.py
  View file @d9a9c73
... ... @@ -4,7 +4,7 @@
4 4 # into a column separated file function table. The table is automatically
5 5 # named "bkg_XXX.txt" where "XXX" is for example "kb_E_50h_v3".
6 6 #
7   -# Copyright (C) 2011 Jurgen Knodlseder
  7 +# Copyright (C) 2011-2014 Juergen Knoedlseder
8 8 #
9 9 # This program is free software: you can redistribute it and/or modify
10 10 # it under the terms of the GNU General Public License as published by
... ... @@ -29,78 +29,78 @@ import sys
29 29 # Generate file function from model #
30 30 # ================================= #
31 31 def make_file_function(irfname, filename):
32   - """
33   - Generate file function from model.
34   -
35   - Parameters:
36   - irfname - IRF filename
37   - filename - File function filename
38   - """
39   - # Open performance file and file function
40   - irf = open(irfname, "r")
41   - file = open(filename, "w")
42   -
43   - # Read data
44   - for row in irf:
45   -
46   - # Skip header and trailer
47   - if row.find("log(E)") != -1:
48   - continue
49   - if row.find("--------") != -1:
50   - break
51   -
52   - # Extract data
53   - fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",")
54   - logE = float(fields[0])
55   - aeff = float(fields[1])
56   - r68 = float(fields[2])
57   - r80 = float(fields[3])
58   - eres = float(fields[4])
59   - bgd = float(fields[5])
60   - sens = float(fields[6])
61   -
62   - # Compute energy (in MeV)
63   - energy = math.pow(10.0, logE)*1.0e6
64   - emin = math.pow(10.0, logE-0.1)*1.0e6
65   - emax = math.pow(10.0, logE+0.1)*1.0e6
66   - ewidth = emax-emin
67   -
68   - # Compute solid angle of r80 region in sr
69   - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
70   -
71   - # Compute background rate per steradian and MeV
72   - bkg_rate = bgd/omega/ewidth
  32 + """
  33 + Generate file function from model.
  34 +
  35 + Parameters:
  36 + irfname - IRF filename
  37 + filename - File function filename
  38 + """
  39 + # Open performance file and file function
  40 + irf = open(irfname, "r")
  41 + file = open(filename, "w")
  42 +
  43 + # Read data
  44 + for row in irf:
  45 +
  46 + # Skip header and trailer
  47 + if row.find("log(E)") != -1:
  48 + continue
  49 + if row.find("--------") != -1:
  50 + break
  51 +
  52 + # Extract data
  53 + fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",")
  54 + logE = float(fields[0])
  55 + aeff = float(fields[1])
  56 + r68 = float(fields[2])
  57 + r80 = float(fields[3])
  58 + eres = float(fields[4])
  59 + bgd = float(fields[5])
  60 + sens = float(fields[6])
  61 +
  62 + # Compute energy (in MeV)
  63 + energy = math.pow(10.0, logE)*1.0e6
  64 + emin = math.pow(10.0, logE-0.1)*1.0e6
  65 + emax = math.pow(10.0, logE+0.1)*1.0e6
  66 + ewidth = emax-emin
  67 +
  68 + # Compute solid angle of r80 region in sr
  69 + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
  70 +
  71 + # Compute background rate per steradian and MeV
  72 + bkg_rate = bgd/omega/ewidth
73 73  
74   - # Write
75   - file.write(str(energy)+" "+str(bkg_rate)+"\n")
76   -
77   - # Close files
78   - irf.close()
79   - file.close()
80   -
81   - # Return
82   - return
  74 + # Write
  75 + file.write(str(energy)+" "+str(bkg_rate)+"\n")
  76 +
  77 + # Close files
  78 + irf.close()
  79 + file.close()
  80 +
  81 + # Return
  82 + return
83 83  
84 84  
85 85 # ================ #
86 86 # Main entry point #
87 87 # ================ #
88 88 if __name__ == '__main__':
89   - """
90   - Main entry point
91   - """
92   - # Get IRF from command line argument
93   - if (len(sys.argv) > 1):
94   - irf = sys.argv[1]
95   - else:
96   - irf = "kb_E_50h_v3"
  89 + """
  90 + Main entry point
  91 + """
  92 + # Get IRF from command line argument
  93 + if (len(sys.argv) > 1):
  94 + irf = sys.argv[1]
  95 + else:
  96 + irf = "kb_E_50h_v3"
97 97  
98   - # Set parameters
99   - caldb = "/usr/local/gamma/share/caldb/cta"
100   -
101   - # Set filenames
102   - irfname = caldb+"/"+irf+".dat"
103   - filename = "bkg_"+irf+".txt"
104   -
105   - # Read
106   - make_file_function(irfname, filename)
  98 + # Set parameters
  99 + caldb = "/usr/local/gamma/share/caldb/cta"
  100 +
  101 + # Set filenames
  102 + irfname = caldb+"/"+irf+".dat"
  103 + filename = "bkg_"+irf+".txt"
  104 +
  105 + # Read
  106 + make_file_function(irfname, filename)
... ...
scripts/cta_make_nodes.py
  View file @d9a9c73
... ... @@ -4,7 +4,7 @@
4 4 # logarithmically spaced energy bins. The intensity scale is set to the
5 5 # HESS Crab intensity (assuming a power law).
6 6 #
7   -# Copyright (C) 2012 Juergen Knoedlseder
  7 +# Copyright (C) 2012-2014 Juergen Knoedlseder
8 8 #
9 9 # This program is free software: you can redistribute it and/or modify
10 10 # it under the terms of the GNU General Public License as published by
... ... @@ -29,91 +29,91 @@ import math
29 29 # Write one node #
30 30 # ============== #
31 31 def write_node(file, energy, scale):
32   - """
33   - Writes one node to XML file.
34   - """
35   - # Convert to strings
36   - str_energy = str(energy)
37   - str_scale = str(scale)
  32 + """
  33 + Writes one node to XML file.
  34 + """
  35 + # Convert to strings
  36 + str_energy = str(energy)
  37 + str_scale = str(scale)
38 38  
39   - # Write start tag
40   - file.write(' <node>\n')
  39 + # Write start tag
  40 + file.write(' <node>\n')
41 41  
42   - # Write energy
43   - file.write(' <parameter scale="1e6" name="Energy"')
44   - file.write(' min="'+str_energy+'" max="'+str_energy+'"')
45   - file.write(' value="'+str_energy+'"')
46   - file.write(' free="0"/>\n')
  42 + # Write energy
  43 + file.write(' <parameter scale="1e6" name="Energy"')
  44 + file.write(' min="'+str_energy+'" max="'+str_energy+'"')
  45 + file.write(' value="'+str_energy+'"')
  46 + file.write(' free="0"/>\n')
47 47  
48   - # Write intensity
49   - file.write(' <parameter scale="'+str_scale+'" name="Intensity"')
50   - file.write(' min="1e-5" max="1e5"')
51   - file.write(' value="1.0"')
52   - file.write(' free="1"/>\n')
  48 + # Write intensity
  49 + file.write(' <parameter scale="'+str_scale+'" name="Intensity"')
  50 + file.write(' min="1e-5" max="1e5"')
  51 + file.write(' value="1.0"')
  52 + file.write(' free="1"/>\n')
53 53  
54   - # Write end tag
55   - file.write(' </node>\n')
  54 + # Write end tag
  55 + file.write(' </node>\n')
56 56  
57   - # Return
58   - return
59   -
  57 + # Return
  58 + return
  59 +
60 60  
61 61 # ============ #
62 62 # Create nodes #
63 63 # ============ #
64 64 def create_nodes(emin, emax, enumbins):
65   - """
66   - Create nodes (energies in TeV).
67   - """
68   - # Open file
69   - file = open("nodes.xml", "w")
  65 + """
  66 + Create nodes (energies in TeV).
  67 + """
  68 + # Open file
  69 + file = open("nodes.xml", "w")
70 70  
71   - # Set node boundaries
72   - elogmin = math.log10(float(emin))
73   - elogmax = math.log10(float(emax))
74   - elogbin = (elogmax - elogmin)/(float(enumbins)-1.0)
75   -
76   - # Fill arrays
77   - for i in range(int(enumbins)):
78   -
79   - # Compute energy
80   - energy = math.pow(10.0, i*elogbin+elogmin)
81   -
82   - # Compute scale (HESS Crab spectrum)
83   - scale = 3.45e-17 * math.pow(energy, -2.63)
84   -
85   - # Write node
86   - write_node(file, energy, scale)
87   -
88   - # Debug
89   - #sys.stdout.write(energy)
  71 + # Set node boundaries
  72 + elogmin = math.log10(float(emin))
  73 + elogmax = math.log10(float(emax))
  74 + elogbin = (elogmax - elogmin)/(float(enumbins)-1.0)
  75 +
  76 + # Fill arrays
  77 + for i in range(int(enumbins)):
  78 +
  79 + # Compute energy
  80 + energy = math.pow(10.0, i*elogbin+elogmin)
  81 +
  82 + # Compute scale (HESS Crab spectrum)
  83 + scale = 3.45e-17 * math.pow(energy, -2.63)
  84 +
  85 + # Write node
  86 + write_node(file, energy, scale)
  87 +
  88 + # Debug
  89 + #sys.stdout.write(energy)
90 90 #sys.stdout.write(math.pow(10.0, int(math.log10(scale))-1.0)±"\n")
91   -
92   - # Close file
93   - file.close()
94   -
95   - # Return
96   - return
97   -
  91 +
  92 + # Close file
  93 + file.close()
  94 +
  95 + # Return
  96 + return
  97 +
98 98  
99 99 # ======================= #
100 100 # Main script entry point #
101 101 # ======================= #
102 102 if __name__ == '__main__':
103   - """
104   - Compute nodes.
105   - """
106   - # Check command line
107   - usage = "Usage: cta_make_nodes emin emax enumbins"
108   - if len(sys.argv) < 3:
109   - sys.stdout.write(usage+"\n")
110   - sys.exit()
111   -
112   - # Extract parameters
113   - emin = sys.argv[1]
114   - emax = sys.argv[2]
115   - enumbins = sys.argv[3]
  103 + """
  104 + Compute nodes.
  105 + """
  106 + # Check command line
  107 + usage = "Usage: cta_make_nodes emin emax enumbins"
  108 + if len(sys.argv) < 3:
  109 + sys.stdout.write(usage+"\n")
  110 + sys.exit()
  111 +
  112 + # Extract parameters
  113 + emin = sys.argv[1]
  114 + emax = sys.argv[2]
  115 + enumbins = sys.argv[3]
116 116  
117   - # Create nodes
118   - create_nodes(emin, emax, enumbins)
119   -
  117 + # Create nodes
  118 + create_nodes(emin, emax, enumbins)
  119 +
... ...
scripts/cta_root2irf.py
  View file @d9a9c73
... ... @@ -11,7 +11,7 @@
11 11 # area (i.e. the area before applying any theta cut). This is required by
12 12 # GammaLib as the maximum likelihood analysis does not imply any theta cut.
13 13 #
14   -# Copyright (C) 2011 Jurgen Knodlseder
  14 +# Copyright (C) 2011-2014 Juergen Knoedlseder
15 15 #
16 16 # This program is free software: you can redistribute it and/or modify
17 17 # it under the terms of the GNU General Public License as published by
... ...