From 0417bf5ce489bf4301f19c13311383079169cc49 Mon Sep 17 00:00:00 2001
From: Jürgen Knödlseder <jurgen.knodlseder@irap.omp.eu>
Date: Fri, 20 Jul 2012 23:23:39 +0200
Subject: [PATCH] Add script to create the node section of a NodeFunction spectral model.

---
 scripts/Makefile.am       |   2 +-
 scripts/cta_make_nodes.py | 118 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 119 insertions(+), 1 deletion(-)
 create mode 100755 scripts/cta_make_nodes.py

diff --git a/scripts/Makefile.am b/scripts/Makefile.am
index 93207a3..7360197 100644
--- a/scripts/Makefile.am
+++ b/scripts/Makefile.am
@@ -1,7 +1,7 @@
 ## Process this file with automake to produce Makefile.in
 
 # Scripts to be distributed in bindir
-binscripts = cta_kb_bkgrates.py cta_root2irf.py
+binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py
 
 # Scripts to be distributed in pythondir
 pythonscripts = cspull.py cssens.py cstsdist.py obsutils.py
diff --git a/scripts/cta_make_nodes.py b/scripts/cta_make_nodes.py
new file mode 100755
index 0000000..4fa9384
--- /dev/null
+++ b/scripts/cta_make_nodes.py
@@ -0,0 +1,118 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This Python script creates the node section of the NodeFunction using
+# logarithmically spaced energy bins. The intensity scale is set to the
+# HESS Crab intensity (assuming a power law).
+#
+# Copyright (C) 2012 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+from gammalib import *
+import sys
+import math
+
+
+# ============== #
+# Write one node #
+# ============== #
+def write_node(file, energy, scale):
+	"""
+	Writes one node to XML file.
+	"""
+	# Convert to strings
+	str_energy = str(energy)
+	str_scale  = str(scale)
+
+	# Write start tag
+	file.write('      <node>\n')
+
+	# Write energy
+	file.write('        <parameter scale="1e6"   name="Energy"')
+	file.write('    min="'+str_energy+'"   max="'+str_energy+'"')
+	file.write(' value="'+str_energy+'"')
+	file.write(' free="0"/>\n')
+
+	# Write intensity
+	file.write('        <parameter scale="'+str_scale+'"   name="Intensity"')
+	file.write(' min="1e-5"   max="1e5"')
+	file.write(' value="1.0"')
+	file.write(' free="1"/>\n')
+
+	# Write end tag
+	file.write('      </node>\n')
+
+	# Return
+	return
+	
+
+# ============ #
+# Create nodes #
+# ============ #
+def create_nodes(emin, emax, enumbins):
+	"""
+	Create nodes (energies in TeV).
+	"""
+	# Open file
+	file = open("nodes.xml", "w")
+
+	# Set node boundaries
+	elogmin = math.log10(float(emin))
+	elogmax = math.log10(float(emax))
+	elogbin = (elogmax - elogmin)/(float(enumbins)-1.0)
+	
+	# Fill arrays
+	for i in range(int(enumbins)):
+		
+		# Compute energy
+		energy = math.pow(10.0, i*elogbin+elogmin)
+		
+		# Compute scale (HESS Crab spectrum)
+		scale = 3.45e-17 * math.pow(energy, -2.63)
+		
+		# Write node
+		write_node(file, energy, scale)
+		
+		# Debug
+		#print energy, math.pow(10.0, int(math.log10(scale))-1.0)
+	
+	# Close file
+	file.close()
+	
+	# Return
+	return
+	
+
+# ======================= #
+# Main script entry point #
+# ======================= #
+if __name__ == '__main__':
+	"""
+	Compute nodes.
+	"""
+	# Check command line
+	usage = "Usage: cta_make_nodes emin emax enumbins"
+	if len(sys.argv) < 3:
+		print usage
+		sys.exit()
+	
+	# Extract parameters
+	emin     = sys.argv[1]
+	emax     = sys.argv[2]
+	enumbins = sys.argv[3]
+
+	# Create nodes
+	create_nodes(emin, emax, enumbins)
+	
--
libgit2 0.21.4