Thierry Reposeur / ctools

Commit 3b0fd743eeae240ef10d0365cc06501d2eebc6c5

Authored by Jürgen Knödlseder
2 parents 00b1324e 55f9af02

Merge branch 'integration' into devel

Inline Side-by-side
Showing 5 changed files with 79 additions and 39 deletions
cscripts/csspec.py
  View file @3b0fd74
... ... @@ -122,6 +122,7 @@ class csspec(ctools.cscript):
122 122 pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels"))
123 123 pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels"))
124 124 pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)"))
  125 + pars.append(gammalib.GApplicationPar("anumbins","i","h","200","","","Angular separation bins for PSF cube in binned mode"))
125 126 pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin"))
126 127 pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin"))
127 128 pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters"))
... ... @@ -167,6 +168,7 @@ class csspec(ctools.cscript):
167 168 self.m_coordsys = self["coordsys"].string()
168 169 self.m_proj = self["proj"].string()
169 170 self.m_ebins = self["nebins"].integer()
  171 + self.m_anumbins = self["anumbins"].integer()
170 172  
171 173 # Read other parameters
172 174 self.m_outfile = self["outfile"].filename()
... ... @@ -254,9 +256,10 @@ class csspec(ctools.cscript):
254 256 if par.is_free() and self.logExplicit():
255 257 self.log(" Fixing \""+par.name()+"\"\n")
256 258 par.fix()
257   - if par.is_fixed() and self.logExplicit():
258   - self.log(" Freeing \""+par.name()+"\"\n")
259   - model.spectral()[0].free()
  259 + normpar = model.spectral()[0]
  260 + if normpar.is_fixed() and self.logExplicit():
  261 + self.log(" Freeing \""+normpar.name()+"\"\n")
  262 + normpar.free()
260 263 if self.m_calc_ts:
261 264 model.tscalc(True)
262 265  
... ... @@ -308,7 +311,7 @@ class csspec(ctools.cscript):
308 311 # Log information
309 312 if self.logTerse():
310 313 self.log("\n")
311   - self.log.header2("Energy bin "+str(i))
  314 + self.log.header2("Energy bin "+str(i+1))
312 315  
313 316 # Get energy boundaries
314 317 emin = self.m_ebounds.emin(i)
... ... @@ -323,33 +326,36 @@ class csspec(ctools.cscript):
323 326 energy_low[i] = (elogmean - emin).TeV()
324 327 energy_high[i] = (emax - elogmean).TeV()
325 328  
326   - # Log information
327   - if self.logExplicit():
328   - self.log.header3("Selecting events")
329   -
330   - # Select events
331   - select = ctools.ctselect(self.obs)
332   - select["emin"] = emin.TeV()
333   - select["emax"] = emax.TeV()
334   - select["tmin"] = "UNDEFINED"
335   - select["tmax"] = "UNDEFINED"
336   - select["rad"] = "UNDEFINED"
337   - select["ra"] = "UNDEFINED"
338   - select["dec"] = "UNDEFINED"
339   - select.run()
340   -
341   - # Retrieve observation
342   - obs = select.obs()
  329 + # unbinned analysis
  330 + if not self.m_binned:
  331 +
  332 + # Log information
  333 + if self.logExplicit():
  334 + self.log.header3("Selecting events")
  335 +
  336 + # Select events
  337 + select = ctools.ctselect(self.obs)
  338 + select["emin"] = emin.TeV()
  339 + select["emax"] = emax.TeV()
  340 + select["tmin"] = "UNDEFINED"
  341 + select["tmax"] = "UNDEFINED"
  342 + select["rad"] = "UNDEFINED"
  343 + select["ra"] = "UNDEFINED"
  344 + select["dec"] = "UNDEFINED"
  345 + select.run()
  346 +
  347 + # Retrieve observation
  348 + obs = select.obs()
343 349  
344 350 # Binned analysis
345   - if self.m_binned:
  351 + else:
346 352  
347 353 # Header
348 354 if self.logTerse():
349 355 self.log.header3("Binning events")
350 356  
351 357 # Bin events
352   - bin = ctools.ctbin(select.obs())
  358 + bin = ctools.ctbin(self.obs)
353 359 bin["usepnt"] = False
354 360 bin["ebinalg"] = "LOG"
355 361 bin["xref"] = self.m_xref
... ... @@ -369,7 +375,7 @@ class csspec(ctools.cscript):
369 375 self.log.header3("Creating exposure cube")
370 376  
371 377 # Create exposure cube
372   - expcube = ctools.ctexpcube(select.obs())
  378 + expcube = ctools.ctexpcube(self.obs)
