Nathan Kelley-Hoskins / ctools

Commit fc54dc9f699260b11776f70a301301ccda88d2d1

Authored by Jürgen Knödlseder
1 parent ea7f03e9

Implement simple zenith angle selection (assuming 20 and 40 deg zenith angles)

Inline Side-by-side
Showing 1 changed file with 109 additions and 38 deletions
examples/make_pointings.py
  View file @fc54dc9
... ... @@ -101,7 +101,7 @@ def get_positions(xmin, xmax, ymin, ymax, step):
101 101 # ====================== #
102 102 def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5,
103 103 separation=3.0, hours=100.0,
104   - site=None, lst=True, autodec=0.0):
  104 + site=None, caldb='prod2', lst=True, autodec=0.0):
105 105 """
106 106 Setup pointing patch on the sky
107 107  
... ... @@ -170,26 +170,19 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5,
170 170 else:
171 171 obs_site = site
172 172  
173   - # Set site dependent Prod2 50 hour IRF and add up exposure for
174   - # a given site
175   - caldb = 'prod2'
176   - if obs_site == 'North':
177   - if lst:
178   - irf = 'North_50h'
179   - else:
180   - irf = 'North_50h'
181   - exposure_north += duration
182   - else:
183   - if lst:
184   - irf = 'South_50h'
185   - else:
186   - irf = 'South_50h'
187   - exposure_south += duration
  173 + # Set IRF
  174 + irf = set_irf(obs_site, obs, caldb, lst=lst)
188 175  
189 176 # Set IRF information
190 177 obs['caldb'] = caldb
191 178 obs['irf'] = irf
192 179  
  180 + # Add up exposure for a given site
  181 + if obs_site == 'North':
  182 + exposure_north += duration
  183 + else:
  184 + exposure_south += duration
  185 +
193 186 # Append observation
194 187 obsdef.append(obs)
195 188  
... ... @@ -206,10 +199,77 @@ def set_patch(lmin=-30.0, lmax=30.0, bmin=-1.5, bmax=+1.5,
206 199 return obsdef
207 200  
208 201  
  202 +# ================================== #
  203 +# Setup Instrument Response Function #
  204 +# ================================== #
  205 +def set_irf(site, obs, caldb, lst=True):
  206 + """
  207 + Setup Instrument Response Function
  208 +
  209 + Parameters
  210 + ----------
  211 + site : str
  212 + Array site (either 'South' or 'North')
  213 + obs : dict
  214 + Dictionary with pointing information
  215 + caldb : str
  216 + Calibration database
  217 + lst : bool, optional
  218 + Use LSTs
  219 +
  220 + Returns
  221 + -------
  222 + irf : str
  223 + IRF name
  224 + """
  225 + # Handle 'prod2'
  226 + if caldb == 'prod2':
  227 + if site == 'North':
  228 + irf = 'North_50h'
  229 + else:
  230 + irf = 'South_50h'
  231 +
  232 + # Handle 'prod3'
  233 + elif caldb == 'prod3':
  234 +
  235 + # Compute Right Ascension and Declination of pointing
  236 + pnt = gammalib.GSkyDir()
  237 + pnt.lb_deg(obs['lon'], obs['lat'])
  238 + ra = pnt.ra_deg()
  239 + dec = pnt.dec_deg()
  240 +
  241 + # Set geographic latitude of array
  242 + if site == 'North':
  243 + geolat = +28.7569
  244 + else:
  245 + geolat = -24.58
  246 +
  247 + # Compute best possible zenith angle
  248 + zenith = abs(dec - geolat)
  249 + print(dec, zenith)
  250 +
  251 + # Set IRF
  252 + if site == 'North':
  253 + irf = 'North_50h'
  254 + else:
  255 + irf = 'South_50h'
  256 + if zenith < 30.0:
  257 + irf += '_z20'
  258 + else:
  259 + irf += '_z40'
  260 +
  261 + # ... otherwise return an empty string
  262 + else:
  263 + irf = ''
  264 +
  265 + # Return IRF
  266 + return irf
  267 +
  268 +
209 269 # =========================== #
210 270 # Setup Galactic plane survey #
211 271 # =========================== #
212   -def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True):
  272 +def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, caldb='prod2', lst=True):
213 273 """
214 274 Setup Galactic plane survey
215 275  
... ... @@ -221,6 +281,8 @@ def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True):
221 281 Minimum Galactic Latitude (deg)
222 282 bmax : float, optional
223 283 Maximum Galactic Latitude (deg)
  284 + caldb : str, optional
  285 + Calibration database name
