Marianne Lemoine / ctools

Commit 0417bf5ce489bf4301f19c13311383079169cc49

Authored by Jürgen Knödlseder
1 parent c0ab3a39

Add script to create the node section of a NodeFunction spectral model. 

The script cta_make_nodes.py can be used to create the node section of a NodeFunction XML model. It initialises the energy values with a logarithmic binning, and sets initial intensity values using a Crab spectrum based on a power law function.
Inline Side-by-side
Showing 2 changed files with 119 additions and 1 deletions
scripts/Makefile.am
  View file @0417bf5
1 1 ## Process this file with automake to produce Makefile.in
2 2  
3 3 # Scripts to be distributed in bindir
4   -binscripts = cta_kb_bkgrates.py cta_root2irf.py
  4 +binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py
5 5  
6 6 # Scripts to be distributed in pythondir
7 7 pythonscripts = cspull.py cssens.py cstsdist.py obsutils.py
... ...
scripts/cta_make_nodes.py 0 → 100755
  View file @0417bf5
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This Python script creates the node section of the NodeFunction using
  4 +# logarithmically spaced energy bins. The intensity scale is set to the
  5 +# HESS Crab intensity (assuming a power law).
  6 +#
  7 +# Copyright (C) 2012 Juergen Knoedlseder
  8 +#
  9 +# This program is free software: you can redistribute it and/or modify
  10 +# it under the terms of the GNU General Public License as published by
  11 +# the Free Software Foundation, either version 3 of the License, or
  12 +# (at your option) any later version.
  13 +#
  14 +# This program is distributed in the hope that it will be useful,
  15 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  17 +# GNU General Public License for more details.
  18 +#
  19 +# You should have received a copy of the GNU General Public License
  20 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  21 +#
  22 +# ==========================================================================
  23 +from gammalib import *
  24 +import sys
  25 +import math
  26 +
  27 +
  28 +# ============== #
  29 +# Write one node #
  30 +# ============== #
  31 +def write_node(file, energy, scale):
  32 + """
  33 + Writes one node to XML file.
  34 + """
  35 + # Convert to strings
  36 + str_energy = str(energy)
  37 + str_scale = str(scale)
  38 +
  39 + # Write start tag
  40 + file.write(' <node>\n')
  41 +
  42 + # Write energy
  43 + file.write(' <parameter scale="1e6" name="Energy"')
  44 + file.write(' min="'+str_energy+'" max="'+str_energy+'"')
  45 + file.write(' value="'+str_energy+'"')
  46 + file.write(' free="0"/>\n')
  47 +
  48 + # Write intensity
  49 + file.write(' <parameter scale="'+str_scale+'" name="Intensity"')
  50 + file.write(' min="1e-5" max="1e5"')
  51 + file.write(' value="1.0"')
  52 + file.write(' free="1"/>\n')
  53 +
  54 + # Write end tag
  55 + file.write(' </node>\n')
  56 +
  57 + # Return
  58 + return
  59 +
  60 +
  61 +# ============ #
  62 +# Create nodes #
  63 +# ============ #
  64 +def create_nodes(emin, emax, enumbins):
  65 + """
  66 + Create nodes (energies in TeV).
  67 + """
  68 + # Open file
  69 + file = open("nodes.xml", "w")
  70 +
  71 + # Set node boundaries
  72 + elogmin = math.log10(float(emin))
  73 + elogmax = math.log10(float(emax))
  74 + elogbin = (elogmax - elogmin)/(float(enumbins)-1.0)
  75 +
  76 + # Fill arrays
  77 + for i in range(int(enumbins)):
  78 +
  79 + # Compute energy
  80 + energy = math.pow(10.0, i*elogbin+elogmin)
  81 +
  82 + # Compute scale (HESS Crab spectrum)
  83 + scale = 3.45e-17 * math.pow(energy, -2.63)
  84 +
  85 + # Write node
  86 + write_node(file, energy, scale)
  87 +
  88 + # Debug
  89 + #print energy, math.pow(10.0, int(math.log10(scale))-1.0)
  90 +
  91 + # Close file
  92 + file.close()
  93 +
  94 + # Return
  95 + return
  96 +
  97 +
  98 +# ======================= #
  99 +# Main script entry point #
  100 +# ======================= #
  101 +if __name__ == '__main__':
  102 + """
  103 + Compute nodes.
  104 + """
  105 + # Check command line
  106 + usage = "Usage: cta_make_nodes emin emax enumbins"
  107 + if len(sys.argv) < 3:
  108 + print usage
  109 + sys.exit()
  110 +
  111 + # Extract parameters
  112 + emin = sys.argv[1]
  113 + emax = sys.argv[2]
  114 + enumbins = sys.argv[3]
  115 +
  116 + # Create nodes
  117 + create_nodes(emin, emax, enumbins)
  118 +
... ...