In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import aplpy
import pylab
from astropy.io import fits, ascii
from astropy.coordinates import ICRS, angles, SkyCoord
import astropy.units as u
from astropy.wcs import wcs
import ROOT
from ROOT import TRolke, TFile, gROOT, gSystem, TGraph, TFeldmanCousins
import rootnotes
import rootprint
import scipy.ndimage as ndimage
import gammapy.spectrum
import os
from gammapy import stats
# import fit
pylab.rcParams['figure.figsize'] = (10.0, 8.0)
plt.rc('font', family='serif', size=14)
In [2]:
%%rootprint
rootDir=os.path.expandvars("$VEGAS")
ROOT.gROOT.Reset()
gSystem.Load("libTreePlayer.so")
gSystem.Load("libPhysics.so")
gSystem.Load(rootDir + "/common/lib/libSP24sharedLite.so")
gSystem.Load(rootDir + "/resultsExtractor/lib/libStage6shared.so")
gROOT.ProcessLine(".L " + rootDir + "/resultsExtractor/include/VAEffectiveAreaCommon.h")
gROOT.ProcessLine(".L " + rootDir + "/common/include/VACommon.h")
gROOT.ProcessLine(".include " + rootDir + "/common/include/")
gROOT.ProcessLine(".include " + rootDir + "/resultsExtractor/src/")
gROOT.ProcessLine(".include " + rootDir + "/resultsExtractor/include/")
gROOT.ProcessLine(".include " + rootDir + "/cfitsio/include/")
In [3]:
cuts = 'Hard'
theta = 'Point'
if cuts == 'Hard':
clean = 'ED'
else:
clean = 'Thresh'
fitsF = "rootdir/St6/All" + cuts + theta + clean + "M31.fits"
save=False
plot=False
livetime = 54.69*60.*60. #seconds
#UL
index = -2.5
diffUL = 3e-9 # TeV-1 m-2 s-1 @ 1TeV (3e-9 = 1% Crab)
#PSF
sigma=0.065
rolkeUL=0.95
if theta == 'Point':
nTestReg = 8 #Number of test regions
sepDist = 0.1 #Test region radius
ULpos = ascii.read("PointLocations.txt")
else:
nTestReg = 2
sepDist = 0.2
ULpos = ascii.read("ExtLocations.txt")
#Plotting options
pngRes=200 #dpi
#cx1 = cubehelix.cmap(reverse=False, start=1.0, rot=-1)
cx1 = 'Blues'
cx1 = 'default'
In [4]:
def standard_setup(sp, exc=True):
###
# Plot formatting
###
sp.ticks.show()
sp.ticks.set_color('black')
sp.tick_labels.set_xformat('dd')
sp.tick_labels.set_yformat('dd')
sp.ticks.set_xspacing(2) # degrees
sp.ticks.set_yspacing(2) # degrees
sp.set_frame_color('black')
sp.set_tick_labels_font(size='10')
sp.set_axis_labels_font(size='12')
###
# IRIS Contours
###
contours = [1500, 2500, 3500 ]
sp.show_contour('M31_IRIS_cropped_ds9.fits',
levels=contours, smooth=1,
linewidths=1., colors="k", zorder=2) #cmap=mpl.cm.gray
###
# Exclusion Regions
###
excColor="black"
excLW=1.
#sp.show_circles(10.6847, 41.2687, 0.3, color=excColor, linewidth=excLW, zorder=4) # M31 center
#sp.show_circles(11, 41.5, 0.3, color=excColor, linewidth=excLW, zorder=4)
#sp.show_circles(11.2, 41.75, 0.3, color=excColor, linewidth=excLW, zorder=4)
#sp.show_circles(10.4, 41, 0.3, color=excColor, linewidth=excLW, zorder=4)
#sp.show_circles(10.25, 40.7, 0.25, color=excColor, linewidth=excLW, zorder=4)
sp.show_circles(12.4535, 41.079, 0.3, color=excColor, linewidth=excLW, zorder=4) # nu Andromedae
if exc:
sp.add_label(12.4535, 41.45, r'$\nu$' + '-Andromedae',
size=12, weight='demi', color='black')
###
# Colorbar
###
sp.add_colorbar()
sp.colorbar.set_location('right')
sp.colorbar.set_width(0.2) # arbitrary units, default is 0.2
sp.colorbar.set_pad(0.05) # arbitrary units, default is 0.05
sp.colorbar.set_font(size=10)#, weight='bold', variant='normal')
def sumInRegion(data, header, raIn, decIn, radius):
'''
This is the total number of counts from the fits file (passed as data and header)
within the region defined by raIn, decIn and its radius.
It returns the total counts, that is all'''
wcs_transformation = wcs.WCS(header)
#init = wcs_transformation.wcs_world2pix(10.6847, 41.2687, 0) # center of fit pt, this case Crab Nebula
target = SkyCoord(raIn, decIn, unit='deg', frame='icrs')
y, x = np.mgrid[:data.shape[0], :data.shape[1]]
ra, dec = wcs_transformation.wcs_pix2world(x, y, 0)
Point = SkyCoord(ra, dec, unit='deg', frame='icrs')
sep = target.separation(Point)
data2 = data[(sep.deg<radius)]
total = np.sum(data2)
return total
def sumInRegionPoint(data, xIn, yIn, radius):
'''
This is the total number of counts
within the region defined by xIn, yIn and its radius.
It returns the total counts, that is all'''
y, x = np.mgrid[:data.shape[0], :data.shape[1]]
sep = np.sqrt((x-xIn)**2 + (y-yIn)**2)
data2 = data[(sep<radius)]
total = np.sum(data2)
return total
def pointInEllipse(x,y,xp,yp,d,D,angle):
#tests if a point[xp,yp] is within
#boundaries defined by the ellipse
#of center[x,y], diameter d D, and tilted at angle
angle = np.deg2rad(angle)
cosa=np.cos(angle)
sina=np.sin(angle)
dd=d/2*d/2
DD=D/2*D/2
a =np.power(cosa*(xp-x)+sina*(yp-y),2)
b =np.power(sina*(xp-x)-cosa*(yp-y),2)
ellipse=(((a/dd)+(b/DD))<1)
return ellipse
In [5]:
if plot:
fig = plt.figure(figsize=(12, 12))
fig1 = aplpy.FITSFigure(fitsF, figure=fig, subplot=(2,3,1), hdu=5)
fig1.show_colorscale()
standard_setup(fig1)
fig1.set_title("Acceptance")
fig1 = aplpy.FITSFigure(fitsF, figure=fig, subplot=(2,3,2), hdu=7)
fig1.show_colorscale()
standard_setup(fig1)
fig1.set_title("Acceptance")
fig1 = aplpy.FITSFigure(fitsF, figure=fig, subplot=(2,3,3), hdu=9)
fig1.show_colorscale()
standard_setup(fig1)
fig1.set_title("Acceptance")
fig1 = aplpy.FITSFigure(fitsF, figure=fig, subplot=(2,3,4), hdu=11)
fig1.show_colorscale()
standard_setup(fig1)
fig1.set_title("Acceptance")
fig1 = aplpy.FITSFigure(fitsF, figure=fig, subplot=(2,3,5), hdu=13)
fig1.show_colorscale()
standard_setup(fig1)
fig1.set_title("Acceptance")
fig1 = aplpy.FITSFigure(fitsF, figure=fig, subplot=(2,3,6), hdu=15)
fig1.show_colorscale()
standard_setup(fig1)
fig1.set_title("Acceptance")
This is done by reading in the on counts and integral acceptance maps and using them to determin the total counts and acceptance within each of the test regions. Since we are dealing with circular regions we can simply use the standard caclulation for separation from astropy
In [6]:
onCounts = 0
onAcc = 0
nbinsOn = 0
onC = []
onA = []
totCounts = 0
grAcc = []
ptAcc = np.zeros([6, nTestReg])
for group in range(6):
#counts setup
extname = 'RawOnMap'+str(group)
onData, onHeader = fits.getdata(fitsF, header=True, extname=extname)
wcs_transformation_on = wcs.WCS(onHeader)
yOn, xOn = np.mgrid[:onData.shape[0], :onData.shape[1]]
raOn, decOn = wcs_transformation_on.all_pix2world(xOn, yOn, 0)
onPos = SkyCoord(raOn, decOn, unit='deg', frame='icrs')
totCounts += np.nansum(onData)
#acceptance setup
accData = fits.getdata(fitsF, header=False, extname='AcceptanceMap'+str(group))
#bin area reading (getting this from the root file as easier)
fName = "rootdir/St6/All" + cuts + theta + clean + "s6.root"
f = TFile(fName, "read")
RBM = f.Get("RingBackgroundModelAnalysis/SkyMapOn")
for i in range(accData.shape[0]):
for j in range(accData.shape[1]):
accData[i][j] = accData[i][j] * RBM.GetBinArea(i,j)
gAcc = 0
for i in range(nTestReg):
Pt = SkyCoord(ULpos['col1'][i], ULpos['col2'][i], unit='deg', frame='icrs')
onSep = Pt.separation(onPos)
cnts = np.nansum(onData[onSep.deg<sepDist])
accSep = Pt.separation(onPos)
acc = np.nansum(accData[onSep.deg<sepDist])*np.nansum(onData)
gAcc += acc
ptAcc[group, i] = acc
if group == 0:
onC.append(cnts)
onA.append(acc)
else:
onC[i] += cnts
onA[i] += acc
grAcc.append(gAcc)
#print i+1, ULpos['col1'][i], ULpos['col2'][i], counts, acc
print onC
print onA
onCounts = np.sum(onC)
onAcc = np.sum(onA)
print "Total", onCounts, onAcc, totCounts
print grAcc
In [7]:
onCountsE = 0
onAccE = 0
nbinsOnE = 0
onCE = 0
onAE = 0
totCountsE = 0
grAccE = 0
inElcent1E = SkyCoord(11.5, 42.00, unit='deg', frame='icrs')
inElcent2E = SkyCoord(10.0, 40.55, unit='deg', frame='icrs')
inElDistE = 1.8
for group in range(6):
#counts setup
extname = 'RawOnMap' + str(group)
onData, onHeader = fits.getdata(fitsF, header=True, extname=extname)
wcs_transformation_on = wcs.WCS(onHeader)
yOn, xOn = np.mgrid[:onData.shape[0], :onData.shape[1]]
raOn, decOn = wcs_transformation_on.all_pix2world(xOn, yOn, 0)
onPos = SkyCoord(raOn, decOn, unit='deg', frame='icrs')
totCountsE += np.nansum(onData)
#acceptance setup
accData = fits.getdata(fitsF, header=False, extname='AcceptanceMap'+str(group))
#bin area reading (getting this from the root file as easier)
fName = "rootdir/St6/All" + cuts + theta + clean + "s6.root"
f = TFile(fName, "read")
RBM = f.Get("RingBackgroundModelAnalysis/SkyMapOn")
for i in range(accData.shape[0]):
for j in range(accData.shape[1]):
accData[i][j] = accData[i][j] * RBM.GetBinArea(i,j)
gAccE = 0
onSep = ((inElcent1E.separation(onPos).deg + inElcent2E.separation(onPos).deg) > inElDistE)
onCE += np.nansum(onData[onSep < inElDistE])
accSep = Pt.separation(onPos)
onAE += np.nansum(accData[onSep < inElDistE])*np.nansum(onData)
gAccE += acc
print onCE
print onAE
print grAcc
The background counts are harder, we are trying to use an ellipse in camera coordinates (that is a flat coord scheme) whereas the data is saved in a spherical coord scheme, to overcome this we cheat a bit by defining the eclipses and then saving them as fits files with a header copied from the data, this means that the projections etc are correct. We can now check this ellipse is okay by drawing it (as a contour) over the skymap and checking that the correct regions are included/excluded. The background will be integrated between the two ellipses (less the bit cut out for nuAndromadae)
In [8]:
inElcent1 = SkyCoord(11.5, 42.00, unit='deg', frame='icrs')
inElcent2 = SkyCoord(10.0, 40.55, unit='deg', frame='icrs')
outElcent1 = inElcent1
outElcent2 = inElcent2
inElDist = 2.2
outElDist = 3.9
inEl = ((inElcent1.separation(onPos).deg + inElcent2.separation(onPos).deg) > inElDist)
outEl = ((outElcent1.separation(onPos).deg + outElcent2.separation(onPos).deg) < outElDist)
El = np.logical_and(inEl, outEl)
nuAndrom = SkyCoord(12.4535, 41.0790, unit='deg', frame='icrs')
sepnuAn = (nuAndrom.separation(onPos).deg > 0.4)
El = np.logical_and(El, sepnuAn)
! rm BGReg.fits
hdu = fits.PrimaryHDU(El*1., header=onHeader)
hdu.writeto('BGReg.fits')
In [9]:
if plot:
fig = plt.figure(figsize=(8, 8))
fig1 = aplpy.FITSFigure(fitsF, figure=fig, hdu=1)
fig1.show_colorscale(vmin=-5,vmax=5,cmap=cx1)
standard_setup(fig1)
fig1.set_title("Significance")
fig1.show_contour("BGReg.fits", colors='k')
for i in range(nTestReg):
fig1.show_circles(ULpos['col1'][i], ULpos['col2'][i], sepDist,
color='purple', linewidth=2, zorder=5)
fig1.add_label(ULpos['col1'][i], ULpos['col2'][i], ULpos['col3'][i],
size=16, weight='bold', color='purple')
In [10]:
if plot:
fig = plt.figure(figsize=(10, 10))
fig1 = aplpy.FITSFigure("M31_IRIS_smoothed.fits", figure=fig)
fig1.show_colorscale(cmap='Blues',vmin=0, vmax=7e3)
fig1.recenter(10.6847, 41.2687, width=4, height=4)
fig1.ticks.show()
fig1.ticks.set_color('black')
fig1.tick_labels.set_xformat('dd.dd')
fig1.tick_labels.set_yformat('dd.dd')
fig1.ticks.set_xspacing(1) # degrees
fig1.set_frame_color('black')
fig1.set_tick_labels_font(size='14')
fig1.set_axis_labels_font(size='16')
fig1.show_grid()
fig1.set_grid_color('k')
fig1.add_label(12.4535, 41.09, r'$\nu$' + '-Andromedae', size=10, weight='demi', color='black')
#fig1.show_contour("BGReg.fits", colors='k')
fig1.show_contour("BGReg.fits", lw=0.5, filled=True, hatches=[None,'/'],
colors='none')
fig1.show_contour("BGReg.fits", linewidths=1., filled=False,
colors='k', levels=1)
plt.savefig("Plots/M31BgReg.pdf")
In [11]:
bgC = []
bgA = []
ptAlpha = np.empty([6, nTestReg])
for group in range(6):
#counts setup
extname = 'RawOnMap'+str(group)
onData, onHeader = fits.getdata(fitsF, header=True, extname=extname)
wcs_transformation_on = wcs.WCS(onHeader)
yOn, xOn = np.mgrid[:onData.shape[0], :onData.shape[1]]
raOn, decOn = wcs_transformation_on.all_pix2world(xOn, yOn, 0)
onPos = SkyCoord(raOn, decOn, unit='deg', frame='icrs')
#acceptance setup
accData = fits.getdata(fitsF, header=False, extname='AcceptanceMap'+str(group))
#bin area reading (getting this from the root file as easier)
fName = "rootdir/St6/All" + cuts + theta + clean + "s6.root"
f = TFile(fName, "read")
RBM = f.Get("RingBackgroundModelAnalysis/SkyMapOn")
for i in range(accData.shape[0]):
for j in range(accData.shape[1]):
accData[i][j] = accData[i][j] * RBM.GetBinArea(i,j)
bgC.append(np.nansum(onData[El]))
bgA.append(np.nansum(accData[El]))#*np.nansum(onData))
ptAlpha[group,] = ptAcc[group,] / np.nansum(accData[El])/ totCounts
#ptAlpha = ptAcc[group, :] / np.nansum(accData[El])/ totCounts
grAlpha = np.array(grAcc) / np.array(bgA)/ totCounts
ptAlpha = np.sum(ptAlpha, axis=0)
alpha = np.sum(grAlpha)
bgCounts = np.sum(bgC)
excess = onCounts - bgCounts * alpha
print "Total:", onCounts, bgCounts, excess, alpha
print "Point Alpha:", ptAlpha
In [19]:
stats.significance_on_off(onCE, bgCounts, )
In [12]:
if plot:
bins = np.linspace(-4.5, 4.5, 100)
sigData, sigHeader = fits.getdata(fitsF, header=True, extname="SignificanceMap")
fig = plt.figure(figsize=(3, 3))
fig, ax = plt.subplots(1)
hist = plt.hist(sigData[(~np.isnan(sigData)) & El], bins=bins, histtype="step")
plt.semilogy()
hist, bins2 = np.histogram(sigData[(~np.isnan(sigData)) & El], bins = bins)
(xf, yf), params, err, chi = fit.fit(fit.gaus, (bins2[0:-1] + bins2[1:])/2, hist)
plt.plot(xf, yf, 'r-', label='Fit')
textstr1 = '$\mu = %.2f $' % (params[1])
textstr2 = '$ %.3f$\n$\sigma = %.2f$' % (err[1], params[2])
textstr3 = '$ %.3f$' % (err[2])
textstr = textstr1 + u"\u00B1" + textstr2 + u"\u00B1" + textstr3
#textstr = textstr1 + textstr2 + textstr3
props = dict(boxstyle='square', alpha=0.5, fc="white")
ax.text(0.95, 0.95, textstr, transform=ax.transAxes, fontsize=14,
verticalalignment='top', horizontalalignment='right',
bbox=props)
plt.ylim(ymin=1e0)
In [13]:
def RUL(on, off, alpha):
rolke = TRolke(rolkeUL)
rolke.SetBounding(True)
rolke.SetPoissonBkgKnownEff(int(on), int(off), 1./(alpha), 1.)
return rolke.GetUpperLimit()
def FCUL(on, off, alpha):
fc = TFeldmanCousins(rolkeUL)
return fc.CalculateUpperLimit((on), (off) * alpha)
In [14]:
if False:
for i in range(nTestReg):
ULCounts = RUL((onC[i]), (bgCounts), ptAlpha[i])
excess = onC[i] - bgCounts * ptAlpha[i]
print "Point", i
print 'On = {0:.0f}, Off = {1:.0f}, alpha = {2:.4f}'.format(onC[i], bgCounts, ptAlpha[i])
print 'Excess = {0:.2f}'.format(excess)
print 'Signif = {0:.3f}'.format(gammapy.stats.significance_on_off(onC[i], bgCounts, ptAlpha[i]))
print 'ULCount = {0:0.3f}'.format(ULCounts)
print ''
excess = onCounts - bgCounts * alpha
ULCounts = RUL((onCounts), (bgCounts), (alpha))
print 'On = {0:.0f}, Off = {1:.0f}, alpha = {2:.4f}'.format(onCounts, bgCounts, alpha)
print 'Excess = {0:.2f}'.format(excess)
print 'Signif = {0:.3f}'.format(gammapy.stats.significance_on_off(onCounts, bgCounts, alpha))
print 'ULCount = {0:0.3f}'.format(ULCounts)
print ''
Needto sum all of the Effective Areas from each of the test positions and then work out the expected flux given a test spectrum
What we need to correct each EA by the difference in the flux distribtuion in its test region. To do this we use the following relation:
$\frac{Frac\; Region\; Flux\; in\; thetaSq}{Frac\; Point\; Source\; Flux\; in\; thetaSq}$
For the bottom bit I take the point source and convolve with PSF, work out the fraction of counts that remain within the thetaSq
For the top bit, I take the model, convolve with the PSF and work out the fraction of counts before to after
Logic: think, not all counts fall within thetaSq, thus the EA is slightly under estimate, since Flux * EA = counts. Thus we need to undo this for the point source and redo this for the extended source. The smoothing factor is put in such that 68\% of the flux falls within a 0.1deg region (this is the standard quoted number) for a point source. I would like to check this with hard cuts etc for sims, but that should be a secondary effect
In [15]:
pointData = np.copy(onData)
pointData.fill(0)
pointData[pointData.shape[0]/2., pointData.shape[1]/2.] = 1000
pointData1 = ndimage.gaussian_filter(pointData, sigma=(-sigma/onHeader['CDELT1'], sigma/onHeader['CDELT2']),
order=0)
wcs_transformation = wcs.WCS(onHeader)
initPos = wcs_transformation.wcs_pix2world(pointData.shape[0]/2., pointData.shape[1]/2., 0)
pointSourceCor = sumInRegion(pointData1, onHeader, initPos[0], initPos[1], sepDist)/np.sum(pointData1)
In [16]:
IRISdata, IRISheader = fits.getdata("M31_IRIS_cropped_ds9.fits", header=True)
IRISdata2 = ndimage.gaussian_filter(IRISdata,
sigma=(-sigma/IRISheader['CDELT1'], sigma/IRISheader['CDELT2']),
order=0)
M31total = np.sum(IRISdata2)
In [17]:
%%rootprint
nPts = 100
En = np.linspace(-1, 2, num=nPts)
Sp1 = (10**En)**index
EA = np.empty([nPts])
EA1 = np.empty([nPts])
minSafeE = 0 #this is the minimum safe energy, I will quote the spectrum here
decorE = 0
EstCounts1 = 0
for j in range(nTestReg):
fName = "rootdir/St6/All" + cuts + theta + clean + str(j+1) + "s6.root"
f = TFile(fName, "read")
UL = f.Get("UpperLimit/VAUpperLimit")
g = UL.GetEffectiveArea()
if UL.GetEnergy() > minSafeE:
minSafeE = UL.GetEnergy()
if UL.GetEdecorr() > decorE:
decorE = UL.GetEdecorr()
# Weight EAs by expected flux from that region
irisReg = sumInRegion(IRISdata, IRISheader, ULpos['col1'][j], ULpos['col2'][j], sepDist)
regW = irisReg / M31total
# Correct for PSF effects
M31RegCor = sumInRegion(IRISdata2, IRISheader, ULpos['col1'][j], ULpos['col2'][j], sepDist) / irisReg
print j, regW, M31RegCor
for i, xval in np.ndenumerate(En):
EA[i] = g.Eval(xval) / pointSourceCor * M31RegCor * regW
EA1[i] += g.Eval(xval) / pointSourceCor * M31RegCor * regW
Fl1 = Sp1 * EA
EstCounts1 += np.trapz(Fl1, 10**En)*livetime
FluxULReg1 = ULCounts / EstCounts1
print FluxULReg1
In [18]:
FluxULM31 = FluxULReg1
FluxULM31_eMin = FluxULM31 * minSafeE **index
FluxULM31_eDec = FluxULM31 * decorE **index
intULM31_eMin = gammapy.spectrum.powerlaw.power_law_integral_flux(FluxULM31, index, 1, minSafeE, 30)
intULM31_eDec = gammapy.spectrum.powerlaw.power_law_integral_flux(FluxULM31, index, 1, decorE, 30)
intULM31_eMin_pcCrab = intULM31_eMin /(gammapy.spectrum.crab_integral_flux(minSafeE, 30, 'hess_pl')[0] *1e2)
intULM31_eDec_pcCrab = intULM31_eDec /(gammapy.spectrum.crab_integral_flux(decorE, 30, 'hess_pl')[0] *1e2)
print 'On = {0:.0f}, Off = {1:.0f}, alpha = {2:.4f}'.format(onCounts, bgCounts, alpha)
print 'Excess = {0:.2f}'.format(excess)
print 'Signif = {0:.3f}'.format(gammapy.stats.significance_on_off(onCounts, bgCounts, alpha))
print 'ULCount = {0:0.3f}'.format(ULCounts)
print ''
print "Differential UL @ 1TeV = {0:.3e}".format(FluxULM31)
print "Differential UL @ min safe E ({0:.1f}GeV) = {1:.3e}".format(minSafeE*1e3, FluxULM31_eMin)
print "Differential UL @ decorrelation ({0:.1f}GeV) = {1:.3e}".format(decorE*1e3, FluxULM31_eDec)
print "Differential UL units = TeV-1 m-2 s-1"
print ""
print "Integral UL between min safe energy and 30TeV = {0:.3e}".format(intULM31_eMin)
print "Integral UL between decorrel energy and 30TeV = {0:.3e}".format(intULM31_eDec)
print "Integral UL units = m-2 s-1"
print ""
print "Integral UL between min safe energy and 30TeV = {0:.3f} %Crab".format(intULM31_eMin_pcCrab)
print "Integral UL between decorrel energy and 30TeV = {0:.3f} %Crab".format(intULM31_eDec_pcCrab)
In [265]:
print 244. / 5.4, 244. / 0.67, 244. / 17.3
print 382. / 7.9, 382. / 3.90, 382. / 25.1
print 65. / 2.4, 65. / 0.30, 65. / 7.75
print 137. / 6.0, 137. / 3.00, 137. / 19.1
In [ ]:
In [ ]: