Veronika Vodeb / ctools

Commit 337b228b2362edceadfabfe37c36710d33af4ce4

Authored by Jürgen Knödlseder
1 parent 420ad0e2

Add spectrum() function to obsutils and create make_spectrum.py example script t… 

…o demonstrate its usage.
Inline Side-by-side
Showing 2 changed files with 170 additions and 68 deletions
examples/make_spectrum.py 0 → 100755
  View file @337b228
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script shows how to make a spectrum using obsutils.
  4 +#
  5 +# Copyright (C) 2014 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import ctools
  23 +try:
  24 + import matplotlib.pyplot as plt
  25 + has_matplotlib = True
  26 +except:
  27 + has_matplotlib = False
  28 +
  29 +
  30 +# =============== #
  31 +# Make a spectrum #
  32 +# =============== #
  33 +def make_spectrum():
  34 + """
  35 + """
  36 + # Set script parameters
  37 + model_name = "${CTOOLS}/share/models/crab.xml"
  38 + caldb = "dummy"
  39 + irf = "cta_dummy_irf"
  40 + ra = 83.63
  41 + dec = 22.01
  42 + rad_sim = 3.0
  43 + tstart = 0.0
  44 + tstop = 1800.0
  45 + emin = 0.1
  46 + emax = 100.0
  47 +
  48 + # Simulate events
  49 + sim = ctools.ctobssim()
  50 + sim["infile"].filename(model_name)
  51 + sim["caldb"].string(caldb)
  52 + sim["irf"].string(irf)
  53 + sim["ra"].real(ra)
  54 + sim["dec"].real(dec)
  55 + sim["rad"].real(rad_sim)
  56 + sim["tmin"].real(tstart)
  57 + sim["tmax"].real(tstop)
  58 + sim["emin"].real(emin)
  59 + sim["emax"].real(emax)
  60 + sim.run()
  61 +
  62 + # Generate an energy binning
  63 + e_min = gammalib.GEnergy(emin, "TeV")
  64 + e_max = gammalib.GEnergy(emax, "TeV")
  65 + ebounds = gammalib.GEbounds(10, e_min, e_max)
  66 +
  67 + # Generate spectral points
  68 + spectrum = ctools.obsutils.spectrum(sim.obs(), "Crab", ebounds)
  69 +
  70 + # Return spectrum
  71 + return spectrum
  72 +
  73 +
  74 +# ============= #
  75 +# Plot spectrum #
  76 +# ============= #
  77 +def plot_spectrum(spectrum):
  78 + """
  79 + Plot spectrum.
  80 + """
  81 + # Create figure
  82 + plt.figure(1)
  83 + plt.title("Crab spectrum")
  84 +
  85 + # Plot spectrum
  86 + plt.loglog(spectrum['energy']['value'], \
  87 + spectrum['flux']['value'], 'ro', label='Crab')
  88 + plt.errorbar(spectrum['energy']['value'], \
  89 + spectrum['flux']['value'], \
  90 + spectrum['flux']['ed_value'], fmt=None, ecolor='r')
  91 +
  92 + # Put labels
  93 + plt.xlabel("Energy ("+spectrum['energy']['unit']+")")
  94 + plt.ylabel("Flux ("+spectrum['flux']['unit']+")")
  95 + plt.legend(loc="lower left")
  96 +
  97 + # Show spectrum
  98 + plt.show()
  99 +
  100 + # Return
  101 + return
  102 +
  103 +
  104 +#==========================#
  105 +# Main routine entry point #
  106 +#==========================#
  107 +if __name__ == '__main__':
  108 + """
  109 + Generate a Crab spectrum.
  110 + """
  111 + # Generate spectrum
  112 + spectrum = make_spectrum()
  113 +
  114 + # Plot spectrum
  115 + if has_matplotlib:
  116 + plot_spectrum(spectrum)
  117 +
0 118 \ No newline at end of file
... ...
pyext/ctools/obsutils.py
  View file @337b228
