In [1]:
import scipy.io
import numpy
import os

percorsoIN4096 = "/home/protoss/Documenti/TESI/DATI/cands4096/candCoinTOT.mat"
percorsoIN8192 = "/home/protoss/Documenti/TESI/DATI/candCoinc8192.mat"

candy8192 = scipy.io.loadmat(percorsoIN8192)
candyH8 = candy8192['candy'][0,0]['cand1']
candyL8 = candy8192['candy'][0,0]['cand2']



candy4096 = scipy.io.loadmat(percorsoIN4096)
candyH4 = candy4096['candy'][0,0]['cand1']
candyL4 = candy4096['candy'][0,0]['cand2']

print(candyH8.shape, candyH4.shape, candyL8.shape, candyL8.shape, )

candyH = numpy.concatenate([candyH8,candyH4],1)
candyL = numpy.concatenate([candyL8,candyL4],1)

print(candyH.shape,candyL.shape)


freqH = candyH[0]
freqL = candyL[0]

critRatH = candyH[5]
critRatL = candyL[5]

freqH8 = candyH8[0]
freqL8 = candyL8[0]

critRatH8 = candyH8[5]
critRatL8 = candyL8[5]

freqH4 = candyH4[0]
freqL4 = candyL4[0]

critRatH4 = candyH4[5]
critRatL4 = candyL4[5]



#ampl1 = candyH[4]
#ampl2 = candyL[4]

#spindown1 = candyH[3]
#spindown2 = candyL[3]


(9, 38575) (9, 89217) (9, 38575) (9, 38575)
(9, 127792) (9, 127792)

In [4]:
percorso = '/home/protoss/Documenti/TESI/DATI/lineeH8.txt'
noteH8 = numpy.loadtxt(percorso)

percorso = '/home/protoss/Documenti/TESI/DATI/lineeH4.txt'
noteH4 = numpy.loadtxt(percorso)

percorso = '/home/protoss/Documenti/TESI/DATI/lineeL8.txt'
noteL8 = numpy.loadtxt(percorso)
percorso = '/home/protoss/Documenti/TESI/DATI/lineeL4.txt'
noteL4 = numpy.loadtxt(percorso)

print(noteH8,noteH4,noteL8,noteL4)


[ 11.111175  11.11123   11.394784  11.395279  19.07328   28.61      29.8019
  35.7048    35.7065    35.7624    35.7628    35.9       36.7       37.3
  60.        99.979     99.9987  ] [ 299.6   302.22  303.31  331.9   500.06  500.21  501.09  501.2   501.26
  501.45  501.61  501.68  501.74  501.81  502.62  502.74  503.    503.11
  504.8   504.87  505.58  505.71  505.8   506.92  507.16  507.19  507.39
  508.    508.14  508.21  508.28  508.64] [ 22.7         23.3         23.9         31.4         32.6         33.8
  60.          99.99925625] [ 306.0596  306.2186  307.3476  307.5051  315.0471  331.3     499.6042
  503.1825  513.332 ]

In [ ]:
#pulizia righe note

for i in freqH8:

MINCR per 8192 e 4096 separatam e poi ensicurv


In [2]:
from matplotlib import pyplot
%matplotlib qt

#a = pyplot.scatter(freq1,ampl1, s = 2, label='H')
#b = pyplot.scatter(freq2,ampl2, s = 2, label = 'L')

a = pyplot.scatter(freqH,critRatH, s = 0.5, label='H')
b = pyplot.scatter(freqL,critRatL, s = 0.5, c='C3',label = 'L')
pyplot.semilogy()
pyplot.legend()
pyplot.xlabel('$\\nu$ ($Hz$)')
pyplot.ylabel('$CR$')
pyplot.show()

#a = pylab.plot(freq1,critRat1,  label='H')
#b = pylab.plot(freq2,critRat2, c='C3',label = 'L')
#pylab.loglog()
#pylab.legend()
#pylab.xlabel('$\\nu$ ($Hz$)')
#pylab.ylabel('$CR$')
#pyplot.show()



# vedere num cand in funzione di frequenza, ogni hertz TODO nb deve fluttuare al massimo in un fattore 2


