Jürgen Knödlseder / ctools

Commit bd73d6fa801425555759f4adaca89797495b365a

Authored by Jürgen Knödlseder
1 parent 58e558e8

Remove cast and adapt to Python 3.

Inline Side-by-side
Showing 1 changed file with 140 additions and 140 deletions
examples/make_unbinned_analysis.py
  View file @bd73d6f
... ... @@ -29,6 +29,11 @@ from math import *
29 29 import os
30 30 import glob
31 31 import sys
  32 +try:
  33 + import matplotlib.pyplot as plt
  34 + has_matplotlib = True
  35 +except:
  36 + has_matplotlib = False
32 37  
33 38  
34 39 # ============================= #
... ... @@ -78,8 +83,8 @@ def pipeline_v1():
78 83 sim["emin"].real(emin)
79 84 sim["emax"].real(emax)
80 85 sim.execute()
81   - print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
82   -
  86 + print("Simulated events ("+str(sim.celapse())+" CPU seconds)")
  87 +
83 88 # Update timing
84 89 wall_seconds += sim.telapse()
85 90 cpu_seconds += sim.celapse()
... ... @@ -97,7 +102,7 @@ def pipeline_v1():
97 102 select["emin"].real(emin)
98 103 select["emax"].real(emax)
99 104 select.execute()
100   - print "Select event ("+str(select.celapse())+" CPU seconds)"
  105 + print("Select event ("+str(select.celapse())+" CPU seconds)")
101 106  
102 107 # Update timing
103 108 wall_seconds += select.telapse()
... ... @@ -112,15 +117,15 @@ def pipeline_v1():
112 117 like["caldb"].string(caldb)
113 118 like["irf"].string(irf)
114 119 like.execute()
115   - print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
  120 + print("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)")
116 121  
117 122 # Update timing
118 123 wall_seconds += like.telapse()
119 124 cpu_seconds += like.celapse()
120 125  
121 126 # Show total times
122   - print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
123   - print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
  127 + print("Total wall time elapsed: "+str(wall_seconds)+" seconds")
  128 + print("Total CPU time used ...: "+str(cpu_seconds)+" seconds")
124 129  
125 130 # Return
126 131 return
... ... @@ -167,7 +172,7 @@ def pipeline_v2(show_data):
167 172 sim["emin"].real(emin)
168 173 sim["emax"].real(emax)
169 174 sim.run()
170   - print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
  175 + print("Simulated events ("+str(sim.celapse())+" CPU seconds)")
171 176  
172 177 # Update timing
173 178 wall_seconds += sim.telapse()
... ... @@ -183,7 +188,7 @@ def pipeline_v2(show_data):
183 188 select["emin"].real(emin)
184 189 select["emax"].real(emax)
185 190 select.run()
186   - print "Select event ("+str(select.celapse())+" CPU seconds)"
  191 + print("Select event ("+str(select.celapse())+" CPU seconds)")
187 192  
188 193 # Update timing
189 194 wall_seconds += select.telapse()
... ... @@ -192,22 +197,25 @@ def pipeline_v2(show_data):
192 197 # Perform maximum likelihood fitting
193 198 like = ctlike(select.obs())
194 199 like.run()
195   - print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
  200 + print("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)")
196 201  
197 202 # Update timing
198 203 wall_seconds += like.telapse()
199 204 cpu_seconds += like.celapse()
200 205  
201 206 # Show total times
202   - print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
203   - print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
  207 + print("Total wall time elapsed: "+str(wall_seconds)+" seconds")
  208 + print("Total CPU time used ...: "+str(cpu_seconds)+" seconds")
204 209  
205 210 # Show model fitting results
206 211 #print like.obs().models()
207 212  
208 213 # Optionally plot counts
209 214 if show_data:
210   - plot_counts(like.obs())
  215 + if has_matplotlib:
  216 + plot_counts(like.obs())
  217 + else:
  218 + sys.stdout.write("Matplotlib is not (correctly) installed on your system. No counts spectra are shown.\n")
211 219  
212 220 # Return
213 221 return
... ... @@ -220,143 +228,135 @@ def plot_counts(observations):
220 228 """
221 229 Plot counts.
222 230 """
223   - # Only proceed if matplotlib is available
224   - try:
225   - # Import matplotlib
226   - import matplotlib.pyplot as plt
227   -
228   - # Set legend fontsize
229   - params = {'legend.fontsize': 10}
230   - plt.rcParams.update(params)
231   -
232   - # Set plot styles
233   - styles = ['b-', 'g-', 'y-', 'n-']
234   -
235   - # Dump header
236   - print ""
237   - print "Make plots (using matplotlib):"
238   - print "=============================="
239   -
240   - # Create figure 1
241   - plt.figure(1,figsize=(12,6))
242   - plt.subplots_adjust(hspace=.7)
  231 + # Set legend fontsize
  232 + params = {'legend.fontsize': 10}
  233 + plt.rcParams.update(params)
  234 +
  235 + # Set plot styles
  236 + styles = ['b-', 'g-', 'y-', 'n-']
