Simon Bonnefoy / ctools

Commit c60d955477d73c48049063074eac914068f878a7

Authored by Jürgen Knödlseder
1 parent 113f32e9

Add pointing pattern definition and observation list definition to obsutils. Enh… 

…ance cspull so that observation patterns can be simulated.
Inline Side-by-side
Showing 2 changed files with 90 additions and 50 deletions
scripts/cspull.py
  View file @c60d955
@@ -18,8 +18,7 @@ @@ -18,8 +18,7 @@
18 # along with this program. If not, see <http://www.gnu.org/licenses/>. 18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 # 19 #
20 # ========================================================================== 20 # ==========================================================================
21 -from ctools import *  
22 -from gammalib import * 21 +import gammalib
23 import obsutils 22 import obsutils
24 import sys 23 import sys
25 import csv 24 import csv
@@ -29,7 +28,7 @@ import math @@ -29,7 +28,7 @@ import math
29 # ============ # 28 # ============ #
30 # cspull class # 29 # cspull class #
31 # ============ # 30 # ============ #
32 -class cspull(GApplication): 31 +class cspull(gammalib.GApplication):
33 """ 32 """
34 This class implements the pull distribution generation script. It derives 33 This class implements the pull distribution generation script. It derives
35 from the GammaLib::GApplication class which provides support for parameter 34 from the GammaLib::GApplication class which provides support for parameter
@@ -42,7 +41,7 @@ class cspull(GApplication): @@ -42,7 +41,7 @@ class cspull(GApplication):
42 """ 41 """
43 # Set name 42 # Set name
44 self.name = "cspull" 43 self.name = "cspull"
45 - self.version = "0.2.0" 44 + self.version = "0.3.0"
46 45
47 # Initialise some members 46 # Initialise some members
48 self.obs = None 47 self.obs = None
@@ -54,9 +53,9 @@ class cspull(GApplication): @@ -54,9 +53,9 @@ class cspull(GApplication):
54 53
55 # Initialise application 54 # Initialise application
56 if len(argv) == 0: 55 if len(argv) == 0:
57 - GApplication.__init__(self, self.name, self.version) 56 + gammalib.GApplication.__init__(self, self.name, self.version)
58 elif len(argv) ==1: 57 elif len(argv) ==1:
59 - GApplication.__init__(self, self.name, self.version, *argv) 58 + gammalib.GApplication.__init__(self, self.name, self.version, *argv)
60 else: 59 else:
61 raise TypeError("Invalid number of arguments given.") 60 raise TypeError("Invalid number of arguments given.")
62 61
@@ -87,29 +86,31 @@ class cspull(GApplication): @@ -87,29 +86,31 @@ class cspull(GApplication):
87 parfile = self.name+".par" 86 parfile = self.name+".par"
88 87
89 try: 88 try:
90 - pars = GApplicationPars(parfile) 89 + pars = gammalib.GApplicationPars(parfile)
91 except: 90 except:
92 # Signal if parfile was not found 91 # Signal if parfile was not found
93 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") 92 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
94 93
95 # Create default parfile 94 # Create default parfile
96 - pars = GApplicationPars()  
97 - pars.append(GApplicationPar("srcmdl","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))  
98 - pars.append(GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))  
99 - pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))  
100 - pars.append(GApplicationPar("caldb","s","h","","","","Calibration database"))  
101 - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))  
102 - pars.append(GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))  
103 - pars.append(GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)"))  
104 - pars.append(GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))  
105 - pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))  
106 - pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))  
107 - pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))  
108 - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))  
109 - pars.append(GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))  
110 - pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))  
111 - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))  
112 - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) 95 + pars = gammalib.GApplicationPars()
  96 + pars.append(gammalib.GApplicationPar("srcmdl","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
  97 + pars.append(gammalib.GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))
  98 + pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials"))
  99 + pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
  100 + pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  101 + pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
  102 + pars.append(gammalib.GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)"))
  103 + pars.append(gammalib.GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))
  104 + pars.append(gammalib.GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
  105 + pars.append(gammalib.GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  106 + pars.append(gammalib.GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
  107 + pars.append(gammalib.GApplicationPar("duration","r","a","1800.0","","","Effective exposure time (s)"))
  108 + pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))
  109 + pars.append(gammalib.GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  110 + pars.append(gammalib.GApplicationPar("pattern","s","h","single","","","Observation pattern (single/four)"))
  111 + pars.append(gammalib.GApplicationPar("offset","r","h","1.5","","","Observation pattern offset (deg)"))
  112 + pars.append(gammalib.GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  113 + pars.append(gammalib.GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
113 pars.append_standard() 114 pars.append_standard()
114 pars.save(parfile) 115 pars.save(parfile)
115 116
@@ -136,6 +137,8 @@ class cspull(GApplication): @@ -136,6 +137,8 @@ class cspull(GApplication):
136 self.m_duration = self["duration"].real() 137 self.m_duration = self["duration"].real()
137 self.m_deadc = self["deadc"].real() 138 self.m_deadc = self["deadc"].real()
138 self.m_rad = self["rad"].real() 139 self.m_rad = self["rad"].real()
  140 + self.m_pattern = self["pattern"].string()
  141 + self.m_offset = self["offset"].real()
139 self.m_npix = self["npix"].integer() 142 self.m_npix = self["npix"].integer()
140 self.m_binsz = self["binsz"].real() 143 self.m_binsz = self["binsz"].real()
141 144
@@ -143,13 +146,12 @@ class cspull(GApplication): @@ -143,13 +146,12 @@ class cspull(GApplication):
143 self.m_log = False # Logging in client tools 146 self.m_log = False # Logging in client tools
144 self.m_debug = False # Debugging in client tools 147 self.m_debug = False # Debugging in client tools
145 148
146 - # Setup observation  
147 - self.obs = self.set_obs(ra=self.m_ra, dec=self.m_dec, \  
148 - emin=self.m_emin, emax=self.m_emax) 149 + # Setup observations
  150 + self.obs = self.set_obs()
149 151
150 # Load source model 152 # Load source model
151 if self.m_srcmdl != None: 153 if self.m_srcmdl != None:
152 - self.model = GModels(self.m_srcmdl) 154 + self.model = gammalib.GModels(self.m_srcmdl)
153 155
154 # Append source model to observation 156 # Append source model to observation
155 self.obs.models(self.model) 157 self.obs.models(self.model)
@@ -225,34 +227,28 @@ class cspull(GApplication): @@ -225,34 +227,28 @@ class cspull(GApplication):
225 # Return 227 # Return
226 return 228 return
227 229
228 - def set_obs(self, ra=83.6331, dec=22.0145, emin=0.1, emax=100.0): 230 + def set_obs(self):
229 """ 231 """
230 - Returns an observation container with a single CTA observation. 232 + Returns an observation container with a set of CTA observations.
231 233
232 Keywords: 234 Keywords:
233 - ra - Right Ascension of pointing [deg] (default: 83.6331)  
234 - dec - Declination of pointing [deg] (default: 22.0145)  
235 - emin - Minimum energy [TeV] (default: 0.1)  
236 - emax - Maximum energy [TeV] (default: 100.0)  
237 """ 235 """
238 - # Allocate observation container  
239 - obs = GObservations()  
240 -  
241 - # Set single pointing  
242 - pntdir = GSkyDir()  
243 - pntdir.radec_deg(ra, dec)  
244 -  
245 - # Create CTA observation  
246 - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \  
247 - duration=self.m_duration, deadc=self.m_deadc, \  
248 - emin=emin, emax=emax, rad=self.m_rad) 236 + # Setup observation definition list
  237 + obsdeflist = obsutils.set_obs_patterns(self.m_pattern, \
  238 + ra=self.m_ra, dec=self.m_dec, \
  239 + offset=self.m_offset)
249 240
250 - # Append observation to container  
251 - obs.append(run) 241 + # Create list of observations
  242 + obs = obsutils.set_obs_list(obsdeflist, \
  243 + tstart=0.0, duration=self.m_duration, \
  244 + deadc=self.m_deadc, \
  245 + emin=self.m_emin, emax=self.m_emax, \
  246 + rad=self.m_rad, \
  247 + irf=self.m_irf, caldb=self.m_caldb)
252 248
253 # Return observation container 249 # Return observation container
254 return obs 250 return obs
255 - 251 +
256 def trial(self, seed): 252 def trial(self, seed):
257 """ 253 """
258 Create the pull for a single trial. 254 Create the pull for a single trial.
scripts/obsutils.py
  View file @c60d955
@@ -342,8 +342,8 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ @@ -342,8 +342,8 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
342 gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) 342 gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration))
343 343
344 # Set energy boundaries 344 # Set energy boundaries
345 - ebounds = gammalib.GEbounds()  
346 - ebounds.append(gammalib.GEnergy(emin, "TeV"), gammalib.GEnergy(emax, "TeV")) 345 + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \
  346 + gammalib.GEnergy(emax, "TeV"))
