About the fast analysis test page

In order to make the first step towards a automated fast supernova signal analysis (based on hitspool data), a simulated signal was combined with background gathered from hitspool data. The goal is to provide meaningful statemeents about the neutrino lightcurve, the signal onset time and supernova internal parameters like accretion time constant to the scientific comunity.

Setting up the environment

In [16]:
import getopt, sys
import datetime
import glob
import re
import json
import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
In [20]:
%matplotlib inline
In [3]:
sns.set_palette(palette='pastel')
sns.set_style(style='whitegrid')
sns.set_context('poster')
In [2]:
import snfluxnew

The above imported module snfluxnew contains several implementation of the Pagliaroli flux parameterization on which all fits to the lightcurve are based on. to be found in the HsAnaSuite SVN project by David Heereman in his sandbox

Import hitspool background data and simulated (Lawrence-Livermore) signal

In [24]:
hsfile = "/data/ana/Supernova/Hitspool/SNALERT_20150221_112506/HitTables/SNALERT_20150221_112506_hitspool_subthreshold_1msSCALERS.txt"
# sgfile = "$SVN/snfit-mcmc/Livermore_signal_10kpc.txt"
sgfile = "/home/hitspool/HsAnaSuite/Livermore_signal_10kpc.txt"
xdata, ydata = snfluxnew.getdata(hsfile, sgfile)

Plotting the hit rates

Basic Hit rate <-> luminosity relation:

$$R(t) = \epsilon_\mathrm{dead time}\frac{n_\mathrm{target} \,L_\mathrm{SN}^\nu(t)}{4\pi d^2\overline{E_{\nu}}(t)} \int_0^{\infty} \,dE_\mathrm{e} \!\! \int_0^{\infty} \,dE_{\nu} \times\, \frac{d\sigma}{dE_\mathrm{e}}(E_\mathrm{e},E_{\nu}) \, \,V_\mathrm{e^+}^\mathrm{eff} \,f(E_{\nu},\overline{E_{\nu}}, \alpha_\nu,t)$$

Approximate expression for $\frac{d\sigma}{dE_e}(E_e,E_v)$:

$$\frac{d\sigma}{dE_e}\approx 10^{-43}\mathrm{cm}^2\cdot p_e\cdot E_e\cdot E_\nu^{-0.07056+0.02018 \ln E_\nu -0.001953 (ln E_\nu)^3}\cdot \delta(E_e-E_\nu+\Delta) =\sigma(E_\nu)\cdot \delta(E_e-E_\nu+\Delta)$$

Note that this is only valid in the coarse approximation that $E_e=E_\nu-\Delta$ with $\Delta=m_n-m_p\approx 1.293$ MeV; $p_e=\sqrt{E_e^2-m_e^2}$.

$V_\mathrm{e^+}^\mathrm{eff} =29.0\pm 3.8$ m$^3$ $\cdot E_e$ and $\epsilon_\mathrm{dead time}\approx 0.87/(1+R(t)\cdot \tau)$, with artificial deadtime $\tau=250\times 10^{-6}$ s and target density $n_\mathrm{target}=2\cdot (916.6\pm 1.6)$kg/m$^3 \cdot N_A /(0.01801528$ kg/mol) ; $N_A=6.022140857\times 10^{23}$.

$$R(t)= \epsilon_\mathrm{dead time}\frac{n_\mathrm{target} \,L_\mathrm{SN}^\nu(t)}{4\pi d^2\overline{E_{\nu}}(t)} \int_0^{\infty} \,dE_{\nu} \times\, \sigma(E_{\nu}) \, \,V_\mathrm{e^+}^\mathrm{eff} \,f(E_{\nu},\overline{E_{\nu}}, \alpha_\nu,t)$$

$$R(t)= 29.0\cdot 10^{-43}\epsilon_\mathrm{dead time}\frac{n_\mathrm{target} \,L_\mathrm{SN}^\nu(t)}{4\pi d^2\overline{E_{\nu}}(t)} \int_0^{\infty} \,dE_{\nu} \times\, \sqrt{(E_\nu-\Delta)^2-m_e^2}\cdot E_\nu^{-0.07056+0.02018 \ln E_\nu -0.001953 (ln E_\nu)^3}\cdot (E_\nu-\Delta)^2\,f(E_{\nu},\overline{E_{\nu}}, \alpha_\nu,t)$$