373 379 expcube["incube"] = "NONE"
374 380 expcube["usepnt"] = False
375 381 expcube["ebinalg"] = "LOG"
... ... @@ -383,7 +389,7 @@ class csspec(ctools.cscript):
383 389 expcube["emin"] = emin.TeV()
384 390 expcube["emax"] = emax.TeV()
385 391 expcube["coordsys"] = self.m_coordsys
386   - expcube["proj"] = self.m_proj
  392 + expcube["proj"] = self.m_proj
387 393 expcube.run()
388 394  
389 395 # Header
... ... @@ -391,21 +397,22 @@ class csspec(ctools.cscript):
391 397 self.log.header3("Creating PSF cube")
392 398  
393 399 # Create psf cube
394   - psfcube = ctools.ctpsfcube(select.obs())
  400 + psfcube = ctools.ctpsfcube(self.obs)
395 401 psfcube["incube"] = "NONE"
396 402 psfcube["usepnt"] = False
397 403 psfcube["ebinalg"] = "LOG"
398 404 psfcube["edisp"] = self.m_edisp
399 405 psfcube["xref"] = self.m_xref
400 406 psfcube["yref"] = self.m_yref
401   - psfcube["binsz"] = self.m_binsz
402   - psfcube["nxpix"] = self.m_nxpix
403   - psfcube["nypix"] = self.m_nypix
  407 + psfcube["binsz"] = self.m_binsz*20.0 # Slow Psf variation
  408 + psfcube["nxpix"] = self.m_nxpix/20
  409 + psfcube["nypix"] = self.m_nypix/20
404 410 psfcube["enumbins"] = self.m_ebins
  411 + psfcube["anumbins"] = self.m_anumbins
405 412 psfcube["emin"] = emin.TeV()
406 413 psfcube["emax"] = emax.TeV()
407 414 psfcube["coordsys"] = self.m_coordsys
408   - psfcube["proj"] = self.m_proj
  415 + psfcube["proj"] = self.m_proj
409 416 psfcube.run()
410 417  
411 418 # Header
... ... @@ -413,7 +420,7 @@ class csspec(ctools.cscript):
413 420 self.log.header3("Creating background cube")
414 421  
415 422 # Create background cube
416   - bkgcube = ctools.ctbkgcube(select.obs())
  423 + bkgcube = ctools.ctbkgcube(self.obs)
417 424 bkgcube["incube"] = "NONE"
418 425 bkgcube["usepnt"] = False
419 426 bkgcube["ebinalg"] = "LOG"
... ... @@ -426,7 +433,7 @@ class csspec(ctools.cscript):
426 433 bkgcube["emin"] = emin.TeV()
427 434 bkgcube["emax"] = emax.TeV()
428 435 bkgcube["coordsys"] = self.m_coordsys
429   - bkgcube["proj"] = self.m_proj
  436 + bkgcube["proj"] = self.m_proj
430 437 bkgcube.run()
431 438  
432 439 # Set new binned observation
... ... @@ -462,7 +469,7 @@ class csspec(ctools.cscript):
462 469  
463 470 # Log information
464 471 if self.logTerse():
465   - self.log("No event in this bin. Bin is skipped\n")
  472 + self.log("No event in this bin. Likelihood is zero. Bin is skipped\n")
466 473  
467 474 # Set all values to 0
468 475 flux[i] = 0.0
... ... @@ -505,8 +512,9 @@ class csspec(ctools.cscript):
505 512  
506 513 # Compute Npred value
507 514 Npred = 0.0
508   - for observation in like.obs():
509   - Npred += observation.npred(source)
  515 + if not self.m_binned:
  516 + for observation in like.obs():
  517 + Npred += observation.npred(source)
510 518  
511 519 # Get differential flux
512 520 fitted_flux = source.spectral().eval(elogmean,gammalib.GTime())
... ...
doc/source/reference_manual/csspec.rst
  View file @3b0fd74
... ... @@ -81,6 +81,9 @@ General parameters
81 81  
82 82 ``binsz [real]``
83 83 Pixel size (in degrees/pixel).
  84 +
  85 +``(anumbins = 200) [integer]``
  86 + Number of angular separation bins.
