Commit e55c60c7e2b977411c50d11c79f1b35eb9f74cd9
1 parent
891c1bfe
changed model query logic
Showing
1 changed file
with
5 additions
and
66 deletions
cscripts/csresspec.py
| @@ -21,7 +21,6 @@ | @@ -21,7 +21,6 @@ | ||
| 21 | import gammalib | 21 | import gammalib |
| 22 | import ctools | 22 | import ctools |
| 23 | from cscripts import obsutils | 23 | from cscripts import obsutils |
| 24 | -import math | ||
| 25 | 24 | ||
| 26 | 25 | ||
| 27 | # =============== # | 26 | # =============== # |
| @@ -130,17 +129,10 @@ class csresspec(ctools.csobservation): | @@ -130,17 +129,10 @@ class csresspec(ctools.csobservation): | ||
| 130 | self['nypix'].real() | 129 | self['nypix'].real() |
| 131 | self['binsz'].real() | 130 | self['binsz'].real() |
| 132 | 131 | ||
| 133 | - # If we are not using a precomputed model query input XML model file | ||
| 134 | - if not self._use_maps: | ||
| 135 | - modxml = self['inmodel'].filename() | ||
| 136 | - # If None check whether models are provided in observation | ||
| 137 | - # container, otherwise throw exception and stop here | ||
| 138 | - if modxml == 'NONE' and self.obs().models().is_empty(): | ||
| 139 | - msg = 'No model provided. Please specify an input XML model file' | ||
| 140 | - raise RuntimeError(msg) | ||
| 141 | - # Otherwise set the input model | ||
| 142 | - else: | ||
| 143 | - self.obs().models(modxml) | 132 | + # If we are not using a precomputed model and no models are available |
| 133 | + # in the observation container query input XML model file | ||
| 134 | + if not self._use_maps and self.obs().models().size() == 0: | ||
| 135 | + self.obs().models(self['inmodel'].filename()) | ||
| 144 | 136 | ||
| 145 | # Unless all observations are On/Off query for mask definition | 137 | # Unless all observations are On/Off query for mask definition |
| 146 | if n_onoff == n_cta: | 138 | if n_onoff == n_cta: |
| @@ -230,7 +222,7 @@ class csresspec(ctools.csobservation): | @@ -230,7 +222,7 @@ class csresspec(ctools.csobservation): | ||
| 230 | cntcube['proj'] = 'TAN' | 222 | cntcube['proj'] = 'TAN' |
| 231 | cntcube['coordsys'] = 'CEL' | 223 | cntcube['coordsys'] = 'CEL' |
| 232 | cntcube['ebinalg'] = self['ebinalg'].string() | 224 | cntcube['ebinalg'] = self['ebinalg'].string() |
| 233 | - if self['ebinalg'].string() == 'FILE' | 225 | + if self['ebinalg'].string() == 'FILE': |
| 234 | cntcube['ebinfile'] = self['ebinfile'].filename() | 226 | cntcube['ebinfile'] = self['ebinfile'].filename() |
| 235 | else: | 227 | else: |
| 236 | cntcube['enumbins'] = self['enumbins'].integer() | 228 | cntcube['enumbins'] = self['enumbins'].integer() |
| @@ -427,59 +419,6 @@ def run(self): | @@ -427,59 +419,6 @@ def run(self): | ||
| 427 | counts = cntcube.counts().total_counts() | 419 | counts = cntcube.counts().total_counts() |
| 428 | model = modcube.counts().total_counts() | 420 | model = modcube.counts().total_counts() |
| 429 | 421 | ||
| 430 | - # Get residual algorithm type | ||
| 431 | - algorithm = self['algorithm'].string() | ||
| 432 | - | ||
| 433 | - # Subtract | ||
| 434 | - if algorithm == 'SUB': | ||
| 435 | - residuals = counts - model | ||
| 436 | - | ||
| 437 | - # Subtract and divide by model map | ||
| 438 | - elif algorithm == 'SUBDIV': | ||
| 439 | - residuals = counts - model | ||
| 440 | - residuals /= model | ||
| 441 | - | ||
| 442 | - # Subtract and divide by sqrt of model map | ||
| 443 | - elif algorithm == 'SUBDIVSQRT': | ||
| 444 | - residuals = counts - model | ||
| 445 | - residuals /= model.sqrt() | ||
| 446 | - | ||
| 447 | - # Calculate significance from Li&Ma derivation | ||
| 448 | - elif algorithm == 'SIGNIFICANCE': | ||
| 449 | - residuals = counts.copy() | ||
| 450 | - | ||
| 451 | - # Compute sign array | ||
| 452 | - sign = (counts - model).sign() ##not implemented!! | ||
| 453 | - | ||
| 454 | - # Loop over energy bins | ||
| 455 | - for i in range(counts.size()): | ||
| 456 | - | ||
| 457 | - # If the model value > 0.0 do the computation as normal ... | ||
| 458 | - model_val = model[i] | ||
| 459 | - if model_val > 0.0: | ||
| 460 | - | ||
| 461 | - # If the data value is also > 0 then compute the | ||
| 462 | - # significance^2 and fill it in the residuals ... | ||
| 463 | - data_val = counts[i] | ||
| 464 | - if data_val > 0.0: | ||
| 465 | - log_val = math.log(data_val / model_val) | ||
| 466 | - residuals[i] = ( | ||
| 467 | - data_val * log_val) + model_val - data_val | ||
| 468 | - | ||
| 469 | - # ... otherwise compute the reduced value of the above | ||
| 470 | - # expression. This is necessary to avoid computing log(0). | ||
| 471 | - else: | ||
| 472 | - residuals[i] = model_val | ||
| 473 | - | ||
| 474 | - # ... otherwise hard-code the significance to 0 | ||
| 475 | - else: | ||
| 476 | - residuals[i] = 0.0 | ||
| 477 | - | ||
| 478 | - # Compute significance | ||
| 479 | - residuals *= 2.0 | ||
| 480 | - residuals = residuals.sqrt() | ||
| 481 | - residuals = residuals * sign | ||
| 482 | - | ||
| 483 | ## Calculate models of individual components if requested | 422 | ## Calculate models of individual components if requested |
| 484 | 423 | ||
| 485 | # otherwise, if On/Off | 424 | # otherwise, if On/Off |