GammaLib

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= 18000

emin  = 0.03   # TeV

emax  = 30.0 # TeV

for i in range(len(DECJ2000)):

    ra=RAJ2000[i]

    dec=DECJ2000[i]

    source=Source_Name[i]

    if(25>dec>0):

        irf   = 'North_z20_50h'

    if(45>dec>25):

        irf   = 'North_z40_50h'

    if(65>dec>45):

        irf   = 'North_z60_50h'

    if(-25<dec<0):

        irf   = 'South_z20_50h'

    if(-45<dec<-25):

        irf   = 'South_z40_50h'

    if(-65<dec<-45):

        irf   = 'South_z60_50h'

    print dec,irf

    for j in range(10):

        s1=random.randint(1,1000)

        evfile = '/Volumes/mac_ext/egal_pop/PL_events_file/%s/events_%f.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

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

        obssim['outevents'] = evfile

        obssim.execute()

    '''

    like = ctools.ctlike()

    like['inobs']   = evfile

    like['caldb']   = caldb

    like['irf']     = irf

    like['inmodel'] = '%s.xml'%(source_name)

    like['outmodel'] = 'results_%s_0.03_%s.xml'%(source_name,s)

    like.execute()

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

    import re

    import sys

    import os

    import glob

    orig_stdout = sys.stdout

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

    sys.stdout = f 

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

    print s1

    sys.stdout = orig_stdout

    f.close()

    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

butterfly = 'butterfly_%s_0.03_%s.txt'%(source_name,s) # text file where the butterfly plot is stored

bttflmk = ctools.ctbutterfly()

bttflmk['inobs'] = 'events.fits'

bttflmk['inmodel'] = 'results_%s_0.03_%s.xml'%(source_name,s)

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

bttflmk['emin'] = emin

bttflmk['emax'] = emax

bttflmk['caldb']     = caldb

bttflmk['irf']       = irf

bttflmk['statistic'] = 'WSTAT'

bttflmk['outfile'] = butterfly

bttflmk.execute()

sed = "sed.fits"

sedmkr = cscripts.csspec()

sedmkr['inobs'] = 'events.fits'

sedmkr['inmodel'] = 'results_%s_0.03_%s.xml'%(source_name,s)

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

sedmkr['outfile'] = sed

sedmkr['statistic'] = 'WSTAT'

sedmkr['method'] = "AUTO"

sedmkr['emin'] = emin

sedmkr['emax'] = emax

sedmkr['caldb']     = caldb

sedmkr['irf']       = irf

sedmkr['enumbins'] = 10

sedmkr['ebinalg'] = 'LOG'

sedmkr.execute()

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()

'''