84 87  
85 88 ``(calc_ts = yes) [boolean]``
86 89 Compute TS for each spectral point?
... ...
src/ctbin/ctbin.cpp
  View file @3b0fd74
... ... @@ -485,6 +485,9 @@ void ctbin::fill_cube(GCTAObservation* obs)
485 485 // Get the RoI
486 486 const GCTARoi& roi = events->roi();
487 487  
  488 + // Get the ebounds
  489 + const GEbounds& obs_ebounds = events->ebounds();
  490 +
488 491 // Check for RoI sanity
489 492 if (!roi.is_valid()) {
490 493 std::string msg = "No RoI information found in input observation "
... ... @@ -526,7 +529,9 @@ void ctbin::fill_cube(GCTAObservation* obs)
526 529  
527 530 // Determine energy bin. Skip if we are outside the energy range
528 531 int index = m_ebounds.index(event->energy());
529   - if (index == -1) {
  532 + if (index == -1 ||
  533 + !obs_ebounds.contains(m_ebounds.emin(index),
  534 + m_ebounds.emax(index))) {
530 535 num_outside_ebds++;
531 536 continue;
532 537 }
... ...
src/ctmodel/ctmodel.cpp
  View file @3b0fd74
... ... @@ -589,6 +589,10 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
589 589 const GEbounds& ebounds = obs->events()->ebounds();
590 590 const GGti& gti = obs->events()->gti();
591 591  
  592 + // Get the ebounds
  593 + const GEbounds& obs_ebounds = obs->ebounds();
  594 + const GEbounds& cube_ebounds = m_cube.ebounds();
  595 +
592 596 // Initialise statistics
593 597 double sum = 0.0;
594 598 int num_outside_ebds = 0;
... ... @@ -603,7 +607,10 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
603 607 GCTAEventBin* bin = m_cube[i];
604 608  
605 609 // Skip bin if it is outside the energy range of the observation
606   - if (!ebounds.contains(bin->energy())) {
  610 + int index = cube_ebounds.index(bin->energy());
  611 + if (index == -1 ||
  612 + !obs_ebounds.contains(cube_ebounds.emin(index),
  613 + cube_ebounds.emax(index))) {
607 614 num_outside_ebds++;
608 615 continue;
609 616 }
... ... @@ -613,7 +620,6 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
613 620  
614 621 // Compute model value for cube bin
615 622 double model = m_obs.models().eval(*bin, *obs) * bin->size();
616   - //double model = m_obs.models().eval(*bin, *ptr) * bin->size();
617 623  
618 624 // Add model to actual value
619 625 value += model;
... ...
src/ctselect/ctselect.cpp
  View file @3b0fd74
... ... @@ -631,7 +631,13 @@ void ctselect::select_events(GCTAObservation* obs, const std::string& filename)
631 631 GCTAEventList* list =
632 632 static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
633 633  
634   - // Get energy boundaries
  634 + // Get existing Roi and energy bounds for possible later use
  635 + // (will be empty if unavailable)
  636 + GCTARoi old_roi = list->roi();
  637 + GEbounds old_ebounds = list->ebounds();
  638 +
  639 + // Determine new energy boundaries for selection
  640 + // taking into account previous existing ones
635 641 double emin = 0.0;
636 642 double emax = 0.0;
637 643 GEbounds ebounds = set_ebounds(obs, list->ebounds());
... ... @@ -888,6 +894,12 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
888 894 GCTAInstDir instdir;
889 895 instdir.dir().radec_deg(ra, dec);
890 896 list->roi(GCTARoi(instdir, rad));
  897 + } // endif: Roi selection was performed
  898 +
  899 + else if (old_roi.is_valid()) {
  900 + // Restore old Roi information in case no selection was performed
  901 + // and RoI was existing before
  902 + list->roi(old_roi);
891 903 }
892 904  
893 905 // Set event list GTI (in any case as any event list has a GTI)
... ... @@ -898,6 +910,12 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
898 910 GEbounds ebounds;
899 911 ebounds.append(GEnergy(emin, "TeV"), GEnergy(emax, "TeV"));
900 912 list->ebounds(ebounds);
  913 + } //endif: energy selection was performed
  914 +
  915 + else if (old_ebounds.size() > 0) {
  916 + // Restore old Ebounds in case no energy selection was performed
  917 + // and observation already had valid Ebounds
  918 + list->ebounds(old_ebounds);
901 919 }
902 920  
903 921 // Recompute ontime and livetime.
... ...