Thierry Reposeur / ctools

Commit 24bcb6e3e149ffa0f64318fef5df00748661ff37

Authored by knodlseder
1 parent 0645558a

initial version

Inline Side-by-side
Showing 1 changed file with 150 additions and 0 deletions
examples/ts_distributions.py 0 → 100755
  View file @24bcb6e
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script generates TS distributions for the standard energy bands that
  4 +# are also used for sensitivity computation. The script makes use of the
  5 +# "processing" Python module, if available, that can be used for parallel
  6 +# computing on multiple cores/CPUs.
  7 +#
  8 +# This script provides an illustration of how a ctools or cscript can be
  9 +# used from within Python. It further illustrates the possibility of
  10 +# parallel computing using the "processing" Python module.
  11 +#
  12 +# Copyright (C) 2011 Jurgen Knodlseder
  13 +#
  14 +# This program is free software: you can redistribute it and/or modify
  15 +# it under the terms of the GNU General Public License as published by
  16 +# the Free Software Foundation, either version 3 of the License, or
  17 +# (at your option) any later version.
  18 +#
  19 +# This program is distributed in the hope that it will be useful,
  20 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  22 +# GNU General Public License for more details.
  23 +#
  24 +# You should have received a copy of the GNU General Public License
  25 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  26 +#
  27 +# ==========================================================================
  28 +from ctatools import *
  29 +from gammalib import *
  30 +from cstsdist import *
  31 +import sys
  32 +import time
  33 +
  34 +# Try importing processing module
  35 +has_processing = False
  36 +try:
  37 + import processing
  38 + has_processing = True
  39 +except:
  40 + pass
  41 +
  42 +
  43 +# ====================== #
  44 +# Create TS distribution #
  45 +# ====================== #
  46 +def create_ts(loge, emin, emax, ntrials=3, duration=180000.0, \
  47 + enumbins=0, log=False):
  48 + """
  49 + Create TS distribution.
  50 +
  51 + Parameters:
  52 + loge - Logarithm of mean energy in TeV
  53 + emin - Minimum energy (TeV)
  54 + emax - Maximum energy (TeV)
  55 + Keywords:
  56 + log - Create log file(s)
  57 + """
  58 + # Generate output filename
  59 + outfile = "ts_"+str(loge)+".dat"
  60 +
  61 + # Setup cstsdist tool
  62 + tsdist = cstsdist()
  63 + tsdist["outfile"] = outfile
  64 + tsdist["ntrials"] = ntrials
  65 + tsdist["caldb"] = "$GAMMALIB/share/caldb/cta"
  66 + tsdist["irf"] = "kb_E_50h_v3"
  67 + tsdist["type"] = "point"
  68 + tsdist["index"] = -2.48
  69 + tsdist["offset"] = 0.0
  70 + tsdist["bkg"] = "$GAMMALIB/share/models/bkg_kb_E_50h_v3.txt"
  71 + tsdist["emin"] = emin
  72 + tsdist["emax"] = emax
  73 + tsdist["enumbins"] = enumbins
  74 + tsdist["duration"] = duration
  75 + tsdist["rad"] = 5.0
  76 + tsdist["npix"] = 200
  77 + tsdist["binsz"] = 0.05
  78 +
  79 + # Optionally open the log file
  80 + if log:
  81 + tsdist.logFileOpen()
  82 +
  83 + # Run tool
  84 + tsdist.run()
  85 +
  86 + # Return
  87 + return
  88 +
  89 +
  90 +# ======================== #
  91 +# Main routine entry point #
  92 +# ======================== #
  93 +if __name__ == '__main__':
  94 + """
  95 + Create TS distribution in a number of energy bands.
  96 + """
  97 + # Get input arguments
  98 + usage = "ts_distributions <max_threads>"
  99 + if len(sys.argv) < 1:
  100 + print usage
  101 + sys.exit()
  102 +
  103 + # Set maximum number of threads (default: 1)
  104 + if len(sys.argv) > 1:
  105 + max_threads = int(sys.argv[1])
  106 + else:
  107 + max_threads = 1
  108 +
  109 + # Loop over energy bands. The energy bands are those that are also
  110 + # used for sensitivity computation.
  111 + for ieng in range(21):
  112 +
  113 + # Set energies
  114 + loge = -1.7 + ieng * 0.2
  115 + emean = pow(10.0, loge)
  116 + emin = pow(10.0, loge-0.1)
  117 + emax = pow(10.0, loge+0.1)
  118 +
  119 + # Processing support?
  120 + if has_processing:
  121 +
  122 + # Wait until one thread has finished
  123 + while len(processing.activeChildren()) >= max_threads:
  124 + time.sleep(60)
  125 +
  126 + # Set arguments
  127 + args = (loge, emin, emax)
  128 + kwargs = {}
  129 +
  130 + # Generate pull distribution
  131 + p = processing.Process(target=create_ts, args=args, kwargs=kwargs)
  132 + p.start()
  133 + print "Process emin=%.4f emax=%.4f started." % (emin, emax)
  134 +
  135 + # Wait a short time to allow process to start
  136 + time.sleep(1)
  137 +
  138 + # ... no
  139 + else:
  140 + create_ts(loge, emin, emax)
  141 +
  142 + # Processing support
  143 + if has_processing:
  144 +
  145 + # Wait until all threads finished
  146 + while len(processing.activeChildren()) > 0:
  147 + time.sleep(60)
  148 +
  149 + # Close file
  150 + file.close()
... ...