1: analyza kompletni matice

In [7]:
%matplotlib inline 

#stazeni souboru
import sys,os
import pyfits
import numpy as np
resp=pyfits.open("/home/data/CALDB/data/swift/xrt/cpf/rmf/swxwt0s6_20010101v014.rmf")

popis formatu resp. matice http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html

organizovana do "channel subsets" (zapis ridke resp. opakujici se matice)

  • F_chan - zacatky podskupin
  • N_chan - delky podskupin
In [8]:
resp[1].columns
Out[8]:
ColDefs(
    name = 'ENERG_LO'; format = 'E'; unit = 'keV'
    name = 'ENERG_HI'; format = 'E'; unit = 'keV'
    name = 'N_GRP'; format = 'I'
    name = 'F_CHAN'; format = 'PI(218)'
    name = 'N_CHAN'; format = 'PI(218)'
    name = 'MATRIX'; format = 'PE(1021)'
)
In [261]:
# soucet hodnot v radcich - normalizace rozdeleni prandepodobnosti
sums=np.array([sum(a) for a in resp[1].data['MATRIX']])
from matplotlib import pyplot as pl
pl.plot(resp[1].data['ENERG_LO'],sums)
Out[261]:
[<matplotlib.lines.Line2D at 0x7f0ca52db1d0>]
In [266]:
# pocty nenulovych prvku ridke matice
lens=np.array([len(a) for a in resp[1].data['MATRIX']])
pl.plot(resp[1].data['ENERG_LO'],lens)
Out[266]:
[<matplotlib.lines.Line2D at 0x7f0ca51be810>]
In [274]:
#rekonstrukce jednoho radku [zvolena energie]
q=resp[1].data[810]
smat1=np.concatenate([np.arange(a,a+b) for a,b in np.r_[[q['F_CHAN'],q['N_CHAN']]].T])
emat=np.zeros(1024)
emat[smat1]=q['MATRIX']
pl.semilogy(emat)
print("energie %.2f keV"%q["ENERG_LO"])

energie 4.05 keV
In [15]:
#rekonstrukce cele matice
gmat=[]
for i in range(0,2400,10):
    q=resp[1].data[i]
    if q['N_GRP']==0: continue
    if q['N_GRP']>1:
        smat=np.concatenate([np.arange(a,a+b) for a,b in np.r_[[q['F_CHAN'],q['N_CHAN']]].T])
    else:
        smat=np.arange(q['F_CHAN'][0],q['F_CHAN'][0]+q['N_CHAN'][0])
    emat[:]=0
    emat[smat]=np.log(q['MATRIX'])
    gmat.append(emat.copy())
gmat=np.array(gmat)
pl.imshow(gmat,aspect='auto',extent=[0,1024,resp[1].data[0]['ENERG_LO'],resp[1].data[-1]['ENERG_HI']],origin='lower')
pl.colorbar()
Out[15]:
<matplotlib.colorbar.Colorbar at 0x7f0cba691c50>
In [16]:
#vertikalni rez - jeden kanal prekladan do energii
b=np.exp(gmat[:,500])
b[gmat[:,500]==0]=0
pl.plot(resp[1].data['ENERG_LO'][10:-10:10],b)
pl.xlabel('keV')
sum(b)
Out[16]:
0.13014538817520865

ridka matice prevedena na hustou bude zjednodusena (zprumerovana) na 46 energetickych a 64 instrumentalnich binu (kanalu)

In [39]:
femat=np.exp(gmat)
femat[femat==1]=0
eemat=femat[:230].reshape(23*2,5,64,16).mean(1).sum(2)
pl.imshow(np.log(eemat))

/opt/rh/python33/root/usr/lib/python3.3/site-packages/ipython-3.0.0_rc1-py3.3.egg/IPython/kernel/__main__.py:4: RuntimeWarning: divide by zero encountered in log
Out[39]:
<matplotlib.image.AxesImage at 0x7f0cb832ad50>
In [45]:
import pickle
#pickle.dump(eemat,open("rebin.rsp","wb"))
eemat=pickle.load(open("rebin.rsp","rb"))
In [40]:
# puvodne 1024 kanalu x 2400 energ. binu
# analogie prvniho obrazku po rebinovani
eemat[eemat==1]=0
pl.plot(eemat.sum(1))
Out[40]:
[<matplotlib.lines.Line2D at 0x7f0cb827ba90>]

