ctools

import numpy

import os, sys

import random

import matplotlib.pyplot as plt

from matplotlib import rc

import gammalib

import ctools

import cscripts

data=numpy.genfromtxt('pl_4lac_with_redshift_v2.dat', dtype='string')

Source_Name=[]

RAJ2000=[]

DECJ2000=[]

Redshift=[]

Pivot_Energy=[]

PL_Flux_Density=[]#units ph / (cm2 MeV s)

PL_Index=[]

for i in range(1,len(data)):

    if(0<float(data[i][3])<4):

        Source_Name.append(data[i][0])

        RAJ2000.append(float(data[i][1]))

        DECJ2000.append(float(data[i][2]))

        Redshift.append(float(data[i][3]))

        Pivot_Energy.append(float(data[i][4]))

        PL_Flux_Density.append(float(data[i][5]))

        PL_Index.append(float(data[i][6]))

caldb = 'prod3b-v2'

tmin = 0.0

tmax= 72000

emin  = 1.0   # TeV still need to do 0.03 and 0.05

emax  = 199 # TeV

irf_used=[]

decs=[]

cta_north_lat=28.71

cta_south_lat=-24.63

alt_cul_north=[]

alt_cul_south=[]

north_irf=[]

south_irf=[]

for i in range(len(DECJ2000)):

    ra=RAJ2000[i]

    dec=DECJ2000[i]

    decs.append(dec)

    source=Source_Name[i]

    alt_cul_north.append(90-numpy.abs(cta_north_lat-dec))

    alt_cul_south.append(90-numpy.abs(cta_south_lat-dec))

#print numpy.min(alt_cul_north)

for i in range(len(decs)):

    if(alt_cul_north[i]<25):

        north_irf.append('North_z20_50h')

    if(25<alt_cul_north[i]<45):

        north_irf.append('North_z40_50h')

    if(45<alt_cul_north[i]<65):

        north_irf.append('North_z60_50h')

    if(65<alt_cul_north[i]<90):

        north_irf.append('North_z60_50h')

    if(alt_cul_south[i]<25):

        south_irf.append('South_z20_50h')

    if(25<alt_cul_south[i]<45):

        south_irf.append('South_z40_50h')

    if(45<alt_cul_south[i]<65):

        south_irf.append('South_z60_50h')

    if(65<alt_cul_south[i]<90):

        south_irf.append('South_z60_50h')

#print len(decs), len(north_irf),len(south_irf)

#y_order=['North_z20_50h','North_z40_50h','North_z60_50h']

'''

plt.figure(figsize=(16,12))

plt.scatter(north_irf,decs)

plt.xticks(fontsize=14, rotation=0)

plt.yticks(fontsize=14, rotation=0)

plt.ylabel('Declination')

plt.xlabel('North IRF used')

plt.savefig('north_irf.png',dpi=100)

plt.figure(figsize=(16,12))

plt.scatter(south_irf,decs)

plt.ylabel('Declination')

plt.xticks(fontsize=14, rotation=0)

plt.yticks(fontsize=14, rotation=0)

plt.xlabel('South IRF used')

plt.savefig('south_irf.png',dpi=100)

for i in range(125,len(Source_Name)):

    if(Source_Name[i]=='4FGL_J0157.7-4614'):

        print i

'''