or in the more naive approximation $\sigma \approx 9.52\times 10^{-44}\frac{p_e E_e}{\mathrm {MeV}^2}$ cm$^2$ $$R(t)= 29.0\cdot 9.52\cdot 10^{-44}\epsilon_\mathrm{dead time}\frac{n_\mathrm{target} \,L_\mathrm{SN}^\nu(t)}{4\pi d^2\overline{E_{\nu}}(t)} \int_0^{\infty} \,dE_{\nu} \times\, \sqrt{(E_\nu-\Delta)^2-m_e^2}\cdot (E_\nu-\Delta)^2\,f(E_{\nu},\overline{E_{\nu}}, \alpha_\nu,t)$$

for $E_e\approx E_\nu\approx p_\nu$:

$$R(t) \approx 0.276\cdot 10^{-41}\epsilon_\mathrm{dead time}\frac{n_\mathrm{target} \,L_\mathrm{SN}^\nu(t)}{4\pi d^2\overline{E_{\nu}}(t)} \int_{0}^{\infty} \,dE_{\nu} \times\, E_{\nu}^3\,f(E_{\nu},\overline{E_{\nu}}, \alpha_\nu,t)\\ \approx 0.276 \cdot 10^{-41}\epsilon_\mathrm{dead time}\frac{n_\mathrm{target} \,L_\mathrm{SN}^\nu(t)}{4\pi d^2} \overline{E_{\nu}^3}/\overline{E_{\nu}}(t)$$

solving for the luminosity $L_\mathrm{SN}^\nu(t)$: $$L_\mathrm{SN}^\nu(t)\approx 4.165\times 10^{41}\frac{R(t)\cdot (1+R(t)\cdot 0.0025)\cdot 4\pi \cdot d^2 \cdot \overline{E_{\nu}}}{n_\mathrm{target}\cdot \overline{E_{\nu}^3}(t)} \approx 89676135863999.17\frac{R(t)\cdot (1+R(t)\cdot 0.0025)\cdot d^2\cdot \overline{E_{\nu}}}{\overline{E_{\nu}^3}(t)} $$

Which leads to our final expression for the relation ofthe hit rate and the luminosity:

$$ R(t) \cdot (1+R(t)\cdot 0.00025 )= \frac{L_\mathrm{SN}^\nu(t)}{d^2} \frac{\overline{E_{\nu}^3}(t)}{\overline{E_{\nu}}} \cdot 1.11 \times 10^{-14}$$

import math NA=6.022140857E23 rho=916.6 mmol=0.01891528 me=0.511 Delta=1.293 ntarget=rho/mmolNA2 print "ntarget = ", ntarget factor = 4.165E414math.pi/ntarget print "factor = ", factor def dt(x): return 0.87(1+x0.00025) def Pe(E): return math.sqrt((E-Delta)2-me2) def Ee(E): return E-Delta def Veff(E): return 29.0*(E-Delta)

Following the considerations by Lutz for the relation of the detector rate to the physical supernova luminosity, we plot $R(t) \cdot (1 + 0.00025 R(t))$ versus time

In [21]:
ydata.mean()
Out[21]:
1421.7305555555556
In [25]:
ydatameanaccess = ydata - ydata.mean()

rateexpression = ydatameanaccess * (1 + ydatameanaccess*0.00025)
In [27]:
plt.figure(figsize(14,10))
Out[27]:
<matplotlib.figure.Figure at 0x59e9750>
<matplotlib.figure.Figure at 0x59e9750>
In [28]:
fig,ax = plt.subplots()
plt.step(xdata[1:], rateexpression[1:], color='#3B9C9C', alpha=0.8, where='post')
ax.text(0.01, 1.01, '+ %i' % int(ydata.mean()), transform=ax.transAxes, fontsize=14)
ax.text(0.7, .7, 'TEST PAGE\nSimulated data*', color='red', transform=ax.transAxes, fontsize=24)
ax.text(0.7, .55, '___________________\n*simulation:\n signal free hitspool data\n 10 kpc distant \n Lawrence-Livermore signal', 
        color='red', transform=ax.transAxes, fontsize=14)

