#
# 2 data sets
# in Xspec make ipl->wdata
# then split file to four parts as lecr1/2/3/4
#

key="lecr"
#key="pow"
keymo="model2"


up=-10.0
#xgrid <- c(30,50,70,100); 
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(5., leftx, 0.1, 0.1))
par(cex.lab=1.6)
par(cex.axis=1.6)
xlim <- c(Emin,Emax)

cl="black"


ylim <- c(-0.005,0.09)

a <- read.table("../galplane_bkg.dat", col.names=c("x","dx","y","dy"))
b <- read.table("../galplane_gc06.dat", col.names=c("x","dx","y","dy"))
c <- read.table("../galplane+20.dat", col.names=c("x","dx","y","dy"))
d <- read.table("../galplane-20.dat", col.names=c("x","dx","y","dy"))


plot(a$x, a$y, pch=3, ylim=ylim, xlim=xlim, ylab=expression("Crab keV"^"-1"~"FOV"^"-1"), xlab="",type="p",log="x",xaxt = 'n')

grid()

cl="red"
segments(b$x,b$y-b$dy,b$x,b$y+b$dy,col=cl)
segments(b$x-b$dx,b$y,b$x+b$dx,b$y,col=cl)

cl="blue"
bias=1.01
segments(c$x*bias,c$y-c$dy,c$x*bias,c$y+c$dy,col=cl)
segments((c$x-c$dx)*bias,c$y,(c$x+c$dx)*bias,c$y,col=cl)

cl="green"
segments(d$x,d$y-d$dy,d$x,d$y+d$dy,col=cl)
segments(d$x-d$dx,d$y,d$x+d$dx,d$y,col=cl)

# background
cl="black"
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)


mtext(side = 1, text = "Energy, keV", line = 4, cex=1.6)

legend( 85.0,0.08,c("GC","LON+20","LON-20","3C 273/Coma"), lty=c(1,1,1,1), lwd=c(2.5,2.5,2.5,2.5),col=c("red","blue","green","black"),cex=1.7)


axis(1, mgp=c(3, 1.5, 0))

dev.off()