Example: gt_apps python Script Used in Binaries RSP Analysis

The like.py script shown below is one of four scripts (getRuns.py; select.py; like.py; and bin.py) that "currently live in /u/ey/richard/GLAST/mQ/pyMQ/ but will soon be in cvs"; see Richard Dubois' confluence page "Scripts used in Binaries RSP Analysis".

gt_apps

There are three key parts to gt_apps:

gt_apps import command;
gt_apps definitions; and the
.run() command (see last line of each of the definition groups).

#!/usr/bin/env python
"""
Generate test data with dataSubselector.

@author J. Chiang <jchiang@slac.stanford.edu>
"""
#
# $Header: /nfs/slac/g/glast/ground/cvs/workbook/pages/sciTools_pyLikelihood_tutorial/v02/example_scriptBinariesRSPanalysis.html,v 1.1 2009/09/10 17:53:49 chuckp Exp $
#

##from setPaths import *
from gt_apps import expMap, diffResps, expCube
from optparse import OptionParser
from UnbinnedAnalysis import *

def run(clean=False):

##
## Define command line arguments
##
parser = OptionParser(version="%prog $Revision 0.1",
description="run gtlike")
parser.add_option("-s","--scfile",action="store",type="string",dest="scfile",
default="FL-FT2-merged.fits",help="Filespec of spacecraft file")
parser.add_option("-e","--evfile",action="store",type="string",dest="evfile",
default="FL-LS5039-10.fits",help="Filespec of events file")
parser.add_option("--expfile",action="store",type="string",dest="expfile", default="mapfile.fits",help="output file name")

parser.add_option("--expCube",action="store",type="string",dest="expCube", default="",help="exposure cube output file")

parser.add_option("--irfs",action="store",type="string",dest="irfs", default="P6_V3_DIFFUSE",help="irfs to use")

parser.add_option("--sourceModel",action="store",type="string",dest="sourceModel", default="",help="source model file")

parser.add_option("--optimizer",action="store",type="string",dest="optimizer", default="NewMinuit",help="select minimisation package")

parser.add_option("--fluxmdl",action="store",type="string",dest="fluxmdl", default="none",help="output xml file of model")

parser.add_option("--likeOnly",action="store",type="string",dest="likeOnly", default="no",help="only run likelihood (yes/no)")

parser.add_option("--diffrsp",action="store",type="string",dest="diffrsp", default="no",help="force diffuse response (yes/no)")

parser.add_option("--plot",action="store",type="string",dest="plot", default="no",help="run the gui (yes/no)")

parser.add_option("--writeCounts",action="store",type="string",dest="writeCounts", default="counts_spectra.fits",help="file spec for counts_spectra fits file")
parser.add_option("--writeResults",action="store",type="string",dest="writeResults", default="results.txt", help="file spec for results.txt")
parser.add_option("--sourceName",action="store",type="string",dest="sourceName", default="",help="name of target source in xml file")

parser.add_option("--resultsLog",action="store",type="string",dest="resultsLog", default="resultsLog.txt",help="summary results output file")

parser.add_option("--emin",action="store",type="float",dest="emin", default=100.,help="min energy for energy integral")

parser.add_option("--emax",action="store",type="float",dest="emax", default=300000.,help="min energy for energy integral")

parser.add_option("--tol",action="store",type="float",dest="tol", default=0.1,help="tolerance value")

(options, args) = parser.parse_args()

expCube['evfile'] = options.evfile
expCube['scfile'] = options.scfile
expCube['outfile'] = options.expCube
expCube['dcostheta'] = 0.025
if (options.likeOnly == 'no'):
expCube.run()

expMap['evfile'] = options.evfile
expMap['outfile'] = options.expfile
expMap['scfile'] = options.scfile
expMap['irfs'] = options.irfs
expMap['expcube'] = expCube['outfile']
if (options.likeOnly == 'no'):
expMap.run()

diffResps['evfile'] = options.evfile
diffResps['scfile'] = options.scfile
diffResps['irfs'] = options.irfs
diffResps['srcmdl'] = options.sourceModel
diffResps['convert'] = 'yes'

if ((options.diffrsp == 'yes') or (options.likeOnly == 'no')):
diffResps.run()

likeObs = UnbinnedObs(eventFile=options.evfile,
scFile=options.scfile,
irfs=options.irfs,
expCube=expCube['outfile'],
expMap=expMap['outfile'])

like2 = UnbinnedAnalysis(likeObs,
srcModel=options.sourceModel)

like2.tol = options.tol
like2.tolType = 1

try:
like2.fit(optimizer='DRMNGB')
except RuntimeError:
try:
like2.logLike.restoreBestFit()
except RuntimeError:
like2.logLike.restoreBestFit()
pass