plt.suptitle('SN candidate from TEST \n 0.001 s binning')

ax.set_xlabel('Time wrt trigger time [s]')
ax.set_ylabel(r'$R(t) \cdot (1+R(t)\cdot 0.00025 )= \frac{L_\mathrm{SN}^\nu(t)}{d^2} \frac{\overline{E_{\nu}^3}(t)}{\overline{E_{\nu}}} \cdot 1.11 \times 10^{-14}$', fontsize = 20)
plt.show()

Plotting hit rates (i.e. luminosity) with different bin sizes

In [11]:
binsizes = [0.001, 0.002, 0.01, 0.02, 0.05, 0.1, 0.2, 
                0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 
                1.0, 1.5 ,2.0 , 2.5 ,3.0 ,3.5 ,4.0,
                4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 
                8.0, 8.5, 9.0, 9.5, 10]
In [15]:
plt.figure(figsize(14,10))
for binsize in binsizes:
    binfactor = int((binsize * 1e3) / 1) # binsize in multiples     of the basic bin size of 2ms
    binavg = ydata.mean() * binfactor
    ydatarebin = np.array([sum(rateexpression[i:i+binfactor]) for i in range(len(rateexpression))[::binfactor]])
    xdata = np.arange(-30,60,binsize)
    fig,ax = plt.subplots()
    plt.suptitle('SN candidate from TEST \n %s s binning' %binsize)
    plt.step(xdata, ydatarebin, color='#3B9C9C', where='post')
    ax.set_xlabel('Time wrt trigger time [s]')
    ax.set_ylabel(r'$R(t) \cdot (1+R(t)\cdot 0.00025 )= \frac{L_\mathrm{SN}^\nu(t)}{d^2} \frac{\overline{E_{\nu}^3}(t)}{\overline{E_{\nu}}} \cdot 1.11 \times 10^{-14}$', fontsize = 20)    
    ax.text(0.01, 1.01, '+ %i' % int(binavg), transform=ax.transAxes, fontsize=12)
    ax.text(0.7, .7, 'TEST PAGE\nSimulated data*', color='red', transform=ax.transAxes, fontsize=24)
    ax.text(0.7, .55, '___________________\n*simulation:\n signal free hitspool data\n 10 kpc distant \n Lawrence-Livermore signal', color='red', transform=ax.transAxes, fontsize=14)

    #plotnameidx = '_'+str(int(binsize*1000)).zfill(6)+'ms_bin.png'
    #plotname = re.sub('.txt', plotnameidx, infile)
    fig.savefig('/home/hitspool/HsAnaSuite/plots/ratesplottest_'+str(binfactor).zfill(5)+'msbins.png')
    plt.show()
    plt.close(fig)

<matplotlib.figure.Figure at 0x5bcac90>

Zoom in and fits

Signal onset time and rise-time $\tau_r$

In opposition to the highly complicated, at least nine parameter, fit of Pagliaroli et al as described in her papers (here and here), we separate the signal onset fit since we only areintereseted in the correct signal onset time and the rise-time in that part and do not need the rest.

In [29]:
from scipy.optimize import curve_fit
In [30]:
ydata = ydata -ydata.mean()
In [31]:
fitparas, fitcov = curve_fit(snfluxnew.fit_onset, xdata[29900:30200], ydata[29900:30200], p0=[0.006, 0.05,600])
In [32]:
myfit = snfluxnew.fit_onset(xdata, *fitparas)
In [33]:
def statsbox(fitparas, fitcov, t0=30000):
    '''
    Build the string that print the stats.
    '''
    fiterr = np.sqrt(np.diag(fitcov))
    
    string1 = '\n$t_a$ : %.4f $\pm$ %.4f s\n$\\tau_a$ : %.2f $\pm$ %.2f s\n$R_{max}$ : %.2d $\pm$ %.2d' %(fitparas[0], 
                                      fiterr[0],
                                      fitparas[1],
                                      fiterr[1],
                                      fitparas[2],
                                      fiterr[2])
    string2 = 'Onset time: t0 {:+} bins\n'.format((np.argmax(myfit>0.0)-t0))
    string3 = '$t_{0}$ [UTC] : %s\n$t_r$: $t_{0}$ + %.6f $\pm$ %.6f s\n$\\tau_r$ : %.2f $\pm$ %.2f s' % ('2016-02-02 02:22:22.123456', 
                                                                                                         fitparas[0], fiterr[0],
                                                                                                        fitparas[1],
                                                                                                        fiterr[1])
    #return string2 + string1 + string3
    return string3
