Johan Bregeon / ctools

Commit 7a8e2181124fce8bbd8b129db9189cb32ef324c1

Authored by Jürgen Knödlseder
1 parent 8830ee7c

Add testing of diffuse map cubes to Python unit test.

Inline Side-by-side
Showing 2 changed files with 221 additions and 161 deletions
test/data/model_background_cube_const.xml 0 → 100644
  View file @7a8e218
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="CTABackgroundModel" type="CTABackground" instrument="CTA">
  4 + <spectrum type="ConstantValue">
  5 + <parameter name="Value" scale="1.0" min="0.01" max="1000.0" value="1.0" free="1"/>
  6 + </spectrum>
  7 + <spatialModel type="MapCubeFunction" file="data/test_cube.fits">
  8 + <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
  9 + </spatialModel>
  10 + </source>
  11 +</source_library>
test/test_python.py
  View file @7a8e218
@@ -30,178 +30,227 @@ import sys @@ -30,178 +30,227 @@ import sys
30 # Analysis pipeline - version 1 # 30 # Analysis pipeline - version 1 #
31 # ============================= # 31 # ============================= #
32 def pipeline_v1(): 32 def pipeline_v1():
33 - """  
34 - Unbinned analysis pipeline - save intermediate results. 33 + """
  34 + Unbinned analysis pipeline - save intermediate results.
