notebook.community

Edit and run



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"
#percorsoIN8192 = "/home/protoss/Documenti/TESI/DATI/candCoinc10to128.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]



amplH8 = candyH8[4]
amplL8 = candyL8[4]

amplH4 = candyH4[4]
amplL4 = candyL4[4]



freqH = numpy.concatenate([freqH8,freqH4])
freqL = numpy.concatenate([freqL8,freqL4])

amplH = numpy.concatenate([amplH8,amplH4])
amplL = numpy.concatenate([amplL8,amplL4])
#spindown1 = candyH[3]
#spindown2 = candyL[3]









    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-1-55a9b39208e5> in <module>()
      7 #percorsoIN8192 = "/home/protoss/Documenti/TESI/DATI/candCoinc10to128.mat"
      8 
----> 9 candy8192 = scipy.io.loadmat(percorsoIN8192)
     10 candyH8 = candy8192['candy'][0,0]['cand1']
     11 candyL8 = candy8192['candy'][0,0]['cand2']

~/.local/lib/python3.6/site-packages/scipy/io/matlab/mio.py in loadmat(file_name, mdict, appendmat, **kwargs)
    139     """
    140     variable_names = kwargs.pop('variable_names', None)
--> 141     MR, file_opened = mat_reader_factory(file_name, appendmat, **kwargs)
    142     matfile_dict = MR.get_variables(variable_names)
    143     if mdict is not None:

~/.local/lib/python3.6/site-packages/scipy/io/matlab/mio.py in mat_reader_factory(file_name, appendmat, **kwargs)
     62 
     63     """
---> 64     byte_stream, file_opened = _open_file(file_name, appendmat)
     65     mjv, mnv = get_matfile_version(byte_stream)
     66     if mjv == 0:

TypeError: 'NoneType' object is not iterable

MINCR per 8192 e 4096 separatam e poi ensicurv



In [2]:

    
#CR8 = 6.43154681572
#CR4 = 6.57801268577
CR8 = 6.77280188991
CR4 = 6.91232393771

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

print(critRatHfiltromedio8.size, critRatLfiltromedio8.size)

from matplotlib import pyplot
%matplotlib qt

a = pyplot.scatter(freqH,critRatH, s = 0.5, label='H')
b = pyplot.scatter(freqL,critRatL, s = 0.5, c='C3',label = 'L')

a = pyplot.scatter(frequenzefiltromedio, critRatHfiltromedio,  s = 10, label='H')
b = pyplot.scatter(frequenzefiltromedio,critRatLfiltromedio, s = 10,label = 'L')



pyplot.semilogy()
#pyplot.loglog()
pyplot.show()



In [3]:

    
#sensicurv con soglia cr
import numpy
import scipy.io
import numpy
import os
import math
import pylab
import scipy.special as scsp

#CR8 = 6.43154681572
#CR4 = 6.57801268577
CR8 = 6.77280188991
CR4 = 6.91232393771
tFft8 = 8192
tFft4 = 4096
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))



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

h0minHsog = numpy.concatenate([h0minH8,h0minH4])
h0minLsog = numpy.concatenate([h0minL8,h0minL4])

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



In [3]:

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


print(numpy.diff(h0Livingston[:4016,0]),numpy.diff(h0Hanford[:4016,0]))









    



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



In [4]:

    
#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"
#percorsoIN8192 = "/home/protoss/Documenti/TESI/DATI/candCoinc10to128.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(candy18.shape, candy14.shape, candy28.shape, candy28.shape, )

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

#candy1 = candy18
#candy2 = candy28

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

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


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

critRatH4 = candyH4[5]
critRatL4 = candyL4[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')
step = numpy.diff(h0Hanford[:4016,0])[0]
print(step)
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(freqsSens,freqsCand,critrats):  
    inizi = numpy.zeros(freqsSens.size)
    for i in numpy.arange(1,freqsSens.size):
        inizi[i] = numpy.where(freqsCand<(numpy.amin(freqsCand)+step*i))[0][-1]

    iniziUnici, quanti = numpy.unique(inizi, return_counts=True)
    iniziUnici = iniziUnici.astype(numpy.int64)

    minCR = numpy.zeros(freqsSens.size)

    inizio = iniziUnici[0]
    fine = iniziUnici[1]
    estrsx = 0
    estrdx = quanti[1]

    minCR[estrsx:estrdx] = numpy.amin(critrats[inizio:fine])

    for i in numpy.arange(1,iniziUnici.size-1):
        inizio = iniziUnici[0+i]
        fine = iniziUnici[1+i]

        estrsx = estrsx+quanti[i]
        estrdx = estrdx+quanti[i+1]

        minCR[estrsx:estrdx] = numpy.amin(critrats[inizio:fine])
        
    return minCR
    
    
    
minCRH8 = shrinka(h0Hanford[:944,0],freqH8,critRatH8)
minCRL8 = shrinka(h0Livingston[:944,0],freqL8,critRatL8)
#print(minCRH8)

minCRH4 = shrinka(h0Hanford[944:4016,0],freqH4,critRatH4)
minCRL4 = shrinka(h0Livingston[944:4016,0],freqL4,critRatL4)

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)/(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
0.125
(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]),)
7.89628814773e-26 6.97229363851e-26
5.04531748254e-22 1.12484295465e-21



