#!/bin/env python

import string
import os
import subprocess
import sys
import shutil

from pylab import *

from astropy.io import fits


filename_sepe = input("Enter the SE+PE filename: ")

filename_pe   = input("Enter the PE filename: ")

filename_rmf  = input("Enter the spectral response filename: ")

ext_pos = input("Enter extension number: ")

print("---------------")
print("SE+PE:     ",filename_sepe)
print("PE:        ",filename_pe)
print("response:  ",filename_rmf)
print("extension: ",ext_pos)
print("---------------")

## energy threshold: 400 keV after 2012 May; 650 keV before.
energy_th = 400.

## PSD efficiency:
pe_eff = 0.85


## read filename_pe and correct for PE efficiency:
ftemp = fits.open(filename_pe)
canaletemp_pe = ftemp[ext_pos].data['channel']
ratetemp_pe = ftemp[ext_pos].data['rate']/pe_eff
eratetemp_pe = ftemp[ext_pos].data['stat_err']/pe_eff
ftemp.close()


## Search channel point where filename_sepe and filename_pe must be glued:
frmf = fits.open(filename_rmf)
canale = frmf[3].data['channel']
e_min = frmf[3].data['E_MIN']
e_max = frmf[3].data['E_MAX']
frmf.close()

canale_th=0
k=0
while (e_min[k] < energy_th):
    canale_th = canale_th +1
    k = k +1 


shutil.copy(filename_sepe,'spectrum.fits')
##shutil.copy(filename_rmf,'.')

i=0
fres = fits.open('spectrum.fits', mode='update')
for i in range(len(canale)):
    if fres[ext_pos].data['channel'][i] >= canale_th :
        fres[ext_pos].data['rate'][i] = ratetemp_pe[i]
        fres[ext_pos].data['stat_err'][i] = eratetemp_pe[i]


fres.flush()
fres.close()