In [1]:
%pylab inline


Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

In [13]:
rcParams['font.size'] = 20.

In [2]:
from VOT import *

In [3]:
from VOTPlot import Plot, PlotSpectrum

This initial GRB is the nominal 080916C.


In [4]:
baseline = VOT("custom", eMin=0.001, eMax=100000, Nbins=10000, redshift=4.35, eblModel='Dominguez', instrument='VERITAS', zenith=20, 
        spectralModel='Band', N0=1e-9,alpha=-0.65, beta=-2.22, Ep = (2 - 0.65)*(265./100000.))


WARNING - Using maximum redshift in EBL table.
WARNING - Spectrum starts below minimum energy in EBL table.
WARNING - Spectrum continues beyond maximum energy in EBL table.
WARNING - Using 659.0 counts/hr as the rate from the Crab Nebula.
WARNING - Time to detection assumes a Crab Nebula Spectrum.
INFO - Using Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_soft-1.summary.csv
INFO - Using EffectiveArea_Azimuth_0_Zenith_20_Noise_4.07
INFO - Safe energy range: 177.83 to 25,118.86 GeV
INFO - dNdE at 1 GeV: 8.52e-06 s^-1 cm^-2 GeV^-1
INFO - dNdE at 400 GeV: 1.42e-11 s^-1 cm^-2 GeV^-1
INFO - dNdE at 1 TeV: 1.86e-12 s^-1 cm^-2 GeV^-1
INFO - tau at min safe E: 5.37
INFO - tau at max safe E: 3,150.77
INFO - Predicted counts/hour: 9.71e+00
INFO - This is approximately 1.4731% of the Crab Nebula's Flux
INFO - This will take approximately 12.3949 hours to detect at a 5 sigma level

You can see that if this GRB went off in the FoV of VERITAS at a good zenith angle it could be detected in ~12.5 hours. This is fine for a non-transient source but not that great for a GRB. What's really killing us here is the EBL. This GRB was at a redshift of 4.35. What if it was at a redshift of 1.0?


In [5]:
baseline_z1 = VOT("custom", eMin=0.001, eMax=100000, Nbins=10000, redshift=1.0, eblModel='Dominguez', instrument='VERITAS', zenith=20, 
        spectralModel='Band', N0=1e-9,alpha=-0.65, beta=-2.22, Ep = (2 - 0.65)*(265./100000.))


WARNING - Spectrum starts below minimum energy in EBL table.
WARNING - Spectrum continues beyond maximum energy in EBL table.
WARNING - Using 659.0 counts/hr as the rate from the Crab Nebula.
WARNING - Time to detection assumes a Crab Nebula Spectrum.
INFO - Using Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_soft-1.summary.csv
INFO - Using EffectiveArea_Azimuth_0_Zenith_20_Noise_4.07
INFO - Safe energy range: 177.83 to 25,118.86 GeV
INFO - dNdE at 1 GeV: 8.52e-06 s^-1 cm^-2 GeV^-1
INFO - dNdE at 400 GeV: 1.42e-11 s^-1 cm^-2 GeV^-1
INFO - dNdE at 1 TeV: 1.86e-12 s^-1 cm^-2 GeV^-1
INFO - tau at min safe E: 2.19
INFO - tau at max safe E: 1,313.03
INFO - Predicted counts/hour: 4.15e+02
INFO - This is approximately 63.0119% of the Crab Nebula's Flux
INFO - This will take approximately 0.0543 hours to detect at a 5 sigma level

This is much more promising. So if you had a burst like 080916C go off at a redshift of 1 and it was in a good location for a IACT then you would nail the detection in about 10 minutes.


In [23]:
fig2 = pyplot.figure(figsize=(16,8))
fig2ax1 = fig2.add_subplot(111)
fig2ax1.set_ylabel(r'Differential Flux [cm$^{-2}$ s$^{-1}$ GeV$^{-1}$]')
fig2ax1.set_xlabel('Energy [GeV]')
fig2ax1.loglog(baseline.VS.EBins, baseline.VS.dNdE)
fig2ax1.loglog(baseline.VS.EBins, baseline.VS.dNdE_absorbed)
fig2ax1.loglog(baseline_z1.VS.EBins, baseline_z1.VS.dNdE_absorbed)
fig2ax1.set_ylim((1e-16,1e-6))
fig2ax1.set_xlim([1,1000])
fig2ax2 = fig2ax1.twinx()
fig2ax2.loglog((10**baseline.VR.EACurve[0:,0])[::2],(baseline.VR.EACurve[0:,1])[::2],'ko')
fig2ax2.set_ylabel(r'Effective Area [cm$^2$]')
fig2ax2.axvline(10**baseline.VR.EASummary['minSafeE'],color='k')
fig2ax2.axvline(10**baseline.VR.EASummary['maxSafeE'],color='k')
show()