243 237  
244   - # Create subplot 1
245   - plt.subplot(121)
246   - plt.title("Spectrum (summed over all pixels)")
  238 + # Dump header
  239 + print("")
  240 + print("Make plots (using matplotlib):")
  241 + print("==============================")
  242 +
  243 + # Create figure 1
  244 + plt.figure(1,figsize=(12,6))
  245 + plt.subplots_adjust(hspace=.7)
  246 +
  247 + # Create subplot 1
  248 + plt.subplot(121)
  249 + plt.title("Spectrum (summed over all pixels)")
  250 +
  251 + # Loop over observations
  252 + for obs in observations:
  253 +
  254 + # Get event list
  255 + list = obs.events()
247 256  
248   - # Loop over observations
249   - for obs in observations:
  257 + # Create energy axis
  258 + ebounds = GEbounds()
  259 + emin = GEnergy()
  260 + emax = GEnergy()
  261 + emin.TeV(0.1)
  262 + emax.TeV(100.0)
  263 + ebounds.setlog(emin, emax, 10)
  264 + energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
  265 +
  266 + # Create spectrum
  267 + print("Extract data:")
  268 + counts = [0.0 for i in range(ebounds.size())]
  269 + print(list.size())
  270 + for atom in list:
  271 + index = ebounds.index(atom.energy())
  272 + counts[index] = counts[index] + 1.0
  273 +
  274 + # Create error bars
  275 + error = [sqrt(c) for c in counts]
  276 +
  277 + # Plot spectrum
  278 + plt.loglog(energy, counts, 'ro', label='data')
  279 + #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
  280 +
  281 + # Extract models
  282 + #print("Extract models:")
  283 + #sum_model = [0.0 for i in range(ebounds.size())]
  284 + #for k, m in enumerate(observations.models()):
  285 + # print("- "+m.name())
  286 + # model = [0.0 for i in range(ebounds.size())]
  287 + # for atom in list:
  288 + # index = ebounds.index(atom.energy())
  289 + # prob = m.eval(atom, obs)
  290 + # model[index] = model[index] + prob * atom.size()
  291 + # for i in range(ebounds.size()):
  292 + # sum_model[i] = sum_model[i] + model[i]
  293 + # plt.loglog(energy, model, styles[k], label=m.name())
  294 + #plt.loglog(energy, sum_model, 'r-', label='total')
  295 +
  296 + # Put labels
  297 + plt.xlabel("Energy (TeV)")
  298 + plt.ylabel("Counts")
  299 + plt.legend(loc="lower left")
  300 +
  301 + # Create subplot 2
  302 + plt.subplot(122)
  303 + plt.title("Offset (summed over all energies)")
  304 +
  305 + # Set Crab direction
  306 + crab = GSkyDir()
  307 + crab.radec_deg(83.63, 22.01)
  308 +
  309 + # Loop over observations
  310 + for obs in observations:
  311 +
  312 + # Get event cube
  313 + list = obs.events()
250 314  
251   - # Get event list
252   - list = cast_GCTAEventList(obs.events())
253   -
254   - # Create energy axis
255   - ebounds = GEbounds()
256   - emin = GEnergy()
257   - emax = GEnergy()
258   - emin.TeV(0.1)
259   - emax.TeV(100.0)
260   - ebounds.setlog(emin, emax, 10)
261   - energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
262   -
263   - # Create spectrum
264   - print "Extract data:"
265   - counts = [0.0 for i in range(ebounds.size())]
266   - print list.size()
267   - for atom in list:
268   - index = ebounds.index(atom.energy())
  315 + # Create offset histogram
  316 + print "Extract data:"
  317 + nx = 30
  318 + doffset2 = 0.01
  319 + offset2 = [(i+0.5)*doffset2 for i in range(nx)]
  320 + counts = [0.0 for i in range(nx)]
  321 + for atom in list:
  322 + off = GCTAInstDir(atom.dir()).dist_deg(crab)
  323 + off2 = off*off
  324 + index = int(off2/doffset2)
  325 + if index < nx:
