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  = 0.03  # TeV

emax  = 199 # TeV

irf_used=[]

decs=[]

cta_north_lat=28.71

cta_south_lat=-24.63

alt_cul_north=[]

alt_cul_south=[]

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

zenith_n=[]

zenith_s=[]

for i in range(len(alt_cul_north)):

    zenith_n.append(90-alt_cul_north[i])

for i in range(len(alt_cul_south)):

    zenith_s.append(90-alt_cul_south[i])

'''

final_ra=[]

final_dec=[]

final_ra_rad=[]

final_dec_rad=[]

ra2=[]

for i in range(len(RAJ2000)):

    ra2.append(RAJ2000[i])

for i in range(len(ra2)):

    if(ra2[i]>180):

        ra2[i]=ra2[i]-360

for i in range(len(ra2)):

    final_ra_rad.append((ra2[i])*0.0174533)

    final_dec_rad.append((DECJ2000[i])*0.0174533)

for i in range(len(final_ra)):

    final_ra_rad.append((final_ra[i])*0.0174533)

    final_dec_rad.append((final_dec[i])*0.0174533)

fig = plt.figure(figsize=(20,16))

ax = fig.add_subplot(111, projection="mollweide")

cm = plt.cm.get_cmap('jet')

plt.grid(True)

for i in range(len(final_ra_rad)):

    if(0<zenith_n[i]<90):

        SC =ax.scatter(final_ra_rad[i], final_dec_rad[i], c=zenith_n[i],cmap=cm,vmin=0, vmax=100)

cbar = plt.colorbar(SC, shrink=1., pad=0.05)

cbar.ax.tick_params(labelsize=20)

cbar.set_label('Zenith Angle', fontsize=20)

plt.title('CTA N')

plt.savefig('Plots/CTA_N_zenith.png',dpi=100)

#plt.show()

fig = plt.figure(figsize=(20,16))

ax = fig.add_subplot(111, projection="mollweide")

cm = plt.cm.get_cmap('jet')

plt.grid(True)

for i in range(len(final_ra_rad)):

    if(0<zenith_s[i]<90):

        SC =ax.scatter(final_ra_rad[i], final_dec_rad[i], c=zenith_s[i],cmap=cm,vmin=0, vmax=100)

cbar = plt.colorbar(SC, shrink=1., pad=0.05)

cbar.ax.tick_params(labelsize=20)

cbar.set_label('Zenith Angle', fontsize=20)

plt.title('CTA S')

plt.savefig('Plots/CTA_S_zenith.png',dpi=100)

'''

north_irf=[]

south_irf=[]

for i in range(len(zenith_n)):

    if(0<zenith_n[i]<25):

        north_irf.append('North_z20_50h')

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

        north_irf.append('North_z40_50h')

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

        north_irf.append('North_z60_50h')

    else:

        north_irf.append('North_z20_50h')

for i in range(len(zenith_s)):

    if(0<zenith_s[i]<25):

        south_irf.append('South_z20_50h')

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

        south_irf.append('South_z40_50h')

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

        south_irf.append('South_z60_50h')

    else:

        south_irf.append('South_z20_50h')

for i in range(826,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,irf_n,irf_s

    for j in range(1):

    #s1=random.randint(1,1000)

        evfile = 'PL_events_file_20hr/%s_events_n30GeV_%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_30

        obssim['nthreads']=0

        obssim['chatter']=2

    #obssim['clobber']='yes'

        obssim['mode']='ql'

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

        obssim['inmodel']   = '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 = 'PL_events_file_20hr/%s_events_s30GeV_%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_30

        obssim['inmodel']   = '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 = 'PL_events_file_20hr/%s_events_n30GeV_%i.fits'%(source,j)

        like['inobs']   = evfile

        like['caldb']   = caldb

        like['irf']     = irf_n_30

        like['inmodel'] = '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'] = 'Fit_results_20hr/results_n30_%s_%s.xml'%(source,j)

        like.execute()

        import re

        import sys

        import os

        import glob

        orig_stdout = sys.stdout

        f = open('Fit_results_20hr/results_n30_%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 = 'PL_events_file_20hr/%s_events_s30GeV_%i.fits'%(source,j)

        like['inobs']   = evfile

        like['caldb']   = caldb

        like['irf']     = irf_s_30

        like['inmodel'] = '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'] = 'Fit_results_20hr/results_s30_%s_%s.xml'%(source,j)

        like.execute()

        orig_stdout = sys.stdout

        f = open('Fit_results_20hr/results_s30_%s_%s.txt'%(source,j),'w')

        sys.stdout = f

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

        #print s1

        sys.stdout = orig_stdout

        f.close()