Lili Yang / ctools

Commit 937b71e2f6e9ab41e56148712ca3eab461ef1647

Authored by knodlseder
1 parent 9aea0130

initial version

Inline Side-by-side
Showing 1 changed file with 105 additions and 0 deletions
scripts/cta_kb_bkgrates.py 0 → 100755
  View file @937b71e
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script converts the background rates in Konrad's performance files
  4 +# into a column separated file function table. The table is automatically
  5 +# named "bkg_XXX.txt" where "XXX" is for example "kb_E_50h_v3".
  6 +#
  7 +# Copyright (C) 2011 Jurgen Knodlseder
  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 math
  25 +import sys
  26 +
  27 +
  28 +# ================================= #
  29 +# Generate file function from model #
  30 +# ================================= #
  31 +def make_file_function(irfname, filename):
  32 + """
  33 + Generate file function from model.
  34 +
  35 + Parameters:
  36 + irfname - IRF filename
  37 + filename - File function filename
  38 + """
  39 + # Open performance file and file function
  40 + irf = open(irfname, "r")
  41 + file = open(filename, "w")
  42 +
  43 + # Read data
  44 + for row in irf:
  45 +
  46 + # Skip header and trailer
  47 + if row.find("log(E)") != -1:
  48 + continue
  49 + if row.find("--------") != -1:
  50 + break
  51 +
  52 + # Extract data
  53 + fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",")
  54 + logE = float(fields[0])
  55 + aeff = float(fields[1])
  56 + r68 = float(fields[2])
  57 + r80 = float(fields[3])
  58 + eres = float(fields[4])
  59 + bgd = float(fields[5])
  60 + sens = float(fields[6])
  61 +
  62 + # Compute energy (in MeV)
  63 + energy = math.pow(10.0, logE)*1.0e6
  64 +
  65 + # Compute solid angle of r68 region in sr
  66 + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r68)))
  67 +
  68 + # Compute background rate per steradian and MeV
  69 + bkg_rate = bgd/omega/energy
  70 +
  71 + # Write
  72 + file.write(str(energy)+" "+str(bkg_rate)+"\n")
  73 +
  74 + print logE, energy, bgd, omega, bkg_rate
  75 +
  76 + # Close files
  77 + irf.close()
  78 + file.close()
  79 +
  80 + # Return
  81 + return
  82 +
  83 +
  84 +# ================ #
  85 +# Main entry point #
  86 +# ================ #
  87 +if __name__ == '__main__':
  88 + """
  89 + Main entry point
  90 + """
  91 + # Get IRF from command line argument
  92 + if (len(sys.argv) > 1):
  93 + irf = sys.argv[1]
  94 + else:
  95 + irf = "kb_E_50h_v3"
  96 +
  97 + # Set parameters
  98 + caldb = "/usr/local/gamma/share/caldb/cta"
  99 +
  100 + # Set filenames
  101 + irfname = caldb+"/"+irf+".dat"
  102 + filename = "bkg_"+irf+".txt"
  103 +
  104 + # Read
  105 + make_file_function(irfname, filename)
... ...