347 347
348 # Allocate event list 348 # Allocate event list
349 events = gammalib.GCTAEventList() 349 events = gammalib.GCTAEventList()
@@ -465,9 +465,12 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ @@ -465,9 +465,12 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
465 obs_cta = set_obs(pntdir, tstart=obs_start, duration=obs_duration, \ 465 obs_cta = set_obs(pntdir, tstart=obs_start, duration=obs_duration, \
466 deadc=obs_deadc, emin=obs_emin, emax=obs_emax, \ 466 deadc=obs_deadc, emin=obs_emin, emax=obs_emax, \
467 rad=obs_rad, irf=obs_irf, caldb=obs_caldb, id=id) 467 rad=obs_rad, irf=obs_irf, caldb=obs_caldb, id=id)
  468 + print obs_cta.__repr__()
468 469
469 # Append to container 470 # Append to container
470 obs.append(obs_cta) 471 obs.append(obs_cta)
  472 + for r in obs:
  473 + print r.__repr__()
471 474
472 # Update start time and identifier 475 # Update start time and identifier
473 obs_start += obs_duration 476 obs_start += obs_duration
@@ -475,3 +478,44 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ @@ -475,3 +478,44 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
475 478
476 # Return observation container 479 # Return observation container
477 return obs 480 return obs
  481 +
  482 +
  483 +# ============================ #
  484 +# Set CTA observation patterns #
  485 +# ============================ #
  486 +def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
  487 + """
  488 + Sets a number of standard patterns.
  489 +
  490 + Parameters:
  491 + pattern - Observation pattern. Possible options are:
  492 + - "single": single pointing
  493 + - "four": four pointings 'offset' around pattern centre
  494 +
  495 + Keywords:
  496 + ra - Right Ascension of pattern centre [deg] (default: 83.6331)
  497 + dec - Declination of pattern centre [deg] (default: 22.0145)
  498 + offset - Offset from pattern centre [deg] (default: 1.5)
  499 + """
  500 + # Initialise observation definition list
  501 + obsdeflist = []
  502 +
  503 + # Add patterns
  504 + if pattern == "single":
  505 + obsdef = {'ra': ra, 'dec': dec}
  506 + obsdeflist.append(obsdef)
  507 + elif pattern == "four":
  508 + # Set pattern centre
  509 + centre = gammalib.GSkyDir()
  510 + centre.radec_deg(ra, dec)
  511 +
  512 + # Append pointings
  513 + for phi in [0.0, 90.0, 180.0, 270.0]:
  514 + pntdir = centre.copy()
  515 + pntdir.rotate_deg(phi, offset)
  516 + obsdeflist.append({'ra': pntdir.ra_deg(), 'dec': pntdir.dec_deg()})
  517 + else:
  518 + print("Warning: Observation pattern '"+str(pattern)+"' not recognized.")
  519 +
  520 + # Return observation definition list
  521 + return obsdeflist