#se c'è molto disturbo ci si aspetta ci siano molto pochi secondari, per questo il fattore 2

In [22]:
CR8 = 6.43154681572
CR4 = 6.57801268577

critRatMedio8 = (critRatH8+critRatL8)/2
filtromedio8 = numpy.where(critRatMedio8 > CR8)

critRatHfiltromedio8 = critRatH8[filtromedio8]
critRatLfiltromedio8 = critRatL8[filtromedio8]

frequenzefiltromedio8 = freqH8[filtromedio8]

critRatMedio4 = (critRatH4+critRatL4)/2
filtromedio4 = numpy.where(critRatMedio4 > CR4)

critRatHfiltromedio4 = critRatH4[filtromedio4]
critRatLfiltromedio4 = critRatL4[filtromedio4]

frequenzefiltromedio4 = freqH4[filtromedio4]


frequenzefiltromedio = numpy.concatenate((frequenzefiltromedio8,frequenzefiltromedio4))
critRatHfiltromedio = numpy.concatenate((critRatHfiltromedio8,critRatHfiltromedio4))
critRatLfiltromedio = numpy.concatenate((critRatLfiltromedio8,critRatLfiltromedio4))

In [25]:



351 351

In [67]:
#sensicurv con soglia cr
import numpy

tFft8 = 4096
tFft4 = 8192
tObs = 9*30*24*60*60
Ntempi8 = tObs/(tFft8/2)*0.6
Ntempi8 = numpy.int(Ntempi8)
Ntempi4 = tObs/(tFft4/2)*0.6
Ntempi4 = numpy.int(Ntempi4)

h0Livingston = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2017-07-20_C01_L1_O2_Sensitivity_strain_asd.txt')
h0Hanford = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2016-12-08_C01_H1_O2_Sensitivity_strain_asd.txt')

print(numpy.where(h0Livingston[:,0]>128))
print(numpy.where(h0Hanford[:,0]>128))

print(numpy.where(h0Livingston[:,0]>512))
print(numpy.where(h0Hanford[:,0]>512))


tFft8 = 8192
tFft4 = 4096

def constizza(CRs,tFft,N):
    theta = 2.5
    sogliaCR = CRs
    gamma = 0.9545


    #p0 è prob di avere picco in peakmap
    p0 = math.exp(-theta)-math.exp(-2*theta)+1/3*math.exp(-3*theta)
    p1 = math.exp(-theta)-2*math.exp(-2*theta)+math.exp(-3*theta)

    probs = p0*(1-p0)/(N*math.pow(p1,2))
    confs = sogliaCR - math.sqrt(2)*scsp.erfcinv(2*gamma)

    const0min = 4.02*math.pow(N,-1/4)*math.pow(theta,-1/2)*math.pow(probs, 1/4)*numpy.power(confs, 1/2)*math.pow(tFft,-1/2)

    return const0min

#lambda0min = (2 / theta**2) * numpy.power(probs, 1/2)*confs
#lambda0min

const0minH8 = constizza(CR8,tFft8,Ntempi8)
const0minL8 = constizza(CR8,tFft8,Ntempi8)

h0minH8 = const0minH8*h0Hanford[:944,1]
h0minL8 = const0minL8*h0Livingston[:944,1]

const0minH4 = constizza(CR4,tFft4,Ntempi4)
const0minL4 = constizza(CR4,tFft4,Ntempi4)

h0minH4 = const0minH4*h0Hanford[944:4016,1]
h0minL4 = const0minL4*h0Livingston[944:4016,1]

print(numpy.amin(h0minH), numpy.amin(h0minL))
print(numpy.amax(h0minH), numpy.amax(h0minL))

h0minH = numpy.concatenate([h0minH8,h0minH4])
h0minL = numpy.concatenate([h0minL8,h0minL4])


from matplotlib import pyplot
%matplotlib qt


#a = pyplot.scatter(h0Hanford[:4016,0],h0minH, s = 1, label='H')
#a = pyplot.scatter(h0Livingston[:4016,0],h0minL, s = 1, label='L')


numpy.save('sensiH',h0minH)
numpy.save('sensiL',h0minL)