224 286 lst : bool, optional
225 287 Use LSTs
226 288  
... ... @@ -235,27 +297,27 @@ def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True):
235 297 # Add inner region South
236 298 obsdef.extend(set_patch(lmin=-60.0, lmax=60.0, bmin=bmin, bmax=bmax,
237 299 separation=separation, hours=780,
238   - site='South', lst=lst))
  300 + site='South', caldb=caldb, lst=lst))
239 301  
240 302 # Add Vela & Carina region
241 303 obsdef.extend(set_patch(lmin=240.0, lmax=300.0, bmin=bmin, bmax=bmax,
242 304 separation=separation, hours=180,
243   - site='South', lst=lst))
  305 + site='South', caldb=caldb, lst=lst))
244 306  
245 307 # Add 210-240 region
246 308 obsdef.extend(set_patch(lmin=210.0, lmax=240.0, bmin=bmin, bmax=bmax,
247 309 separation=separation, hours=60,
248   - site='South', lst=lst))
  310 + site='South', caldb=caldb, lst=lst))
249 311  
250 312 # Add Cygnus, Perseus
251 313 obsdef.extend(set_patch(lmin=60.0, lmax=150.0, bmin=bmin, bmax=bmax,
252 314 separation=separation, hours=450,
253   - site='North', lst=lst))
  315 + site='North', caldb=caldb, lst=lst))
254 316  
255 317 # Add Anticentre
256 318 obsdef.extend(set_patch(lmin=150.0, lmax=210.0, bmin=bmin, bmax=bmax,
257 319 separation=separation, hours=150,
258   - site='North', lst=lst))
  320 + site='North', caldb=caldb, lst=lst))
259 321  
260 322 # Return observation definition
261 323 return obsdef
... ... @@ -264,7 +326,7 @@ def set_gps(separation=3.0, bmin=-1.3, bmax=1.3, lst=True):
264 326 # ========================== #
265 327 # Setup Extragalactic survey #
266 328 # ========================== #
267   -def set_extgal(separation=3.0, lst=True):
  329 +def set_extgal(separation=3.0, caldb='prod2', lst=True):
268 330 """
269 331 Setup Extragalactic survey.
270 332  
... ... @@ -272,6 +334,8 @@ def set_extgal(separation=3.0, lst=True):
272 334 ----------
273 335 separation : float, optional
274 336 Pointing separation (deg)
  337 + caldb : str, optional
  338 + Calibration database name
275 339 lst : bool, optional
276 340 Use LSTs
277 341  
... ... @@ -286,7 +350,8 @@ def set_extgal(separation=3.0, lst=True):
286 350 # Set patch
287 351 obsdef.extend(set_patch(lmin=-90.0, lmax=90.0, bmin=+5.0, bmax=+88.0,
288 352 separation=separation, hours=1000,
289   - site='Automatic', lst=lst, autodec=10.0))
  353 + site='Automatic', caldb=caldb, lst=lst,
  354 + autodec=10.0))
290 355  
291 356 # Return observation definition
292 357 return obsdef
... ... @@ -295,12 +360,14 @@ def set_extgal(separation=3.0, lst=True):
295 360 # ========================= #
296 361 # Setup Galactic centre KSP #
297 362 # ========================= #
298   -def set_gc(lst=True):
  363 +def set_gc(caldb='prod2', lst=True):
299 364 """
300 365 Setup Galactic centre KSP.
301 366  
302 367 Parameters
303 368 ----------
  369 + caldb : str, optional
  370 + Calibration database name