2:now for something rebinned

kalibrace pozorovani AGN 3C 273 pomoci PN kamery druzice XMM

nastroji specificke analyzy dat (XMM SAS) se vytvori kalibracni matice pro konkretni polohu zdroje a natoceni dalekohledu (attitude)

  • spec.rmf
  • spec.arf
In [276]:
myspec="/home/limu/data/agn/3C273/0159960101/pn/agn_grpspec.pha"
sp1=pyfits.open(myspec)[1].data['COUNTS']
#sp1=sp1.reshape(len(sp1)//2,2).mean(1)
pl.semilogy(sp1)
pl.xlabel("energy[keV]")
pl.ylabel("counts")
Out[276]:
<matplotlib.text.Text at 0x7f0ca4fcc550>
In [77]:
ene=pyfits.open(myspec)[1].header
In [72]:
#data obsahuji take regrupovani, zde ho nepouzijeme
mydat.columns
Out[72]:
ColDefs(
    name = 'CHANNEL'; format = 'I'
    name = 'COUNTS'; format = 'J'; unit = 'count'
    name = 'GROUPING'; format = 'I'
    name = 'QUALITY'; format = 'I'
)
In [94]:
resp2=pyfits.open("/home/limu/data/agn/3C273/0159960101/pn/spec.rmf")
lens=np.array([len(a) for a in resp2[1].data['MATRIX']])
pl.plot(resp2[1].data['ENERG_LO'],lens)
pl.ylim(0,900)
resp2[1].columns
Out[94]:
ColDefs(
    name = 'ENERG_LO'; format = 'E'; unit = 'keV'
    name = 'ENERG_HI'; format = 'E'; unit = 'keV'
    name = 'N_GRP'; format = 'I'
    name = 'F_CHAN'; format = '3I'
    name = 'N_CHAN'; format = '3I'
    name = 'MATRIX'; format = '1PE(800)'
)
In [43]:
# new response matrix
#resp2=pyfits.open("/home/limu/data/agn/3C273/0159960101/pn/spec.rmf")
print(resp2[1].columns)
llen=800
emat=np.zeros(llen)
gmat2=[]
for i in range(0,2000,5):
    q=resp2[1].data[i]
    if q['N_GRP']==0: continue
    if q['N_GRP']>1:
        smat=np.concatenate([np.arange(a,a+b) for a,b in np.r_[[q['F_CHAN'],q['N_CHAN']]].T])
    else:
        smat=np.arange(q['F_CHAN'][0],q['F_CHAN'][0]+q['N_CHAN'][0])
    emat[:]=0
    emat[smat]=np.log(q['MATRIX'])
    gmat2.append(emat.copy())
gmat2=np.array(gmat2)
pl.imshow(gmat2,aspect='auto',extent=[0,llen,resp2[1].data[0]['ENERG_LO'],resp2[1].data[i]['ENERG_HI']],origin='lower')
pl.colorbar()
gmat2.shape

ColDefs(
    name = 'ENERG_LO'; format = 'E'; unit = 'keV'
    name = 'ENERG_HI'; format = 'E'; unit = 'keV'
    name = 'N_GRP'; format = 'I'
    name = 'F_CHAN'; format = '3I'
    name = 'N_CHAN'; format = '3I'
    name = 'MATRIX'; format = '1PE(800)'
)
Out[43]:
(400, 800)

diagonalni tvar je tady vyrazne deformovan - duvodem je ve skutecnosti nepravidelne rozdeleni energet. binu, jak je videt z ARF matice nize (srovnani s polohami binu u "velke" matice v uvodu)

In [44]:
arf2=pyfits.open("/home/limu/data/agn/3C273/0159960101/pn/spec.arf")
pl.plot(resp[1].data['ENERG_LO'])
pl.plot(resp2[1].data['ENERG_LO'])
pl.xlabel("channel")
pl.ylabel("energy")
pl.legend(["fullmat","rebinned"],loc=2)
arf2[1].data['ENERG_LO'].shape,resp2[1].data['ENERG_LO'].shape
Out[44]:
((2067,), (2067,))