This plot shows the effect of the EBL. The blue line is the raw band model, the red is the spectrum after accounting for a redshift of 4.35 and the green is for a redshift of 1.0. The black curve is the VERITAS EA at a zenith angle of 20 degrees. The vertical black lines are the low and high energy ranges for the VERITAS analysis.

Another issue is that most GRBs for occur at high zenith angles due to solid angle arguments and high zenith means high energy threashold for the Cherenkov technique. What if this GRB occurred at a zenith angle of 45 degrees.


In [16]:
baseline_zenith50 = VOT("custom", eMin=0.001, eMax=100000, Nbins=10000, redshift=4.35, eblModel='Dominguez', instrument='VERITAS', zenith=50, 
        spectralModel='Band', N0=1e-9,alpha=-0.65, beta=-2.22, Ep = (2 - 0.65)*(265./100000.))


WARNING - Using maximum redshift in EBL table.
WARNING - Spectrum starts below minimum energy in EBL table.
WARNING - Spectrum continues beyond maximum energy in EBL table.
WARNING - Using 659.0 counts/hr as the rate from the Crab Nebula.
WARNING - Time to detection assumes a Crab Nebula Spectrum.
INFO - Using Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_soft-1.summary.csv
INFO - Using EffectiveArea_Azimuth_0_Zenith_50_Noise_4.07
INFO - Safe energy range: 630.96 to 56,234.13 GeV
INFO - dNdE at 1 GeV: 8.52e-06 s^-1 cm^-2 GeV^-1
INFO - dNdE at 400 GeV: 1.42e-11 s^-1 cm^-2 GeV^-1
INFO - dNdE at 1 TeV: 1.86e-12 s^-1 cm^-2 GeV^-1
INFO - tau at min safe E: 19.33
INFO - tau at max safe E: 3,708.95
INFO - Predicted counts/hour: 3.78e-06
INFO - This is approximately 0.0000% of the Crab Nebula's Flux
INFO - This will take approximately 241.6449 hours to detect at a 5 sigma level

This thing is not detectable if it happened at a zenith angle of 50. In fact, it'll only be detectable at a zenith angle of 30 and if it was nearby.


In [17]:
baseline_zenith30_z1 = VOT("custom", eMin=0.001, eMax=100000, Nbins=10000, redshift=1.0, eblModel='Dominguez', instrument='VERITAS', zenith=30, 
        spectralModel='Band', N0=1e-9,alpha=-0.65, beta=-2.22, Ep = (2 - 0.65)*(265./100000.))


WARNING - Spectrum starts below minimum energy in EBL table.
WARNING - Spectrum continues beyond maximum energy in EBL table.
WARNING - Using 659.0 counts/hr as the rate from the Crab Nebula.
WARNING - Time to detection assumes a Crab Nebula Spectrum.
INFO - Using Effective_Areas/VERITAS/ea_Nov2010_na_ATM21_vegasv240rc1_7sam_050off_soft-1.summary.csv
INFO - Using EffectiveArea_Azimuth_0_Zenith_30_Noise_4.07
INFO - Safe energy range: 223.87 to 31,622.78 GeV
INFO - dNdE at 1 GeV: 8.52e-06 s^-1 cm^-2 GeV^-1
INFO - dNdE at 400 GeV: 1.42e-11 s^-1 cm^-2 GeV^-1
INFO - dNdE at 1 TeV: 1.86e-12 s^-1 cm^-2 GeV^-1
INFO - tau at min safe E: 3.24
INFO - tau at max safe E: 1,668.79
INFO - Predicted counts/hour: 1.15e+02
INFO - This is approximately 17.3907% of the Crab Nebula's Flux
INFO - This will take approximately 0.2131 hours to detect at a 5 sigma level