a = pylab.plot(h0Hanford[:4016,0],h0minH,label='H',linewidth = 0.5)
b = pylab.plot(h0Livingston[:4016,0],h0minL,'C3',label='L', linewidth = 0.5)

pyplot.loglog()
pyplot.ylim(7*1e-27,2*1e-22)
#pyplot.semilogy()
pyplot.legend()
pyplot.xlabel('$\\nu$ ($Hz$)')
pyplot.ylabel('$h_{0_{min}} \;(1/\sqrt{Hz})$')
pyplot.show()
# vedere num cand in funzione di frequenza, ogni hertz TODO nb deve fluttuare al massimo in un fattore 2


#se c'è molto disturbo ci si aspetta ci siano molto pochi secondari, per questo il fattore 2


(array([  944,   945,   946, ..., 39918, 39919, 39920]),)
(array([  944,   945,   946, ..., 39918, 39919, 39920]),)
(array([ 4016,  4017,  4018, ..., 39918, 39919, 39920]),)
(array([ 4016,  4017,  4018, ..., 39918, 39919, 39920]),)
1.12706978679e-26 1.0007394402e-26
4.38830931319e-23 1.38230408679e-22
Warning: Cannot change to a different GUI toolkit: qt. Using notebook instead.

In [87]:
%matplotlib notebook
pfa8 = scsp.erfc(CR8/math.sqrt(2))/2

pfa4 = scsp.erfc(CR4/math.sqrt(2))/2

pfa8 = numpy.ones(944)*pfa8
pfa4 = numpy.ones(4016-944)*pfa4

pfa = numpy.concatenate([pfa8,pfa4])

#a = pylab.plot(h0Hanford[:4016,0],pfaH,label='H',linewidth = 0.7)
#b = pylab.plot(h0Livingston[:4016,0],pfaL,label='L', linewidth = 0.7)
#a = pyplot.scatter(h0Hanford[:4016,0],pfaH,label='H',s = 0.8)
#b = pyplot.scatter(h0Livingston[:4016,0],pfaL,color = 'C3',label='L', s= 0.8)
print(pfa.size,h0Hanford[:4016,0].size)
#a = pylab.plot(h0Hanford[:4016,0],pfa,'C1',label='$CR_{thr}$',linewidth = 1)

#pyplot.semilogy()

#pyplot.ylim(1e-7,1)

#pyplot.legend(loc='upper left')
#pyplot.xlabel('$\\nu$ ($Hz$)')
#pyplot.ylabel('$P_{fa}$')
#pyplot.show()


4016 4016

In [ ]:
import pylab
%matplotlib qt

import numpy
h0Livingston = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2017-07-20_C01_L1_O2_Sensitivity_strain_asd.txt')
h0Hanford = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2016-12-08_C01_H1_O2_Sensitivity_strain_asd.txt')

print(numpy.where(h0Livingston[:,0]>128))
print(numpy.where(h0Hanford[:,0]>128))

print(numpy.where(h0Livingston[:,0]>512))
print(numpy.where(h0Hanford[:,0]>512))


#b = pyplot.scatter(h0Hanford[:4015,0],h0Hanford[:4015,1], s = 5,label = 'H')
#a = pyplot.scatter(h0Livingston[:4015,0],h0Livingston[:4015,1], s = 5, label = 'L')

#a = pylab.plot(h0Hanford[:4016,0],h0Hanford[:4016,1],label='H',linewidth = 0.7)
#b = pylab.plot(h0Livingston[:4016,0],h0Livingston[:4016,1],label='L', linewidth = 0.7)

b = pylab.plot(h0Hanford[:,0],h0Hanford[:,1], label = 'Hanford',linewidth = 0.5)
#a = pylab.plot(h0Livingston[:,0],h0Livingston[:,1],linewidth = 0.1)
#a = pylab.plot(h0Livingston[:,0],h0Livingston[:,1],linewidth = 0.1)
a = pylab.plot(h0Livingston[:,0],h0Livingston[:,1],'C3', label = 'Livingston',linewidth = 0.5)


