Commit eede74f2af6373dcff6622ed125481c054205bf4
1 parent
d6f6ec21
Remove new code before it gets corrected
Showing
1 changed file
with
0 additions
and
94 deletions
pyext/ctools/obsutils.py
| @@ -516,97 +516,3 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): | @@ -516,97 +516,3 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): | ||
| 516 | 516 | ||
| 517 | # Return observation definition list | 517 | # Return observation definition list |
| 518 | return obsdeflist | 518 | return obsdeflist |
| 519 | - | ||
| 520 | - | ||
| 521 | -# ================================================= # | ||
| 522 | -# Make spectrum for one Source in a CTA Observation # | ||
| 523 | -# ================================================= # | ||
| 524 | -def getlogbins(nbins=12,emin=0.1,emax=100): | ||
| 525 | - """ | ||
| 526 | - Returns list with logarithmiscenergy bins in | ||
| 527 | - format [[emin1,emax1],[emin2,emax2],..] | ||
| 528 | - | ||
| 529 | - Parameters: | ||
| 530 | - nbins Number of bins | ||
| 531 | - emin Minimum energy | ||
| 532 | - emax Maximum energy | ||
| 533 | - """ | ||
| 534 | - lgemin=numpy.log10(emin) | ||
| 535 | - lgemax=numpy.log10(emax) | ||
| 536 | - lgbins=[] | ||
| 537 | - lgedges=numpy.arange(lgemin,lgemax,(lgemax-lgemin)/nbins) | ||
| 538 | - lgedges=numpy.append(lgedges,[lgemax]) | ||
| 539 | - for ii in range(len(lgedges)-1): | ||
| 540 | - lgbins.append([lgedges[ii],lgedges[ii+1]]) | ||
| 541 | - ebins=numpy.power(10,lgbins) | ||
| 542 | - return ebins | ||
| 543 | - | ||
| 544 | -#Please note that specpoints can make spectra for any given | ||
| 545 | -#spectral model, as long as it has a "Prefactor" parameter | ||
| 546 | -#which scales the funtion | ||
| 547 | -def specpoints(obs,source,ebins): | ||
| 548 | - """ | ||
| 549 | - Returns spectral points in the given binning for | ||
| 550 | - the passed observations and source. | ||
| 551 | - | ||
| 552 | - Parameters: | ||
| 553 | - obs - Pointing direction | ||
| 554 | - source - Source name for which spectrum should be done | ||
| 555 | - ebins - Binning passed as a list in units of TeV as | ||
| 556 | - [[emin1,emax1],[emin2,emax2],..] | ||
| 557 | - Keywords: | ||
| 558 | - """ | ||
| 559 | - #Fix all parameters except prefactor | ||
| 560 | - for model in obs.models(): | ||
| 561 | - for par in range(model.size()): | ||
| 562 | - model[par].fix() | ||
| 563 | - | ||
| 564 | - opt_spectral = obs.models()[source].spectral() | ||
| 565 | - opt_spectral["Prefactor"].free() | ||
| 566 | - | ||
| 567 | - #Define spectrum dictionary, units are in format readable by astropy | ||
| 568 | - pars = {"value":[],"eu_value":[],"ed_value":[],"unit":[]} | ||
| 569 | - spec = {"energy":pars.copy(),"flux":pars.copy(),"TS":[]} | ||
| 570 | - | ||
| 571 | - #Loop over energy bins, fit and write result to spectrum | ||
| 572 | - for ebin in ebins: | ||
| 573 | - #Get logarithmic bin center and flux value of all energy fit | ||
| 574 | - ecenter = numpy.power(10,(numpy.log10(ebin[1])+numpy.log10(ebin[0]))/2.) | ||
| 575 | - allflux = opt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) | ||
| 576 | - | ||
| 577 | - #Clone observations and reset energy thresholds | ||
| 578 | - select = ctools.ctselect(obs) | ||
| 579 | - select["ra"].real(-1.0) | ||
| 580 | - select["dec"].real(-1.0) | ||
| 581 | - select["rad"].real(-1.0) | ||
| 582 | - select["tmin"].real(0) | ||
| 583 | - select["tmax"].real(0) | ||
| 584 | - select["emin"].real(ebin[0]) | ||
| 585 | - select["emax"].real(ebin[1]) | ||
| 586 | - select.run() | ||
| 587 | - obsclone = select.obs() | ||
| 588 | - | ||
| 589 | - #Refit and write out to spectrum | ||
| 590 | - binlike=fit(obsclone,log=False, debug=False, edisp=False) | ||
| 591 | - binopt_spectral = binlike.obs().models()[source].spectral() | ||
| 592 | - binflux=binopt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0)) | ||
| 593 | - if binopt_spectral["Prefactor"].value()>0: | ||
| 594 | - relerr=binopt_spectral["Prefactor"].error()/binopt_spectral["Prefactor"].value() | ||
| 595 | - else: | ||
| 596 | - relerr=0 | ||
| 597 | - | ||
| 598 | - #Store and convert to TeV from MeV | ||
| 599 | - spec["energy"]["value"].append(ecenter) | ||
| 600 | - spec["energy"]["ed_value"].append(ecenter-ebin[0]) | ||
| 601 | - spec["energy"]["eu_value"].append(ebin[1]-ecenter) | ||
| 602 | - spec["energy"]["unit"].append("1e12 eV") | ||
| 603 | - | ||
| 604 | - spec["flux"]["value"].append(binflux*1e6) | ||
| 605 | - spec["flux"]["ed_value"].append(relerr*binflux*1e6) | ||
| 606 | - spec["flux"]["eu_value"].append(relerr*binflux*1e6) | ||
| 607 | - spec["flux"]["unit"].append("cm-2 1e-12 eV-1") | ||
| 608 | - | ||
| 609 | - spec["TS"].append(binlike.obs().models()[source].ts()) | ||
| 610 | - | ||
| 611 | - #Return spectrum | ||
| 612 | - return spec |