In [1]:

    
from matplotlib import pyplot
import matplotlib as mpl
%matplotlib qt
mpl.rcParams['font.size'] = 14

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

#h = pylab.plot(h0Hanford[:4016,0],h0Hanford[:4016,1],'--', color = 'gray', label = 'H strain',linewidth = 0.5)
h = pylab.plot(h0Hanford[:4016,0],h0minH,label='Hanford',linewidth = 0.5)
#h = pylab.plot(h0Hanford[:4016,0],h0minHsog,label='$CR_{thr}(T_{FFT})$', color = 'C3' ,linewidth = 0.5)

#l = pylab.plot(h0Livingston[:4016,0],h0Livingston[:4016,1],'--', color = 'gray', label = 'L strain',linewidth = 0.5)
l = pylab.plot(h0Livingston[:4016,0],h0minL,label='Livingston', color = 'C3', linewidth = 0.5)
#l = pylab.plot(h0Livingston[:4016,0],h0minLsog,label='$CR_{thr}(T_{FFT})$', color = 'C3' ,linewidth = 0.5)





pyplot.loglog()

#pyplot.semilogy()
pyplot.ylim(6*1e-26,2e-21)
#pyplot.semilogy()
pyplot.legend()
pyplot.xlabel('$\\nu$ ($Hz$)', fontsize = 14)
pyplot.ylabel('$h_{0_{min}}$', fontsize = 14)
#pyplot.title('Search sensitivity: Hanford')
pyplot.title('Search sensitivity', fontsize = 16)
pyplot.show()









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-22afef3e0b00> in <module>()
      9 
     10 #h = pylab.plot(h0Hanford[:4016,0],h0Hanford[:4016,1],'--', color = 'gray', label = 'H strain',linewidth = 0.5)
---> 11 h = pylab.plot(h0Hanford[:4016,0],h0minH,label='Hanford',linewidth = 0.5)
     12 #h = pylab.plot(h0Hanford[:4016,0],h0minHsog,label='$CR_{thr}(T_{FFT})$', color = 'C3' ,linewidth = 0.5)
     13 

NameError: name 'pylab' is not defined

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

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



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

pfamin = numpy.stack((pfaH,pfaL),1)
print(minCRmin)
pfamin = numpy.amin(pfamin,1)

#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.scatter(h0Hanford[:4016,0],pfamin,label='$CR = CR_{min}(\\nu)$',s = 4)
#b = pylab.plot(h0Livingston[:4016,0],pfaL,'C3',label='L', linewidth = 0.3)
a = pylab.plot(h0Hanford[:4016,0],pfa,'--',color = 'C3',label='$CR = CR_{thr}$',linewidth = 1.5)

pyplot.semilogy()

pyplot.ylim(1e-12,1)

pyplot.legend(loc = 'lower right')
pyplot.xlabel('$\\nu$ ($Hz$)')
pyplot.ylabel('$P_{fa}$')
pyplot.title('False alarm probability')
pyplot.show()
#da dire quanto è grande lo spazio dei parametri









    



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






    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-6-d64c01e91663> in <module>()
     37 
     38 pfamin = numpy.stack((pfaH,pfaL),1)
---> 39 print(minCRmin)
     40 pfamin = numpy.amin(pfamin,1)
     41 

NameError: name 'minCRmin' is not defined



In [7]:

    
from matplotlib import pyplot
%matplotlib qt



a = pyplot.scatter(h0Hanford[:4016,0],minCRH, s = 1, label='H, $CR = CR_{min}(\\nu)$')
a = pyplot.scatter(h0Livingston[:4016,0],minCRL, s=1, color = 'C3',label='L, $CR = CR_{min}(\\nu)$')
a = pyplot.scatter(frequenzefiltromedio, critRatHfiltromedio,  s = 5, color = 'c',label='H, $CR>CR_{thr}(T_{FFT})$')
b = pyplot.scatter(frequenzefiltromedio,critRatLfiltromedio, s = 5, color = 'r',label = 'L, $CR>CR_{thr}(T_{FFT})$')

