diff --git a/README.md b/README.md index 1e21990..67a639f 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ The preprocessed data will be continuously updated until INTEGRAL stops register # INSTALL -## Step 1. (optionally) Creation of Python virtual environment +## Step 1. Creation of Python virtual environment ``` mkdir ~/work; cd ~/work python3 -m venv venv diff --git a/figures/cutoffpl_gb.r b/figures/cutoffpl_gb.r new file mode 100644 index 0000000..fd468d9 --- /dev/null +++ b/figures/cutoffpl_gb.r @@ -0,0 +1,98 @@ +# +# 2 data sets +# in Xspec make ipl->wdata +# then split file to four parts as lecr1/2/3/4 +# + + +library(sfsmisc) + + +up=2.0 +#xgrid <- c(30,50,70,100); +Emin=25.0 +Emax=180.0 +postscript('cutoffpl_gb.eps', horizontal = FALSE, onefile = FALSE, paper = "special",width = 9.0, height = 6.0) +par(oma=c(1,1,1,1)) +par(lwd=2) + +layout(matrix(c(1,2), 2, 1, byrow = TRUE), heights = c(1.5,1), TRUE) +leftx=5.3 +par(mar=c(0.1, leftx, 1, 1)) +par(cex.lab=1.6) +par(cex.axis=1.6) +xlim <- c(Emin,Emax) +ylim <- c(0.2,2) + +cl="black" +a <- read.table('./cutoffpl_gb_eeuf.dat', col.names=c("x","dx","y","dy","total","cutoffpl","pow")) +plot(a$x, a$y, pch=3, bg="white", col="white", ylim=ylim, xlim=xlim, ylab=expression("keV"^"2"~"(Phot. keV"^"-1"~"cm"^"-2"~"s"^"-1"*")"),type="p",xaxt = 'n',xlab="",log="xy") + +#axis(side=2, at=c(1e-6,1e-5, 1e-4, 1e-3, 1e-2), labels=expression('-6','-5','-4','-3','-2')) + +#abline(v=xgrid, col="lightgray", lty="dotted") +grid() + +################################################################################### +## +## Plot model components +## +################################################################################### + +lines(a$x,a$cutoffpl,pch=1,bg=cl,col="blue",lwd=3,lty="longdash") +lines(a$x,a$pow,pch=1,bg=cl,col="red",lwd=3,lty="dotdash") +lines(a$x,a$total,pch=1,bg=cl,col="black") + + +# +# Plot upper limits +# +upx=a$x[(a$y/a$dy)=up] +py=a$y[(a$y/a$dy)>=up] +pdx=a$dx[(a$y/a$dy)>=up] +pdy=a$dy[(a$y/a$dy)>=up] + +segments(px,py-pdy,px,py+pdy,col=cl) +segments(px-pdx,py,px+pdx,py,col=cl) + +text(100, 1.5, "GB",cex=2.5) + +################################################################################### + +par(mar=c(4.5, leftx, 0.1, 1)) +ylim <- c(-3.1,3.1) + +cl="black" +a <- read.table("./cutoffpl_gb_delchi.dat", col.names=c("x","dx","y","dy")) +plot(a$x, a$y, pch=3, ylim=ylim, xlim=xlim, ylab=expression(Delta~chi), xlab="",type="p",log="x",xaxt = 'n') +grid() +segments(a$x,a$y-a$dy,a$x,a$y+a$dy,col=cl) +segments(a$x-a$dx,a$y,a$x+a$dx,a$y,col=cl) + +#abline(v=xgrid, col="lightgray", lty="dotted") +#grid() +#abline(h=0, col = "black",lty=2) +#dev.off() + +#abline(v=xgrid, col="lightgray", lty="dotted") +abline(h=0, col = "black",lty=2) + +mtext(side = 1, text = "Energy, keV", line = 4, cex=1.6) + +#axis(1, mgp=c(3, 1.5, 0)) +### Log axis using sfsmisc ### +atx=c(8,30,50,80,100,150) +eaxis(1, at = atx, labels = pretty10exp(atx, sub10=c(1,100), drop.1=TRUE), las=0) +abline(v=atx, col="lightgray", lty="dotted") +### + +dev.off() diff --git a/figures/cutoffpl_l+20.r b/figures/cutoffpl_l+20.r new file mode 100644 index 0000000..57358c0 --- /dev/null +++ b/figures/cutoffpl_l+20.r @@ -0,0 +1,99 @@ +# +# 2 data sets +# in Xspec make ipl->wdata +# then split file to four parts as lecr1/2/3/4 +# + +library(sfsmisc) + + +up=2.0 +#xgrid <- c(30,50,70,100); +Emin=25.0 +Emax=180.0 +postscript('cutoffpl_l+20.eps', horizontal = FALSE, onefile = FALSE, paper = "special",width = 9.0, height = 6.0) +par(oma=c(1,1,1,1)) +par(lwd=2) + +layout(matrix(c(1,2), 2, 1, byrow = TRUE), heights = c(1.5,1), TRUE) +leftx=5.3 +par(mar=c(0.1, leftx, 1, 1)) +par(cex.lab=1.6) +par(cex.axis=1.6) +xlim <- c(Emin,Emax) +ylim <- c(0.05,1.0) + +cl="black" +a <- read.table('./cutoffpl_l+20_eeuf.dat', col.names=c("x","dx","y","dy","total","cutoffpl","pow")) +plot(a$x, a$y, pch=3, bg="white", col="white", ylim=ylim, xlim=xlim,ylab=expression("keV"^"2"~"(Phot. keV"^"-1"~"cm"^"-2"~"s"^"-1"*")"),type="p",xaxt = 'n',xlab="",log="xy") + +#axis(side=2, at=c(1e-6,1e-5, 1e-4, 1e-3, 1e-2), labels=expression('-6','-5','-4','-3','-2')) + +#abline(v=xgrid, col="lightgray", lty="dotted") +grid() + +################################################################################### +## +## Plot model components +## +################################################################################### + +lines(a$x,a$cutoffpl,pch=1,bg=cl,col="blue",lwd=3,lty="longdash") +lines(a$x,a$pow,pch=1,bg=cl,col="red",lwd=3,lty="dotdash") +lines(a$x,a$total,pch=1,bg=cl,col="black") + + +# +# Plot upper limits +# +upx=a$x[(a$y/a$dy)=up] +py=a$y[(a$y/a$dy)>=up] +pdx=a$dx[(a$y/a$dy)>=up] +pdy=a$dy[(a$y/a$dy)>=up] + +segments(px,py-pdy,px,py+pdy,col=cl) +segments(px-pdx,py,px+pdx,py,col=cl) + +#legend( 35.0,1e-2,c("FPMA","FPMB"), lty=c(1,1), lwd=c(2.5,2.5),col=c("black","red"),cex=1.7) +text(100, 0.75, "L+20",cex=2.5) + + +################################################################################### + +par(mar=c(4.5, leftx, 0.1, 1)) +ylim <- c(-3.1,3.1) + +cl="black" +a <- read.table("./cutoffpl_l+20_delchi.dat", col.names=c("x","dx","y","dy")) +plot(a$x, a$y, pch=3, ylim=ylim, xlim=xlim, ylab=expression(Delta~chi), xlab="",type="p",log="x",xaxt = 'n') +grid() +segments(a$x,a$y-a$dy,a$x,a$y+a$dy,col=cl) +segments(a$x-a$dx,a$y,a$x+a$dx,a$y,col=cl) + +#abline(v=xgrid, col="lightgray", lty="dotted") +#grid() +#abline(h=0, col = "black",lty=2) +#dev.off() + +#abline(v=xgrid, col="lightgray", lty="dotted") +abline(h=0, col = "black",lty=2) + +mtext(side = 1, text = "Energy, keV", line = 4, cex=1.6) + +#axis(1, mgp=c(3, 1.5, 0)) +### Log axis using sfsmisc ### +atx=c(8,30,50,80,100,150) +eaxis(1, at = atx, labels = pretty10exp(atx, sub10=c(1,100), drop.1=TRUE), las=0) +abline(v=atx, col="lightgray", lty="dotted") +### + +dev.off() diff --git a/figures/cutoffpl_l-20.r b/figures/cutoffpl_l-20.r new file mode 100644 index 0000000..be90a46 --- /dev/null +++ b/figures/cutoffpl_l-20.r @@ -0,0 +1,99 @@ +# +# 2 data sets +# in Xspec make ipl->wdata +# then split file to four parts as lecr1/2/3/4 +# + + +library(sfsmisc) + + +up=2.0 +#xgrid <- c(30,50,70,100); +Emin=25.0 +Emax=180.0 +postscript('cutoffpl_l-20.eps', horizontal = FALSE, onefile = FALSE, paper = "special",width = 9.0, height = 6.0) +par(oma=c(1,1,1,1)) +par(lwd=2) + +layout(matrix(c(1,2), 2, 1, byrow = TRUE), heights = c(1.5,1), TRUE) +leftx=5.3 +par(mar=c(0.1, leftx, 1, 1)) +par(cex.lab=1.6) +par(cex.axis=1.6) +xlim <- c(Emin,Emax) +ylim <- c(0.05,1) + +cl="black" +a <- read.table('./cutoffpl_l-20_eeuf.dat', col.names=c("x","dx","y","dy","total","cutoffpl","pow")) +plot(a$x, a$y, pch=3, bg="white", col="white", ylim=ylim, xlim=xlim, ylab=expression("keV"^"2"~"(Phot. keV"^"-1"~"cm"^"-2"~"s"^"-1"*")"),type="p",xaxt = 'n',xlab="",log="xy") + +#axis(side=2, at=c(1e-6,1e-5, 1e-4, 1e-3, 1e-2), labels=expression('-6','-5','-4','-3','-2')) + +#abline(v=xgrid, col="lightgray", lty="dotted") +grid() +################################################################################### +## +## Plot model components +## +################################################################################### + +lines(a$x,a$cutoffpl,pch=1,bg=cl,col="blue",lwd=3,lty="longdash") +lines(a$x,a$pow,pch=1,bg=cl,col="red",lwd=3,lty="dotdash") +lines(a$x,a$total,pch=1,bg=cl,col="black") + +# +# Plot upper limits +# +upx=a$x[(a$y/a$dy)=up] +py=a$y[(a$y/a$dy)>=up] +pdx=a$dx[(a$y/a$dy)>=up] +pdy=a$dy[(a$y/a$dy)>=up] + +segments(px,py-pdy,px,py+pdy,col=cl) +segments(px-pdx,py,px+pdx,py,col=cl) + +#legend( 35.0,1e-2,c("FPMA","FPMB"), lty=c(1,1), lwd=c(2.5,2.5),col=c("black","red"),cex=1.7) +text(100, 0.75, "L-20",cex=2.5) + + +################################################################################### + +par(mar=c(4.5, leftx, 0.1, 1)) +ylim <- c(-3.1, 3.1) + +cl="black" +a <- read.table("./cutoffpl_l-20_delchi.dat", col.names=c("x","dx","y","dy")) +plot(a$x, a$y, pch=3, ylim=ylim, xlim=xlim, ylab=expression(Delta~chi), xlab="",type="p",log="x",xaxt = 'n') +grid() +segments(a$x,a$y-a$dy,a$x,a$y+a$dy,col=cl) +segments(a$x-a$dx,a$y,a$x+a$dx,a$y,col=cl) + +#abline(v=xgrid, col="lightgray", lty="dotted") +#grid() +#abline(h=0, col = "black",lty=2) +#dev.off() + +#abline(v=xgrid, col="lightgray", lty="dotted") + +abline(h=0, col = "black",lty=2) + +mtext(side = 1, text = "Energy, keV", line = 4, cex=1.6) + +#axis(1, mgp=c(3, 1.5, 0)) +### Log axis using sfsmisc ### +atx=c(8,30,50,80,100,150) +eaxis(1, at = atx, labels = pretty10exp(atx, sub10=c(1,100), drop.1=TRUE), las=0) +abline(v=atx, col="lightgray", lty="dotted") +### + +dev.off() diff --git a/figures/cutoffpl_scox1.r b/figures/cutoffpl_scox1.r new file mode 100644 index 0000000..15bc55e --- /dev/null +++ b/figures/cutoffpl_scox1.r @@ -0,0 +1,98 @@ +# +# 2 data sets +# in Xspec make ipl->wdata +# then split file to four parts as lecr1/2/3/4 +# + + +library(sfsmisc) + + +up=2.0 +#xgrid <- c(30,50,70,100); +Emin=25.0 +Emax=180.0 +postscript('cutoffpl_scox1.eps', horizontal = FALSE, onefile = FALSE, paper = "special",width = 9.0, height = 6.0) +par(oma=c(1,1,1,1)) +par(lwd=2) + +layout(matrix(c(1,2), 2, 1, byrow = TRUE), heights = c(1.5,1), TRUE) +leftx=5.3 +par(mar=c(0.1, leftx, 1, 1)) +par(cex.lab=1.6) +par(cex.axis=1.6) +xlim <- c(Emin,Emax) +ylim <- c(0.01,10) + +cl="black" +a <- read.table('./cutoffpl_scox1_eeuf.dat', col.names=c("x","dx","y","dy","total","cutoffpl","pow")) +plot(a$x, a$y, pch=3, bg="white", col="white", ylim=ylim, xlim=xlim, ylab=expression("keV"^"2"~"(Phot. keV"^"-1"~"cm"^"-2"~"s"^"-1"*")"),type="p",xaxt = 'n',xlab="",log="xy") + +#axis(side=2, at=c(1e-6,1e-5, 1e-4, 1e-3, 1e-2), labels=expression('-6','-5','-4','-3','-2')) + +#abline(v=xgrid, col="lightgray", lty="dotted") +grid() + +################################################################################### +## +## Plot model components +## +################################################################################### + +lines(a$x,a$cutoffpl,pch=1,bg=cl,col="blue",lwd=3,lty="longdash") +lines(a$x,a$pow,pch=1,bg=cl,col="red",lwd=3,lty="dotdash") +lines(a$x,a$total,pch=1,bg=cl,col="black") + + +# +# Plot upper limits +# +upx=a$x[(a$y/a$dy)=up] +py=a$y[(a$y/a$dy)>=up] +pdx=a$dx[(a$y/a$dy)>=up] +pdy=a$dy[(a$y/a$dy)>=up] + +segments(px,py-pdy,px,py+pdy,col=cl) +segments(px-pdx,py,px+pdx,py,col=cl) + +text(100, 1.5, "Sco X-1",cex=2.5) + +################################################################################### + +par(mar=c(4.5, leftx, 0.1, 1)) +ylim <- c(-3.1,3.1) + +cl="black" +a <- read.table("./cutoffpl_scox1_delchi.dat", col.names=c("x","dx","y","dy")) +plot(a$x, a$y, pch=3, ylim=ylim, xlim=xlim, ylab=expression(Delta~chi), xlab="",type="p",log="x",xaxt = 'n') +grid() +segments(a$x,a$y-a$dy,a$x,a$y+a$dy,col=cl) +segments(a$x-a$dx,a$y,a$x+a$dx,a$y,col=cl) + +#abline(v=xgrid, col="lightgray", lty="dotted") +#grid() +#abline(h=0, col = "black",lty=2) +#dev.off() + +#abline(v=xgrid, col="lightgray", lty="dotted") +abline(h=0, col = "black",lty=2) + +mtext(side = 1, text = "Energy, keV", line = 4, cex=1.6) + +#axis(1, mgp=c(3, 1.5, 0)) +### Log axis using sfsmisc ### +atx=c(8,30,50,80,100,150) +eaxis(1, at = atx, labels = pretty10exp(atx, sub10=c(1,100), drop.1=TRUE), las=0) +abline(v=atx, col="lightgray", lty="dotted") +### + +dev.off() diff --git a/figures/cutoffpl_solo_gb.r b/figures/cutoffpl_solo_gb.r new file mode 100644 index 0000000..cb88668 --- /dev/null +++ b/figures/cutoffpl_solo_gb.r @@ -0,0 +1,98 @@ +# +# 2 data sets +# in Xspec make ipl->wdata +# then split file to four parts as lecr1/2/3/4 +# + + +library(sfsmisc) + + +up=2.0 +#xgrid <- c(30,50,70,100); +Emin=25.0 +Emax=180.0 +postscript('cutoffpl_solo_gb.eps', horizontal = FALSE, onefile = FALSE, paper = "special",width = 9.0, height = 6.0) +par(oma=c(1,1,1,1)) +par(lwd=2) + +layout(matrix(c(1,2), 2, 1, byrow = TRUE), heights = c(1.5,1), TRUE) +leftx=5.3 +par(mar=c(0.1, leftx, 1, 1)) +par(cex.lab=1.6) +par(cex.axis=1.6) +xlim <- c(Emin,Emax) +ylim <- c(0.2,2) + +cl="black" +a <- read.table('./cutoffpl_solo_gb_eeuf.dat', col.names=c("x","dx","y","dy","cutoffpl")) +plot(a$x, a$y, pch=3, bg="white", col="white", ylim=ylim, xlim=xlim, ylab=expression("keV"^"2"~"(Phot. keV"^"-1"~"cm"^"-2"~"s"^"-1"*")"),type="p",xaxt = 'n',xlab="",log="xy") + +#axis(side=2, at=c(1e-6,1e-5, 1e-4, 1e-3, 1e-2), labels=expression('-6','-5','-4','-3','-2')) + +#abline(v=xgrid, col="lightgray", lty="dotted") +grid() + +################################################################################### +## +## Plot model components +## +################################################################################### + +lines(a$x,a$cutoffpl,pch=1,bg=cl,col="blue", lwd=3, lty="longdash") +#lines(a$x,a$pow,pch=1,bg=cl,col="red") +#lines(a$x,a$total,pch=1,bg=cl,col="black") + + +# +# Plot upper limits +# +upx=a$x[(a$y/a$dy)=up] +py=a$y[(a$y/a$dy)>=up] +pdx=a$dx[(a$y/a$dy)>=up] +pdy=a$dy[(a$y/a$dy)>=up] + +segments(px,py-pdy,px,py+pdy,col=cl) +segments(px-pdx,py,px+pdx,py,col=cl) + +text(100, 1.5, "GB",cex=2.5) + +################################################################################### + +par(mar=c(4.5, leftx, 0.1, 1)) +ylim <- c(-4.5,4.5) + +cl="black" +a <- read.table("./cutoffpl_solo_gb_delchi.dat", col.names=c("x","dx","y","dy")) +plot(a$x, a$y, pch=3, ylim=ylim, xlim=xlim, ylab=expression(Delta~chi), xlab="",type="p",log="x",xaxt = 'n') +grid() +segments(a$x,a$y-a$dy,a$x,a$y+a$dy,col=cl) +segments(a$x-a$dx,a$y,a$x+a$dx,a$y,col=cl) + +#abline(v=xgrid, col="lightgray", lty="dotted") +#grid() +#abline(h=0, col = "black",lty=2) +#dev.off() + +#abline(v=xgrid, col="lightgray", lty="dotted") +abline(h=0, col = "black",lty=2) + +mtext(side = 1, text = "Energy, keV", line = 4, cex=1.6) + +#axis(1, mgp=c(3, 1.5, 0)) +### Log axis using sfsmisc ### +atx=c(8,30,50,80,100,150) +eaxis(1, at = atx, labels = pretty10exp(atx, sub10=c(1,100), drop.1=TRUE), las=0) +abline(v=atx, col="lightgray", lty="dotted") +### + +dev.off() diff --git a/figures/galplane.r b/figures/galplane.r new file mode 100644 index 0000000..77f2037 --- /dev/null +++ b/figures/galplane.r @@ -0,0 +1,83 @@ +# +# 2 data sets +# in Xspec make ipl->wdata +# then split file +# + + + +library(sfsmisc) + + + + +up=-10.0 + +Emin=25.0 +Emax=180.0 +postscript('galplane.eps', horizontal = FALSE, onefile = FALSE, paper = "special",width = 9.0, height = 6.0) +par(oma=c(1,1,1,1)) +par(lwd=2) + +#layout(matrix(c(1,2), 2, 1, byrow = TRUE), heights = c(1.5,1), TRUE) + +leftx=4.5 +par(mar=c(4., leftx, 0.1, 0.1)) +par(cex.lab=1.6) +par(cex.axis=1.6) +xlim <- c(Emin,Emax) + +cl="black" + +sc=1000.0 + +ylim <- c(-0.005,0.09)*sc + +a <- read.table("./galplane_1.dat", col.names=c("x","dx","y","dy")) +b <- read.table("./galplane_3.dat", col.names=c("x","dx","y","dy")) +c <- read.table("./galplane_2.dat", col.names=c("x","dx","y","dy")) +d <- read.table("./galplane_4.dat", col.names=c("x","dx","y","dy")) + + +plot(a$x, a$y*sc, pch=0, ylim=ylim, xlim=xlim, ylab=expression("mCrab keV"^"-1"~"FOV"^"-1"), xaxt="no",log="x",xlab="",type="n") + +### Log axis using sfsmisc ### +atx=c(8,30,50,80,100,150) +eaxis(1, at = atx, labels = pretty10exp(atx, sub10=c(1,100), drop.1=TRUE), las=0) +abline(v=atx, col="lightgray", lty="dotted") +### + +grid() + +cex=1.5 + +cl="black" +segments(a$x,(a$y-a$dy)*sc,a$x,(a$y+a$dy)*sc,col=cl) +segments(a$x-a$dx,(a$y)*sc,a$x+a$dx,a$y*sc,col=cl) +points(a$x,a$y*sc,pch=21,col=cl,bg="white",cex=cex) + +# L-20 +cl="blue" +bias=1.01 +segments(c$x*bias,(c$y-c$dy)*sc,c$x*bias,(c$y+c$dy)*sc,col=cl) +segments((c$x-c$dx)*bias,c$y*sc,(c$x+c$dx)*bias,c$y*sc,col=cl) +points(c$x*bias,c$y*sc,pch=22,col=cl,bg="white",cex=cex) + +# L+20 +cl="red" +segments(b$x,(b$y-b$dy)*sc,b$x,(b$y+b$dy)*sc,col=cl) +segments(b$x-b$dx,b$y*sc,b$x+b$dx,b$y*sc,col=cl) +points(b$x,b$y*sc,pch=23,col=cl,bg="white",cex=cex) + +# BKG +cl="magenta" +segments(d$x,(d$y-d$dy)*sc,d$x,(d$y+d$dy)*sc,col=cl) +segments(d$x-d$dx,(d$y)*sc,d$x+d$dx,d$y*sc,col=cl) +points(d$x,d$y*sc,pch=24,col=cl,bg="white", ,cex=cex) + + +mtext(side = 1, text = "Energy, keV", line = 3, cex=1.6) + +legend( 85.0,0.09*sc,c("GB","L+20","L-20","3C 273/Coma"), pch=c(21,22,23,24),bg=c("white","white","white","white"), col=c("black","blue","red","magenta"),cex=1.7) + +dev.off() diff --git a/ridge/ridge/config.py b/ridge/ridge/config.py index ca08571..3e2b061 100644 --- a/ridge/ridge/config.py +++ b/ridge/ridge/config.py @@ -68,6 +68,10 @@ grxe_err_cut=90 crab_ra = 83.6287 crab_dec = 22.014 +""" Sco X-1 coordinates """ +sco_ra = 244.979455 +sco_dec = -15.640283 + bands={ 'A01':' 25.00 60.00', 'A02':' 60.00 80.00', diff --git a/ridge/ridge/utils.py b/ridge/ridge/utils.py index 32d47a8..72de1df 100644 --- a/ridge/ridge/utils.py +++ b/ridge/ridge/utils.py @@ -520,3 +520,176 @@ def get_spec(df, grxe_err_cut=None, skey=None, enkey=None, sigma=3, n_bins = 60, hdus[1].add_checksum() return sg_mean, sg_sem, skew_val, skew_err + + + + +def get_spec_src(df, grxe_err_cut=None, skey=None, enkey=None, sigma=3, n_bins = 60, plotme=False, gaussfit=False, fout=None, nscw_min=10, bootstrap=False): + src = np.array(df['SRC']) + + + rev = np.array(df['REV']) + phase = np.array(df['PHASE']) + clean = np.array(df['CLEAN']) + model = np.array(df['MODEL']) + crab = np.array(df['CRAB']) + resid = np.array(df['RESID']) + if(plotme): + print("get_spec_src: Initial size: {}".format(len(src))) + + """ + perc = np.percentile(src, grxe_err_cut, axis=0, keepdims=False) + print("{} {}: {}% cut of SCO: {:.2f} mCrab".format(skey,enkey,grxe_err_cut,perc)) + idx = np.where(src < perc) + + # in some cases grxe_err is just one number + if (len(src[idx])>nscw_min): + src=src[idx] + + if(plotme): + print("get_spec_scr: {} percentage error cut size: {}".format(grxe_err_cut, len(src))) + + clean=clean[idx] + model=model[idx] + crab=crab[idx] + rev=rev[idx] + phase=phase[idx] + resid=resid[idx] + + + idx = filtered_src.mask==False + src=src[idx] + + clean=clean[idx] + model=model[idx] + crab=crab[idx] + rev=rev[idx] + phase=phase[idx] + resid=resid[idx] + + + """ + filtered_data = sigma_clip(src, sigma=sigma, maxiters=10, return_bounds=True) + filtered_src = filtered_data[0] + filtered_min = filtered_data[1] + filtered_max = filtered_data[2] + + skew_val = skew(src) + skew_err = np.sqrt(6/len(src)) + print("{} SKEW: {:.5f} {:.5f}".format(enkey,skew_val, skew_err)) + + #print("length orig: {} taken: {} filtered: {}".format(len(src),len(src[filtered_src.mask==False]),len(src[filtered_src.mask==True]))) + + sg_mean, sg_med, sg_std = sigma_clipped_stats(src, sigma=sigma, maxiters=10) + sg_sem = sem(src) + + if(plotme): + print("{}: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f} N={}".format(enkey, sg_mean, sg_med, sg_std, sg_sem, len(src))) + + if(gaussfit==True): + k=1.1 + n, bins, patches = plt.hist(src, n_bins, density=True, facecolor='green', alpha=0.75, + range=[filtered_min*k, filtered_max*k]) + # add a 'best fit' line + area = np.sum(n * np.diff(bins)) + xdata = bins[:-1]+np.diff(bins)/2 + ydata = n + (mu, sg) = norm.fit(src) + + if(plotme==True): + print("Initial Gaiss fit: mu={:.2f} sigma={:.2f}".format(mu,sg)) + + y_peak = norm.pdf(mu, mu, sg) + plt.xlim(filtered_min*k, filtered_max*k) + plt.ylim(y_peak*area*0.001,y_peak*area*10) + fraction=0.1 + + params=[ + y_peak*area, # height + mu, # mu + sg, # sigma + 0.0, # const + 0.0, + ] + + #pfit, perr = fit_curvefit(params, xdata, ydata, const_gaussian_ff) + pfit, perr = fit_leastsq(params, xdata, ydata, const_gaussian_ff) + #print(pfit) + #print(perr) + + #pfit, perr = fit_leastsq(params, xdata, ydata, double_gaussian2_ff) + #pfit, perr = fit_curvefit(params, xdata, ydata, double_gaussian_ff) + #plt.plot(bins, double_gaussian_ff(bins, pfit), c='r' ) + + sg_mean = pfit[1] + sg_sem = perr[1] + if(plotme): + print("Const-Gauss fit: {:.2f} {:.2f}".format(sg_mean,sg_sem)) + + #sg_mean = mu + #sg_sem = sg/np.sqrt(len(grxe)) + + #sg_sem = perr[1] + #print("Narrow Gaiss fit: mu={:.2f} sigma={:.2f}".format(sg_mean,sg_sem)) + #print(" Wide Gaiss fit: mu={:.2f} sigma={:.2f}".format(pfit[4],pfit[5])) + + if(plotme): + plt.yscale("log") + k=1.1 + n, bins, patches = plt.hist(src, n_bins, density=True, facecolor='green', alpha=0.75, + range=[filtered_min*k, filtered_max*k]) + #plt.hist(grxe, bins=n_bins, range=[filtered_min*k, filtered_max*k]) + #plt.hist(grxe[filtered_grxe.mask], bins=n_bins, range=[filtered_min*k, filtered_max*k]) + plt.axvline(sg_mean, color="black") + plt.axvline(sg_mean+sg_sem, color="black", linestyle="dashed") + plt.axvline(sg_mean-sg_sem, color="black", linestyle="dashed") + plt.axvline(sg_mean+sg_std, color="blue", linestyle="dashed") + plt.axvline(sg_mean-sg_std, color="blue", linestyle="dashed") + plt.xlabel("Flux, mCrab") + plt.ylabel('Probability') + area = np.sum(n * np.diff(bins)) + xdata = bins[:-1]+np.diff(bins)/2 + ydata = n + #plt.plot(xdata, ydata, c='m' ) + + if(gaussfit==True): + # add a 'best fit' line + (mu, sg) = norm.fit(src) + print("Gauss fit: {:.2f} {:.2f}".format(mu,sg)) + y = norm.pdf(bins, mu, sg) + l = plt.plot(bins, y*area, 'r--', linewidth=2) + plt.plot(bins, const_gaussian_ff(bins, pfit), c='black' ) + #plt.plot(bins, double_gaussian2_ff(bins, pfit), c='black' ) + #plt.plot(bins, one_gaussian_ff(bins, pfit[:3]), c='b' ) + #plt.plot(bins, one_gaussian_ff(bins, [pfit[3],pfit[1],pfit[4]]), c='m' ) + + + plt.title("[DATA] {}: {}".format(skey, enkey)) + plt.grid(True) + plt.show() + + if(fout != None): + coldefs = fits.ColDefs([ + fits.Column(name='REV', format='J', unit='', array=rev), + fits.Column(name='CLEAN', format='D', unit='cts/s', array=clean), + fits.Column(name='PHASE', format='D', unit='', array=phase), + fits.Column(name='MODEL', format='D', unit='cts/s', array=model), + fits.Column(name='CRAB', format='D', unit='cts/s', array=crab), + fits.Column(name='RESID', format='D', unit='cts/s', array=resid), + fits.Column(name='SRC', format='D', unit='cts/s', array=src), + ]) + + + hdu = fits.BinTableHDU.from_columns(coldefs, name='SRC') + hdu.header['MISSION'] = ('INTEGRAL', '') + hdu.header['INSTITUT'] = ('IKI', 'Affiliation') + hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person') + hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail') + #hdu.add_checksum() + print(hdu.columns) + hdu.writeto(fout, overwrite=True) + + with fits.open(fout, mode='update') as hdus: + hdus[1].add_checksum() + + return sg_mean, sg_sem, skew_val, skew_err diff --git a/scripts/00_scox1_lightcurve.py b/scripts/00_scox1_lightcurve.py new file mode 100755 index 0000000..1492459 --- /dev/null +++ b/scripts/00_scox1_lightcurve.py @@ -0,0 +1,209 @@ +#!/usr/bin/env python + +__author__ = "Roman Krivonos" +__copyright__ = "Space Research Institute (IKI)" + +import numpy as np +import pandas as pd +from astropy.io import fits +from astropy import wcs +import matplotlib.pyplot as plt +import math, sys, os +import pickle +from numpy.polynomial import Polynomial +from astropy.table import Table, Column + +from sklearn.linear_model import LinearRegression +from sklearn.linear_model import HuberRegressor +from sklearn.linear_model import RANSACRegressor +from sklearn.linear_model import TheilSenRegressor +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import RepeatedKFold + + +from astropy.coordinates import SkyCoord # High-level coordinates +from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames +from astropy.coordinates import Angle, Latitude, Longitude # Angles +import astropy.units as u + +#from statsmodels.robust.scale import huber +from astropy.stats import sigma_clip +from numpy import absolute +from numpy import arange + +from ridge.utils import * +from ridge.config import * + +if not os.path.exists(proddir): + os.makedirs(proddir) + +#enkey = sys.argv[1] +enkey='A01' +inkey="ALL" + +fn="detcnts.{}.{}.resid.fits".format(enkey,inkey) +dat = Table.read(proddir+fn, unit_parse_strict='silent') +df = dat.to_pandas() + +with open(proddir+"detcnts.{}.ignored_scw.pkl".format(enkey), 'rb') as fp: + ignored_scw = pickle.load(fp) + +sco_crd = SkyCoord(sco_ra, sco_dec, frame=FK5(), unit="deg") + +plotme=False + +npix = 50 +sw = 30.0 # deg +pix = sw/npix # pixel size in degrees + + +crabmodel={} +rota_arr=[] +a0=[] +b0=[] +a_full0=[] +b_full0=[] +b_est0=[] +err0=[] +rev0=[] + +totx=[] +toty=[] + +for i,rec in df.iterrows(): + obsid = rec['OBSID']#.decode("utf-8") + + if (obsid in ignored_scw): + print("Skip ScW",obsid) + continue + +# accumulate full data set +for rev in range(revmin,revmax): + + + df0 = df.query('SRC > 0.0 & REV == {} & PHASE > {} & PHASE < {} & SCO_SEP < {}'.format(rev,phmin,phmax,crab_sep_max)) + nobs=len(df0) + if not (nobs> crab_nmax): + continue + + print(rev,nobs) + + for n in df0['SCO_SEP'].values: + totx.append(n) + + for n in df0['SRC'].values: + toty.append(n) + + +x = np.array(totx) +y = np.array(toty) +x = x.reshape((-1, 1)) +model = LinearRegression() +#model = TheilSenRegressor() +results = evaluate_model(x, y, model) +a_full,b_full,err_full = plot_best_fit(x, y, model) +if(plotme): + plot_ab(x, y, a_full, b_full, err_full, title="REGRESSION") + + + +# go over orbits +poly_x=[] +poly_y=[] + +ntotal=0 +nrev=0 +for rev in range(revmin,revmax): + df0 = df.query('SRC > 0.0 & REV == {} & PHASE > {} & PHASE < {} & SCO_SEP < {}'.format(rev,phmin,phmax,crab_sep_max)) + nobs=len(df0) + if not (nobs> crab_nmax): + continue + #cen_ra = np.array(df0['RA'].values) + #cen_dec = np.array(df0['DEC'].values) + print("*** Orbit ",rev) + x = np.array(df0['SCO_SEP'].values) + y = np.array(df0['SRC'].values) + + #rota_deg = np.array(df0['ROTA'].values) + #rota = np.array(df0['ROTA'].values) * np.pi / 180. # in radians + #detx = np.cos(rota)*x + #dety = np.sin(rota)*x + + x = x.reshape((-1, 1)) + model = LinearRegression() + #model = TheilSenRegressor() + results = evaluate_model(x, y, model) + a,b,err = plot_best_fit(x, y, model) + + b_est = np.mean(y - a_full*x) + + a_full0.append(a_full) + b_full0.append(b_full) + b_est0.append(b_est) + + + a0.append(a) + b0.append(b) + err0.append(err) + rev0.append(rev) + + poly_x.append(rev) + poly_y.append(b_est) + + + crabmodel[rev]={'a':a_full, 'b':b_est, 'err':err} + if(plotme): + plot_ab(x, y, a_full, b_est, err, title="REGRESSION rev {}".format(rev)) + + + + print("ax+b: a={:.2e}, b={:.2e}, std={:.2e}\n".format(a,b,err)) + ntotal=ntotal+nobs + nrev=nrev+1 + +print("Orbits: {}, obs: {}".format(nrev,ntotal)) + +npoly=4 + +z = np.polyfit(poly_x, poly_y, npoly) + +p = np.poly1d(z) +poly_z=[] +for t in poly_x: + poly_z.append(p(t)) +plt.scatter(poly_x, poly_y) +plt.plot(poly_x, poly_z, color='r') +plt.title("Crab detector count rate evolution") +plt.ylabel("Crab count rate cts/s/pix") +plt.xlabel("INTEGRAL orbit") +#plt.show() +plt.savefig(proddir+fn.replace(".fits",".sco_rate.png")) + +indices = sorted( + range(len(rev0)), + key=lambda index: rev0[index] +) + +coldefs = fits.ColDefs([ + fits.Column(name='REV', format='J', unit='', array=[rev0[index] for index in indices]), + fits.Column(name='A', format='D', unit='', array=[a0[index] for index in indices]), + fits.Column(name='B', format='D', unit='', array=[b0[index] for index in indices]), + fits.Column(name='ERR', format='D', unit='', array=[err0[index] for index in indices]), + fits.Column(name='A_FULL', format='D', unit='', array=[a_full0[index] for index in indices]), + fits.Column(name='B_FULL', format='D', unit='', array=[b_full0[index] for index in indices]), + fits.Column(name='B_EST', format='D', unit='', array=[b_est0[index] for index in indices]), + fits.Column(name='B_POLY', format='D', unit='', array=[poly_z[index] for index in indices]), +]) + +fout = fn.replace(".fits",".scox1.fits") +hdu = fits.BinTableHDU.from_columns(coldefs, name='GRXE') +hdu.header['MISSION'] = ('INTEGRAL', '') +hdu.header['TELESCOP'] = ('IBIS', '') +hdu.header['INSTITUT'] = ('IKI', 'Affiliation') +hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person') +hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail') +#hdu.add_checksum() +print(hdu.columns) +hdu.writeto(proddir+fout, overwrite=True) + + diff --git a/scripts/02_grxe_resid.py b/scripts/02_grxe_resid.py index 866b319..56c8651 100755 --- a/scripts/02_grxe_resid.py +++ b/scripts/02_grxe_resid.py @@ -36,6 +36,8 @@ from numpy import arange from ridge.utils import * from ridge.config import * +sco_crd = SkyCoord(sco_ra, sco_dec, frame=FK5(), unit="deg") + plotme=False enkey = sys.argv[1] outkey = sys.argv[2] @@ -71,6 +73,8 @@ clean = data.field('clean') phase = data.field('phase') texp = data.field('exposure') +src = data.field('src') # for Sco X-1 testing + obsid0=[] rev0=[] phase0=[] @@ -90,8 +94,10 @@ lat0=[] base0=[] c0=[] texp0=[] +src0=[] +sco_sep0=[] -hdulist = fits.open(datadir+modelrxte) +hdulist = fits.open(modelsdir+modelrxte) w = wcs.WCS(hdulist[0].header) smap = hdulist[0].data sx=int(hdulist[0].header['NAXIS1']) @@ -191,6 +197,12 @@ for i, row in df.iterrows(): rev0.append(orbit) lon0.append(row['LON']) lat0.append(row['LAT']) + src0.append(1000*(float(row['SRC'])/p(orbit))) + + ra=float(row['RA']) + dec=float(row['DEC']) + sc = SkyCoord(ra, dec, frame=FK5(), unit="deg") + sco_sep0.append(sco_crd.separation(sc).deg) lon=row['LON'] lat=row['LAT'] @@ -276,6 +288,7 @@ coldefs = fits.ColDefs([ fits.Column(name='TEXP', format='D', unit='', array=[texp0[index] for index in indices]), fits.Column(name='PHASE', format='D', unit='', array=[phase0[index] for index in indices]), fits.Column(name='CLEAN', format='D', unit='cts/s', array=[clean0[index] for index in indices]), + fits.Column(name='SRC', format='D', unit='cts/s', array=[src0[index] for index in indices]), fits.Column(name='MODEL', format='D', unit='cts/s', array=[model0[index] for index in indices]), fits.Column(name='MODEL_ERR', format='D', unit='', array=[model_err0[index] for index in indices]), fits.Column(name='RESID', format='D', unit='cts/s', array=[resid0[index] for index in indices]), @@ -288,6 +301,7 @@ coldefs = fits.ColDefs([ fits.Column(name='C', format='D', unit='', array=[c0[index] for index in indices]), fits.Column(name='DS9X', format='D', unit='', array=[ds9x[index] for index in indices]), fits.Column(name='DS9Y', format='D', unit='', array=[ds9y[index] for index in indices]), + fits.Column(name='SCO_SEP', format='D', unit='', array=[sco_sep0[index] for index in indices]), ]) fout = fn.replace(".fits",".{}.resid.fits".format(outkey)) diff --git a/scripts/02_grxe_resid.sh b/scripts/02_grxe_resid.sh index 0e10147..2f0e7e0 100755 --- a/scripts/02_grxe_resid.sh +++ b/scripts/02_grxe_resid.sh @@ -1,14 +1,14 @@ for band in A01 A02 A03 do ./02_grxe_resid.py $band ALL -./02_grxe_resid.py $band BKG -./02_grxe_resid.py $band GAL +#./02_grxe_resid.py $band BKG +#./02_grxe_resid.py $band GAL done for band in B01 B02 B03 B04 B05 B06 B07 B08 B09 B10 B11 B12 B13 B14 B15 B16 B17 B18 B19 B20 B21 do ./02_grxe_resid.py $band ALL -./02_grxe_resid.py $band BKG -./02_grxe_resid.py $band GAL +#./02_grxe_resid.py $band BKG +#./02_grxe_resid.py $band GAL done diff --git a/scripts/03_scox1_spec.py b/scripts/03_scox1_spec.py new file mode 100755 index 0000000..c0624b8 --- /dev/null +++ b/scripts/03_scox1_spec.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python + +__author__ = "Roman Krivonos" +__copyright__ = "Space Research Institute (IKI)" + +import numpy as np +import pandas as pd +from astropy.io import fits +from astropy.table import Table, Column +from astropy import units as u + +import matplotlib.pyplot as plt +import math, sys, os +import pickle + +from sklearn.linear_model import LinearRegression +from sklearn.linear_model import HuberRegressor +from sklearn.linear_model import RANSACRegressor +from sklearn.linear_model import TheilSenRegressor +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import RepeatedKFold + +#from statsmodels.robust.scale import huber +from astropy.stats import sigma_clip +from astropy.stats import sigma_clipped_stats + +from scipy.stats import norm +from scipy.stats import describe +from scipy.stats import sem +import subprocess + +from numpy import absolute +from numpy import arange + +from ridge.utils import * +from ridge.config import * + +inkey="ALL" + +skey="SCOX1" + + +plotme=False + + +ebands0={'B01':[0.0,0.0], + 'B02':[0.0,0.0], + 'B03':[0.0,0.0], + 'B04':[0.0,0.0], + 'B05':[0.0,0.0], + 'B06':[0.0,0.0], + 'B07':[0.0,0.0], + 'B08':[0.0,0.0], + 'B09':[0.0,0.0], + 'B10':[0.0,0.0], + 'B11':[0.0,0.0], + 'B12':[0.0,0.0], + 'B13':[0.0,0.0], + 'B14':[0.0,0.0], + 'B15':[0.0,0.0], + 'B16':[0.0,0.0], + 'B17':[0.0,0.0], + 'B18':[0.0,0.0], + 'B19':[0.0,0.0], + 'B20':[0.0,0.0], + 'B21':[0.0,0.0], + } + +""" +ebands0={'A01':[0.0,0.0], + 'A02':[0.0,0.0], + 'A03':[0.0,0.0], + } +""" + +mcrab=u.def_unit('mCrab') +ctss=u.def_unit('cts/s') +u.add_enabled_units([mcrab,ctss]) + +if not os.path.exists(specdir): + os.makedirs(specdir) + +fitsdir = "{}fits/".format(specdir) +if not os.path.exists(fitsdir): + os.makedirs(fitsdir) + +for enkey in ebands0.keys(): + + fn="detcnts.{}.{}.resid.fits".format(enkey,inkey) + dat = Table.read(proddir+fn, unit_parse_strict='silent') + df = dat.to_pandas() + + #crab_sep_max=2.0 + query = 'PHASE > {} & PHASE < {} & SCO_SEP < {}'.format(phmin,phmax,crab_sep_max) + df = df.query(query) + print("{}: {} N={}".format(enkey, query, df.shape[0])) + + t = Table.from_pandas(df) + t.write("{}fits/SCOX1.{}.fits".format(specdir,enkey),overwrite=True) + + if not (df.shape[0]>0): + print("continue") + continue + + sg_mean,sg_sem,skew_val,skew_err = get_spec_src(df, sigma=3, grxe_err_cut=grxe_err_cut, enkey=enkey, plotme=True, gaussfit=True) + ebands0[enkey]=[sg_mean,sg_sem] + + + + +fspec="{}{}.spec".format(specdir,skey) +with open(fspec, 'w') as fp: + for enkey in ebands0.keys(): + flux=ebands0[enkey][0] + err=ebands0[enkey][1] + print("[DATA] {}: {} {:.6f} {:.6f}".format(skey,enkey,flux,err)) + fp.write("0 {} {:.6f} {:.6f} 0.0\n".format(bands[enkey],flux,err)) +subprocess.run(["perl", "do_pha_igr_v3_mCrab.pl", fspec]) + + +try: + for remfile in ["cols","cols1","cols2","header",]: + os.remove(remfile) +except OSError: + pass diff --git a/scripts/_00_scan_detcnts.py b/scripts/_00_scan_detcnts.py index 6b3af0b..e494a49 100644 --- a/scripts/_00_scan_detcnts.py +++ b/scripts/_00_scan_detcnts.py @@ -3,7 +3,7 @@ __author__ = "Roman Krivonos" __copyright__ = "Space Research Institute (IKI)" -""" Работает только в архиве """ +""" TO BE RUN IN ARCHIVE ONLY """ import re import sys