Andreas Specovius / ctools

Commit 2549d836386ffe608eb0454915fe3ac336a3ce85

Authored by Luigi Tibaldo
1 parent a9a8d129

handling exception due to empty observations containers returned from On/Off ana… 

…lysis when using on-axis observations
Inline Side-by-side
Showing 2 changed files with 105 additions and 69 deletions
cscripts/cslightcrv.py
  View file @2549d83
... ... @@ -488,21 +488,82 @@ class cslightcrv(ctools.csobservation):
488 488 self._log_header3(gammalib.EXPLICIT, 'Fitting the data')
489 489  
490 490 # Do maximum likelihood model fitting
491   - like = ctools.ctlike(obs)
492   - like['edisp'] = self['edisp'].boolean()
493   - like.run()
  491 + if obs.size() > 0:
  492 + like = ctools.ctlike(obs)
  493 + like['edisp'] = self['edisp'].boolean()
  494 + like.run()
494 495  
495   - # Skip bin if no event was present
496   - if like.obs().logL() == 0.0:
  496 + # Skip bin if no event was present
  497 + if like.obs().logL() == 0.0:
497 498  
498   - # Signal skipping of bin
  499 + # Signal skipping of bin
  500 + self._log_value(gammalib.TERSE, 'Warning',
  501 + 'No event in this time bin, skip bin.')
  502 +
  503 + # Set all results to 0
  504 + for par in pars:
  505 + result['values'][par] = 0.0
  506 + result['values']['e_'+par] = 0.0
  507 +
  508 + # Append result
  509 + results.append(result)
  510 +
  511 + # Continue with next time bin
  512 + continue
  513 +
  514 + # Retrieve model fitting results for source of interest
  515 + source = like.obs().models()[self._srcname]
  516 +
  517 + # Extract parameter values
  518 + for par in pars:
  519 + result['values'][par] = source[par].value()
  520 + result['values']['e_'+par] = source[par].error()
  521 +
  522 + # Calculate upper limit (-1 if not computed)
  523 + #ul_diff, ul_flux, ul_eflux = self._compute_ulimit(like.obs())
  524 + ul_diff, ul_flux, ul_eflux = self._compute_ulimit(obs)
  525 + if ul_diff > 0.0:
  526 + result['ul_diff'] = ul_diff
  527 + result['ul_flux'] = ul_flux
  528 + result['ul_eflux'] = ul_eflux
  529 +
  530 + # Extract Test Statistic value
  531 + if self['calc_ts'].boolean():
  532 + result['ts'] = source.ts()
  533 +
  534 + # Append result to list of dictionaries
  535 + results.append(result)
  536 +
  537 + # Log results for this time bin
  538 + self._log.header3('Results')
  539 + pars = self._get_free_par_names()
  540 + for par in pars:
  541 + value = source[par].value()
  542 + error = source[par].error()
  543 + unit = source[par].unit()
  544 + self._log_value(gammalib.NORMAL, par,
  545 + str(value)+' +/- '+str(error)+' '+unit)
  546 + if ul_diff > 0.0:
  547 + self._log_value(gammalib.NORMAL, 'Upper flux limit',
  548 + str(result['ul_diff'])+' ph/cm2/s/MeV')
  549 + self._log_value(gammalib.NORMAL, 'Upper flux limit',
  550 + str(result['ul_flux'])+' ph/cm2/s')
  551 + self._log_value(gammalib.NORMAL, 'Upper flux limit',
  552 + str(result['ul_eflux'])+' erg/cm2/s')
  553 + if self['calc_ts'].boolean():
  554 + self._log_value(gammalib.NORMAL, 'Test Statistic', result['ts'])
  555 +
  556 + # Otherwise, if observations size is 0, signal bin is skipped and
  557 + # fill results table with zeros
  558 + else:
499 559 self._log_value(gammalib.TERSE, 'Warning',
500   - 'No event in this time bin, skip bin.')
  560 + 'No observations available in this time bin',
  561 + 'Skip bin.')
501 562  
502 563 # Set all results to 0
503 564 for par in pars:
504   - result['values'][par] = 0.0
505   - result['values']['e_'+par] = 0.0
  565 + result['values'][par] = 0.0
  566 + result['values']['e_' + par] = 0.0