for i in range(len(Source_Name)):

    irf_n=north_irf[i]

    irf_s=south_irf[i]

    irf_n_30='North_z20_50h'

    irf_s_30='South_z20_50h'

    ra=RAJ2000[i]

    dec=DECJ2000[i]

    source=Source_Name[i]

    print source,ra,dec

    for j in range(1):

        #s1=random.randint(1,1000)

        evfile = '/Volumes/mac_ext/egal_pop/PL_events_file_50hr/%s_events_n1000GeV_%i.fits'%(source,j)

        obssim = ctools.ctobssim()

        obssim['ra']        = ra

        obssim['dec']       = dec

        obssim['rad']       = 3.0

        obssim['tmin']      = tmin#'2020-01-01T00:00:00'

        obssim['tmax']      = tmax

        #obssim['seed'] = s1 # Any integer value you like

        obssim['emin']      = emin

        obssim['emax']      = emax

        obssim['caldb']     = caldb

        obssim['irf']       = irf_n

        obssim['nthreads']=0

        obssim['chatter']=2

        #obssim['clobber']='yes'

        obssim['mode']='ql'

        obssim['logfile']='ctobssim_n_%s.log'%(source)

        obssim['inmodel']   = '/Volumes/mac_ext/egal_pop/PL_input_model/input_loglaw_pl_%s.xml'%(source)

        obssim['outevents'] = evfile

        obssim.execute()

    for j in range(1):

        #s1=random.randint(1,1000)

        evfile = '/Volumes/mac_ext/egal_pop/PL_events_file_50hr/%s_events_s1000GeV_%i.fits'%(source,j)

        obssim = ctools.ctobssim()

        obssim['ra']        = ra

        obssim['dec']       = dec

        obssim['rad']       = 3.0

        obssim['tmin']      = tmin#'2020-01-01T00:00:00'

        obssim['tmax']      = tmax

        #obssim['seed'] = s1 # Any integer value you like

        obssim['emin']      = emin

        obssim['emax']      = emax

        obssim['caldb']     = caldb

        obssim['irf']       = irf_s

        obssim['inmodel']   = '/Volumes/mac_ext/egal_pop/PL_input_model/input_loglaw_pl_%s.xml'%(source)

        obssim['nthreads']=0

        obssim['chatter']=2

        #obssim['clobber']='yes'

        obssim['mode']='ql'

        obssim['logfile']='ctobssim_s_%s.log'%(source)

        obssim['outevents'] = evfile

        obssim.execute()

        like = ctools.ctlike()

        evfile = '/Volumes/mac_ext/egal_pop/PL_events_file_50hr/%s_events_n1000GeV_%i.fits'%(source,j)

        like['inobs']   = evfile

        like['caldb']   = caldb

        like['irf']     = irf_n

        like['inmodel'] = '/Volumes/mac_ext/egal_pop/PL_input_model/input_loglaw_pl_%s.xml'%(source)

        like['outcovmat']='NONE'

        like['statistic']='DEFAULT'

        like['like_accuracy']= 0.005

        like['max_iter']=50

        like['logfile']='ctlike_n_%s.log'%(source)

        like['outmodel'] = '/Volumes/mac_ext/egal_pop/Fit_results_50hr/results_n1000_%s_%s.xml'%(source,j)

        like.execute()

        import re

        import sys

        import os

        import glob

        orig_stdout = sys.stdout

        f = open('/Volumes/mac_ext/egal_pop/Fit_results_50hr/results_n1000_%s_%s.txt'%(source,j),'w')

        sys.stdout = f

        print(like.obs().models()[0])

        #print s1

        sys.stdout = orig_stdout

        f.close()

        like = ctools.ctlike()

        evfile = '/Volumes/mac_ext/egal_pop/PL_events_file_50hr/%s_events_s1000GeV_%i.fits'%(source,j)

        like['inobs']   = evfile

        like['caldb']   = caldb

        like['irf']     = irf_s

        like['inmodel'] = '/Volumes/mac_ext/egal_pop/PL_input_model/input_loglaw_pl_%s.xml'%(source)

        like['outcovmat']='NONE'

        like['statistic']='DEFAULT'

        like['like_accuracy']= 0.005

        like['max_iter']=50

        like['logfile']='ctlike_s_%s.log'%(source)

        like['outmodel'] = '/Volumes/mac_ext/egal_pop/Fit_results_50hr/results_s1000_%s_%s.xml'%(source,j)

        like.execute()

        orig_stdout = sys.stdout

        f = open('/Volumes/mac_ext/egal_pop/Fit_results_50hr/results_s1000_%s_%s.txt'%(source,j),'w')

        sys.stdout = f

        print(like.obs().models()[0])

        #print s1

        sys.stdout = orig_stdout

        f.close()

