Johan Bregeon / ctools

Commit 8d8bb720c913ea606da6d4d5cefaa743fdf5be73

Authored by Jürgen Knödlseder
1 parent f497c92a

Use input model for sensitivity computation and simplify the parameter interface (#1426)

Inline Side-by-side
Showing 6 changed files with 292 additions and 247 deletions
ChangeLog
  View file @8d8bb72
1   -2015-02-16 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2015-02-20 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.0.0 released
4 4 ========================
... ... @@ -6,6 +6,7 @@
6 6 Adapt to GammaLib 1.0.0 interface
7 7 Do not import automatically Python scripts
8 8 Support INDEF/NONE/UNDEF/UNDEFINED parameters in ctselect
  9 + Modify cssens parameter interface (#1426)
9 10  
10 11  
11 12 2015-01-08 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 1.0.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -16 February 2015
  5 +20 February 2015
6 6  
7 7  
8 8 1. Introduction
... ... @@ -99,7 +99,10 @@ None
99 99  
100 100 18. Scripts
101 101 -----------
102   -None
  102 +The user parameter interface of cssens has been modified. The script
  103 +now reads a source model that is used for sensitivity estimation.
  104 +This provides more flexibility in the choice of the model for
  105 +sensitivity estimation.
103 106  
104 107  
105 108 19. caldb
... ...
doc/source/reference_manual/cssens.rst
  View file @8d8bb72
... ... @@ -12,19 +12,38 @@ Synopsis
12 12 Computes the CTA sensitivity for a number of energy bins using :doc:`ctlike`.
13 13 Crab spectra are fitted in narrow energy bins to simulated data, and the
14 14 flux level is determined that leads to a particular significance.
  15 +The significance is estimated using the Test Statistic value.
  16 +The simplified assumption is made that the significance (in Gaussian
  17 +sigma) is the square root of the Test Statistic.
15 18  
16 19  
17 20 General parameters
18 21 ------------------
19 22  
20   -``outfile [file]``
21   - ASCII file containing the sensitivity values.
22   -
  23 +``(inobs = NONE) [file]``
  24 + Input event list, counts cube or observation definition XML file.
  25 +
  26 +``inmodel [file]``
  27 + Input source model XML file.
  28 +
  29 +``srcname [string]``
  30 + Name of the source in the source model XML file which should be used
  31 + for sensitivity computation.
  32 +
  33 +``(offset = 0.0) [real]``
  34 + Offset angle of source in field of view (in degrees).
  35 +
23 36 ``caldb [string]``
24 37 Calibration database.
25 38  
26 39 ``irf [string]``
27 40 Instrumental response function.
  41 +
  42 +``(deadc = 0.95) [real]``
  43 + Average deadtime correction factor.
  44 +
  45 +``outfile [file]``
  46 + ASCII file containing the sensitivity values.
28 47  
29 48 ``type <point|gauss|shell|disk> [string]``
30 49 Source model type.
... ... @@ -32,14 +51,20 @@ General parameters
32 51 ``offset [real]``
33 52 Source offset angle (in degrees).
34 53  
35   -``bkg [string]``
36   - Background model file function (none=power law for subarray E).
37   -
38 54 ``duration [real]``
39 55 Effective exposure time (in seconds).
40 56  
41 57 ``rad [real]``
42   - Radius of ROI (in degrees).
  58 + Radius of Region of Interest (RoI) (in degrees).
  59 +
  60 +``emin [real]``
  61 + Lower energy limit for sensitivity computation (in TeV).
  62 +
  63 +``emax [real]``
  64 + Upper energy limit for sensitivity computation (in TeV).
  65 +
  66 +``bins [integer]``
  67 + Number of energy bins for sensitivity computation.
43 68  
44 69 ``(enumbins = 0) [integer]``
45 70 Number of energy bins (0=unbinned).
... ... @@ -53,18 +78,6 @@ General parameters
53 78 ``(sigma = 5.0) [real]``
54 79 Significance threshold.
55 80  
56   -``(ts_use = yes) [boolean]``
57   - Use TS for signifiance estimation?
58   -
59   -``(index = -2.48) [real]``
60   - Spectral index for source model.
61   -
62   -``(radius = 0.1) [real]``
63   - Extended source model radius (in degrees).
64   -
65   -``(width = 0.05) [real]``
66   - Extended source model width (in degrees).
67   -
68 81 ``(max_iter = 50) [integer]``
69 82 Maximum number of iterations.
70 83  
... ... @@ -96,6 +109,9 @@ Standard parameters
96 109 ``(mode = ql) [string]``
97 110 Mode of automatic parameters (default is "ql", i.e. "query and learn").
98 111  
  112 +``(logfile = cssens.log) [filename]``
  113 + Log filename.
  114 +
99 115  
100 116 Related tools
101 117 -------------
... ...
doc/source/reference_manual/cttsmap.rst
  View file @8d8bb72
... ... @@ -3,13 +3,13 @@
3 3 cttsmap
4 4 =======
5 5  
6   -Generate a Test Statistics map.
  6 +Generate a Test Statistic map.
7 7  
8 8  
9 9 Synopsis
10 10 --------
11 11  
12   -This tool generates a test statistics map.
  12 +This tool generates a Test Statistic map.
13 13  
14 14  
15 15 General parameters
... ...
scripts/cssens.py
  View file @8d8bb72
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This script computes the CTA sensitivity for a number of energy bins using
4   -# ctlike. Crab spectra are fitted in narrow energy bins to simulated data,
5   -# and the flux level is determined that leads to a particular significance.
6   -#
7   -# The significance can be either determined using the Test Statistic (which
8   -# is defined as twice the likelihood difference between fitting with and
9   -# without the test source) or using the source model scaling factor (or
10   -# Prefactor) divided by its statistical error. The latter has been turned
11   -# out to be more reliable, and is used as the default parameter.
12   -#
13   -# As background model, a GCTAModelRadialGauss model is used for the spatial
14   -# component. For the spectral component, either a file function is used
15   -# (specified by the "bkg" parameter), or if left blank, a power law adjusted
16   -# to Konrad's configuration "E" performance file is employed.
17   -#
18   -# The source is modelled as a point source
  3 +# This script computes the CTA sensitivity.
19 4 #
20 5 # Copyright (C) 2011-2015 Juergen Knoedlseder
21 6 #
... ... @@ -46,10 +31,12 @@ import math
46 31 # ============ #
47 32 class cssens(ctools.cscript):
48 33 """
49   - This class implements the sensitivity computation script. It derives
50   - from the ctools.cscript 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.
  34 + This class computes the CTA sensitivity for a number of energy bins using
  35 + ctlike. Spectra are fitted in narrow energy bins to simulated data,
  36 + and the flux level is determined that leads to a particular significance.
  37 +
  38 + The significance is determined using the Test Statistic, defined as twice
  39 + the likelihood difference between fitting with and without the test source.
53 40 """
54 41 def __init__(self, *argv):
55 42 """
... ... @@ -58,6 +45,11 @@ class cssens(ctools.cscript):
58 45 # Set name
59 46 self.name = "cssens"
60 47 self.version = "1.0.0"
  48 +
  49 + # Initialise some parameters
  50 + self.obs = None
  51 + self.m_ra = None
  52 + self.m_dec = None
61 53  
62 54 # Make sure that parfile exists
63 55 file = self.parfile()
... ... @@ -101,29 +93,31 @@ class cssens(ctools.cscript):
101 93 pars = gammalib.GApplicationPars(parfile)
102 94 except:
103 95 # Signal if parfile was not found
104   - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  96 + print("Parfile "+parfile+" not found. Create default parfile.")
105 97  
106 98 # Create default parfile
107 99 pars = gammalib.GApplicationPars()
108   - pars.append(gammalib.GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
109   - pars.append(gammalib.GApplicationPar("caldb","s","a","$CTOOLS/share/caldb/data/cta/dummy","","","Calibration database"))
  100 + pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file"))
  101 + pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
  102 + pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
  103 + pars.append(gammalib.GApplicationPar("offset","r","h","0.0","","","Source offset angle (deg)"))
  104 + pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
110 105 pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
111   - pars.append(gammalib.GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
112   - pars.append(gammalib.GApplicationPar("offset","r","a","0.0","","","Source offset angle (deg)"))
113   - pars.append(gammalib.GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
  106 + pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","0","1","Deadtime correction factor"))
  107 + pars.append(gammalib.GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name"))
114 108 pars.append(gammalib.GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
115 109 pars.append(gammalib.GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)"))
  110 + pars.append(gammalib.GApplicationPar("emin","r","a","0.020","","","Lower energy limit (TeV)"))
  111 + pars.append(gammalib.GApplicationPar("emax","r","a","200.0","","","Upper energy limit (TeV)"))
  112 + pars.append(gammalib.GApplicationPar("bins","i","a","21","","","Number of energy bins for differential sensitivity computation"))
116 113 pars.append(gammalib.GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)"))
117 114 pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
118 115 pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
119 116 pars.append(gammalib.GApplicationPar("sigma","r","h","5.0","","","Significance threshold"))
120   - pars.append(gammalib.GApplicationPar("ts_use","b","h","yes","","","Use TS?"))
121   - pars.append(gammalib.GApplicationPar("index","r","h","-2.48","","","Assumed spectral index"))
122   - pars.append(gammalib.GApplicationPar("radius","r","h","0.1","","","Extended source model radius"))
123   - pars.append(gammalib.GApplicationPar("width","r","h","0.05","","","Extended source model width"))
124 117 pars.append(gammalib.GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations"))
125 118 pars.append(gammalib.GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
126 119 pars.append_standard()
  120 + pars.append(gammalib.GApplicationPar("logfile","f","h","cssens.log","","","Log filename"))
127 121 pars.save(parfile)
128 122  
129 123 # Return
... ... @@ -131,31 +125,48 @@ class cssens(ctools.cscript):
131 125  
132 126 def get_parameters(self):
133 127 """
134   - Get parameters from parfile.
  128 + Get user parameters from parfile.
135 129 """
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()
  130 + # Set observation if not done before
  131 + if self.obs == None or self.obs.size() == 0:
  132 + self.obs = self.set_obs()
  133 +
  134 + # Set models if we have none
  135 + if self.obs.models().size() == 0:
  136 + self.obs.models(self["inmodel"].filename())
  137 +
  138 + # Get source name
  139 + self.m_srcname = self["srcname"].string()
  140 +
  141 + # Read further parameters
  142 + self.m_outfile = self["outfile"].filename()
  143 + self.m_emin = self["emin"].real()
  144 + self.m_emax = self["emax"].real()
  145 + self.m_bins = self["bins"].integer()
  146 +
  147 + # Read parameters for binned if requested
  148 + self.m_enumbins = self["enumbins"].integer()
  149 + if not self.m_enumbins == 0:
  150 + self.m_npix = self["npix"].integer()
  151 + self.m_binsz = self["binsz"].real()
  152 + else:
  153 + # Set dummy values (required by obsutils)
  154 + self.m_npix = 0
  155 + self.m_binsz = 0.0
  156 +
  157 + # Read remaining parameters
148 158 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 159 self.m_max_iter = self["max_iter"].integer()
154 160 self.m_num_avg = self["num_avg"].integer()
155   -
  161 +
156 162 # Set some fixed parameters
157 163 self.m_log = False # Logging in client tools
158 164 self.m_debug = self["debug"].boolean() # Debugging in client tools
  165 +
  166 + # Derive some parameters
  167 + self.m_ebounds = gammalib.GEbounds(self.m_bins, \
  168 + gammalib.GEnergy(self.m_emin, "TeV"), \
  169 + gammalib.GEnergy(self.m_emax, "TeV"))
159 170  
160 171 # Return
161 172 return
... ... @@ -190,17 +201,8 @@ class cssens(ctools.cscript):
190 201 colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens']
191 202 results = []
192 203  
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 = gammalib.GModels()
197   - full_model = gammalib.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 + # Initialise models
  205 + full_model, bkg_model = self.set_models()
204 206  
205 207 # Write models into logger
206 208 if self.logTerse():
... ... @@ -218,21 +220,19 @@ class cssens(ctools.cscript):
218 220 self.log("\n")
219 221 self.log.header1("Sensitivity determination")
220 222  
221   - # Loop over energy bands. The number of energy band is still
222   - # fixed ...
223   - for ieng in range(21):
  223 + # Loop over energy bins
  224 + for ieng in range(self.m_ebounds.size()):
224 225  
225 226 # 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)
  227 + emean = self.m_ebounds.elogmean(ieng).TeV()
  228 + emin = self.m_ebounds.emin(ieng).TeV()
  229 + emax = self.m_ebounds.emax(ieng).TeV()
230 230  
231   - # Setup observation(s)
232   - obs = self.set_obs(emin=emin, emax=emax)
  231 + # Set energy boundaries
  232 + self.set_obs_ebounds(emin, emax)
233 233  
234 234 # Determine sensitivity
235   - result = self.get_sensitivity(obs, bkg_model, full_model)
  235 + result = self.get_sensitivity(self.obs, bkg_model, full_model)
236 236  
237 237 # Write results
238 238 if ieng == 0:
... ... @@ -252,10 +252,10 @@ class cssens(ctools.cscript):
252 252  
253 253 # Return
254 254 return
255   -
  255 +
256 256 def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
257 257 """
258   - Returns an observation container with a single CTA observation.
  258 + Returns an observation container.
259 259  
260 260 Keywords:
261 261 lpnt - Galactic longitude of pointing [deg] (default: 0.0)
... ... @@ -263,134 +263,125 @@ class cssens(ctools.cscript):
263 263 emin - Minimum energy [TeV] (default: 0.1)
264 264 emax - Maximum energy [TeV] (default: 100.0)
265 265 """
266   - # Allocate observation container
267   - obs = gammalib.GObservations()
  266 + # If an observation was provided on input then load it from XML file
  267 + filename = self["inobs"].filename()
  268 + if filename != "NONE" and filename != "":
  269 + obs = self.get_observations()
  270 +
  271 + # ... otherwise allocate a single observation
  272 + else:
  273 +
  274 + # Read relevant user parameters
  275 + caldb = self["caldb"].string()
  276 + irf = self["irf"].string()
  277 + deadc = self["deadc"].real()
  278 + duration = self["duration"].real()
  279 + rad = self["rad"].real()
  280 +
  281 + # Allocate observation container
  282 + obs = gammalib.GObservations()
268 283  
269   - # Set single pointing
270   - pntdir = gammalib.GSkyDir()
271   - pntdir.lb_deg(lpnt, bpnt)
  284 + # Set single pointing
  285 + pntdir = gammalib.GSkyDir()
  286 + pntdir.lb_deg(lpnt, bpnt)
272 287  
273   - # Create CTA observation
274   - run = obsutils.set_obs(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
275   - duration=self.m_duration, \
276   - emin=emin, emax=emax, rad=self.m_roi)
  288 + # Create CTA observation
  289 + run = obsutils.set_obs(pntdir, caldb=caldb, irf=irf, \
  290 + duration=duration, deadc=deadc, \
  291 + emin=emin, emax=emax, rad=rad)
277 292  
278   - # Append observation to container
279   - obs.append(run)
  293 + # Append observation to container
  294 + obs.append(run)
  295 +
  296 + # Set source position
  297 + offset = self["offset"].real()
  298 + pntdir.lb_deg(lpnt, bpnt+offset)
  299 + self.m_ra = pntdir.ra_deg()
  300 + self.m_dec = pntdir.dec_deg()
280 301  
281 302 # Return observation container
282 303 return obs
283   -
284   - def set_bkg_model(self, fitidx=False):
  304 +
  305 + def set_obs_ebounds(self, emin, emax):
285 306 """
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)
  307 + Set energy boundaries for observation in container.
291 308 """
292   - # Define radial component
293   - radial = gammalib.GCTAModelRadialGauss(3.0)
  309 + # Loop over all observations in container
  310 + for obs in self.obs:
294 311  
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 = gammalib.GModelSpectralFunc(self.m_bkg, 1.0)
299   - else:
300   - spectrum = gammalib.GModelSpectralPlaw(61.8e-6, -1.85, gammalib.GEnergy(1.0, "TeV"))
301   - if fitidx:
302   - spectrum["Index"].free()
303   - else:
304   - spectrum["Index"].fix()
  312 + # Set energy boundaries
  313 + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \
  314 + gammalib.GEnergy(emax, "TeV"))
  315 + obs.events().ebounds(ebounds)
305 316  
306   - # Create background model
307   - model = gammalib.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   - 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).
  317 + # Return
  318 + return
320 319  
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)
  320 + def set_models(self, fitpos=False, fitspec=False):
330 321 """
331   - # Set source location
332   - location = gammalib.GSkyDir()
333   - location.lb_deg(l, b)
334   -
335   - # Set source spectrum
336   - spectrum = gammalib.GModelSpectralPlaw(flux*5.7e-16, index, gammalib.GEnergy(0.3, "TeV"))
337   - if fitidx:
338   - spectrum["Index"].free()
  322 + Set full and background model.
  323 + """
  324 + # Retrieve full model from observation container
  325 + full_model = self.obs.models().copy()
  326 +
  327 + # Get source model
  328 + model = full_model[self.m_srcname]
  329 +
  330 + # Check that model has a Prefactor
  331 + if not model.has_par("Prefactor"):
  332 + msg = "Model \""+self.m_srcname+"\" has no parameter \"Prefactor\"."+ \
  333 + " Only spectral models with a \"Prefactor\" parameter are supported."
  334 + raise gammalib.GException.invalid_value("cssens", msg)
  335 +
  336 + # Set source position
  337 + if self.m_ra != None and self.m_dec != None:
  338 + if model.has_par("RA") and model.has_par("DEC"):
  339 + model["RA"].value(self.m_ra)
  340 + model["DEC"].value(self.m_dec)
  341 +
  342 + # Fit or fix spatial parameters
  343 + if fitpos:
  344 + if model.has_par("RA"):
  345 + model["RA"].free()
  346 + if model.has_par("DEC"):
  347 + model["DEC"].free()
  348 + if model.has_par("Sigma"):
  349 + model["Sigma"].free()
  350 + if model.has_par("Radius"):
  351 + model["Radius"].free()
  352 + if model.has_par("Width"):
  353 + model["Width"].free()
339 354 else:
340   - spectrum["Index"].fix()
341   -
342   - # Set source
343   - if type == "point":
344   - spatial = gammalib.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 = gammalib.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 = gammalib.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 = gammalib.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()
  355 + if model.has_par("RA"):
  356 + model["RA"].fix()
  357 + if model.has_par("DEC"):
  358 + model["DEC"].fix()
  359 + if model.has_par("Sigma"):
  360 + model["Sigma"].fix()
  361 + if model.has_par("Radius"):
  362 + model["Radius"].fix()
  363 + if model.has_par("Width"):
  364 + model["Width"].fix()
  365 +
  366 + # Fit or fix spectral parameters
  367 + if fitspec:
  368 + if model.has_par("Index"):
  369 + model["Index"].free()
  370 + if model.has_par("Cutoff"):
  371 + model["Cutoff"].free()
383 372 else:
384   - self.log("ERROR: Unknown source type '"+type+"'.\n")
385   - return None
386   - source = gammalib.GModelSky(spatial, spectrum)
387   -
388   - # Set source name
389   - source.name("Test")
390   -
391   - # Return source
392   - return source
393   -
  373 + if model.has_par("Index"):
  374 + model["Index"].fix()
  375 + if model.has_par("Cutoff"):
  376 + model["Cutoff"].fix()
  377 +
  378 + # Create background model
  379 + bkg_model = full_model.copy()
  380 + bkg_model.remove(self.m_srcname)
  381 +
  382 + # Return models
  383 + return full_model, bkg_model
  384 +
394 385 def get_sensitivity(self, obs, bkg_model, full_model):
395 386 """
396 387 Determine sensitivity for a given observations.
... ... @@ -413,10 +404,25 @@ class cssens(ctools.cscript):
413 404 erg_width = (emax.TeV()-emin.TeV()) * tev2erg
414 405 break
415 406  
  407 + # Compute Crab unit (this is the factor with which the Prefactor needs
  408 + # to be multiplied to get 1 Crab
  409 + crab_flux = self.get_crab_flux(emin, emax)
  410 + src_flux = full_model[self.m_srcname].spectral().flux(emin, emax)
  411 + crab_unit = crab_flux/src_flux
  412 +
416 413 # Write header
417 414 if self.logTerse():
418 415 self.log("\n")
419 416 self.log.header2("Energies: "+str(emin)+" - "+str(emax))
  417 + self.log.parformat("Crab flux")
  418 + self.log(crab_flux)
  419 + self.log(" ph/cm2/s\n")
  420 + self.log.parformat("Source model flux")
  421 + self.log(src_flux)
  422 + self.log(" ph/cm2/s\n")
  423 + self.log.parformat("Crab unit factor")
  424 + self.log(crab_unit)
  425 + self.log("\n")
420 426  
421 427 # Initialise loop
422 428 flux_value = []
... ... @@ -424,7 +430,7 @@ class cssens(ctools.cscript):
424 430 eflux_value = []
425 431 diffSens_value = []
426 432 iter = 0
427   - test_flux = 0.1 # This is the initial test flux in Crab units
  433 + test_cflux = 0.1 # This is the initial test flux in Crab units (100 mCrab)
428 434  
429 435 # Loop until we break
430 436 while True:
... ... @@ -437,13 +443,14 @@ class cssens(ctools.cscript):
437 443 self.log.header2("Iteration "+str(iter))
438 444  
439 445 # Set test flux
440   - src_model = full_model.copy()
441   - crab_flux = src_model["Test"]['Prefactor'].value()
442   - src_model["Test"]['Prefactor'].value(crab_flux*test_flux)
  446 + src_model = full_model.copy()
  447 + crab_prefactor = src_model[self.m_srcname]['Prefactor'].value() * crab_unit
  448 + src_model[self.m_srcname]['Prefactor'].value(crab_prefactor * test_cflux)
443 449 obs.models(src_model)
444 450  
  451 +
445 452 # Simulate events
446   - sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \
  453 + sim = obsutils.sim(obs, nbins=self.m_enumbins, seed=iter, \
447 454 binsz=self.m_binsz, npix=self.m_npix, \
448 455 log=self.m_log, debug=self.m_debug)
449 456  
... ... @@ -544,22 +551,24 @@ class cssens(ctools.cscript):
544 551 continue
545 552  
546 553 # Get fitted Crab, photon and energy fluxes
547   - cflux = result_all["Test"]['Prefactor'].value() / crab_flux
548   - cflux_err = result_all["Test"]['Prefactor'].error() / crab_flux
549   - pflux = result_all["Test"].spectral().flux(emin, emax)
550   - eflux = result_all["Test"].spectral().eflux(emin, emax)
  554 + cflux = result_all[self.m_srcname]['Prefactor'].value() / crab_prefactor
  555 + cflux_err = result_all[self.m_srcname]['Prefactor'].error() / crab_prefactor
  556 + pflux = result_all[self.m_srcname].spectral().flux(emin, emax)
  557 + eflux = result_all[self.m_srcname].spectral().eflux(emin, emax)
551 558  
552 559 # Compute differential sensitivity
553 560 diffSens = pflux / erg_width * erg_mean*erg_mean
554 561  
555 562 # Compute flux correction factor based on average TS
556 563 correct = 1.0
557   - if self.m_use_ts:
558   - if ts > 0:
559   - correct = math.sqrt(self.m_ts_thres/ts)
560   - else:
561   - if cflux > 0 and cflux_err > 0:
562   - correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err)
  564 + if ts > 0:
  565 + correct = math.sqrt(self.m_ts_thres/ts)
  566 + #if self.m_use_ts:
  567 + # if ts > 0:
  568 + # correct = math.sqrt(self.m_ts_thres/ts)
  569 + #else:
  570 + # if cflux > 0 and cflux_err > 0:
  571 + # correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err)
