Iftach Sadeh / ctools

Commit af530c3337fce06cff37f9d26a28523a2a121a5d

Authored by knodlseder
1 parent 8bcf40c9

small improvements

Inline Side-by-side
Showing 1 changed file with 65 additions and 775 deletions
scripts/obsutils.py
  View file @af530c3
1 1 # ==========================================================================
2 2 # This script provides a number of functions that are useful for handling
3   -# CTA data.
  3 +# CTA observations.
4 4 #
5 5 # Copyright (C) 2011 Jurgen Knodlseder
6 6 #
... ... @@ -20,7 +20,7 @@
20 20 # ==========================================================================
21 21 from ctatools import *
22 22 from gammalib import *
23   -from math import *
  23 +#from math import *
24 24 #import os
25 25 #import glob
26 26 #import sys
... ... @@ -32,7 +32,7 @@ from math import *
32 32 # ===================== #
33 33 # Simulate observations #
34 34 # ===================== #
35   -def sim_obs(obs, log=False, seed=0, nbins=0, binsz=0.05, npix=200):
  35 +def sim(obs, log=False, debug=True, seed=0, nbins=0, binsz=0.05, npix=200):
36 36 """
37 37 Simulate events for all observations in the container.
38 38  
... ... @@ -40,6 +40,7 @@ def sim_obs(obs, log=False, seed=0, nbins=0, binsz=0.05, npix=200):
40 40 obs - Observation container
41 41 Keywords:
42 42 log - Create log file(s)
  43 + debug - Create screen dump
43 44 seed - Seed value for simulations (default: 0)
44 45 nbins - Number of energy bins (default: 0=unbinned)
45 46 binsz - Pixel size for binned simulation (deg/pixel)
... ... @@ -52,6 +53,10 @@ def sim_obs(obs, log=False, seed=0, nbins=0, binsz=0.05, npix=200):
52 53 # Optionally open the log file
53 54 if log:
54 55 sim.logFileOpen()
  56 +
  57 + # Optionally switch-on debugging model
  58 + if debug:
  59 + sim["debug"].boolean(True)
55 60  
56 61 # Run ctobssim application. This will loop over all observations in the
57 62 # container and simulation the events for each observation. Note that
... ... @@ -92,6 +97,10 @@ def sim_obs(obs, log=False, seed=0, nbins=0, binsz=0.05, npix=200):
92 97 if log:
93 98 bin.logFileOpen()
94 99  
  100 + # Optionally switch-on debugging model
  101 + if debug:
  102 + bin["debug"].boolean(True)
  103 +
95 104 # Run ctbin application. This will loop over all observations in
96 105 # the container and bin the events in counts maps
97 106 bin.run()
... ... @@ -115,14 +124,15 @@ def sim_obs(obs, log=False, seed=0, nbins=0, binsz=0.05, npix=200):
115 124 # ================ #
116 125 # Fit observations #
117 126 # ================ #
118   -def fit_obs(obs, log=False):
  127 +def fit(obs, log=False, debug=False):
119 128 """
120 129 Perform maximum likelihood fitting of observations in the container.
121 130  
122 131 Parameters:
123   - obs - Observation container
  132 + obs - Observation container
124 133 Keywords:
125   - log - Create log file(s)
  134 + log - Create log file(s)
  135 + debug - Create screen dump
126 136 """
127 137 # Allocate ctlike application
128 138 like = ctlike(obs)
... ... @@ -130,6 +140,9 @@ def fit_obs(obs, log=False):
130 140 # Optionally open the log file
131 141 if log:
132 142 like.logFileOpen()
  143 +
  144 + # Optionally switch-on debugging model
  145 + if debug:
133 146 like["debug"].boolean(True)
134 147  
135 148 # Run ctlike application.
... ... @@ -142,25 +155,27 @@ def fit_obs(obs, log=False):
142 155 # ================= #
143 156 # Create counts map #
144 157 # ================= #
145   -def make_counts_map(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
146   - nxpix=600, nypix=200, pixsize=0.05, \
147   - outname="survey_cntmap.fits"):
  158 +def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
  159 + binsz=0.05, nxpix=200, nypix=200, \
  160 + outname="cntmap.fits"):
148 161 """
149 162 Creates a counts map by combining the events of all observations.
150 163 The counts map will be a summed map over all energies.
151 164  
152 165 Parameters:
153   - obs - Observation container
  166 + obs - Observation container
154 167 Keywords:
155   - xval - Reference longitude value [deg]
156   - yval - Reference latitude value [deg]
157   - nxpix - Number of pixels in longitude axis
158   - nypix - Number of pixels in latitude axis
159   - pixsize - Pixel size [deg]
160   - outname - Counts map FITS filename
  168 + proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN)
  169 + coord - Coordinate type (GAL, CEL) (default: GAL)
  170 + xval - Reference longitude value [deg] (default: 0.0)
  171 + yval - Reference latitude value [deg] (default: 0.0)
  172 + binsz - Pixel size [deg/pixel] (default: 0.05)
  173 + nxpix - Number of pixels in X direction (default: 200)
  174 + nypix - Number of pixels in Y direction (default: 200)
  175 + outname - Counts map FITS filename (default: cntmap.fits)