... ... @@ -522,109 +522,94 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
522 522  
523 523 # Return observation definition list
524 524 return obsdeflist
525   -
526   -# ================================================= #
527   -# Make spectrum for one Source in a CTA Observation #
528   -# ================================================= #
529 525  
530   -def getlogbins(nbins=12,emin=0.1,emax=100):
531   - """
532   - Returns list with logarithmiscenergy bins in
533   - format [[emin1,emax1],[emin2,emax2],..]
534   -
535   - Parameters:
536   - nbins Number of bins
537   - emin Minimum energy
538   - emax Maximum energy
539   - """
540   - lgemin=math.log10(emin)
541   - lgemax=math.log10(emax)
542   - lgsize=(lgemax-lgemin)/nbins
543   - ebins=[]
544   - for ii in range(0,nbins):
545   - ebins.append([math.pow(10,lgemin),math.pow(10,lgemin+lgsize)])
546   - lgemin=lgemin+lgsize
547   - return ebins
548   -
549   -# Please note that specpoints can make spectra for any given
550   -# spectral model, as long as it has a "Prefactor" or "Integral"
551   -# parameter which scales the function
552   -def specpoints(obs,source,ebins):
  526 +
  527 +# ================================================================= #
  528 +# Generate a spectrum by fitting an observation in energy intervals #
  529 +# ================================================================= #
  530 +def spectrum(obs, source, ebounds):
553 531 """
554   - Returns spectral points in the given binning for
555   - the passed observations and source.
  532 + Returns spectral points in the given binning for the passed observations
  533 + and source.
  534 +
  535 + Please note that specpoints can make spectra for any given spectral
  536 + model, as long as it has a "Prefactor" or "Integral" parameter which
  537 + scales the function
556 538  
557 539 Parameters:
558   - obs - Pointing direction
559   - source - Source name for which spectrum should be done
560   - ebins - Binning passed as a list in units of TeV as
561   - [[emin1,emax1],[emin2,emax2],..]
  540 + obs - Observation container
  541 + source - Source name for which spectrum should be computed
  542 + ebounds - GEbounds energy boundaries
562 543 Keywords:
  544 + None
563 545 """
564   - #Check if parameters Prefactor or Integral exist to scale mdoel
565   - normpar = "Prefactor"
  546 + # Check if parameters Prefactor or Integral exist to scale mdoel
  547 + normpar = "Prefactor"
566 548 srcmodel = obs.models()[source]
567 549 if not srcmodel.has_par(normpar):
568 550 if srcmodel.has_par("Integral"):
569 551 normpar = "Integral"
570 552 else:
571   - print "Error: your spectral model of ",source,"does not\
572   - have a scaling parameter named 'Prefactor' or 'Integral'"
  553 + print("Error: your spectral model of "+source+"does not"+\
  554 + "have a scaling parameter named 'Prefactor'"+\
  555 + " or 'Integral'")
573 556 return
574 557  
575   - #Fix all parameters except prefactor
  558 + # Fix all parameters except prefactor
576 559 opt_spectral = srcmodel.spectral()
577 560 for model in obs.models():
578 561 for par in range(model.size()):
579 562 model[par].fix()
580 563 opt_spectral[normpar].free()
581 564  
582   - #Define spectrum dictionary, units are in format readable by astropy
  565 + # Define spectrum dictionary, units are in format readable by astropy
583 566 spec = {"energy":{"value":[],"eu_value":[],"ed_value":[],"unit":"none"},
584 567 "flux" :{"value":[],"eu_value":[],"ed_value":[],"unit":"none"},
585 568 "TS" :{"value":[],"eu_value":[],"ed_value":[],"unit":"none"}}
586 569  
587   - #Loop over energy bins, fit and write result to spectrum
588   - for ebin in ebins:
589   - #Get logarithmic bin center and flux value of all energy fit
590   - ecenter = math.pow(10,(math.log10(ebin[1])+math.log10(ebin[0]))/2.)
591   - allflux = opt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0))
  570 + # Loop over energy bins, fit and write result to spectrum
  571 + for i in range(len(ebounds)):