563 572  
564 573 # Compute extrapolated fluxes
565 574 flux = correct * cflux
... ... @@ -631,8 +640,8 @@ class cssens(ctools.cscript):
631 640  
632 641 # Compare average flux to last average
633 642 if iter > self.m_num_avg:
634   - if test_flux > 0:
635   - ratio = flux/test_flux
  643 + if test_cflux > 0:
  644 + ratio = flux/test_cflux
636 645  
637 646 # We have 2 convergence criteria:
638 647 # 1. The average flux does not change
... ... @@ -648,7 +657,7 @@ class cssens(ctools.cscript):
648 657 break
649 658  
650 659 # Use average for next iteration
651   - test_flux = flux
  660 + test_cflux = flux
652 661  
653 662 # Exit loop if number of trials exhausted
654 663 if (iter >= self.m_max_iter):
... ... @@ -712,6 +721,19 @@ class cssens(ctools.cscript):
712 721 # Return result
713 722 return result
714 723  
  724 + def get_crab_flux(self, emin, emax):
  725 + """
  726 + Determine the Crab photon flux in a given energy interval.
  727 + """
  728 + # Set Crab spectral model
  729 + crab = gammalib.GModelSpectralPlaw(5.7e-16, -2.48, gammalib.GEnergy(0.3, "TeV"))
  730 +
  731 + # Determine flux
  732 + flux = crab.flux(emin, emax)
  733 +
  734 + # Return flux
  735 + return flux
  736 +
