Fabio Iocco / ctools

Commit 1566d76aa200e1f7e420c497798e46fa4e73a427

Authored by Jürgen Knödlseder
2 parents 299899f4 7955dfc9

Merge branch '1364-write-csspec' into integration

Inline Side-by-side
Showing 2 changed files with 266 additions and 41 deletions
scripts/csspec.py
  View file @1566d76
... ... @@ -21,10 +21,11 @@
21 21 import gammalib
22 22 import ctools
23 23 import sys
  24 +import tempfile
24 25  
25   -# ============== #
26   -# cstsdist class #
27   -# ============== #
  26 +# ============ #
  27 +# csspec class #
  28 +# ============ #
28 29 class csspec(ctools.cscript):
29 30 """
30 31 This class implements the creation of spectral points. It derives from
... ... @@ -42,14 +43,14 @@ class csspec(ctools.cscript):
42 43 self.version = "1.0.0"
43 44  
44 45 # Initialise some members
45   - self.obs = None
  46 + self.obs = None
46 47  
47 48 # Initialise some members
48 49 if isinstance(argv[0],gammalib.GObservations):
49 50 self.obs = argv[0]
50   - argv = argv[1:]
  51 + argv = argv[1:]
51 52 else:
52   - self.obs = gammalib.GObservations()
  53 + self.obs = gammalib.GObservations()
53 54 self.obs.clear()
54 55 self.m_outfile = ""
55 56  
... ... @@ -104,12 +105,22 @@ class csspec(ctools.cscript):
104 105 pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.fits","","","Output file name"))
105 106 pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
106 107 pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
  108 + pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)"))
107 109 pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
108 110 pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
109   - pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit (TeV)"))
110   - pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit (TeV)"))
111   - pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of energy bins"))
112   - pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm"))
  111 + pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit for spectral points(TeV)"))
  112 + pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit for spectral points(TeV)"))
  113 + pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of spectral points"))
  114 + pars.append(gammalib.GApplicationPar("ebinalg","s", "h","LOG","FILE|LIN|LOG","", "Algorithm for defining energy bins"))
  115 + pars.append(gammalib.GApplicationPar("binned","b","a","no","yes|no","","Use binned analysis in each energy bin"))
  116 + pars.append(gammalib.GApplicationPar("nebins","i","h","5","","","Number of energy bins per spectral point"))
  117 + pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
  118 + pars.append(gammalib.GApplicationPar("proj","s","a","CAR","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method"))
  119 + pars.append(gammalib.GApplicationPar("xref","r","a","83.63","0","360","First coordinate of image center in degrees (RA or galactic l)"))
  120 + pars.append(gammalib.GApplicationPar("yref","r","a","22.01","-90","90","Second coordinate of image center in degrees (DEC or galactic b)"))
  121 + pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels"))
  122 + pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels"))
  123 + pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)"))
113 124 pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin"))
114 125 pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin"))
115 126 pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters"))
... ... @@ -125,7 +136,6 @@ class csspec(ctools.cscript):
125 136 """
126 137 Get parameters from parfile and setup the observation.
127 138 """
128   - # Get parameters
129 139 # Set observation if not done before
130 140 if self.obs == None or self.obs.size() == 0:
131 141 self.require_inobs("csspec::get_parameters()")
... ... @@ -138,9 +148,21 @@ class csspec(ctools.cscript):
138 148 # Get source name
139 149 self.m_srcname = self["srcname"].string()
140 150  
141   - # Get ebounds
  151 + # Get ebounds
142 152 self.m_ebounds = self.create_ebounds()
143 153  
  154 + # Get binning flag
  155 + self.m_binned = self["binned"].boolean()
  156 + if self.m_binned:
  157 + self.m_xref = self["xref"].real()
  158 + self.m_yref = self["yref"].real()
  159 + self.m_nxpix = self["nxpix"].integer()
  160 + self.m_nypix = self["nypix"].integer()
  161 + self.m_binsz = self["binsz"].real()
  162 + self.m_coordsys = self["coordsys"].string()
  163 + self.m_proj = self["proj"].string()
  164 + self.m_ebins = self["nebins"].integer()
  165 +
144 166 # Read other parameters
145 167 self.m_outfile = self["outfile"].filename()
146 168  
... ... @@ -210,21 +232,39 @@ class csspec(ctools.cscript):
210 232 self.log("\n")
211 233 self.log.header1("Adjust model parameters")
212 234  
  235 + # Adjust model parameters dependent on input user parameters
213 236 for model in self.obs.models():
  237 +
  238 + # Set TS flag for all models to false.
  239 + # Source of interest will be set to true later
214 240 model.tscalc(False)
  241 +
  242 + # Log model name
  243 + if self.logExplicit():
  244 + self.log.header3(model.name())
  245 +
  246 + # Deal with the source of interest
215 247 if model.name() == self.m_srcname:
216 248 for par in model:
  249 + if par.is_free() and self.logExplicit():
  250 + self.log(" Fixing \""+par.name()+"\"\n")
217 251 par.fix()
  252 + if par.is_fixed() and self.logExplicit():
  253 + self.log(" Freeing \""+par.name()+"\"\n")
218 254 model.spectral()[0].free()
219 255 if self.m_calc_ts:
220 256 model.tscalc(True)
221 257  
222   - elif self.m_fix_bkg and not model.classname() == "GModelData":
  258 + elif self.m_fix_bkg and not model.classname() == "GModelSky":
223 259 for par in model:
  260 + if par.is_free() and self.logExplicit():
  261 + self.log(" Fixing \""+par.name()+"\"\n")
224 262 par.fix()
225 263  
226 264 elif self.m_fix_srcs and model.classname() == "GModelSky":
227 265 for par in model:
  266 + if par.is_free() and self.logExplicit():
  267 + self.log(" Fixing \""+par.name()+"\"\n")
228 268 par.fix()
229 269  
230 270 # Write header
... ... @@ -237,14 +277,16 @@ class csspec(ctools.cscript):
237 277 table.extname("SPECTRUM")
238 278  
239 279 # Add Header for compatibility with gammalib.GMWLSpectrum
240   - table.card("Instrument", "CTA", "Name of Instrument")
241   - table.card("Telescope", "CTA", "Name of Telescope")
  280 + table.card("INSTRUME", "CTA", "Name of Instrument")
  281 + table.card("TELESCOP", "CTA", "Name of Telescope")
242 282  
243 283 # Create FITS table columns
244 284 energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size())
245 285 energy.unit("TeV")
246   - energy_err = gammalib.GFitsTableDoubleCol("e_Energy", self.m_ebounds.size())
247   - energy_err.unit("TeV")
  286 + energy_low = gammalib.GFitsTableDoubleCol("ed_Energy", self.m_ebounds.size())
  287 + energy_low.unit("TeV")
  288 + energy_high = gammalib.GFitsTableDoubleCol("eu_Energy", self.m_ebounds.size())
  289 + energy_high.unit("TeV")
248 290 flux = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size())
249 291 flux.unit("erg/cm2/s")
250 292 flux_err = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size())
... ... @@ -252,16 +294,31 @@ class csspec(ctools.cscript):
252 294 TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size())
253 295 ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size())
254 296 ulim_values.unit("erg/cm2/s")
255   -
  297 +
256 298 # Loop over energy bins
257 299 for i in range(self.m_ebounds.size()):
258 300  
  301 + # Log information
  302 + if self.logTerse():
  303 + self.log("\n")
  304 + self.log.header2("Energy bin "+str(i))
  305 +
259 306 # Get energy boundaries
260 307 emin = self.m_ebounds.emin(i)
261 308 emax = self.m_ebounds.emax(i)
262   - emean = self.m_ebounds.emean(i)
263 309 elogmean = self.m_ebounds.elogmean(i)
264   - elogmean2 = elogmean.MeV() * elogmean.MeV()
  310 + elogmean2 = elogmean.MeV() * elogmean.MeV()
  311 +
  312 + # Store energy as TeV
  313 + energy[i] = elogmean.TeV()
  314 +
  315 + # Store energy errors
  316 + energy_low[i] = (elogmean - emin).TeV()
  317 + energy_high[i] = (emax - elogmean).TeV()
  318 +
  319 + # Log information
  320 + if self.logExplicit():
  321 + self.log.header3("Selecting events")
265 322  
266 323 # Select events
267 324 select = ctools.ctselect(self.obs)
... ... @@ -272,23 +329,163 @@ class csspec(ctools.cscript):
272 329 select["rad"].value("UNDEFINED")
273 330 select["ra"].value("UNDEFINED")
274 331 select["dec"].value("UNDEFINED")
275   - select.run()
276   -
277   - # likelihood
278   - like = ctools.ctlike(select.obs())
  332 + select.run()
  333 +
  334 + # Retrieve observation
  335 + obs = select.obs()
  336 +
  337 + # Binned analysis
  338 + if self.m_binned:
  339 +
  340 + # Header
  341 + if self.logTerse():
  342 + self.log.header3("Binning events")
  343 +
  344 + # Bin events
  345 + bin = ctools.ctbin(select.obs())
  346 + bin["usepnt"].boolean(False)
  347 + bin["ebinalg"].string("LOG")
  348 + bin["xref"].real(self.m_xref)
  349 + bin["yref"].real(self.m_yref)
  350 + bin["binsz"].real(self.m_binsz)
  351 + bin["nxpix"].integer(self.m_nxpix)
  352 + bin["nypix"].integer(self.m_nypix)
  353 + bin["enumbins"].integer(self.m_ebins)
  354 + bin["emin"].real(emin.TeV())
  355 + bin["emax"].real(emax.TeV())
  356 + bin["coordsys"].string(self.m_coordsys)
  357 + bin["proj"].string(self.m_proj)
  358 + bin.run()
  359 +
  360 + # Header
  361 + if self.logTerse():
  362 + self.log.header3("Creating exposure cube")
  363 +
  364 + # Create exposure cube
  365 + expcube = ctools.ctexpcube(select.obs())
  366 + expcube["incube"].filename("NONE")
  367 + expcube["usepnt"].boolean(False)
  368 + expcube["ebinalg"].string("LOG")
  369 + expcube["xref"].real(self.m_xref)
  370 + expcube["yref"].real(self.m_yref)
  371 + expcube["binsz"].real(self.m_binsz)
  372 + expcube["nxpix"].integer(self.m_nxpix)
  373 + expcube["nypix"].integer(self.m_nypix)
  374 + expcube["enumbins"].integer(self.m_ebins)
  375 + expcube["emin"].real(emin.TeV())
  376 + expcube["emax"].real(emax.TeV())
  377 + expcube["coordsys"].string(self.m_coordsys)
  378 + expcube["proj"].string(self.m_proj)
  379 + expcube.run()
  380 +
  381 + # Header
  382 + if self.logTerse():
  383 + self.log.header3("Creating PSF cube")
  384 +
  385 + # Create psf cube
  386 + psfcube = ctools.ctpsfcube(select.obs())
  387 + psfcube["incube"].filename("NONE")
  388 + psfcube["usepnt"].boolean(False)
  389 + psfcube["ebinalg"].string("LOG")
  390 + psfcube["xref"].real(self.m_xref)
  391 + psfcube["yref"].real(self.m_yref)
  392 + psfcube["binsz"].real(self.m_binsz)
  393 + psfcube["nxpix"].integer(self.m_nxpix)
  394 + psfcube["nypix"].integer(self.m_nypix)
  395 + psfcube["enumbins"].integer(self.m_ebins)
  396 + psfcube["emin"].real(emin.TeV())
  397 + psfcube["emax"].real(emax.TeV())
  398 + psfcube["coordsys"].string(self.m_coordsys)
  399 + psfcube["proj"].string(self.m_proj)
  400 + psfcube.run()
  401 +
  402 + # Header
  403 + if self.logTerse():
  404 + self.log.header3("Creating background cube")
  405 +
  406 + # Create background cube
  407 + bkgcube = ctools.ctbkgcube(select.obs())
  408 + bkgcube["incube"].filename("NONE")
  409 + bkgcube["usepnt"].boolean(False)
  410 + bkgcube["ebinalg"].string("LOG")
  411 + bkgcube["xref"].real(self.m_xref)
  412 + bkgcube["yref"].real(self.m_yref)
  413 + bkgcube["binsz"].real(self.m_binsz)
  414 + bkgcube["nxpix"].integer(self.m_nxpix)
  415 + bkgcube["nypix"].integer(self.m_nypix)
  416 + bkgcube["enumbins"].integer(self.m_ebins)
  417 + bkgcube["emin"].real(emin.TeV())
  418 + bkgcube["emax"].real(emax.TeV())
  419 + bkgcube["coordsys"].string(self.m_coordsys)
  420 + bkgcube["proj"].string(self.m_proj)
  421 + bkgcube.run()
  422 +
  423 + # Set new binned observation
  424 + obs = bin.obs()
  425 +
  426 + # Set precomputed binned response
  427 + obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.cube())
  428 +
  429 + # Get new models
  430 + models = bkgcube.models()
  431 +
  432 + # Fix background models if required
  433 + if self.m_fix_bkg:
  434 + for model in models:
  435 + if not model.classname() == "GModelSky":
  436 + for par in model:
  437 + par.fix()
  438 +
  439 + # Set new models to binned observation
  440 + obs.models(models)
  441 +
  442 + # Header
  443 + if self.logTerse():
  444 + self.log.header3("Performing fit")
  445 +
  446 + # Likelihood
  447 + like = ctools.ctlike(obs)
279 448 like.run()
280 449  
  450 + # Skip bin if no event was present
  451 + if like.obs().logL() == 0.0:
  452 +
  453 + # Log information
  454 + if self.logTerse():
  455 + self.log("No event in this bin. Bin is skipped\n")
  456 +
  457 + # Set all values to 0
  458 + flux[i] = 0.0
  459 + flux_err[i] = 0.0
  460 + TSvalues[i] = 0.0
  461 + ulim_values[i] = 0.0
  462 + continue
  463 +
281 464 # Get results
282 465 fitted_models = like.obs().models()
283 466 source = fitted_models[self.m_srcname]
284   -
  467 +
285 468 # Calculate Upper Limit
286 469 ulimit_value = -1.0
287   - #if self.m_calc_ulimit:
288   - # ulimit = ctools.ctulimit(like.obs())
289   - # ulimit["eref"].real(elogmean.TeV())
290   - # ulimit.run()
291   - # ulimit_value.ulimit.diff_ulimit()
  470 + if self.m_calc_ulimit:
  471 +
  472 + # Logging information
  473 + if self.logTerse():
  474 + self.log.header3("Computing upper limit")
  475 +
  476 + # Create upper limit object
  477 + ulimit = ctools.ctulimit(like.obs())
  478 + ulimit["srcname"].string(self.m_srcname)
  479 + ulimit["eref"].real(elogmean.TeV())
  480 +
  481 + # Try to run upper limit and catch exceptions
  482 + try:
  483 + ulimit.run()
  484 + ulimit_value = ulimit.diff_ulimit()
  485 + except:
  486 + if self.logTerse():
  487 + self.log("Upper limit calculation failed\n")
  488 + ulimit_value = -1.0
292 489  
293 490 # Get TS value
294 491 TS = -1.0
... ... @@ -306,19 +503,12 @@ class csspec(ctools.cscript):
306 503 # Set values for storage
307 504 TSvalues[i] = TS
308 505  
309   - # Store energy as TeV
310   - energy[i] = elogmean.TeV()
311   -
312   - # Store energy error
313   - energy_err[i] = (emax - emean).TeV()
314   -
315 506 # Convert fluxes to nuFnu
316 507 flux[i] = fitted_flux * elogmean2 * gammalib.MeV2erg
317 508 flux_err[i] = e_flux * elogmean2 * gammalib.MeV2erg
318 509 if ulimit_value > 0.0:
319 510 ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg
320   -
321   -
  511 +
322 512 # Log information
323 513 if self.logExplicit():
324 514 self.log("Bin "+str(i)+" ["+str(emin.TeV())+"-"+str(emax.TeV())+"] TeV: ")
... ... @@ -329,10 +519,11 @@ class csspec(ctools.cscript):
329 519 if self.m_calc_ulimit and ulim_values[i] > 0.0:
330 520 self.log(", UL = "+str(ulim_values[i])+" [erg/cm2/s]")
331 521 self.log("\n")
332   -
  522 +
333 523 # Append filles columns to fits table
334 524 table.append(energy)
335   - table.append(energy_err)
  525 + table.append(energy_low)
  526 + table.append(energy_high)
336 527 table.append(flux)
337 528 table.append(flux_err)
338 529 table.append(TSvalues)
... ...
test/test_cscripts.sh.in
  View file @1566d76
... ... @@ -191,7 +191,7 @@ $ECHO " ok"
191 191  
192 192 #
193 193 # Test csspec
194   -# =============
  194 +# ===========
195 195 $ECHO -n "Test csspec: "
196 196 chmod +x @abs_top_srcdir@/scripts/csspec.py
197 197 @abs_top_srcdir@/scripts/csspec.py \
... ... @@ -203,6 +203,7 @@ chmod +x @abs_top_srcdir@/scripts/csspec.py
203 203 srcname="Crab" \
204 204 emin=0.1 \
205 205 emax=100.0 \
  206 + binned=no \
206 207 outfile="spectrum.fits" > csspec.log
207 208 $ECHO -n "."
208 209 if [ -s "spectrum.fits" ]
... ... @@ -214,4 +215,37 @@ else
214 215 fi
215 216 $ECHO " ok"
216 217  
  218 +#
  219 +# Test csspec Run2
  220 +# ================
  221 +$ECHO -n "Test csspec: "
  222 +chmod +x @abs_top_srcdir@/scripts/csspec.py
  223 +@abs_top_srcdir@/scripts/csspec.py \
  224 + inobs="data/crab_events.fits" \
  225 + inmodel="data/crab.xml" \
  226 + caldb="irf" \
  227 + irf="cta_dummy_irf" \
  228 + enumbins=3 \
  229 + srcname="Crab" \
  230 + emin=0.1 \
  231 + emax=100.0 \
  232 + binned=yes \
  233 + nxpix=10 \
  234 + nypix=10 \
  235 + binsz=0.1 \
  236 + nebins=2 \
  237 + coordsys="CEL" \
  238 + xref=83.63 \
  239 + yref=22.01 \
  240 + proj="CAR" \
  241 + outfile="spectrum.fits" > csspec.log
  242 +$ECHO -n "."
  243 +if [ -s "spectrum.fits" ]
  244 +then
  245 + $ECHO -n "."
  246 +else
  247 + $ECHO " spectrum.fits file is not valid"
  248 + exit 1
  249 +fi
  250 +$ECHO " ok"
217 251  
... ...