35 35
36 - This function implements an analysis pipeline that successively calls  
37 - ctobssim, ctselect and ctlike by saving the intermediate results as FITS  
38 - files on disk. This replicates the way how the analysis would be done  
39 - in a ftools-like approach.  
40 - """  
41 - # Set script parameters  
42 - model_name = "data/crab.xml"  
43 - events_name = "events.fits"  
44 - selected_events_name = "selected_events.fits"  
45 - result_name = "results.xml"  
46 - caldb = "irf"  
47 - irf = "cta_dummy_irf"  
48 - ra = 83.63  
49 - dec = 22.01  
50 - rad_sim = 10.0  
51 - tstart = 0.0  
52 - tstop = 1800.0  
53 - emin = 0.1  
54 - emax = 100.0  
55 - rad_select = 3.0  
56 -  
57 - # Initialise timing  
58 - wall_seconds = 0.0  
59 - cpu_seconds = 0.0  
60 -  
61 - # Simulate events  
62 - sim = ctobssim()  
63 - sim.logFileOpen()  
64 - sim["infile"].filename(model_name)  
65 - sim["outfile"].filename(events_name)  
66 - sim["caldb"].string(caldb)  
67 - sim["irf"].string(irf)  
68 - sim["ra"].real(ra)  
69 - sim["dec"].real(dec)  
70 - sim["rad"].real(rad_sim)  
71 - sim["tmin"].real(tstart)  
72 - sim["tmax"].real(tstop)  
73 - sim["emin"].real(emin)  
74 - sim["emax"].real(emax)  
75 - sim.execute()  
76 - sys.stdout.write(".")  
77 -  
78 - # Select events  
79 - select = ctselect()  
80 - select.logFileOpen()  
81 - select["infile"].filename(events_name)  
82 - select["outfile"].filename(selected_events_name)  
83 - select["ra"].real(ra)  
84 - select["dec"].real(dec)  
85 - select["rad"].real(rad_select)  
86 - select["tmin"].real(tstart)  
87 - select["tmax"].real(tstop)  
88 - select["emin"].real(emin)  
89 - select["emax"].real(emax)  
90 - select.execute()  
91 - sys.stdout.write(".")  
92 -  
93 - # Perform maximum likelihood fitting  
94 - like = ctlike()  
95 - like.logFileOpen()  
96 - like["infile"].filename(selected_events_name)  
97 - like["srcmdl"].filename(model_name)  
98 - like["outmdl"].filename(result_name)  
99 - like["caldb"].string(caldb)  
100 - like["irf"].string(irf)  
101 - like.execute()  
102 - sys.stdout.write(".")  
103 -  
104 - # Return  
105 - return 36 + This function implements an analysis pipeline that successively calls
  37 + ctobssim, ctselect and ctlike by saving the intermediate results as FITS
  38 + files on disk. This replicates the way how the analysis would be done
  39 + in a ftools-like approach.
  40 + """
  41 + # Set script parameters
  42 + model_name = "data/crab.xml"
  43 + events_name = "events.fits"
  44 + selected_events_name = "selected_events.fits"
  45 + result_name = "results.xml"
  46 + caldb = "irf"
  47 + irf = "cta_dummy_irf"
  48 + ra = 83.63
  49 + dec = 22.01
  50 + rad_sim = 10.0
  51 + tstart = 0.0
  52 + tstop = 1800.0
  53 + emin = 0.1
  54 + emax = 100.0
  55 + rad_select = 3.0
  56 +
  57 + # Initialise timing
  58 + wall_seconds = 0.0
  59 + cpu_seconds = 0.0
  60 +
  61 + # Simulate events
  62 + sim = ctobssim()
  63 + sim.logFileOpen()
  64 + sim["infile"].filename(model_name)
  65 + sim["outfile"].filename(events_name)
  66 + sim["caldb"].string(caldb)
  67 + sim["irf"].string(irf)
  68 + sim["ra"].real(ra)
  69 + sim["dec"].real(dec)
  70 + sim["rad"].real(rad_sim)
  71 + sim["tmin"].real(tstart)
  72 + sim["tmax"].real(tstop)
  73 + sim["emin"].real(emin)
  74 + sim["emax"].real(emax)
  75 + sim.execute()
  76 + sys.stdout.write(".")
  77 +
  78 + # Select events
  79 + select = ctselect()
  80 + select.logFileOpen()
  81 + select["infile"].filename(events_name)
  82 + select["outfile"].filename(selected_events_name)
  83 + select["ra"].real(ra)
  84 + select["dec"].real(dec)
  85 + select["rad"].real(rad_select)
  86 + select["tmin"].real(tstart)
  87 + select["tmax"].real(tstop)
  88 + select["emin"].real(emin)
  89 + select["emax"].real(emax)
  90 + select.execute()
  91 + sys.stdout.write(".")
  92 +
  93 + # Perform maximum likelihood fitting
  94 + like = ctlike()
  95 + like.logFileOpen()
  96 + like["infile"].filename(selected_events_name)
  97 + like["srcmdl"].filename(model_name)
  98 + like["outmdl"].filename(result_name)
  99 + like["caldb"].string(caldb)
  100 + like["irf"].string(irf)
  101 + like.execute()
  102 + sys.stdout.write(".")
  103 +
  104 + # Return
  105 + return
106 106
107 107
108 # ============================= # 108 # ============================= #
109 # Analysis pipeline - version 2 # 109 # Analysis pipeline - version 2 #
110 # ============================= # 110 # ============================= #
111 -def pipeline_v2():  
112 - """  
113 - Unbinned analysis pipeline - keep intermediate results in memory. 111 +def pipeline_v2(model_name="data/crab.xml", \
  112 + ra_sim=83.63, dec_sim=22.01, rad_sim=10.0, \
  113 + ra_select=83.63, dec_select=22.01, rad_select=3.0, \
  114 + tstart=0.0, tstop=1800.0, \
  115 + emin=0.1, emax=100.0, \
  116 + ref={}):
  117 + """
  118 + Unbinned analysis pipeline - keep intermediate results in memory.
