Dmitry Malyshev / ctools

Commit 3ef43b0ee3f5d4a5b942c4d1775ac7cf9caf4eaf

Authored by knodlseder
1 parent 910cb4f3

renamed from sensivitity_pull.py

Inline Side-by-side
Showing 1 changed file with 209 additions and 0 deletions
examples/make_pull_at_sensitivity_limit.py 0 → 100755
  View file @3ef43b0
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script generates pull distributions for a source at the 5 sigma
  4 +# sensitivity limit. The script makes use of the "processing" Python module,
  5 +# if available, that can be used for parallel computing on multiple
  6 +# cores/CPUs.
  7 +#
  8 +# Required 3rd party modules:
  9 +# - processing (optional)
  10 +#
  11 +# Copyright (C) 2011 Jurgen Knodlseder
  12 +#
  13 +# This program is free software: you can redistribute it and/or modify
  14 +# it under the terms of the GNU General Public License as published by
  15 +# the Free Software Foundation, either version 3 of the License, or
  16 +# (at your option) any later version.
  17 +#
  18 +# This program is distributed in the hope that it will be useful,
  19 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  21 +# GNU General Public License for more details.
  22 +#
  23 +# You should have received a copy of the GNU General Public License
  24 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  25 +#
  26 +# ==========================================================================
  27 +from ctatools import *
  28 +from gammalib import *
  29 +from cspull import *
  30 +import sys
  31 +import time
  32 +
  33 +# Try importing processing module
  34 +has_processing = False
  35 +try:
  36 + import processing
  37 + has_processing = True
  38 +except:
  39 + pass
  40 +
  41 +
  42 +# ======================== #
  43 +# Create pull distribution #
  44 +# ======================== #
  45 +def create_pull(loge, emin, emax, models, ntrials=100, duration=180000.0, \
  46 + enumbins=0, log=False):
  47 + """
  48 + Create pull distribution
  49 +
  50 + Parameters:
  51 + emin - Minimum energy (TeV)
  52 + emax - Maximum energy (TeV)
  53 + Keywords:
  54 + log - Create log file(s)
  55 + """
  56 + # Generate output filename
  57 + outfile = "pull_"+loge+".dat"
  58 +
  59 + # Setup cspull tool
  60 + pull = cspull()
  61 + pull.models(models)
  62 + pull["outfile"] = outfile
  63 + pull["ntrials"] = ntrials
  64 + pull["caldb"] = "$GAMMALIB/share/caldb/cta"
  65 + pull["irf"] = "kb_E_50h_v3"
  66 + pull["ra"] = 83.6331
  67 + pull["dec"] = 22.0145
  68 + pull["emin"] = emin
  69 + pull["emax"] = emax
  70 + pull["enumbins"] = enumbins
  71 + pull["duration"] = duration
  72 +
  73 + # Optionally open the log file
  74 + if log:
  75 + pull.logFileOpen()
  76 +
  77 + # Run tool
  78 + pull.run()
  79 +
  80 + # Return
  81 + return
  82 +
  83 +
  84 +# ============ #
  85 +# Create model #
  86 +# ============ #
  87 +def create_model(flux, index=-2.48, fitidx=False):
  88 + """
  89 + Add standard CTA background model to observations container.
  90 + We use a simple power law here, scaled to Konrad's E configuration
  91 + performance table.
  92 +
  93 + Parameters:
  94 + flux - Flux in Crab units
  95 + """
  96 + # Define background model
  97 + bgd_radial = GCTAModelRadialGauss(3.0)
  98 + bgd_spectrum = GModelSpectralPlaw(61.8, -1.85)
  99 + bgd_spectrum["Prefactor"].scale(1.0e-6)
  100 + bgd_spectrum["PivotEnergy"].value(1.0)
  101 + bgd_spectrum["PivotEnergy"].scale(1.0e6)
  102 + if fitidx:
  103 + bgd_spectrum["Index"].free()
  104 + else:
  105 + bgd_spectrum["Index"].fix()
  106 + bgd_model = GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum)
  107 + bgd_model.name("Background")
  108 + bgd_model.instruments("CTA")
  109 +
  110 + # Define source spectrum
  111 + location = GSkyDir()
  112 + location.radec_deg(83.6331, 22.0145)
  113 + src_spatial = GModelSpatialPtsrc(location)
  114 + src_spectrum = GModelSpectralPlaw(flux, index)
  115 + src_spectrum["Prefactor"].scale(5.7e-16)
  116 + src_spectrum["PivotEnergy"].value(0.3)
  117 + src_spectrum["PivotEnergy"].scale(1.0e6)
  118 + if fitidx:
  119 + src_spectrum["Index"].free()
  120 + else:
  121 + src_spectrum["Index"].fix()
  122 + src_model = GModelPointSource(src_spatial, src_spectrum)
  123 + src_model.name("Test")
  124 +
  125 + # Add models to container
  126 + models = GModels()
  127 + models.append(bgd_model)
  128 + models.append(src_model)
  129 +
  130 + # Return model container
  131 + return models
  132 +
  133 +
  134 +# ======================== #
  135 +# Main routine entry point #
  136 +# ======================== #
  137 +if __name__ == '__main__':
  138 + """
  139 + Create pull distributions for the energy bands and fluxes given in the
  140 + file specified in the argument. The file needs to be created using the
  141 + cssens script.
  142 + """
  143 + # Get input arguments
  144 + usage = "make_pull_at_sensitivity_limit filename [max_threads]"
  145 + if len(sys.argv) < 2 or len(sys.argv) > 3:
  146 + print usage
  147 + sys.exit()
  148 +
  149 + # Get sensitivity file
  150 + filename = sys.argv[1]
  151 +
  152 + # Set maximum number of threads
  153 + if len(sys.argv) > 2:
  154 + max_threads = int(sys.argv[2])
  155 + else:
  156 + max_threads = 1
  157 +
  158 + # Load sensitivity file
  159 + file = open(filename, 'r')
  160 + reader = csv.DictReader(file)
  161 +
  162 + # Loop over rows
  163 + for row in reader:
  164 +
  165 + # Get parameters
  166 + loge = float(row['loge'])
  167 + if loge < 0:
  168 + loge = "m"+str(abs(loge))
  169 + else:
  170 + loge = "p"+str(abs(loge))
  171 + emin = float(row['emin'])
  172 + emax = float(row['emax'])
  173 + flux = float(row['crab'])
  174 +
  175 + # Generate model
  176 + models = create_model(flux)
  177 +
  178 + # Processing support?
  179 + if has_processing:
  180 +
  181 + # Wait until one thread has finished
  182 + while len(processing.activeChildren()) >= max_threads:
  183 + time.sleep(60)
  184 +
  185 + # Set arguments
  186 + args = (loge, emin, emax, models)
  187 + kwargs = {}
  188 +
  189 + # Generate pull distribution
  190 + p = processing.Process(target=create_pull, args=args, kwargs=kwargs)
  191 + p.start()
  192 + print "Process emin=%.4f emax=%.4f started." % (emin, emax)
  193 +
  194 + # Wait some time
  195 + time.sleep(1)
  196 +
  197 + # ... no
  198 + else:
  199 + create_pull(loge, emin, emax, models)
  200 +
  201 + # Processing support
  202 + if has_processing:
  203 +
  204 + # Wait until all threads finished
  205 + while len(processing.activeChildren()) > 0:
  206 + time.sleep(60)
  207 +
  208 + # Close file
  209 + file.close()
... ...