ctools

import gammalib

import ctools

import cscripts

ra         = 350.8524

dec        = 58.8159

max_offset = 3.

srcrad     = 0.2 #[deg]

stacking   = True

emin       = 0.1 #TeV

emax       = 50.

enumbins   = 30

obs = 'obs.xml'

obsselect = cscripts.csobsselect()

obsselect['inobs']     = '$CTADATA/obs/obs_gps_baseline.xml'

obsselect['outobs']    = obs

obsselect['pntselect'] = 'CIRCLE'

obsselect['coordsys']  = 'CEL'

obsselect['ra']        = ra

obsselect['dec']       = dec

obsselect['rad']       = max_offset

obsselect['tmin']      = 'UNDEF'

obsselect.execute()

obs_selected = "obs_selected.xml"

select = ctools.ctselect()

select['inobs']  = obs

select['ra']     = 'UNDEF'

select['tmin']   = 'UNDEF'

select['emin']   = emin

select['emax']   = emax

select['expr']   = 'MC_ID==78'

select['outobs'] = obs_selected

select.execute()

model_sky = 'model_sky.xml'

modelselect = cscripts.csmodelselect()

modelselect['inobs']    = obs_selected

modelselect['inmodel']  = '$CTADATA/models/model_galactic_bright.xml'

modelselect['outmodel'] = model_sky

modelselect.execute()

modelfile = model_sky

model = gammalib.GModels(modelfile)

for source in model:

    if source.name() == 'Cassiopeia A':

        for par in source:

            if par.name() == 'RA' or par.name() == 'DEC':

                par.fix()

model['Cassiopeia A'].tscalc(True)

model.save(modelfile)

onoff_obs = 'onoff_obs.xml' #output file

phatool = cscripts.csphagen()

phatool['inobs']     = obs_selected

phatool['outobs']    = onoff_obs

phatool['ebinalg']   = 'LOG'

phatool['emin']      = emin

phatool['emax']      = emax

phatool['enumbins']  = enumbins

phatool['coordsys']  = 'CEL'

phatool['ra']        = ra

phatool['dec']       = dec

phatool['rad']       = srcrad #deg

phatool['bkgmethod'] = 'REFLECTED'

phatool['maxoffset'] = max_offset

phatool['stack']     = stacking # spectra from all observations are stacked into a single one

phatool['debug']     = True

phatool.execute()

outmodelfile_onoff = 'bestfit_model_onoff.xml'

like = ctools.ctlike()

like['inobs']     = onoff_obs

like['inmodel']   = modelfile

like['statistic'] = 'WSTAT' #background from likelihood profiling

like['outmodel']  = outmodelfile_onoff

like.execute()

sed_onoff = "spectrum_issue.fits"

sedmkr = cscripts.csspec()

sedmkr['inobs'] = onoff_obs

sedmkr['inmodel'] = outmodelfile_onoff

sedmkr['srcname'] = 'Cassiopeia A'

sedmkr['outfile'] = sed_onoff

sedmkr['statistic'] = 'WSTAT'

sedmkr['method'] = "AUTO"

sedmkr['emin'] = emin

sedmkr['emax'] = emax

sedmkr['enumbins'] = enumbins

sedmkr['ebinalg'] = 'LOG'

sedmkr['debug']   = True

sedmkr.execute()

outmodelfile_onoff_bkg = 'bestfit_model_onoff_bkg.xml'

bkg = gammalib.GCTAModelIrfBackground( gammalib.GModelSpectralPlaw( 0, 0, gammalib.GEnergy( 0, 'TeV' ) ) )

bkg.instruments( 'CTAOnOff' )

model_bkg = gammalib.GModels(outmodelfile_onoff)

model_bkg.append(bkg)

model_bkg.save(outmodelfile_onoff_bkg)

sed_onoff = "spectrum_works.fits"

sedmkr = cscripts.csspec()

sedmkr['inobs'] = onoff_obs

sedmkr['inmodel'] = outmodelfile_onoff_bkg

sedmkr['srcname'] = 'Cassiopeia A'

sedmkr['outfile'] = sed_onoff

sedmkr['statistic'] = 'WSTAT'

sedmkr['method'] = "AUTO"

sedmkr['emin'] = emin

sedmkr['emax'] = emax

sedmkr['enumbins'] = enumbins

sedmkr['ebinalg'] = 'LOG'

sedmkr['debug']   = True

sedmkr.execute()