161 176 """
162 177 # Allocate counts map
163   - cntmap = GSkymap(proj, coord, xval, yval, -pixsize, pixsize, nxpix, nypix, 1)
  178 + map = GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
164 179  
165 180 # Fill all observations
166 181 for run in obs:
... ... @@ -173,43 +188,45 @@ def make_counts_map(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
173 188  
174 189 # Determine sky pixel
175 190 skydir = event.dir().skydir()
176   - pixel = cntmap.dir2pix(skydir)
  191 + pixel = map.dir2pix(skydir)
177 192  
178 193 # Set pixel
179   - cntmap[pixel] += 1.0
  194 + map[pixel] += 1.0
180 195  
181 196 # Save sky map. The clobber flag is set to True, so any existing FITS
182 197 # file will be overwritten.
183   - cntmap.save(outname, True)
  198 + map.save(outname, True)
184 199  
185 200 # Return counts map
186   - return cntmap
  201 + return map
187 202  
188 203  
189 204 # ================ #
190 205 # Create model map #
191 206 # ================ #
192   -def make_model_map(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
193   - nxpix=600, nypix=200, pixsize=0.05, \
194   - outname="survey_modmap.fits"):
  207 +def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
  208 + binsz=0.05, nxpix=200, nypix=200, \
  209 + outname="modmap.fits"):
195 210 """
196 211 Make model map for a given reference energy by combining all observations.
197 212 The model map will be evaluated for a given reference energy 'eref' and will
198 213 be given in units of [counts/(sr MeV s)].
199 214  
200 215 Parameters:
201   - obs - Observation container
  216 + obs - Observation container
202 217 Keywords:
203   - eref - Reference energy for which model is created [TeV]
204   - xval - Reference longitude value [deg]
205   - yval - Reference latitude value [deg]
206   - nxpix - Number of pixels in longitude axis
207   - nypix - Number of pixels in latitude axis
208   - pixsize - Pixel size [deg]
209   - outname - Counts map FITS filename
  218 + eref - Reference energy for which model is created [TeV] (default: 0.1)
  219 + proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN)
  220 + coord - Coordinate type (GAL, CEL) (default: GAL)
  221 + xval - Reference longitude value [deg] (default: 0.0)
  222 + yval - Reference latitude value [deg] (default: 0.0)
  223 + binsz - Pixel size [deg/pixel] (default: 0.05)
  224 + nxpix - Number of pixels in X direction (default: 200)
  225 + nypix - Number of pixels in Y direction (default: 200)
  226 + outname - Model map FITS filename (default: modmap.fits)
