Eirik Hatlen / ctools

Commit 329b821ae01122da81141bf6d8d66fe458737ea9

Authored by Josh Cardenzana
1 parent d403041f

(#2925) Add 'get_butterfly_file()' and 'get_butterfly_csv()' methods 

Allows extracting butterfly parameters for plotting
Inline Side-by-side
Showing 1 changed file with 59 additions and 20 deletions
examples/show_butterfly.py
  View file @329b821
... ... @@ -30,28 +30,45 @@ except (ImportError, RuntimeError):
30 30 sys.exit()
31 31  
32 32  
33   -# ============== #
34   -# Plot butterfly #
35   -# ============== #
36   -def plot_butterfly(filename, plotfile):
  33 +def get_butterfly_file(filename):
37 34 """
38   - Plot butterfly diagram
  35 + Get butterfly diagram values
39 36  
40 37 Parameters
41 38 ----------
42 39 filename : str
43 40 Butterfly CSV file
44   - plotfile : str
45   - Plot filename
  41 +
  42 + Returns
  43 + -------
  44 + See `get_butterfly_csv`
46 45 """
47 46 # Read butterfly file
48 47 csv = gammalib.GCsv(filename)
49 48  
  49 + return get_butterfly_csv(csv)
  50 +
  51 +
  52 +def get_butterfly_csv(csv):
  53 + """
  54 + Get butterfly diagram values from GCsv object
  55 +
  56 + Parameters
  57 + ----------
  58 + csv : `~gammalib.GCsv`
  59 + GCsv object containing butterfly plot information
  60 +
  61 + Returns
  62 + -------
  63 + 'dict' defining butterfly plot and best fit spectrum
  64 + """
50 65 # Initialise arrays to be filled
51   - butterfly_x = []
52   - butterfly_y = []
53   - line_x = []
54   - line_y = []
  66 + btrfly = {
  67 + 'butterfly_x' : [],
  68 + 'butterfly_y' : [],
  69 + 'line_x' : [],
  70 + 'line_y' : []
  71 + }
55 72  
56 73 # Loop over rows of the file
57 74 nrows = csv.nrows()
... ... @@ -61,30 +78,52 @@ def plot_butterfly(filename, plotfile):
61 78 conv = csv.real(row,0) * csv.real(row,0) * gammalib.MeV2erg
62 79  
63 80 # Compute upper edge of confidence band
64   - butterfly_x.append(csv.real(row,0)/1.0e6) # TeV
65   - butterfly_y.append(csv.real(row,2)*conv)
  81 + btrfly['butterfly_x'].append(csv.real(row,0)/1.0e6) # TeV
  82 + btrfly['butterfly_y'].append(csv.real(row,2)*conv)
66 83  
67 84 # Set line values
68   - line_x.append(csv.real(row,0)/1.0e6) # TeV
69   - line_y.append(csv.real(row,1)*conv)
  85 + btrfly['line_x'].append(csv.real(row,0)/1.0e6) # TeV
  86 + btrfly['line_y'].append(csv.real(row,1)*conv)
70 87  
71 88 # Loop over the rows backwards to compute the lower edge of the
72 89 # confidence band
73 90 for row in range(nrows):
74 91 index = nrows - 1 - row
75 92 conv = csv.real(index,0) * csv.real(index,0) * gammalib.MeV2erg
76   - butterfly_x.append(csv.real(index,0)/1.0e6)
  93 + btrfly['butterfly_x'].append(csv.real(index,0)/1.0e6)
77 94 low_error = max(csv.real(index,3)*conv, 1e-26)
78   - butterfly_y.append(low_error)
  95 + btrfly['butterfly_y'].append(low_error)
  96 +
  97 + return btrfly
  98 +
  99 +
  100 +# ============== #
  101 +# Plot butterfly #
  102 +# ============== #
  103 +def plot_butterfly(filename, plotfile):
  104 + """
  105 + Plot butterfly diagram
  106 +
  107 + Parameters
  108 + ----------
  109 + filename : str
  110 + Butterfly CSV file
  111 + plotfile : str
  112 + Plot filename
  113 + """
  114 + # Read butterfly file
  115 + btrfly = get_butterfly_file(filename)
79 116  
80 117 # Plot the butterfly and spectral line
81 118 plt.figure()
82 119 plt.loglog()
83 120 plt.grid()
84   - plt.plot(line_x,line_y,color='black',ls='-')
85   - plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5)
86 121 plt.xlabel('Energy (TeV)')
87 122 plt.ylabel(r'E$^2$ $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)')
  123 +
  124 + # Plot the butterfly and spectral line
  125 + plt.plot(btrfly['line_x'],btrfly['line_y'],color='black',ls='-')
  126 + plt.fill(btrfly['butterfly_x'],btrfly['butterfly_y'],color='green',alpha=0.5)
88 127  
89 128 # Show spectrum or save it into file
90 129 if len(plotfile) > 0:
... ... @@ -116,7 +155,7 @@ def show_butterfly():
116 155 plotfile = options[0]['value']
117 156  
118 157 # Plot butterfly diagram
119   - plot_butterfly(args[0], plotfile)
  158 + plot_butterfly(filename=args[0], plotfile=plotfile)
120 159  
121 160 # Return
122 161 return
... ...