This plot shows why...


In [25]:
fig2 = pyplot.figure(figsize=(16,8))
fig2ax1 = fig2.add_subplot(111)
fig2ax1.set_ylabel(r'Differential Flux [cm$^{-2}$ s$^{-1}$ GeV$^{-1}$]')
fig2ax1.set_xlabel('Energy [GeV]')
fig2ax1.loglog(baseline.VS.EBins, baseline.VS.dNdE,'k')
fig2ax1.loglog(baseline.VS.EBins, baseline.VS.dNdE_absorbed,'k')
fig2ax1.loglog(baseline_z1.VS.EBins, baseline_z1.VS.dNdE_absorbed,'k')

fig2ax1.set_ylim((1e-16,1e-6))
fig2ax1.set_xlim([1,1000])
fig2ax2 = fig2ax1.twinx()
fig2ax2.set_ylabel(r'Effective Area [cm$^2$]')

fig2ax2.loglog((10**baseline.VR.EACurve[0:,0])[::2],(baseline.VR.EACurve[0:,1])[::2],'bo')
fig2ax2.axvline(10**baseline.VR.EASummary['minSafeE'],color='b')
fig2ax2.axvline(10**baseline.VR.EASummary['maxSafeE'],color='b')

fig2ax2.loglog((10**baseline_zenith30_z1.VR.EACurve[0:,0])[::2],(baseline_zenith30_z1.VR.EACurve[0:,1])[::2],'ro')
fig2ax2.axvline(10**baseline_zenith30_z1.VR.EASummary['minSafeE'],color='r')
fig2ax2.axvline(10**baseline_zenith30_z1.VR.EASummary['maxSafeE'],color='r')


show()


The black curves are what were shown before (band model, absorbed at z = 4.35, and absorbed at z = 1.0). The blue curve and vertical lines are the EA at a zenith angle of 20 while the red is the same at a zenith angle of 30. The high energy EA (above ~300 GeV) is the same but the low energy threshold changes by about 50 GeV which is huge at these redshifts.

Things we will do:

  • Add additional components
  • Comprehensive instrument analysis (HAWC, VERITAS, HESS, HESS2, CTA…)
  • Statistical test for solid angle
  • Light curves
  • Different EBL models
  • ...

In [6]:
baseline_hawc = VOT("custom", eMin=0.001, eMax=100000, Nbins=10000, redshift=4.35, eblModel='Dominguez', instrument='HAWC', zenith=20, 
        spectralModel='Band', N0=1e-9,alpha=-0.65, beta=-2.22, Ep = (2 - 0.65)*(265./100000.))


WARNING - Using maximum redshift in EBL table.
WARNING - Spectrum starts below minimum energy in EBL table.
WARNING - Spectrum continues beyond maximum energy in EBL table.
CRITICAL - Unsupported instrument.
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-6-e4dedadb88b7> in <module>()
      1 baseline_hawc = VOT("custom", eMin=0.001, eMax=100000, Nbins=10000, redshift=4.35, eblModel='Dominguez', instrument='HAWC', zenith=20, 
----> 2         spectralModel='Band', N0=1e-9,alpha=-0.65, beta=-2.22, Ep = (2 - 0.65)*(265./100000.))

/Users/jsperki1/src/VHEObserversTools/VOT.pyc in __init__(self, source, input, output, jsonString, dataOut, sparseData, **pars)
     40                    "custom" : self.calculateRateAndTime}
     41 
---> 42         sources[source](**pars)
     43         self.Print(output, dataOut, sparseData)
     44 

/Users/jsperki1/src/VHEObserversTools/VOT.pyc in calculateRateAndTime(self, eMin, eMax, Nbins, redshift, eblModel, instrument, zenith, spectralModel, **spectralPars)
     96 
     97         self.VR = VOTResponse(instrument,zenith=zenith, azimuth=0, noise=4.07)
---> 98         if(not self.VR.EASummary):
     99             return
    100         self.EACurve_interpolated = self.VR.interpolateEA(self.VS.EBins, self.VR.EACurve)

AttributeError: VOTResponse instance has no attribute 'EASummary'

In [ ]: