Thierry Reposeur / ctools

Commit 50d6ad8850eb912ac8f5d7405bdb7faa952b18df

Authored by Michael Mayer
1 parent d43f11d9

csspec now supports binned analysis in each energy bin

Inline Side-by-side
Showing 2 changed files with 247 additions and 29 deletions
scripts/csspec.py
  View file @50d6ad8
... ... @@ -21,6 +21,7 @@
21 21 import gammalib
22 22 import ctools
23 23 import sys
  24 +import tempfile
24 25  
25 26 # ============== #
26 27 # cstsdist class #
... ... @@ -106,10 +107,19 @@ class csspec(ctools.cscript):
106 107 pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
107 108 pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
108 109 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"))
  110 + pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit for spectral points(TeV)"))
  111 + pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit for spectral points(TeV)"))
  112 + pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of spectral points"))
  113 + pars.append(gammalib.GApplicationPar("ebinalg","s", "h","LOG","FILE|LIN|LOG","", "Algorithm for defining energy bins"))
  114 + pars.append(gammalib.GApplicationPar("binned","b","a","no","yes|no","","Use binned analysis in each energy bin"))
  115 + pars.append(gammalib.GApplicationPar("nebins","i","h","5","","","Number of energy bins per spectral point"))
  116 + pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
  117 + pars.append(gammalib.GApplicationPar("proj","s","a","CAR","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method"))
  118 + pars.append(gammalib.GApplicationPar("xref","r","a","83.63","0","360","First coordinate of image center in degrees (RA or galactic l)"))
  119 + pars.append(gammalib.GApplicationPar("yref","r","a","22.01","-90","90","Second coordinate of image center in degrees (DEC or galactic b)"))
  120 + pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels"))
  121 + pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels"))
  122 + pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)"))
113 123 pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin"))
114 124 pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin"))
115 125 pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel 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
... ... @@ -243,8 +283,10 @@ class csspec(ctools.cscript):
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,156 @@ 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()
  332 + select.run()
276 333  
  334 + obs = select.obs()
  335 +
  336 + if self.m_binned:
  337 +
  338 + if self.logTerse():
  339 + self.log.header3("Binning events")
  340 +
  341 + # bin events
  342 + bin = ctools.ctbin(select.obs())
  343 + bin["usepnt"].boolean(False)
  344 + bin["ebinalg"].string("LOG")
  345 + bin["xref"].real(self.m_xref)
  346 + bin["yref"].real(self.m_yref)
  347 + bin["binsz"].real(self.m_binsz)
  348 + bin["nxpix"].integer(self.m_nxpix)
  349 + bin["nypix"].integer(self.m_nypix)
  350 + bin["enumbins"].integer(self.m_ebins)
  351 + bin["emin"].real(emin.TeV())
  352 + bin["emax"].real(emax.TeV())
  353 + bin["coordsys"].string(self.m_coordsys)
  354 + bin["proj"].string(self.m_proj)
  355 + bin.run()
  356 +
  357 + if self.logTerse():
  358 + self.log.header3("Creating exposure cube")
  359 +
  360 + # create exposure cube
  361 + expcube = ctools.ctexpcube(select.obs())
  362 + expcube["incube"].filename("NONE")
  363 + expcube["usepnt"].boolean(False)
  364 + expcube["ebinalg"].string("LOG")
  365 + expcube["xref"].real(self.m_xref)
  366 + expcube["yref"].real(self.m_yref)
  367 + expcube["binsz"].real(self.m_binsz)
  368 + expcube["nxpix"].integer(self.m_nxpix)
  369 + expcube["nypix"].integer(self.m_nypix)
  370 + expcube["enumbins"].integer(self.m_ebins)
  371 + expcube["emin"].real(emin.TeV())
  372 + expcube["emax"].real(emax.TeV())
  373 + expcube["coordsys"].string(self.m_coordsys)
  374 + expcube["proj"].string(self.m_proj)
  375 + expcube.run()
  376 +
  377 + if self.logTerse():
  378 + self.log.header3("Creating PSF cube")
  379 +
  380 + # create psf cube
  381 + psfcube = ctools.ctpsfcube(select.obs())
  382 + psfcube["incube"].filename("NONE")
  383 + psfcube["usepnt"].boolean(False)
  384 + psfcube["ebinalg"].string("LOG")
  385 + psfcube["xref"].real(self.m_xref)
  386 + psfcube["yref"].real(self.m_yref)
  387 + psfcube["binsz"].real(self.m_binsz)
  388 + psfcube["nxpix"].integer(self.m_nxpix)
  389 + psfcube["nypix"].integer(self.m_nypix)
  390 + psfcube["enumbins"].integer(self.m_ebins)
  391 + psfcube["emin"].real(emin.TeV())
  392 + psfcube["emax"].real(emax.TeV())
  393 + psfcube["coordsys"].string(self.m_coordsys)
  394 + psfcube["proj"].string(self.m_proj)
  395 + psfcube.run()
  396 +
  397 + if self.logTerse():
  398 + self.log.header3("Creating background cube")
  399 +
  400 + # create background cube
  401 + bkgcube = ctools.ctbkgcube(select.obs())
  402 + bkgcube["incube"].filename("NONE")
  403 + bkgcube["usepnt"].boolean(False)
  404 + bkgcube["ebinalg"].string("LOG")
  405 + bkgcube["xref"].real(self.m_xref)
  406 + bkgcube["yref"].real(self.m_yref)
  407 + bkgcube["binsz"].real(self.m_binsz)
  408 + bkgcube["nxpix"].integer(self.m_nxpix)
  409 + bkgcube["nypix"].integer(self.m_nypix)
  410 + bkgcube["enumbins"].integer(self.m_ebins)
  411 + bkgcube["emin"].real(emin.TeV())
  412 + bkgcube["emax"].real(emax.TeV())
  413 + bkgcube["coordsys"].string(self.m_coordsys)
  414 + bkgcube["proj"].string(self.m_proj)
  415 + bkgcube.run()
  416 +
  417 + # Set new binned observation
  418 + obs = bin.obs()
  419 +
  420 + # Set precomputed binned response
  421 + obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.cube())
  422 +
  423 + # Get new models
  424 + models = bkgcube.models()
  425 +
  426 + # Fix background models if required
  427 + if self.m_fix_bkg:
  428 + for model in models:
  429 + if not model.classname() == "GModelSky":
  430 + for par in model:
  431 + par.fix()
  432 +
  433 + # Set new models to binned observation
  434 + obs.models(models)
  435 +
  436 + if self.logTerse():
  437 + self.log.header3("Performing fit")
  438 +