304 371 lst : bool, optional
305 372 Use LSTs
306 373  
... ... @@ -315,12 +382,12 @@ def set_gc(lst=True):
315 382 # Central wobble
316 383 obsdef.extend(set_patch(lmin=-0.5, lmax=0.5, bmin=-0.5, bmax=0.5,
317 384 separation=0.1, hours=525,
318   - site='South', lst=lst))
  385 + site='South', caldb=caldb, lst=lst))
319 386  
320 387 # Extended region
321 388 obsdef.extend(set_patch(lmin=-10.0, lmax=10.0, bmin=-10.0, bmax=10.0,
322 389 separation=1.5, hours=300,
323   - site='South', lst=lst))
  390 + site='South', caldb=caldb, lst=lst))
324 391  
325 392 # Return observation definition
326 393 return obsdef
... ... @@ -329,7 +396,7 @@ def set_gc(lst=True):
329 396 # ============= #
330 397 # Setup LMC KSP #
331 398 # ============= #
332   -def set_lmc(hours=250.0, lst=True):
  399 +def set_lmc(hours=250.0, caldb='prod2', lst=True):
333 400 """
334 401 Setup LMC KSP.
335 402  
... ... @@ -337,6 +404,8 @@ def set_lmc(hours=250.0, lst=True):
337 404 ----------
338 405 hours : float, optional
339 406 Total observation duration (h)
  407 + caldb : str, optional
  408 + Calibration database name
340 409 lst : bool, optional
341 410 Use LSTs
342 411  
... ... @@ -372,12 +441,10 @@ def set_lmc(hours=250.0, lst=True):
372 441 # Set positions and duration
373 442 obs = {'lon': lon, 'lat': lat, 'duration': duration}
374 443  
  444 + # Set IRF
  445 + irf = set_irf(obs_site, obs, caldb, lst=lst)
  446 +
375 447 # Add IRF
376   - caldb = 'prod2'
377   - if lst:
378   - irf = 'South_50h'
379   - else:
380   - irf = 'South_50h'
381 448 obs['caldb'] = caldb
382 449 obs['irf'] = irf
383 450  
... ... @@ -467,27 +534,31 @@ def make_pointings():
467 534 print(' lmc LMC survey')
468 535 sys.exit()
469 536  
  537 + # Set calibration database
  538 + caldb = 'prod3'
  539 + lst = True
  540 +
470 541 # Galactic plane survey
471 542 if obsname == 'gps':
472   - obsdef = set_gps(lst=True)
  543 + obsdef = set_gps(caldb=caldb, lst=lst)
473 544 write_obsdef('gps.dat', obsdef)
474 545 elif obsname == 'gps3':
475   - obsdef = set_gps(separation=1.5, lst=True)
  546 + obsdef = set_gps(separation=1.5, caldb=caldb, lst=lst)
476 547 write_obsdef('gps3.dat', obsdef)
477 548  
478 549 # Extragalactic survey
479 550 elif obsname == 'extgal':
480   - obsdef = set_extgal(lst=True)
  551 + obsdef = set_extgal(caldb=caldb, lst=lst)
481 552 write_obsdef('extgal.dat', obsdef)
482 553  
483 554 # Galactic centre
484 555 elif obsname == 'gc':
485   - obsdef = set_gc(lst=True)
  556 + obsdef = set_gc(caldb=caldb, lst=lst)
486 557 write_obsdef('gc.dat', obsdef)
487 558  
488 559 # LMC
489 560 elif obsname == 'lmc':
490   - obsdef = set_lmc(lst=True)
  561 + obsdef = set_lmc(caldb=caldb, lst=lst)
491 562 write_obsdef('lmc.dat', obsdef)
492 563  
493 564 # Invalid pattern
... ...