Michael Mayer / ctools

Commit b3ecfbf9714058dc1652c21d3e87fe2889d02239

Authored by Jürgen Knödlseder
1 parent 44e28016

Correctly implement flux computation in sensitivity estimation

Inline Side-by-side
Showing 1 changed file with 11 additions and 2 deletions
examples/show_response.py
  View file @b3ecfbf
@@ -167,6 +167,7 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig @@ -167,6 +167,7 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig
167 """ 167 """
168 # Set constants 168 # Set constants
169 r68_to_sigma = 0.6624305 169 r68_to_sigma = 0.6624305
  170 + TeV2erg = 1.0e6 * gammalib.MeV2erg
170 171
171 # Set to figure 3 172 # Set to figure 3
172 plt.figure(3) 173 plt.figure(3)
@@ -189,8 +190,16 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig @@ -189,8 +190,16 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig
189 if Noff > 0: 190 if Noff > 0:
190 Non = Non_lima(sigma, Noff) 191 Non = Non_lima(sigma, Noff)
191 src_counts = Non - bgd_counts 192 src_counts = Non - bgd_counts
192 - flux[i] = src_counts / (0.68*aeff * duration) * E[i]*1.0e6 * gammalib.MeV2erg  
193 - #print logE[i], flux[i] 193 + if src_counts < 0.05*bgd_counts:
  194 + src_counts = 0.05*bgd_counts
  195 + if src_counts < 10:
  196 + src_counts = 10.0
  197 + emin = gammalib.GEnergy(math.pow(10.0, logE[i]-0.1), "TeV")
  198 + emax = gammalib.GEnergy(math.pow(10.0, logE[i]+0.1), "TeV")
  199 + epivot = gammalib.GEnergy(math.pow(10.0, logE[i]), "TeV")
  200 + plaw = gammalib.GModelSpectralPlaw(1.0e-6, -2.6, epivot)
  201 + conv = TeV2erg*E[i]*E[i]/plaw.flux(emin, emax)
  202 + flux[i] = conv * src_counts / (duration * aeff*0.68)
194 203
195 # Plot data 204 # Plot data
196 plt.loglog(E, flux, color+'-', label=name) 205 plt.loglog(E, flux, color+'-', label=name)