zkusime model

model 1: $p_0\ x^{p_1}\ e^{-x/p_2}$

In [96]:
# defining absorbed PL - 3 parameters
model1=lambda x,p:p[0]*pow(x,p[1])*np.exp(-x/p[2]) 
x1=resp2[1].data['ENERG_LO'][2:2000:5][:400]
pl.loglog(x1,model1(x1,[2,-1.3,3.]))
Out[96]:
[<matplotlib.lines.Line2D at 0x7f0cb3444110>]
In [287]:
smat=np.exp(gmat2)
smat[gmat2==0]=0 #pri logaritmovani jsme upravili nulove body
pl.imshow(smat,origin="lower")
pl.colorbar()
#sum(gmat2==1)
Out[287]:
<matplotlib.colorbar.Colorbar at 0x7f0ca35d04d0>

v nezlogaritmovanem tvaru vidime vyrazne dominujici prvky kolem "diagonaly"
vyhledanim maximalniho prvku pro dany kanal muzeme priblizne prepocitavat instrumentalni kanaly na energii (fce argmax o nekolik poli nize)...

In [277]:
#aplikace response matice na model
pix1=np.dot(model1(x1,[2,-1.3,10]),smat)
pl.semilogy(pix1)
Out[277]:
[<matplotlib.lines.Line2D at 0x7f0ca472c0d0>]
In [279]:
llen=800
ematff=np.zeros(llen)
gmat2ff=[]
for i in range(0,2000):
    q=resp2[1].data[i]
    if q['N_GRP']==0: continue
    if q['N_GRP']>1:
        smat=np.concatenate([np.arange(a,a+b) for a,b in np.r_[[q['F_CHAN'],q['N_CHAN']]].T])
    else:
        smat=np.arange(q['F_CHAN'][0],q['F_CHAN'][0]+q['N_CHAN'][0])
    ematff[:]=0
    ematff[smat]=np.log(q['MATRIX'])
    gmat2ff.append(ematff.copy())
gmat2ff=np.array(gmat2ff)
In [286]:
#spocteno jeste jednou bez prumerovani po 5ti binech
x1b=(resp2[1].data['ENERG_LO'][:2000]+resp2[1].data['ENERG_HI'][:2000])/2.
smatff=np.exp(gmat2ff)
smatff[gmat2ff==0]=0 
pix1b=np.dot(model1(x1b,[2,-1.3,10]),smatff)
pl.semilogy(pix1b)
Out[286]:
[<matplotlib.lines.Line2D at 0x7f0ca458ce10>]

muzeme pouzit uz rebinovanou matici z analyzy

  • RMF matice 800 x 2067 (redukujeme na 800 x 400)
  • ARF matice 2067 bodu # efektivni plocha det.
In [80]:
pl.plot(resp[1].data['ENERG_LO'])
pl.plot(resp2[1].data['ENERG_LO'])
Out[80]:
[<matplotlib.lines.Line2D at 0x7f0759975310>]
In [122]:
x1=resp2[1].data['ENERG_LO'][2:2000:5][:400]
#ene2=resp2[1].data['ENERG_LO'][2:2000:5] #kazdy 5. prvek
pl.plot(x1[smat.argmax(0)])
Out[122]:
[<matplotlib.lines.Line2D at 0x7f0cb25b5a10>]
In [125]:
np.polyfit(np.r_[:800],x1[smat.argmax(0)],1)
Out[125]:
array([ 0.01501085,  0.02567662])
  • spektrum spoctene ve 400 bodech konvertujeme na odezvu v 800 kanalech (vyse diskutovana fce argmax)
  • transformace (zhruba) E[eV]=15*ch+25.6 (fit predchoziho obrazku)
  • pouziti ARF matice koriguje ucinnost detekce v kazdem spekt. bode