506 567  
507 568 # Append result
508 569 results.append(result)
... ... @@ -510,48 +571,6 @@ class cslightcrv(ctools.csobservation):
510 571 # Continue with next time bin
511 572 continue
512 573  
513   - # Retrieve model fitting results for source of interest
514   - source = like.obs().models()[self._srcname]
515   -
516   - # Extract parameter values
517   - for par in pars:
518   - result['values'][par] = source[par].value()
519   - result['values']['e_'+par] = source[par].error()
520   -
521   - # Calculate upper limit (-1 if not computed)
522   - #ul_diff, ul_flux, ul_eflux = self._compute_ulimit(like.obs())
523   - ul_diff, ul_flux, ul_eflux = self._compute_ulimit(obs)
524   - if ul_diff > 0.0:
525   - result['ul_diff'] = ul_diff
526   - result['ul_flux'] = ul_flux
527   - result['ul_eflux'] = ul_eflux
528   -
529   - # Extract Test Statistic value
530   - if self['calc_ts'].boolean():
531   - result['ts'] = source.ts()
532   -
533   - # Append result to list of dictionaries
534   - results.append(result)
535   -
536   - # Log results for this time bin
537   - self._log.header3('Results')
538   - pars = self._get_free_par_names()
539   - for par in pars:
540   - value = source[par].value()
541   - error = source[par].error()
542   - unit = source[par].unit()
543   - self._log_value(gammalib.NORMAL, par,
544   - str(value)+' +/- '+str(error)+' '+unit)
545   - if ul_diff > 0.0:
546   - self._log_value(gammalib.NORMAL, 'Upper flux limit',
547   - str(result['ul_diff'])+' ph/cm2/s/MeV')
548   - self._log_value(gammalib.NORMAL, 'Upper flux limit',
549   - str(result['ul_flux'])+' ph/cm2/s')
550   - self._log_value(gammalib.NORMAL, 'Upper flux limit',
551   - str(result['ul_eflux'])+' erg/cm2/s')
552   - if self['calc_ts'].boolean():
553   - self._log_value(gammalib.NORMAL, 'Test Statistic', result['ts'])
554   -
555 574 # Create FITS table from results
556 575 table = self._create_fits_table(results)
557 576  
... ...
cscripts/csphasecrv.py
  View file @2549d83
... ... @@ -242,11 +242,17 @@ class csphasecrv(ctools.csobservation):
242 242  
243 243 # Store fit results
244 244 phname = str(phmin)+'-'+str(phmax)
245   - source = self._fitmodels[phname][self._srcname]
246   - for par in pars:
247   - result['values'][par] = source[par].value()
248   - result['values']['e_'+par] = source[par].error()
249   -
  245 + # If the model contains the source of interest fill results
  246 + try:
  247 + source = self._fitmodels[phname][self._srcname]
  248 + for par in pars:
  249 + result['values'][par] = source[par].value()
  250 + result['values']['e_'+par] = source[par].error()
  251 + # Otherwise fills with zeros
  252 + except:
  253 + for par in pars:
  254 + result['values'][par] = 0.
  255 + result['values']['e_'+par] = 0.
250 256 # Append result to list of dictionaries
251 257 results.append(result)
252 258  
... ... @@ -374,19 +380,30 @@ class csphasecrv(ctools.csobservation):
374 380 self._log_header3(gammalib.EXPLICIT, 'Fitting the data')
375 381  
376 382 # Do maximum likelihood model fitting
377   - like = ctools.ctlike(obs)
378   - like['edisp'] = self['edisp'].boolean()
379   - like.run()
380   -
381   - # Renormalize models to phase selection
382   - # TODO move the scaling from the temporal to the spectral component
383   - for model in like.obs().models():
384   - scaled_norm = model['Normalization'].value()/(phbin[1]-phbin[0])
385   - model['Normalization'].value(scaled_norm)
386   -
387   - # Store fit model
388   - self._fitmodels[str(phbin[0])+'-'+str(phbin[1])] = \
389   - like.obs().models().copy()
  383 + # The On/Off analysis can produce empty observation containers,
  384 + # e.g., when on-axis observations are used. To avoid ctlike asking
  385 + # for a new observations contained (or hang, in on interactive mode)
  386 + # we'll run ctlike only if the size is >0
  387 + if obs.size() > 0:
  388 + like = ctools.ctlike(obs)
  389 + like['edisp'] = self['edisp'].boolean()
  390 + like.run()
  391 +
  392 + # Renormalize models to phase selection
  393 + # TODO move the scaling from the temporal to the spectral component
  394 + for model in like.obs().models():
  395 + scaled_norm = model['Normalization'].value()/(phbin[1]-phbin[0])
  396 + model['Normalization'].value(scaled_norm)
  397 + # Store fit model
  398 + self._fitmodels[str(phbin[0])+'-'+str(phbin[1])] = \
  399 + like.obs().models().copy()
  400 + else:
  401 + self._log_string(gammalib.TERSE,
  402 + 'PHASE %f - %f: no observations available'
  403 + ' for fitting' %(phbin[0], phbin[1]))
  404 + # empty models container
  405 + self._fitmodels[str(phbin[0]) + '-' + str(phbin[1])] = \
  406 + gammalib.GModels()
390 407  
391 408 # Create FITS file
392 409 self._create_fits()
... ...