In [34]:
fig, ax = plt.subplots()
plt.step(xdata[1:], ydata[1:], lw=3)

plt.plot(xdata[1:30150], myfit[1:30150], lw=4, color='red')

props = dict(boxstyle='round',facecolor='white', alpha=0.8)
plt.title('TEST SN \nat 2016-02-02 02:22:22.123456')
# place a text box in upper left in axes coords
ax.text(0.2, 0.25, statsbox(fitparas, fitcov), transform=ax.transAxes, fontsize=24,
        verticalalignment='top', bbox=props)
ax.set_xlabel('Time post bounce [s]')
ax.set_ylabel('Events per bin [1.0 ms]')
plt.xlim(-0.01, 0.1)
#fig.savefig('/home/david/workspace/dheereman/HsAnaHTML/hs_sn/ratesplottest_onsetfit.png')
plt.show()

Fitting the rest

The remaining part of the original Pagliaroli model is now fitted. in order to have converging fit, it was found that holding a couple of parameters fixed at reasonable values was the only way to success. All lot of different techniques for the fit (eg Markov Chain Monte Carlos) were tested but didnt not reveal an inprovement to this method in terms of speed.

In [35]:
fig, ax = plt.subplots()
plt.step(xdata[1:], ydata[1:], color='#3B9C9C', alpha=0.8, where='post')
plt.xlim(-10,20)
Out[35]:
(-10, 20)
In [37]:
import snflux

The above imported module snfluxnew contains several implementation of the Pagliaroli flux parameterization on which all fits to the lightcurve are based on. to be found in the HsAnaSuite SVN project by David Heereman in his sandbox

In [38]:
cppmodel = snflux.FluxPara()

cppfluxmodel = cppmodel.fluxallparams(xdata, 12., 12., 0.1, 0.33, 2.17, .66, 17.4, 4.47, 4.63)
In [39]:
# cpp coded model with fixed parameters par0, par3, par4, par6 used to fit the data:
popts1, perr1 = curve_fit(cppmodel.fluxfixedparams0346, xdata, ydata, p0=[ 9.83e-3, 3.02e+3, 2.74e-1, 1.06, 5.56])
In [40]:
def cppstatsbox(params, pcov):
    perrs = np.sqrt(np.diag(pcov))
    par0=12.
    par3=0.2
    par4=2.5
    par6=16.
    paramstring = r'''
    $E_{\nu}$ = %.2f  
    $M_{a}$ = %.1f
    $T_{a}$ = %.1f 
    $\tau_{a}$ = %.1f $\pm$ %.2f
    $R_{c}$ = %.1f 
    $T_{c}$ = %.1f $\pm$ %.1f 
    $\tau_{c}$ = %.1f $\pm$ %.2f''' %(par0, par3, par4, 
                                      params[2],perrs[2], par6, 
                                      params[3],perrs[3], params[4], perrs[4])
    return paramstring
In [41]:
fig, ax = plt.subplots()
plt.step(xdata,ydata,label='signal')
plt.plot(xdata, cppmodel.fluxfixedparams0346(xdata, *popts1), 'k-', lw=4, label='fit')
props = dict(boxstyle='round', facecolor='white', alpha=0.8)
ax.text(0.35, 0.7, cppstatsbox(popts1, perr1), transform=ax.transAxes, fontsize=18,
        verticalalignment='top', bbox=props)
ax.set_xlabel('Time post bounce [s]')
ax.set_ylabel('Events per bin [1.0 ms]')
ax.set_xlim(-3,15)
ax.set_ylim(-200,800)
plt.legend()
plt.show()