210 227 """
211 228 # Allocate model map
212   - modmap = GSkymap(proj, coord, xval, yval, -pixsize, pixsize, nxpix, nypix, 1)
  229 + map = GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
213 230  
214 231 # Set reference energy, time and direction. The time is not initialised and is
215 232 # in fact not used (as the IRF is assumed to be time independent for now).
... ... @@ -220,10 +237,10 @@ def make_model_map(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
220 237 energy.TeV(eref)
221 238  
222 239 # Loop over all map pixels
223   - for pixel in range(modmap.npix()):
  240 + for pixel in range(map.npix()):
224 241  
225 242 # Get sky direction
226   - skydir = modmap.pix2dir(pixel)
  243 + skydir = map.pix2dir(pixel)
227 244 instdir.skydir(skydir)
228 245  
229 246 # Create event atom for map pixel
... ... @@ -240,21 +257,21 @@ def make_model_map(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \
240 257 value += obs.models().eval(atom, run)
241 258  
242 259 # Set map value
243   - modmap[pixel] = value
  260 + map[pixel] = value
244 261  
245 262 # Save sky map
246   - modmap.save(outname, True)
  263 + map.save(outname, True)
247 264  
248 265 # Return model map
249   - return modmap
  266 + return map
250 267  
251 268  
252 269 # ======================= #
253 270 # Set one CTA observation #
254 271 # ======================= #
255   -def set_one_obs(pntdir, tstart=0.0, duration=1800.0, \
256   - emin=0.1, emax=100.0, rad=5.0, \
257   - irf="kb_E_50h_v3", caldb="/usr/local/gamma/share/caldb/cta"):
  272 +def set(pntdir, tstart=0.0, duration=1800.0, \
  273 + emin=0.1, emax=100.0, rad=5.0, \
  274 + irf="kb_E_50h_v3", caldb="$GAMMALIB/share/caldb/cta"):
258 275 """
259 276 Returns a single CTA observation. By looping over this function we can
260 277 add CTA observations to the observation container.
... ... @@ -262,13 +279,13 @@ def set_one_obs(pntdir, tstart=0.0, duration=1800.0, \
262 279 Parameters:
263 280 pntdir - Pointing direction
264 281 Keywords:
265   - tstart - Start time [seconds]
266   - duration - Duration of observation [seconds]
267   - emin - Minimum event energy [TeV]
268   - emax - Maximum event energy [TeV]
269   - rad - ROI radius used for analysis [deg]
270   - irf - Instrument response function
271   - caldb - Calibration database path
  282 + tstart - Start time [seconds] (default: 0.0)
  283 + duration - Duration of observation [seconds] (default: 1800.0)
  284 + emin - Minimum event energy [TeV] (default: 0.1)
  285 + emax - Maximum event energy [TeV] (default: 100.0)
  286 + rad - ROI radius used for analysis [deg] (default: 5.0)
  287 + irf - Instrument response function (default: kb_E_50h_v3)
  288 + caldb - Calibration database path (default: $GAMMALIB/share/caldb/cta)
272 289 """
273 290 # Allocate CTA observation
274 291 obs = GCTAObservation()
... ... @@ -310,733 +327,6 @@ def set_one_obs(pntdir, tstart=0.0, duration=1800.0, \
310 327  
311 328 # Set instrument response
312 329 obs.response(irf, caldb)
313   - #print obs
314 330  
315 331 # Return observation
316 332 return obs
317   -
318   -
319   -# ================== #
320   -# Setup source model #
321   -# ================== #
322   -def add_source_model(models, source):
323   - """
324   - This function adds one source model to the observations container.
325   -
326   - Parameters:
327   - obs - Observations container
328   - source - Source to be added to container
329   - """
330   - # Append source to model container
331   - models.append(source)
332   -
333   - # Return model container
334   - return models
335   -
336   -
337   -# ======================================== #
338   -# Add CTA background model to observations #
339   -# ======================================== #
340   -def add_background_model(models):
341   - """
342   - Add standard CTA background model to observations container.
343   - We use a simple power law here, scaled to Konrad's E configuration
344   - performance table.
345   -
346   - Parameters:
347   - obs - Observation container
348   - """
349   - # Define background model
350   - bgd_radial = GCTAModelRadialGauss(3.0)
351   - bgd_spectrum = GModelSpectralPlaw(61.8, -1.85)
352   - bgd_spectrum["Prefactor"].scale(1.0e-6)
353   - bgd_spectrum["PivotEnergy"].value(1.0)
354   - bgd_spectrum["PivotEnergy"].scale(1.0e6)
355   - bgd_model = GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum)
356   - bgd_model.name("Background")
357   - bgd_model.instruments("CTA")
358   -
359   - # Add background model to container
360   - models.append(bgd_model)
361   -
362   - # Return model container
363   - return models
364   -
365   -
366   -# ==================== #
367   -# Set Crab-like source #
368   -# ==================== #
369   -def crab_like(name, l, b, flux=1.0, type="point", sigma=1.0, radius=1.0, width=0.1, \
370   - fitpos=True):
371   - """
372   - Returns a single source with Crab-like spectrum. The source flux
373   - can be scaled in Crab units.
374   -
375   - Parameters:
376   - name - Unique source name
377   - l - Galactic longitude of source location [deg]
378   - b - Galactic latitude of source location [deg]
379   - Keywords:
380   - flux - Source flux in Crab units
381   - type - Source type ("point", "gauss", "disk", "shell")
382   - sigma - Gaussian sigma (for type="gauss")
383   - radius - Disk or shell inner radius [deg] (for type="disk" and type="shell")
384   - width - Shell width [deg] (for type="shell")
385   - fitpos - Specifies if position should be fitted or not
386   - """
387   - # Set source location
388   - location = GSkyDir()
389   - location.lb_deg(l, b)
390   -
391   - # Set source spectrum
392   - spectrum = crab_spec(flux=flux)
393   -
394   - # Set source
395   - if type == "point":
396   - spatial = GModelSpatialPtsrc(location)
397   - if fitpos:
398   - spatial[0].free()
399   - spatial[1].free()
400   - source = GModelPointSource(spatial, spectrum)
401   - elif type == "gauss":
402   - radial = GModelRadialGauss(location, sigma)
403   - if fitpos:
404   - radial[0].free()
405   - radial[1].free()
406   - source = GModelExtendedSource(radial, spectrum)
407   - elif type == "disk":
408   - radial = GModelRadialDisk(location, radius)
409   - if fitpos:
410   - radial[0].free()
411   - radial[1].free()
412   - source = GModelExtendedSource(radial, spectrum)
413   - elif type == "shell":
414   - radial = GModelRadialShell(location, radius, width)
415   - if fitpos:
416   - radial[0].free()
417   - radial[1].free()
418   - source = GModelExtendedSource(radial, spectrum)
419   - else:
420   - print "ERROR: Unknown source type '"+type+"'."
421   - return None
422   -
423   - # Set source name
424   - source.name(name)
425   -
426   - # Return source
427   - return source
428   -
429   -
430   -# ================= #
431   -# Set Crab spectrum #
432   -# ================= #
433   -def crab_spec(flux=1.0):
434   - """
435   - Set Crab spectrum based on MAGIC observations
436   - (Albert et al. 2008, ApJ, 674, 1037)
437   -
438   - Keywords:
439   - flux - Source flux in Crab units
440   - """
441   - # Set parameters
442   - spectrum = GModelSpectralPlaw(flux, -2.48)
443   - spectrum["Prefactor"].scale(5.7e-16)
444   - spectrum["PivotEnergy"].value(0.3)
445   - spectrum["PivotEnergy"].scale(1.0e6)
446   -
447   - # Return spectrum
448   - return spectrum
449   -
450   -
451   -# ========================= #
452   -# Set any point-like source #
453   -# ========================= #
454   -def any_pointlike(name, l, b, flux=1.0, index=-2.48, fitpos=True, fitidx=True):
455   - """
456   - Returns a single point-source with power-law spectrum.
457   - The source flux can be scaled in Crab units.
458   -
459   - Parameters:
460   - name - Unique source name
461   - l - Galactic longitude of source location [deg]
462   - b - Galactic latitude of source location [deg]
463   - Keywords:
464   - flux - Source flux in Crab units
465   - fitpos - Specifies if position should be fitted or not
466   - fitidx - Specifies if index should be fitted or not
467   - """
468   - # Set source location
469   - location = GSkyDir()
470   - location.lb_deg(l, b)
471   -
472   -
473   - # Set source spatial model
474   - spatial = GModelSpatialPtsrc(location)
475   -
476   - # Set source spectral model
477   - spectrum = GModelSpectralPlaw(flux, index)
478   - spectrum["Prefactor"].scale(5.7e-16)
479   - spectrum["PivotEnergy"].value(0.3)
480   - spectrum["PivotEnergy"].scale(1.0e6)
481   -
482   - # Set fit properties
483   - if fitpos:
484   - spatial[0].free()
485   - spatial[1].free()
486   - else:
487   - spatial[0].fix()
488   - spatial[1].fix()
489   -
490   - if fitidx:
491   - spectrum["Index"].free()
492   - else:
493   - spectrum["Index"].fix()
494   -
495   - # Assign spatial and spectral models to source object
496   - source = GModelPointSource(spatial, spectrum)
497   -
498   - # Set source name
499   - source.name(name)
500   -
501   - # Return source
502   - return source
503   -
504   -
505   -# =============================== #
506   -# Setup single observation survey #
507   -# =============================== #
508   -def survey_single(caldb, irf, lpnt=0.0, bpnt=0.0):
509   - """
510   - This is an example of how to create a pointing strategy. Here, we
511   - only have a single CTA pointing, so this is a very basic pointing
512   - strategy.
513   -
514   - Parameters:
515   - caldb - Calibration database
516   - irf - Response function
517   - lpnt - Galactic longitude of pointing [deg]
518   - bpnt - Galactic latitude of pointing [deg]
519   - """
520   - # Allocate observation container
521   - obs = GObservations()
522   -
523   - # Set single pointing at galactic centre
524   - pntdir = GSkyDir()
525   - pntdir.lb_deg(lpnt, bpnt)
526   - run = set_one_obs(pntdir, caldb=caldb, irf=irf)
527   - obs.append(run)
528   -
529   - # Return observation container
530   - return obs
531   -
532   -
533   -# =========================== #
534   -# Setup Galactic plane survey #
535   -# =========================== #
536   -def survey_gplane(caldb, irf, bpnt=0.0, lrange=10, lstep=2):
537   - """
538   - Creates a single observation survey for test purposes.
539   -
540   - Keywords:
541   - bpnt - Latitude value for scan [deg]
542   - lrange - Longitude range (integer deg)
543   - lstep - Longitude step size (integer deg)
544   - """
545   - # Allocate observation container
546   - obs = GObservations()
547   -
548   - # Loop over longitudes
549   - for l in range(-lrange,lrange+lstep,lstep):
550   -
551   - # Set pointing
552   - pntdir = GSkyDir()
553   - pntdir.lb_deg(l, bpnt)
554   - run = set_one_obs(pntdir, caldb=caldb, irf=irf)
555   - obs.append(run)
556   -
557   - # Return observation container
558   - return obs
559   -
560   -
561   -# =========================== #
562   -# Setup Galactic plane survey #
563   -# =========================== #
564   -def survey_square(caldb, irf, lrange=2.0, lstep=2.0, duration=1800.0, \
565   - emin=0.1, emax=100.0):
566   - """
567   - Creates a single observation survey over a square grid.
568   - For the moment, the grid is centred on (l,b)=(0,0)
569   -
570   - Arguments:
571   - lrange - Longitude/latitude range (float deg)
572   - lstep - Longitude/latitude step size (float deg)
573   - duration/emin/emax - See set_one_obs function
574   - """
575   - # Allocate observation container
576   - obs = GObservations()
577   -
578   - # Compute number of points along each axis
579   - num_pnt= int(round(2.0*lrange/lstep)+1.0)
580   -
581   - # Loop over longitudes and latitudes
582   - for i in range(0,num_pnt):
583   - for j in range(0,num_pnt):
584   -
585   - # Set pointing
586   - l= -lrange+i*lstep
587   - b= -lrange+j*lstep
588   - pntdir = GSkyDir()
589   - pntdir.lb_deg(l, b)
590   - run = set_one_obs(pntdir, caldb=caldb, irf=irf, duration=duration, \
591   - emin=emin, emax=emax)
592   - obs.append(run)
593   -
594   - # Return observation container
595   - return obs
596   -
597   -
598   -# =========================== #
599   -# Setup Galactic plane survey #
600   -# =========================== #
601   -def survey_galplane(caldb, irf, lmin=-2.0, lstep=2.0, nstep=3, duration=1800.0, \
602   - emin=0.1, emax=100.0):
603   - """
604   - Creates a single observation survey along the galactic plane
605   -
606   - Arguments:
607   - lmin - Starting longitude(float deg)
608   - lstep - Longitude step size (float deg)
609   - nstep - Number of steps/pointings
610   - duration/emin/emax - See set_one_obs function
611   - """
612   - # Allocate observation container
613   - obs = GObservations()
614   -
615   - # Loop over longitudes
616   - for i in range(0,nstep):
617   -
618   - # Set pointing
619   - l= lmin+i*lstep
620   - b= 0.0
621   - pntdir = GSkyDir()
622   - pntdir.lb_deg(l, b)
623   - run = set_one_obs(pntdir, caldb=caldb, irf=irf, duration=duration, \
624   - emin=emin, emax=emax)
625   - obs.append(run)
626   -
627   - # Return observation container
628   - return obs
629   -
630   -# ===================================== #
631   -# Setup set of point source Crab models #
632   -# ===================================== #
633   -def add_point_sources(obs):
634   - """
635   - Add point-source model made of Crab-like sources.
636   - """
637   - # Add point-source Crabs
638   - obs = add_source_model(obs, crab_like("Crab", 0.0, 0.0))
639   - obs = add_source_model(obs, crab_like("Crab 1deg off", 1.0, 0.0))
640   - obs = add_source_model(obs, crab_like("Crab 2deg off", 2.0, 0.0))
641   - obs = add_source_model(obs, crab_like("Crab 3deg off", 3.0, 0.0))
642   - obs = add_source_model(obs, crab_like("Crab 4deg off", 4.0, 0.0))
643   - obs = add_source_model(obs, crab_like("Crab 5deg off", 5.0, 0.0))
644   -
645   - # Return
646   - return obs
647   -
648   -
649   -# ================= #
650   -# Setup disk models #
651   -# ================= #
652   -def add_disk_sources(obs):
653   - """
654   - Add disk-source model made of Crab-like sources.
655   - """
656   - # Add disk Crabs
657   - obs = add_source_model(obs, crab_like("Disk 0deg off", 0.0, 0.0, type="disk", radius=0.3))
658   - obs = add_source_model(obs, crab_like("Disk 1deg off", -1.0, 0.0, type="disk", radius=0.3))
659   - obs = add_source_model(obs, crab_like("Disk 2deg off", -2.0, 0.0, type="disk", radius=0.3))
660   - obs = add_source_model(obs, crab_like("Disk 3deg off", -3.0, 0.0, type="disk", radius=0.3))
661   - #obs = add_source_model(obs, crab_like("Disk 4deg off", -4.0, 0.0, type="disk", radius=0.3))
662   - #obs = add_source_model(obs, crab_like("Disk 5deg off", -5.0, 0.0, type="disk", radius=0.3))
663   -
664   - # Return
665   - return obs
666   -
667   -
668   -# ================== #
669   -# Setup Gauss models #
670   -# ================== #
671   -def add_gauss_sources(obs):
672   - """
673   - Add Gauss-source model made of Crab-like sources.
674   - """
675   - # Add Gauss Crabs
676   - obs = add_source_model(obs, crab_like("Gauss 0deg off", 0.0, 0.0, type="gauss", sigma=0.3))
677   - obs = add_source_model(obs, crab_like("Gauss 1deg off", -1.0, 0.0, type="gauss", sigma=0.3))
678   - obs = add_source_model(obs, crab_like("Gauss 2deg off", -2.0, 0.0, type="gauss", sigma=0.3))
679   - obs = add_source_model(obs, crab_like("Gauss 3deg off", -3.0, 0.0, type="gauss", sigma=0.3))
680   - #obs = add_source_model(obs, crab_like("Gauss 4deg off", -4.0, 0.0, type="gauss", sigma=0.3))
681   - #obs = add_source_model(obs, crab_like("Gauss 5deg off", -5.0, 0.0, type="gauss", sigma=0.3))
682   -
683   - # Return
684   - return obs
685   -
686   -
687   -# ================== #
688   -# Setup Shell models #
689   -# ================== #
690   -def add_shell_sources(obs):
691   - """
692   - Add shell-source model made of Crab-like sources.
693   - """
694   - # Add Shell Crabs
695   - obs = add_source_model(obs, crab_like("Shell 0deg off", 0.0, 0.0, type="shell", radius=0.3, width=0.1))
696   - obs = add_source_model(obs, crab_like("Shell 1deg off", -1.0, 0.0, type="shell", radius=0.3, width=0.1))
697   - obs = add_source_model(obs, crab_like("Shell 2deg off", -2.0, 0.0, type="shell", radius=0.3, width=0.1))
698   - #obs = add_source_model(obs, crab_like("Shell 3deg off", -3.0, 0.0, type="shell", radius=0.3, width=0.1))
699   - #obs = add_source_model(obs, crab_like("Shell 4deg off", -4.0, 0.0, type="shell", radius=0.3, width=0.1))
700   - #obs = add_source_model(obs, crab_like("Shell 5deg off", -5.0, 0.0, type="shell", radius=0.3, width=0.1))
701   -
702   - # Return
703   - return obs
704   -
705   -
706   -# ==================== #
707   -# Setup big disk model #
708   -# ==================== #
709   -def add_big_disk_source(obs):
710   - """
711   - Add big disk-source model.
712   - """
713   - # Add one big offset disk Crab
714   - obs = add_source_model(obs, crab_like("Big Disk", 1.0, 1.0, type="disk", radius=1.0))
715   -
716   - # Return
717   - return obs
718   -
719   -
720   -# ===================== #
721   -# Setup big Gauss model #
722   -# ===================== #
723   -def add_big_gauss_source(obs):
724   - """
725   - Add big gauss-source model.
726   - """
727   - # Add one big offset Gauss Crab
728   - obs = add_source_model(obs, crab_like("Big Gauss", -1.0, 1.0, type="gauss", sigma=1.0))
729   -
730   - # Return
731   - return obs
732   -
733   -
734   -# ===================== #
735   -# Setup big shell model #
736   -# ===================== #
737   -def add_big_shell_source(obs):
738   - """
739   - Add big shell-source model.
740   - """
741   - # Add one big offset Shell Crab
742   - obs = add_source_model(obs, crab_like("Big Shell", -1.0, -1.0, type="shell", radius=1.0, width=0.1))
743   -
744   - # Return
745   - return obs
746   -
747   -
748   -# ================================================= #
749   -# Setup set of point source for the survey analysis #
750   -# ================================================= #
751   -def add_survey_sources(obs):
752   - """
753   - Add point-source model made of Crab-like sources.
754   - """
755   - # Add point-source Crabs
756   - obs = add_source_model(obs, crab_like("Source 1", 0.0, 0.0, flux=0.1, type="point"))
757   - obs = add_source_model(obs, crab_like("Source 2", 1.0, 0.0, flux=0.2, type="point"))
758   - obs = add_source_model(obs, crab_like("Source 3", -1.0, 0.0, flux=0.5, type="point"))
759   - obs = add_source_model(obs, crab_like("Source 4", 0.0, 1.0, flux=1.0, type="point"))
760   - obs = add_source_model(obs, crab_like("Source 5", 0.0, -1.0, flux=2.0, type="point"))
761   - obs = add_source_model(obs, crab_like("Source 6", 0.0, 4.0, flux=5.0, type="point"))
762   -
763   - # Return
764   - return obs
765   -
766   -
767   -# ================================================= #
768   -# Setup set of point source for the survey analysis #
769   -# ================================================= #
770   -def add_random_sources(model):
771   - """
772   - Add point-sources uniformly distributed over the sky, with a parameterised distribution in flux.
773   -
774   - Arguments
775   - - model: a source model container to hold the source model over the loop
776   -
777   - Notes
778   - - Contrary to other add_source functions, this one works on a model container and not on an observation container
779   - - This is because we need to reuse the same random source model in a loop over pointing strategies
780   - """
781   - # Define source population characteristics
782   - flux_min=0.001
783   - flux_max=10.0
784   - flux_index=-2.0
785   - field=4.0
786   - l_min=0.0-0.5*field
787   - b_min=0.0-0.5*field
788   - num_sources=4
789   -
790   - # Loop over sources to be generated
791   - for i in range(0,num_sources):
792   -
793   - # Get random flux
794   - flux=math.log10(flux_min)+(math.log10(flux_max)-math.log10(flux_min))*random.random()
795   - flux=math.pow(10.0,flux)
796   -
797   - # Get random longitude/latitude
798   - l=l_min+field*random.random()
799   - b=b_min+field*random.random()
800   -
801   - # Add point-source Crabs
802   - sourcename="Source "+str(i)
803   - model.append(crab_like(sourcename, l, b, flux=flux, type="point"))
804   -
805   - # Return
806   - return model
807   -
808   -
809   -# ================================================= #
810   -# Setup a test point source for the survey analysis #
811   -# ================================================= #
812   -def add_test_sources(models, l=0.0, b=0.0, flux=1.0, index=-2.5):
813   - """
814   - Add a test point-source with specified position, flux, and spectral index
815   -
816   - Arguments
817   - - l: longitude
818   - - b: latitude
819   - - flux: flux in Crab units
820   - - index: spectral index
821   - """
822   -
823   - # Set source name
824   - sourcename="Test"
825   -
826   - # Add it to container
827   - models = add_source_model(models, any_pointlike(sourcename, l, b, \
828   - flux=flux, index=index, fitpos=False, fitidx=False))
829   -
830   - # Return
831   - return models
832   -
833   -
834   -#==============================#
835   -# Modified routine entry point #
836   -#==============================#
837   -if __name__ == '__main__':
838   - """
839   - CTA survey simulation and analysis script.
840   -
841   - Notes:
842   - - Based on the initial script
843   - - Added functions: survey_square, add_survey_sources, add_random_sources
844   - - Added import: math, random
845   -
846   - """
847   - # Set job parameters
848   - irf = "kb_E_50h_v3" # IRF
849   - survey_time = 50.0*3600.0 # Total survey time
850   - pnt_step = 1.0 # Spacing between pointings
851   - pnt_mintime = 3600.0 # Minimum duration of a pointing
852   - pnt_fov = 1.0 # Typical FoV size (useful to know where to probe sensitivity)
853   - pnt_maxnum = 1 # Maximum number of pointings
854   - emin = 0.1 # minimum energy (TeV)
855   - emax = 100.0 # maximum energy (TeV)
856   - #emin = 1.0 # minimum energy (TeV)
857   - #emax = 5.0 # maximum energy (TeV)
858   - nbins = 20 # number of logarithmic energy bins
859   - test_thres = 25.0 # Reference TS
860   - index = -3.0 # Spectral index
861   -
862   - # Initialise response database
863   - caldb = "/usr/local/gamma/share/caldb/cta"
864   -
865   - # Initialise flags
866   - need_help = False
867   -
868   - # Test for command line arguments
869   - if (len(sys.argv) > 1):
870   - surveyrange= float(sys.argv[1])
871   -
872   - if (len(sys.argv) > 2):
873   - surveystep= float(sys.argv[2])
874   -
875   - if (len(sys.argv) > 3):
876   - irf=sys.argv[3]
877   -
878   - # Print help if needed and exit
879   - if need_help:
880   - print "Usage: example_survey.py [OPTIONS]"
881   - print " galscan Simulate Galactic plane scan"
882   - print " extended Simulate also extended sources"
883   - print " -h Display this usage message"
884   - sys.exit()
885   -
886   - # Dump header
887   - print "***********************"
888   - print "* Simulate CTA survey *"
889   - print "***********************"
890   -
891   - # Remove any existing result files
892   - list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
893   - for files in list:
894   - for file in files:
895   - os.remove(file)
896   -
897   - # Set number of observations
898   - first = True
899   - results_file = 'sensitivity.dat'
900   - for i in range(1,pnt_maxnum+1):
901   -
902   - # Initialise models
903   - bkg_model = GModels()
904   - full_model = GModels()
905   -
906   - # Setup models
907   - bkg_model = add_background_model(bkg_model)
908   - full_model = add_background_model(full_model)
909   - full_model = add_test_sources(full_model, 0.0, 0.0, flux=1.0, index=index)
910   -
911   - # Compute duration of pointings and starting longitude
912   - pnt_duration = survey_time/float(i)
913   - pnt_lmin = 0.0-float(i-1)*pnt_step/2.0
914   -
915   - # Setup observations for survey
916   - survey = survey_galplane(caldb, irf, lmin=pnt_lmin, lstep=pnt_step, \
917   - nstep=i, duration=pnt_duration, \
918   - emin=emin, emax=emax)
919   -
920   - print ""
921   - print "Simulating an array of "+str(i)+" observations along the plane"
922   -
923   - # Look for smallest flux for detection
924   - test_num = 0
925   - test_max = 50
926   - test_flux = 0.001
927   -
928   - # Loop
929   - flux_value = []
930   - while True:
931   -
932   - # Update test counter
933   - test_num = test_num + 1
934   -
935   - # Set test flux
936   - full_model[1]['Prefactor'].value(test_flux)
937   - survey.models(full_model)
938   -
939   - # Simulate events
940   - obs = sim_obs(survey, nbins=nbins, seed=test_num)
941   -
942   - # Determine number of events in simulation
943   - nevents = 0.0
944   - for x in obs:
945   - nevents += x.events().number()
946   -
947   - # Fit background only
948   - obs.models(bkg_model)
949   - like = fit_obs(obs,log=False)
950   -
951   - # Store results
952   - LogL_bgm = like.opt().value()
953   - Pre_bgm = like.obs().models()[0]['Prefactor'].value()
954   - ePre_bgm = like.obs().models()[0]['Prefactor'].error()
955   - Inx_bgm = like.obs().models()[0]['Index'].value()
956   - eInx_bgm = like.obs().models()[0]['Index'].error()
957   - npred_bgm = like.obs().npred()
958   -
959   - # Fit background and test source
960   - obs.models(full_model)
961   - like = fit_obs(obs,log=False)
962   -
963   - # Store results
964   - LogL_all = like.opt().value()
965   - test_ts = 2.0*(LogL_bgm-LogL_all)
966   - Pre_all = like.obs().models()[0]['Prefactor'].value()
967   - ePre_all = like.obs().models()[0]['Prefactor'].error()
968   - Inx_all = like.obs().models()[0]['Index'].value()
969   - eInx_all = like.obs().models()[0]['Index'].error()
970   - SPre_all = like.obs().models()[1]['Prefactor'].value()
971   - eSPre_all = like.obs().models()[1]['Prefactor'].error()
972   - SInx_all = like.obs().models()[1]['Index'].value()
973   - eSInx_all = like.obs().models()[1]['Index'].error()
974   - npred_all = like.obs().npred()
975   - if test_ts > 0:
976   - flux = math.sqrt(test_thres/test_ts)*SPre_all
977   - else:
978   - flux = 0.0
979   - flux_value.append(flux)
980   -
981   - # Compute sliding average
982   - flux = 0.0
983   - num = 0.0
984   - for k in range(3):
985   - inx = test_num - k - 1
986   - if inx >= 0:
987   - flux += flux_value[inx]
988   - num += 1.0
989   - flux /= num
990   -
991   - # Print results
992   - print "TS=%9.2f %8.5f lnL(BGD)=%8.1f lnL(ALL)=%8.1f #%5d" \
993   - " BGD: #%8.2f %5.2f+/-%4.2f %5.2f+/-%5.3f" \
994   - " ALL: #%8.2f %5.2f+/-%4.2f %5.2f+/-%5.3f" \
995   - " %8.5f+/-%7.5f %5.2f+/-%5.3f Flux=%8.5f" % \
996   - (test_ts, test_flux, LogL_bgm, LogL_all, nevents, \
997   - npred_bgm, Pre_bgm, ePre_bgm, Inx_bgm, eInx_bgm, \
998   - npred_all, Pre_all, ePre_all, Inx_all, eInx_all, \
999   - SPre_all, eSPre_all, SInx_all, eSInx_all, \
1000   - flux)
1001   -
1002   - # Compare average flux to last average
1003   - if test_num > 3:
1004   - ratio = flux/test_flux
1005   - if ratio >= 0.99 and ratio <= 1.01:
1006   - print "Converged ("+str(ratio)+")"
1007   - break
1008   -
1009   - # Use average for next iteration
1010   - test_flux = flux
1011   -
1012   - # Exit loop if number of trials exhausted
1013   - if (test_num >= test_max):
1014   - print "Test ended after "+str(test_max)+" iterations."
1015   - break
1016   -
1017   - # Print results
1018   - print "Sensitivity of "+str(i)+" observations: "+str(flux*1000.0)+" mCrab"
1019   -
1020   - # Write to CSV file
1021   - result = {}
1022   - colnames = ['emin', 'emax', 'observations', 'sensitivity']
1023   - result['emin'] = emin
1024   - result['emax'] = emax
1025   - result['observations'] = i
1026   - result['sensitivity'] = flux
1027   - if first:
1028   - f = open(results_file, 'w')
1029   - writer = csv.DictWriter(f, colnames)
1030   - writer.writerow(dict((_,_) for _ in colnames))
1031   - first = False
1032   - else:
1033   - f = open(results_file, 'a')
1034   - writer = csv.DictWriter(f, colnames)
1035   - writer.writerow(result)
1036   - f.close()
1037   -
1038   - # Delete some objects
1039   - del obs
1040   - del like
1041   -
1042   -
... ...