Thierry Reposeur / ctools

Commit ce3d128823e212fb7e0d8dea28595b1399a63e4d

Authored by Michael Mayer
1 parent 7a6df46b

adapted cstsdist and test cases

Inline Side-by-side
Showing 2 changed files with 55 additions and 66 deletions
scripts/cstsdist.py
  View file @ce3d128
... ... @@ -20,16 +20,16 @@
20 20 #
21 21 # ==========================================================================
22 22 import gammalib
  23 +import ctools
23 24 from ctools import obsutils
24 25 import sys
25 26 import csv
26   -import math
27 27  
28 28  
29 29 # ============== #
30 30 # cstsdist class #
31 31 # ============== #
32   -class cstsdist(gammalib.GApplication):
  32 +class cstsdist(ctools.cscript):
33 33 """
34 34 This class implements the TS distribution generation script. It derives
35 35 from the GammaLib::GApplication class which provides support for parameter
... ... @@ -54,9 +54,9 @@ class cstsdist(gammalib.GApplication):
54 54  
55 55 # Initialise application
56 56 if len(argv) == 0:
57   - gammalib.GApplication.__init__(self, self.name, self.version)
  57 + ctools.cscript.__init__(self, self.name, self.version)
58 58 elif len(argv) ==1:
59   - gammalib.GApplication.__init__(self, self.name, self.version, *argv)
  59 + ctools.cscript.__init__(self, self.name, self.version, *argv)
60 60 else:
61 61 raise TypeError("Invalid number of arguments given.")
62 62  
... ... @@ -94,21 +94,28 @@ class cstsdist(gammalib.GApplication):
94 94  
95 95 # Create default parfile
96 96 pars = gammalib.GApplicationPars()
  97 + pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file"))
97 98 pars.append(gammalib.GApplicationPar("outfile","f","h","ts.dat","","","Output file name"))
98   - pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
99   - pars.append(gammalib.GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  99 + pars.append(gammalib.GApplicationPar("expcube","s","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
  100 + pars.append(gammalib.GApplicationPar("psfcube","s","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
  101 + pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
100 102 pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
101   - pars.append(gammalib.GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
102   - pars.append(gammalib.GApplicationPar("index","r","h","-2.48","","","Spectral index"))
103   - pars.append(gammalib.GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
104   - pars.append(gammalib.GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)"))
105   - pars.append(gammalib.GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
106   - pars.append(gammalib.GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  103 + pars.append(gammalib.GApplicationPar("ra","r","h","83.6331","0","360","RA of pointing (deg)"))
  104 + pars.append(gammalib.GApplicationPar("dec","r","h","22.0145","-90","90","Dec of pointing (deg)"))
  105 + pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit (TeV)"))
  106 + pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit (TeV)"))
  107 + pars.append(gammalib.GApplicationPar("tmin","r","h","0.0","","","Start time (MET in s)"))
  108 + pars.append(gammalib.GApplicationPar("tmax","r","a","1800.0","","","Observation duration (in s)"))
107 109 pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
108   - pars.append(gammalib.GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
109   - pars.append(gammalib.GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
110 110 pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
111 111 pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  112 + pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))
  113 + pars.append(gammalib.GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  114 + pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
  115 + pars.append(gammalib.GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
  116 + pars.append(gammalib.GApplicationPar("index","r","h","-2.48","","","Spectral index"))
  117 + pars.append(gammalib.GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
  118 + pars.append(gammalib.GApplicationPar("bkg","s","a","$CTOOLS/share/models/bkg_dummy.txt","","","Background model file function"))
112 119 pars.append_standard()
113 120 pars.save(parfile)
114 121  
... ... @@ -119,38 +126,45 @@ class cstsdist(gammalib.GApplication):
119 126 """
120 127 Get parameters from parfile and setup the observation.
121 128 """
122   - # Get parameters
  129 +
  130 + # Set observation if not done before
  131 + if self.obs == None or self.obs.size() == 0:
  132 + self.obs = self.get_observations()
  133 +
  134 + # Get number of energy bins
  135 + self.m_enumbins = self["enumbins"].integer()
  136 +
  137 + # Read parameters for binned if requested
  138 + if not self.m_enumbins == 0:
  139 + self.m_npix = self["npix"].integer()
  140 + self.m_binsz = self["binsz"].real()
  141 + else:
  142 + # Set dummy values (required by obsutils)
  143 + self.m_npix = 0
  144 + self.m_binsz = 0.0
  145 +
  146 + # Get other parameters
123 147 self.m_outfile = self["outfile"].filename()
124 148 self.m_ntrials = self["ntrials"].integer()
125   - self.m_caldb = self["caldb"].string()
126   - self.m_irf = self["irf"].string()
127 149 self.m_type = self["type"].string()
128 150 self.m_index = self["index"].real()
129 151 self.m_offset = self["offset"].real()
130 152 self.m_bkg = self["bkg"].string()
131   - self.m_emin = self["emin"].real()
132   - self.m_emax = self["emax"].real()
133   - self.m_enumbins = self["enumbins"].integer()
134   - self.m_duration = self["duration"].real()
135   - self.m_rad = self["rad"].real()
136   - self.m_npix = self["npix"].integer()
137   - self.m_binsz = self["binsz"].real()
138 153  
139 154 # Set some fixed parameters
140 155 self.m_log = False # Logging in client tools
141 156 self.m_debug = False # Debugging in client tools
142   -
143   - # Setup observation
144   - self.obs = self.set_obs(emin=self.m_emin, emax=self.m_emax)
145   -
146   - # Initialise models. Note that we centre the point source at the Galactic
147   - # center as our observation is also centred at the Galactic centre, so
148   - # we're onaxis.
  157 +
  158 + pnt = gammalib.GSkyDir()
  159 + pnt.radec_deg(self["ra"].real(),self["dec"].real())
  160 +
  161 + # Initialise models. Note that we centre the point source at the
  162 + # center of our observation, so we're onaxis.
149 163 self.bkg_model = gammalib.GModels()
150 164 self.full_model = gammalib.GModels()
151 165 self.bkg_model.append(self.set_bkg_model())
152 166 self.full_model.append(self.set_bkg_model())
153   - self.full_model.append(self.set_src_model(0.0, self.m_offset, \
  167 + self.full_model.append(self.set_src_model(pnt.l_deg(), pnt.b_deg()+self.m_offset, \
154 168 flux=0.010, \
155 169 type=self.m_type, \
156 170 index=self.m_index))
... ... @@ -166,7 +180,7 @@ class cstsdist(gammalib.GApplication):
166 180 Set model.
167 181 """
168 182 # Copy models
169   - self.model = models.copy()
  183 + self.obs.models(models.copy())