#c = pyplot.scatter(frequenzefiltroboth,critRat1filtroboth*1e-21, s = 10)
#d = pyplot.scatter(frequenzefiltroboth,critRat2filtroboth*1e-21, s = 10)
#c = pyplot.scatter(frequenzefiltromedio,critRat1filtromedio*1e-21, s = 10)
#d = pyplot.scatter(frequenzefiltromedio,critRat2filtromedio*1e-21, s = 10)


ymin = numpy.amin(h0Livingston[:,1])
ymax = numpy.amax(h0Livingston[:,1])
pylab.ylim((ymin,ymax))
#pyplot.semilogy()
pylab.ylabel('Strain ($1/\sqrt{Hz}$)')
pylab.xlabel('$\\nu$ ($Hz$)')
pylab.legend()
pylab.loglog()
#pylab.ylim(4*1e-24,1.5*1e-19)
pylab.ylim(5*1e-24,3*1e-19)
#pylab.show()

In [84]:
#sensicurv senza soglia cr

import scipy.io
import numpy
import os
import math
import pylab
import scipy.special as scsp

tFft8 = 4096
tFft4 = 8192
tObs = 9*30*24*60*60
Ntempi8 = tObs/(tFft8/2)*0.6
Ntempi8 = numpy.int(Ntempi8)
Ntempi4 = tObs/(tFft4/2)*0.6
Ntempi4 = numpy.int(Ntempi4)

print(Ntempi8,Ntempi4)


percorsoIN4096 = "/home/protoss/Documenti/TESI/DATI/cands4096/candCoinTOT.mat"
percorsoIN8192 = "/home/protoss/Documenti/TESI/DATI/candCoinc8192.mat"

candy8192 = scipy.io.loadmat(percorsoIN8192)
candy18 = candy8192['candy'][0,0]['cand1']
candy28 = candy8192['candy'][0,0]['cand2']



candy4096 = scipy.io.loadmat(percorsoIN4096)
candy14 = candy4096['candy'][0,0]['cand1']
candy24 = candy4096['candy'][0,0]['cand2']

#print(candy18.shape, candy14.shape, candy28.shape, candy28.shape, )

#candy1 = numpy.concatenate([candy18,candy14],1)
#candy2 = numpy.concatenate([candy28,candy24],1)

#candy1 = candy18
#candy2 = candy28

freq1 = candy1[0]
freq2 = candy2[0]

critRat18 = candy18[5]
critRat28 = candy28[5]

critRat14 = candy14[5]
critRat24 = candy24[5]



ampl1 = candy1[4]
ampl2 = candy2[4]

spindown1 = candy1[3]
spindown2 = candy2[3]

coo1=candy1[6:8]
coo2=candy2[6:8]


h0Livingston = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2017-07-20_C01_L1_O2_Sensitivity_strain_asd.txt')
h0Hanford = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2016-12-08_C01_H1_O2_Sensitivity_strain_asd.txt')

print(numpy.where(h0Livingston[:,0]>128))
print(numpy.where(h0Hanford[:,0]>128))

print(numpy.where(h0Livingston[:,0]>512))
print(numpy.where(h0Hanford[:,0]>512))

def shrinka(sens,critrats):
    size = sens.size

    rapportoInd = numpy.int(critrats.size/size)

    minCR = numpy.zeros(size)
    for i in numpy.arange(size-1):
        inizio = i*rapportoInd
        fine = (i+1)*rapportoInd
        minCR[i] = numpy.amin(critrats[inizio:fine])

    inizio = (size-2)*rapportoInd
    minCR[size-1] =numpy.amin(critrats[inizio:])
    return minCR
    
    
    
minCRH8 = shrinka(h0Hanford[:944,0],critRat18)
minCRL8 = shrinka(h0Livingston[:944,0],critRat28)

minCRH4 = shrinka(h0Hanford[944:4016,0],critRat14)
minCRL4 = shrinka(h0Livingston[944:4016,0],critRat24)

minCRH = numpy.concatenate((minCRH8,minCRH4))
minCRL = numpy.concatenate((minCRL8,minCRL4))