277 439 # likelihood
278   - like = ctools.ctlike(select.obs())
  440 + like = ctools.ctlike(obs)
279 441 like.run()
280 442  
  443 + # Skip bin if no event was present
  444 + if like.obs().logL() == 0.0:
  445 +
  446 + # Log information
  447 + if self.logTerse():
  448 + self.log("No event in this bin. Bin is skipped\n")
  449 +
  450 + flux[i] = 0.0
  451 + flux_err[i] = 0.0
  452 + TSvalues[i] = 0.0
  453 + ulim_values[i] = 0.0
  454 + continue
  455 +
281 456 # Get results
282 457 fitted_models = like.obs().models()
283 458 source = fitted_models[self.m_srcname]
284   -
  459 +
285 460 # Calculate Upper Limit
286 461 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()
  462 + if self.m_calc_ulimit:
  463 +
  464 + # Logging information
  465 + if self.logTerse():
  466 + self.log.header3("Computing upper limit")
  467 +
  468 + # Create upper limit object
  469 + ulimit = ctools.ctulimit(like.obs())
  470 + ulimit["srcname"].string(self.m_srcname)
  471 + ulimit["eref"].real(elogmean.TeV())
  472 +
  473 + # Try to run upper limit and catch exceptions
  474 + try:
  475 + ulimit.run()
  476 + ulimit_value = ulimit.diff_ulimit()
  477 + except:
  478 +
  479 + if self.logTerse():
  480 + self.log("Upper limit calculation failed\n")
  481 + ulimit_value = -1.0
292 482  
293 483 # Get TS value
294 484 TS = -1.0
... ... @@ -306,19 +496,12 @@ class csspec(ctools.cscript):
306 496 # Set values for storage
307 497 TSvalues[i] = TS
308 498  
309   - # Store energy as TeV
310   - energy[i] = elogmean.TeV()
311   -
312   - # Store energy error
313   - energy_err[i] = (emax - emean).TeV()
314   -
315 499 # Convert fluxes to nuFnu
316 500 flux[i] = fitted_flux * elogmean2 * gammalib.MeV2erg
317 501 flux_err[i] = e_flux * elogmean2 * gammalib.MeV2erg
318 502 if ulimit_value > 0.0:
319 503 ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg
320   -
321   -
  504 +
322 505 # Log information
323 506 if self.logExplicit():
324 507 self.log("Bin "+str(i)+" ["+str(emin.TeV())+"-"+str(emax.TeV())+"] TeV: ")
... ... @@ -329,10 +512,11 @@ class csspec(ctools.cscript):
329 512 if self.m_calc_ulimit and ulim_values[i] > 0.0:
330 513 self.log(", UL = "+str(ulim_values[i])+" [erg/cm2/s]")
331 514 self.log("\n")
332   -
  515 +
333 516 # Append filles columns to fits table
334 517 table.append(energy)
335   - table.append(energy_err)
  518 + table.append(energy_low)
  519 + table.append(energy_high)
336 520 table.append(flux)
337 521 table.append(flux_err)
338 522 table.append(TSvalues)
... ...
test/test_cscripts.sh.in
  View file @50d6ad8
... ... @@ -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  
... ...