269 326 counts[index] = counts[index] + 1.0
270 327  
271   - # Create error bars
272   - error = [sqrt(c) for c in counts]
273   -
274   - # Plot spectrum
275   - plt.loglog(energy, counts, 'ro', label='data')
276   - #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
277   -
278   - # Extract models
279   - print "Extract models:"
280   - sum_model = [0.0 for i in range(ebounds.size())]
281   - for k, m in enumerate(observations.models()):
282   - print "- "+m.name()
283   - model = [0.0 for i in range(ebounds.size())]
284   - for atom in list:
285   - index = ebounds.index(atom.energy())
286   - prob = m.eval(atom, obs)
287   - model[index] = model[index] + prob * atom.size()
288   - for i in range(ebounds.size()):
289   - sum_model[i] = sum_model[i] + model[i]
290   - #plt.loglog(energy, model, styles[k], label=m.name())
291   - #plt.loglog(energy, sum_model, 'r-', label='total')
292   -
293   - # Put labels
294   - plt.xlabel("Energy (TeV)")
295   - plt.ylabel("Counts")
296   - plt.legend(loc="lower left")
297   -
298   - # Create subplot 2
299   - plt.subplot(122)
300   - plt.title("Offset (summed over all energies)")
301   -
302   - # Set Crab direction
303   - crab = GSkyDir()
304   - crab.radec_deg(83.63, 22.01)
305   -
306   - # Loop over observations
307   - for obs in observations:
308   -
309   - # Get event cube
310   - list = cast_GCTAEventList(obs.events())
311   -
312   - # Create offset histogram
313   - print "Extract data:"
314   - nx = 30
315   - doffset2 = 0.01
316   - offset2 = [(i+0.5)*doffset2 for i in range(nx)]
317   - counts = [0.0 for i in range(nx)]
  328 + # Create error bars
  329 + error = [sqrt(c) for c in counts]
  330 +
  331 + # Plot distribution
  332 + plt.semilogy(offset2, counts, 'ro', label='data')
  333 +
  334 + # Extract models
  335 + print "Extract models:"
  336 + sum_model = [0.0 for i in range(nx)]
  337 + for k, m in enumerate(observations.models()):
  338 + print "- "+m.name()
  339 + model = [0.0 for i in range(nx)]
318 340 for atom in list:
319   - off = atom.dir().dist_deg(crab)
  341 + off = GCTAInstDir(atom.dir()).dist_deg(crab)
320 342 off2 = off*off
321 343 index = int(off2/doffset2)
322 344 if index < nx:
323   - counts[index] = counts[index] + 1.0
324   -
325   - # Create error bars
326   - error = [sqrt(c) for c in counts]
327   -
328   - # Plot distribution
329   - plt.semilogy(offset2, counts, 'ro', label='data')
330   -
331   - # Extract models
332   - print "Extract models:"
333   - sum_model = [0.0 for i in range(nx)]
334   - for k, m in enumerate(observations.models()):
335   - print "- "+m.name()
336   - model = [0.0 for i in range(nx)]
337   - for atom in list:
338   - off = atom.dir().dist_deg(crab)
339   - off2 = off*off
340   - index = int(off2/doffset2)
341   - if index < nx:
342   - prob = m.eval(atom, obs)
343   - model[index] = model[index] + prob * atom.size()
344   - for i in range(nx):
345   - sum_model[i] = sum_model[i] + model[i]
346   - #plt.plot(offset2, model, styles[k], label=m.name())
347   - #plt.plot(offset2, sum_model, 'r-', label='total')
348   - #plt.ylim(ymin=0.1)
349   -
350   - # Put labels
351   - plt.xlabel("Offset (deg^2)")
352   - plt.ylabel("Counts")
353   - plt.legend(loc="upper right")
354   -
355   - # Show counts spectra
356   - plt.show()
357   -
358   - except ImportError:
359   - print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown."
  345 + prob = m.eval(atom, obs)
  346 + model[index] = model[index] + prob * atom.size()
  347 + for i in range(nx):
  348 + sum_model[i] = sum_model[i] + model[i]
  349 + #plt.plot(offset2, model, styles[k], label=m.name())
  350 + #plt.plot(offset2, sum_model, 'r-', label='total')
  351 + #plt.ylim(ymin=0.1)
  352 +
  353 + # Put labels
  354 + plt.xlabel("Offset (deg^2)")
  355 + plt.ylabel("Counts")
  356 + plt.legend(loc="upper right")
  357 +
  358 + # Show counts spectra
  359 + plt.show()
360 360  
361 361 # Return
362 362 return
... ...