import math
def constizza(CRs,tFft,N):
    theta = 2.5
    sogliaCR = CRs
    gamma = 0.9545


    #p0 è prob di avere picco in peakmap
    p0 = math.exp(-theta)-math.exp(-2*theta)+1/3*math.exp(-3*theta)
    p1 = math.exp(-theta)-2*math.exp(-2*theta)+math.exp(-3*theta)

    probs = p0*(1-p0)/(N*math.pow(p1,2))
    confs = sogliaCR - math.sqrt(2)*scsp.erfcinv(2*gamma)

    const0min = 4.02*math.pow(N,-1/4)*math.pow(theta,-1/2)*math.pow(probs, 1/4)*numpy.power(confs, 1/2)*math.pow(tFft,-1/2)

    return const0min

#lambda0min = (2 / theta**2) * numpy.power(probs, 1/2)*confs
#lambda0min

const0minH8 = constizza(minCRH8,tFft8,Ntempi8)
const0minL8 = constizza(minCRL8,tFft8,Ntempi8)

h0minH8 = const0minH8*h0Hanford[:944,1]
h0minL8 = const0minL8*h0Livingston[:944,1]

const0minH4 = constizza(minCRH4,tFft4,Ntempi4)
const0minL4 = constizza(minCRL4,tFft4,Ntempi4)

h0minH4 = const0minH4*h0Hanford[944:4016,1]
h0minL4 = const0minL4*h0Livingston[944:4016,1]



h0minH = numpy.concatenate([h0minH8,h0minH4])
h0minL = numpy.concatenate([h0minL8,h0minL4])


print(numpy.amin(h0minH), numpy.amin(h0minL))
print(numpy.amax(h0minH), numpy.amax(h0minL))

from matplotlib import pyplot
%matplotlib qt


#a = pyplot.scatter(h0Hanford[:4016,0],h0minH, s = 1, label='H')
#a = pyplot.scatter(h0Livingston[:4016,0],h0minL, s = 1, label='L')


numpy.save('sensiH',h0minH)
numpy.save('sensiL',h0minL)

#a = pylab.plot(h0Hanford[:4016,0],h0Hanford[:4016,1],'--',label='H',linewidth = 0.3)
#b = pylab.plot(h0Livingston[:4016,0],h0Livingston[:4016,1],'--',label='L', linewidth = 0.3)
a = pylab.plot(h0Hanford[:4016,0],h0minH,label='H',linewidth = 0.5)
b = pylab.plot(h0Livingston[:4016,0],h0minL,'C3',label='L', linewidth = 0.5)

pyplot.loglog()
pyplot.ylim(7*1e-27,2*1e-22)
#pyplot.semilogy()
pyplot.legend()
pyplot.xlabel('$\\nu$ ($Hz$)')
pyplot.ylabel('$h_{0_{min}} \;(1/\sqrt{Hz})$')
pyplot.show()
# vedere num cand in funzione di frequenza, ogni hertz TODO nb deve fluttuare al massimo in un fattore 2


#se c'è molto disturbo ci si aspetta ci siano molto pochi secondari, per questo il fattore 2


6834 3417
(array([  944,   945,   946, ..., 39918, 39919, 39920]),)
(array([  944,   945,   946, ..., 39918, 39919, 39920]),)
(array([ 4016,  4017,  4018, ..., 39918, 39919, 39920]),)
(array([ 4016,  4017,  4018, ..., 39918, 39919, 39920]),)
9.70955156002e-27 8.52419517708e-27
4.56240405396e-23 1.26277804719e-22
Warning: Cannot change to a different GUI toolkit: qt. Using notebook instead.
$$h_{0,min} \approx \frac{4.02}{N^{1/4}\theta_{thr}^{1/2}} \left(\frac{p_0(1-p_0)}{p_1^2}\right)^{1/4}\sqrt{CR_{thr}-\sqrt{2}erfc^{-1}(2\Gamma)} \cdot \sqrt{\frac{S_n(\nu)}{T_{FFT}}}$$

Poniamo

  • $\Gamma = 0.9545$
  • $\theta_{thr} = 2.5$ (TODO CONTROL)
  • $T_{FFT} = 8192$ (o $4096$)
  • $N \sim 3400$ (TODO CONTROL)
  • $p_0 = 0.0755$, $p_1 = 0.0692$

In [89]:
%matplotlib notebook
pfaH8 = scsp.erfc(minCRH8/math.sqrt(2))/2
pfaL8 = scsp.erfc(minCRL8/math.sqrt(2))/2