170 184  
171 185 # Return
172 186 return
... ... @@ -235,35 +249,7 @@ class cstsdist(gammalib.GApplication):
235 249  
236 250 # Return
237 251 return
238   -
239   - def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
240   - """
241   - Returns an observation container with a single CTA observation.
242   -
243   - Keywords:
244   - lpnt - Galactic longitude of pointing [deg] (default: 0.0)
245   - bpnt - Galactic latitude of pointing [deg] (default: 0.0)
246   - emin - Minimum energy [TeV] (default: 0.1)
247   - emax - Maximum energy [TeV] (default: 100.0)
248   - """
249   - # Allocate observation container
250   - obs = gammalib.GObservations()
251   -
252   - # Set single pointing
253   - pntdir = gammalib.GSkyDir()
254   - pntdir.lb_deg(lpnt, bpnt)
255   -
256   - # Create CTA observation
257   - run = obsutils.set_obs(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
258   - duration=self.m_duration, \
259   - emin=emin, emax=emax, rad=self.m_rad)
260   -
261   - # Append observation to container
262   - obs.append(run)
263   -
264   - # Return observation container
265   - return obs
266   -
  252 +
267 253 def set_bkg_model(self, fitsigma=False):
268 254 """
269 255 Setup CTA background model.
... ... @@ -343,10 +329,13 @@ class cstsdist(gammalib.GApplication):
343 329 self.log("ERROR: Unknown source type '"+type+"'.\n")
344 330 return None
345 331 source = gammalib.GModelSky(spatial, spectrum)
346   -
  332 +
347 333 # Set source name
348 334 source.name("Test")
349   - source.tscalc(True)
  335 +
  336 + # Turn TS value computation off
  337 + # since we compute it by hand later on
  338 + source.tscalc(False)
350 339  
351 340 # Return source
352 341 return source
... ...
test/test_cscripts.sh.in
  View file @ce3d128
... ... @@ -62,7 +62,7 @@ chmod +x @abs_top_srcdir@/scripts/cspull.py
62 62 emin=0.1 \
63 63 emax=100.0 \
64 64 enumbins=0 \
65   - duration=1800.0 \
  65 + tmax=1800.0 \
66 66 deadc=0.95 \
67 67 rad=5.0 \
68 68 npix=200 \
... ... @@ -118,7 +118,7 @@ chmod +x @abs_top_srcdir@/scripts/cstsdist.py
118 118 emin=0.1 \
119 119 emax=100.0 \
120 120 enumbins=0 \
121   - duration=1800.0 \
  121 + tmax=1800.0 \
122 122 rad=5.0 \
123 123 npix=200 \
124 124 binsz=0.05 > cstsdist.log
... ...