Commit d0732e80658ef853ba65a2f8a9cd69f8629b581c
1 parent
a660dbce
Removed spectrum() function form obsutils.
Showing
1 changed file
with
0 additions
and
87 deletions
pyext/ctools/obsutils.py
| ... | ... | @@ -571,90 +571,3 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): |
| 571 | 571 | |
| 572 | 572 | # Return observation definition list |
| 573 | 573 | return obsdeflist |
| 574 | - | |
| 575 | - | |
| 576 | -# ================================================================= # | |
| 577 | -# Generate a spectrum by fitting an observation in energy intervals # | |
| 578 | -# ================================================================= # | |
| 579 | -def spectrum(obs, source, ebounds): | |
| 580 | - """ | |
| 581 | - Returns spectral points in the given binning for the passed observations | |
| 582 | - and source. | |
| 583 | - | |
| 584 | - Please note that specpoints can make spectra for any given spectral | |
| 585 | - model, as long as it has a "Prefactor" or "Integral" parameter which | |
| 586 | - scales the function | |
| 587 | - | |
| 588 | - Parameters: | |
| 589 | - obs - Observation container | |
| 590 | - source - Source name for which spectrum should be computed | |
| 591 | - ebounds - GEbounds energy boundaries | |
| 592 | - Keywords: | |
| 593 | - None | |
| 594 | - """ | |
| 595 | - # Check if parameters Prefactor or Integral exist to scale mdoel | |
| 596 | - normpar = "Prefactor" | |
| 597 | - srcmodel = obs.models()[source] | |
| 598 | - if not srcmodel.has_par(normpar): | |
| 599 | - if srcmodel.has_par("Integral"): | |
| 600 | - normpar = "Integral" | |
| 601 | - else: | |
| 602 | - print("Error: your spectral model of "+source+"does not"+\ | |
| 603 | - "have a scaling parameter named 'Prefactor'"+\ | |
| 604 | - " or 'Integral'") | |
| 605 | - return | |
| 606 | - | |
| 607 | - # Fix all parameters except prefactor | |
| 608 | - opt_spectral = srcmodel.spectral() | |
| 609 | - for model in obs.models(): | |
| 610 | - for par in range(model.size()): | |
| 611 | - model[par].fix() | |
| 612 | - opt_spectral[normpar].free() | |
| 613 | - | |
| 614 | - # Define spectrum dictionary, units are in format readable by astropy | |
| 615 | - spec = {"energy":{"value":[],"eu_value":[],"ed_value":[],"unit":"none"}, | |
| 616 | - "flux" :{"value":[],"eu_value":[],"ed_value":[],"unit":"none"}, | |
| 617 | - "TS" :{"value":[],"eu_value":[],"ed_value":[],"unit":"none"}} | |
| 618 | - | |
| 619 | - # Loop over energy bins, fit and write result to spectrum | |
| 620 | - for i in range(len(ebounds)): | |
| 621 | - | |
| 622 | - # Get logarithmic bin center and flux value of all energy fit | |
| 623 | - ecenter = ebounds.elogmean(i) | |
| 624 | - allflux = opt_spectral.eval(ecenter, gammalib.GTime(0.0)) | |
| 625 | - | |
| 626 | - # Clone observations and reset energy thresholds | |
| 627 | - select = ctools.ctselect(obs) | |
| 628 | - select["ra"] = "UNDEF" | |
| 629 | - select["dec"] = "UNDEF" | |
| 630 | - select["rad"] = "UNDEF" | |
| 631 | - select["tmin"] = 0 | |
| 632 | - select["tmax"] = 0 | |
| 633 | - select["emin"] = ebounds.emin(i).TeV() | |
| 634 | - select["emax"] = ebounds.emax(i).TeV() | |
| 635 | - select.run() | |
| 636 | - | |
| 637 | - # Refit and write out to spectrum | |
| 638 | - binlike = fit(select.obs(), log=False, debug=False, edisp=False) | |
| 639 | - binopt_spectral = binlike.obs().models()[source].spectral() | |
| 640 | - binflux = binopt_spectral.eval(ecenter, gammalib.GTime(0.0)) | |
| 641 | - if binopt_spectral[normpar].value() > 0.0: | |
| 642 | - relerr = binopt_spectral[normpar].error()/binopt_spectral[normpar].value() | |
| 643 | - else: | |
| 644 | - relerr = 0.0 | |
| 645 | - | |
| 646 | - # Store and convert to TeV from MeV | |
| 647 | - spec["energy"]["value"].append(ecenter.TeV()) | |
| 648 | - spec["energy"]["ed_value"].append(ecenter.TeV()-ebounds.emin(i).TeV()) | |
| 649 | - spec["energy"]["eu_value"].append(ebounds.emax(i).TeV()-ecenter.TeV()) | |
| 650 | - spec["energy"]["unit"] = "TeV" | |
| 651 | - | |
| 652 | - spec["flux"]["value"].append(binflux*1.0e6) | |
| 653 | - spec["flux"]["ed_value"].append(relerr*binflux*1.0e6) | |
| 654 | - spec["flux"]["eu_value"].append(relerr*binflux*1.0e6) | |
| 655 | - spec["flux"]["unit"] = "ph/cm2/s/TeV" | |
| 656 | - | |
| 657 | - spec["TS"]["value"].append(binlike.obs().models()[source].ts()) | |
| 658 | - | |
| 659 | - # Return spectrum | |
| 660 | - return spec | ... | ... |