In [290]:
conv2=0.0150*np.r_[:800]+0.0256
wid=resp2[1].data['ENERG_HI'][2:2000:5][:400]-x1
widb=(resp2[1].data['ENERG_HI'][:2000:]-x1b)*2
x1=x1+wid/2.
pix2=np.dot(model1(x1,[2,-1.3,10])*wid,smat)
pl.semilogy(conv2,pix2)
pl.xlabel("energy")
Out[290]:
<matplotlib.text.Text at 0x7f0ca35a9290>
In [167]:
pl.plot(x1,wid)
Out[167]:
[<matplotlib.lines.Line2D at 0x7f0ca6fd8e90>]
In [285]:
pl.plot(arf2[1].data['ENERG_LO'],arf2[1].data['SPECRESP'])
pl.xlabel("energy")
pl.title("ARF matice")
Out[285]:
<matplotlib.text.Text at 0x7f0ca3893a50>
In [284]:
pl.semilogx(arf2[1].data['ENERG_LO'],arf2[1].data['SPECRESP'])
pl.xlabel("energy")
Out[284]:
<matplotlib.text.Text at 0x7f0ca45bc790>
In [178]:
arpl=arf2[1].data['SPECRESP'][:2000].reshape(400,5).mean(1) #rebinovani  ARF matice po 5 energ. binech
pix2b=np.dot(model1(x1,[1,-1.3,10])*wid*arpl,smat)
pl.semilogy(conv2,sp1)
norm=sp1[10:].sum()/pix2b[10:].sum()
pix2b=np.dot(model1(x1,[norm,-1.3,10])*wid*arpl,smat)
pl.semilogy(conv2,pix2b)
pl.xlabel("energy")
norm
Out[178]:
5697.9899146928847
In [302]:
arplb=arf2[1].data['SPECRESP'][:2000]
pix2bb=np.dot(model1(x1b,[1,-1.3,10])*widb*arplb,smatff)
normb=sp1[10:].sum()/pix2bb[10:].sum()
pix2bb=np.dot(model1(x1b,[normb,-1.3,10])*widb*arplb,smatff)
#pl.semilogy(conv2,pix2bb)
pl.semilogy(conv2,abs(pix2bb-pix2b))
Out[302]:
[<matplotlib.lines.Line2D at 0x7f0ca2fa6990>]

bez prumerovani nejsou problemy s "vlnkami" ve vyslednem spektru (ale je to lehce vypocetne narocnejsi)

In [304]:
pl.semilogy(conv2,pix2bb)
Out[304]:
[<matplotlib.lines.Line2D at 0x7f0ca3712110>]

3: fitovani

In [307]:
sel=(conv2>0.7)*(conv2<10)
#dif=lambda p:((np.dot(model1(x1,p)*wid*arpl,smat)[sel]-sp1[sel])**2/sp1[sel]).sum()
dif=lambda p:((np.dot(model1(x1b,p)*widb*arplb,smatff)[sel]-sp1[sel])**2/sp1[sel]).sum()
dif([normb,-1.3,10])
Out[307]:
13362.291510615321
In [308]:
from scipy import optimize as op
bst1=op.fmin(dif,[normb,-1.3,10])
bst1,dif(bst1)

Warning: Maximum number of function evaluations has been exceeded.
Out[308]:
(array([  1.13876051e+03,  -1.71023323e+00,   9.55943477e+06]),
 3516.8899426683133)
In [313]:
pl.semilogy(conv2[sel],sp1[sel])
val1=np.dot(model1(x1b,bst1)*widb*arplb,smatff)
pl.semilogy(conv2[sel],val1[sel])
pl.grid()
In [314]:
pl.plot(conv2[sel],(sp1[sel]-val1[sel])/np.sqrt(sp1[sel]))
pl.ylim(-6,6)
pl.ylabel('residuum')
pl.xlabel('energy')
pl.grid()

model2 : broken power-law

In [325]:
model2=lambda x,p:p[0]*np.concatenate([pow(x[x<p[2]]/p[2],p[1]),pow(x[x>=p[2]]/p[2],p[3])])
pl.loglog(x1,model2(x1,[5,-1.,3,-2.]))
pl.grid()
In [326]:
cst0=[bst1[0],bst1[1]+.5,2.,bst1[1]]
#difarr=lambda p:(np.dot(model2(x1,p)*wid*arpl,smat)[sel]-sp1[sel])/np.sqrt(sp1[sel])
difarr=lambda p:(np.dot(model2(x1b,p)*widb*arplb,smatff)[sel]-sp1[sel])/np.sqrt(sp1[sel])
#cst1=op.fmin(lambda p:sum(difarr(p)**2),cst0)
zmin=lambda p:sum(difarr(p)**2)
#zmin(cst0)
cst1=op.fmin(zmin,cst0)