715 737  
716 738 # ======================== #
717 739 # Main routine entry point #
... ...
test/test_cscripts.sh.in
  View file @8d8bb72
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # This script tests all cscripts that are shipped with ctools
4 4 #
5   -# Copyright (C) 2011-2013 Juergen Knoedlseder
  5 +# Copyright (C) 2011-2015 Juergen Knoedlseder
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -66,6 +66,7 @@ chmod +x @abs_top_srcdir@/scripts/cspull.py
66 66 deadc=0.95 \
67 67 rad=5.0 \
68 68 npix=200 \
  69 + debug=yes \
69 70 binsz=0.05 > cspull.log
70 71 $ECHO -n "."
71 72 if [ -s "pull.dat" ]
... ... @@ -81,25 +82,27 @@ $ECHO &quot; ok&quot;
81 82 #
82 83 # Test cssens
83 84 # ===========
84   -#$ECHO -n "Test cssens: "
85   -#chmod +x @abs_top_srcdir@/scripts/cssens.py
86   -#@abs_top_srcdir@/scripts/cssens.py \
87   -# duration=3600.0 \
88   -# caldb="irf" \
89   -# irf="cta_dummy_irf" \
90   -# type="point" \
91   -# offset="0.0" \
92   -# bkg="data/bkg_dummy.txt" \
93   -# rad=5.0 > cssens.log
94   -#echo @ECHO_N@ "."
95   -#if [ -s "sensitivity.dat" ]
96   -#then
97   -# $ECHO -n "."
98   -#else
99   -# echo " sensitivity.dat file is not valid"
100   -# exit 1
101   -#fi
102   -#echo " ok"
  85 +$ECHO -n "Test cssens: "
  86 +chmod +x @abs_top_srcdir@/scripts/cssens.py
  87 +@abs_top_srcdir@/scripts/cssens.py \
  88 + inmodel="data/crab.xml" \
  89 + srcname="Crab" \
  90 + caldb="irf" \
  91 + irf="cta_dummy_irf" \
  92 + duration=1800.0 \
  93 + rad=3.0 \
  94 + emin=1.0 \
  95 + emax=10.0 \
  96 + bins=1 > cssens.log
  97 +$ECHO -n "."
  98 +if [ -s "sensitivity.dat" ]
  99 +then
  100 + $ECHO -n "."
  101 +else
  102 + $ECHO " sensitivity.dat file is not valid"
  103 + exit 1
  104 +fi
  105 +$ECHO " ok"
103 106  
104 107  
105 108 #
... ...