notebook.community

Edit and run



In [11]:

    
import scipy.io
import numpy
import os

percorsoIN4096 = "DATI/cands4096/candCoinTOT.mat"
percorsoIN8192 = "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)

print(candy1.shape,candy2.shape)


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

critRat1 = candy1[5]
critRat2 = candy2[5]

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

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

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









    



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



In [ ]:

    
MINCR per 8192 e 4096 separatam e poi ensicurv



In [9]:

    
from matplotlib import pyplot
%matplotlib notebook


a = pyplot.scatter(coo1[0],coo1[1], s = 5, label='H')
b = pyplot.scatter(coo2[0],coo2[1], s = 5, label = 'L')
#a = pyplot.scatter(freq1,spindown1*1e9, s = 0.1, label='H')
#b = pyplot.scatter(freq2,spindown2*1e9, s = 0.1, label = 'L')



#pyplot.loglog()
pyplot.legend();

pyplot.show()



In [10]:

    
from matplotlib import pyplot
%matplotlib notebook


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

pyplot.show()



In [14]:

    
from matplotlib import pyplot
%matplotlib qt


a = pyplot.scatter(freq1,critRat1, s = 0.5, label='H')
b = pyplot.scatter(freq2,critRat2, s = 0.5, c='C3',label = 'L')
pyplot.loglog()
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 [8]:

    
sogliaCR8 = 6.43154681572
sogliaCR4 = 6.57801268577


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

critRatMedio = (critRat1+critRat2)/2
filtromedio = numpy.where(critRatMedio > sogliaCR)

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

critRat1filtromedio = critRat1[filtromedio]
critRat2filtromedio = critRat2[filtromedio]

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]
frequenzefiltromedio = freq1[filtromedio]
frequenzefiltroboth = freq1[nonzeri1]
spindowni1filtroboth = spindown1[nonzeri1]
spindowni2filtroboth = spindown2[nonzeri2]

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









    



6.1548125 6.0705



In [9]:

    
print(critRat1filtromedio.size, critRat2filtromedio.size)

from matplotlib import pyplot
%matplotlib notebook

pyplot.figure(figsize=(10, 8))
a = pyplot.scatter(frequenzefiltromedio, critRat1filtromedio, s = 10)
b = pyplot.scatter(frequenzefiltromedio,critRat2filtromedio, s = 10)
#pyplot.semilogx()
pyplot.loglog()
pyplot.show()



In [10]:

    
from matplotlib import pyplot
%matplotlib notebook

pyplot.figure(figsize=(10, 8))
a = pyplot.scatter(frequenzefiltro1,critRat1filtro1, s = 10)
b = pyplot.scatter(frequenzefiltro1,critRat2filtro1, s = 10)
#pyplot.semilogx()
pyplot.loglog()

pyplot.show()



In [11]:

    
from matplotlib import pyplot
%matplotlib notebook

pyplot.figure(figsize=(10, 8))
a = pyplot.scatter(frequenzefiltro2,critRat1filtro2, s = 10)
b = pyplot.scatter(frequenzefiltro2,critRat2filtro2, s = 10)
#pyplot.semilogx()
pyplot.loglog()
pyplot.show()



In [12]:

    
from matplotlib import pyplot
%matplotlib notebook

pyplot.figure(figsize=(10, 8))
a = pyplot.scatter(frequenzefiltroboth,critRat1filtroboth, s = 10)
b = pyplot.scatter(frequenzefiltroboth,critRat2filtroboth, s = 10)
#pyplot.semilogx()
pyplot.loglog()

pyplot.show()



In [1]:

    
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))









    



(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]),)



In [4]:

    
import pylab
%matplotlib qt

#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()









    Out[4]:





(5e-24, 3e-19)



In [89]:

    
#sensicurv

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


#percorsoIN4096 = "DATI/cands4096/candCoinTOT.mat"
percorsoIN8192 = "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))









    



(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]),)

$$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 [91]:

    
tFft8 = 4096
tFft4 = 8192
tObs = 9*30*24*60*60
Ntempi = tObs/(tFft/2)*0.6
Ntempi = numpy.int(Ntempi)
print(Ntempi)



In [90]:

    
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)



In [61]:

    
from matplotlib import pyplot
%matplotlib qt


a = pyplot.scatter(h0Hanford[:944,0],minCRH, s = 1, label='H')
a = pyplot.scatter(h0Livingston[:944,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



In [92]:

    
#VARIO N

def constizza(CRs,t):
    N = Ntempi#numpy.linspace(1,Ntempi,10*Ntempi)
    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)
const0minL8 = constizza(minCRL8,tFft8)

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



In [93]:

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

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



In [94]:

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

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









    



9.70955156002e-27 8.52419517708e-27
4.56240405396e-23 1.26277804719e-22



In [95]:

    
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.7)
b = pylab.plot(h0Livingston[:4016,0],h0minL,label='L', linewidth = 0.7)

pyplot.loglog()
pyplot.ylim(7*1e-27,2*1e-20)
#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



In [107]:

    
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,label='L', s= 0.8)


#pyplot.loglog()

#pyplot.ylim(0,0.05)

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



In [ ]: