diff --git a/.gitignore b/.gitignore index 50d115b..1c6ae15 100644 --- a/.gitignore +++ b/.gitignore @@ -72,6 +72,7 @@ src/ctulimit/ctulimit src/cterror/cterror test/pfiles/ test/caldb/ +test/prod3 test/csroot2caldb_* test/iactdata_* test/*.fits @@ -80,7 +81,9 @@ test/*.xml test/*.reg test/*.trs test/*.txt +test/*.eps test/test_cscripts.sh +test/test_examples.sh test/reports/*.xml test/butterfly.txt test/.coverage diff --git a/ChangeLog b/ChangeLog index bddbb13..a161fc6 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2016-06-21 Juergen Knoedlseder +2016-06-25 Juergen Knoedlseder * Version 1.1.0 released ======================== + Add unit tests for all example scripts Handle zero parameter errors in ctulimit (#1795) Compute counts cube weights in ctbin (#1789) Adapt to new stacked cube interface (#1788) diff --git a/NEWS b/NEWS index 23da3f5..db7bb1f 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.1.0 -21 June 2016 +25 June 2016 Introduction diff --git a/README.md b/README.md index 41cc4f9..dee57cc 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ctools information ================== -* Version: 1.1.0.dev1 (21 June 2016) +* Version: 1.1.0.dev1 (25 June 2016) * Author: Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu) * GammaLib dependency: 1.1.0.dev1 diff --git a/configure.ac b/configure.ac index c278b20..c191943 100644 --- a/configure.ac +++ b/configure.ac @@ -1025,6 +1025,7 @@ AC_CONFIG_FILES([Makefile test/test_cscripts.sh test/test_python_ctools.sh test/test_python_cscripts.sh + test/test_examples.sh pyext/Makefile pyext/setup.py pyext/ctools/__init__.py @@ -1038,6 +1039,7 @@ AC_CONFIG_FILES([Makefile [if test -f test/test_cscripts.sh; then chmod +x test/test_cscripts.sh; fi; if test -f test/test_python_ctools.sh; then chmod +x test/test_python_ctools.sh; fi; if test -f test/test_python_cscripts.sh; then chmod +x test/test_python_cscripts.sh; fi; + if test -f test/test_examples.sh; then chmod +x test/test_examples.sh; fi; if test -f testconf.sh; then chmod +x testconf.sh; fi]) AC_OUTPUT diff --git a/examples/Makefile.am b/examples/Makefile.am index e867429..6bb953d 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -1,12 +1,9 @@ ## Process this file with automake to produce Makefile.in # Scripts to be distributed -scripts = README \ - benchmark_cta_analysis.py \ - check_models.py \ - make_pull_at_sensitivity_limit.py \ +scripts = README.md \ + make_pointings.py \ make_spectrum.py \ - make_survey.py \ make_ts_distributions.py \ pipeline_binned_disk.py \ pipeline_binned_mem.py \ @@ -15,10 +12,13 @@ scripts = README \ pipeline_unbinned_disk.py \ pipeline_unbinned_mem.py \ show_butterfly.py \ + show_model.py \ show_pha.py \ + show_pull_evolution.py \ show_pull_histogram.py \ show_response.py \ show_sensitivity.py \ + show_spectrum.py \ show_ts_distribution.py # Set directories diff --git a/examples/README b/examples/README deleted file mode 100644 index a7fc211..0000000 --- a/examples/README +++ /dev/null @@ -1,86 +0,0 @@ -This folder contains Python example scripts that illustrate the -usage of the ctools. - -benchmark_cta_analysis.py - This script peforms a benchmark of the various CTA analysis types - (unbinned, binned, stacked), and if matplotlib is installed, creates - a plot of the benchmark results. - -check_models.py - Checks a variety of models using an unbinned analysis. - -make_pointings.py - This script creates pointings for CTA surveys. For the moment the - Galactic Plane Survey (gps), the Extragalactic Survey (extgal), - the Galactic Centre Survey (gc) and the LMC survey (lmc) are - implemented. The script outputs and ASCII file that can be used by - csobsdef. - -make_pull_at_sensitivity_limit.py - This script creates pull distributions as function of energy - for a test source that is at the 5 sigma sensitivity limit for - 50h of effective observing time. It takes as input file the - output file generated using cspull. - If the Python processing module is installed the script can make - use of multiple cores and runs in parallel. - -make_spectrum.py - This script shows how to make a spectrum using obsutils. - -make_survey.py - This script simulates a CTA galactic plane survey. - -make_ts_distributions.py - This script creates TS distributions as function of energy for - an off observation of 50h effective observing time. - If the Python processing module is installed the script can make - use of multiple cores and runs in parallel. - -benchmark_cta_analysis.py - Performs performance benchmarks of the various CTA analysis types. - -pipeline_binned_disk.py - Example pipeline for simulating events and performing a binned - analysis by saving intermediate results to disk. - -pipeline_binned_mem.py - Example pipeline for simulating events and performing a binned - analysis in memory. No files will be written to disk. - -pipeline_stacked_disk.py - Example pipeline for simulating events and performing a stacked - analysis by saving intermediate results to disk. - -pipeline_stacked_mem.py - Example pipeline for simulating events and performing a stacked - analysis in memory. No files will be written to disk. - -pipeline_unbinned_disk.py - Example pipeline for simulating events and performing an unbinned - analysis by saving intermediate results to disk. - -pipeline_unbinned_mem.py - Example pipeline for simulating events and performing an unbinned - analysis in memory. No files will be written to disk. - -show_butterfly.py - Displays a butterfly generated using the ctbutterfly tool. - -show_model.py - Displays a spectral model. - -show_pha.py - Displays a spectrum in PHA format. - -show_pull_histogram.py - Displays pull histograms generated by the cspull script. - -show_response.py - Displays one or several CTA response functions. - -show_sensitivity.py - Displays sensitivity information generated by the cssens script. - -show_ts_distribution.py - Displays Test Statistic histograms generated using the cstsdist - script. diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000..5aabb5b --- /dev/null +++ b/examples/README.md @@ -0,0 +1,74 @@ +ctools examples +=============== + +This folder contains Python example scripts that illustrate the usage of +the ctools. + +`generate_prod3_irfs/py` +> Generates Prod3 Instrument Response Functions from FITS files. + +`make_pointings.py` +> Creates pointings for CTA surveys. For the moment the Galactic Plane +> Survey (gps), the Extragalactic Survey (extgal), the Galactic Centre +> Survey (gc) and the LMC survey (lmc) are implemented. The script outputs +> and ASCII file that can be used by csobsdef. + +`make_spectrum.py` +> Simulates events and generates a spectrum from the using csspec. + +`make_ts_distributions.py` +> Creates TS distributions as function of energy for an OFF observaition +> of 50 hours effective observing time. If the Python processing module is +> installed the script can make use of multiple cores and runs in parallel. + +`pipeline_binned_disk.py` +> Example pipeline for simulating events and performing a binned analysis +> by saving intermediate results to disk. + +`pipeline_binned_mem.py` +> Example pipeline for simulating events and performing a binned analysis +> in memory. No files will be written to disk. + +`pipeline_stacked_disk.py` +> Example pipeline for simulating events and performing a stacked analysis +> by saving intermediate results to disk. + +`pipeline_stacked_mem.py` +> Example pipeline for simulating events and performing a stacked analysis +> in memory. No files will be written to disk. + +`pipeline_unbinned_disk.py` +> Example pipeline for simulating events and performing an unbinned analysis +> by saving intermediate results to disk. + +`pipeline_unbinned_mem.py` +> Example pipeline for simulating events and performing an unbinned analysis +> in memory. No files will be written to disk. + +`show_butterfly.py` +> Displays a butterfly generated using the ctbutterfly tool. + +`show_model.py` +> Displays a spectral model. + +`show_pha.py` +> Displays a spectrum in PHA format. + +`show_pull_evolution.py` +> Displays the evolution of the mean and rms pull as function of number of +> trials. + +`show_pull_histogram.py` +> Displays pull histograms generated by the cspull script. + +`show_response.py` +> Displays one or several CTA response functions. + +`show_sensitivity.py` +> Displays the sensitivity generated by the cssens script. + +`show_spectrum.py` +> Displays spectrum generated by the csspec script. + +`show_ts_distribution.py` +> Displays Test Statistic histograms generated using the cstsdist script. diff --git a/examples/generate_prod3_irfs.py b/examples/generate_prod3_irfs.py index fbb263d..174d318 100755 --- a/examples/generate_prod3_irfs.py +++ b/examples/generate_prod3_irfs.py @@ -27,15 +27,30 @@ import cscripts # ======================= # def set_prod3_desy(): """ - Set Prod3 DESY database. + Set Prod3 DESY database + + This function sets the calibration database. By adding dictionaries to + the db list a number of calibrations can be defined for which IRFs + should be generated. All calibrations will be added to the calibration + database. """ # Set database attributes - path = "/project-data/cta/performance/prod3" - psftype = "Gauss" + #path = '/project-data/cta/performance/prod3' + #file = 'DESY.d20160304.V5.ID0NIM2.prod3-paranalp05-NN.S.3HB1-ND-3.180000s.root' + path = 'data' + file = 'irf.root' + psftype = 'Gauss' + oversample = 3 + #norm1d = True # Will normalise the on-axis 2D IRFs on the 1D IRFs + norm1d = False + bgdinfill = False + bgdethres = 10.0 # Set database content - db = [{'inst': "prod3", 'id': "South_50h", 'psftype': psftype, 'path': path, - 'file': "DESY.d20160304.V5.ID0NIM2.prod3-paranalp05-NN.S.3HB1-ND-3.180000s.root"} + db = [{'inst': 'prod3', 'id': 'South_50h', 'psftype': psftype, + 'oversample': oversample, 'norm1d': norm1d, + 'bgdinfill': bgdinfill, 'bgdethres': bgdethres, + 'path': path, 'file': file} ] # Return database @@ -66,23 +81,43 @@ if __name__ == '__main__': bgdscale = entry['bgdscale'] else: bgdscale = 1.0 + if entry.has_key('bgdethres'): + bgdethres = entry['bgdethres'] + else: + bgdethres = 1000.0 + if entry.has_key('bgdinfill'): + bgdinfill = entry['bgdinfill'] + else: + bgdinfill = False + if entry.has_key('oversample'): + oversample = entry['oversample'] + else: + oversample = 1 + if entry.has_key('norm1d'): + norm1d = entry['norm1d'] + else: + norm1d = False # Set filename - filename = entry['path']+"/"+entry['file'] + filename = entry['path']+'/'+entry['file'] # Set-up csroot2caldb caldb = cscripts.csroot2caldb() - caldb["infile"] = filename - caldb["outdir"] = "prod3" - caldb["inst"] = entry['inst'] - caldb["id"] = entry['id'] - caldb["version"] = entry['inst'] - caldb["analysis"] = "DESY" - caldb["psftype"] = entry['psftype'] - caldb["rebin"] = rebin - caldb["eascale"] = eascale - caldb["bgdscale"] = bgdscale - caldb["logfile"] = "generate_prod3_irfs.log" + caldb['infile'] = filename + caldb['outdir'] = 'prod3' + caldb['inst'] = entry['inst'] + caldb['id'] = entry['id'] + caldb['version'] = entry['inst'] + caldb['analysis'] = 'DESY' + caldb['psftype'] = entry['psftype'] + caldb['norm1d'] = norm1d + caldb['rebin'] = rebin + caldb['eascale'] = eascale + caldb['bgdscale'] = bgdscale + caldb['bgdoversample'] = oversample + caldb['bgdethres'] = bgdethres + caldb['bgdinfill'] = bgdinfill + caldb['logfile'] = 'generate_prod3_irfs.log' # Add CALDB entry caldb.execute() diff --git a/examples/make_pointings.py b/examples/make_pointings.py index 9005284..67b0ab6 100755 --- a/examples/make_pointings.py +++ b/examples/make_pointings.py @@ -1,13 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script generates pointing patterns for CTA observations. +# Generates pointing patterns for CTA observations # -# Usage: -# ./make_pointings.py -# -# ========================================================================== -# -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -33,10 +28,30 @@ import gammalib # ============================================================ # def get_positions(xmin, xmax, ymin, ymax, step): """ + Get pointing positions for an interleaved pattern + Get positions for a patch of interleaved pointings on the sky. The function does respect the latitude dependence of the pointing distance. The step parameter is for a latitude of zero. + + Parameters + ---------- + xmin : float + Minimum longitude (deg) + xmax : float + Maximum longitude (deg) + ymin : float + Minimum latitude (deg) + ymax : float + Maximum latitude (deg) + step : float + Pointing step size (deg) + + Returns + ------- + positions : list of dict + List of dictionaries with x and y positions """ # Initialise positions positions = [] @@ -46,7 +61,7 @@ def get_positions(xmin, xmax, ymin, ymax, step): ny = int((ymax-ymin)/ystep+1.5) rescale = (ymax-ymin)/(ystep*(ny-1)) ystep *= rescale - print("Number of rows ....: %d" % (ny)) + print('Number of rows ....: %d' % (ny)) # Loop over rows y = ymin @@ -63,7 +78,6 @@ def get_positions(xmin, xmax, ymin, ymax, step): x = xmin else: x = xmin + 0.5 * xstep - #print(row, y, nx) # Append pointings for pnt in range(nx): @@ -86,7 +100,33 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5, separation=3.0, hours=100.0, site=None, lst=True, autodec=0.0): """ - Setup pointing patch on the sky. + Setup pointing patch on the sky + + Parameters + ---------- + lmin : float, optional + Minimum Galactic Longitude (deg) + lmax : float, optional + Maximum Galactic Longitude (deg) + bmin : float, optional + Minimum Galactic Latitude (deg) + bmax : float, optional + Maximum Galactic Latitude (deg) + separation : float, optional + Pointing separation (deg) + hours : float, optional + Total observation duration (h) + site : str, optional + Array site (one of 'Automatic', 'South', 'North') + lst : bool, optional + Use LSTs + autodec : float, optional + Declination for switching automatically between South and North site + + Returns + ------- + obsdef : list of dict + List of pointing definitions """ # Initialise observation definition obsdef = [] @@ -112,30 +152,30 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5, if site != None: # Handle automatic site switching - if site == "Automatic": + if site == 'Automatic': pos = gammalib.GSkyDir() pos.lb_deg(pnt['x'], pnt['y']) dec = pos.dec_deg() if (dec >= autodec): - obs_site = "North" + obs_site = 'North' else: - obs_site = "South" + obs_site = 'South' else: obs_site = site # Set site IRF - caldb = "prod2" - if obs_site == "North": + caldb = 'prod2' + if obs_site == 'North': if lst: - irf = "North_50h" + irf = 'North_50h' else: - irf = "North_50h" + irf = 'North_50h' exposure_north += duration else: if lst: - irf = "South_50h" + irf = 'South_50h' else: - irf = "South_50h" + irf = 'South_50h' exposure_south += duration # Set IRF information @@ -146,11 +186,11 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5, obsdef.append(obs) # Dump statistics - print("Number of pointings: %d (%.2f sec)" % (n_pnt,duration)) - print("South array .......: %.2f hours (%.2f%%)" % + print('Number of pointings: %d (%.2f s)' % (n_pnt,duration)) + print('South array .......: %.2f hours (%.2f%%)' % (exposure_south/3600.0, 100.0*exposure_south/(exposure_south+exposure_north))) - print("North array .......: %.2f hours (%.2f%%)" % + print('North array .......: %.2f hours (%.2f%%)' % (exposure_north/3600.0, 100.0*exposure_north/(exposure_south+exposure_north))) @@ -163,7 +203,23 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5, # =========================== # def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True): """ - Setup Galactic plane survey. + Setup Galactic plane survey + + Parameters + ---------- + separation : float, optional + Pointing separation (deg) + bmin : float, optional + Minimum Galactic Latitude (deg) + bmax : float, optional + Maximum Galactic Latitude (deg) + lst : bool, optional + Use LSTs + + Returns + ------- + obsdef : list of dict + List of pointing definitions """ # Initialise observation definition obsdef = [] @@ -171,27 +227,27 @@ def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True): # Add inner region South obsdef.extend(set_patch(lmin=-60.0, lmax=60.0, bmin=bmin, bmax=bmax, separation=separation, hours=780, - site="South", lst=lst)) + site='South', lst=lst)) # Add Vela & Carina region obsdef.extend(set_patch(lmin=240.0, lmax=300.0, bmin=bmin, bmax=bmax, separation=separation, hours=180, - site="South", lst=lst)) + site='South', lst=lst)) # Add 210-240 region obsdef.extend(set_patch(lmin=210.0, lmax=240.0, bmin=bmin, bmax=bmax, separation=separation, hours=60, - site="South", lst=lst)) + site='South', lst=lst)) # Add Cygnus, Perseus obsdef.extend(set_patch(lmin=60.0, lmax=150.0, bmin=bmin, bmax=bmax, separation=separation, hours=450, - site="North", lst=lst)) + site='North', lst=lst)) # Add Anticentre obsdef.extend(set_patch(lmin=150.0, lmax=210.0, bmin=bmin, bmax=bmax, separation=separation, hours=150, - site="North", lst=lst)) + site='North', lst=lst)) # Return observation definition return obsdef @@ -203,6 +259,18 @@ def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True): def set_extgal(separation=3.0, lst=True): """ Setup Extragalactic survey. + + Parameters + ---------- + separation : float, optional + Pointing separation (deg) + lst : bool, optional + Use LSTs + + Returns + ------- + obsdef : list of dict + List of pointing definitions """ # Initialise observation definition obsdef = [] @@ -210,7 +278,7 @@ def set_extgal(separation=3.0, lst=True): # Set patch obsdef.extend(set_patch(lmin=-90.0, lmax=90.0, bmin=+5.0, bmax=+88.0, separation=separation, hours=1000, - site="Automatic", lst=lst, autodec=10.0)) + site='Automatic', lst=lst, autodec=10.0)) # Return observation definition return obsdef @@ -222,6 +290,16 @@ def set_extgal(separation=3.0, lst=True): def set_gc(lst=True): """ Setup Galactic centre KSP. + + Parameters + ---------- + lst : bool, optional + Use LSTs + + Returns + ------- + obsdef : list of dict + List of pointing definitions """ # Initialise observation definition obsdef = [] @@ -229,12 +307,12 @@ def set_gc(lst=True): # Central wobble obsdef.extend(set_patch(lmin=-0.5, lmax=0.5, bmin=-0.5, bmax=0.5, separation=0.1, hours=525, - site="South", lst=lst)) + site='South', lst=lst)) # Extended region obsdef.extend(set_patch(lmin=-10.0, lmax=10.0, bmin=-10.0, bmax=10.0, separation=1.5, hours=300, - site="South", lst=lst)) + site='South', lst=lst)) # Return observation definition return obsdef @@ -246,6 +324,18 @@ def set_gc(lst=True): def set_lmc(hours=250.0, lst=True): """ Setup LMC KSP. + + Parameters + ---------- + hours : float, optional + Total observation duration (h) + lst : bool, optional + Use LSTs + + Returns + ------- + obsdef : list of dict + List of pointing definitions """ # Initialise observation definition obsdef = [] @@ -275,11 +365,11 @@ def set_lmc(hours=250.0, lst=True): obs = {'lon': lon, 'lat': lat, 'duration': duration} # Add IRF - caldb = "prod2" + caldb = 'prod2' if lst: - irf = "South_50h" + irf = 'South_50h' else: - irf = "South_50h" + irf = 'South_50h' obs['caldb'] = caldb obs['irf'] = irf @@ -287,7 +377,7 @@ def set_lmc(hours=250.0, lst=True): obsdef.append(obs) # Dump statistics - print("Number of pointings: %d (%.2f sec)" % (n_pnt,duration)) + print('Number of pointings: %d (%.2f s)' % (n_pnt,duration)) # Return observation definition return obsdef @@ -298,31 +388,38 @@ def set_lmc(hours=250.0, lst=True): # ======================================= # def write_obsdef(filename, obsdef): """ - Write observation definition file. + Write observation definition file + + Parameters + ---------- + filename : str + Observation definition file name + obsdef : list of dict + List of pointing definitions """ # Open file file = open(filename, 'w') # Write header - file.write("ra,dec,duration,caldb,irf\n") + file.write('ra,dec,duration,caldb,irf\n') # Loop over pointings for obs in obsdef: # If we have lon,lat then convert into RA,Dec - if "lon" in obs and "lat" in obs: - lon = obs["lon"] - lat = obs["lat"] + if 'lon' in obs and 'lat' in obs: + lon = obs['lon'] + lat = obs['lat'] dir = gammalib.GSkyDir() dir.lb_deg(lon,lat) ra = dir.ra_deg() dec = dir.dec_deg() else: - ra = obs["ra"] - dec = obs["dec"] + ra = obs['ra'] + dec = obs['dec'] # Write information - file.write("%8.4f,%8.4f,%.4f,%s,%s\n" % + file.write('%8.4f,%8.4f,%.4f,%s,%s\n' % (ra, dec, obs['duration'], obs['caldb'], obs['irf'])) # Close file @@ -336,15 +433,13 @@ def write_obsdef(filename, obsdef): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - CTA pointing pattern generation. - """ + # Initialise flags need_help = False # Test for command line arguments if (len(sys.argv) > 1): - if sys.argv[1] == "-h": + if sys.argv[1] == '-h': need_help = True else: obsname = sys.argv[1] @@ -353,37 +448,37 @@ if __name__ == '__main__': # Print help if needed and exit if need_help: - print("Usage: make_pointing.py [OPTIONS]") - print(" -h Display this usage message") - print(" gps Galactic plane survey (2 row scheme)") - print(" gps3 Galactic plane survey (3 row scheme)") - print(" extgal Extragalactic survey") - print(" gc Galactic centre survey") - print(" lmc LMC survey") + print('Usage: make_pointing.py [OPTIONS]') + print(' -h Display this usage message') + print(' gps Galactic plane survey (2 row scheme)') + print(' gps3 Galactic plane survey (3 row scheme)') + print(' extgal Extragalactic survey') + print(' gc Galactic centre survey') + print(' lmc LMC survey') sys.exit() # Galactic plane survey - if obsname == "gps": + if obsname == 'gps': obsdef = set_gps(lst=True) - write_obsdef("gps.dat", obsdef) - elif obsname == "gps3": + write_obsdef('gps.dat', obsdef) + elif obsname == 'gps3': obsdef = set_gps(separation=1.5, lst=True) - write_obsdef("gps3.dat", obsdef) + write_obsdef('gps3.dat', obsdef) # Extragalactic survey - elif obsname == "extgal": + elif obsname == 'extgal': obsdef = set_extgal(lst=True) - write_obsdef("extgal.dat", obsdef) + write_obsdef('extgal.dat', obsdef) # Galactic centre - elif obsname == "gc": + elif obsname == 'gc': obsdef = set_gc(lst=True) - write_obsdef("gc.dat", obsdef) + write_obsdef('gc.dat', obsdef) # LMC - elif obsname == "lmc": + elif obsname == 'lmc': obsdef = set_lmc(lst=True) - write_obsdef("lmc.dat", obsdef) + write_obsdef('lmc.dat', obsdef) # Invalid pattern else: diff --git a/examples/make_spectrum.py b/examples/make_spectrum.py index b2ced97..b258186 100755 --- a/examples/make_spectrum.py +++ b/examples/make_spectrum.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This script shows how to make a spectrum using obsutils. +# Makes spectrum and optionally plot it # # Copyright (C) 2014-2016 Juergen Knoedlseder # @@ -18,11 +18,13 @@ # along with this program. If not, see . # # ========================================================================== -#import gammalib +import sys import ctools import cscripts try: import matplotlib.pyplot as plt + plt.figure() + plt.close() has_matplotlib = True except: has_matplotlib = False @@ -33,12 +35,17 @@ except: # =============== # def make_spectrum(): """ - Make a spectrum. + Simulate events and generate a source spectrum + + Returns + ------- + spectrum : `~gammalib.GFits` + Spectrum FITS file """ # Set script parameters - model_name = "${CTOOLS}/share/models/crab.xml" - caldb = "prod2" - irf = "South_50h" + model_name = 'data/crab.xml' + caldb = 'prod2' + irf = 'South_0.5h' ra = 83.63 dec = 22.01 rad_sim = 3.0 @@ -50,36 +57,32 @@ def make_spectrum(): # Simulate events sim = ctools.ctobssim() - sim["inmodel"] = model_name - sim["caldb"] = caldb - sim["irf"] = irf - sim["ra"] = ra - sim["dec"] = dec - sim["rad"] = rad_sim - sim["tmin"] = tstart - sim["tmax"] = tstop - sim["emin"] = emin - sim["emax"] = emax + sim['inmodel'] = model_name + sim['caldb'] = caldb + sim['irf'] = irf + sim['ra'] = ra + sim['dec'] = dec + sim['rad'] = rad_sim + sim['tmin'] = tstart + sim['tmax'] = tstop + sim['emin'] = emin + sim['emax'] = emax sim.run() - # Generate an energy binning - #e_min = gammalib.GEnergy(emin, "TeV") - #e_max = gammalib.GEnergy(emax, "TeV") - #ebounds = gammalib.GEbounds(10, e_min, e_max) - - # Setup csspec run + # Create spectrum spec = cscripts.csspec(sim.obs()) - spec["srcname"] = "Crab" - spec["outfile"] = "spectrum.fits" - spec["expcube"] = "NONE" - spec["psfcube"] = "NONE" - spec["bkgcube"] = "NONE" - spec["edisp"] = False - spec["emin"] = emin - spec["emax"] = emax - spec["enumbins"] = enumbins - spec["binned"] = False + spec['srcname'] = 'Crab' + spec['outfile'] = 'example_spectrum.fits' + spec['expcube'] = 'NONE' + spec['psfcube'] = 'NONE' + spec['bkgcube'] = 'NONE' + spec['edisp'] = False + spec['emin'] = emin + spec['emax'] = emax + spec['enumbins'] = enumbins + spec['ebinalg'] = 'LOG' spec.run() + spec.save() # Get copy of spectrum spectrum = spec.spectrum().copy() @@ -91,19 +94,26 @@ def make_spectrum(): # ============= # # Plot spectrum # # ============= # -def plot_spectrum(spectrum): +def plot_spectrum(spectrum, file=False): """ - Plot spectrum. + Plot spectrum + + Parameters + ---------- + spectrum : `~gammalib.GFits` + Spectrum FITS file + file : bool, optional + Plot into a file """ - # Read spectrum file + # Extract columns from spectrum file table = spectrum.table(1) - c_energy = table["Energy"] - c_ed = table["ed_Energy"] - c_eu = table["eu_Energy"] - c_flux = table["Flux"] - c_eflux = table["e_Flux"] - c_ts = table["TS"] - c_upper = table["UpperLimit"] + c_energy = table['Energy'] + c_ed = table['ed_Energy'] + c_eu = table['eu_Energy'] + c_flux = table['Flux'] + c_eflux = table['e_Flux'] + c_ts = table['TS'] + c_upper = table['UpperLimit'] # Initialise arrays to be filled energies = [] @@ -115,7 +125,7 @@ def plot_spectrum(spectrum): ul_ed_engs = [] ul_eu_engs = [] ul_flux = [] - + # Loop over rows of the file nrows = table.nrows() for row in range(nrows): @@ -146,33 +156,43 @@ def plot_spectrum(spectrum): # Create figure plt.figure() - plt.title("Crab spectrum") + plt.title('Crab spectrum') # Plot the spectrum plt.loglog() plt.grid() - plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], fmt='ro') - plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], yerr=1.0e-11, uplims=True, fmt='ro') - plt.xlabel("Energy (TeV)") - plt.ylabel(r"E$^{2}$dN/dE (erg cm$^{-2}$ s$^{-1}$)") - - # Show spectrum - plt.show() + plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], + fmt='ro') + plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], + yerr=1.0e-11, uplims=True, fmt='ro') + plt.xlabel('Energy (TeV)') + plt.ylabel(r'E$^{2}$dN/dE (erg cm$^{-2}$ s$^{-1}$)') + + # Show spectrum or save it into file + if file: + plt.savefig('example_spectrum.eps') + else: + plt.show() # Return return -#==========================# +# ======================== # # Main routine entry point # -#==========================# +# ======================== # if __name__ == '__main__': - """ - Generate a Crab spectrum. - """ + + # Initialise option + file = False + + # Get any command line arguments + if (len(sys.argv) > 1): + file = True + # Generate spectrum spectrum = make_spectrum() - # Plot spectrum + # If matplotlib is installed then plot spectrum if has_matplotlib: - plot_spectrum(spectrum) + plot_spectrum(spectrum, file=file) diff --git a/examples/make_ts_distributions.py b/examples/make_ts_distributions.py index 44b7c91..a1412ca 100755 --- a/examples/make_ts_distributions.py +++ b/examples/make_ts_distributions.py @@ -1,18 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script generates TS distributions for the standard energy bands that -# are also used for sensitivity computation. The script makes use of the -# "processing" Python module, if available, that can be used for parallel -# computing on multiple cores/CPUs. +# Generate TS distributions as function of energy for an OFF observation # -# This script provides an illustration of how a ctools or cscript can be -# used from within Python. It further illustrates the possibility of -# parallel computing using the "processing" Python module. -# -# Required 3rd party modules: -# - processing (optional) -# -# Copyright (C) 2011-2015 Jurgen Knodlseder +# Copyright (C) 2011-2016 Jurgen Knodlseder # # 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 @@ -30,15 +20,14 @@ # ========================================================================== import sys import time -from cscripts import cstsdist +import cscripts -# Try importing processing module -has_processing = False +# Optionally import processing module try: import processing has_processing = True except: - pass + has_processing = False # ====================== # @@ -49,35 +38,44 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, """ Create TS distribution. - Parameters: - loge - Logarithm of mean energy in TeV - emin - Minimum energy (TeV) - emax - Maximum energy (TeV) - Keywords: - log - Create log file(s) + Parameters + ---------- + loge : float + Logarithm of mean energy (TeV) + emin : float + Minimum energy (TeV) + emax : float + Maximum energy (TeV) + ntrials : int, optional + Number of MC samples + durations : float, optional + Observation duration (s) + enumbins : int, optional + Number of energy bins + log : bool, optional + Create log file(s) """ # Generate output filename - outfile = "ts_"+str(loge)+".dat" + outfile = 'ts_'+str(loge)+'.dat' # Setup cstsdist tool - tsdist = cstsdist() - tsdist["inmodel"] = "$CTOOLS/share/models/crab.xml" - tsdist["srcname"] = "Crab" - tsdist["outfile"] = outfile - tsdist["ntrials"] = ntrials - tsdist["caldb"] = "prod2" - tsdist["irf"] = "South_50h" - tsdist["ra"] = 83.63 - tsdist["dec"] = 22.01 - tsdist["emin"] = emin - tsdist["emax"] = emax - tsdist["enumbins"] = enumbins - tsdist["tmin"] = 0 - tsdist["tmax"] = duration - tsdist["rad"] = 5.0 - tsdist["npix"] = 200 - tsdist["binsz"] = 0.05 - #tsdist["debug"].boolean(True) + tsdist = cscripts.cstsdist() + tsdist['inmodel'] = 'data/crab.xml' + tsdist['srcname'] = 'Crab' + tsdist['outfile'] = outfile + tsdist['ntrials'] = ntrials + tsdist['caldb'] = 'prod2' + tsdist['irf'] = 'South_50h' + tsdist['ra'] = 83.63 + tsdist['dec'] = 22.01 + tsdist['emin'] = emin + tsdist['emax'] = emax + tsdist['enumbins'] = enumbins + tsdist['tmin'] = 0 + tsdist['tmax'] = duration + tsdist['rad'] = 5.0 + tsdist['npix'] = 200 + tsdist['binsz'] = 0.05 # Optionally open the log file if log: @@ -94,11 +92,9 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Create TS distribution in a number of energy bands. - """ + # Get input arguments - usage = "make_ts_distributions [-n ntrials] [-e enumbins] [-m max_threads]" + usage = 'make_ts_distributions.py [-n ntrials] [-e enumbins] [-m max_threads]' if len(sys.argv) < 1: print(usage) sys.exit() @@ -113,7 +109,7 @@ if __name__ == '__main__': pars = [{'option': '-n', 'value': ntrials}, {'option': '-e', 'value': enumbins}, {'option': '-d', 'value': duration}, - {'option': '-m', 'value': max_threads}] + {'option': '-m', 'value': max_threads}] # Gather parameters from command line i = 1 @@ -123,7 +119,7 @@ if __name__ == '__main__': for par in pars: if sys.argv[i] == par['option']: if len(sys.argv) > i+1: - i += 1 + i += 1 try: par['value'] = int(sys.argv[i]) except: @@ -152,11 +148,12 @@ if __name__ == '__main__': emin = pow(10.0, loge-0.1) emax = pow(10.0, loge+0.1) if loge < 0: - loge = "m"+str(abs(loge)) + loge = 'm'+str(abs(loge)) else: - loge = "p"+str(abs(loge)) + loge = 'p'+str(abs(loge)) - # Processing support? + # If processing support is available then start max_threads processes + # in parallel if has_processing: # Wait until one thread has finished @@ -171,19 +168,17 @@ if __name__ == '__main__': # Generate pull distribution p = processing.Process(target=create_ts, args=args, kwargs=kwargs) p.start() - print("Process emin=%.4f emax=%.4f started." % (emin, emax)) + print('Process emin=%.4f emax=%.4f started.' % (emin, emax)) # Wait a short time to allow process to start time.sleep(1) - # ... no + # ... otherwise simply start a single process else: create_ts(loge, emin, emax, ntrials=ntrials, enumbins=enumbins, - duration=duration) + duration=duration) - # Processing support + # If processing support is available then wait until all threads finished if has_processing: - - # Wait until all threads finished while len(processing.activeChildren()) > 0: time.sleep(10) diff --git a/examples/pipeline_binned_disk.py b/examples/pipeline_binned_disk.py index 89c3e4d..f89ca25 100755 --- a/examples/pipeline_binned_disk.py +++ b/examples/pipeline_binned_disk.py @@ -1,15 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script illustrates how to perform a binned CTA analysis based on -# simulated CTA data. You may use and adapt this script to implement your -# own pipeline. +# Perform binned analysis of simulated CTA data. # -# Usage: -# ./pipeline_binned_disk.py -# -# ========================================================================== -# -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -33,25 +26,41 @@ from cscripts import obsutils # ================== # # Setup observations # # ================== # -def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, - emin=0.1, emax=100.0, rad=5.0, duration=1800.0, - deadc=0.95, - caldb="prod2", irf="South_50h"): +def setup_observations(pattern='four', ra=83.63, dec=22.01, offset=1.5, + emin=0.1, emax=100.0, rad=5.0, duration=180.0, + deadc=0.95, caldb='prod2', irf='South_0.5h'): """ - Returns an observation container. - - Keywords: - pattern - Pointing pattern, either "single" or "four" - ra - RA of pattern centre [deg] (default: 83.6331) - dec - DEC of pattern centre [deg] (default: 22.0145) - offset - Offset between observations of pattern [deg] (default: 1.5) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - duration - Duration of one CTA observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - caldb - Calibration database path (default: "dummy") - irf - Instrument response function (default: cta_dummy_irf) + Returns an observation container + + Parameters + ---------- + pattern : str, optional + Pointing pattern, either 'single' or 'four' + ra : float, optional + Right Ascension of pattern centre (deg) + dec : float, optional + Declination of pattern centre (deg) + offset : float, optional + Offset between observations of pattern (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + rad : float, optional + ROI radius used for analysis (deg) + duration : float, optional + Duration of one CTA observation (s) + deadc : float, optional + Deadtime correction factor + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Set list of observations obs_def_list = obsutils.set_obs_patterns(pattern, @@ -75,12 +84,21 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, # =========== # # Setup model # # =========== # -def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): +def setup_model(obs, model='data/crab.xml'): """ - Setup model for analysis. - - Keywords: - model - Model Xml file + Setup model for analysis + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + model : str, optional + Model definition XML file + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Append model from file to observation container obs.models(gammalib.GModels(model)) @@ -94,56 +112,71 @@ def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): # ================================ # def run_pipeline(obs, emin=0.1, emax=100.0, enumbins=20, nxpix=200, nypix=200, binsz=0.02, - coordsys="CEL", proj="CAR", - model="${CTOOLS}/share/models/crab.xml", - caldb="prod2", irf="South_50h", + coordsys='CEL', proj='CAR', + model='data/crab.xml', + caldb='prod2', irf='South_0.5h', debug=False): """ - Simulation and binned analysis pipeline. - - Keywords: - emin - Minimum energy of cube [TeV] (default: 0.1) - emax - Maximum energy of cube [TeV] (default: 100.0) - enumbins - Number of energy bins in cube (default: 20) - nxpix - Number of RA pixels in cube (default: 200) - nypix - Number of DEC pixels in cube (default: 200) - binsz - Spatial cube bin size [deg] (default: 0.02) - coordsys - Cube coordinate system (CEL or GAL) - proj - Cube World Coordinate System (WCS) projection - model - Model Xml file - caldb - Calibration database path (default: "dummy") - irf - Instrument response function (default: cta_dummy_irf) - debug - Enable debugging (default: False) + Simulation and binned analysis pipeline + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + enumbins : int, optional + Number of energy bins + nxpix : int, optional + Number of pixels in X axis + nypix : int, optional + Number of pixels in Y axis + binsz : float, optional + Pixel size (deg) + coordsys : str, optional + Coordinate system + proj : str, optional + Coordinate projection + model : str, optional + Model definition XML file + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + debug : bool, optional + Debug function """ # Simulate events sim = ctools.ctobssim(obs) - sim["debug"] = debug - sim["outevents"] = "obs.xml" + sim['debug'] = debug + sim['outevents'] = 'obs.xml' sim.execute() # Bin events by looping over all observations in the container - sim_obs = gammalib.GObservations("obs.xml") + sim_obs = gammalib.GObservations('obs.xml') obs = gammalib.GObservations() for run in sim_obs: # Get event filename and set counts cube filename - eventfile = run.eventfile() - cubefile = "cube_"+eventfile + eventfile = run.eventfile().url() + cubefile = 'cube_'+eventfile # Bin events for that observation bin = ctools.ctbin() - bin["inobs"] = eventfile - bin["outcube"] = cubefile - bin["ebinalg"] = "LOG" - bin["emin"] = emin - bin["emax"] = emax - bin["enumbins"] = enumbins - bin["nxpix"] = nxpix - bin["nypix"] = nypix - bin["binsz"] = binsz - bin["coordsys"] = coordsys - bin["usepnt"] = True - bin["proj"] = proj + bin['inobs'] = eventfile + bin['outcube'] = cubefile + bin['ebinalg'] = 'LOG' + bin['emin'] = emin + bin['emax'] = emax + bin['enumbins'] = enumbins + bin['nxpix'] = nxpix + bin['nypix'] = nypix + bin['binsz'] = binsz + bin['coordsys'] = coordsys + bin['usepnt'] = True + bin['proj'] = proj bin.execute() # Set observation ID @@ -156,19 +189,19 @@ def run_pipeline(obs, emin=0.1, emax=100.0, # Save XML file xml = gammalib.GXml() obs.write(xml) - xml.save("obs_cube.xml") + xml.save('obs_cube.xml') # Perform maximum likelihood fitting like = ctools.ctlike() - like["inobs"] = "obs_cube.xml" - like["inmodel"] = model - like["outmodel"] = "fit_results.xml" - like["expcube"] = "NONE" - like["psfcube"] = "NONE" - like["bkgcube"] = "NONE" - like["caldb"] = caldb - like["irf"] = irf - like["debug"] = True # Switch this always on for results in console + like['inobs'] = 'obs_cube.xml' + like['inmodel'] = model + like['outmodel'] = 'fit_results.xml' + like['expcube'] = 'NONE' + like['psfcube'] = 'NONE' + like['bkgcube'] = 'NONE' + like['caldb'] = caldb + like['irf'] = irf + like['debug'] = True # Switch this always on for results in console like.execute() # Return @@ -179,12 +212,11 @@ def run_pipeline(obs, emin=0.1, emax=100.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header - print("********************************************") - print("* CTA binned analysis pipeline *") - print("********************************************") + print('********************************************') + print('* CTA binned analysis pipeline *') + print('********************************************') # Setup observations obs = setup_observations() diff --git a/examples/pipeline_binned_mem.py b/examples/pipeline_binned_mem.py index f85a5b0..fa89ba4 100755 --- a/examples/pipeline_binned_mem.py +++ b/examples/pipeline_binned_mem.py @@ -1,15 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script illustrates how to perform a binned CTA analysis based on -# simulated CTA data. You may use and adapt this script to implement your -# own pipeline. +# Perform binned in-memory analysis of simulated CTA data. # -# Usage: -# ./pipeline_binned_mem.py -# -# ========================================================================== -# -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -33,25 +26,41 @@ from cscripts import obsutils # ================== # # Setup observations # # ================== # -def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, - emin=0.1, emax=100.0, rad=5.0, duration=1800.0, - deadc=0.95, - caldb="prod2", irf="South_50h"): +def setup_observations(pattern='four', ra=83.63, dec=22.01, offset=1.5, + emin=0.1, emax=100.0, rad=5.0, duration=180.0, + deadc=0.95, caldb='prod2', irf='South_0.5h'): """ - Returns an observation container. - - Keywords: - pattern - Pointing pattern, either "single" or "four" - ra - RA of pattern centre [deg] (default: 83.6331) - dec - DEC of pattern centre [deg] (default: 22.0145) - offset - Offset between observations of pattern [deg] (default: 1.5) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - duration - Duration of one CTA observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - caldb - Calibration database path (default: "dummy") - irf - Instrument response function (default: cta_dummy_irf) + Returns an observation container + + Parameters + ---------- + pattern : str, optional + Pointing pattern, either 'single' or 'four' + ra : float, optional + Right Ascension of pattern centre (deg) + dec : float, optional + Declination of pattern centre (deg) + offset : float, optional + Offset between observations of pattern (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + rad : float, optional + ROI radius used for analysis (deg) + duration : float, optional + Duration of one CTA observation (s) + deadc : float, optional + Deadtime correction factor + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Set list of observations obs_def_list = obsutils.set_obs_patterns(pattern, @@ -75,12 +84,21 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, # =========== # # Setup model # # =========== # -def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): +def setup_model(obs, model='data/crab.xml'): """ - Setup model for analysis. - - Keywords: - model - Model Xml file + Setup model for analysis + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + model : str, optional + Model definition XML file + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Append model from file to observation container obs.models(gammalib.GModels(model)) @@ -94,24 +112,36 @@ def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): # ================================ # def run_pipeline(obs, emin=0.1, emax=100.0, enumbins=20, nxpix=200, nypix=200, binsz=0.02, - coordsys="CEL", proj="CAR", debug=False): + coordsys='CEL', proj='CAR', debug=False): """ - Simulation and binned analysis pipeline. - - Keywords: - emin - Minimum energy of cube [TeV] (default: 0.1) - emax - Maximum energy of cube [TeV] (default: 100.0) - enumbins - Number of energy bins in cube (default: 20) - nxpix - Number of RA pixels in cube (default: 200) - nypix - Number of DEC pixels in cube (default: 200) - binsz - Spatial cube bin size [deg] (default: 0.02) - coordsys - Cube coordinate system (CEL or GAL) - proj - Cube World Coordinate System (WCS) projection - debug - Enable debugging (default: False) + Simulation and binned analysis pipeline + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + enumbins : int, optional + Number of energy bins + nxpix : int, optional + Number of pixels in X axis + nypix : int, optional + Number of pixels in Y axis + binsz : float, optional + Pixel size (deg) + coordsys : str, optional + Coordinate system + proj : str, optional + Coordinate projection + debug : bool, optional + Debug function """ # Simulate events sim = ctools.ctobssim(obs) - sim["debug"] = debug + sim['debug'] = debug sim.run() # Bin events by looping over all observations in the container @@ -125,16 +155,16 @@ def run_pipeline(obs, emin=0.1, emax=100.0, # Bin events for that observation bin = ctools.ctbin(container) - bin["ebinalg"] = "LOG" - bin["emin"] = emin - bin["emax"] = emax - bin["enumbins"] = enumbins - bin["nxpix"] = nxpix - bin["nypix"] = nypix - bin["binsz"] = binsz - bin["coordsys"] = coordsys - bin["usepnt"] = True - bin["proj"] = proj + bin['ebinalg'] = 'LOG' + bin['emin'] = emin + bin['emax'] = emax + bin['enumbins'] = enumbins + bin['nxpix'] = nxpix + bin['nypix'] = nypix + bin['binsz'] = binsz + bin['coordsys'] = coordsys + bin['usepnt'] = True + bin['proj'] = proj bin.run() # Append result to observations @@ -142,7 +172,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0, # Perform maximum likelihood fitting like = ctools.ctlike(obs) - like["debug"] = True # Switch this always on for results in console + like['debug'] = True # Switch this always on for results in console like.run() # Return @@ -153,12 +183,11 @@ def run_pipeline(obs, emin=0.1, emax=100.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header - print("********************************************") - print("* CTA binned analysis pipeline *") - print("********************************************") + print('********************************************') + print('* CTA binned analysis pipeline *') + print('********************************************') # Setup observations obs = setup_observations() diff --git a/examples/pipeline_stacked_disk.py b/examples/pipeline_stacked_disk.py index 3277caf..4db9d11 100755 --- a/examples/pipeline_stacked_disk.py +++ b/examples/pipeline_stacked_disk.py @@ -1,15 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script illustrates how to perform a stacked CTA analysis based on -# simulated CTA data. You may use and adapt this script to implement your -# own pipeline. +# Perform stacked analysis of simulated CTA data. # -# Usage: -# ./pipeline_stacked_disk.py -# -# ========================================================================== -# -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -33,25 +26,41 @@ from cscripts import obsutils # ================== # # Setup observations # # ================== # -def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, - emin=0.1, emax=100.0, rad=5.0, duration=1800.0, - deadc=0.95, - caldb="prod2", irf="South_50h"): +def setup_observations(pattern='four', ra=83.63, dec=22.01, offset=1.5, + emin=0.1, emax=100.0, rad=5.0, duration=180.0, + deadc=0.95, caldb='prod2', irf='South_0.5h'): """ - Returns an observation container. - - Keywords: - pattern - Pointing pattern, either "single" or "four" - ra - RA of pattern centre [deg] (default: 83.6331) - dec - DEC of pattern centre [deg] (default: 22.0145) - offset - Offset between observations of pattern [deg] (default: 1.5) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - duration - Duration of one CTA observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - caldb - Calibration database path (default: "dummy") - irf - Instrument response function (default: cta_dummy_irf) + Returns an observation container + + Parameters + ---------- + pattern : str, optional + Pointing pattern, either 'single' or 'four' + ra : float, optional + Right Ascension of pattern centre (deg) + dec : float, optional + Declination of pattern centre (deg) + offset : float, optional + Offset between observations of pattern (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + rad : float, optional + ROI radius used for analysis (deg) + duration : float, optional + Duration of one CTA observation (s) + deadc : float, optional + Deadtime correction factor + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Set list of observations obs_def_list = obsutils.set_obs_patterns(pattern, @@ -75,12 +84,21 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, # =========== # # Setup model # # =========== # -def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): +def setup_model(obs, model='data/crab.xml'): """ - Setup model for analysis. + Setup model for analysis - Keywords: - model - Model Xml file + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + model : str, optional + Model definition XML file + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Append model from file to observation container obs.models(gammalib.GModels(model)) @@ -94,126 +112,146 @@ def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): # ================================ # def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, enumbins=20, nxpix=200, nypix=200, binsz=0.02, - coordsys="CEL", proj="CAR", - model="${CTOOLS}/share/models/crab.xml", - caldb="prod2", irf="South_50h", + coordsys='CEL', proj='CAR', + model='data/crab.xml', + caldb='prod2', irf='South_0.5h', debug=False): """ - Simulation and stacked analysis pipeline. - - Keywords: - ra - RA of cube centre [deg] (default: 83.6331) - dec - DEC of cube centre [deg] (default: 22.0145) - emin - Minimum energy of cube [TeV] (default: 0.1) - emax - Maximum energy of cube [TeV] (default: 100.0) - enumbins - Number of energy bins in cube (default: 20) - nxpix - Number of RA pixels in cube (default: 200) - nypix - Number of DEC pixels in cube (default: 200) - binsz - Spatial cube bin size [deg] (default: 0.02) - coordsys - Cube coordinate system (CEL or GAL) - proj - Cube World Coordinate System (WCS) projection - debug - Enable debugging (default: False) + Simulation and stacked analysis pipeline + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + ra : float, optional + Right Ascension of counts cube centre (deg) + dec : float, optional + Declination of Region of counts cube centre (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + enumbins : int, optional + Number of energy bins + nxpix : int, optional + Number of pixels in X axis + nypix : int, optional + Number of pixels in Y axis + binsz : float, optional + Pixel size (deg) + coordsys : str, optional + Coordinate system + proj : str, optional + Coordinate projection + model : str, optional + Model definition XML file + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + debug : bool, optional + Debug function """ # Simulate events sim = ctools.ctobssim(obs) - sim["debug"] = debug - sim["outevents"] = "obs.xml" + sim['debug'] = debug + sim['outevents'] = 'obs.xml' sim.execute() # Bin events into counts map bin = ctools.ctbin() - bin["inobs"] = "obs.xml" - bin["outcube"] = "cntcube.fits" - bin["ebinalg"] = "LOG" - bin["emin"] = emin - bin["emax"] = emax - bin["enumbins"] = enumbins - bin["nxpix"] = nxpix - bin["nypix"] = nypix - bin["binsz"] = binsz - bin["coordsys"] = coordsys - bin["proj"] = proj - bin["xref"] = ra - bin["yref"] = dec - bin["debug"] = debug + bin['inobs'] = 'obs.xml' + bin['outcube'] = 'cntcube.fits' + bin['ebinalg'] = 'LOG' + bin['emin'] = emin + bin['emax'] = emax + bin['enumbins'] = enumbins + bin['nxpix'] = nxpix + bin['nypix'] = nypix + bin['binsz'] = binsz + bin['coordsys'] = coordsys + bin['proj'] = proj + bin['xref'] = ra + bin['yref'] = dec + bin['debug'] = debug bin.execute() # Create exposure cube expcube = ctools.ctexpcube() - expcube["inobs"] = "obs.xml" - expcube["incube"] = "cntcube.fits" - expcube["outcube"] = "expcube.fits" - expcube["caldb"] = caldb - expcube["irf"] = irf - expcube["ebinalg"] = "LOG" - expcube["emin"] = emin - expcube["emax"] = emax - expcube["enumbins"] = enumbins - expcube["nxpix"] = nxpix - expcube["nypix"] = nypix - expcube["binsz"] = binsz - expcube["coordsys"] = coordsys - expcube["proj"] = proj - expcube["xref"] = ra - expcube["yref"] = dec - expcube["debug"] = debug + expcube['inobs'] = 'obs.xml' + expcube['incube'] = 'cntcube.fits' + expcube['outcube'] = 'expcube.fits' + expcube['caldb'] = caldb + expcube['irf'] = irf + expcube['ebinalg'] = 'LOG' + expcube['emin'] = emin + expcube['emax'] = emax + expcube['enumbins'] = enumbins + expcube['nxpix'] = nxpix + expcube['nypix'] = nypix + expcube['binsz'] = binsz + expcube['coordsys'] = coordsys + expcube['proj'] = proj + expcube['xref'] = ra + expcube['yref'] = dec + expcube['debug'] = debug expcube.execute() # Create PSF cube psfcube = ctools.ctpsfcube() - psfcube["inobs"] = "obs.xml" - psfcube["incube"] = "NONE" - psfcube["outcube"] = "psfcube.fits" - psfcube["caldb"] = caldb - psfcube["irf"] = irf - psfcube["ebinalg"] = "LOG" - psfcube["emin"] = emin - psfcube["emax"] = emax - psfcube["enumbins"] = enumbins - psfcube["nxpix"] = 10 - psfcube["nypix"] = 10 - psfcube["binsz"] = 1.0 - psfcube["coordsys"] = coordsys - psfcube["proj"] = proj - psfcube["xref"] = ra - psfcube["yref"] = dec - psfcube["debug"] = debug + psfcube['inobs'] = 'obs.xml' + psfcube['incube'] = 'NONE' + psfcube['outcube'] = 'psfcube.fits' + psfcube['caldb'] = caldb + psfcube['irf'] = irf + psfcube['ebinalg'] = 'LOG' + psfcube['emin'] = emin + psfcube['emax'] = emax + psfcube['enumbins'] = enumbins + psfcube['nxpix'] = 10 + psfcube['nypix'] = 10 + psfcube['binsz'] = 1.0 + psfcube['coordsys'] = coordsys + psfcube['proj'] = proj + psfcube['xref'] = ra + psfcube['yref'] = dec + psfcube['debug'] = debug psfcube.execute() # Create background cube bkgcube = ctools.ctbkgcube() - bkgcube["inobs"] = "obs.xml" - bkgcube["inmodel"] = model - bkgcube["incube"] = "cntcube.fits" - bkgcube["outcube"] = "bkgcube.fits" - bkgcube["outmodel"] = "model_bkg.xml" - bkgcube["caldb"] = caldb - bkgcube["irf"] = irf - bkgcube["ebinalg"] = "LOG" - bkgcube["emin"] = emin - bkgcube["emax"] = emax - bkgcube["enumbins"] = enumbins - bkgcube["nxpix"] = 10 - bkgcube["nypix"] = 10 - bkgcube["binsz"] = 1.0 - bkgcube["coordsys"] = coordsys - bkgcube["proj"] = proj - bkgcube["xref"] = ra - bkgcube["yref"] = dec - bkgcube["debug"] = debug + bkgcube['inobs'] = 'obs.xml' + bkgcube['inmodel'] = model + bkgcube['incube'] = 'NONE' + bkgcube['outcube'] = 'bkgcube.fits' + bkgcube['outmodel'] = 'model_bkg.xml' + bkgcube['caldb'] = caldb + bkgcube['irf'] = irf + bkgcube['ebinalg'] = 'LOG' + bkgcube['emin'] = emin + bkgcube['emax'] = emax + bkgcube['enumbins'] = enumbins + bkgcube['nxpix'] = 10 + bkgcube['nypix'] = 10 + bkgcube['binsz'] = 1.0 + bkgcube['coordsys'] = coordsys + bkgcube['proj'] = proj + bkgcube['xref'] = ra + bkgcube['yref'] = dec + bkgcube['debug'] = debug bkgcube.execute() # Perform maximum likelihood fitting like = ctools.ctlike() - like["inobs"] = "cntcube.fits" - like["inmodel"] = "model_bkg.xml" - like["outmodel"] = "fit_results.xml" - like["expcube"] = "expcube.fits" - like["psfcube"] = "psfcube.fits" - like["bkgcube"] = "bkgcube.fits" - like["caldb"] = caldb - like["irf"] = irf - like["debug"] = True # Switch this always on for results in console + like['inobs'] = 'cntcube.fits' + like['inmodel'] = 'model_bkg.xml' + like['outmodel'] = 'fit_results.xml' + like['expcube'] = 'expcube.fits' + like['psfcube'] = 'psfcube.fits' + like['bkgcube'] = 'bkgcube.fits' + like['caldb'] = caldb + like['irf'] = irf + like['debug'] = True # Switch this always on for results in console like.execute() # Return @@ -224,12 +262,11 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header - print("********************************************") - print("* CTA stacked analysis pipeline *") - print("********************************************") + print('********************************************') + print('* CTA stacked analysis pipeline *') + print('********************************************') # Setup observations obs = setup_observations() diff --git a/examples/pipeline_stacked_mem.py b/examples/pipeline_stacked_mem.py index 555055a..a0c6986 100755 --- a/examples/pipeline_stacked_mem.py +++ b/examples/pipeline_stacked_mem.py @@ -1,15 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script illustrates how to perform a stacked CTA analysis based on -# simulated CTA data. You may use and adapt this script to implement your -# own pipeline. +# Perform stacked in-memory analysis of simulated CTA data # -# Usage: -# ./pipeline_stacked_mem.py -# -# ========================================================================== -# -# Copyright (C) 2014-2015 Juergen Knoedlseder +# Copyright (C) 2014-2016 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 @@ -33,25 +26,41 @@ from cscripts import obsutils # ================== # # Setup observations # # ================== # -def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, - emin=0.1, emax=100.0, rad=5.0, duration=1800.0, - deadc=0.95, - caldb="prod2", irf="South_50h"): +def setup_observations(pattern='four', ra=83.63, dec=22.01, offset=1.5, + emin=0.1, emax=100.0, rad=5.0, duration=180.0, + deadc=0.95, caldb='prod2', irf='South_0.5h'): """ - Returns an observation container. - - Keywords: - pattern - Pointing pattern, either "single" or "four" - ra - RA of pattern centre [deg] (default: 83.6331) - dec - DEC of pattern centre [deg] (default: 22.0145) - offset - Offset between observations of pattern [deg] (default: 1.5) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - duration - Duration of one CTA observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - caldb - Calibration database path (default: "dummy") - irf - Instrument response function (default: cta_dummy_irf) + Returns an observation container + + Parameters + ---------- + pattern : str, optional + Pointing pattern, either 'single' or 'four' + ra : float, optional + Right Ascension of pattern centre (deg) + dec : float, optional + Declination of pattern centre (deg) + offset : float, optional + Offset between observations of pattern (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + rad : float, optional + ROI radius used for analysis (deg) + duration : float, optional + Duration of one CTA observation (s) + deadc : float, optional + Deadtime correction factor + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Set list of observations obs_def_list = obsutils.set_obs_patterns(pattern, @@ -75,12 +84,21 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, # =========== # # Setup model # # =========== # -def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): +def setup_model(obs, model='data/crab.xml'): """ - Setup model for analysis. - - Keywords: - model - Model Xml file + Setup model for analysis + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + model : str, optional + Model definition XML file + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Append model from file to observation container obs.models(gammalib.GModels(model)) @@ -94,93 +112,107 @@ def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): # ================================ # def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, enumbins=20, nxpix=200, nypix=200, binsz=0.02, - coordsys="CEL", proj="CAR", debug=False): + coordsys='CEL', proj='CAR', debug=False): """ - Simulation and stacked analysis pipeline. - - Keywords: - ra - RA of cube centre [deg] (default: 83.6331) - dec - DEC of cube centre [deg] (default: 22.0145) - emin - Minimum energy of cube [TeV] (default: 0.1) - emax - Maximum energy of cube [TeV] (default: 100.0) - enumbins - Number of energy bins in cube (default: 20) - nxpix - Number of RA pixels in cube (default: 200) - nypix - Number of DEC pixels in cube (default: 200) - binsz - Spatial cube bin size [deg] (default: 0.02) - coordsys - Cube coordinate system (CEL or GAL) - proj - Cube World Coordinate System (WCS) projection - debug - Enable debugging (default: False) + Simulation and stacked analysis pipeline + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + ra : float, optional + Right Ascension of counts cube centre (deg) + dec : float, optional + Declination of Region of counts cube centre (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + enumbins : int, optional + Number of energy bins + nxpix : int, optional + Number of pixels in X axis + nypix : int, optional + Number of pixels in Y axis + binsz : float, optional + Pixel size (deg) + coordsys : str, optional + Coordinate system + proj : str, optional + Coordinate projection + debug : bool, optional + Debug function """ # Simulate events sim = ctools.ctobssim(obs) - sim["debug"] = debug + sim['debug'] = debug sim.run() # Bin events into counts map bin = ctools.ctbin(sim.obs()) - bin["ebinalg"] = "LOG" - bin["emin"] = emin - bin["emax"] = emax - bin["enumbins"] = enumbins - bin["nxpix"] = nxpix - bin["nypix"] = nypix - bin["binsz"] = binsz - bin["coordsys"] = coordsys - bin["proj"] = proj - bin["xref"] = ra - bin["yref"] = dec - bin["debug"] = debug + bin['ebinalg'] = 'LOG' + bin['emin'] = emin + bin['emax'] = emax + bin['enumbins'] = enumbins + bin['nxpix'] = nxpix + bin['nypix'] = nypix + bin['binsz'] = binsz + bin['coordsys'] = coordsys + bin['proj'] = proj + bin['xref'] = ra + bin['yref'] = dec + bin['debug'] = debug bin.run() # Create exposure cube expcube = ctools.ctexpcube(sim.obs()) - expcube["incube"] = "NONE" - expcube["ebinalg"] = "LOG" - expcube["emin"] = emin - expcube["emax"] = emax - expcube["enumbins"] = enumbins - expcube["nxpix"] = nxpix - expcube["nypix"] = nypix - expcube["binsz"] = binsz - expcube["coordsys"] = coordsys - expcube["proj"] = proj - expcube["xref"] = ra - expcube["yref"] = dec - expcube["debug"] = debug + expcube['incube'] = 'NONE' + expcube['ebinalg'] = 'LOG' + expcube['emin'] = emin + expcube['emax'] = emax + expcube['enumbins'] = enumbins + expcube['nxpix'] = nxpix + expcube['nypix'] = nypix + expcube['binsz'] = binsz + expcube['coordsys'] = coordsys + expcube['proj'] = proj + expcube['xref'] = ra + expcube['yref'] = dec + expcube['debug'] = debug expcube.run() # Create PSF cube psfcube = ctools.ctpsfcube(sim.obs()) - psfcube["incube"] = "NONE" - psfcube["ebinalg"] = "LOG" - psfcube["emin"] = emin - psfcube["emax"] = emax - psfcube["enumbins"] = enumbins - psfcube["nxpix"] = 10 - psfcube["nypix"] = 10 - psfcube["binsz"] = 1.0 - psfcube["coordsys"] = coordsys - psfcube["proj"] = proj - psfcube["xref"] = ra - psfcube["yref"] = dec - psfcube["debug"] = debug + psfcube['incube'] = 'NONE' + psfcube['ebinalg'] = 'LOG' + psfcube['emin'] = emin + psfcube['emax'] = emax + psfcube['enumbins'] = enumbins + psfcube['nxpix'] = 10 + psfcube['nypix'] = 10 + psfcube['binsz'] = 1.0 + psfcube['coordsys'] = coordsys + psfcube['proj'] = proj + psfcube['xref'] = ra + psfcube['yref'] = dec + psfcube['debug'] = debug psfcube.run() # Create background cube bkgcube = ctools.ctbkgcube(sim.obs()) - bkgcube["incube"] = "NONE" - bkgcube["ebinalg"] = "LOG" - bkgcube["emin"] = emin - bkgcube["emax"] = emax - bkgcube["enumbins"] = enumbins - bkgcube["nxpix"] = 10 - bkgcube["nypix"] = 10 - bkgcube["binsz"] = 1.0 - bkgcube["coordsys"] = coordsys - bkgcube["proj"] = proj - bkgcube["xref"] = ra - bkgcube["yref"] = dec - bkgcube["debug"] = debug + bkgcube['incube'] = 'NONE' + bkgcube['ebinalg'] = 'LOG' + bkgcube['emin'] = emin + bkgcube['emax'] = emax + bkgcube['enumbins'] = enumbins + bkgcube['nxpix'] = 10 + bkgcube['nypix'] = 10 + bkgcube['binsz'] = 1.0 + bkgcube['coordsys'] = coordsys + bkgcube['proj'] = proj + bkgcube['xref'] = ra + bkgcube['yref'] = dec + bkgcube['debug'] = debug bkgcube.run() # Attach background model to observation container @@ -192,7 +224,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, # Perform maximum likelihood fitting like = ctools.ctlike(bin.obs()) - like["debug"] = True # Switch this always on for results in console + like['debug'] = True # Switch this always on for results in console like.run() # Return @@ -203,12 +235,11 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header - print("********************************************") - print("* CTA stacked analysis pipeline *") - print("********************************************") + print('********************************************') + print('* CTA stacked analysis pipeline *') + print('********************************************') # Setup observations obs = setup_observations() diff --git a/examples/pipeline_unbinned_disk.py b/examples/pipeline_unbinned_disk.py index 55fae66..1f93093 100755 --- a/examples/pipeline_unbinned_disk.py +++ b/examples/pipeline_unbinned_disk.py @@ -1,15 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script illustrates how to perform an unbinned CTA analysis based on -# simulated CTA data. You may use and adapt this script to implement your -# own pipeline. +# Perform unbinned analysis of simulated CTA data. # -# Usage: -# ./pipeline_unbinned_disk.py -# -# ========================================================================== -# -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -33,25 +26,41 @@ from cscripts import obsutils # ================== # # Setup observations # # ================== # -def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, - emin=0.1, emax=100.0, rad=5.0, duration=1800.0, - deadc=0.95, - caldb="prod2", irf="South_50h"): +def setup_observations(pattern='four', ra=83.63, dec=22.01, offset=1.5, + emin=0.1, emax=100.0, rad=5.0, duration=180.0, + deadc=0.95, caldb='prod2', irf='South_0.5h'): """ - Returns an observation container. - - Keywords: - pattern - Pointing pattern, either "single" or "four" - ra - RA of pattern centre [deg] (default: 83.6331) - dec - DEC of pattern centre [deg] (default: 22.0145) - offset - Offset between observations of pattern [deg] (default: 1.5) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - duration - Duration of one CTA observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - caldb - Calibration database path (default: "dummy") - irf - Instrument response function (default: cta_dummy_irf) + Returns an observation container + + Parameters + ---------- + pattern : str, optional + Pointing pattern, either 'single' or 'four' + ra : float, optional + Right Ascension of pattern centre (deg) + dec : float, optional + Declination of pattern centre (deg) + offset : float, optional + Offset between observations of pattern (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + rad : float, optional + ROI radius used for analysis (deg) + duration : float, optional + Duration of one CTA observation (s) + deadc : float, optional + Deadtime correction factor + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Set list of observations obs_def_list = obsutils.set_obs_patterns(pattern, @@ -75,12 +84,21 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, # =========== # # Setup model # # =========== # -def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): +def setup_model(obs, model='data/crab.xml'): """ - Setup model for analysis. - - Keywords: - model - Model Xml file + Setup model for analysis + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + model : str, optional + Model definition XML file + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Append model from file to observation container obs.models(gammalib.GModels(model)) @@ -95,55 +113,69 @@ def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, emin=0.1, emax=100.0, tmin=0.0, tmax=0.0, - model="${CTOOLS}/share/models/crab.xml", - caldb="prod2", irf="South_50h", + model='data/crab.xml', + caldb='prod2', irf='South_0.5h', debug=False): """ - Simulation and unbinned analysis pipeline. - - Keywords: - ra - RA of cube centre [deg] (default: 83.63) - dec - DEC of cube centre [deg] (default: 22.01) - rad - Selection radius [deg] (default: 3.0) - emin - Minimum energy of cube [TeV] (default: 0.1) - emax - Maximum energy of cube [TeV] (default: 100.0) - tmin - Start time [MET] (default: 0.0) - tmax - Stop time [MET] (default: 0.0) - model - Model Xml file - caldb - Calibration database path (default: "dummy") - irf - Instrument response function (default: cta_dummy_irf) - debug - Enable debugging (default: False) + Simulation and unbinned analysis pipeline + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + ra : float, optional + Right Ascension of Region of Interest centre (deg) + dec : float, optional + Declination of Region of Interest centre (deg) + rad : float, optional + Radius of Region of Interest (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + tmin : float, optional + Start time (s) + tmax : float, optional + Stop time (s) + model : str, optional + Model definition XML file + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + debug : bool, optional + Debug function """ # Get model # Simulate events sim = ctools.ctobssim(obs) - sim["debug"] = debug - sim["outevents"] = "obs.xml" + sim['debug'] = debug + sim['outevents'] = 'obs.xml' sim.execute() # Select events select = ctools.ctselect() - select["inobs"] = "obs.xml" - select["outobs"] = "obs_selected.xml" - select["ra"] = ra - select["dec"] = dec - select["rad"] = rad - select["emin"] = emin - select["emax"] = emax - select["tmin"] = tmin - select["tmax"] = tmax - select["debug"] = debug + select['inobs'] = 'obs.xml' + select['outobs'] = 'obs_selected.xml' + select['ra'] = ra + select['dec'] = dec + select['rad'] = rad + select['emin'] = emin + select['emax'] = emax + select['tmin'] = tmin + select['tmax'] = tmax + select['debug'] = debug select.execute() # Perform maximum likelihood fitting like = ctools.ctlike() - like["inobs"] = "obs_selected.xml" - like["inmodel"] = model - like["outmodel"] = "fit_results.xml" - like["caldb"] = caldb - like["irf"] = irf - like["debug"] = True # Switch this always on for results in console + like['inobs'] = 'obs_selected.xml' + like['inmodel'] = model + like['outmodel'] = 'fit_results.xml' + like['caldb'] = caldb + like['irf'] = irf + like['debug'] = True # Switch this always on for results in console like.execute() # Return @@ -154,12 +186,11 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header - print("********************************************") - print("* CTA unbinned analysis pipeline *") - print("********************************************") + print('********************************************') + print('* CTA unbinned analysis pipeline *') + print('********************************************') # Setup observations obs = setup_observations() diff --git a/examples/pipeline_unbinned_mem.py b/examples/pipeline_unbinned_mem.py index 392c21c..925a515 100755 --- a/examples/pipeline_unbinned_mem.py +++ b/examples/pipeline_unbinned_mem.py @@ -1,15 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script illustrates how to perform an unbinned CTA analysis based on -# simulated CTA data. You may use and adapt this script to implement your -# own pipeline. +# Perform unbinned in-memory analysis of simulated CTA data # -# Usage: -# ./pipeline_unbinned_mem.py -# -# ========================================================================== -# -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -33,25 +26,41 @@ from cscripts import obsutils # ================== # # Setup observations # # ================== # -def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, - emin=0.1, emax=100.0, rad=5.0, duration=1800.0, - deadc=0.95, - caldb="prod2", irf="South_50h"): +def setup_observations(pattern='four', ra=83.63, dec=22.01, offset=1.5, + emin=0.1, emax=100.0, rad=5.0, duration=180.0, + deadc=0.95, caldb='prod2', irf='South_0.5h'): """ - Returns an observation container. - - Keywords: - pattern - Pointing pattern, either "single" or "four" - ra - RA of pattern centre [deg] (default: 83.6331) - dec - DEC of pattern centre [deg] (default: 22.0145) - offset - Offset between observations of pattern [deg] (default: 1.5) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - duration - Duration of one CTA observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - caldb - Calibration database path (default: "dummy") - irf - Instrument response function (default: cta_dummy_irf) + Returns an observation container + + Parameters + ---------- + pattern : str, optional + Pointing pattern, either 'single' or 'four' + ra : float, optional + Right Ascension of pattern centre (deg) + dec : float, optional + Declination of pattern centre (deg) + offset : float, optional + Offset between observations of pattern (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + rad : float, optional + ROI radius used for analysis (deg) + duration : float, optional + Duration of one CTA observation (s) + deadc : float, optional + Deadtime correction factor + caldb : str, optional + Calibration database path + irf : str, optional + Instrument response function + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Set list of observations obs_def_list = obsutils.set_obs_patterns(pattern, @@ -75,12 +84,21 @@ def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, # =========== # # Setup model # # =========== # -def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"): +def setup_model(obs, model='data/crab.xml'): """ - Setup model for analysis. - - Keywords: - model - Model Xml file + Setup model for analysis + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + model : str, optional + Model definition XML file + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ # Append model from file to observation container obs.models(gammalib.GModels(model)) @@ -97,38 +115,49 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, tmin=0.0, tmax=0.0, debug=False): """ - Simulation and unbinned analysis pipeline. - - Keywords: - ra - RA of cube centre [deg] (default: 83.63) - dec - DEC of cube centre [deg] (default: 22.01) - rad - Selection radius [deg] (default: 3.0) - emin - Minimum energy of cube [TeV] (default: 0.1) - emax - Maximum energy of cube [TeV] (default: 100.0) - tmin - Start time [MET] (default: 0.0) - tmax - Stop time [MET] (default: 0.0) - debug - Enable debugging (default: False) + Simulation and binned analysis pipeline + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + ra : float, optional + Right Ascension of Region of Interest centre (deg) + dec : float, optional + Declination of Region of Interest centre (deg) + rad : float, optional + Radius of Region of Interest (deg) + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + tmin : float, optional + Start time (s) + tmax : float, optional + Stop time (s) + debug : bool, optional + Debug function """ # Simulate events sim = ctools.ctobssim(obs) - sim["debug"] = debug + sim['debug'] = debug sim.run() # Select events select = ctools.ctselect(sim.obs()) - select["ra"] = ra - select["dec"] = dec - select["rad"] = rad - select["emin"] = emin - select["emax"] = emax - select["tmin"] = tmin - select["tmax"] = tmax - select["debug"] = debug + select['ra'] = ra + select['dec'] = dec + select['rad'] = rad + select['emin'] = emin + select['emax'] = emax + select['tmin'] = tmin + select['tmax'] = tmax + select['debug'] = debug select.run() # Perform maximum likelihood fitting like = ctools.ctlike(select.obs()) - like["debug"] = True # Switch this always on for results in console + like['debug'] = True # Switch this always on for results in console like.run() # Return @@ -139,12 +168,11 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header - print("********************************************") - print("* CTA unbinned analysis pipeline *") - print("********************************************") + print('********************************************') + print('* CTA unbinned analysis pipeline *') + print('********************************************') # Setup observations obs = setup_observations() diff --git a/examples/show_butterfly.py b/examples/show_butterfly.py index 1916d62..a3539e7 100755 --- a/examples/show_butterfly.py +++ b/examples/show_butterfly.py @@ -1,8 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script shows how to plot a butterfly created with ctbutterfly +# Shows butterfly diagram created with ctbutterfly # -# Copyright (C) 2014-2015 Michael Mayer +# Copyright (C) 2014-2016 Michael Mayer # # 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 @@ -18,23 +18,30 @@ # along with this program. If not, see . # # ========================================================================== -import gammalib import sys +import gammalib try: import matplotlib.pyplot as plt + plt.figure() + plt.close() except: - sys.exit("This script needs matplotlib") + print('This script needs the "matplotlib" module') + sys.exit() # ============ # # Script entry # # ============ # -if __name__ == "__main__": - """ - """ +if __name__ == '__main__': + # Check for given butterfly file - if not len(sys.argv) == 2: - sys.exit("Usage: show_butterfly.py butterfly.txt") + if len(sys.argv) < 2: + sys.exit('Usage: show_butterfly.py butterfly.txt [file]') + + # Check if plotting in file is requested + plotfile = '' + if len(sys.argv) == 3: + plotfile = sys.argv[2] # Read given butterfly file filename = sys.argv[1] @@ -75,6 +82,11 @@ if __name__ == "__main__": plt.grid() plt.plot(line_x,line_y,color='black',ls='-') plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5) - plt.xlabel("Energy (MeV)") - plt.ylabel(r"E dN/dE (MeV$^{-1}$ s$^{-1}$ cm$^{-2}$)") - plt.show() + plt.xlabel('Energy (MeV)') + plt.ylabel(r'E $\times$ dN/dE (MeV$^{-1}$ s$^{-1}$ cm$^{-2}$)') + + # Show spectrum or save it into file + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() diff --git a/examples/show_model.py b/examples/show_model.py index 4fee000..a160016 100755 --- a/examples/show_model.py +++ b/examples/show_model.py @@ -1,8 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script shows a model +# Shows the photon spectrum of a model # -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -18,27 +18,44 @@ # along with this program. If not, see . # # ========================================================================== -import gammalib import sys import math +import gammalib try: import matplotlib.pyplot as plt + plt.figure() + plt.close() except: - sys.exit("This script needs matplotlib") + print('This script needs the "matplotlib" module') + sys.exit() # ======================= # # Plot spectral component # # ======================= # -def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100): +def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100, plotfile=''): """ + Plot spectral model component + + Parameters + ---------- + model : `~gammalib.GModel` + Model + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + enumbins : integer, optional + Number of energy bins + plotfile : str, optional + Name of plot file """ # Get spectral component spectrum = model.spectral() # Setup energy axis - e_min = gammalib.GEnergy(emin, "TeV") - e_max = gammalib.GEnergy(emax, "TeV") + e_min = gammalib.GEnergy(emin, 'TeV') + e_max = gammalib.GEnergy(emax, 'TeV') ebounds = gammalib.GEbounds(enumbins, e_min, e_max) # Setup model plot @@ -48,7 +65,7 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100): max = 0.0 for i in range(enumbins): energy = ebounds.elogmean(i) - value = spectrum.eval(energy, gammalib.GTime()) + value = spectrum.eval(energy) if value > max: max = value if value < min: @@ -57,22 +74,27 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100): y.append(value) # Set minimum and maximum - max = math.pow(10.0, math.log10(max)+1) - if min < 0.0001*max: - min = 0.0001*max - min = math.pow(10.0, math.log10(min)-1) + #max = math.pow(10.0, math.log10(max)+1) + #if min < 0.0001*max: + # min = 0.0001*max + #min = math.pow(10.0, math.log10(min)-1) # Show spectrum plt.figure() - plt.title(model.name()+" ("+spectrum.type()+")") + plt.title(model.name()+' ('+spectrum.type()+')') plt.loglog() plt.grid() plt.loglog(x, y, color='red') - plt.xlabel("Energy [TeV]") - plt.ylabel(r"dN/dE [ph s$^{-1}$ cm$^{-2}$ MeV$^{-1}$]") - plt.ylim([min,max]) - plt.show() + plt.xlabel('Energy (TeV)') + plt.ylabel(r'dN/dE (ph s$^{-1}$ cm$^{-2}$ MeV$^{-1}$)') + #plt.ylim([min,max]) + + # Show spectrum or save it into file + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() # Return return @@ -81,25 +103,28 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100): # ============= # # Script entry # # ============= # -if __name__ == "__main__": - """ - """ +if __name__ == '__main__': + # Check for model filename if len(sys.argv) < 2: - sys.exit("Usage: show_model.py model.xml [index]") + sys.exit('Usage: show_model.py model.xml [name] [file]') # Get optional model index - if len(sys.argv) == 3: - index = int(sys.argv[2]) + if len(sys.argv) >= 3: + name = sys.argv[2] + else: + name = 0 + if len(sys.argv) == 4: + plotfile = sys.argv[3] else: - index = 0 + plotfile = '' # Read models XML file filename = sys.argv[1] models = gammalib.GModels(filename) # Extract relevant model - model = models[index] + model = models[name] # Plot spectrum - plot_spectrum(model) + plot_spectrum(model, plotfile=plotfile) diff --git a/examples/show_pha.py b/examples/show_pha.py index 547752f..904ea9c 100755 --- a/examples/show_pha.py +++ b/examples/show_pha.py @@ -1,12 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays a spectrum in PHA format. +# Show a spectrum in PHA format. # -# Required 3rd party modules: -# - matplotlib -# - numpy -# -# Copyright (C) 2013 Juergen Knoedlseder +# Copyright (C) 2013-2016 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 @@ -25,17 +21,32 @@ import sys import math import gammalib +try: + import matplotlib.pyplot as plt + plt.figure() + plt.close() +except: + print('This script needs the "matplotlib" module') + sys.exit() # ================= # # Show PHA spectrum # # ================= # -def show_pha(pha): +def show_pha(pha, plotfile=''): """ + Show PHA spectrum + + Parameters + ---------- + spectrum : `~gammalib.GPha` + PHA file + plotfile : str, optional + Plot file name """ # Create figure plt.figure(1) - plt.title("PHA spectrum (" + pha.filename() + ")") + plt.title('PHA spectrum (' + pha.filename() + ')') # Generate energy vector energy = [] @@ -59,51 +70,46 @@ def show_pha(pha): plt.semilogy(energy, counts, 'ro') else: plt.loglog(energy, counts, 'ro') - plt.errorbar(energy, counts, error, fmt=None, ecolor='r') + plt.errorbar(energy, counts, error, fmt='ro', ecolor='r') # Set axes if channels: - plt.xlabel("Channels") + plt.xlabel('Channels') else: - plt.xlabel("Energy (TeV)") - plt.ylabel("Counts") + plt.xlabel('Energy (TeV)') + plt.ylabel('Counts') - # Show plot - plt.show() + # Show spectrum or save it into file + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() # Return - #return + return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Display spectrum in PHA format. - """ + # Print usage information - usage = "Usage: show_pha filename" \ - " [-n bins] [-c column] [-t title] [-p plot]" + usage = 'Usage: show_pha.py filename [file]' if len(sys.argv) < 2: print(usage) sys.exit() + # Check if plotting in file is requested + plotfile = '' + if len(sys.argv) == 3: + plotfile = sys.argv[2] + # Extract parameters filename = sys.argv[1] # Load PHA spectrum pha = gammalib.GPha(filename) - #print(pha) - - # Try importing matplotlib - try: - import matplotlib.pyplot as plt - has_matplotlib = True - except ImportError: - print("Matplotlib is not (correctly) installed on your system.") - has_matplotlib = False # Show PHA spectrum - if has_matplotlib: - show_pha(pha) + show_pha(pha, plotfile=plotfile) diff --git a/examples/show_pull_evolution.py b/examples/show_pull_evolution.py index 8e6b3bb..06bace2 100755 --- a/examples/show_pull_evolution.py +++ b/examples/show_pull_evolution.py @@ -1,19 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays the evolution of the pull distributions as function -# of the number of samples. The pull distributions are generated using the -# cspull script. The pull is defined by +# Shows the evolution of the mean and rms pull # -# (fitted value - simulated value) / estimated error -# -# and the pull histogram allows to assess possible biases in the fit -# parameters and errors. -# -# Required 3rd party modules: -# - matplotlib -# - numpy -# -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -29,10 +18,20 @@ # along with this program. If not, see . # # ========================================================================== -import matplotlib.pyplot as plt -import numpy as np import sys import csv +try: + import matplotlib.pyplot as plt + plt.figure() + plt.close() +except: + print('This script needs the "matplotlib" module') + sys.exit() +try: + import numpy as np +except: + print('This script needs the "numpy" module') + sys.exit() # ====================== # @@ -40,7 +39,14 @@ import csv # ====================== # def read_pull(filename, parname): """ - Read pull distribution. + Read pull distribution + + Parameters + ---------- + filename : str + Pull distribution ASCII file + parname : str + Parameter """ # Initialise list values = [] @@ -58,9 +64,9 @@ def read_pull(filename, parname): try: index = row.index(parname) except: - sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n') + print('ERROR: Parameter "'+parname+'" not found in list:') for p in row: - sys.stdout.write(' "'+p+'"\n') + print(' "'+p+'"') raise NameError(parname) # Handle data rows @@ -81,15 +87,18 @@ def read_pull(filename, parname): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Display pull evolution generated by cspull using matplotlib. - """ + # Print usage information - usage = "Usage: show_pull_evolution filename parname" + usage = 'Usage: show_pull_evolution.py filename parname [file]' if len(sys.argv) < 3: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() + # Check if plotting in file is requested + plotfile = '' + if len(sys.argv) == 4: + plotfile = sys.argv[3] + # Extract parameters filename = sys.argv[1] parname = sys.argv[2] @@ -108,17 +117,20 @@ if __name__ == '__main__': rms.append(np.std(values[0:i])) # Plot mean and rms evolution - plt.semilogx(samples, mean, 'r-', label="Mean") - plt.semilogx(samples, rms, 'b-', label="Std. deviation") + plt.semilogx(samples, mean, 'r-', label='Mean') + plt.semilogx(samples, rms, 'b-', label='Std. deviation') # Set plot - plt.xlabel("Samples") - plt.ylabel("Value") + plt.xlabel('Samples') + plt.ylabel('Value') plt.title(parname) plt.grid(True) # Set legend - plt.legend(loc="upper left") + plt.legend(loc='upper left') - # Show histogram - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() diff --git a/examples/show_pull_histogram.py b/examples/show_pull_histogram.py index 3cf3707..7a8521e 100755 --- a/examples/show_pull_histogram.py +++ b/examples/show_pull_histogram.py @@ -1,18 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays pull histograms that have been generated using the -# cspull script. The pull is defined by +# Shows the pull histogram # -# (fitted value - simulated value) / estimated error -# -# and the pull histogram allows to assess possible biases in the fit -# parameters and errors. -# -# Required 3rd party modules: -# - matplotlib -# - numpy -# -# Copyright (C) 2011 Jurgen Knodlseder +# Copyright (C) 2011-2016 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 @@ -28,11 +18,21 @@ # along with this program. If not, see . # # ========================================================================== -import matplotlib.mlab as mlab -import matplotlib.pyplot as plt -import numpy as np import sys import csv +try: + import matplotlib.pyplot as plt + import matplotlib.mlab as mlab + plt.figure() + plt.close() +except: + print('This script needs the "matplotlib" module') + sys.exit() +try: + import numpy as np +except: + print('This script needs the "numpy" module') + sys.exit() # ====================== # @@ -40,7 +40,14 @@ import csv # ====================== # def read_pull(filename, parname): """ - Read pull distribution. + Read pull distribution + + Parameters + ---------- + filename : str + Pull distribution ASCII file + parname : str + Parameter """ # Initialise list values = [] @@ -58,9 +65,9 @@ def read_pull(filename, parname): try: index = row.index(parname) except: - sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n') + print('ERROR: Parameter "'+parname+'" not found in list:') for p in row: - sys.stdout.write(' "'+p+'"\n') + print(' "'+p+'"') raise NameError(parname) # Handle data rows @@ -81,13 +88,11 @@ def read_pull(filename, parname): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Display pull histogram generated by cspull using matplotlib. - """ + # Print usage information - usage = "Usage: show_pull_histogram filename parname [bins]" + usage = 'Usage: show_pull_histogram.py filename parname [bins] [file]' if len(sys.argv) < 3: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() # Extract parameters @@ -97,13 +102,17 @@ if __name__ == '__main__': nbins = int(sys.argv[3]) else: nbins = 50 + if len(sys.argv) == 5: + plotfile = sys.argv[4] + else: + plotfile = '' # Read values from CSV file values = read_pull(filename, parname) # Create histogram n, bins, patches = plt.hist(values, nbins, range=[-4.0,4.0], - normed=True, facecolor='green') + normed=True, facecolor='green') # Create expected distribution y = mlab.normpdf(bins, 0.0, 1.0) @@ -118,5 +127,8 @@ if __name__ == '__main__': plt.title(name) plt.grid(True) - # Show histogram - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() diff --git a/examples/show_response.py b/examples/show_response.py index 3ba9414..468f736 100755 --- a/examples/show_response.py +++ b/examples/show_response.py @@ -1,12 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays one or several CTA response functions. +# Shows one or several CTA response functions. # -# Required 3rd party modules: -# - matplotlib -# - numpy -# -# Copyright (C) 2014-2015 Juergen Knoedlseder +# Copyright (C) 2014-2016 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 @@ -25,7 +21,13 @@ import sys import math import gammalib -import matplotlib.pyplot as plt +try: + import matplotlib.pyplot as plt + plt.figure() + plt.close() +except: + print('This script needs the "matplotlib" module') + sys.exit() # ================================== # @@ -35,11 +37,19 @@ def sigma_lima(Non, Noff, alpha=0.2): """ Compute Eq. (17) of Li & Ma (1983). - Parameters: - Non - Number of on counts - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure + Parameters + ---------- + Non : float + Number of On counts + Noff : float + Number of Off counts + alpha : float, optional + Ratio of on-to-off exposure + + Returns + ------- + sigma : float + Detection significance in Gaussian sigma """ # Compute sensitivity alpha1 = alpha + 1.0 @@ -61,11 +71,19 @@ def Non_lima(sigma, Noff, alpha=0.2): """ Solve Eq. (17) of Li & Ma (1983) for Non - Parameters: - sigma - Required significance - Noff - Number of off counts - Keywords: - alpha - Ratio of on-to-off exposure + Parameters + ---------- + sigma : float + Required detection significance in Gaussian sigma + Noff : float + Number of Off counts + alpha : float, optional + Ratio of on-to-off exposure + + Returns + ------- + Non : float + Number of On counts """ # Set initial guess for Non S = 25.0 @@ -89,9 +107,18 @@ def Non_lima(sigma, Noff, alpha=0.2): # ======================= # # Show one effective area # # ======================= # -def show_one_effective_area(rsp, name, color="r"): +def show_one_effective_area(rsp, name, color='r'): """ - Show one effective area. + Show one effective area + + Parameters + ---------- + rsp : `~gammalib.GCTAResponse` + Response function + name : str + Name of the response function + color : str, optional + Color for plot """ # Set to figure 1 plt.figure(1) @@ -111,11 +138,11 @@ def show_one_effective_area(rsp, name, color="r"): plt.loglog(E, aeff, color+'o') # Set axes - plt.xlabel("Energy [TeV]") - plt.ylabel("Effective area [cm2]") + plt.xlabel('Energy (TeV)') + plt.ylabel(r'Effective area (cm$^{2}$)') # Set legend - plt.legend(loc="lower right") + plt.legend(loc='lower right') # Return return @@ -124,9 +151,18 @@ def show_one_effective_area(rsp, name, color="r"): # ======================== # # Show one background rate # # ======================== # -def show_one_background_rate(rsp, name, color="r"): +def show_one_background_rate(rsp, name, color='r'): """ - Show one background rate. + Show one background rate + + Parameters + ---------- + rsp : `~gammalib.GCTAResponse` + Response function + name : str + Name of the response function + color : str, optional + Color for plot """ # Set to figure 2 plt.figure(2) @@ -148,11 +184,11 @@ def show_one_background_rate(rsp, name, color="r"): plt.loglog(E, bgrate, color+'o') # Set axes - plt.xlabel("Energy [TeV]") - plt.ylabel("Background rate [events/s/MeV/sr]") + plt.xlabel('Energy (TeV)') + plt.ylabel(r'Background rate (counts s$^{-1}$ MeV$^{-1}$ sr$^{-1}$)') # Set legend - plt.legend(loc="lower left") + plt.legend(loc='lower left') # Return return @@ -161,9 +197,25 @@ def show_one_background_rate(rsp, name, color="r"): # ==================== # # Show one sensitivity # # ==================== # -def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sigma=5.0): +def show_one_sensitivity(rsp, name, color='r', duration=180000.0, alpha=0.2, + sigma=5.0): """ - Show one sensitivity. + Show one sensitivity + + Parameters + ---------- + rsp : `~gammalib.GCTAResponse` + Response function + name : str + Name of the response function + color : str, optional + Color for plot + duration : float, optional + Duration of observation (s) + alpha : float, optional + Ratio of on-to-off exposure + sigma : float, optional + Required detection significance in Gaussian sigma """ # Set constants r68_to_sigma = 0.6624305 @@ -194,9 +246,9 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig src_counts = 0.05*bgd_counts if src_counts < 10: src_counts = 10.0 - emin = gammalib.GEnergy(math.pow(10.0, logE[i]-0.1), "TeV") - emax = gammalib.GEnergy(math.pow(10.0, logE[i]+0.1), "TeV") - epivot = gammalib.GEnergy(math.pow(10.0, logE[i]), "TeV") + emin = gammalib.GEnergy(math.pow(10.0, logE[i]-0.1), 'TeV') + emax = gammalib.GEnergy(math.pow(10.0, logE[i]+0.1), 'TeV') + epivot = gammalib.GEnergy(math.pow(10.0, logE[i]), 'TeV') plaw = gammalib.GModelSpectralPlaw(1.0e-6, -2.6, epivot) conv = TeV2erg*E[i]*E[i]/plaw.flux(emin, emax) flux[i] = conv * src_counts / (duration * aeff*0.68) @@ -206,11 +258,11 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig plt.loglog(E, flux, color+'o') # Set axes - plt.xlabel("Energy [TeV]") - plt.ylabel("Sensitivity [erg/cm2/s]") + plt.xlabel('Energy (TeV)') + plt.ylabel(r'Sensitivity (erg cm$^{-2}$ s$^{-1}$)') # Set legend - plt.legend(loc="upper right") + plt.legend(loc='upper right') # Return return @@ -219,9 +271,22 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig # ==================== # # Show one sensitivity # # ==================== # -def show_one_response(rspname, dbname, name, rootdir=None, color="r"): +def show_one_response(rspname, dbname, name, rootdir=None, color='r'): """ Show one response. + + Parameters + ---------- + rspname : str + Response name + dbname : str + Database name + name : str + Name of the response function + rootdir : str, optional + Response root directory + color : str, optional + Color for plot """ # Set-up calibration database caldb = gammalib.GCaldb() @@ -230,7 +295,7 @@ def show_one_response(rspname, dbname, name, rootdir=None, color="r"): if gammalib.dir_exists(dbname): caldb.rootdir(dbname) else: - caldb.open("cta", dbname) + caldb.open('cta', dbname) # Load response function rsp = gammalib.GCTAResponseIrf(rspname, caldb) @@ -252,31 +317,36 @@ def show_one_response(rspname, dbname, name, rootdir=None, color="r"): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Display response information. - """ + + # Print usage information + usage = 'Usage: show_response.py [file]' + if len(sys.argv) < 1: + print(usage) + sys.exit() + + # Check if plotting in file is requested + plotfile = '' + if len(sys.argv) == 2: + plotfile = sys.argv[1] + # Create figures plt.figure(1) - plt.title("Effective area") + plt.title('Effective area') plt.figure(2) - plt.title("Background rate") + plt.title('Background rate') plt.figure(3) - plt.title("Sensitivity") + plt.title('Sensitivity') # Set response dictionary rsps = [ - {'dbname': "dummy", - 'rspname': "cta_dummy_irf", - 'name': "Prod1 (Array E, MPIK)", - 'color': "b"}, - {'dbname': "e", - 'rspname': "IFAE20120510_50h", - 'name': "Prod1 (Array E, IFAE)", - 'color': "g"}, - {'dbname': "aar", - 'rspname': "DESY20140105_50h", - 'name': "Prod2 (Aar, DESY)", - 'color': "r"} + {'dbname': 'irf', + 'rspname': 'cta_dummy_irf.dat', + 'name': 'Prod1 (Array E, MPIK)', + 'color': 'b'}, + {'dbname': 'irf', + 'rspname': 'irf_file.fits', + 'name': 'Prod2 (Aar, IFAE)', + 'color': 'r'} ] # Loop over all responses @@ -295,5 +365,8 @@ if __name__ == '__main__': # Show response show_one_response(rspname, dbname, name, rootdir=rootdir, color=color) - # Show plot - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() diff --git a/examples/show_sensitivity.py b/examples/show_sensitivity.py index 28c6ccf..f72dcd3 100755 --- a/examples/show_sensitivity.py +++ b/examples/show_sensitivity.py @@ -1,8 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays a differential sensitivity plot generated by cssens. +# Display differential sensitivity plot generated by cssens # -# Copyright (C) 2015 Juergen Knoedlseder +# Copyright (C) 2015-2016 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 @@ -23,8 +23,11 @@ import math import gammalib try: import matplotlib.pyplot as plt + plt.figure() + plt.close() except: - sys.exit("This script needs matplotlib") + print('This script needs the "matplotlib" module') + sys.exit() # ===================== # @@ -32,7 +35,17 @@ except: # ===================== # def read_sensitivity(filename): """ - Read sensitivity information from CSV file. + Read sensitivity information from CSV file + + Parameters + ---------- + filename : str + Name of CSV file + + Returns + ------- + sensitivity : dict + Dictionary with sensitivity information """ # Read filename csv = gammalib.GCsv(filename,',') @@ -40,39 +53,39 @@ def read_sensitivity(filename): # Create dictionnary sensitivity = {} for column in range(csv.ncols()): - name = csv[0,column] + name = csv[0,column].rstrip('\r') values = [] for row in range(csv.nrows()-1): values.append(float(csv[row+1,column])) sensitivity[name] = values # Check where file contains differential or integral sensitivity - mode = "Integral" - if sensitivity.has_key("emax"): + mode = 'Integral' + if sensitivity.has_key('emax'): emax_ref = -1.0 - for value in sensitivity["emax"]: + for value in sensitivity['emax']: if emax_ref < 0.0: emax_ref = value elif emax_ref != value: - mode = "Differential" + mode = 'Differential' break # Add mode - sensitivity["mode"] = mode + sensitivity['mode'] = mode # Add linear energy values - if mode == "Differential": - if sensitivity.has_key("loge"): + if mode == 'Differential': + if sensitivity.has_key('loge'): name = 'energy' values = [] - for value in sensitivity["loge"]: + for value in sensitivity['loge']: values.append(math.pow(10.0, value)) sensitivity[name] = values else: - if sensitivity.has_key("emin"): + if sensitivity.has_key('emin'): name = 'energy' values = [] - for value in sensitivity["emin"]: + for value in sensitivity['emin']: values.append(value) sensitivity[name] = values @@ -83,9 +96,18 @@ def read_sensitivity(filename): # ===================== # # Show sensitivity data # # ===================== # -def show_sensitivity(sensitivity, filename): +def show_sensitivity(sensitivity, filename, plotfile=''): """ - Plot sensitivity. + Plot sensitivity data + + Parameters + ---------- + sensitivity : dict + Dictionary with sensitivity information + filename : str + Name of CSV file + plotfile : str, optional + Plot file name """ # Build title title = sensitivity['mode'] + ' sensitivity ('+filename+')' @@ -102,14 +124,17 @@ def show_sensitivity(sensitivity, filename): plt.plot(sensitivity['energy'], sensitivity['sensitivity'], 'ro-') # Set labels - plt.xlabel("Energy (TeV)") - if sensitivity['mode'] == "Differential": - plt.ylabel(r"E$\times$ F(E) (erg cm$^{-2}$ s$^{-1}$)") + plt.xlabel('Energy (TeV)') + if sensitivity['mode'] == 'Differential': + plt.ylabel(r'E $\times$ F(E) (erg cm$^{-2}$ s$^{-1}$)') else: - plt.ylabel(r"E$\times$ F($>$E) (erg cm$^{-2}$ s$^{-1}$)") + plt.ylabel(r'E $\times$ F($>$E) (erg cm$^{-2}$ s$^{-1}$)') - # Show plot - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() # Return return @@ -119,12 +144,11 @@ def show_sensitivity(sensitivity, filename): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Print usage information - usage = "Usage: show_sensitivity [filename]" + usage = 'Usage: show_sensitivity.py [filename] [file]' if len(sys.argv) < 1: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() # Extract parameters @@ -132,9 +156,12 @@ if __name__ == '__main__': filename = 'sensitivity.dat' else: filename = sys.argv[1] + plotfile = '' + if len(sys.argv) == 3: + plotfile = sys.argv[2] # Read sensitivity data sensitivity = read_sensitivity(filename) # Show sensitivity data - show_sensitivity(sensitivity, filename) + show_sensitivity(sensitivity, filename, plotfile=plotfile) diff --git a/examples/show_spectrum.py b/examples/show_spectrum.py index e8945e5..970bdc1 100755 --- a/examples/show_spectrum.py +++ b/examples/show_spectrum.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This script shows a spectrum created with csspec +# Display spectrum generated by csspec # # Copyright (C) 2015-2016 Juergen Knoedlseder # @@ -18,37 +18,41 @@ # along with this program. If not, see . # # ========================================================================== - -import gammalib import sys +import gammalib try: import matplotlib.pyplot as plt + plt.figure() + plt.close() except: - sys.exit("This script needs matplotlib") + print('This script needs the "matplotlib" module') + sys.exit() -# ============ # -# Script entry # -# ============ # -if __name__ == "__main__": - """ +# ============= # +# Show spectrum # +# ============= # +def show_spectrum(filename, plotfile=''): """ + Plot sensitivity data - # Check for given spectrum file - if not len(sys.argv) == 2: - sys.exit("Usage: show_spectrum.py spectrum.fits") - + Parameters + ---------- + filename : str + Name of spectrum FITS file + plotfile : str, optional + Plot file name + """ # Read spectrum file - filename = sys.argv[1] fits = gammalib.GFits(filename) table = fits.table(1) - c_energy = table["Energy"] - c_ed = table["ed_Energy"] - c_eu = table["eu_Energy"] - c_flux = table["Flux"] - c_eflux = table["e_Flux"] - c_ts = table["TS"] - c_upper = table["UpperLimit"] + c_energy = table['Energy'] + c_ed = table['ed_Energy'] + c_eu = table['eu_Energy'] + c_flux = table['Flux'] + c_eflux = table['e_Flux'] + c_ts = table['TS'] + c_upper = table['UpperLimit'] # Initialise arrays to be filled energies = [] @@ -65,25 +69,22 @@ if __name__ == "__main__": nrows = table.nrows() for row in range(nrows): - # Get TS + # Get Test Statistic, flux and flux error ts = c_ts.real(row) flx = c_flux.real(row) e_flx = c_eflux.real(row) - # Switch + # If Test Statistic is larger than 9 and flux error is smaller than + # flux then append flux plots ... if ts > 9.0 and e_flx < flx: - - # Add information energies.append(c_energy.real(row)) flux.append(c_flux.real(row)) ed_engs.append(c_ed.real(row)) eu_engs.append(c_eu.real(row)) e_flux.append(c_eflux.real(row)) - # + # ... otherwise append upper limit else: - - # Add information ul_energies.append(c_energy.real(row)) ul_flux.append(c_upper.real(row)) ul_ed_engs.append(c_ed.real(row)) @@ -93,8 +94,39 @@ if __name__ == "__main__": plt.figure() plt.loglog() plt.grid() - plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], fmt='ro') - plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], yerr=1.0e-11, uplims=True, fmt='ro') - plt.xlabel("Energy (TeV)") - plt.ylabel(r"E dN/dE (erg cm$^{-2}$ s$^{-1}$)") - plt.show() + plt.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs], + fmt='ro') + plt.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs], + yerr=1.0e-11, uplims=True, fmt='ro') + plt.xlabel('Energy (TeV)') + plt.ylabel(r'E $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)') + + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Print usage information + usage = 'Usage: show_spectrum.py filename [file]' + if len(sys.argv) < 2: + print(usage) + sys.exit() + + # Extract parameters + filename = sys.argv[1] + plotfile = '' + if len(sys.argv) == 3: + plotfile = sys.argv[2] + + # Show spectrum + show_spectrum(filename, plotfile=plotfile) diff --git a/examples/show_ts_distribution.py b/examples/show_ts_distribution.py index 027b4c8..13f580f 100755 --- a/examples/show_ts_distribution.py +++ b/examples/show_ts_distribution.py @@ -1,12 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script displays the TS distribution values generated using cstsdist. +# Display TS distribution generated by cstsdist # -# Required 3rd party modules: -# - matplotlib -# - numpy -# -# Copyright (C) 2011 Jurgen Knodlseder +# Copyright (C) 2011-2016 Jurgen Knodlseder # # 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 @@ -22,20 +18,42 @@ # along with this program. If not, see . # # ========================================================================== -import matplotlib.mlab as mlab -import matplotlib.pyplot as plt -import numpy as np import sys import csv import math +try: + import matplotlib.pyplot as plt + import matplotlib.mlab as mlab + plt.figure() + plt.close() +except: + print('This script needs the "matplotlib" module') + sys.exit() +try: + import numpy as np +except: + print('This script needs the "numpy" module') + sys.exit() # ==================== # # Read TS distribution # # ==================== # -def read_ts(filename, tsname="TS"): +def read_ts(filename, tsname='TS'): """ - Read TS distribution values from CSV file. + Read TS distribution values from CSV file + + Parameters + ---------- + filename : str + Name of CSV file + tsname : str, optional + Column name of TS + + Returns + ------- + a : `numpy.array` + Array of TS values """ # Initialise list values = [] @@ -53,9 +71,9 @@ def read_ts(filename, tsname="TS"): try: index = row.index(tsname) except: - sys.stdout.write('ERROR: Column "TS" not found in file\n') - sys.stdout.write(row+"\n") - raise NameError("TS") + print('ERROR: Column "'+tsname+'" not found in file') + print(row) + raise NameError(tsname) # Handle data rows else: @@ -76,7 +94,17 @@ def read_ts(filename, tsname="TS"): # =========== # def erf(x): """ - Compute error function. + Compute error function + + Parameters + ---------- + x : float + Value + + Returns + ------- + erf(-x) : float + Error function """ # save the sign of x sign = 1 @@ -103,9 +131,20 @@ def erf(x): # =============================== # # Display cumulative distribution # # =============================== # -def dist_cdf(values, nbins, title): +def dist_cdf(values, nbins, title, expect='', plotfile=''): """ - Displays cumulative TS distribution and overplots expectation. + Display cumulative TS distribution and overplots expectation + + Parameters + ---------- + values : float + TS values + nbins : int + Number of bins + title : str + Plot title + plotfile : str, optional + Plot file name """ # Set range and adapt number of bins. We make sure that we have a # histogram bin centred on 0 that will capture the TS=0 part. @@ -142,19 +181,23 @@ def dist_cdf(values, nbins, title): hist_expect.append(norm) # Show expected distribution (half) - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label="simulations") - l = plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label="expectation") - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2) + plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label='simulations') + if expect == 'show': + plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label='expectation') + plt.semilogy(hist_x, hist_y, 'r-', linewidth=2) # Set plot plt.xlabel('TS') plt.ylabel('Number of trials > TS') plt.title(title) plt.grid(True) - plt.legend(loc="upper right") + plt.legend(loc='upper right') - # Show histogram - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() # Return return @@ -163,10 +206,20 @@ def dist_cdf(values, nbins, title): # ======================================== # # Display probability density distribution # # ======================================== # -def dist_pdf(values, nbins, title): +def dist_pdf(values, nbins, title, expect='', plotfile=''): """ - Display probability density TS distribution and overplots - expectation. + Display probability density TS distribution and overplots expectation + + Parameters + ---------- + values : float + TS values + nbins : int + Number of bins + title : str + Plot title + plotfile : str, optional + Plot file name """ # Set range and adapt number of bins. We make sure that we have a # histogram bin centred on 0 that will capture the TS=0 part. @@ -180,9 +233,9 @@ def dist_pdf(values, nbins, title): # Create histogram n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value], align='mid', facecolor='green', - log=True, label="simulations") + log=True, label='simulations') - # Create expected distribution (onyl for positive TS). We compute the + # Create expected distribution (only for positive TS). We compute the # full and half of the distribution, as for positively constrained # quantities half of the values have TS=0 x = [] @@ -202,47 +255,48 @@ def dist_pdf(values, nbins, title): y_half.append(0.5*y) # Show expected distribution (half) - l = plt.semilogy(x, y_half, 'ro', linewidth=2, label="expectation") + if expect == 'show': + plt.semilogy(x, y_half, 'ro', linewidth=2, label='expectation') # Set plot plt.xlabel('TS') plt.ylabel('Number of trials') plt.title(title) plt.grid(True) - plt.legend(loc="upper right") + plt.legend(loc='upper right') - # Show histogram - plt.show() + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() + + # Return + return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Display TS distribution generated using cstsdist using matplotlib. - """ + # Print usage information - usage = "Usage: show_ts_distribution filename" \ - " [-n bins] [-c column] [-t title] [-p plot]" + usage = 'Usage: show_ts_distribution.py filename' \ + ' [-n bins] [-c column] [-t title] [-p plot] [-f file]' if len(sys.argv) < 2: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() # Extract parameters filename = sys.argv[1] - # Set default parameters - nbins = 30 - tsname = "TS" - title = "TS distribution" - plot = "pdf" - # Parameter dictionnary - pars = [{'option': '-n', 'value': nbins}, - {'option': '-c', 'value': tsname}, - {'option': '-t', 'value': title}, - {'option': '-p', 'value': plot}] + pars = [{'option': '-n', 'value': 30}, + {'option': '-c', 'value': 'TS'}, + {'option': '-t', 'value': 'TS distribution'}, + {'option': '-p', 'value': 'pdf'}, + {'option': '-e', 'value': ''}, + {'option': '-f', 'value': ''}] # Gather parameters from command line i = 2 @@ -256,27 +310,29 @@ if __name__ == '__main__': try: par['value'] = sys.argv[i] except: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() else: - sys.stdout.write(usage+"\n") + print(usage) sys.exit() # Next item i += 1 # Recover parameters - nbins = int(pars[0]['value']) - tsname = pars[1]['value'] - title = pars[2]['value'] - plot = pars[3]['value'] + nbins = int(pars[0]['value']) + tsname = pars[1]['value'] + title = pars[2]['value'] + plot = pars[3]['value'] + expect = pars[4]['value'] + plotfile = pars[5]['value'] # Read values from CSV file values = read_ts(filename, tsname=tsname) - sys.stdout.write(len(values), "values read.\n") + print('%d values read.' % len(values)) # Show histogram - if plot == "pdf": - dist_pdf(values, nbins, title) + if plot == 'pdf': + dist_pdf(values, nbins, title, expect=expect, plotfile=plotfile) else: - dist_cdf(values, nbins, title) + dist_cdf(values, nbins, title, expect=expect, plotfile=plotfile) diff --git a/sonar-project.properties b/sonar-project.properties new file mode 100644 index 0000000..5c75e43 --- /dev/null +++ b/sonar-project.properties @@ -0,0 +1,42 @@ +# Required metadata +sonar.projectKey=ctools +sonar.projectName=ctools +sonar.projectVersion=1.1.0.dev1 + +# Description +sonar.projectDescription="Cherenkov Telescope Array Science Analysis Software" + +# Links +sonar.links.ci=https://cta-jenkins.irap.omp.eu +sonar.links.homepage=http://cta.irap.omp.eu/ctools +sonar.links.issue=https://cta-redmine.irap.omp.eu/projects/ctools/issues +sonar.links.scm=https://cta-gitlab.irap.omp.eu/ctools/ctools + +# Path to source directories +sonar.sources=. + +# Path to tests +#sonar.tests=test + +# Path to include directories +sonar.cxx.includeDirectories=pyext,/home/jenkins/gcc/gcc492/include/c++/4.9.2,/home/jenkins/gcc/gcc492/include/c++/4.9.2/tr1,/home/jenkins/gcc/gcc492/include/c++/4.9.2/ext,/home/jenkins/gcc/gcc492/include/c++/4.9.2/x86_64-unknown-linux-gnu,/home/jenkins/gcc/gcc492/lib/gcc/x86_64-unknown-linux-gnu/4.9.2/include,/usr/include,/usr/include/bits,/usr/include/gnu,/usr/include/linux +sonar.cxx.missingIncludeWarnings=False + +# Exclusions +sonar.exclusions=pyext/ctools/**,pyext/build/**,pyext/setup.py,doc/**,test/dev/** +sonar.coverage.exclusions=pyext/ctools/**,pyext/build/**,examples/**,doc/**,test/dev/** +#sonar.test.exclusions= + +# C++ configuration +sonar.cxx.cppcheck.reportPath=cppcheck-result.xml +sonar.cxx.cppncss.reportPath=cppncss-raw-report.xml +#sonar.cxx.valgrind.reportPath=test/valgrind_test_*.xml +sonar.cxx.xunit.reportPath=test/reports/*.xml +#sonar.cxx.xunit.provideDetails=true +#sonar.cxx.xunit.xsltURL=cppunit-1.x-to-junit-1.0.xsl +sonar.cxx.coverage.reportPath=coverage-report-cxx.xml + +# Python configuration +sonar.python.coverage.forceZeroCoverage=True +sonar.python.coverage.reportPath=coverage-report-python.xml +#sonar.python.xunit.reportPath=test/reports/*.xml diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index ec19f93..2c78676 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -52,6 +52,7 @@ #define G_SETUP_OBSERVATION "ctool::setup_observations(GObservations&)" #define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)" #define G_CREATE_EBOUNDS "ctool::create_ebounds()" +#define G_SET_OBS_RESPONSE "ctool::set_obs_response(GCTAObservation*)" #define G_PROVIDE_HELP "ctool::provide_help()" /* __ Debug definitions __________________________________________________ */ @@ -931,6 +932,9 @@ void ctool::set_response(GObservations& obs) * * @param[in,out] obs CTA observation * + * @exception GException::invalid_value + * Energy dispersion requested but not energy dispersion cube given + * * Set the response for one CTA observation. If the CTA observation contains * a counts cube the method first attempts to set the response information * using the following user parameters @@ -964,13 +968,23 @@ void ctool::set_obs_response(GCTAObservation* obs) // response information for a stacked analysis is provided, and if so, // set that response information for the CTA observation if (dynamic_cast(obs->events()) != NULL) { - if (has_par("expcube") && has_par("psfcube") && - has_par("edispcube") && has_par("bkgcube")) { + + // Check if response cube parameters are provided. We need all four + // response cube parameters. + if (has_par("expcube") && has_par("psfcube") && + has_par("bkgcube") && has_par("edispcube")) { + + // Check if we need to query the energy dispersion cube. We do + // not need to query the cube is an "edisp" parameter exists and + // if this parameter is set to "no". + bool query_edisp = true; + if (has_par("edisp")) { + query_edisp = (*this)["edisp"].boolean(); + } // Get filenames GFilename expcube((*this)["expcube"].filename()); GFilename psfcube((*this)["psfcube"].filename()); - GFilename edispcube((*this)["edispcube"].filename()); GFilename bkgcube((*this)["bkgcube"].filename()); // Extract stacked response information if available @@ -981,21 +995,53 @@ void ctool::set_obs_response(GCTAObservation* obs) !psfcube.is_empty() && !bkgcube.is_empty()) { - // Get exposure, PSF and background cubes + // Load exposure, PSF and background cubes GCTACubeExposure exposure(expcube); GCTACubePsf psf(psfcube); GCTACubeBackground background(bkgcube); - // If an energy dispersion cube is available, get also - // this cube, and set all information as stacked response - // for this CTA observation - if ((gammalib::toupper(edispcube.url()) != "NONE") && - !edispcube.is_empty()) { - GCTACubeEdisp edisp(edispcube); - obs->response(exposure, psf, edisp, background); - } + // If querying of energy dispersion cube is requested then + // query it now + if (query_edisp) { + + // Get filename + GFilename edispcube((*this)["edispcube"].filename()); + + // If filename is valid then use energy dispersion ... + if ((gammalib::toupper(edispcube.url()) != "NONE") && + !edispcube.is_empty()) { + + // Load energy dispersion cube + GCTACubeEdisp edisp(edispcube); + + // Set response with all four cubes + obs->response(exposure, psf, edisp, background); + + } // endif: energy dispersion cube was provided + + // ... otherwise work without energy dispersion + else { + + // If energy dispersion was requested but no cube + // was provided then throw an exception + if (has_par("edisp") && (*this)["edisp"].boolean()) { + std::string msg = "Energy dispersion requested but " + "no energy dispersion cube was " + "specified."; + throw GException::invalid_value(G_SET_OBS_RESPONSE, + msg); + } + + // Set response without energy dispersion + obs->response(exposure, psf, background); + + } // endelse: no energy dispersion was available + + } // endif: energy dispersion needed + + // ... otherwise work without energy dispersion else { - obs->response(exposure, psf, background); + obs->response(exposure, psf, background); } // Signal that response is available diff --git a/test/Makefile.am b/test/Makefile.am index 3d9f41e..2213dab 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -9,6 +9,7 @@ if WITH_PYTHON TEST_CSCRIPTS = test_cscripts.sh TEST_PYTHON_CTOOLS = test_python_ctools.sh TEST_PYTHON_CSCRIPTS = test_python_cscripts.sh + TEST_PYTHON_EXAMPLES = test_examples.sh TEST_PYTHON_ENV = PYTHONPATH=$(top_builddir):$(top_builddir)/pyext:$(top_builddir)/pyext/build/ctools/ctools:$(pythondir):$(pyexecdir):$(PYTHONPATH) endif @@ -18,7 +19,8 @@ TESTS_ENVIRONMENT = @RUNSHARED@=$(top_builddir)/src/.libs$(TEST_ENV_DIR):$(LIBDI # Tests (these are the tests to be executed) TESTS = $(TEST_PYTHON_CTOOLS) \ - $(TEST_PYTHON_CSCRIPTS) + $(TEST_PYTHON_CSCRIPTS) \ + $(TEST_PYTHON_EXAMPLES) # Add Valgrind rule valgrind: diff --git a/test/README b/test/README deleted file mode 100644 index c7493dd..0000000 --- a/test/README +++ /dev/null @@ -1,21 +0,0 @@ -This folder contains files used for unit testing and files that serve -as example scripts for ctools. - -test_cscripts.sh - Unit test for cscripts. - -test_ctools.sh - Unit test for ctools. - -test_python.py - Unit test for ctools Python interface. - -example_models.py - Example Python pipeline for unbinned analysis that is executed - on a number of models (point source, disk source, Gaussian - source, Shell source). - -example_survey.py - This script illustrates how a survey can be simulated and analysed. - The script can be easily extended to adapt specific survey - strategies and simulated source models. diff --git a/test/README.md b/test/README.md new file mode 100644 index 0000000..c93c30e --- /dev/null +++ b/test/README.md @@ -0,0 +1,4 @@ +ctools unit test +================ + +This folder contains files used for ctools unit testing. diff --git a/test/data/butterfly.txt b/test/data/butterfly.txt new file mode 100644 index 0000000..b6f559a --- /dev/null +++ b/test/data/butterfly.txt @@ -0,0 +1,100 @@ +103514 7.97876e-15 7.61983e-15 8.34788e-15 +110917 6.72257e-15 6.4272e-15 7.02571e-15 +118850 5.66416e-15 5.42109e-15 5.91312e-15 +127350 4.77239e-15 4.57237e-15 4.97683e-15 +136458 4.02101e-15 3.8564e-15 4.18893e-15 +146218 3.38794e-15 3.25246e-15 3.52586e-15 +156675 2.85454e-15 2.743e-15 2.96785e-15 +167880 2.40511e-15 2.31326e-15 2.49824e-15 +179887 2.02645e-15 1.95077e-15 2.10302e-15 +192752 1.7074e-15 1.64502e-15 1.7704e-15 +206538 1.43859e-15 1.38713e-15 1.49045e-15 +221309 1.21209e-15 1.16962e-15 1.25483e-15 +237137 1.02126e-15 9.86165e-16 1.05651e-15 +254097 8.6047e-16 8.31443e-16 8.89582e-16 +272270 7.24996e-16 7.00958e-16 7.49072e-16 +291743 6.10851e-16 5.90916e-16 6.30794e-16 +312608 5.14678e-16 4.9812e-16 5.31227e-16 +334965 4.33646e-16 4.1987e-16 4.47406e-16 +358922 3.65372e-16 3.53889e-16 3.76836e-16 +384592 3.07847e-16 2.98257e-16 3.17421e-16 +412098 2.59379e-16 2.51353e-16 2.67393e-16 +441570 2.18542e-16 2.1181e-16 2.25268e-16 +473151 1.84135e-16 1.78475e-16 1.89793e-16 +506991 1.55144e-16 1.50375e-16 1.59917e-16 +543250 1.30718e-16 1.26691e-16 1.34754e-16 +582103 1.10138e-16 1.06729e-16 1.13559e-16 +623735 9.27973e-17 8.99058e-17 9.57049e-17 +668344 7.81872e-17 7.57292e-17 8.06636e-17 +716143 6.58773e-17 6.37836e-17 6.79909e-17 +767361 5.55054e-17 5.37188e-17 5.7313e-17 +822243 4.67666e-17 4.52393e-17 4.83151e-17 +881049 3.94036e-17 3.80959e-17 4.07323e-17 +944061 3.31998e-17 3.20787e-17 3.43416e-17 +1.01158e+06 2.79728e-17 2.70103e-17 2.89551e-17 +1.08393e+06 2.35687e-17 2.27416e-17 2.44147e-17 +1.16145e+06 1.9858e-17 1.91465e-17 2.05873e-17 +1.24451e+06 1.67315e-17 1.6119e-17 1.73606e-17 +1.33352e+06 1.40973e-17 1.35697e-17 1.46403e-17 +1.42889e+06 1.18778e-17 1.1423e-17 1.23467e-17 +1.53109e+06 1.00077e-17 9.61563e-18 1.04128e-17 +1.64059e+06 8.4321e-18 8.09391e-18 8.78208e-18 +1.75792e+06 7.10453e-18 6.81279e-18 7.40699e-18 +1.88365e+06 5.98598e-18 5.73426e-18 6.24739e-18 +2.01837e+06 5.04354e-18 4.82635e-18 5.26946e-18 +2.16272e+06 4.24948e-18 4.06207e-18 4.44474e-18 +2.31739e+06 3.58043e-18 3.41874e-18 3.74917e-18 +2.48313e+06 3.01672e-18 2.87723e-18 3.16252e-18 +2.66073e+06 2.54177e-18 2.42144e-18 2.66772e-18 +2.85102e+06 2.14159e-18 2.03782e-18 2.25037e-18 +3.05492e+06 1.80441e-18 1.71493e-18 1.89835e-18 +3.27341e+06 1.52032e-18 1.44319e-18 1.60142e-18 +3.50752e+06 1.28096e-18 1.21448e-18 1.35096e-18 +3.75837e+06 1.07928e-18 1.02201e-18 1.13968e-18 +4.02717e+06 9.09359e-19 8.60023e-19 9.61461e-19 +4.31519e+06 7.66188e-19 7.23704e-19 8.11119e-19 +4.62381e+06 6.45558e-19 6.08984e-19 6.84293e-19 +4.9545e+06 5.43921e-19 5.12443e-19 5.77303e-19 +5.30884e+06 4.58285e-19 4.31203e-19 4.87047e-19 +5.68853e+06 3.86132e-19 3.62839e-19 4.10905e-19 +6.09537e+06 3.25339e-19 3.0531e-19 3.46669e-19 +6.53131e+06 2.74117e-19 2.569e-19 2.92478e-19 +6.99842e+06 2.30959e-19 2.16165e-19 2.4676e-19 +7.49894e+06 1.94597e-19 1.81887e-19 2.0819e-19 +8.03526e+06 1.63959e-19 1.53044e-19 1.7565e-19 +8.60994e+06 1.38145e-19 1.28773e-19 1.48196e-19 +9.22571e+06 1.16395e-19 1.08351e-19 1.25035e-19 +9.88553e+06 9.80699e-20 9.11675e-20 1.05493e-19 +1.05925e+07 8.26296e-20 7.67083e-20 8.9007e-20 +1.13501e+07 6.96202e-20 6.45422e-20 7.50971e-20 +1.21619e+07 5.86591e-20 5.43052e-20 6.33617e-20 +1.30317e+07 4.94237e-20 4.56918e-20 5.34601e-20 +1.39637e+07 4.16424e-20 3.84443e-20 4.51063e-20 +1.49624e+07 3.50861e-20 3.23463e-20 3.80579e-20 +1.60325e+07 2.95621e-20 2.72155e-20 3.2111e-20 +1.71791e+07 2.49078e-20 2.28983e-20 2.70936e-20 +1.84077e+07 2.09863e-20 1.9266e-20 2.28601e-20 +1.97242e+07 1.76822e-20 1.62098e-20 1.92882e-20 +2.11349e+07 1.48983e-20 1.36383e-20 1.62745e-20 +2.26464e+07 1.25527e-20 1.14748e-20 1.37317e-20 +2.42661e+07 1.05763e-20 9.65444e-21 1.15863e-20 +2.60016e+07 8.91119e-21 8.12282e-21 9.77606e-21 +2.78612e+07 7.5082e-21 6.83419e-21 8.24867e-21 +2.98538e+07 6.32609e-21 5.74997e-21 6.95995e-21 +3.1989e+07 5.3301e-21 4.83773e-21 5.87259e-21 +3.42768e+07 4.49092e-21 4.07022e-21 4.95511e-21 +3.67282e+07 3.78386e-21 3.42447e-21 4.18097e-21 +3.9355e+07 3.18813e-21 2.88116e-21 3.5278e-21 +4.21697e+07 2.68618e-21 2.42404e-21 2.97667e-21 +4.51856e+07 2.26327e-21 2.03945e-21 2.51165e-21 +4.84172e+07 1.90693e-21 1.71587e-21 2.11927e-21 +5.188e+07 1.6067e-21 1.44362e-21 1.7882e-21 +5.55904e+07 1.35374e-21 1.21457e-21 1.50885e-21 +5.95662e+07 1.14061e-21 1.02186e-21 1.27314e-21 +6.38263e+07 9.61027e-22 8.59731e-22 1.07426e-21 +6.83912e+07 8.09721e-22 7.23319e-22 9.06442e-22 +7.32825e+07 6.82237e-22 6.0855e-22 7.64845e-22 +7.85236e+07 5.74825e-22 5.11991e-22 6.45368e-22 +8.41395e+07 4.84324e-22 4.30753e-22 5.44554e-22 +9.01571e+07 4.08071e-22 3.62405e-22 4.59489e-22 +9.66051e+07 3.43824e-22 3.04901e-22 3.87712e-22 diff --git a/test/data/obs_stacked.xml b/test/data/obs_stacked.xml index 151c335..c54bcea 100644 --- a/test/data/obs_stacked.xml +++ b/test/data/obs_stacked.xml @@ -1,9 +1,9 @@ - - - - + + + + diff --git a/test/data/obs_unbinned.xml b/test/data/obs_unbinned.xml index 4adfdbb..fe07707 100644 --- a/test/data/obs_unbinned.xml +++ b/test/data/obs_unbinned.xml @@ -1,7 +1,7 @@ - + diff --git a/test/data/pha.fits b/test/data/pha.fits new file mode 100644 index 0000000..d3e5718 --- /dev/null +++ b/test/data/pha.fits diff --git a/test/data/pull.dat b/test/data/pull.dat new file mode 100644 index 0000000..3a1ecc2 --- /dev/null +++ b/test/data/pull.dat @@ -0,0 +1,101 @@ +LogL,Sim_Events,Npred_Events,Crab_Prefactor,Pull_Crab_Prefactor,Unc_Crab_Prefactor,Crab_Index,Pull_Crab_Index,Unc_Crab_Index,Background_Prefactor,Pull_Background_Prefactor,Unc_Background_Prefactor,Background_Index,Pull_Background_Index,Unc_Background_Index +104702.14348660898,15105.0,15104.999075727912,5.67748323458501e-16,-0.08501991498574554,2.6484107186845484e-17,-2.512994429094607,-1.1338279073158557,0.02910003262551158,0.9958218078266116,-0.17872001258319373,0.023378423675096273,-0.012229059990362811,-0.9236754329409402,0.013239563979119857 +103330.80595191609,14818.0,14817.997797150003,5.732404092766631e-16,0.12283753926158997,2.637963358873928e-17,-2.484339613991382,-0.1554901325428316,0.027909256493730588,0.9787908971847574,-0.9101986145091182,0.023301620632196762,-0.008200032682593691,-0.6103278351854464,0.013435455848252658 +104268.96595026793,14979.0,14978.999805565223,5.572825610413261e-16,-0.48658024310290293,2.6136365253088216e-17,-2.4950417511405485,-0.516821157597015,0.02910436409083146,1.0063046580826853,0.2716574371824787,0.02320811882816348,-0.0008417569678340521,-0.06468934227489305,0.013012297516599597 +104821.06173253663,15107.0,15106.996064070843,5.572524365485819e-16,-0.4875817324560644,2.614446482070941e-17,-2.488742729886503,-0.3182574459681672,0.027470621653193025,0.9996316452115123,-0.015557106670414662,0.023677589688849032,-0.010871023195874276,-0.8107102927475592,0.013409257651129044 +103881.48989070744,14862.0,14861.992397431473,5.301740735337788e-16,-1.5134244536513575,2.631510702112374e-17,-2.465301587043766,0.4827690936075433,0.030446052058549868,1.0115707161034393,0.4942920042462879,0.023408665331503164,0.0048734492051994045,0.3732789541610912,0.013055783485441917 +104709.9140731044,15024.0,15023.989651557333,5.226837712898727e-16,-1.7954633953701884,2.6353212675979732e-17,-2.4343195831835844,1.5525972082532988,0.029421936722278975,1.017947066365573,0.7598967774868103,0.02361776875134134,0.0009326118036593681,0.07132952806775084,0.013074694715119194 +103713.58340281199,14892.0,14891.999561522098,5.591199213471972e-16,-0.41862333613425523,2.599013889974217e-17,-2.470115751370293,0.35799682935892585,0.027609877571840587,0.9794463837010542,-0.8909844728935377,0.023068433765401575,-0.012917295469889406,-0.9740485577904611,0.013261449202482362 +104795.70589373197,15029.0,15028.996906992426,5.435055239409328e-16,-1.0152303535510279,2.609700937959154e-17,-2.4655187514911407,0.49983528642007014,0.028972041194964666,1.0209631987169243,0.8832401674398217,0.023734426365241833,0.00466053101676357,0.3544111893768436,0.013150067369368664 +104193.29654018967,14931.0,14930.996770551044,5.485461677632551e-16,-0.8193085256981937,2.6185291088558448e-17,-2.479088507285966,0.031657689546980165,0.028792142669831974,1.0185759089559494,0.7867635712209643,0.023610535153682555,0.008184752146355114,0.6243663808647385,0.013108893106991681 +103628.40847429863,14848.0,14847.998122677556,6.014015930608307e-16,1.1849643924993984,2.6500030937298085e-17,-2.473150703559013,0.25987554591515094,0.026356063695286987,0.9902956776180022,-0.4187159470974005,0.02317638592289015,0.001885900145097897,0.14262284727070834,0.013222987629172232 +105660.09206884849,15205.0,15204.992131202547,6.109854020247515e-16,1.5187149825791957,2.6986895167878697e-17,-2.4922552592625755,-0.44934245174528803,0.0272737623943052,1.00323842606026,0.13976115495783228,0.023171145524928653,-0.005994371257142875,-0.46075233874827176,0.013009963820103039 +103835.59625385483,14858.0,14857.999995507591,5.297499346626201e-16,-1.54328568208427,2.6080761199715524e-17,-2.4319384773924724,1.701399801117752,0.028248223948276648,1.0130018036654878,0.5601780025624995,0.023210128934038485,0.006713522884626453,0.5203716357596322,0.012901400505479386 +103482.23089354827,14807.0,14806.999524388077,5.736758548183581e-16,0.13761738296308404,2.6710686827580283e-17,-2.4701072439260923,0.3421196454712564,0.028916071336040444,0.9891275161084382,-0.4697332187670207,0.02314608262132377,-0.0006918586515145906,-0.05255682430971379,0.013164011726384277 +102857.43398505163,14731.0,14730.999998518842,5.258007078860843e-16,-1.7041864102309696,2.5935714455043313e-17,-2.4853332222867714,-0.17651127204646544,0.03021462722996824,0.9985936056978045,-0.060149237976978034,0.023381747624695796,0.0019325678652796843,0.14603466473388227,0.013233624145345124 +103817.91519576444,14910.0,14909.994926336014,5.329622298577745e-16,-1.4389390287441313,2.5739638304585517e-17,-2.4916943762429473,-0.39969910169305867,0.02925794977624896,1.0079906020378802,0.3407215111753502,0.023452003398070874,0.00020664358393950301,0.015738123620475775,0.013130128401752638 +103967.78649231851,14870.0,14869.990058174044,5.260603968355675e-16,-1.6789644755360136,2.6170656857051498e-17,-2.4301686389439427,1.7352508788074,0.02871709311007826,1.0126332673580076,0.536452771161907,0.02354963575012417,0.005447900085284499,0.4147492279428597,0.013135407417889301 +105284.40661958465,15058.0,15057.999998318428,4.927243299991711e-16,-2.9598366118818897,2.6108086402680312e-17,-2.4504493609494307,0.9673695847716703,0.03054741384860079,1.0639477867734564,2.6531725316838397,0.0241023853555697,0.02273846071513058,1.7729798337272424,0.012824996811908856 +104830.07226830821,15032.0,15031.998093973076,5.743212958004981e-16,0.16282839716227676,2.653895681470984e-17,-2.460255056292079,0.7331750878331816,0.02693073460293763,1.009527339928993,0.40649631413313464,0.023437703117457166,0.0017910275983860455,0.13653669474155755,0.013117554967740933 +104698.0498891591,15040.0,15039.998968232076,5.570480545515681e-16,-0.4911050777527014,2.6373063597102404e-17,-2.4920533487352197,-0.4183571629577559,0.02881114464493293,1.0144387152840497,0.6180934305102588,0.0233600853387651,0.0013844182971720093,0.10655399897505911,0.012992645142263101 +104925.21327728468,15098.0,15097.992620766281,5.87969306777613e-16,0.6774288392926258,2.6525748145556783e-17,-2.5084885893716833,-1.0402720874605131,0.027385709676426055,1.0166071274045811,0.716918888388743,0.023164583432730087,0.003874758455068439,0.30154828940328254,0.012849545466618257 +103716.61350507324,14871.0,14870.999265696919,5.610566829273711e-16,-0.33697545311885546,2.653996601192925e-17,-2.454976233044967,0.886385961633361,0.028231231132001506,0.9787988077612177,-0.922234954456611,0.022988927210283613,-0.011704744624020038,-0.8844478536183279,0.013233956729201437 +103344.09380500649,14825.0,14824.999999989699,5.185881816208363e-16,-1.9630789475358086,2.6189378905876087e-17,-2.418893532303184,2.1636019917057876,0.028242933742466933,0.9802133791662438,-0.8562764603125603,0.02310774819914305,-0.013591256173913914,-1.026686644454609,0.01323797893673175 +105469.3219667073,15214.0,15213.997120280635,5.796556830349523e-16,0.36143547823531713,2.6714818041924188e-17,-2.481092011680722,-0.03892934595915295,0.028051118091422435,0.9960380118113645,-0.16810935916980949,0.023567921549408868,-0.014903616203581289,-1.1156393733398904,0.013358811601426653 +104469.9971681874,14961.0,14960.999082266428,5.947706132978403e-16,0.9305884258012391,2.6618226286784883e-17,-2.475757601829857,0.1572808660063672,0.02697339020228572,0.9946878034334904,-0.22520963190307067,0.023587785840331894,-0.0017152659481122583,-0.12790022377679958,0.01341096909342076 +104826.19001967587,15063.0,15062.99608195958,5.287050347495291e-16,-1.5838597523756837,2.6072362271047778e-17,-2.4579789012405966,0.7746752179149065,0.028426233665606017,1.0185925724596263,0.7958686132990104,0.02336135908483304,-0.0002351029057570468,-0.018165082808064857,0.012942572750214319 +103320.95695538132,14818.0,14817.998687527173,5.478498323029024e-16,-0.8523687875773152,2.598660112843289e-17,-2.490995186934019,-0.37945854402407125,0.028975989886582468,0.9953315815622235,-0.19898386538883064,0.02346129133964744,-0.0013328500645812436,-0.09997748749545347,0.01333150190078397 +104825.19976994366,15116.0,15115.990754728477,5.315233940698427e-16,-1.460682378789519,2.6341528102805798e-17,-2.464718006751812,0.5076150170072764,0.030105479026773865,0.9928387367913682,-0.3049904020195074,0.023480290400003188,-0.018504566807966617,-1.3859269563155745,0.01335176195516118 +103770.47872478547,14858.0,14857.991542917502,5.241689052795333e-16,-1.7674486406832426,2.5930651485720014e-17,-2.470999774704249,0.31047198034973555,0.028988848802434315,1.0142624255798223,0.613387112964016,0.02325191592451823,0.005780894385199397,0.4468102900554024,0.012938140669237033 +103629.36892843783,14923.0,14922.999996487131,5.984682923028192e-16,1.0727263172672183,2.653826222455558e-17,-2.517783256274881,-1.3204723419015199,0.028613440112249505,0.9670801291995181,-1.4156504768341795,0.023254236366380916,-0.017944371819906522,-1.3217953180960127,0.013575756831817646 +104202.31840266004,14969.0,14968.993306219852,6.11363636982842e-16,1.5511916803318888,2.6665716111881122e-17,-2.4805177040241033,-0.01902684336053044,0.027209138914617383,0.9743574698978508,-1.1040125515345534,0.02322666537305819,-0.012864906504281222,-0.9561276967352017,0.013455217904689712 +104558.3210911674,14988.0,14987.99592260394,5.727035536558679e-16,0.10076888899713153,2.6829249411937305e-17,-2.4539572754877907,0.9040406708094028,0.028807027552082087,0.991942847147923,-0.3490713839169001,0.02308167676670701,-0.007233321441344711,-0.5529892440587004,0.013080401687843472 +103295.47162555656,14883.0,14883.000000254477,5.800709280699981e-16,0.3789335257832652,2.6577031021948253e-17,-2.498719352789656,-0.6638277678957628,0.028199110816038476,0.9595872183999479,-1.7482932734147603,0.023115562025310562,-0.02304437391034645,-1.6988602430501567,0.013564608392372682 +103258.47720670803,14845.0,14844.992307538843,5.83283553420523e-16,0.4952239600807473,2.6823325386673746e-17,-2.4589665237005023,0.7632059346150839,0.02755937204564027,0.964755489938738,-1.5351536401647652,0.022958294947911025,-0.016523138292911025,-1.2320888966887555,0.01341067055901326 +103657.9671271524,14817.0,14816.999863026358,5.363242117977063e-16,-1.2786443217150414,2.6337103782797478e-17,-2.4323288184790695,1.6135859025493344,0.029543627919414728,1.0151591013135277,0.6351636323121661,0.023866450379636053,0.010915006913822117,0.8190878811605987,0.013325806869900458 +104392.7082185512,15006.0,15005.995853606486,6.005823640352908e-16,1.1354123615360034,2.6935028251690447e-17,-2.4730683275167515,0.25611675736022566,0.027064501966574228,0.9714394796871246,-1.2460905368592647,0.022920100480709322,-0.017794798460961776,-1.3401186899781854,0.013278524203890809 +105182.36065597368,15075.0,15074.989261140847,5.26313996325609e-16,-1.67692009148971,2.6051332974120466e-17,-2.4565236365198113,0.8036379946781819,0.029212610199683024,1.029479540250354,1.2428947693473755,0.023718452259504866,0.005601907621198788,0.43008379976278804,0.013025153759077906 +103071.25151722605,14836.0,14835.999997691571,5.626870367730519e-16,-0.2758504395736606,2.651060929339331e-17,-2.4857663826203016,-0.19466445034259733,0.02962216578401001,0.9595379837552862,-1.7643536985481396,0.022933052640187388,-0.022964614426049572,-1.7082617255518946,0.013443264625407625 +103697.65696037607,14939.0,14938.998760928765,5.389129205397533e-16,-1.1876025861552415,2.6176331899788434e-17,-2.4690490243020475,0.3886240085819805,0.02817884499187428,0.9885163863773131,-0.49579679121988196,0.0231619361521725,-0.012105401262481892,-0.9167134357594129,0.01320521854515384 +104637.09745185435,15053.0,15052.996846537539,5.406589368501109e-16,-1.1101425795758761,2.6429995290423552e-17,-2.4389502165668464,1.4305536253574138,0.028695032961730167,0.9973647045748609,-0.11185075823051355,0.02356081860176562,-0.010884674103109395,-0.8187292691775412,0.013294595066869482 +104582.78060549527,15024.0,15023.999395677984,5.827997309903262e-16,0.48195052823758067,2.6558184378659893e-17,-2.480763172768796,-0.027230974562203435,0.028025907301000918,0.9972408492496609,-0.11684592426392378,0.023613581455411925,-0.004849155581888444,-0.36150282462815947,0.013413880201008851 +104660.06320954843,15025.0,15024.998351527545,5.950507262965721e-16,0.9335275835620825,2.6834478956674745e-17,-2.468672941074628,0.4151788322094137,0.027282361350393892,0.9942447528437467,-0.24704233712761414,0.023296602611399077,-0.004813342700068887,-0.3639786391517183,0.013224245003186924 +105487.22379340562,15160.0,15159.998753761023,5.897989793185405e-16,0.7404510509478629,2.6739079231767564e-17,-2.4830216972591455,-0.10999404551766875,0.027471462158922832,1.008007621452472,0.3392549160552585,0.023603553179367128,-0.003582415027088905,-0.2709663149118288,0.013220886988312945 +104309.2333224919,14961.0,14960.997133833138,5.993271060876569e-16,1.1023581510449902,2.6603972637981608e-17,-2.4663916652348883,0.5047187705067517,0.026962212543529037,1.0006059944066537,0.02571854685831731,0.023562544571124427,0.002877938002773262,0.2161085252295913,0.013317096119719353 +104536.81282712151,14994.0,14993.997986296345,5.649645367970281e-16,-0.1890146539982087,2.6640596887369527e-17,-2.445895378373151,1.2164245981656017,0.02803677406579858,0.9871622217101877,-0.5497357570786212,0.02335263465129891,-0.011401907900413608,-0.8561347270640302,0.01331788974325867 +104885.66713609539,15058.0,15057.999302468945,5.85376938187508e-16,0.5763874122125322,2.667812978163725e-17,-2.4783224283833145,0.061953861563189286,0.027077757130190266,1.011221208894034,0.4738519298700737,0.023680834004645206,0.0027051157314355177,0.20430259331033007,0.01324073124870481 +103965.1770144278,14841.0,14840.999947621616,5.441964953783418e-16,-0.9786223699914871,2.636717227492218e-17,-2.4240347973833023,2.0660578382534425,0.027087916698405833,1.0048749802449943,0.2084331684438968,0.02338869711279427,0.004696774739682412,0.3576192508777786,0.013133450529170757 +104657.89204800228,15025.0,15024.997235329069,5.291370559902183e-16,-1.5604400748894864,2.618680759828323e-17,-2.4474869070028125,1.119624494026985,0.029039283412107895,1.0070782890882277,0.3024366368726101,0.02340420513011186,-0.005298128388872043,-0.4048299005085303,0.013087295138567476 +105104.3320267707,15071.0,15070.991729175797,5.95019096110726e-16,0.9335257688335624,2.6800648622680576e-17,-2.48966742734985,-0.3400391003696072,0.028430340332455842,1.0279217001184597,1.1706788386728464,0.02385086259021603,0.01347745991263468,1.0255680015129696,0.013141459067318844 +103117.18933059537,14779.0,14778.999954219906,5.741061709636474e-16,0.1552774697240199,2.6444087290612467e-17,-2.5085237338167152,-0.9578959261781123,0.029777487341992842,0.9880764025321682,-0.5036977763199469,0.02367212647819584,-0.0007966132391824033,-0.05871465634982645,0.013567536433086148 +103861.51126322581,14909.0,14908.999988128835,5.141729840097847e-16,-2.147589053088375,2.5995204208147898e-17,-2.4563869001257395,0.8147506117659684,0.0289819971083903,1.001915565907135,0.08215253976479471,0.023317184260150284,-0.005365587573949675,-0.4086698661866238,0.01312939371825232 +104401.21304313006,15015.0,15014.998970160144,5.994262650853293e-16,1.0978905429865868,2.6802549009376762e-17,-2.4967114366338805,-0.5986746885080364,0.02791405241388654,0.9947189099510139,-0.22577026819208507,0.023391432766040707,-0.004090920537099004,-0.30811253302204894,0.013277358428020362 +103896.10903963698,14907.0,14906.999883600616,5.868618849352218e-16,0.6376991535457983,2.6441755240640774e-17,-2.479523442602133,0.017226036087355272,0.027664948305589025,0.9909518757304067,-0.3908429511409154,0.023150281316781655,-0.002695578042956474,-0.20441807152789476,0.013186593645115359 +103877.93353330657,14893.0,14892.999790511029,5.558775494680619e-16,-0.5356629666367859,2.6364433256618963e-17,-2.4743403482955824,0.2062639870749244,0.02743887473852516,1.0021831242943509,0.09496376921188177,0.022989023208207993,0.0008862487268151458,0.06861458581377454,0.01291633136459492 +103396.93676252602,14770.0,14769.999940468431,5.867378273072098e-16,0.6352366873526242,2.6348961954583328e-17,-2.5155847264284814,-1.2077020193866677,0.029464823157746442,1.015575822278027,0.6678015218032541,0.023324029325311725,0.018553263639159025,1.4254318720711299,0.013015889431601825 +104160.03334950031,14964.0,14963.999998356765,6.10907582479636e-16,1.520382351595521,2.690611505501016e-17,-2.470844076924422,0.34594540723430006,0.026466381354145058,0.9927990773010172,-0.3077142659053865,0.023401328754763813,-0.0006925360964855318,-0.05198559811791703,0.01332169142143325 +104864.6850193089,15060.0,15059.999218909168,5.609746644889459e-16,-0.3411527664813018,2.6455407658400925e-17,-2.4647004340264353,0.5481903216886338,0.027909223071352622,1.0027414961956898,0.11570915352832659,0.02369299326883945,-0.005698317662346879,-0.4266719822022961,0.013355265637398148 +104923.17543779277,15101.0,15100.995028939238,6.127117410688088e-16,1.584235014851103,2.69604829260911e-17,-2.503434209444339,-0.8389182406124437,0.02793384183330101,1.0016273884843798,0.06972600471156862,0.02333976385298056,-0.0021331039675606374,-0.1622460782952347,0.013147337611939605 +104453.17589728598,14909.0,14908.97953769627,5.536892055352507e-16,-0.6201645362100462,2.630075328787414e-17,-2.4540288821633847,0.9320828943896502,0.02786352801123104,1.0282733858014168,1.2048710613023608,0.023465901630051383,0.016505203933255115,1.281884889389383,0.01287573016101099 +104688.55914493323,15031.0,15030.991877826402,6.006938103496527e-16,1.1617511752543976,2.6420296362456006e-17,-2.499085904547056,-0.7049249572708577,0.027075087000675516,1.017186344545262,0.733530932235582,0.02342961119973933,0.009323142088312936,0.7143190125704435,0.013051790480508207 +104926.93916741286,14995.0,14994.999995528324,5.463460924576351e-16,-0.9060463828756001,2.610672918012471e-17,-2.5100838929434004,-1.0195942133745202,0.029505750963251018,1.0595024692788106,2.4864521555702046,0.023930671316362114,0.029147843606924958,2.2708541683724275,0.01283562987570262 +104418.96129561374,14975.0,14974.999993368268,5.131487853632615e-16,-2.168120600530305,2.6221426346317257e-17,-2.4229966833284244,1.9386863136486103,0.029403063440570344,1.021692983447443,0.9054165163310347,0.023959120533109186,0.004508823226188785,0.3403330078401741,0.013248268966923719 +104597.41903067115,15024.0,15023.998585086974,5.5052416135314e-16,-0.7400989417831648,2.6315182399714976e-17,-2.4745524720532823,0.18891505066199857,0.028835859967897517,0.9960930724767406,-0.1660999857570656,0.023521540386966532,-0.009836850928478167,-0.7370514625111277,0.013346219943671375 +105023.41205975645,15104.0,15103.99832493881,5.797871057178e-16,0.36796529273395895,2.659790450638014e-17,-2.4868954513434605,-0.2490035791048457,0.02769217763153948,0.9920885937026824,-0.3372923372634315,0.023455636026319258,-0.01239104177720175,-0.9265919584774807,0.013372705929330483 +103289.81254198536,14773.0,14772.998796750477,5.619254179684634e-16,-0.3038067169570706,2.6578023397282623e-17,-2.482157660822906,-0.07419044378169867,0.029082732396842604,1.0184661779445874,0.779855959686005,0.02367895983255992,0.017436881008395338,1.3229099579545722,0.0131807013043847 +102660.54003349427,14770.0,14769.999999705335,5.364790752662568e-16,-1.2975231162634542,2.5834549160307202e-17,-2.466043384426776,0.48372644345193777,0.028852289888531828,0.958831803268213,-1.8022271411880775,0.02284295680102107,-0.023476141574471938,-1.7556113526783363,0.01337206070048309 +105006.72279247809,15048.0,15047.99465532865,5.993539114702324e-16,1.1005580354163806,2.6671843306406585e-17,-2.466000437205987,0.5172157125141147,0.027067164541391957,1.018356543752431,0.7674507997149594,0.023918854158792646,0.00983138430351044,0.737755349439672,0.013326076606530082 +103102.28361685216,14716.0,14715.999442685392,5.674552108899239e-16,-0.09631102227788653,2.6422615500161828e-17,-2.454995788646379,0.9125078449265976,0.027401639879197092,0.9984852826342058,-0.06491044261688289,0.0233354958728976,0.009111631303170683,0.6896669209855639,0.013211640323635889 +104531.59108373195,14950.0,14949.998318397973,5.833625674801249e-16,0.49877785891182885,2.679061879225682e-17,-2.4544331622678186,0.9389281566958414,0.027229812579221184,1.01219836264865,0.514842214521027,0.023693400239136235,0.008740556642838928,0.6596348798548487,0.013250598035025481 +103466.32364242684,14842.0,14841.999992364486,5.559773024581159e-16,-0.5338636010041037,2.626644243119372e-17,-2.5146215085708925,-1.1450405559233547,0.030236054427761224,1.009723793948652,0.4124474334233883,0.023575838181226716,0.00721903046759916,0.545000471514331,0.013245916003596215 +104316.6824251344,14911.0,14910.975299452331,5.863638892368771e-16,0.6201956888974903,2.638504189857696e-17,-2.4780783592346936,0.07238488048354566,0.026547543526624376,1.0340466137193263,1.440548699949866,0.023634476030217625,0.023923346759196155,1.8463064644490832,0.01295740832838096 +105529.5151767736,15190.0,15189.99957462132,5.834717114215463e-16,0.5010886028147278,2.6884888911606908e-17,-2.4889834457321363,-0.3136735758566292,0.028639472443934322,1.0016659371683858,0.07092028990347908,0.02349027578219314,-0.009923624425733822,-0.7499001090737418,0.013233261744675887 +105457.94736951304,15122.0,15122.00000033066,6.027761674204073e-16,1.228476892803507,2.668032879772676e-17,-2.484405248255628,-0.16595468883164577,0.026544885755513475,1.033269220492325,1.3891723116963275,0.0239489516255185,0.015526648583225076,1.1821004068806855,0.01313479675063867 +104470.87925587903,15024.0,15023.99640239301,6.03790838846066e-16,1.2589189731428367,2.684115464691797e-17,-2.483869890747045,-0.14371975264646822,0.026926644916821823,0.9864947020215399,-0.5740287753518823,0.023527214241448623,-0.008837024981601296,-0.655230584788196,0.013486893296438342 +105162.29126301195,15112.0,15111.999999317562,6.177356170407144e-16,1.78390671313331,2.6759032122744834e-17,-2.506380869101536,-0.9667329798802878,0.027288682242746005,1.0253507279309722,1.0574779091807105,0.023972820340628042,0.0127067997195174,0.9573832730291858,0.013272427122434208 +103777.57790654837,14877.0,14876.995407234004,6.038262246294493e-16,1.268964516464256,2.6656556736275105e-17,-2.4620622880379215,0.6649687898115889,0.02697526897038421,0.9726737228224851,-1.1837182765463132,0.0230851189163385,-0.010148044300178573,-0.7591884031873629,0.013366964323444887 +104806.70057053298,15096.0,15095.997686243283,5.439260391738729e-16,-0.9824410031651785,2.6539976183937103e-17,-2.485384750079234,-0.18354769817720248,0.0293370613345176,1.0113235065508546,0.4749012032977535,0.02384392052962426,-0.0047982174071041815,-0.3597467862069802,0.013337763090796671 +105784.79092823583,15259.0,15258.997768580222,5.827961339350463e-16,0.47726383852805365,2.6811446629837503e-17,-2.5047295496039585,-0.8455934842339454,0.029245198863330865,1.0148334072555911,0.6252990634905237,0.02372210054623874,-0.005345331494155609,-0.40347832963608515,0.013248125367666706 +105162.46385599423,15152.0,15151.990215977432,6.104114704168049e-16,1.4901748418897318,2.711861003206714e-17,-2.474470457622988,0.20105659418141558,0.027502417413986915,0.9773365909881715,-0.9704918298581112,0.02335249851113323,-0.01955450863830202,-1.4522701704146073,0.013464787087597772 +103698.33193109326,14764.0,14763.999994387923,5.411019601513098e-16,-1.1049254973794107,2.615383563618418e-17,-2.4445338644781205,1.2428914744224622,0.028535182879389883,1.0311968370972786,1.3211549589914875,0.023613306588269496,0.023708273144118674,1.8300034868884438,0.01295531583080744 +104821.39827733487,15059.0,15058.999246115178,5.68594851432458e-16,-0.053139051470666006,2.6442861297923464e-17,-2.466495670942182,0.4820149547355415,0.02801641095394487,1.0120587400167682,0.5137233646888311,0.023473216998943287,0.0011567486920515345,0.08829982954016835,0.013100236977528022 +106505.23587817501,15264.0,15263.999999850623,6.218963100694398e-16,1.9102781400790376,2.7166886842610554e-17,-2.471373714397668,0.3260255898816594,0.026458921845561856,1.0244822100573021,1.0264179091613137,0.023852087769304933,0.006396545662894279,0.48529211649252724,0.013180815112195989 +103716.65615689153,15034.0,15033.999998061632,5.577294742966784e-16,-0.46570958681866953,2.634802042007193e-17,-2.482904112111916,-0.0985414772370263,0.029470961805563557,0.9613012436149456,-1.6877220023130035,0.022929579831286302,-0.03170581928856768,-2.3653482434347395,0.013404292317873425 +104225.27879959819,14900.0,14900.00000050222,5.199590965611082e-16,-1.9135528494747858,2.615078201400422e-17,-2.4556855368707633,0.8170335145337638,0.02975944400899101,1.0386428710542701,1.615496590731249,0.0239201192227701,0.018544221167335707,1.4217032000742367,0.01304366563032804 +102786.34653474073,14809.0,14808.997062384125,5.671608003596035e-16,-0.10806725632187256,2.6272524509552447e-17,-2.5099327217773117,-1.0258010259966965,0.02917985166590013,0.9698004119754605,-1.3141855751175404,0.02297969829857434,-0.014914062983158742,-1.1175505167317183,0.013345314381648704 +103693.80219318194,14925.0,14925.000000354028,6.262827577806683e-16,2.0904413060784135,2.6923864170218004e-17,-2.5411901412740177,-2.2024997893381664,0.027782132634121534,0.9887409816798305,-0.4825792434259602,0.023330921239460494,-0.000780752700560303,-0.05835466951287918,0.013379438305069774 +104450.08142122456,15026.0,15025.994475005335,5.564571488415918e-16,-0.5143813331337013,2.632842657003666e-17,-2.4774149392783564,0.09322806764920688,0.027728352488978595,0.9801244763987085,-0.8536014915849469,0.02328431217287018,-0.019296953745379125,-1.43996889663141,0.013400951777862312 +104737.74337394757,15065.0,15064.99996910083,5.554067901763217e-16,-0.5523022868838343,2.6422504795363372e-17,-2.4843390111399652,-0.14983936941816967,0.028957750935643706,1.0003050186229425,0.013048813191263154,0.02337520037046332,-0.008653486101637216,-0.657271863832087,0.013165763784843708 +104812.57801906214,15006.0,15005.99998335691,5.178722422960302e-16,-1.9922068346204702,2.6165836196370893e-17,-2.427920517426411,1.7570265877928588,0.029640691231093017,1.0336105370795834,1.4217166073837435,0.023640813439911784,0.0108890268433584,0.8436609848926396,0.012906874963222445 +106339.57096486905,15207.0,15206.999997505669,5.492624962630607e-16,-0.7780861160226986,2.6651939046210067e-17,-2.445831650210193,1.190829912702277,0.028692888401057177,1.0621195059244974,2.582425705182309,0.024054711738594643,0.02255630002191606,1.7588965807729846,0.012824119546587107 +103728.05730247138,14836.0,14835.999999972239,5.209387563455267e-16,-1.8911654411555279,2.5942333011593222e-17,-2.4645567717301,0.5223658413199294,0.029564008685709207,1.012433963029872,0.5286900139121922,0.02351843746369132,0.005379995631329906,0.4100847184331302,0.01311922973351953 +103639.17031866539,14849.0,14848.995763098657,5.526901347779302e-16,-0.6572633079617131,2.6336271951268132e-17,-2.4790756256950353,0.03172639935343961,0.029135808783939735,1.020468619932421,0.8708978450847256,0.023502894223408958,0.013820223566015184,1.0582294652684516,0.01305976068480504 +104906.65266527233,15057.0,15056.999676679454,5.898564205186925e-16,0.7409408234691087,2.6798928996413046e-17,-2.4884855686348057,-0.300606433374182,0.028228167107266427,0.9987704521192002,-0.052390350346280826,0.023468976112450243,-0.00464859897558333,-0.35058822561408415,0.013259426974311319 +103900.09432355365,14941.0,14940.999999124993,5.626494348423733e-16,-0.2802725807392193,2.622648686589148e-17,-2.5374766771077253,-1.956661148119976,0.029374875237314736,1.0177793643158666,0.7655517169966314,0.023224249807207744,0.008053383365417891,0.624695729612739,0.01289168947962929 +105421.80459221234,15086.0,15085.991362965473,5.516294154681729e-16,-0.6898096398915645,2.6631382731495095e-17,-2.432532234555485,1.6716428368356826,0.028395877635182285,1.0247181424780918,1.0377693753611312,0.023818531424180996,0.006247790610896622,0.47575967680589687,0.013132240741465004 +104763.55636490823,14971.0,14970.992992053427,5.754182702542043e-16,0.20500866554129127,2.6429469407540822e-17,-2.471106898042058,0.3212152695157941,0.02768580077574692,1.0213216150467694,0.912626363652129,0.023362918162308007,0.01187716446399299,0.918512218279798,0.012930872586797705 +104010.23755771558,14983.0,14982.986816532728,5.82399109476417e-16,0.4663708032910838,2.658637588142092e-17,-2.487410895062763,-0.26203121450146666,0.02828248946165792,0.9697278266558659,-1.3016693628373415,0.02325642302754035,-0.02063547654962738,-1.5264684883010589,0.01351844254092295 +103802.99200127616,14915.0,14914.999158519993,5.98955655549965e-16,1.0774187851228139,2.687502385311054e-17,-2.492051222393842,-0.42355797140608237,0.028452356483424656,0.9793257627282156,-0.8930739375796235,0.023149524806215908,-0.00909363354268816,-0.6818341388097926,0.013337017061894227 +104450.88752347854,15017.0,15016.999255893796,5.809988580730963e-16,0.41359203147894474,2.65934960926737e-17,-2.4601309159680382,0.7212149320641501,0.02754946292514427,0.9867690971578991,-0.5699875949468021,0.02321261543128109,-0.010889108810203578,-0.821966945763293,0.013247623722985307 +103991.60794535361,14949.0,14948.999999199335,5.116779889999498e-16,-2.2436615470269636,2.5994121563179454e-17,-2.4490704069391716,1.0231506789081568,0.0302297537385545,0.9920210767039928,-0.34551318825890465,0.02309296306810833,-0.01362458177359732,-1.0422571545953772,0.013072188292040678 +104522.21161729087,14952.0,14951.996157985375,5.778641034812695e-16,0.29445980626519813,2.670688261672947e-17,-2.463321751237089,0.6036096703921505,0.027630850831258973,1.0145116280579276,0.6221761632371783,0.02332398589882278,0.00909011274137726,0.7004224247969091,0.012978043562801382 diff --git a/test/data/sensitivity.dat b/test/data/sensitivity.dat new file mode 100644 index 0000000..8872926 --- /dev/null +++ b/test/data/sensitivity.dat @@ -0,0 +1,11 @@ +loge,emin,emax,crab_flux,photon_flux,energy_flux,sensitivity +-1.51402150412,0.02,0.0468734582318,0.144989163852,6.60543533444e-10,3.05652923293e-11,3.56378721075e-11 +-1.14412450369,0.0468734582318,0.109856054331,0.0561134683318,7.24745374549e-11,7.85976500012e-12,9.16416230711e-12 +-0.774227503252,0.109856054331,0.257466658709,0.0315037806261,1.15354211075e-11,2.93193651747e-12,3.41851723554e-12 +-0.404330502818,0.257466658709,0.603417633655,0.0244879762385,2.54200331885e-12,1.51423828719e-12,1.76553948308e-12 +-0.0344335023848,0.603417633655,1.41421356237,0.022256802649,6.54995675049e-13,9.14436478628e-13,1.0661952755e-12 +0.335463498049,1.41421356237,3.31445401734,0.0250217133244,2.08758967789e-13,6.83057479879e-13,7.96416891675e-13 +0.705360498482,3.31445401734,7.76799609716,0.042333217441,1.00129446033e-13,7.67840339358e-13,8.95270214276e-13 +1.07525749892,7.76799609716,18.2056420303,0.0929194773185,6.23074149979e-14,1.11981387734e-12,1.30565686449e-12 +1.44515449935,18.2056420303,42.6680700645,0.14972962577,2.8463818009e-14,1.19893673197e-12,1.39791085453e-12 +1.81505149978,42.6680700645,100.0,0.382388032371,2.06082989569e-14,2.03442769003e-12,2.37205890421e-12 diff --git a/test/data/spectrum.fits b/test/data/spectrum.fits new file mode 100644 index 0000000..80da2e5 --- /dev/null +++ b/test/data/spectrum.fits diff --git a/test/data/ts.dat b/test/data/ts.dat new file mode 100644 index 0000000..573ca7f --- /dev/null +++ b/test/data/ts.dat @@ -0,0 +1,101 @@ +TS,LogL_bgm,LogL_all,Nevents,Npred_bkg,Npred_all,Crab_Prefactor,Unc_Crab_Prefactor,CTABackgroundModel_Prefactor,Unc_CTABackgroundModel_Prefactor,CTABackgroundModel_Index,Unc_CTABackgroundModel_Index +21858.52363364247,164117.78313692164,153188.5213201004,22947.0,22946.99999940957,22946.999182570653,5.793209910482599e-16,9.843787551860997e-18,0.9961007850466849,0.01197837072408249,-0.002988471615299611,0.007462197925324776 +22254.447332882788,165680.6922036114,154553.46853717,23099.0,23098.999999613512,23098.998093109363,5.834991758798705e-16,9.870580046002392e-18,1.0126678419515491,0.011967610062684665,0.004782477384744285,0.007316347479857281 +21669.561687524547,165499.35987704716,154664.57903328488,23134.0,23133.99999969615,23133.999315025623,5.759131518978944e-16,9.8078414532141e-18,1.0105363690914626,0.011967661087729033,-0.00017543331960004447,0.0073220367405137776 +20586.996292978816,164330.2704312009,154036.77228471148,22897.0,22896.999999836524,22896.997968569693,5.555266953291272e-16,9.65303499566792e-18,1.0126646529020689,0.011972560615151874,0.005782431464096715,0.007323711254845693 +21217.911190702056,163471.7261244653,152862.77052911426,22869.0,22868.999999925014,22868.99897334564,5.70225985492044e-16,9.781422622377391e-18,0.9910189865469351,0.011911721710682954,-0.006124951668038474,0.007435839123654169 +21499.99581397162,164613.52000714888,153863.52210016307,22972.0,22971.999999709296,22971.999752804753,5.6910222413269675e-16,9.758624936302264e-18,0.996047801086967,0.011904474629350275,-0.006618987374187345,0.007377903239793396 +21218.138624030922,164860.51429341873,154251.44498140327,23028.0,23027.999999721134,23027.999861910866,5.631722425238875e-16,9.709038606307067e-18,1.0013335385708708,0.011894934859389805,-0.006251748911501181,0.007319737866598979 +21516.115901010286,164574.57805417257,153816.52010366743,22977.0,22976.999999663327,22976.99992616064,5.737703984639025e-16,9.8058501757969e-18,1.0017304453830707,0.011926394544762287,-0.0012345022004066025,0.007362686667357548 +21814.264059167763,163757.54297349413,152850.41094391025,22865.0,22864.999999773583,22864.999780063565,5.781984731779256e-16,9.833138684097321e-18,0.9972279507003987,0.01184624889933737,0.0009075599364746032,0.007331742276660691 +21722.886111398868,165089.47672604493,154228.0336703455,23109.0,23108.999999674452,23108.9997904007,5.709992445993981e-16,9.76964807969154e-18,0.9967680654306393,0.011876107200647554,-0.011021721057770708,0.007332464028069347 +21598.863410390506,165201.15784713515,154401.7261419399,23121.0,23120.999999727013,23120.99956818323,5.796621137864304e-16,9.85685177987898e-18,1.0005077278386079,0.011874257987756785,-0.006437475225365512,0.007308543239518941 +21794.114624339214,164895.88104651406,153998.82373434445,23001.0,23000.99999958087,23000.999482800085,5.798901936139489e-16,9.862163118374043e-18,1.005069551247296,0.011900250356816858,0.0019410618022905287,0.0073140103821301514 +21227.425172272255,166038.31113511338,155424.59854897726,23106.0,23105.999999409567,23105.996163342123,5.683306232669589e-16,9.767724432554512e-18,1.0208070375514615,0.012000389287969418,0.006896248020332655,0.007278861707651359 +21573.883430126647,163893.0064978525,153106.06478278918,22893.0,22892.999999855037,22892.997657277247,5.746330138649936e-16,9.809856798510093e-18,0.9887958563489442,0.011919769035438853,-0.007702982005473525,0.007460539914015142 +21216.534766604716,164024.07526321124,153415.80787990888,22891.0,22890.99999976823,22890.999903703054,5.634003761042508e-16,9.713723962581571e-18,0.9988843468354479,0.011884925398741884,-0.002623626233923857,0.007345635611935417 +21597.08751892927,165975.3560595174,155176.81230005276,23150.0,23149.99999944865,23149.998888937902,5.718580414582826e-16,9.78000524693473e-18,1.0132558638441092,0.011992693047890666,0.0002528999590185808,0.007321935314776597 +21387.703713911586,163364.79937473897,152670.94751778318,22794.0,22793.999999815976,22793.999973829385,5.6892940405553605e-16,9.758446104823037e-18,0.9943975957184543,0.011824486448349799,-0.0008002832127617221,0.007333535186744982 +20912.907039812184,162038.81793470908,151582.36441480299,22666.0,22665.99999997628,22665.990260466566,5.601818185889693e-16,9.694320601561969e-18,0.9751815451597996,0.011793838389151424,-0.01294140924200319,0.007455055791860874 +21911.95291238872,163189.5450744877,152233.56861829333,22810.0,22809.999999897966,22809.998057095352,5.791145848008091e-16,9.84239549367057e-18,0.9854579489782141,0.011857156203560955,-0.005834836728800691,0.007437257754506665 +21121.265116654686,164055.15655683712,153494.52399850977,22912.0,22911.999999787822,22911.99984896773,5.654325050779409e-16,9.741375490234794e-18,0.9999753099585894,0.01194854296456228,-0.0021040061041819437,0.00739916979858624 +21195.306229041773,164281.77518035137,153684.12206583048,23034.0,23033.999999960743,23033.998176414847,5.70433524603058e-16,9.77671925031924e-18,0.9886794379131255,0.011878644954228475,-0.01446464700955386,0.007402981305713397 +20846.777400515275,162770.057188395,152346.66848813737,22731.0,22730.999999907664,22730.99981226707,5.624363678525482e-16,9.712847646970035e-18,0.9909653476365609,0.011763408213370439,-0.002595590722067562,0.007301229915792391 +20897.083739825,164019.2197518424,153570.6778819299,23003.0,23002.999999907483,23002.999556953055,5.649662099748417e-16,9.736547752010038e-18,0.9935802602295315,0.01188380653681554,-0.010805113941356031,0.007371306800476196 +20602.686458923854,164911.49560107855,154610.15237161663,22935.0,22934.99999955501,22934.99341642103,5.530693214901318e-16,9.625543333068213e-18,1.0206021943272432,0.011985188441572748,0.00972446272247772,0.00727243068453318 +21836.23271708103,163824.30785453672,152906.1914959962,22755.0,22754.999999495434,22754.999710735854,5.75568768540384e-16,9.811525606728942e-18,1.0043798388913345,0.01188200956164694,0.010276896104008231,0.007325165400710747 +21420.40131879633,164409.64030012718,153699.43964072902,23001.0,23000.999999780248,23000.99928771896,5.684570967501337e-16,9.762208897347998e-18,0.9913746713865513,0.011878388317921098,-0.011560419259827103,0.007384451973090443 +21727.81684923329,164232.132445953,153368.22402133635,22985.0,22984.99999976957,22984.998900142164,5.738813127112364e-16,9.809036890551794e-18,0.9911339994790322,0.01189255974895351,-0.009747052205701144,0.00740674910573882 +20868.72095527401,164301.13159811447,153866.77112047747,22992.0,22991.999999802247,22991.99969008748,5.605800164321631e-16,9.703156512461817e-18,0.9973271760042328,0.011871321613801066,-0.008569072284481264,0.0073309043881014035 +22492.009746829106,164682.28638778496,153436.2815143704,22951.0,22950.9999996003,22950.99665538983,5.879507535878784e-16,9.910190021652269e-18,1.001162832181959,0.011846315465073663,0.002987242954331668,0.0073024800531300005 +22653.864335719903,165519.4792556914,154192.54708783145,23076.0,23075.999999249838,23075.993670520074,5.912861048051989e-16,9.940531531808007e-18,1.0024002934241467,0.011941973255326052,-0.00023952533210942385,0.007373252448619404 +21449.306879230775,166768.10502219904,156043.45158258366,23297.0,23296.99999917116,23296.992123045784,5.7341518582148665e-16,9.807974620497261e-18,1.024760753033504,0.012028256916479971,0.0035902945030410723,0.007258563117465229 +21486.346823231317,164876.94234162793,154133.76893001227,23003.0,23002.999999741518,23002.999859281495,5.757455861467064e-16,9.822563593405131e-18,1.004168442832293,0.011942314003147688,0.00011042172295081477,0.0073583821510989885 +20652.239219280833,163809.36895235302,153483.2493427126,22957.0,22956.999999943808,22956.995521902107,5.528577248898855e-16,9.622644416274187e-18,0.9853639781514777,0.011834892876956514,-0.018391417942683855,0.0073803794078678014 +21787.521341275133,164528.73527216076,153634.9746015232,22991.0,22990.99999973,22990.998358757795,5.720325280614957e-16,9.785480546850936e-18,0.9855686593361147,0.011839563015884493,-0.014781546186952868,0.007392004202286942 +21402.38363453705,164305.6816226261,153604.4898053576,22994.0,22993.99999989152,22993.99965618869,5.712540995621592e-16,9.778682761426494e-18,0.9949474763637698,0.011898956019282798,-0.007806261712544803,0.007380964683300804 +20141.589426982973,163856.60168552134,153785.80697202985,22929.0,22928.999999913434,22928.999999743348,5.438078945642788e-16,9.558024400005279e-18,0.9932852802680611,0.011875856351872056,-0.01340284403577944,0.007359275829012058 +21943.395299219526,163065.74164094013,152094.04399133037,22784.0,22783.999999305994,22783.99675835668,5.79676765978459e-16,9.848370996477712e-18,0.9837834848183966,0.01187815045838842,-0.005954920583608014,0.007472494212571003 +20834.792521972733,164064.23789739315,153646.84163640678,22895.0,22894.99999975629,22894.999530192188,5.589302026667169e-16,9.677416099688684e-18,0.9986526033294748,0.011854084906733467,-0.004096128598782811,0.007311666493293262 +21305.324882587476,163265.06469459928,152612.40225330554,22828.0,22827.999999891683,22827.99808304598,5.713178115397861e-16,9.79512015984894e-18,0.987396852431259,0.011894410275773492,-0.007029367080902169,0.007451877733083311 +20847.034227716562,163814.44659629313,153390.92948243485,22897.0,22896.9999996484,22896.99953025752,5.607597257232211e-16,9.688770805376503e-18,1.0047595349843508,0.011799861066908465,0.0010146013322702991,0.007215496921685529 +21216.92292083666,164176.47242941064,153568.0109689923,22999.0,22998.999999831543,22998.9975955996,5.642001984733995e-16,9.72351332200732e-18,0.9899066275912457,0.011903062056425938,-0.013687904197762802,0.007417248945029567 +20816.64541355765,163181.06644880498,152772.74374202616,22777.0,22776.99999988913,22776.994775080595,5.516451209091623e-16,9.618999946029458e-18,0.9869057530226409,0.011856380365721023,-0.010358640899461943,0.007409856539614678 +21679.522176688828,163799.51890170656,152959.75781336214,22837.0,22836.99999978841,22836.99973823716,5.717636386575965e-16,9.776297868466177e-18,0.9908230482181739,0.011847216840415472,-0.0046007921454711475,0.007378504976736998 +21106.35513021017,164821.97117069492,154268.79360558983,22917.0,22916.999999453958,22916.99841731256,5.647103860714512e-16,9.727705266059672e-18,1.01541521494636,0.011955291437595085,0.00943775316656054,0.007291252547122076 +20837.813223739155,163866.034636364,153447.12802449442,22791.0,22790.999999713706,22790.998309820432,5.572360798379893e-16,9.654736904973028e-18,1.0091814893739082,0.011893607890048716,0.007811373787627806,0.007283734723217876 +20780.997344649804,164108.8221487588,153718.3234764339,22982.0,22981.999999861746,22981.99986802358,5.630081040612143e-16,9.720705407249754e-18,0.9960740698856181,0.011864470017298944,-0.008531796083847434,0.007334719933560889 +21835.693036562065,164481.7513907408,153563.90487245977,22863.0,22862.999999581152,22862.999821485086,5.719539799575279e-16,9.779590491901723e-18,1.0043623557350347,0.011875625884065816,0.0049348547201782614,0.007305172253576122 +21653.10258389666,165542.21112609853,154715.6598341502,23262.0,23261.999999841548,23261.99499680352,5.749522283549226e-16,9.819619041947758e-18,0.9902674093404891,0.011918436489673211,-0.02105969680201884,0.007407329162175303 +21846.154768515087,164624.27857300622,153701.20118874867,23012.0,23011.99999979313,23011.998528939348,5.833817292347365e-16,9.884966321550659e-18,0.9935369248121374,0.011882112655646851,-0.006556706861582095,0.007382216624463548 +22013.17601596308,166235.52968944985,155228.9416814683,23227.0,23226.999999677897,23226.998673585054,5.797707114245376e-16,9.848651501252573e-18,1.006518546940583,0.011988963639015245,-0.00598127326508908,0.007364840449676914 +21651.91926104814,163154.12158821282,152328.16195768875,22773.0,22772.99999985566,22772.998844177993,5.714396502849746e-16,9.769436981298296e-18,0.9851000302656024,0.011819450107849325,-0.006582743622206519,0.007400947537603796 +21989.977012775722,165301.46496142616,154306.4764550383,23054.0,23053.999999486277,23053.99866206013,5.831451508774712e-16,9.873936836478257e-18,1.0070990373467372,0.011907048508928606,0.0022074936263132647,0.007300838328239046 +21438.695088208653,163698.03314021966,152978.68559611533,22850.0,22849.999999890217,22849.999942469298,5.723101116737367e-16,9.79415042715683e-18,0.9959990059596918,0.011909255151422262,-0.0009492657684907494,0.007402029667490372 +20838.374803555547,163909.36890479646,153490.18150301868,22825.0,22824.999999715288,22824.99939814383,5.603315148685522e-16,9.69540670910627e-18,1.0087802770334235,0.011925882318992111,0.006917357818017513,0.007318975346866552 +21270.20557162707,165625.21774982946,154990.11496401593,23098.0,23097.999999612013,23097.99846084131,5.695134364689924e-16,9.766410999308169e-18,1.0165267626116348,0.012010348326125827,0.0042397629393144945,0.007319288376668166 +21527.705370273616,164243.86907532276,153480.01639018595,22965.0,22964.999999743457,22964.99993603194,5.701748498974925e-16,9.771290624802224e-18,0.995702026500132,0.011874291328625755,-0.006337078712909612,0.007353017018362593 +21055.645831402973,164272.05201691244,153744.22910121095,22913.0,22912.99999970991,22912.999921720784,5.630081490231017e-16,9.7255367108198e-18,0.9979621989156474,0.011868129706833142,-0.0043189322848428005,0.0073321006440042955 +21156.85824884911,162399.58053200433,151821.15140757977,22706.0,22705.99999991407,22705.997403443518,5.654797040066863e-16,9.740490054150095e-18,0.9784619159402815,0.011744743142957838,-0.010617863436450973,0.007382073296970704 +21201.688190487737,164789.7753571481,154188.93126190422,23026.0,23025.99999906124,23025.99972065303,5.667600645790645e-16,9.742584040372483e-18,1.0120410866982648,0.01207095002314141,0.002979507608702241,0.007418035938329006 +21579.612940530176,164202.44370278623,153412.63723252114,22874.0,22873.999999804666,22874.00001950805,5.716335760382397e-16,9.78228682145465e-18,1.002323437487339,0.011916380859040334,0.0028258681377384813,0.007357128369306409 +21460.62560378993,164011.29522604778,153280.9824241528,22951.0,22950.99999985298,22950.999919398004,5.718007013749168e-16,9.78668674780352e-18,0.9941518980348495,0.011860449976759028,-0.006570586126134151,0.007354999480418171 +20888.288302358415,164286.87337022196,153842.72921904275,22976.0,22975.99999977232,22975.998047214835,5.544737396797197e-16,9.641570891212174e-18,1.0007779663248277,0.011916634775209257,-0.006801761329085501,0.0073460207848977796 +20408.074956183205,164422.41163892217,154218.37416083057,22955.0,22954.999999835327,22954.995893393418,5.513465666103151e-16,9.625616231421954e-18,1.0058950024352238,0.011974149384797315,-0.0028091604581978514,0.007363395502675995 +21802.37343549641,164500.592398738,153599.4056809898,22961.0,22960.999999757954,22960.9989212055,5.832468881001641e-16,9.890309949834711e-18,1.0034751072565908,0.011921220309959827,0.0031744106933537897,0.0073533366536048545 +20876.450997727632,164582.49714246154,154144.27164359772,22993.0,22992.999999842265,22992.9998257974,5.648090586661093e-16,9.736675647490638e-18,1.0025225514093496,0.011933600093706552,-0.0035293846222004244,0.007354639908367345 +21706.519023373083,164473.15900975114,153619.8994980646,22978.0,22977.999999763233,22977.999793557632,5.759940859191368e-16,9.813811252134857e-18,0.9984861958740264,0.011920224946766407,-0.0032243711331649107,0.007381393296644344 +21656.740131683822,165712.62332411483,154884.25325827292,23262.0,23261.99999977213,23261.998757088804,5.797595166213485e-16,9.856541718951383e-18,0.9985566042942451,0.011918622406145106,-0.013482108055460828,0.007351626660765964 +21108.65355328942,164601.01057088916,154046.68379424445,23016.0,23015.99999983359,23015.999758528997,5.653188521209414e-16,9.736939795430797e-18,0.9959645629245217,0.011889610707927175,-0.00938547147416546,0.007358851980664596 +20647.59149535757,162669.5612982715,152345.7655505927,22734.0,22733.99999995898,22733.995853601802,5.5190264735296685e-16,9.620164290713148e-18,0.9903309324372139,0.01187206487521323,-0.005898045252815988,0.007403049951329795 +22290.401570641203,163883.56355593016,152738.36277060956,22782.0,22781.999999546068,22781.99916149079,5.824528487203314e-16,9.864372038097642e-18,0.9966719584506485,0.011871832280126494,0.004939494731895152,0.007375380023497896 +21624.4816204338,164856.28447951959,154044.0436693027,22996.0,22995.99999963064,22995.999887018945,5.729990754258132e-16,9.793654933564995e-18,1.0001246408759745,0.011920784556020848,-0.003433668355429103,0.007365033404665828 +21599.218091811636,162947.49765919705,152147.88861329123,22666.0,22665.999999721185,22665.9995259345,5.76124788921756e-16,9.820444219997191e-18,0.9997434910652248,0.011837575389850108,0.010453639932065167,0.0073277935176205644 +21352.83466984972,163983.0417331816,153306.62439825674,22833.0,22832.999999748183,22832.999677072614,5.611046191249927e-16,9.692013011636335e-18,0.9935935257673363,0.011855675126920649,-0.004991283544181841,0.007359215729362822 +22138.833458752255,164579.4216104209,153510.00488104476,22972.0,22971.99999954863,22971.999029847222,5.82889920739159e-16,9.860552383071182e-18,1.0016245840621145,0.011965267887535018,0.0012043877283222457,0.007408104140671816 +21460.73646459746,164025.76170268297,153295.39347038424,22919.0,22918.99999972083,22918.99998951889,5.739817859219702e-16,9.810075323312415e-18,0.9966380414703272,0.01185241980257522,-0.0028036440323137035,0.007332677333553906 +20272.70840003481,163913.44129798363,153777.08709796623,22905.0,22904.99999986819,22904.993807344785,5.491116393214731e-16,9.603700699255587e-18,0.9993535632641539,0.01192103384820017,-0.006429848120336535,0.007364414061681116 +22040.583799355198,164196.0694411831,153175.7775415055,22948.0,22947.999999762113,22947.998520762158,5.830729768890433e-16,9.878618700156061e-18,0.9889198731570649,0.01183479859649009,-0.00766945927540272,0.00737644490898403 +22189.423341836373,162999.79537715373,151905.08370623554,22762.0,22761.999999802018,22761.99812426319,5.840469854404995e-16,9.88689591899704e-18,0.9874701389124014,0.011840024598898219,-0.001045075500813165,0.007415569662489257 +21404.67148317158,163675.38150747615,152973.04576589036,22873.0,22872.99999985997,22872.999797008175,5.707784704943953e-16,9.776315626974606e-18,0.9932817805516895,0.011877174992673342,-0.004381770811720796,0.007385324275316955 +21680.787235579046,164409.20965859844,153568.81604080892,23000.0,22999.999999595042,22999.996041215003,5.756363373649316e-16,9.818959282992611e-18,0.9911211088710422,0.011963680831576883,-0.009911884410269314,0.007474654715810191 +21579.654693948105,164174.39769977125,153384.5703527972,22910.0,22909.999999767482,22909.999826840383,5.71987722050035e-16,9.791834857909584e-18,0.9969717770530057,0.011930237932024544,-0.0026897435561697674,0.007407652576797861 +21297.579771767836,163385.8273468829,152737.03746099898,22844.0,22843.999999879547,22843.9991035317,5.688808155458493e-16,9.7608380955326e-18,0.987985497494533,0.011852158347747057,-0.007831002182125445,0.007402005571818617 +22413.117900692392,165160.81083219431,153954.25188184812,23075.0,23074.99999963741,23074.9980466868,5.849262030225733e-16,9.88285352545655e-18,0.9941410946295832,0.011801812232478627,-0.008211194148192122,0.007291859455219944 +21122.470503843215,163334.07738528133,152772.84213335972,22879.0,22878.999999898744,22878.99709745025,5.633723709743721e-16,9.722100961324131e-18,0.984984844375672,0.011857253919565828,-0.01296529513990941,0.007421939723453924 +21282.707195805793,163941.15171670247,153299.79811879958,22909.0,22908.9999998743,22908.99996153768,5.729219609420647e-16,9.797812724663524e-18,0.9954654665694013,0.011861939156832611,-0.003583079120707077,0.007352470265065281 +20724.34760198067,163318.72421718622,152956.55041619588,22790.0,22789.999999879707,22789.99813018993,5.550590116159096e-16,9.645783754938671e-18,0.9933982950272806,0.01187452330418287,-0.00495332364849543,0.0073808910696461375 +21757.13207008707,165179.7837367039,154301.21770166035,23030.0,23029.999999515534,23029.999095339273,5.763186525066147e-16,9.828374022376591e-18,1.0136834107892345,0.011971605614733427,0.006510480574850554,0.007314863606454032 +21671.1096591925,166527.63624669318,155692.08141709692,23264.0,23263.99999942336,23263.992478841574,5.84129987251449e-16,9.896922008537031e-18,1.0234213270128842,0.012014512018187535,0.006575992519346971,0.007267333319003582 +21221.41757199116,163903.57706571487,153292.8682797193,22906.0,22905.999999861677,22905.998783516217,5.646919165369301e-16,9.724700045763518e-18,0.9900591679883312,0.011885002043500589,-0.009753133051711622,0.007408431197796472 +21681.41252249159,165556.59745661478,154715.89119536898,23196.0,23195.999999720592,23195.99929739878,5.773631747722166e-16,9.836520659176885e-18,0.9977101559045379,0.01186422662665964,-0.012136729410717443,0.0073081040770472065 +21153.477707419137,166277.3899099234,155700.65105621383,23255.0,23254.999999606127,23254.99604384548,5.653184419662256e-16,9.730321230819675e-18,1.019156144918438,0.01201068815115771,-0.0010363684582943512,0.007280823299749718 +22226.132660409086,164336.5451073256,153223.47877712105,22935.0,22934.999999669544,22934.9976589201,5.8647466247030845e-16,9.907864241886007e-18,0.9988210866347508,0.011904203335735244,0.001435599985409938,0.0073741953380605406 +21424.358322370157,164683.5971306915,153971.4179695064,22940.0,22939.99999953921,22939.999271899742,5.679278962446394e-16,9.75563702408002e-18,1.0098384925060064,0.011899202395887487,0.0050329170142471,0.007276269049906902 +21455.353563566052,164947.2565606395,154219.57977885648,22987.0,22986.99999956718,22986.998933887335,5.706645322939917e-16,9.770262988270668e-18,1.0142244697962193,0.011960235218104793,0.007217238595620098,0.007300805026314625 +21452.114580445574,164273.43465029492,153547.37736007213,22845.0,22844.9999995031,22844.999775048265,5.677909293087916e-16,9.751768271419505e-18,1.0082800787052053,0.011922505792565435,0.007632916428577924,0.007322032069579832 +21019.69961927872,162935.32861177463,152425.47880213527,22657.0,22656.999999781452,22656.999043407268,5.573701521579324e-16,9.661153343333429e-18,0.994168966337365,0.011892058320532648,0.0016080034919236333,0.007408684998744538 +21441.518593108456,164053.07922474376,153332.31992818954,22960.0,22959.99999984039,22959.999974187285,5.727077800111475e-16,9.793315534540368e-18,0.9933895462403782,0.011810199688400778,-0.0072931345463456265,0.007310683986440049 +21781.365669648803,164462.0895137939,153571.4066789695,22877.0,22876.999999410484,22876.999794261017,5.738537591812932e-16,9.802523721098818e-18,1.0067930568779322,0.0119066595609444,0.006739517264159702,0.007318677995153281 +21013.571768569178,162221.01090832267,151714.22502403808,22714.0,22713.999999970183,22714.000000122807,5.578526953798294e-16,9.670657005596775e-18,0.96868333949408,0.011813808188693846,-0.02057929881999145,0.007517120515436912 +20965.61130349373,164112.7074564671,153629.90180472023,22893.0,22892.999999919695,22892.99927046633,5.5999754133615e-16,9.687391521031491e-18,0.9954549917387471,0.011910503281561615,-0.0062261048953279885,0.007392491292287055 diff --git a/test/dev/README.md b/test/dev/README.md new file mode 100644 index 0000000..0de5415 --- /dev/null +++ b/test/dev/README.md @@ -0,0 +1,30 @@ +Developer test scripts +====================== + +This folder contains scripts that are used for various testing purposes +during code development. None of the scripts actually needs to work, they +may actually be outdated, but could still be useful for further testing +if needed. + +`benchmark_cta_analysis.py` +> This script peforms a benchmark of the various CTA analysis types +> (unbinned, binned, stacked), and if matplotlib is installed, creates +> a plot of the benchmark results. + +`check_models.py` +> Checks a variety of models using an unbinned analysis. + +`make_pull_at_sensitivity_limit.py` +> This script creates pull distributions as function of energy +> for a test source that is at the 5 sigma sensitivity limit for +> 50h of effective observing time. It takes as input file the +> output file generated using cspull. +> If the Python processing module is installed the script can make +> use of multiple cores and runs in parallel. + +`make_survey.py` +> This script simulates a CTA galactic plane survey. + +`simulate_events.py` +> Simulates events and compares the simulation to the model. + diff --git a/examples/benchmark_cta_analysis.py b/test/dev/benchmark_cta_analysis.py index 921ec75..5fa01d4 100755 --- a/examples/benchmark_cta_analysis.py +++ b/test/dev/benchmark_cta_analysis.py @@ -299,8 +299,7 @@ def stacked_pipeline(model_name, duration): # Main routine entry point # #==========================# if __name__ == '__main__': - """ - """ + # Initialise parameters title = "CTA analysis benchmark" model_name = "${CTOOLS}/share/models/crab.xml" diff --git a/examples/check_models.py b/test/dev/check_models.py index 936d61d..2a4902b 100755 --- a/examples/check_models.py +++ b/test/dev/check_models.py @@ -1,14 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script performs an unbinned maximum likelihood analysis for a -# variety of models. This allows checking the models. +# Performs unbinned maximum likelihood analysis for various models # -# Usage: -# ./check_models.py -# -# ========================================================================== -# -# Copyright (C) 2011-2015 Juergen Knoedlseder +# Copyright (C) 2011-2016 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 @@ -33,26 +27,26 @@ import ctools # ================= # def pipeline(model_name): """ - Unbinned analysis pipeline - keep intermediate results in memory. + Unbinned in-memory analysis pipeline This function implements an analysis pipeline that successively calls ctobssim, ctselect and ctlike without saving the intermediate results as FITS files on disk. All data is only hold in memory. """ # Set script parameters - caldb = "prod2" - irf = "South_50h" + caldb = 'prod2' + irf = 'South_0.5h' ra = 83.63 dec = 22.01 - rad_sim = 10.0 + rad_sim = 5.0 tstart = 0.0 - tstop = 1800.0 + tstop = 180.0 emin = 0.1 emax = 100.0 rad_select = 3.0 # Write model name - print("*** Model: "+model_name+" ************************************") + print('*** Model: '+model_name+' ************************************') # Initialise timing wall_seconds = 0.0 @@ -60,27 +54,27 @@ def pipeline(model_name): # Simulate events sim = ctools.ctobssim() - sim["inmodel"] = model_name - sim["caldb"] = caldb - sim["irf"] = irf - sim["ra"] = ra - sim["dec"] = dec - sim["rad"] = rad_sim - sim["tmin"] = tstart - sim["tmax"] = tstop - sim["emin"] = emin - sim["emax"] = emax + sim['inmodel'] = model_name + sim['caldb'] = caldb + sim['irf'] = irf + sim['ra'] = ra + sim['dec'] = dec + sim['rad'] = rad_sim + sim['tmin'] = tstart + sim['tmax'] = tstop + sim['emin'] = emin + sim['emax'] = emax sim.run() # Select events select = ctools.ctselect(sim.obs()) - select["ra"] = ra - select["dec"] = dec - select["rad"] = rad_select - select["tmin"] = tstart - select["tmax"] = tstop - select["emin"] = emin - select["emax"] = emax + select['ra'] = ra + select['dec'] = dec + select['rad'] = rad_select + select['tmin'] = tstart + select['tmax'] = tstop + select['emin'] = emin + select['emax'] = emax select.run() # Perform maximum likelihood fitting @@ -98,39 +92,37 @@ def pipeline(model_name): # Main routine entry point # #==========================# if __name__ == '__main__': - """ - Perform unbinned analyses for various models. - """ + # Initialise flags need_help = False # Test for command line arguments print(sys.argv[0]) if (len(sys.argv) > 1): - if sys.argv[1] == "-h": + if sys.argv[1] == '-h': need_help = True else: need_help = True # Print help if needed and exit if need_help: - print("Usage: check_models.py [OPTIONS]") - print(" -h Display this usage message") + print('Usage: check_models.py [OPTIONS]') + print(' -h Display this usage message') sys.exit() # Dump header - print("*******************************************") - print("* Check models using an unbinned analysis *") - print("*******************************************") + print('*******************************************') + print('* Check models using an unbinned analysis *') + print('*******************************************') # Perform analysis for Crab model - pipeline("$CTOOLS/share/models/crab.xml") + pipeline('$CTOOLS/share/models/crab.xml') # Perform analysis for disk model - pipeline("$CTOOLS/share/models/disk.xml") + pipeline('$CTOOLS/share/models/disk.xml') # Perform analysis for Gaussian model - pipeline("$CTOOLS/share/models/gauss.xml") + pipeline('$CTOOLS/share/models/gauss.xml') # Perform analysis for shell model - pipeline("$CTOOLS/share/models/shell.xml") + pipeline('$CTOOLS/share/models/shell.xml') diff --git a/cscripts/cta_kb_bkgrates.py b/test/dev/cta_kb_bkgrates.py index 559e73d..559e73d 100755 --- a/cscripts/cta_kb_bkgrates.py +++ b/test/dev/cta_kb_bkgrates.py diff --git a/cscripts/cta_kb_root2caldb.py b/test/dev/cta_kb_root2caldb.py index fca22eb..fca22eb 100755 --- a/cscripts/cta_kb_root2caldb.py +++ b/test/dev/cta_kb_root2caldb.py diff --git a/cscripts/cta_make_nodes.py b/test/dev/cta_make_nodes.py index 9d8a3dd..9d8a3dd 100755 --- a/cscripts/cta_make_nodes.py +++ b/test/dev/cta_make_nodes.py diff --git a/cscripts/cta_root2irf.py b/test/dev/cta_root2irf.py index e0f6654..e0f6654 100755 --- a/cscripts/cta_root2irf.py +++ b/test/dev/cta_root2irf.py diff --git a/examples/make_pull_at_sensitivity_limit.py b/test/dev/make_pull_at_sensitivity_limit.py index 950cb55..319b903 100755 --- a/examples/make_pull_at_sensitivity_limit.py +++ b/test/dev/make_pull_at_sensitivity_limit.py @@ -135,11 +135,7 @@ def create_model(flux, index=-2.48, fitidx=False): # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Create pull distributions for the energy bands and fluxes given in the - file specified in the argument. The file needs to be created using the - cssens script. - """ + # Get input arguments usage = """ make_pull_at_sensitivity_limit filename [max_threads] diff --git a/examples/make_survey.py b/test/dev/make_survey.py index 5af2d9c..b923eb6 100755 --- a/examples/make_survey.py +++ b/test/dev/make_survey.py @@ -303,9 +303,7 @@ def survey_gplane(lrange=10, lstep=2): # Main routine entry point # #==========================# if __name__ == '__main__': - """ - CTA survey simulation and analysis script. - """ + # Initialise flags need_help = False diff --git a/examples/simulate_events.py b/test/dev/simulate_events.py index 8313e42..01fcac6 100755 --- a/examples/simulate_events.py +++ b/test/dev/simulate_events.py @@ -239,17 +239,15 @@ def show_events(events, xmlname, duration, emin, emax, ebins=30): # Main routine entry point # #==========================# if __name__ == '__main__': - """ - Simulate events. - """ + # Dump header - print("") - print("*******************") - print("* Simulate events *") - print("*******************") + print('') + print('*******************') + print('* Simulate events *') + print('*******************') # Set XML names - xmlname = "model_bgd.xml" + xmlname = 'model_bgd.xml' # Set simulation parameters duration = 360000.0 # Observation duration (s) diff --git a/test/test_cscaldb.py b/test/test_cscaldb.py index 83d8a43..f9dbed2 100755 --- a/test/test_cscaldb.py +++ b/test/test_cscaldb.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test cscaldb on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): cscaldb = '../cscripts/cscaldb.py' else: cscaldb = 'cscaldb' diff --git a/test/test_csfindobs.py b/test/test_csfindobs.py index c9da2c5..775494c 100755 --- a/test/test_csfindobs.py +++ b/test/test_csfindobs.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csfindobs on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): csfindobs = '../cscripts/csfindobs.py' else: csfindobs = 'csfindobs' diff --git a/test/test_csiactcopy.py b/test/test_csiactcopy.py index ec80654..c969b84 100755 --- a/test/test_csiactcopy.py +++ b/test/test_csiactcopy.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csiactcopy on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): csiactcopy = "../cscripts/csiactcopy.py" else: csiactcopy = "csiactcopy" diff --git a/test/test_csiactdata.py b/test/test_csiactdata.py index cb10874..14f6d2b 100755 --- a/test/test_csiactdata.py +++ b/test/test_csiactdata.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csiactdata on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): csiactdata = "../cscripts/csiactdata.py" else: csiactdata = "csiactdata" diff --git a/test/test_csiactobs.py b/test/test_csiactobs.py index 06dfd4a..122d659 100755 --- a/test/test_csiactobs.py +++ b/test/test_csiactobs.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csiactobs on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): csiactobs = "../cscripts/csiactobs.py" else: csiactobs = "csiactobs" diff --git a/test/test_cslightcrv.py b/test/test_cslightcrv.py index a8f8e9d..6f64efe 100755 --- a/test/test_cslightcrv.py +++ b/test/test_cslightcrv.py @@ -70,7 +70,7 @@ class Test(gammalib.GPythonTestSuite): Test cslightcrv on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): cslightcrv = '../cscripts/cslightcrv.py' else: cslightcrv = 'cslightcrv' diff --git a/test/test_csmodelinfo.py b/test/test_csmodelinfo.py index bb1ec1f..79818a6 100755 --- a/test/test_csmodelinfo.py +++ b/test/test_csmodelinfo.py @@ -69,7 +69,7 @@ class Test(gammalib.GPythonTestSuite): Test csmodelinfo on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): csmodelinfo = "../cscripts/csmodelinfo.py" else: csmodelinfo = "csmodelinfo" diff --git a/test/test_csmodelmerge.py b/test/test_csmodelmerge.py index c99b5a6..9d505f5 100755 --- a/test/test_csmodelmerge.py +++ b/test/test_csmodelmerge.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csmodelmerge on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): csmodelmerge = "../cscripts/csmodelmerge.py" else: csmodelmerge = "csmodelmerge" diff --git a/test/test_csobs2caldb.py b/test/test_csobs2caldb.py index 3a807e4..b294c12 100755 --- a/test/test_csobs2caldb.py +++ b/test/test_csobs2caldb.py @@ -73,7 +73,7 @@ class Test(gammalib.GPythonTestSuite): Test csobs2caldb on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): csobs2caldb = "../cscripts/csobs2caldb.py" else: csobs2caldb = "csobs2caldb" diff --git a/test/test_csobsdef.py b/test/test_csobsdef.py index 2d6f838..a2250c2 100755 --- a/test/test_csobsdef.py +++ b/test/test_csobsdef.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csobsdef on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): csobsdef = '../cscripts/csobsdef.py' else: csobsdef = 'csobsdef' diff --git a/test/test_csobsinfo.py b/test/test_csobsinfo.py index 2e151a9..ca3daf7 100755 --- a/test/test_csobsinfo.py +++ b/test/test_csobsinfo.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csobsinfo on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): csobsinfo = "../cscripts/csobsinfo.py" else: csobsinfo = "csobsinfo" diff --git a/test/test_cspull.py b/test/test_cspull.py index d1de25a..8cca24a 100755 --- a/test/test_cspull.py +++ b/test/test_cspull.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test cspull on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): cspull = "../cscripts/cspull.py" else: cspull = "cspull" diff --git a/test/test_csresmap.py b/test/test_csresmap.py index 432e06a..5cec03f 100755 --- a/test/test_csresmap.py +++ b/test/test_csresmap.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csresmap on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): csresmap = '../cscripts/csresmap.py' else: csresmap = 'csresmap' diff --git a/test/test_csroot2caldb.py b/test/test_csroot2caldb.py index 004aae0..b0d4289 100644 --- a/test/test_csroot2caldb.py +++ b/test/test_csroot2caldb.py @@ -69,7 +69,7 @@ class Test(gammalib.GPythonTestSuite): Test csroot2caldb on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): csroot2caldb = "../cscripts/csroot2caldb.py" else: csroot2caldb = "csroot2caldb" diff --git a/test/test_cssens.py b/test/test_cssens.py index 370a33c..8f4e000 100755 --- a/test/test_cssens.py +++ b/test/test_cssens.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test cssens on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): cssens = '../cscripts/cssens.py' else: cssens = 'cssens' diff --git a/test/test_csspec.py b/test/test_csspec.py index cc2a300..83a78fd 100755 --- a/test/test_csspec.py +++ b/test/test_csspec.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test csspec on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): csspec = '../cscripts/csspec.py' else: csspec = 'csspec' diff --git a/test/test_cstsdist.py b/test/test_cstsdist.py index 9162ba9..10ff345 100755 --- a/test/test_cstsdist.py +++ b/test/test_cstsdist.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test cstsdist on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): cstsdist = '../cscripts/cstsdist.py' else: cstsdist = 'cstsdist' diff --git a/test/test_cstsmapmerge.py b/test/test_cstsmapmerge.py index e05752c..8be2a0c 100755 --- a/test/test_cstsmapmerge.py +++ b/test/test_cstsmapmerge.py @@ -66,7 +66,7 @@ class Test(gammalib.GPythonTestSuite): Test cstsmapmerge on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): cstsmapmerge = '../cscripts/cstsmapmerge.py' else: cstsmapmerge = 'cstsmapmerge' diff --git a/test/test_csworkflow.py b/test/test_csworkflow.py index 1c78f8c..fda26ab 100755 --- a/test/test_csworkflow.py +++ b/test/test_csworkflow.py @@ -64,7 +64,7 @@ class Test(gammalib.GPythonTestSuite): def _test_cmd(self): # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): csworkflow = '../cscripts/csworkflow.py' else: csworkflow = 'csworkflow' diff --git a/test/test_ctbin.py b/test/test_ctbin.py index eab2c70..22bd625 100755 --- a/test/test_ctbin.py +++ b/test/test_ctbin.py @@ -69,7 +69,7 @@ class Test(gammalib.GPythonTestSuite): Test ctbin on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): ctbin = "../src/ctbin/ctbin" else: ctbin = "ctbin" diff --git a/test/test_ctbkgcube.py b/test/test_ctbkgcube.py index 86f8128..5ae5de0 100755 --- a/test/test_ctbkgcube.py +++ b/test/test_ctbkgcube.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctbkgcube on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctbkgcube = '../src/ctbkgcube/ctbkgcube' else: ctbkgcube = 'ctbkgcube' diff --git a/test/test_ctbutterfly.py b/test/test_ctbutterfly.py index d0ef5ae..6d018c4 100755 --- a/test/test_ctbutterfly.py +++ b/test/test_ctbutterfly.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctbutterfly on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctbutterfly = '../src/ctbutterfly/ctbutterfly' else: ctbutterfly = 'ctbutterfly' diff --git a/test/test_ctcubemask.py b/test/test_ctcubemask.py index e9a6e7e..6380f5a 100755 --- a/test/test_ctcubemask.py +++ b/test/test_ctcubemask.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctcubemask on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctcubemask = '../src/ctcubemask/ctcubemask' else: ctcubemask = 'ctcubemask' diff --git a/test/test_ctedispcube.py b/test/test_ctedispcube.py index 8b0dc3c..38d8ce5 100755 --- a/test/test_ctedispcube.py +++ b/test/test_ctedispcube.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctedispcube on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctedispcube = '../src/ctedispcube/ctedispcube' else: ctedispcube = 'ctedispcube' diff --git a/test/test_cterror.py b/test/test_cterror.py index 7693807..173efaa 100755 --- a/test/test_cterror.py +++ b/test/test_cterror.py @@ -69,7 +69,7 @@ class Test(gammalib.GPythonTestSuite): Test cterror on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): cterror = '../src/cterror/cterror' else: cterror = 'cterror' diff --git a/test/test_ctexpcube.py b/test/test_ctexpcube.py index d9c02e1..2d8db47 100755 --- a/test/test_ctexpcube.py +++ b/test/test_ctexpcube.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctexpcube on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctexpcube = '../src/ctexpcube/ctexpcube' else: ctexpcube = 'ctexpcube' diff --git a/test/test_ctlike.py b/test/test_ctlike.py index df82b05..dfdff89 100755 --- a/test/test_ctlike.py +++ b/test/test_ctlike.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctlike on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctlike = '../src/ctlike/ctlike' else: ctlike = 'ctlike' diff --git a/test/test_ctmapcube.py b/test/test_ctmapcube.py index 852c0a7..e83dbe9 100644 --- a/test/test_ctmapcube.py +++ b/test/test_ctmapcube.py @@ -69,7 +69,7 @@ class Test(gammalib.GPythonTestSuite): Test ctmapcube on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): ctmapcube = "../src/ctmapcube/ctmapcube" else: ctmapcube = "ctmapcube" diff --git a/test/test_ctmodel.py b/test/test_ctmodel.py index 3bacae3..33c6475 100755 --- a/test/test_ctmodel.py +++ b/test/test_ctmodel.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctmodel on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctmodel = '../src/ctmodel/ctmodel' else: ctmodel = 'ctmodel' diff --git a/test/test_ctobssim.py b/test/test_ctobssim.py index 4413eee..ec65cf6 100755 --- a/test/test_ctobssim.py +++ b/test/test_ctobssim.py @@ -72,7 +72,7 @@ class Test(gammalib.GPythonTestSuite): Test ctobssim on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile("README"): + if os.path.isfile("README.md"): ctobssim = "../src/ctobssim/ctobssim" else: ctobssim = "ctobssim" diff --git a/test/test_ctpsfcube.py b/test/test_ctpsfcube.py index 626f33c..9fb5e79 100755 --- a/test/test_ctpsfcube.py +++ b/test/test_ctpsfcube.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctpsfcube on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctpsfcube = '../src/ctpsfcube/ctpsfcube' else: ctpsfcube = 'ctpsfcube' diff --git a/test/test_ctselect.py b/test/test_ctselect.py index 1ac97bd..1d2c386 100755 --- a/test/test_ctselect.py +++ b/test/test_ctselect.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctselect on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctselect = '../src/ctselect/ctselect' else: ctselect = 'ctselect' diff --git a/test/test_ctskymap.py b/test/test_ctskymap.py index 330bf44..fe2bc50 100755 --- a/test/test_ctskymap.py +++ b/test/test_ctskymap.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctskymap on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctskymap = '../src/ctskymap/ctskymap' else: ctskymap = 'ctskymap' diff --git a/test/test_cttsmap.py b/test/test_cttsmap.py index ef215d7..d068da2 100755 --- a/test/test_cttsmap.py +++ b/test/test_cttsmap.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test cttsmap on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): cttsmap = '../src/cttsmap/cttsmap' else: cttsmap = 'cttsmap' diff --git a/test/test_ctulimit.py b/test/test_ctulimit.py index 9cd8d31..139d75f 100755 --- a/test/test_ctulimit.py +++ b/test/test_ctulimit.py @@ -62,7 +62,7 @@ class Test(gammalib.GPythonTestSuite): Test ctulimit on the command line. """ # Kluge to set the command (installed version has no README file) - if os.path.isfile('README'): + if os.path.isfile('README.md'): ctulimit = '../src/ctulimit/ctulimit' else: ctulimit = 'ctulimit' diff --git a/test/test_examples.py b/test/test_examples.py new file mode 100755 index 0000000..21032ee --- /dev/null +++ b/test/test_examples.py @@ -0,0 +1,425 @@ +#! /usr/bin/env python +# ========================================================================== +# Performs unit tests for the ctools example scripts +# +# Copyright (C) 2016 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 . +# +# ========================================================================== +import os +import sys +import gammalib +import ctools +import cscripts + + +# =============================================== # +# Test class for example executables unit testing # +# =============================================== # +class Test(gammalib.GPythonTestSuite): + """ + Test class for GammaLib example executables + """ + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + gammalib.GPythonTestSuite.__init__(self) + + # Return + return + + # Execute Python script + def _execute_python(self, name, args='', log=''): + """ + Execute Python script + + Parameters + ---------- + name : str + Python script name with .py extension + args : str, optional + String with arguments to be passed to script + log : str, optional + String with name to be used for logging + """ + # Setup command + cmd = '../examples/' + name + '.py' + if len(args) > 0: + cmd += ' ' + args + + # Set logname + if len(log) > 0: + logname = log + else: + logname = name + + # Execute Python script, make sure we catch any exception + try: + rc = os.system(cmd+' > example_'+logname+'.log 2>&1') + except: + pass + + # Check if execution was successful + self.test_assert(rc == 0, 'Check "'+name+'" script from command line') + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions + """ + # Set test name + self.name('examples') + + # Append tests + self.append(self.test_generate_prod3_irfs, 'Test generate_prod3_irfs.py') + self.append(self.test_make_pointings, 'Test make_pointings.py') + self.append(self.test_make_spectrum, 'Test make_spectrum.py') + self.append(self.test_make_ts_distributions, 'Test make_ts_distributions.py') + self.append(self.test_pipeline_binned_disk, 'Test pipeline_binned_disk.py') + self.append(self.test_pipeline_binned_mem, 'Test pipeline_binned_mem.py') + self.append(self.test_pipeline_stacked_disk, 'Test pipeline_stacked_disk.py') + self.append(self.test_pipeline_stacked_mem, 'Test pipeline_stacked_mem.py') + self.append(self.test_pipeline_unbinned_disk, 'Test pipeline_unbinned_disk.py') + self.append(self.test_pipeline_unbinned_mem, 'Test pipeline_unbinned_mem.py') + self.append(self.test_show_butterfly, 'Test show_butterfly.py') + self.append(self.test_show_model, 'Test show_model.py') + self.append(self.test_show_pha, 'Test show_pha.py') + self.append(self.test_show_pull_evolution, 'Test show_pull_evolution.py') + self.append(self.test_show_pull_histogram, 'Test show_pull_histogram.py') + self.append(self.test_show_response, 'Test test_show_response.py') + self.append(self.test_show_sensitivity, 'Test test_show_sensitivity.py') + self.append(self.test_show_spectrum, 'Test test_show_spectrum.py') + self.append(self.test_show_ts_distribution, 'Test show_ts_distribution.py') + + # Return + return + + # Test make_pointings + def test_make_pointings(self): + """ + Test make_pointings + """ + # Set script name + script = 'make_pointings' + + # Execute script + self._execute_python(script, args='-h', log='make_pointings_h') + self._execute_python(script, args='gps', log='make_pointings_gps') + self._execute_python(script, args='gps3', log='make_pointings_gps3') + self._execute_python(script, args='extgal', log='make_pointings_extgal') + self._execute_python(script, args='gc', log='make_pointings_gc') + self._execute_python(script, args='lmc', log='make_pointings_lmc') + + #TODO: Do any testing + + # Return + return + + # Test generate_prod3_irfs + def test_generate_prod3_irfs(self): + """ + Test generate_prod3_irfs + """ + # Execute script + self._execute_python('generate_prod3_irfs') + + #TODO: Do any testing + + # Return + return + + # Test make_spectrum + def test_make_spectrum(self): + """ + Test make_spectrum + """ + # Execute script + self._execute_python('make_spectrum', args='eps') + + #TODO: Do any testing + + # Return + return + + # Test make_ts_distributions + def test_make_ts_distributions(self): + """ + Test make_ts_distributions + """ + # Execute script + self._execute_python('make_ts_distributions', args='-n 2 -e 0 -d 1800') + + #TODO: Do any testing + + # Return + return + + # Test pipeline_binned_disk + def test_pipeline_binned_disk(self): + """ + Test pipeline_binned_disk + """ + # Execute script + self._execute_python('pipeline_binned_disk') + + #TODO: Do any testing + + # Return + return + + # Test pipeline_binned_mem + def test_pipeline_binned_mem(self): + """ + Test pipeline_binned_mem + """ + # Execute script + self._execute_python('pipeline_binned_mem') + + #TODO: Do any testing + + # Return + return + + # Test pipeline_stacked_disk + def test_pipeline_stacked_disk(self): + """ + Test pipeline_stacked_disk + """ + # Execute script + self._execute_python('pipeline_stacked_disk') + + #TODO: Do any testing + + # Return + return + + # Test pipeline_stacked_mem + def test_pipeline_stacked_mem(self): + """ + Test pipeline_stacked_mem + """ + # Execute script + self._execute_python('pipeline_stacked_mem') + + #TODO: Do any testing + + # Return + return + + # Test pipeline_unbinned_disk + def test_pipeline_unbinned_disk(self): + """ + Test pipeline_unbinned_disk + """ + # Execute script + self._execute_python('pipeline_unbinned_disk') + + #TODO: Do any testing + + # Return + return + + # Test pipeline_unbinned_mem + def test_pipeline_unbinned_mem(self): + """ + Test pipeline_unbinned_mem + """ + # Execute script + self._execute_python('pipeline_unbinned_mem') + + #TODO: Do any testing + + # Return + return + + # Test show_butterfly + def test_show_butterfly(self): + """ + Test show_butterfly + """ + # Execute script + self._execute_python('show_butterfly', + args='data/butterfly.txt ' + 'example_show_butterfly.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_model + def test_show_model(self): + """ + Test show_model + """ + # Execute script + self._execute_python('show_model', + args='data/crab.xml Crab example_show_model.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_pha + def test_show_pha(self): + """ + Test show_pha + """ + # Execute script + self._execute_python('show_pha', + args='data/pha.fits example_show_pha.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_pull_evolution + def test_show_pull_evolution(self): + """ + Test show_pull_evolution + """ + # Execute script + self._execute_python('show_pull_evolution', + args='data/pull.dat Pull_Crab_Prefactor ' + 'example_show_pull_evolution.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_pull_histogram + def test_show_pull_histogram(self): + """ + Test show_pull_histogram + """ + # Execute script + self._execute_python('show_pull_histogram', + args='data/pull.dat Pull_Crab_Prefactor 50 ' + 'example_show_pull_histogram.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_response + def test_show_response(self): + """ + Test show_response + """ + # Execute script + self._execute_python('show_response', args='example_show_response.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_sensitivity + def test_show_sensitivity(self): + """ + Test show_sensitivity + """ + # Execute script + self._execute_python('show_sensitivity', + args='data/sensitivity.dat ' + 'example_show_sensitivity.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_spectrum + def test_show_spectrum(self): + """ + Test show_spectrum + """ + # Execute script + self._execute_python('show_spectrum', + args='data/spectrum.fits ' + 'example_show_spectrum.eps') + + #TODO: Do any testing + + # Return + return + + # Test show_ts_distribution + def test_show_ts_distribution(self): + """ + Test show_ts_distribution + """ + # Execute script + self._execute_python('show_ts_distribution', + args='data/ts.dat -n 100 ' + '-f example_show_ts_distribution.eps') + + #TODO: Do any testing + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Set calibration database + os.environ['CALDB'] = '../caldb' + + # Set PFILES environment variable + try: + os.mkdir('pfiles') + except: + pass + os.environ['PFILES'] = 'pfiles' + + # Copy over pfiles + os.system('cp -r ../src/*/*.par pfiles/') + os.system('cp -r ../cscripts/*.par pfiles/') + + # Allocate test suites + suites = gammalib.GTestSuites('Examples testing') + + # Allocate test suite, setup tests and append them to the container + suite = Test() + suite.set() + suites.append(suite) + + # Run test suite + success = suites.run() + + # Save test results + suites.save('reports/examples.xml') + + # Set return code + if success: + rc = 0 + else: + rc = 1 + + # Exit with return code + sys.exit(rc) diff --git a/test/test_examples.sh.in b/test/test_examples.sh.in new file mode 100755 index 0000000..fc5d689 --- /dev/null +++ b/test/test_examples.sh.in @@ -0,0 +1,25 @@ +#! /bin/sh +# ========================================================================== +# This script launches the examples unit tests +# +# Copyright (C) 2016 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 . +# +# ========================================================================== + +# +# Run the examples test +# ===================== +@PYTHON_COVERAGE@./test_examples.py diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py index af6dbc3..bd8d68c 100755 --- a/test/test_python_cscripts.py +++ b/test/test_python_cscripts.py @@ -188,6 +188,10 @@ def test(installed=False, debug=False): else: rc = 1 + # Remove temporary direction + if installed: + os.system('rm -rf %s' % (path)) + # Exit with return code sys.exit(rc) diff --git a/test/test_python_ctools.py b/test/test_python_ctools.py index 67cf118..92d6291 100755 --- a/test/test_python_ctools.py +++ b/test/test_python_ctools.py @@ -158,6 +158,10 @@ def test(installed=False): else: rc = 1 + # Remove temporary direction + if installed: + os.system('rm -rf %s' % (path)) + # Exit with return code sys.exit(rc) diff --git a/testall.sh b/testall.sh index 91f1682..bb429a0 100755 --- a/testall.sh +++ b/testall.sh @@ -1,8 +1,8 @@ #!/bin/sh # ===================================================================== -# Run all test and example scripts that come with the package. +# Run some test scripts that come with the package # -# Copyright (C) 2011-2015 Juergen Knoedlseder +# Copyright (C) 2011-2016 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 @@ -25,11 +25,11 @@ base=$PWD # -# examples +# test/dev # ======== echo -echo "=====> examples" -cd examples +echo "=====> test/dev" +cd test/dev # Create local environment rm -rf ctools @@ -41,18 +41,8 @@ export CTOOLS=$PWD/ctools # Run checkers ./check_models.py -# Run pipelines -./pipeline_binned_disk.py -./pipeline_binned_mem.py -./pipeline_stacked_disk.py -./pipeline_stacked_mem.py -./pipeline_unbinned_disk.py -./pipeline_unbinned_mem.py - # Run makers ./make_survey.py -./make_ts_distributions.py -n 2 -e 0 -d 1800 -#./make_pull_at_sensitivity_limit.py # Remove local environment rm -rf ctools