592 572  
593   - #Clone observations and reset energy thresholds
  573 + # Get logarithmic bin center and flux value of all energy fit
  574 + ecenter = ebounds.elogmean(i)
  575 + allflux = opt_spectral.eval(ecenter, gammalib.GTime(0.0))
  576 +
  577 + # Clone observations and reset energy thresholds
594 578 select = ctools.ctselect(obs)
595 579 select["ra"].real(-1.0)
596 580 select["dec"].real(-1.0)
597 581 select["rad"].real(-1.0)
598   - select["tmin"].real(0)
599   - select["tmax"].real(0)
600   - select["emin"].real(ebin[0])
601   - select["emax"].real(ebin[1])
  582 + select["tmin"].real(0.0)
  583 + select["tmax"].real(0.0)
  584 + select["emin"].real(ebounds.emin(i).TeV())
  585 + select["emax"].real(ebounds.emax(i).TeV())
602 586 select.run()
603   - obsclone = select.obs()
604 587  
605   - #Refit and write out to spectrum
606   - binlike=fit(obsclone,log=False, debug=False, edisp=False,tscalc=True)
  588 + # Refit and write out to spectrum
  589 + binlike = fit(select.obs(), log=False, debug=False, edisp=False, tscalc=True)
607 590 binopt_spectral = binlike.obs().models()[source].spectral()
608   - binflux=binopt_spectral.eval(gammalib.GEnergy(ecenter,"TeV"),gammalib.GTime(0))
609   - if binopt_spectral[normpar].value()>0:
610   - relerr=binopt_spectral[normpar].error()/binopt_spectral[normpar].value()
  591 + binflux = binopt_spectral.eval(ecenter, gammalib.GTime(0.0))
  592 + if binopt_spectral[normpar].value() > 0.0:
  593 + relerr = binopt_spectral[normpar].error()/binopt_spectral[normpar].value()
611 594 else:
612   - relerr=0
  595 + relerr = 0.0
613 596  
614   - #Store and convert to TeV from MeV
615   - spec["energy"]["value"].append(ecenter)
616   - spec["energy"]["ed_value"].append(ecenter-ebin[0])
617   - spec["energy"]["eu_value"].append(ebin[1]-ecenter)
618   - spec["energy"]["unit"]="1e12 eV"
  597 + # Store and convert to TeV from MeV
  598 + spec["energy"]["value"].append(ecenter.TeV())
  599 + spec["energy"]["ed_value"].append(ecenter.TeV()-ebounds.emin(i).TeV())
  600 + spec["energy"]["eu_value"].append(ebounds.emax(i).TeV()-ecenter.TeV())
  601 + spec["energy"]["unit"] = "TeV"
619 602  
620   - spec["flux"]["value"].append(binflux*1e6)
621   - spec["flux"]["ed_value"].append(relerr*binflux*1e6)
622   - spec["flux"]["eu_value"].append(relerr*binflux*1e6)
623   - spec["flux"]["unit"]="cm-2 1e-12 eV-1"
  603 + spec["flux"]["value"].append(binflux*1.0e6)
  604 + spec["flux"]["ed_value"].append(relerr*binflux*1.0e6)
  605 + spec["flux"]["eu_value"].append(relerr*binflux*1.0e6)
  606 + spec["flux"]["unit"] = "ph/cm2/s/TeV"
624 607  
625 608 spec["TS"]["value"].append(binlike.obs().models()[source].ts())
626 609  
627   - print spec["energy"]["value"]
628   - print spec["flux"]["value"]
629   - #Return spectrum
  610 + # Debugging???
  611 + #print(spec["energy"]["value"])
  612 + #print(spec["flux"]["value"])
  613 +
  614 + # Return spectrum
630 615 return spec
... ...