#pyplot.loglog()
pyplot.legend()
pyplot.xlabel('$\\nu$ ($Hz$)')
pyplot.ylabel('$CR$')
pyplot.title('Candidates critical ratio')
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]:

    
from matplotlib import pyplot
%matplotlib qt

minCRmin = numpy.stack((minCRH,minCRL),1)
print(minCRmin)
minCRmin = numpy.amin(minCRmin,1)

pyplot.grid(linestyle = '--')
a = pyplot.scatter(h0Hanford[:4016,0],minCRmin, s = 4, label='$CR = CR_{min}(\\nu)$')
#a = pyplot.scatter(h0Livingston[:4016,0],minCRL, s=1, color = 'C3',label='L, $CR = CR_{min}(\\nu)$')
a = pyplot.scatter(frequenzefiltromedio, (critRatHfiltromedio+critRatLfiltromedio)/2,  s = 5, color = 'C3',label='$CR>CR_{thr}(T_{FFT})$')

#pyplot.loglog()
pyplot.legend()
pyplot.xlabel('$\\nu$ ($Hz$)', fontsize = 14)
pyplot.ylabel('$CR$', fontsize = 14)
pyplot.title('Candidates critical ratio', fontsize = 16)
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









    



[[ 1.92714286  1.7705625 ]
 [ 3.31118182  2.0945    ]
 [ 2.698       1.25906667]
 ..., 
 [ 2.39822222  0.910575  ]
 [ 2.62305556  0.63823656]
 [ 0.          0.        ]]



In [10]:

    
print(minCRmin.size, (critRatHfiltromedio.size))



In [10]:









    Out[10]:





(array([], dtype=int64),)

Stima ellittic

$$h_0 = \frac{4G I}{c^4 r} \omega^2 \epsilon$$$$h_0 = \frac{4G I}{c^4 r} \left(2\pi \nu \right)^2 \epsilon$$

poniamo

$G = 6.674$, $c = 299 792 458$
$r = $
$I = $



In [9]:

    
import math
import numpy

def ellit(nu,h0, I):
    G = 6.674e-11
    c = 299792458
    r = 2.484e+20


    cost = (4*G*I/(math.pow(c,4)*r))*math.pow(2*math.pi,2)
    #h0 = cost*numpy.power(nu,2)*epsilon

    epsilon = h0/(cost*numpy.power(nu,2))
    return epsilon

epsiH = ellit(h0Hanford[:4016,0],h0minH,1e38)
epsiL = ellit(h0Livingston[:4016,0],h0minL,1e38)

epsiHsog = ellit(h0Hanford[:4016,0],h0minH,5e38)
epsiLsog = ellit(h0Livingston[:4016,0],h0minL,5e38)



In [12]:

    
from matplotlib import pyplot
%matplotlib qt

minCRmin = numpy.stack((minCRH,minCRL),1)
print(minCRmin)
minCRmin = numpy.amin(minCRmin,1)

pyplot.grid(linestyle = '--')

h = pylab.plot(h0Hanford[:4016,0],epsiH, color = 'C0',linewidth = 0.2)
h = pylab.plot(h0Hanford[:4016,0],epsiHsog, color = 'C0' ,linewidth = 0.2)
f = pyplot.fill_between(h0Hanford[:4016,0], epsiH, epsiHsog, label='Hanford',color = 'C0', alpha = 0.3)


l = pylab.plot(h0Livingston[:4016,0],epsiL, color = 'C3',linewidth = 0.2)
l = pylab.plot(h0Livingston[:4016,0],epsiLsog, color = 'C3' ,linewidth = 0.2)
f = pyplot.fill_between(h0Livingston[:4016,0], epsiL, epsiLsog, label = 'Livingston', color = 'C3', alpha = 0.3)


pyplot.loglog()
pyplot.legend()
pyplot.title('Minimum detectable ellipticity', fontsize = 16)
pyplot.xlabel('$\\nu$ ($Hz$)', fontsize = 14)
pyplot.ylabel('$\epsilon$', fontsize = 14)
pyplot.show()









    



[[ 1.92714286  1.7705625 ]
 [ 3.31118182  2.0945    ]
 [ 2.698       1.25906667]
 ..., 
 [ 2.39822222  0.910575  ]
 [ 2.62305556  0.63823656]
 [ 0.          0.        ]]

archiviate



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)



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 ]