'''

        butterfly = '/Volumes/mac_ext/egal_pop/PL_obtained_spectra/butterfly_%s_n30GeV.txt'%(source) # text file where the butterfly plot is stored

        bttflmk = ctools.ctbutterfly()

        bttflmk['inobs'] = '/Volumes/mac_ext/egal_pop/PL_events_file/%s_events_n_%i.fits'%(source,j)

        bttflmk['inmodel'] = '/Volumes/mac_ext/egal_pop/Fit_results/results_n30_%s_%s.xml'%(source,j)

        bttflmk['srcname'] = '%s'%(source)

        bttflmk['emin'] = emin

        bttflmk['emax'] = emax

        bttflmk['caldb']     = caldb

        bttflmk['irf']       = irf_n_30

        #bttflmk['statistic'] = 'WSTAT'

        bttflmk['outfile'] = butterfly

        bttflmk.execute()

        butterfly = '/Volumes/mac_ext/egal_pop/PL_obtained_spectra/butterfly_%s_s30GeV.txt'%(source) # text file where the butterfly plot is stored

        bttflmk = ctools.ctbutterfly()

        bttflmk['inobs'] = '/Volumes/mac_ext/egal_pop/PL_events_file/%s_events_s_%i.fits'%(source,j)

        bttflmk['inmodel'] = '/Volumes/mac_ext/egal_pop/Fit_results/results_s30_%s_%s.xml'%(source,j)

        bttflmk['srcname'] = '%s'%(source)

        bttflmk['emin'] = emin

        bttflmk['emax'] = emax

        bttflmk['caldb']     = caldb

        bttflmk['irf']       = irf_s_30

        #bttflmk['statistic'] = 'WSTAT'

        bttflmk['outfile'] = butterfly

        bttflmk.execute()

        sed = "/Volumes/mac_ext/egal_pop/PL_obtained_spectra/sed%s_n30GeV.fits"%(source)

        sedmkr = cscripts.csspec()

        sedmkr['inobs'] = '/Volumes/mac_ext/egal_pop/PL_events_file/%s_events_n_%i.fits'%(source,j)

        sedmkr['inmodel'] = '/Volumes/mac_ext/egal_pop/Fit_results/results_n30_%s_%s.xml'%(source,j)

        sedmkr['srcname'] = '%s'%(source)

        sedmkr['outfile'] = sed

              #sedmkr['statistic'] = 'WSTAT'

        sedmkr['method'] = "AUTO"

        sedmkr['emin'] = emin

        sedmkr['emax'] = emax

        sedmkr['caldb']     = caldb

        sedmkr['irf']       = irf_n_30

        sedmkr['enumbins'] = 10

        sedmkr['ebinalg'] = 'LOG'

        sedmkr.execute()

        sed = "/Volumes/mac_ext/egal_pop/PL_obtained_spectra/sed%s_s30GeV.fits"%(source)

        sedmkr = cscripts.csspec()

        sedmkr['inobs'] = '/Volumes/mac_ext/egal_pop/PL_events_file/%s_events_s_%i.fits'%(source,j)

        sedmkr['inmodel'] = '/Volumes/mac_ext/egal_pop/Fit_results/results_s30_%s_%s.xml'%(source,j)

        sedmkr['srcname'] = '%s'%(source)

        sedmkr['outfile'] = sed

        #sedmkr['statistic'] = 'WSTAT'

        sedmkr['method'] = "AUTO"

        sedmkr['emin'] = emin

        sedmkr['emax'] = emax

        sedmkr['caldb']     = caldb

        sedmkr['irf']       = irf_s_30

        sedmkr['enumbins'] = 10

        sedmkr['ebinalg'] = 'LOG'

        sedmkr.execute()

#print (like.obs().models()[0])

    text_file = open('results_%s_0.03_%s.txt'%(source_name,s))

#print text_file

    lines = text_file.readlines()

    #print lines

    if (numpy.float(re.findall("\d+\.\d+", lines[3])[0])>10):

        ts.append(numpy.float(re.findall("\d+\.\d+", lines[3])[0]))

print ts

plt.hist(ts, bins=10)

plt.xlabel('TS')

plt.ylabel('Number')

plt.xlim(numpy.min(ts),numpy.max(ts)+10)

plt.savefig('ts_hist.png')

#plt.show()

maxt=0

for t in range(len(ts)):

    if (ts[t]==numpy.max(ts)):

        maxt= t

s=maxt

def plot_butterfly(filename,ax,color,label):

    #this snippet is extracted from

    #$CTOOLS/examples/show_butterfly.py

    # Read butterfly file

        csv = gammalib.GCsv(filename)

    # Initialise arrays to be filled

        butterfly_x = []

        butterfly_y = []

        line_x      = []

        line_y      = []

    # Loop over rows of the file

        nrows = csv.nrows()

        for row in range(nrows):

        # Get conversion coefficient

            conv = csv.real(row,0) * csv.real(row,0) * gammalib.MeV2erg

        # Compute upper edge of confidence band

            butterfly_x.append(csv.real(row,0)/1.0e6) # TeV

            butterfly_y.append(csv.real(row,2)*conv)

        # Set line values

            line_x.append(csv.real(row,0)/1.0e6) # TeV

            line_y.append(csv.real(row,1)*conv)

    # Loop over the rows backwards to compute the lower edge of the

    # confidence band

        for row in range(nrows):

            index = nrows - 1 - row

            conv  = csv.real(index,0) * csv.real(index,0) * gammalib.MeV2erg

            butterfly_x.append(csv.real(index,0)/1.0e6)

            low_error = max(csv.real(index,3)*conv, 1e-26)

            butterfly_y.append(low_error)

    # Plot

        ax.plot(line_x,line_y,color=color,ls='--',alpha=0.5,label=label)

            #ax.fill(butterfly_x,butterfly_y,color=color,alpha=0.2)

def plot_sed(filename,ax,color,ts_thresh=9):

    #this snippet is extracted from

    #$CTOOLS/examples/show_spectrum.py

    fits     = gammalib.GFits(filename)

    table    = fits.table(1)

    c_energy = table['Energy']

    c_ed     = table['ed_Energy']

    c_eu     = table['eu_Energy']

    c_flux   = table['Flux']

    c_eflux  = table['e_Flux']

    c_ts     = table['TS']

    c_upper  = table['UpperLimit']

    # Initialise arrays to be filled

    energies    = []

    flux        = []

    ed_engs     = []

    eu_engs     = []

    e_flux      = []

    ul_energies = []

    ul_ed_engs  = []

    ul_eu_engs  = []

    ul_flux     = []

    # Loop over rows of the file

    nrows = table.nrows()

    for row in range(nrows):

        # Get Test Statistic, flux and flux error

        ts    = c_ts.real(row)

        flx   = c_flux.real(row)

        e_flx = c_eflux.real(row)

        # If Test Statistic is larger than threshold and flux error is smaller than

        # flux then append flux plots ...

        if ts > ts_thresh and e_flx < flx:

            energies.append(c_energy.real(row))

            flux.append(c_flux.real(row))

            ed_engs.append(c_ed.real(row))

            eu_engs.append(c_eu.real(row))

            e_flux.append(c_eflux.real(row))

        # ... otherwise append upper limit

        else:

            ul_energies.append(c_energy.real(row))

            ul_flux.append(c_upper.real(row))

            ul_ed_engs.append(c_ed.real(row))

            ul_eu_engs.append(c_eu.real(row))

# Set upper limit errors

    yerr = [0.6 * x for x in ul_flux]

    # Plot

    ax.errorbar(energies, flux, yerr=e_flux, xerr=[ed_engs, eu_engs],

                fmt='o'+color,alpha=0.4)

    ax.errorbar(ul_energies, ul_flux, xerr=[ul_ed_engs, ul_eu_engs],

                            yerr=yerr, uplims=True, fmt='o'+color,alpha=0.4)

fig6 = plt.figure()

ax = plt.subplot()

# define axes properties

ax.set_xscale('log')

ax.set_yscale('log')

ax.set_xlabel('Energy (TeV)')

ax.set_ylabel(r'E$^2$ $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)')

# plot results

plot_butterfly(butterfly,ax,color='b',label='%s'%(source_name))

plot_sed(sed,ax,color='g')

# legend

ax.legend(loc='best')

pylab.xlim([0.03,30])

#pylab.ylim([10**(-10),10**(-9)])

plt.savefig('butterfly_%s_0.03_%s.png'%(source_name,s))

#plt.show()

ulimit=ctools.ctulimit()

ulimit['inobs']   = evfile

ulimit['emin']      = emin

ulimit['emax']      = emax

ulimit['caldb']     = caldb

ulimit['irf']       = irf

ulimit['srcname'] = '%s'%(source_name)

ulimit['inmodel']   = 'input_loglaw.xml'

ulimit.execute()

'''