like2.fit(optimizer=options.optimizer,covar=True)
print like2.model
print options.sourceName + '\n'
print "TS = " + str(like2.Ts(options.sourceName))
print "Npred" + options.sourceName + " " + str(like2.logLike.NpredValue(options.sourceName))
print "NObs " + str(sum(like2._Nobs()))
print options.sourceName + ' flux ' +str(options.emin) + ' - ' + str(options.emax) + ' MeV \n'
prt1 = str(like2.flux(options.sourceName,emin=options.emin,emax=options.emax)) + " +/- " + str(like2.fluxError(options.sourceName, options.emin,options.emax))
print prt1
print options.sourceName + ' flux 0.1-0.3 GeV \n'
prt2 = str(like2.flux(options.sourceName,emin=100.,emax=300.)) + " +/- " + str(like2.fluxError(options.sourceName,emin=100.,emax=300.))
print prt2
print options.sourceName + ' flux 0.3-1 GeV \n'
prt3 = str(like2.flux(options.sourceName,emin=300.,emax=1000.)) + " +/- " + str(like2.fluxError(options.sourceName,emin=300.,emax=1000.))
print prt3
print options.sourceName + ' flux 1-3 GeV \n'
prt4 = str(like2.flux(options.sourceName,emin=1000.,emax=3000.)) + " +/- " + str(like2.fluxError(options.sourceName,emin=1000.,emax=3000.))
print prt4
print options.sourceName + ' flux 3-10 GeV \n'
prt5 = str(like2.flux(options.sourceName,emin=3000.,emax=10000.)) + " +/- " + str(like2.fluxError(options.sourceName,emin=3000.,emax=10000.))
print prt5
print options.sourceName + ' flux 10-300 GeV \n'
prt6 = str(like2.flux(options.sourceName,emin=10000.,emax=300000.)) + " +/- " + str(like2.fluxError(options.sourceName,emin=10000., emax=300000.))
print prt6

print options.sourceName + ' energyFlux 0.1-300 GeV \n'
prt11 = str(like2.energyFlux(options.sourceName,emin=100.,emax=300000)) + " +/- " + str(like2.energyFluxError(options.sourceName, emin=100., emax=300000.))
print prt11
print options.sourceName + ' energyFlux 0.1-0.3 GeV \n'
prt12 = str(like2.energyFlux(options.sourceName,emin=100.,emax=300.)) + " +/- " + str(like2.energyFluxError(options.sourceName, emin=100., emax=300.))
print prt12
print options.sourceName + ' energyFlux 0.3-1 GeV \n'
prt13 = str(like2.energyFlux(options.sourceName,emin=300.,emax=1000.)) + " +/- " + str(like2.energyFluxError(options.sourceName, emin=300., emax=1000.))
print prt13
print options.sourceName + ' energyFlux 1-3 GeV \n'
prt14 = str(like2.energyFlux(options.sourceName,emin=1000.,emax=3000.)) + " +/- " + str(like2.energyFluxError(options.sourceName, emin=1000., emax=3000.))
print prt14
print options.sourceName + ' energyFlux 3-10 GeV \n'
prt15 = str(like2.energyFlux(options.sourceName,emin=3000.,emax=10000.)) + " +/- " + str(like2.energyFluxError(options.sourceName, emin=3000., emax=10000.))
print prt15
print options.sourceName + ' energyFlux 10-300 GeV \n'
prt16 = str(like2.energyFlux(options.sourceName,emin=10000.,emax=300000.)) + " +/- " + str(like2.energyFluxError(options.sourceName, emin=10000., emax=300000.))
print prt16

resultsFile = open(options.writeResults,'w')
print >> resultsFile, like2.model
print >> resultsFile, options.sourceName, ' TS \n'
print >> resultsFile, like2.Ts(options.sourceName)
print >> resultsFile, options.sourceName + ' flux 0.1-300 GeV \n'
print >> resultsFile, prt1
print >> resultsFile, options.sourceName + ' flux 0.1-0.3 GeV \n'
print >> resultsFile, prt2
print >> resultsFile, options.sourceName + ' flux 0.3-1 GeV \n'
print >> resultsFile, prt3
print >> resultsFile, options.sourceName + ' flux 1-3 GeV \n'
print >> resultsFile, prt4
print >> resultsFile, options.sourceName + ' flux 3-10 GeV \n'
print >> resultsFile, prt5
print >> resultsFile, options.sourceName + ' flux 10-300 GeV \n'
print >> resultsFile, prt6

print >> resultsFile, options.sourceName + ' energyFlux 0.1-300 GeV \n'
print >> resultsFile, prt11
print >> resultsFile, options.sourceName + ' energyFlux 0.1-0.3 GeV \n'
print >> resultsFile, prt12
print >> resultsFile, options.sourceName + ' energyFlux 0.3-1 GeV \n'
print >> resultsFile, prt13
print >> resultsFile, options.sourceName + ' energyFlux 1-3 GeV \n'
print >> resultsFile, prt14
print >> resultsFile, options.sourceName + ' energyFlux 3-10 GeV \n'
print >> resultsFile, prt15
print >> resultsFile, options.sourceName + ' energyFlux 10-300 GeV \n'
print >> resultsFile, prt16

resultsFile.close()

like2.writeCountsSpectra(options.writeCounts)
if (options.fluxmdl != "none"):
like2.writeXml(options.fluxmdl)

if __name__ == "__main__":
run()

Owned by: Richard Dubois

Last updated by: Chuck Patterson 09/09/2009