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]
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)
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]:
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
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()
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
Poniamo
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
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)