Optimization terminated successfully.
         Current function value: 3464.833506
         Iterations: 332
         Function evaluations: 596
In [328]:
pl.semilogy(conv2[sel],sp1[sel])
val2=np.dot(model2(x1b,cst1)*widb*arplb,smatff)
pl.semilogy(conv2[sel],val2[sel])
pl.grid()
In [332]:
pl.plot(conv2[sel],(val2[sel]-sp1[sel])/np.sqrt(sp1[sel]))
pl.grid()

robustni postupy?

potrebujeme znat nejistoty predpovedi modelu v jednotlivych binech

zavisi take na nejistotach parametru

T.B.D. (stay tuned)

analyza nejistot

pocitame profil chi2 v zavislosti na 4 parametrech

In [211]:
s0=(difarr(cst1)**2).sum()
s0/sum(sel)
Out[211]:
7.4904398398077561
In [232]:
npt=21
zgrid=np.mgrid[-5:5:1j*npt,-5:5:1j*npt,-5:5:1j*npt,-5:5:1j*npt]
zgrid[:,12,12,12,12]
Out[232]:
array([ 1.,  1.,  1.,  1.])
In [234]:
alpars=np.array([cst1[i]*(1+0.01*zgrid[i].ravel()) for i in range(len(cst1))]).T
alchis=np.array([(difarr(c)**2).sum() for c in alpars]).reshape(npt,npt,npt,npt)
In [237]:
pl.imshow(alchis[10,10]/s0*sum(sel))
pl.colorbar()
Out[237]:
<matplotlib.colorbar.Colorbar at 0x7f0ca6793310>
In [239]:
pl.imshow(alchis[10,:,10,:]/s0*sum(sel))
pl.colorbar()
Out[239]:
<matplotlib.colorbar.Colorbar at 0x7f0ca6b15050>
In [240]:
pl.imshow(alchis.mean(2).mean(0)/s0*sum(sel))
pl.colorbar()
Out[240]:
<matplotlib.colorbar.Colorbar at 0x7f0ca61fe4d0>
In [254]:
xp=np.r_[1-0.05:1.05:npt*1j]
pl.plot(xp,alchis[10,10,10,:])
pl.plot(xp,alchis[10,10,8,:])
pl.plot(xp,alchis[10,14,10,:])
pl.xlabel("relat. p4")
Out[254]:
<matplotlib.text.Text at 0x7f0ca55021d0>
In [249]:
#der2=lambda p:p[0]
print(1/np.polyfit(xp-1,alchis[11,11,11,:]/s0,2)[0]/2.)
print(1/np.polyfit(xp[5:16]-1,alchis[11,11,11,5:16]/s0,2)[0]/2.)
print(1/np.polyfit(xp[5:16]-1,alchis[11,11,8,5:16]/s0,2)[0]/2.)
print(1/np.polyfit(xp[5:16]-1,alchis[11,14,11,5:16]/s0,2)[0]/2.)

0.00261408619715
0.00263046872989
0.00276578452421
0.00263044339036
In [252]:
print(1/np.sqrt(np.polyfit(xp[5:17]-1,alchis[5:17,11,11,11]/s0,2)[0]/2.))
print(1/np.sqrt(np.polyfit(xp[5:17]-1,alchis[11,5:17,11,11]/s0,2)[0]/2.))
print(1/np.sqrt(np.polyfit(xp[5:17]-1,alchis[11,11,5:17,11]/s0,2)[0]/2.))
print(1/np.sqrt(np.polyfit(xp[5:17]-1,alchis[11,11,11,5:17]/s0,2)[0]/2.))

0.131418304118
0.505813946525
0.075717277037
0.103465065099
In [ ]:
pl.plot(xp,alchis[10,:,10,10]/s0,label="p2")
pl.plot(xp,alchis[10,10,10,:]/s0,label="p4")
pl.legend()
Out[ ]:
<matplotlib.legend.Legend at 0x7f0ca533f6d0>
In [260]:
#4D matice zazalohovana
import pickle
pickle.dump(alchis,open("chi2_bknpow.npy","wb"))