pfaH4 = scsp.erfc(minCRH4/math.sqrt(2))/2
pfaL4 = scsp.erfc(minCRL4/math.sqrt(2))/2

pfaH = numpy.concatenate([pfaH8,pfaH4])
pfaL = numpy.concatenate([pfaL8,pfaL4])

#a = pylab.plot(h0Hanford[:4016,0],pfaH,label='H',linewidth = 0.7)
#b = pylab.plot(h0Livingston[:4016,0],pfaL,label='L', linewidth = 0.7)
#a = pyplot.scatter(h0Hanford[:4016,0],pfaH,label='H',s = 0.8)
#b = pyplot.scatter(h0Livingston[:4016,0],pfaL,color = 'C3',label='L', s= 0.8)

a = pylab.plot(h0Hanford[:4016,0],pfaH,label='H',linewidth = 0.3)
b = pylab.plot(h0Livingston[:4016,0],pfaL,'C3',label='L', linewidth = 0.3)
a = pylab.plot(h0Hanford[:4016,0],pfa,'--',color = 'C4',label='$CR_{thr}$',linewidth = 1.5)

pyplot.semilogy()

pyplot.ylim(1e-11,1)

pyplot.legend(loc = 'lower right')
pyplot.xlabel('$\\nu$ ($Hz$)')
pyplot.ylabel('$P_{fa}$')
pyplot.show()



In [41]:
from matplotlib import pyplot
%matplotlib qt


a = pyplot.scatter(h0Hanford[:4016,0],minCRH, s = 1, label='H')
a = pyplot.scatter(h0Livingston[:4016,0],minCRL, s = 1, label='L')
#pyplot.loglog()
pyplot.legend()
pyplot.xlabel('$\\nu$ ($Hz$)')
pyplot.ylabel('$CR$')
pyplot.show()
# vedere num cand in funzione di frequenza, ogni hertz TODO nb deve fluttuare al massimo in un fattore 2


#se c'è molto disturbo ci si aspetta ci siano molto pochi secondari, per questo il fattore 2


Warning: Cannot change to a different GUI toolkit: qt. Using notebook instead.

In [ ]:
CR8 = 6.43154681572
CR4 = 6.57801268577


#filtro1 = numpy.where(critRat1 > sogliaCR)
#filtro2 = numpy.where(critRat2 > sogliaCR)

critRatMedio8 = (critRatH8+critRatL8)/2
filtromedio8 = numpy.where(critRatMedio8 > CR8)

#critRat1filtro1 = critRat1[filtro1]
#critRat2filtro1 = critRat2[filtro1]
#critRat1filtro2 = critRat1[filtro2]
#critRat2filtro2 = critRat2[filtro2]

critRatHfiltromedio8 = critRatH8[filtromedio8]
critRatLfiltromedio8 = critRatL8[filtromedio8]

#annullatore1 = numpy.where(critRat1 <= sogliaCR)
#annullatore2 = numpy.where(critRat2 <= sogliaCR)
#critRat1[annullatore1] = 0
#critRat1[annullatore2] = 0
#critRat2[annullatore1] = 0
#critRat2[annullatore2] = 0

#nonzeri1 = numpy.nonzero(critRat1)
#nonzeri2 = numpy.nonzero(critRat2)

#critRat1filtroboth = critRat1[nonzeri1[0]]#[filtroboth]
#critRat2filtroboth = critRat2[nonzeri2[0]]#[filtroboth]

#print(numpy.amin(critRat1filtroboth), numpy.amin(critRat2filtroboth))


#frequenzefiltro1 = freq1[filtro1]
#frequenzefiltro2 = freq2[filtro2]
frequenzefiltromedio8 = freqH8[filtromedio8]
#frequenzefiltroboth = freq1[nonzeri1]
#spindowni1filtroboth = spindown1[nonzeri1]
#spindowni2filtroboth = spindown2[nonzeri2]

#print(critRat1filtro1.size, critRat1filtro2.size, critRat1filtroboth.size, critRat2filtro1.size, critRat2filtro2.size, critRat2filtroboth.size)