simArp220.py
1 |
#! /usr/bin/env python
|
---|---|
2 |
# ==========================================================================
|
3 |
# This script simulates CTA observations of Arp220
|
4 |
#
|
5 |
# Author: Rolf Buehler (rolf.buehler@desy.de)
|
6 |
# Date: 20.3.2014
|
7 |
# ==========================================================================
|
8 |
import gammalib as gl |
9 |
import obsutils as obsu |
10 |
|
11 |
def simarp220(obstime=100*3600.,write=False): |
12 |
""""Simulates and fits Arp220 & background for one run
|
13 |
@param [float] obstime Observation duration in seconds
|
14 |
@param [bool] write Write out simulated observations adn fit results
|
15 |
"""
|
16 |
|
17 |
#-----------------------------------------------------
|
18 |
#Create one observation (run)
|
19 |
#-----------------------------------------------------
|
20 |
ra_degree=233.7382 #Arp220 position from NED |
21 |
dec_degree=23.503
|
22 |
offset_degree=-0.5 #Observing offset in dec from Arp220 in deg |
23 |
pntdir = gl.GSkyDir() |
24 |
pntdir.radec_deg(ra_degree, dec_degree+offset_degree) |
25 |
|
26 |
#irfs = "data/cta/aar/bcf/DESY20140105_50h/irf_file.fits"
|
27 |
#filename="$CALDB/data/cta/aar/bcf/DESY20140105_50h/irf_file.fits"
|
28 |
#caldbs = "aar"
|
29 |
|
30 |
irfs = "data/cta/i/bcf/IFAE20120510_50h/irf_file.fits"
|
31 |
filename="$CALDB/data/cta/i/bcf/IFAE20120510_50h/irf_file.fits"
|
32 |
caldbs = "i"
|
33 |
|
34 |
ob = obsu.set(pntdir,tstart=0.0,duration=obstime,deadc=0.95, \ |
35 |
emin=0.1, emax=10.0, rad=3.0, irf=irfs,caldb=caldbs) |
36 |
ob.id("run_1")
|
37 |
|
38 |
#----------------------------------------------------
|
39 |
# Create background model
|
40 |
#----------------------------------------------------
|
41 |
fixallpars=True
|
42 |
spec = gl.GModelSpectralPlaw(1,0,gl.GEnergy(1.0,"TeV")) |
43 |
bck = gl.GCTAModelBackground(ob,filename,spec) |
44 |
bck.name("Background")
|
45 |
bck.instruments("CTA")
|
46 |
if fixallpars :
|
47 |
for par in range(0,bck.size()): |
48 |
bck[par].fix() |
49 |
|
50 |
#----------------------------------------------------
|
51 |
# Create Arp220 model
|
52 |
#----------------------------------------------------
|
53 |
arp220dir = gl.GSkyDir() |
54 |
arp220dir.radec_deg(ra_degree, dec_degree) |
55 |
arp220_spatial = gl.GModelSpatialPointSource(arp220dir) |
56 |
|
57 |
arp220_spectrum = gl.GModelSpectralPlaw(1.0e-19, -2.2,gl.GEnergy(1,"TeV") ) |
58 |
arp220 = gl.GModelSky(arp220_spatial, arp220_spectrum) |
59 |
arp220.name("Arp220")
|
60 |
arp220.ids("")
|
61 |
|
62 |
#----------------------------------------------------
|
63 |
# Create GObservations and add models
|
64 |
#----------------------------------------------------
|
65 |
obs = gl.GObservations() |
66 |
obs.append(ob) |
67 |
obs.models().append(bck) |
68 |
obs.models().append(arp220) |
69 |
|
70 |
#----------------------------------------------------
|
71 |
# Simulate, fit and write out
|
72 |
#----------------------------------------------------
|
73 |
obssim = obsu.sim(obs, log=write, debug=False, edisp=False, seed=4, nbins=40, binsz=0.02) |
74 |
if write:
|
75 |
obssim.save("arp220_obs.xml")
|
76 |
obssim.models().save("arp220_sim.xml")
|
77 |
print "Creating counts map.." |
78 |
obsu.cntmap(obssim, proj="TAN", coord="EQU", xval=ra_degree, yval=dec_degree+offset_degree, \ |
79 |
binsz=0.02, nxpix=300, nypix=300, outname="arp220_cntmap.fits") |
80 |
print "Creating model map.." |
81 |
obsu.modmap(obssim, eref=0.1, proj="TAN", coord="EQU", xval=ra_degree, yval=dec_degree+offset_degree, \ |
82 |
binsz=0.02, nxpix=300, nypix=300, outname="arp220_modmap.fits") |
83 |
for ob in obssim: |
84 |
events = ob.events() |
85 |
events.save("arp220_events_run_1.fits")
|
86 |
print "Fitting.." |
87 |
obsu.fit(obssim,log=write, debug=True, edisp=False) |
88 |
|
89 |
if __name__ == '__main__': |
90 |
simarp220() |