114 119
115 - This function implements an analysis pipeline that successively calls  
116 - ctobssim, ctselect and ctlike without saving the intermediate results as  
117 - FITS files on disk. All data is only hold in memory.  
118 - """  
119 - # Set script parameters  
120 - model_name = "data/crab.xml"  
121 - caldb = "irf"  
122 - irf = "cta_dummy_irf"  
123 - ra = 83.63  
124 - dec = 22.01  
125 - rad_sim = 10.0  
126 - tstart = 0.0  
127 - tstop = 1800.0  
128 - emin = 0.1  
129 - emax = 100.0  
130 - rad_select = 3.0  
131 -  
132 - # Initialise timing  
133 - wall_seconds = 0.0  
134 - cpu_seconds = 0.0  
135 -  
136 - # Simulate events  
137 - sim = ctobssim()  
138 - sim["infile"].filename(model_name)  
139 - sim["caldb"].string(caldb)  
140 - sim["irf"].string(irf)  
141 - sim["ra"].real(ra)  
142 - sim["dec"].real(dec)  
143 - sim["rad"].real(rad_sim)  
144 - sim["tmin"].real(tstart)  
145 - sim["tmax"].real(tstop)  
146 - sim["emin"].real(emin)  
147 - sim["emax"].real(emax)  
148 - sim.run()  
149 - sys.stdout.write(".")  
150 -  
151 - # Select events  
152 - select = ctselect(sim.obs())  
153 - select["ra"].real(ra)  
154 - select["dec"].real(dec)  
155 - select["rad"].real(rad_select)  
156 - select["tmin"].real(tstart)  
157 - select["tmax"].real(tstop)  
158 - select["emin"].real(emin)  
159 - select["emax"].real(emax)  
160 - select.run()  
161 - sys.stdout.write(".")  
162 -  
163 - # Perform maximum likelihood fitting  
164 - like = ctlike(select.obs())  
165 - like.run()  
166 - sys.stdout.write(".")  
167 -  
168 - # Return  
169 - return 120 + This function implements an analysis pipeline that successively calls
  121 + ctobssim, ctselect and ctlike without saving the intermediate results as
  122 + FITS files on disk. All data is only hold in memory.
  123 + """
  124 + # Set script parameters
  125 + #model_name = "data/crab.xml"
  126 + caldb = "irf"
  127 + irf = "cta_dummy_irf"
  128 + #ra = 83.63
  129 + #dec = 22.01
  130 + #rad_sim = 10.0
  131 + #tstart = 0.0
  132 + #tstop = 1800.0
  133 + #emin = 0.1
  134 + #emax = 100.0
  135 + #rad_select = 3.0
  136 +
  137 + # Simulate events
  138 + sim = ctobssim()
  139 + sim["infile"].filename(model_name)
  140 + sim["caldb"].string(caldb)
  141 + sim["irf"].string(irf)
  142 + sim["ra"].real(ra_sim)
  143 + sim["dec"].real(dec_sim)
  144 + sim["rad"].real(rad_sim)
  145 + sim["tmin"].real(tstart)
  146 + sim["tmax"].real(tstop)
  147 + sim["emin"].real(emin)
  148 + sim["emax"].real(emax)
  149 + sim.run()
  150 + sys.stdout.write(".")
  151 +
  152 + # Select events
  153 + select = ctselect(sim.obs())
  154 + select["ra"].real(ra_select)
  155 + select["dec"].real(dec_select)
  156 + select["rad"].real(rad_select)
  157 + select["tmin"].real(tstart)
  158 + select["tmax"].real(tstop)
  159 + select["emin"].real(emin)
  160 + select["emax"].real(emax)
  161 + select.run()
  162 + sys.stdout.write(".")
  163 +
  164 + # Perform maximum likelihood fitting
  165 + like = ctlike(select.obs())
  166 + like.run()
  167 + sys.stdout.write(".")
  168 +
  169 + # Analyse results
  170 + errors = 0
  171 + models = like.obs().models()
  172 + #print(models)
  173 + for model in models:
  174 + for key in ref:
  175 + if key == model.name():
  176 + for par in model:
  177 + for value in ref[key]:
  178 + if value == par.name():
  179 + if abs(float(ref[key][value])-par.value()) < 0.0001:
  180 + sys.stdout.write(".")
  181 + else:
  182 + sys.stdout.write("E")
  183 + errors += 1
  184 +
  185 + # Return number of errors
  186 + return errors
170 187
171 188
172 #==========================# 189 #==========================#
173 # Main routine entry point # 190 # Main routine entry point #
174 #==========================# 191 #==========================#
175 if __name__ == '__main__': 192 if __name__ == '__main__':
176 - """  
177 - Test ctools Python modules.  
178 - """  
179 - # Dump header  
180 - sys.stdout.write("\n")  
181 - sys.stdout.write("***********************************\n")  
182 - sys.stdout.write("* ctools Python interface testing *\n")  
183 - sys.stdout.write("***********************************\n")  
184 -  
185 - # Set PFILES environment variable  
186 - try:  
187 - os.mkdir("pfiles")  
188 - except:  
189 - pass  
190 - os.system("cp -r ../src/*/*.par pfiles/")  
191 - os.environ['PFILES'] = "pfiles"  
192 -  
193 - # Remove any existing result files  
194 - list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]  
195 - for files in list:  
196 - for file in files:  
197 - os.remove(file)  
198 -  
199 - # Save intermediate results on disk  
200 - sys.stdout.write("Test executable analysis: ")  
201 - pipeline_v1()  
202 - sys.stdout.write(" ok\n") 193 + """
  194 + Test ctools Python modules.
  195 + """
  196 + # Dump header
  197 + print("")
  198 + print("***********************************")
  199 + print("* ctools Python interface testing *")
  200 + print("***********************************")
  201 +
  202 + # Set PFILES environment variable
  203 + try:
  204 + os.mkdir("pfiles")
  205 + except:
  206 + pass
  207 + os.system("cp -r ../src/*/*.par pfiles/")
  208 + os.environ['PFILES'] = "pfiles"
  209 +
  210 + # Remove any existing result files
  211 + list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
  212 + for files in list:
  213 + for file in files:
  214 + os.remove(file)
  215 +
  216 + # Save intermediate results on disk
  217 + sys.stdout.write("Test executable analysis: ")
  218 + pipeline_v1()
  219 + print(" ok")
  220 +
  221 + # Keep results in memory
  222 + sys.stdout.write("Test in-memory analysis: ")
  223 + err = pipeline_v2()
  224 + if err == 0:
  225 + print(" ok")
  226 + else:
  227 + print(" NOK")
  228 +
  229 + # Test diffuse map cube with constant spectrum
  230 + sys.stdout.write("Test diffuse map cube with constant analysis: ")
  231 + ref = {'CTABackgroundModel': {'Value': 1.03547}}
  232 + err = pipeline_v2(model_name="data/model_background_cube_const.xml", \
  233 + ra_sim=84.17263, dec_sim=22.01444, rad_sim=2.0, \
  234 + ra_select=83.80, dec_select=22.20, rad_select=1.0, \
  235 + tstart=0.0, tstop=1.0, \
  236 + emin=0.3, emax=8.0, \
  237 + ref=ref)
  238 + if err == 0:
  239 + print(" ok")
  240 + else:
  241 + print(" NOK")
  242 +
  243 + # Test diffuse map cube with power law spectrum
  244 + sys.stdout.write("Test diffuse map cube with power law analysis: ")
  245 + ref = {'CTABackgroundModel': {'Prefactor': 1.02396, 'Index': -0.00710087}}
  246 + err = pipeline_v2(model_name="data/model_background_cube.xml", \
  247 + ra_sim=84.17263, dec_sim=22.01444, rad_sim=2.0, \
  248 + ra_select=84.00, dec_select=22.20, rad_select=1.0, \
  249 + tstart=0.0, tstop=1.0, \
  250 + emin=0.3, emax=8.0, \
  251 + ref=ref)
  252 + if err == 0:
  253 + print(" ok")
  254 + else:
  255 + print(" NOK")
203 256
204 - # Keep results in memory  
205 - sys.stdout.write("Test in-memory analysis: ")  
206 - pipeline_v2()  
207 - sys.stdout.write(" ok\n")