Luigi Tibaldo / ctools

Commit 5c53c14a266ac43c11429af62c1fe6dc9e88ae2d

Authored by Christoph Deil
1 parent d46008b1

Fix Python source code formatting (PEP8) 

This command was run:
  autopep8 --select E101,W191 --in-place --recursive .

It fixes only these two issues:
  - W191 indentation contains tabs
  - E101 indentation contains mixed spaces and tabs
Inline Side-by-side
Showing 53 changed files with 1492 additions and 1508 deletions
cscripts/cscaldb.py
  View file @5c53c14
... ... @@ -42,7 +42,7 @@ class cscaldb(ctools.cscript):
42 42 # Set name
43 43 self.name = "cscaldb"
44 44 self.version = "1.0.0"
45   -
  45 +
46 46 # Make sure that parfile exists
47 47 file = self.parfile()
48 48  
... ... @@ -60,7 +60,7 @@ class cscaldb(ctools.cscript):
60 60  
61 61 # Return
62 62 return
63   -
  63 +
64 64 def __del__(self):
65 65 """
66 66 Destructor.
... ... @@ -75,35 +75,35 @@ class cscaldb(ctools.cscript):
75 75 """
76 76 # Set parfile name
77 77 parfile = self.name+".par"
78   -
  78 +
79 79 try:
80 80 pars = gammalib.GApplicationPars(parfile)
81 81 except:
82 82 # Signal if parfile was not found
83 83 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
84   -
  84 +
85 85 # Create default parfile
86 86 pars = gammalib.GApplicationPars()
87 87 pars.append_standard()
88 88 pars.save(parfile)
89   -
  89 +
90 90 # Return
91 91 return
92   -
  92 +
93 93 def get_parameters(self):
94 94 """
95 95 Get parameters from parfile.
96 96 """
97 97 # Return
98 98 return
99   -
  99 +
100 100 def execute(self):
101 101 """
102 102 Execute the script.
103 103 """
104 104 # Run the script
105 105 self.run()
106   -
  106 +
107 107 # Return
108 108 return
109 109  
... ... @@ -117,7 +117,7 @@ class cscaldb(ctools.cscript):
117 117  
118 118 # Get parameters
119 119 self.get_parameters()
120   -
  120 +
121 121 # Write input parameters into logger
122 122 if self.logTerse():
123 123 self.log_parameters()
... ... @@ -156,7 +156,7 @@ class cscaldb(ctools.cscript):
156 156  
157 157 # Loop over instruments
158 158 for instrument in instruments:
159   -
  159 +
160 160 # Write mission into logger
161 161 if self.logTerse():
162 162 self.log.header3("Response functions in database \""+instrument+"\"")
... ... @@ -166,7 +166,7 @@ class cscaldb(ctools.cscript):
166 166 fits = gammalib.GFits(cifname)
167 167 cif = fits["CIF"]
168 168 cals = cif["CAL_CBD"]
169   -
  169 +
170 170 # Extract response names
171 171 names = []
172 172 nrows = cals.length()
... ... @@ -181,13 +181,13 @@ class cscaldb(ctools.cscript):
181 181 if names.count(name) == 0:
182 182 names.append(name)
183 183 names.sort()
184   -
  184 +
185 185 # Print response name
186 186 if self.logTerse():
187 187 for name in names:
188 188 self.log(name+"\n")
189 189 self.log("\n")
190   -
  190 +
191 191 # Return
192 192 return
193 193  
... ... @@ -201,10 +201,9 @@ if __name__ == '__main__':
201 201 """
202 202 # Create instance of application
203 203 app = cscaldb(sys.argv)
204   -
  204 +
205 205 # Open logfile
206 206 app.logFileOpen()
207   -
  207 +
208 208 # Execute application
209 209 app.execute()
210   -
211 210 \ No newline at end of file
... ...
cscripts/cshessobs.py
  View file @5c53c14
... ... @@ -52,12 +52,12 @@ class cshessobs(ctools.cscript):
52 52 except KeyError:
53 53 raise RuntimeError("Environment variable \"HESSFITS\" not set. " \
54 54 "Please set the environment variable \"HESSFITS\" to use this script")
55   -
  55 +
56 56 # Initialise some members
57 57 self.inmodels = None
58 58 self.outobs = "obs.xml"
59 59 self.nodes = 1
60   -
  60 +
61 61 # Make sure that parfile exists
62 62 file = self.parfile()
63 63  
... ... @@ -75,7 +75,7 @@ class cshessobs(ctools.cscript):
75 75  
76 76 # Return
77 77 return
78   -
  78 +
79 79 def __del__(self):
80 80 """
81 81 Destructor.
... ... @@ -88,64 +88,64 @@ class cshessobs(ctools.cscript):
88 88 Check if parfile exists. If parfile does not exist then create a
89 89 default parfile. This kluge avoids shipping the cscript with a parfile.
90 90 """
91   -
  91 +
92 92 # Set parfile name
93 93 parfile = self.name+".par"
94   -
  94 +
95 95 try:
96 96 pars = gammalib.GApplicationPars(parfile)
97 97 except:
98 98 # Signal that parfile was not found
99 99 print("Parfile \""+parfile+"\" not found. Create default parfile.")
100   -
  100 +
101 101 # Check for folders and possible/available configurations
102 102 chains = os.listdir(self.datapath)
103   -
  103 +
104 104 # Initialise folder dictionary
105 105 folders = {}
106   -
  106 +
107 107 # Loop over possible analysis chains
108 108 for chain in chains:
109   -
  109 +
110 110 # Omit hidden files/folders
111 111 if chain[0]=="." or not os.path.isdir(self.datapath+"/"+chain):
112 112 continue
113   -
  113 +
114 114 # Initialise dictionary per chain
115 115 folders[chain] = {}
116   -
  116 +
117 117 # Read available DST versions
118 118 dstversions = os.listdir(self.datapath+"/"+chain)
119   -
  119 +
120 120 # Loop over DST versions
121 121 for dstver in dstversions:
122   -
  122 +
123 123 # Omit hidden files/folders
124 124 if dstver[0] == "." or not os.path.isdir(self.datapath+"/"+chain+"/"+dstver):
125 125 continue
126   -
  126 +
127 127 # Initialise available cut configurations
128 128 folders[chain][dstver] = []
129 129 configs = os.listdir(self.datapath+"/"+chain+"/"+dstver)
130   -
  130 +
131 131 # Loop over configurations
132 132 for config in configs:
133   -
  133 +
134 134 # Omit hidden files/folders
135 135 if config[0] == "." or not os.path.isdir(self.datapath+"/"+chain+"/"+dstver+"/"+config):
136 136 continue
137   -
  137 +
138 138 # Append configuration to corresponding dictionary
139 139 folders[chain][dstver].append(config)
140   -
  140 +
141 141 # Check for appropriate folder structure
142 142 checked_folders = self.check_structure(folders)
143 143 if not len(checked_folders):
144 144 raise KeyError("Appropriate folder structure not found in $HESSFITS. The correct structure can be found on HESS internal documentation")
145   -
  145 +
146 146 # Create default parfile
147 147 pars = gammalib.GApplicationPars()
148   -
  148 +
149 149 # Create parameters for each available config
150 150 avchains = ""
151 151 for chain in checked_folders:
... ... @@ -168,38 +168,38 @@ class cshessobs(ctools.cscript):
168 168 pars.append(gammalib.GApplicationPar("usetrig","b","h","yes","yes|no","","Pre-estimate background normalisation? Recommended for binned analysis"))
169 169 pars.append_standard()
170 170 pars.save(parfile)
171   -
  171 +
172 172 # Return
173 173 return
174   -
  174 +
175 175 def get_parameters(self):
176 176 """
177 177 Get parameters from parfile and setup the observation.
178 178 """
179 179 # Get parameters
180 180 self.runlistfile = self["runlistfile"].filename()
181   -
  181 +
182 182 # Read models if user input
183 183 self.inmodel = self["inmodel"].filename()
184 184 if not (self.inmodel == "NONE" or self.inmodel == ""):
185 185 self.inmodels = gammalib.GModels(self.inmodel)
186   -
  186 +
187 187 # Observation outfile
188 188 self.outobs = self["outobs"].filename()
189   -
  189 +
190 190 # Output model file
191 191 self.outmodel = self["outmodel"].filename()
192   -
  192 +
193 193 # Read hidden parameters
194 194 self.nodes = self["nodes"].integer()
195 195 self.usetrig = self["usetrig"].boolean()
196 196 self.datapath = self["datapath"].string()
197   -
  197 +
198 198 # Read run configuration parameters
199 199 self.chain = self["chain"].string()
200 200 self.dstver = self[self.chain+"_dstver"].string()
201 201 self.config = self[self.chain+"_"+self.dstver+"_config"].string()
202   -
  202 +
203 203 # Build directory name prefix
204 204 self.directory = self.datapath+"/"
205 205 self.directory += self.chain +"/"
... ... @@ -208,22 +208,22 @@ class cshessobs(ctools.cscript):
208 208  
209 209 if not os.path.isdir(self.directory):
210 210 raise gammalib.Exception.runtime_error("Folder structure not found. Please check your HESSFITS environment variable")
211   -
  211 +
212 212 self.m_log = False # Logging in client tools
213 213 self.m_debug = False # Debugging in client tools
214 214 self.m_clobber = self["clobber"].boolean()
215   -
  215 +
216 216 # Return
217 217 return
218   -
  218 +
219 219 def check_structure(self,folders):
220   -
  220 +
221 221 # Initialise valid folders
222 222 newfolders = {}
223   -
  223 +
224 224 # Check if we have found possible chains
225 225 if len(folders):
226   -
  226 +
227 227 # Loop over chains
228 228 for chain in folders:
229 229  
... ... @@ -232,7 +232,7 @@ class cshessobs(ctools.cscript):
232 232  
233 233 # Loop over configs:
234 234 for config in folders[chain][dst]:
235   -
  235 +
236 236 # check if config is valid
237 237 path = self.datapath+"/"+chain+"/"+dst+"/"+config
238 238 if os.path.isdir(path) and len(glob.glob(path+"/run*")):
... ... @@ -241,25 +241,25 @@ class cshessobs(ctools.cscript):
241 241 if not newfolders[chain].has_key(dst):
242 242 newfolders[chain][dst] = []
243 243 newfolders[chain][dst].append(config)
244   -
  244 +
245 245 # Return valid configs
246 246 return newfolders
247   -
  247 +
248 248 # Create instrument background model
249 249 def hess_inst_background(self,run,aefffile,trgrate=0.0):
250   -
  250 +
251 251 # Initialise return spectrum
252 252 spec = 0
253   -
  253 +
254 254 # Initialise default value
255 255 value = 1.0
256   -
  256 +
257 257 if self.usetrig and trgrate > 0.0:
258 258 # Use value from lookup
259 259 value = 0.00451521513379 * trgrate + 0.162316854264
260   -
  260 +
261 261 if self.nodes <= 1:
262   -
  262 +
263 263 # Use a constant scaling function
264 264 spec = gammalib.GModelSpectralConst()
265 265 spec["Value"].min(0.1)
... ... @@ -269,9 +269,9 @@ class cshessobs(ctools.cscript):
269 269 spec["Value"].fix()
270 270 else:
271 271 spec["Value"].free()
272   -
  272 +
273 273 elif self.nodes == 2:
274   -
  274 +
275 275 # User a power law scaling function
276 276 e = gammalib.GEnergy(1.0,"TeV")
277 277 spec = gammalib.GModelSpectralPlaw(1.0,0.0,e)
... ... @@ -285,7 +285,7 @@ class cshessobs(ctools.cscript):
285 285 spec[1].max(2.0)
286 286 spec[1].scale(1.0)
287 287 else:
288   -
  288 +
289 289 # Use several energy-dependent scaling factors
290 290 fits = gammalib.GFits(aefffile)
291 291 emin = fits["EFFECTIVE AREA"].real("LO_THRES")
... ... @@ -301,39 +301,39 @@ class cshessobs(ctools.cscript):
301 301 elif "Intensity" in par.name():
302 302 par.min(0.01)
303 303 par.max(10.0)
304   -
  304 +
305 305 # Create background model instance
306 306 bck = gammalib.GCTAModelIrfBackground(spec)
307   -
  307 +
308 308 # Assign specific run id
309 309 bck.ids(str(int(run)))
310   -
  310 +
311 311 # Assign instrument
312 312 bck.instruments("HESS")
313   -
  313 +
314 314 # Set name (arbitrary)
315 315 bck.name("bkg_"+run)
316   -
  316 +
317 317 # Turn off TS calculation for background model
318 318 bck.tscalc(False)
319   -
  319 +
320 320 # Return model
321 321 return bck
322   -
323   -
  322 +
  323 +
324 324 def execute(self):
325 325 """
326 326 Execute the script.
327 327 """
328 328 # Run the script
329 329 self.run()
330   -
  330 +
331 331 # Save residual map
332 332 self.save(self.outobs, self.m_clobber)
333   -
  333 +
334 334 # Return
335 335 return
336   -
  336 +
337 337 # Get runfolder from runnumber (subject to defined folder structure)
338 338 def get_runfolder(self,runnr):
339 339  
... ... @@ -346,13 +346,13 @@ class cshessobs(ctools.cscript):
346 346 end = str(int_end)
347 347 if len(end) == 5:
348 348 end = "0"+ end
349   -
  349 +
350 350 runfold = runnr
351 351 if len(runfold) == 5:
352 352 runfold = "0"+runfold
353   -
  353 +
354 354 return self.directory+"/run"+start+"-"+end+"/run"+runfold+"/"
355   -
  355 +
356 356 def run(self):
357 357 """
358 358 Run the script.
... ... @@ -360,7 +360,7 @@ class cshessobs(ctools.cscript):
360 360 # Switch screen logging on in debug mode
361 361 if self.logDebug():
362 362 self.log.cout(True)
363   -
  363 +
364 364 # Get parameters
365 365 self.get_parameters()
366 366  
... ... @@ -368,7 +368,7 @@ class cshessobs(ctools.cscript):
368 368 if self.logTerse():
369 369 self.log_parameters()
370 370 self.log("\n")
371   -
  371 +
372 372 # Initialise empty models
373 373 self.models = gammalib.GModels()
374 374  
... ... @@ -376,7 +376,7 @@ class cshessobs(ctools.cscript):
376 376 self.xml = gammalib.GXml()
377 377 self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\""))
378 378 lib = self.xml.element("observation_list", 0)
379   -
  379 +
380 380 # Read runlist from file
381 381 runlist = []
382 382 runfile = open(self.runlistfile)
... ... @@ -390,30 +390,30 @@ class cshessobs(ctools.cscript):
390 390 else:
391 391 runlist.append('')
392 392 runfile.close()
393   -
  393 +
394 394 # Loop over runs
395 395 # Log output
396 396 if self.logTerse():
397 397 self.log("\n")
398 398 self.log.header1("Looping over runs")
399   -
  399 +
400 400 for runnr in runlist:
401   -
  401 +
402 402 # Unify run number string length
403 403 run = runnr
404 404 if len(run) == 5:
405 405 run = "0"+run
406   -
  406 +
407 407 # Get folder in which run files should be located
408 408 runfolder = self.get_runfolder(run)
409   -
  409 +
410 410 # construct filenames
411 411 eventfile = runfolder + "hess_events_"+run+".fits.gz"
412 412 aefffile = runfolder + "hess_aeff_2d_"+run+".fits.gz"
413 413 psffile = runfolder + "hess_psf_king_"+run+".fits.gz"
414 414 edispfile = ""#runfolder + "hess_edisp_2d_"+run+".fits.gz"
415 415 bgfile = runfolder + "hess_bkg_offruns_"+run+".fits.gz"
416   -
  416 +
417 417 # Check for existence of files
418 418 msg = ""
419 419 skip = False
... ... @@ -434,7 +434,7 @@ class cshessobs(ctools.cscript):
434 434 skip = True
435 435 else:
436 436 msg = "Adding run "+str(int(run))
437   -
  437 +
438 438 # Log ouput and skip run if not available
439 439 if self.logTerse():
440 440 self.log(msg+"\n")
... ... @@ -446,11 +446,11 @@ class cshessobs(ctools.cscript):
446 446 self.log(" Point spread function: "+psffile+"\n")
447 447 self.log(" Energy dispersion: "+edispfile+"\n")
448 448 self.log(" Background file: "+bgfile+"\n")
449   -
  449 +
450 450 # Open fits file to determine the observation name
451 451 fits = gammalib.GFits(eventfile)
452 452 object_name = fits["EVENTS"].string("OBJECT")
453   -
  453 +
454 454 # Retrieve trigger rate if available
455 455 trgrate = 0.0
456 456 if fits["EVENTS"].has_card("ZTRGRATE"):
... ... @@ -458,7 +458,7 @@ class cshessobs(ctools.cscript):
458 458  
459 459 # Close FITS file
460 460 fits.close()
461   -
  461 +
462 462 # Append observation to XML and set attributes
463 463 obs = lib.append("observation");
464 464 obs.attribute("name", object_name);
... ... @@ -469,29 +469,29 @@ class cshessobs(ctools.cscript):
469 469 par = gammalib.GXmlElement("parameter name=\"EventList\"")
470 470 par.attribute("file", eventfile);
471 471 obs.append(par)
472   -
  472 +
473 473 # Append effective area
474 474 aeff = gammalib.GXmlElement("parameter name=\"EffectiveArea\"")
475 475 aeff.attribute("file",aefffile)
476   -
  476 +
477 477 # Append PSF
478 478 psf = gammalib.GXmlElement("parameter name=\"PointSpreadFunction\"")
479 479 psf.attribute("file",psffile)
480   -
  480 +
481 481 # Append energy dispersion
482 482 edisp = gammalib.GXmlElement("parameter name=\"EnergyDispersion\"")
483 483 edisp.attribute("file",edispfile)
484   -
  484 +
485 485 # Append background
486 486 bck = gammalib.GXmlElement("parameter name=\"Background\"")
487 487 bck.attribute("file",bgfile)
488   -
  488 +
489 489 # assign events and IRFs to observation
490 490 obs.append(aeff)
491 491 obs.append(psf)
492 492 obs.append(edisp)
493 493 obs.append(bck)
494   -
  494 +
495 495 # Append instrumental background model
496 496 self.models.append(self.hess_inst_background(run,aefffile,trgrate))
497 497  
... ... @@ -500,7 +500,7 @@ class cshessobs(ctools.cscript):
500 500 if self.logTerse():
501 501 self.log("\n")
502 502 self.log.header1("Appending models")
503   -
  503 +
504 504 # Loop over input models
505 505 for model in self.inmodels:
506 506 if self.logTerse():
... ... @@ -509,12 +509,12 @@ class cshessobs(ctools.cscript):
509 509  
510 510 # Return
511 511 return
512   -
  512 +
513 513 def save(self,outfile,clobber):
514   -
  514 +
515 515 # Save observation XML file
516 516 self.xml.save(outfile)
517   -
  517 +
518 518 # Save model XML file
519 519 self.models.save(self.outmodel)
520 520  
... ... @@ -528,10 +528,9 @@ if __name__ == &#39;__main__&#39;:
528 528 """
529 529 # Create instance of application
530 530 app = cshessobs(sys.argv)
531   -
  531 +
532 532 # Open logfile
533 533 app.logFileOpen()
534   -
  534 +
535 535 # Execute application
536 536 app.execute()
537   -
... ...
cscripts/cslightcrv.py
  View file @5c53c14
... ... @@ -36,14 +36,14 @@ class cslightcrv(ctools.cscript):
36 36 """
37 37 Constructor.
38 38 """
39   -
  39 +
40 40 # Set name
41 41 self.name = "cslightcrv"
42 42 self.version = "1.0.0"
43   -
  43 +
44 44 # Initialise some members
45 45 self.obs = None
46   -
  46 +
47 47 # Initialise some members
48 48 if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations):
49 49 self.obs = argv[0]
... ... @@ -52,7 +52,7 @@ class cslightcrv(ctools.cscript):
52 52 self.obs = gammalib.GObservations()
53 53 self.obs.clear()
54 54 self.m_outfile = ""
55   -
  55 +
56 56 # Make sure that parfile exists
57 57 file = self.parfile()
58 58  
... ... @@ -70,7 +70,7 @@ class cslightcrv(ctools.cscript):
70 70  
71 71 # Return
72 72 return
73   -
  73 +
74 74 def __del__(self):
75 75 """
76 76 Destructor.
... ... @@ -78,7 +78,7 @@ class cslightcrv(ctools.cscript):
78 78 # Write separator into logger
79 79 if self.logTerse():
80 80 self.log("\n")
81   -
  81 +
82 82 # Return
83 83 return
84 84  
... ... @@ -89,13 +89,13 @@ class cslightcrv(ctools.cscript):
89 89 """
90 90 # Set parfile name
91 91 parfile = self.name+".par"
92   -
  92 +
93 93 try:
94 94 pars = gammalib.GApplicationPars(parfile)
95 95 except:
96 96 # Signal if parfile was not found
97 97 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
98   -
  98 +
99 99 # Create default parfile
100 100 pars = gammalib.GApplicationPars()
101 101 pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file"))
... ... @@ -129,12 +129,12 @@ class cslightcrv(ctools.cscript):
129 129 pars.append_standard()
130 130 pars.append(gammalib.GApplicationPar("logfile","f","h","cslightcrv.log","","","Log filename"))
131 131 pars.save(parfile)
132   -
  132 +
133 133 # Return
134 134 return
135 135  
136 136 def create_tbounds(self):
137   -
  137 +
138 138 # Create time bin container
139 139 self.m_tbins = gammalib.GGti()
140 140  
... ... @@ -172,7 +172,7 @@ class cslightcrv(ctools.cscript):
172 172  
173 173 # Return
174 174 return
175   -
  175 +
176 176 def get_parameters(self):
177 177 """
178 178 Get parameters from parfile and setup the observation.
... ... @@ -181,14 +181,14 @@ class cslightcrv(ctools.cscript):
181 181 if self.obs == None or self.obs.size() == 0:
182 182 self.require_inobs("cslightcrv::get_parameters()")
183 183 self.obs = self.get_observations()
184   -
  184 +
185 185 # Set models if we have none
186 186 if self.obs.models().size() == 0:
187 187 self.obs.models(self["inmodel"].filename())
188   -
  188 +
189 189 # Get source name
190 190 self.m_srcname = self["srcname"].string()
191   -
  191 +
192 192 # Get time boundaries
193 193 self.create_tbounds()
194 194  
... ... @@ -215,32 +215,32 @@ class cslightcrv(ctools.cscript):
215 215  
216 216 # Read other parameters
217 217 self.m_outfile = self["outfile"].filename()
218   -
  218 +
219 219 # Get other parameeters
220 220 self.m_calc_ulimit = self["calc_ulim"].boolean()
221 221 self.m_calc_ts = self["calc_ts"].boolean()
222 222 self.m_fix_bkg = self["fix_bkg"].boolean()
223 223 self.m_fix_srcs = self["fix_srcs"].boolean()
224   -
  224 +
225 225 # Set some fixed parameters
226 226 self.m_log = False # Logging in client tools
227 227 self.m_chatter = self["chatter"].integer()
228 228 self.m_clobber = self["clobber"].boolean()
229 229 self.m_debug = self["debug"].boolean()
230   -
  230 +
231 231 # Return
232 232 return
233   -
  233 +
234 234 def models(self, models):
235 235 """
236 236 Set model.
237 237 """
238 238 # Copy models
239 239 self.obs.models(models.clone())
240   -
  240 +
241 241 # Return
242 242 return
243   -
  243 +
244 244 def execute(self):
245 245 """
246 246 Execute the script.
... ... @@ -250,7 +250,7 @@ class cslightcrv(ctools.cscript):
250 250  
251 251 # Save residual map
252 252 self.fits.saveto(self.m_outfile, self.m_clobber)
253   -
  253 +
254 254 # Return
255 255 return
256 256  
... ... @@ -264,12 +264,12 @@ class cslightcrv(ctools.cscript):
264 264  
265 265 # Get parameters
266 266 self.get_parameters()
267   -
  267 +
268 268 # Write input parameters into logger
269 269 if self.logTerse():
270 270 self.log_parameters()
271 271 self.log("\n")
272   -
  272 +
273 273 # Write observation into logger
274 274 if self.logTerse():
275 275 self.log("\n")
... ... @@ -284,51 +284,51 @@ class cslightcrv(ctools.cscript):
284 284  
285 285 # Adjust model parameters dependent on input user parameters
286 286 for model in self.obs.models():
287   -
  287 +
288 288 # Set TS flag for all models to false.
289 289 # Source of interest will be set to true later
290 290 model.tscalc(False)
291   -
  291 +
292 292 # Log model name
293 293 if self.logExplicit():
294 294 self.log.header3(model.name())
295   -
  295 +
296 296 # Deal with the source of interest
297 297 if model.name() == self.m_srcname:
298 298 if self.m_calc_ts:
299 299 model.tscalc(True)
300   -
  300 +
301 301 elif self.m_fix_bkg and not model.classname() == "GModelSky":
302 302 for par in model:
303 303 if par.is_free() and self.logExplicit():
304 304 self.log(" Fixing \""+par.name()+"\"\n")
305 305 par.fix()
306   -
  306 +
307 307 elif self.m_fix_srcs and model.classname() == "GModelSky":
308 308 for par in model:
309 309 if par.is_free() and self.logExplicit():
310 310 self.log(" Fixing \""+par.name()+"\"\n")
311 311 par.fix()
312   -
  312 +
313 313 # Write header
314 314 if self.logTerse():
315 315 self.log("\n")
316 316 self.log.header1("Generate lightcurve")
317   -
  317 +
318 318 # Initialise FITS Table with extension "LIGHTCURVE"
319 319 table = gammalib.GFitsBinTable(self.m_tbins.size())
320 320 table.extname("LIGHTCURVE")
321   -
  321 +
322 322 # Add Header for compatibility with gammalib.GMWLSpectrum
323 323 table.card("INSTRUME", "CTA", "Name of Instrument")
324 324 table.card("TELESCOP", "CTA", "Name of Telescope")
325   -
  325 +
326 326 # Create FITS table columns
327 327 MJD = gammalib.GFitsTableDoubleCol("MJD", self.m_tbins.size())
328 328 MJD.unit("days")
329 329 e_MJD = gammalib.GFitsTableDoubleCol("e_MJD", self.m_tbins.size())
330 330 e_MJD.unit("days")
331   -
  331 +
332 332 # Create a FITS column for every free parameter
333 333 columns = []
334 334 for par in self.obs.models()[self.m_srcname]:
... ... @@ -339,7 +339,7 @@ class cslightcrv(ctools.cscript):
339 339 e_col = gammalib.GFitsTableDoubleCol("e_"+par.name(), self.m_tbins.size())
340 340 e_col.unit(par.unit())
341 341 columns.append(e_col)
342   -
  342 +
343 343 # Create TS and upper limit columns
344 344 TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_tbins.size())
345 345 ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_tbins.size())
... ... @@ -356,7 +356,7 @@ class cslightcrv(ctools.cscript):
356 356 # Get time boundaries
357 357 tmin = self.m_tbins.tstart(i)
358 358 tmax = self.m_tbins.tstop(i)
359   -
  359 +
360 360 # Compute time bin center and time width
361 361 tmean = (tmin + tmax)
362 362 tmean *= 0.5
... ... @@ -366,11 +366,11 @@ class cslightcrv(ctools.cscript):
366 366 # Store time as MJD
367 367 MJD[i] = tmean.mjd()
368 368 e_MJD[i] = twidth.days()
369   -
  369 +
370 370 # Log information
371 371 if self.logExplicit():
372 372 self.log.header3("Selecting events")
373   -
  373 +
374 374 # Select events
375 375 select = ctools.ctselect(self.obs)
376 376 select["emin"] = self.m_emin
... ... @@ -384,14 +384,14 @@ class cslightcrv(ctools.cscript):
384 384  
385 385 # Retrieve observation
386 386 obs = select.obs()
387   -
  387 +
388 388 # Binned analysis
389 389 if self.m_binned:
390 390  
391 391 # Header
392 392 if self.logTerse():
393 393 self.log.header3("Binning events")
394   -
  394 +
395 395 # Bin events
396 396 bin = ctools.ctbin(select.obs())
397 397 bin["usepnt"] = False
... ... @@ -407,11 +407,11 @@ class cslightcrv(ctools.cscript):
407 407 bin["coordsys"] = self.m_coordsys
408 408 bin["proj"] = self.m_proj
409 409 bin.run()
410   -
  410 +
411 411 # Header
412 412 if self.logTerse():
413 413 self.log.header3("Creating exposure cube")
414   -
  414 +
415 415 # Create exposure cube
416 416 expcube = ctools.ctexpcube(select.obs())
417 417 expcube["incube"] = "NONE"
... ... @@ -428,11 +428,11 @@ class cslightcrv(ctools.cscript):
428 428 expcube["coordsys"] = self.m_coordsys
429 429 expcube["proj"] = self.m_proj
430 430 expcube.run()
431   -
  431 +
432 432 # Header
433 433 if self.logTerse():
434 434 self.log.header3("Creating PSF cube")
435   -
  435 +
436 436 # Create psf cube
437 437 psfcube = ctools.ctpsfcube(select.obs())
438 438 psfcube["incube"] = "NONE"
... ... @@ -449,11 +449,11 @@ class cslightcrv(ctools.cscript):
449 449 psfcube["coordsys"] = self.m_coordsys
450 450 psfcube["proj"] = self.m_proj
451 451 psfcube.run()
452   -
  452 +
453 453 # Header
454 454 if self.logTerse():
455 455 self.log.header3("Creating background cube")
456   -
  456 +
457 457 # Create background cube
458 458 bkgcube = ctools.ctbkgcube(select.obs())
459 459 bkgcube["incube"] = "NONE"
... ... @@ -470,37 +470,37 @@ class cslightcrv(ctools.cscript):
470 470 bkgcube["coordsys"] = self.m_coordsys
471 471 bkgcube["proj"] = self.m_proj
472 472 bkgcube.run()
473   -
  473 +
474 474 # Set new binned observation
475 475 obs = bin.obs()
476   -
  476 +
477 477 # Set precomputed binned response
478 478 obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube())
479 479  
480 480 # Get new models
481 481 models = bkgcube.models()
482   -
  482 +
483 483 # Fix background models if required
484 484 if self.m_fix_bkg:
485 485 for model in models:
486 486 if not model.classname() == "GModelSky":
487 487 for par in model:
488 488 par.fix()
489   -
  489 +
490 490 # Set new models to binned observation
491 491 obs.models(models)
492   -
  492 +
493 493 # Header
494 494 if self.logTerse():
495 495 self.log.header3("Performing fit")
496   -
  496 +
497 497 # Likelihood
498 498 like = ctools.ctlike(obs)
499 499 like.run()
500   -
  500 +
501 501 # Skip bin if no event was present
502 502 if like.obs().logL() == 0.0:
503   -
  503 +
504 504 # Log information
505 505 if self.logTerse():
506 506 self.log("No event in this time bin. Bin is skipped\n")
... ... @@ -511,7 +511,7 @@ class cslightcrv(ctools.cscript):
511 511 TSvalues[i] = 0.0
512 512 ulim_values[i] = 0.0
513 513 continue
514   -
  514 +
515 515 # Get results
516 516 fitted_models = like.obs().models()
517 517 source = fitted_models[self.m_srcname]
... ... @@ -519,16 +519,16 @@ class cslightcrv(ctools.cscript):
519 519 # Calculate Upper Limit
520 520 ulimit_value = -1.0
521 521 if self.m_calc_ulimit:
522   -
  522 +
523 523 # Logging information
524 524 if self.logTerse():
525 525 self.log.header3("Computing upper limit")
526   -
  526 +
527 527 # Create upper limit object
528 528 ulimit = ctools.ctulimit(like.obs())
529 529 ulimit["srcname"] = self.m_srcname
530 530 ulimit["eref"] = 1.0
531   -
  531 +
532 532 # Try to run upper limit and catch exceptions
533 533 try:
534 534 ulimit.run()
... ... @@ -537,26 +537,26 @@ class cslightcrv(ctools.cscript):
537 537 if self.logTerse():
538 538 self.log("Upper limit calculation failed\n")
539 539 ulimit_value = -1.0
540   -
  540 +
541 541 # Get TS value
542 542 TS = -1.0
543 543 if self.m_calc_ts:
544 544 TS = source.ts()
545   -
  545 +
546 546 # Set values for storage
547 547 TSvalues[i] = TS
548   -
  548 +
549 549 # Set FITS column values
550 550 for col in columns:
551 551 if "e_" == col.name()[:2]:
552 552 col[i] = source.spectral()[col.name()[2:]].error()
553 553 else:
554 554 col[i] = source.spectral()[col.name()].value()
555   -
  555 +
556 556 # Store upper limit value if available
557 557 if ulimit_value > 0.0:
558 558 ulim_values[i] = ulimit_value
559   -
  559 +
560 560 # Log information
561 561 if self.logExplicit():
562 562 self.log.header3("Results of bin "+str(i)+": MJD "+str(tmin.mjd())+"-"+str(tmax.mjd()))
... ... @@ -580,14 +580,14 @@ class cslightcrv(ctools.cscript):
580 580 table.append(col)
581 581 table.append(TSvalues)
582 582 table.append(ulim_values)
583   -
  583 +
584 584 # Create the FITS file now
585 585 self.fits = gammalib.GFits()
586 586 self.fits.append(table)
587   -
  587 +
588 588 # Return
589 589 return
590   -
  590 +
591 591  
592 592 # ======================== #
593 593 # Main routine entry point #
... ... @@ -598,9 +598,9 @@ if __name__ == &#39;__main__&#39;:
598 598 """
599 599 # Create instance of application
600 600 app = cslightcrv(sys.argv)
601   -
  601 +
602 602 # Open logfile
603 603 app.logFileOpen()
604   -
  604 +
605 605 # Execute application
606 606 app.execute()
... ...
cscripts/csobsdef.py
  View file @5c53c14
... ... @@ -33,7 +33,7 @@ class csobsdef(ctools.cscript):
33 33 file with header keywords in the first row followed by a list of
34 34 pointings (one pointing per row). The following header keywords
35 35 are supported (case sensitive, column order irrelevant):
36   -
  36 +
37 37 name - Observation name string
38 38 id - Unique observation identifier string
39 39 ra - Right Ascension of pointing (deg)
... ... @@ -47,14 +47,14 @@ class csobsdef(ctools.cscript):
47 47 deadc - Deadtime correction factor [0-1]
48 48 caldb - Calibration database
49 49 irf - Response function name
50   -
  50 +
51 51 Only the pairs (ra,dec) or (lon,lat) are mandatory header keywords.
52 52 All other keywords are optional and can be specified when calling
53 53 csobsdef as user parameters. The only exception is the "duration"
54 54 keyword that will automatically queried.
55   -
  55 +
56 56 Here some usage examples:
57   -
  57 +
58 58 csobsdef
59 59 Creates minimal observation definition file.
60 60  
... ... @@ -85,7 +85,7 @@ class csobsdef(ctools.cscript):
85 85 self.clobber = True
86 86 self.debug = False
87 87 self.read_ahead = False
88   -
  88 +
89 89 # Make sure that parfile exists
90 90 file = self.parfile()
91 91  
... ... @@ -103,7 +103,7 @@ class csobsdef(ctools.cscript):
103 103  
104 104 # Return
105 105 return
106   -
  106 +
107 107 def __del__(self):
108 108 """
109 109 Destructor.
... ... @@ -125,7 +125,7 @@ class csobsdef(ctools.cscript):
125 125 except:
126 126 # Signal that parfile was not found
127 127 print("Parfile \""+parfile+"\" not found. Create default parfile.")
128   -
  128 +
129 129 # Create default parfile
130 130 pars = gammalib.GApplicationPars()
131 131 pars.append(gammalib.GApplicationPar("inpnt","f","a","NONE","","","Pointing definition file"))
... ... @@ -140,10 +140,10 @@ class csobsdef(ctools.cscript):
140 140 pars.append_standard()
141 141 pars.append(gammalib.GApplicationPar("logfile","f","h","csobsdef.log","","","Log filename"))
142 142 pars.save(parfile)
143   -
  143 +
144 144 # Return
145 145 return
146   -
  146 +
147 147 def get_parameters(self):
148 148 """
149 149 Get parameters from parfile.
... ... @@ -162,20 +162,20 @@ class csobsdef(ctools.cscript):
162 162  
163 163 # Return
164 164 return
165   -
  165 +
166 166 def execute(self):
167 167 """
168 168 Execute the script.
169 169 """
170 170 # Set read ahead flag
171 171 self.read_ahead = True
172   -
  172 +
173 173 # Run the script
174 174 self.run()
175 175  
176 176 # Save observation definition file
177 177 self.obs.save(self.outfile)
178   -
  178 +
179 179 # Return
180 180 return
181 181  
... ... @@ -189,7 +189,7 @@ class csobsdef(ctools.cscript):
189 189  
190 190 # Get parameters
191 191 self.get_parameters()
192   -
  192 +
193 193 # Write input parameters into logger
194 194 if self.logTerse():
195 195 self.log_parameters()
... ... @@ -325,13 +325,13 @@ class csobsdef(ctools.cscript):
325 325 # Create an empty event list
326 326 list = gammalib.GCTAEventList()
327 327 list.gti(gti)
328   -
  328 +
329 329 # Set optional information
330 330 if has_ebounds:
331 331 list.ebounds(ebounds)
332 332 if has_roi:
333 333 list.roi(roi)
334   -
  334 +
335 335 # Attach event list to CTA observation
336 336 obs.events(list)
337 337  
... ... @@ -358,10 +358,10 @@ class csobsdef(ctools.cscript):
358 358 # Create CTA response
359 359 response = gammalib.GCTAResponseIrf()
360 360 response.read(xml)
361   -
  361 +
362 362 # Attach response to observation
363 363 obs.response(response)
364   -
  364 +
365 365 # Return observation
366 366 return obs
367 367  
... ... @@ -375,10 +375,9 @@ if __name__ == &#39;__main__&#39;:
375 375 """
376 376 # Create instance of application
377 377 app = csobsdef(sys.argv)
378   -
  378 +
379 379 # Open logfile
380 380 app.logFileOpen()
381   -
  381 +
382 382 # Execute application
383 383 app.execute()
384   -
385 384 \ No newline at end of file
... ...
cscripts/cspull.py
  View file @5c53c14
... ... @@ -42,7 +42,7 @@ class cspull(ctools.cscript):
42 42 # Set name
43 43 self.name = "cspull"
44 44 self.version = "1.0.0"
45   -
  45 +
46 46 # Initialise some members
47 47 self.obs = None
48 48 self.model = None
... ... @@ -53,7 +53,7 @@ class cspull(ctools.cscript):
53 53 self.m_psfcube = None
54 54 self.m_bckcube = None
55 55 self.m_stackmodels = None
56   -
  56 +
57 57 # Make sure that parfile exists
58 58 file = self.parfile()
59 59  
... ... @@ -71,7 +71,7 @@ class cspull(ctools.cscript):
71 71  
72 72 # Return
73 73 return
74   -
  74 +
75 75 def __del__(self):
76 76 """
77 77 Destructor.
... ... @@ -79,7 +79,7 @@ class cspull(ctools.cscript):
79 79 # Write separator into logger
80 80 if self.logTerse():
81 81 self.log("\n")
82   -
  82 +
83 83 # Return
84 84 return
85 85  
... ... @@ -90,13 +90,13 @@ class cspull(ctools.cscript):
90 90 """
91 91 # Set parfile name
92 92 parfile = self.name+".par"
93   -
  93 +
94 94 try:
95 95 pars = gammalib.GApplicationPars(parfile)
96 96 except:
97 97 # Signal if parfile was not found
98 98 print("Parfile "+parfile+" not found. Create default parfile.")
99   -
  99 +
100 100 # Create default parfile
101 101 pars = gammalib.GApplicationPars()
102 102 pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file"))
... ... @@ -126,29 +126,29 @@ class cspull(ctools.cscript):
126 126 pars.append_standard()
127 127 pars.append(gammalib.GApplicationPar("logfile","f","h","cspull.log","","","Log filename"))
128 128 pars.save(parfile)
129   -
  129 +
130 130 # Return
131 131 return
132   -
  132 +
133 133 def get_parameters(self):
134 134 """
135 135 Get parameters from parfile and setup the observation.
136 136 """
137 137 # Get parameters
138   -
  138 +
139 139 # Set observation if not done before
140 140 if self.obs == None or self.obs.size() == 0:
141 141 self.obs = self.get_observations()
142   -
  142 +
143 143 # Check for requested pattern and use above
144 144 # observation parameters to set wobble pattern
145 145 self.m_pattern = self["pattern"].string()
146 146 if self.m_pattern == "four":
147 147 self.obs = self.set_obs()
148   -
  148 +
149 149 # Get number of energy bins
150 150 self.m_enumbins = self["enumbins"].integer()
151   -
  151 +
152 152 # Read parameters for binned if requested
153 153 if not self.m_enumbins == 0:
154 154 self.m_npix = self["npix"].integer()
... ... @@ -157,11 +157,11 @@ class cspull(ctools.cscript):
157 157 # Set dummy values (required by obsutils)
158 158 self.m_npix = 0
159 159 self.m_binsz = 0.0
160   -
  160 +
161 161 # Set models if we have none
162 162 if self.obs.models().size() == 0:
163 163 self.obs.models(self["inmodel"].filename())
164   -
  164 +
165 165 # Read other parameters
166 166 self.m_outfile = self["outfile"].filename()
167 167 self.m_ntrials = self["ntrials"].integer()
... ... @@ -176,7 +176,7 @@ class cspull(ctools.cscript):
176 176  
177 177 # Return
178 178 return
179   -
  179 +
180 180  
181 181 def models(self, models):
182 182 """
... ... @@ -184,10 +184,10 @@ class cspull(ctools.cscript):
184 184 """
185 185 # Copy models
186 186 self.obs.models(models.copy())
187   -
  187 +
188 188 # Return
189 189 return
190   -
  190 +
191 191  
192 192 def execute(self):
193 193 """
... ... @@ -195,11 +195,11 @@ class cspull(ctools.cscript):
195 195 """
196 196 # Run the script
197 197 self.run()
198   -
  198 +
199 199 # Return
200 200 return
201 201  
202   -
  202 +
203 203 def set_stacked_irf(self):
204 204 """
205 205 For stacked analysis prepare stacked irfs
... ... @@ -227,7 +227,7 @@ class cspull(ctools.cscript):
227 227 # Notify exposure computation
228 228 if self.logTerse():
229 229 self.log("Computed exposure cube\n")
230   -
  230 +
231 231 # Get stacked Psf
232 232 psfcube = ctools.ctpsfcube(self.obs)
233 233 psfcube["incube"] = "NONE"
... ... @@ -246,7 +246,7 @@ class cspull(ctools.cscript):
246 246 # Notify Psf computation
247 247 if self.logTerse():
248 248 self.log("Computed Psf cube\n")
249   -
  249 +
250 250 # Get stacked background
251 251 bkgcube = ctools.ctbkgcube(self.obs)
252 252 bkgcube["incube"] = "NONE"
... ... @@ -265,17 +265,17 @@ class cspull(ctools.cscript):
265 265 # Notify background cube computation
266 266 if self.logTerse():
267 267 self.log("Computed background cube\n")
268   -
  268 +
269 269 # Store results
270 270 self.m_exposure = expcube.expcube().copy()
271 271 self.m_psfcube = psfcube.psfcube().copy()
272 272 self.m_bckcube = bkgcube.bkgcube().copy()
273 273 self.m_stackmodels = bkgcube.models().copy()
274   -
  274 +
275 275 # Return
276 276 return
277   -
278   -
  277 +
  278 +
279 279 def run(self):
280 280 """
281 281 Run the script.
... ... @@ -286,12 +286,12 @@ class cspull(ctools.cscript):
286 286  
287 287 # Get parameters
288 288 self.get_parameters()
289   -
  289 +
290 290 # Write input parameters into logger
291 291 if self.logTerse():
292 292 self.log_parameters()
293 293 self.log("\n")
294   -
  294 +
295 295 # Write observation into logger
296 296 if self.logTerse():
297 297 self.log("\n")
... ... @@ -310,10 +310,10 @@ class cspull(ctools.cscript):
310 310  
311 311 # Loop over trials
312 312 for seed in range(self.m_ntrials):
313   -
  313 +
314 314 # Make a trial
315 315 result = self.trial(seed)
316   -
  316 +
317 317 # Write out result immediately
318 318 if seed == 0:
319 319 file = open(self.m_outfile, 'w')
... ... @@ -324,23 +324,23 @@ class cspull(ctools.cscript):
324 324 writer = csv.DictWriter(file, result['colnames'])
325 325 writer.writerow(result['values'])
326 326 file.close()
327   -
  327 +
328 328 # Return
329 329 return
330   -
  330 +
331 331  
332 332 def set_obs(self):
333 333 """
334 334 Returns an observation container with a set of CTA observations.
335   -
  335 +
336 336 Keywords:
337 337 """
338   -
  338 +
339 339 # Setup observation definition list
340 340 obsdeflist = obsutils.set_obs_patterns(self.m_pattern,
341 341 ra=self["ra"].real(), dec=self["dec"].real(),
342 342 offset=self["offset"].real())
343   -
  343 +
344 344 # Create list of observations
345 345 obs = obsutils.set_obs_list(obsdeflist,
346 346 tstart=self["tmin"].real(), duration=self["tmax"].real()-self["tmin"].real(),
... ... @@ -348,7 +348,7 @@ class cspull(ctools.cscript):
348 348 emin=self["emin"].real(), emax=self["emax"].real(),
349 349 rad=self["rad"].real(),
350 350 irf=self["irf"].string(), caldb=self["caldb"].string())
351   -
  351 +
352 352 # Return observation container
353 353 return obs
354 354  
... ... @@ -356,11 +356,11 @@ class cspull(ctools.cscript):
356 356 def trial(self, seed):
357 357 """
358 358 Create the pull for a single trial.
359   -
  359 +
360 360 Parameters:
361 361 seed - Random number generator seed
362 362 """
363   -
  363 +
364 364 # Write header
365 365 if self.logNormal():
366 366 self.log.header2("Trial "+str(seed+1))
... ... @@ -375,12 +375,12 @@ class cspull(ctools.cscript):
375 375 log=self.m_log,
376 376 debug=self.m_debug,
377 377 chatter=self.m_chatter)
378   -
  378 +
379 379 # If stacked, add stacked responses and model
380 380 if self.obs.size() > 0 and self.m_enumbins > 0:
381 381 obs[0].response(self.m_exposure, self.m_psfcube, self.m_bckcube)
382 382 obs.models(self.m_stackmodels)
383   -
  383 +
384 384 # Determine number of events in simulation
385 385 nevents = 0.0
386 386 for run in obs:
... ... @@ -416,7 +416,7 @@ class cspull(ctools.cscript):
416 416 # Write result header
417 417 if self.logNormal():
418 418 self.log.header3("Pulls")
419   -
  419 +
420 420 # Gather results
421 421 colnames = []
422 422 values = {}
... ... @@ -432,15 +432,15 @@ class cspull(ctools.cscript):
432 432 for k in range(model.size()):
433 433 par = model[k]
434 434 if par.is_free():
435   -
  435 +
436 436 # Set parameter name
437 437 name = model_name+"_"+par.name()
438   -
  438 +
439 439 # Append parameter, Pull_parameter and Unc_parameter
440 440 colnames.append(name)
441 441 colnames.append("Pull_"+name)
442 442 colnames.append("Unc_"+name)
443   -
  443 +
444 444 # Compute pull
445 445 fitted_value = par.value()
446 446 real_value = self.obs.models()[i][k].value()
... ... @@ -449,7 +449,7 @@ class cspull(ctools.cscript):
449 449 pull = (fitted_value - real_value) / error
450 450 else:
451 451 pull = 99.0
452   -
  452 +
453 453 # Store results
454 454 values[name] = fitted_value
455 455 values["Pull_"+name] = pull
... ... @@ -464,10 +464,10 @@ class cspull(ctools.cscript):
464 464 self.log(" +/- ")
465 465 self.log(error)
466 466 self.log(")\n")
467   -
  467 +
468 468 # Bundle together results
469 469 result = {'colnames': colnames, 'values': values}
470   -
  470 +
471 471 # Return
472 472 return result
473 473  
... ... @@ -481,10 +481,9 @@ if __name__ == &#39;__main__&#39;:
481 481 """
482 482 # Create instance of application
483 483 app = cspull(sys.argv)
484   -
  484 +
485 485 # Open logfile
486 486 app.logFileOpen()
487   -
  487 +
488 488 # Execute application
489 489 app.execute()
490   -
... ...
cscripts/csresmap.py
  View file @5c53c14
... ... @@ -36,11 +36,11 @@ class csresmap(ctools.cscript):
36 36 """
37 37 Constructor.
38 38 """
39   -
  39 +
40 40 # Set name
41 41 self.name = "csresmap"
42 42 self.version = "1.0.0"
43   -
  43 +
44 44 # Initialise some members
45 45 self.obs = None
46 46 self.algorithm = "SUB"
... ... @@ -48,7 +48,7 @@ class csresmap(ctools.cscript):
48 48 self.m_modcube = "None"
49 49 self.m_use_maps = False
50 50 self.m_skip_binning = False
51   -
  51 +
52 52 # Initialise some members
53 53 if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations):
54 54 self.obs = argv[0]
... ... @@ -57,7 +57,7 @@ class csresmap(ctools.cscript):
57 57 self.obs = gammalib.GObservations()
58 58 self.obs.clear()
59 59 self.outfile = ""
60   -
  60 +
61 61 # Make sure that parfile exists
62 62 file = self.parfile()
63 63  
... ... @@ -75,7 +75,7 @@ class csresmap(ctools.cscript):
75 75  
76 76 # Return
77 77 return
78   -
  78 +
79 79 def __del__(self):
80 80 """
81 81 Destructor.
... ... @@ -83,7 +83,7 @@ class csresmap(ctools.cscript):
83 83 # Write separator into logger
84 84 if self.logTerse():
85 85 self.log("\n")
86   -
  86 +
87 87 # Return
88 88 return
89 89  
... ... @@ -94,13 +94,13 @@ class csresmap(ctools.cscript):
94 94 """
95 95 # Set parfile name
96 96 parfile = self.name+".par"
97   -
  97 +
98 98 try:
99 99 pars = gammalib.GApplicationPars(parfile)
100 100 except:
101 101 # Signal if parfile was not found
102 102 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
103   -
  103 +
104 104 # Create default parfile
105 105 pars = gammalib.GApplicationPars()
106 106 pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file"))
... ... @@ -128,10 +128,10 @@ class csresmap(ctools.cscript):
128 128 pars.append_standard()
129 129 pars.append(gammalib.GApplicationPar("logfile","f","h","csresmap.log","","","Log filename"))
130 130 pars.save(parfile)
131   -
  131 +
132 132 # Return
133 133 return
134   -
  134 +
135 135 def get_parameters(self):
136 136 """
137 137 Get parameters from parfile and setup the observation.
... ... @@ -153,28 +153,28 @@ class csresmap(ctools.cscript):
153 153 self.m_modcube = self["modcube"].filename()
154 154 if self.m_modcube != "NONE":
155 155 self.m_use_maps = True
156   -
  156 +
157 157 # If not two maps are given, proceed to set up observation
158 158 if not self.m_use_maps:
159   -
  159 +
160 160 # Set observation if not done before
161 161 if self.obs.size() == 0:
162 162 self.require_inobs("csresmap.get_parameters()")
163 163 self.obs = self.get_observations()
164   -
  164 +
165 165 # Check if we have exactly one binned CTA observation
166 166 if self.obs.size() == 1:
167   -
  167 +
168 168 if self.obs[0].classname() == "GCTAObservation":
169 169 if self.obs[0].eventtype() == "CountsCube":
170 170 # Skip ctbin step later on
171 171 self.m_skip_binning = True
172   -
173   -
  172 +
  173 +
174 174 # Set models if we have none
175 175 if self.obs.models().size() == 0:
176 176 self.obs.models(self["inmodel"].filename())
177   -
  177 +
178 178 # Skip query for spatial parameters if a binning is provided in the observation
179 179 if not self.m_skip_binning:
180 180 # Read other parameters
... ... @@ -189,30 +189,30 @@ class csresmap(ctools.cscript):
189 189 self.m_nxpix = self["nxpix"].integer()
190 190 self.m_nypix = self["nypix"].integer()
191 191 self.m_binsz = self["binsz"].real()
192   -
  192 +
193 193 # Read necessary parameters
194 194 self.m_outfile = self["outmap"].filename()
195 195 self.m_algorithm = self["algorithm"].string()
196   -
  196 +
197 197 # Set some fixed parameters
198 198 self.m_log = False # Logging in client tools
199 199 self.m_chatter = self["chatter"].integer()
200 200 self.m_clobber = self["clobber"].boolean()
201 201 self.m_debug = self["debug"].boolean()
202   -
  202 +
203 203 # Return
204 204 return
205   -
  205 +
206 206 def models(self, models):
207 207 """
208 208 Set model.
209 209 """
210 210 # Copy models
211 211 self.obs.models(models.clone())
212   -
  212 +
213 213 # Return
214 214 return
215   -
  215 +
216 216 def execute(self):
217 217 """
218 218 Execute the script.
... ... @@ -222,7 +222,7 @@ class csresmap(ctools.cscript):
222 222  
223 223 # Save residual map
224 224 self.resmap.save(self.m_outfile, self.m_clobber)
225   -
  225 +
226 226 # Return
227 227 return
228 228  
... ... @@ -236,12 +236,12 @@ class csresmap(ctools.cscript):
236 236  
237 237 # Get parameters
238 238 self.get_parameters()
239   -
  239 +
240 240 # Write input parameters into logger
241 241 if self.logTerse():
242 242 self.log_parameters()
243 243 self.log("\n")
244   -
  244 +
245 245 # Write observation into logger
246 246 if self.logTerse():
247 247 self.log("\n")
... ... @@ -253,19 +253,19 @@ class csresmap(ctools.cscript):
253 253 if self.m_use_maps:
254 254 countmap = gammalib.GSkymap(self["inobs"].filename())
255 255 modelmap = gammalib.GSkymap(self.m_modcube)
256   -
  256 +
257 257 else:
258   -
  258 +
259 259 if self.m_skip_binning:
260 260 cta_counts_cube = gammalib.GCTAEventCube(self.obs[0].events().clone())
261 261  
262 262 else:
263   -
  263 +
264 264 # Write header
265 265 if self.logTerse():
266 266 self.log("\n")
267 267 self.log.header1("Generate binned map (ctbin)")
268   -
  268 +
269 269 # Create countsmap
270 270 bin = ctools.ctbin(self.obs)
271 271 bin["nxpix"].integer(self.m_nxpix)
... ... @@ -283,18 +283,18 @@ class csresmap(ctools.cscript):
283 283 bin["clobber"].boolean(self.m_clobber)
284 284 bin["debug"].boolean(self.m_debug)
285 285 bin.run()
286   -
  286 +
287 287 # Retrieve counts cube
288 288 cta_counts_cube = bin.cube()
289   -
  289 +
290 290 # Assign GCTAEventCube to skymap
291 291 countmap = cta_counts_cube.map()
292   -
  292 +
293 293 # Write header
294 294 if self.logTerse():
295 295 self.log("\n")
296 296 self.log.header1("Generate model map (ctmodel)")
297   -
  297 +
298 298 # Create model map
299 299 model = ctools.ctmodel(self.obs)
300 300 model.cube(cta_counts_cube)
... ... @@ -302,11 +302,11 @@ class csresmap(ctools.cscript):
302 302 model["clobber"].boolean(self.m_clobber)
303 303 model["debug"].boolean(self.m_debug)
304 304 model.run()
305   -
  305 +
306 306 # Get model map into GSkymap object
307 307 modelmap = model.cube().map().copy()
308   -
309   -
  308 +
  309 +
310 310 # Store counts map as residual map. Note that we need a
311 311 # special construct here to avoid memory leaks. This seems
312 312 # to be a SWIG feature as SWIG creates a new object when
... ... @@ -315,15 +315,15 @@ class csresmap(ctools.cscript):
315 315 self.resmap = countmap.copy()
316 316 self.resmap.stack_maps()
317 317 modelmap.stack_maps()
318   -
  318 +
319 319 # Continue calculations depending on given algorithm
320 320 if self.m_algorithm == "SUB":
321   -
  321 +
322 322 # Subtract maps
323 323 self.resmap -= modelmap
324   -
  324 +
325 325 elif self.m_algorithm == "SUBDIV":
326   -
  326 +
327 327 # Subtract and divide by model map
328 328 self.resmap -= modelmap
329 329 self.resmap /= modelmap
... ... @@ -331,7 +331,7 @@ class csresmap(ctools.cscript):
331 331 # if pixel != 0.0:
332 332 # pixel = 1.0/pixel
333 333 #self.resmap *= modelmap
334   -
  334 +
335 335 elif self.m_algorithm == "SUBDIVSQRT":
336 336  
337 337 # subtract and divide by sqrt of model map
... ... @@ -341,15 +341,15 @@ class csresmap(ctools.cscript):
341 341 # if pixel != 0.0:
342 342 # pixel = 1.0/math.sqrt(pixel)
343 343 #self.resmap *= modelmap
344   -
  344 +
345 345 else:
346   -
  346 +
347 347 # Raise error if algorithm is unkown
348 348 raise TypeError("Algorithm \""+self.m_algorithm+"\" not known")
349   -
  349 +
350 350 # Return
351 351 return
352   -
  352 +
353 353  
354 354 # ======================== #
355 355 # Main routine entry point #
... ... @@ -360,10 +360,9 @@ if __name__ == &#39;__main__&#39;:
360 360 """
361 361 # Create instance of application
362 362 app = csresmap(sys.argv)
363   -
  363 +
364 364 # Open logfile
365 365 app.logFileOpen()
366   -
  366 +
367 367 # Execute application
368 368 app.execute()
369   -
370 369 \ No newline at end of file
... ...
cscripts/cssens.py
  View file @5c53c14
... ... @@ -50,7 +50,7 @@ class cssens(ctools.cscript):
50 50 self.obs = None
51 51 self.m_ra = None
52 52 self.m_dec = None
53   -
  53 +
54 54 # Make sure that parfile exists
55 55 file = self.parfile()
56 56  
... ... @@ -69,7 +69,7 @@ class cssens(ctools.cscript):
69 69  
70 70 # Return
71 71 return
72   -
  72 +
73 73 def __del__(self):
74 74 """
75 75 Destructor.
... ... @@ -77,7 +77,7 @@ class cssens(ctools.cscript):
77 77 # Write separator into logger
78 78 if self.logTerse():
79 79 self.log("\n")
80   -
  80 +
81 81 # Return
82 82 return
83 83  
... ... @@ -88,13 +88,13 @@ class cssens(ctools.cscript):
88 88 """
89 89 # Set parfile name
90 90 parfile = self.name+".par"
91   -
  91 +
92 92 try:
93 93 pars = gammalib.GApplicationPars(parfile)
94 94 except:
95 95 # Signal if parfile was not found
96 96 print("Parfile "+parfile+" not found. Create default parfile.")
97   -
  97 +
98 98 # Create default parfile
99 99 pars = gammalib.GApplicationPars()
100 100 pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file"))
... ... @@ -120,10 +120,10 @@ class cssens(ctools.cscript):
120 120 pars.append_standard()
121 121 pars.append(gammalib.GApplicationPar("logfile","f","h","cssens.log","","","Log filename"))
122 122 pars.save(parfile)
123   -
  123 +
124 124 # Return
125 125 return
126   -
  126 +
127 127 def get_parameters(self):
128 128 """
129 129 Get user parameters from parfile.
... ... @@ -160,7 +160,7 @@ class cssens(ctools.cscript):
160 160 self.m_max_iter = self["max_iter"].integer()
161 161 self.m_num_avg = self["num_avg"].integer()
162 162 self.m_type = self["type"].string()
163   -
  163 +
164 164 # Set some fixed parameters
165 165 self.m_log = False # Logging in client tools
166 166 self.m_debug = self["debug"].boolean() # Debugging in client tools
... ... @@ -169,17 +169,17 @@ class cssens(ctools.cscript):
169 169 self.m_ebounds = gammalib.GEbounds(self.m_bins,
170 170 gammalib.GEnergy(self.m_emin, "TeV"),
171 171 gammalib.GEnergy(self.m_emax, "TeV"))
172   -
  172 +
173 173 # Return
174 174 return
175   -
  175 +
176 176 def execute(self):
177 177 """
178 178 Execute the script.
179 179 """
180 180 # Run the script
181 181 self.run()
182   -
  182 +
183 183 # Return
184 184 return
185 185  
... ... @@ -193,12 +193,12 @@ class cssens(ctools.cscript):
193 193  
194 194 # Get parameters
195 195 self.get_parameters()
196   -
  196 +
197 197 # Write input parameters into logger
198 198 if self.logTerse():
199 199 self.log_parameters()
200 200 self.log("\n")
201   -
  201 +
202 202 # Initialise script
203 203 colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux',
204 204 'energy_flux', 'sensitivity']
... ... @@ -228,7 +228,7 @@ class cssens(ctools.cscript):
228 228  
229 229 # Loop over energy bins
230 230 for ieng in range(self.m_ebounds.size()):
231   -
  231 +
232 232 # Set energies
233 233 if self.m_type == "Differential":
234 234 emin = self.m_ebounds.emin(ieng)
... ... @@ -240,11 +240,11 @@ class cssens(ctools.cscript):
240 240 msg = "Invalid sensitivity type \""+self.m_type+"\" encountered."+ \
241 241 " Either use \"Differential\" or \"Integral\"."
242 242 raise gammalib.GException.invalid_value("cssens", msg)
243   -
  243 +
244 244 # Determine sensitivity
245 245 result = self.get_sensitivity(self.obs, emin, emax,
246 246 bkg_model, full_model)
247   -
  247 +
248 248 # Write results
249 249 if ieng == 0:
250 250 f = open(self.m_outfile, 'w')
... ... @@ -257,17 +257,17 @@ class cssens(ctools.cscript):
257 257 writer = csv.DictWriter(f, colnames)
258 258 writer.writerow(result)
259 259 f.close()
260   -
  260 +
261 261 # Append results
262 262 results.append(result)
263   -
  263 +
264 264 # Return
265 265 return
266 266  
267 267 def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
268 268 """
269 269 Returns an observation container.
270   -
  270 +
271 271 Keywords:
272 272 lpnt - Galactic longitude of pointing [deg] (default: 0.0)
273 273 bpnt - Galactic latitude of pointing [deg] (default: 0.0)
... ... @@ -291,16 +291,16 @@ class cssens(ctools.cscript):
291 291  
292 292 # Allocate observation container
293 293 obs = gammalib.GObservations()
294   -
  294 +
295 295 # Set single pointing
296 296 pntdir = gammalib.GSkyDir()
297 297 pntdir.lb_deg(lpnt, bpnt)
298   -
  298 +
299 299 # Create CTA observation
300 300 run = obsutils.set_obs(pntdir, caldb=caldb, irf=irf,
301 301 duration=duration, deadc=deadc,
302 302 emin=emin, emax=emax, rad=rad)
303   -
  303 +
304 304 # Append observation to container
305 305 obs.append(run)
306 306  
... ... @@ -309,25 +309,25 @@ class cssens(ctools.cscript):
309 309 pntdir.lb_deg(lpnt, bpnt+offset)
310 310 self.m_ra = pntdir.ra_deg()
311 311 self.m_dec = pntdir.dec_deg()
312   -
  312 +
313 313 # Return observation container
314 314 return obs
315 315  
316 316 def set_obs_ebounds(self, emin, emax):
317 317 """
318 318 Set energy boundaries for observation in container.
319   -
  319 +
320 320 Parameters:
321 321 emin - Minimum energy
322 322 emax - Maximum energy
323 323 """
324 324 # Loop over all observations in container
325 325 for obs in self.obs:
326   -
  326 +
327 327 # Set energy boundaries
328 328 ebounds = gammalib.GEbounds(emin, emax)
329 329 obs.events().ebounds(ebounds)
330   -
  330 +
331 331 # Return
332 332 return
333 333  
... ... @@ -340,7 +340,7 @@ class cssens(ctools.cscript):
340 340  
341 341 # Get source model
342 342 model = full_model[self.m_srcname]
343   -
  343 +
344 344 # Check that model has a Prefactor
345 345 if not model.has_par("Prefactor"):
346 346 msg = "Model \""+self.m_srcname+"\" has no parameter \"Prefactor\"."+ \
... ... @@ -399,7 +399,7 @@ class cssens(ctools.cscript):
399 399 def get_sensitivity(self, obs, emin, emax, bkg_model, full_model):
400 400 """
401 401 Determine sensitivity for given observations.
402   -
  402 +
403 403 Parameters:
404 404 obs - Observation container
405 405 emin - Minimum energy for fitting and flux computation
... ... @@ -412,7 +412,7 @@ class cssens(ctools.cscript):
412 412  
413 413 # Set energy boundaries
414 414 self.set_obs_ebounds(emin, emax)
415   -
  415 +
416 416 # Determine energy boundaries from first observation in the container
417 417 loge = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
418 418 e_mean = math.pow(10.0, loge)
... ... @@ -517,7 +517,7 @@ class cssens(ctools.cscript):
517 517 self.log("\n")
518 518 for par in model:
519 519 self.log(str(par)+"\n")
520   -
  520 +
521 521 # Fit background and test source
522 522 sim.models(src_model)
523 523 like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
... ... @@ -580,7 +580,7 @@ class cssens(ctools.cscript):
580 580 correct = 1.0
581 581 if ts > 0:
582 582 correct = math.sqrt(self.m_ts_thres/ts)
583   -
  583 +
584 584 # Compute extrapolated fluxes
585 585 crab_flux = correct * crab_flux
586 586 photon_flux = correct * photon_flux
... ... @@ -590,7 +590,7 @@ class cssens(ctools.cscript):
590 590 photon_flux_value.append(photon_flux)
591 591 energy_flux_value.append(energy_flux)
592 592 sensitivity_value.append(sensitivity)
593   -
  593 +
594 594 # Write background and test source fit results
595 595 if self.logExplicit():
596 596 self.log.parformat("Photon flux")
... ... @@ -627,7 +627,7 @@ class cssens(ctools.cscript):
627 627 self.log(" mCrab = ")
628 628 self.log(sensitivity)
629 629 self.log(" erg/cm2/s\n")
630   -
  630 +
631 631 # Compute sliding average of extrapolated fitted prefactor,
632 632 # photon and energy flux. This damps out fluctuations and
633 633 # improves convergence
... ... @@ -648,12 +648,12 @@ class cssens(ctools.cscript):
648 648 photon_flux /= num
649 649 energy_flux /= num
650 650 sensitivity /= num
651   -
  651 +
652 652 # Compare average flux to last average
653 653 if iter > self.m_num_avg:
654 654 if test_crab_flux > 0:
655 655 ratio = crab_flux/test_crab_flux
656   -
  656 +
657 657 # We have 2 convergence criteria:
658 658 # 1. The average flux does not change
659 659 # 2. The flux correction factor is small
... ... @@ -666,10 +666,10 @@ class cssens(ctools.cscript):
666 666 if self.logTerse():
667 667 self.log(" Flux is zero.\n")
668 668 break
669   -
  669 +
670 670 # Use average for next iteration
671 671 test_crab_flux = crab_flux
672   -
  672 +
673 673 # Exit loop if number of trials exhausted
674 674 if (iter >= self.m_max_iter):
675 675 if self.logTerse():
... ... @@ -723,13 +723,13 @@ class cssens(ctools.cscript):
723 723 self.log("\n")
724 724 for par in model:
725 725 self.log(str(par)+"\n")
726   -
  726 +
727 727 # Store result
728 728 result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
729 729 'crab_flux': crab_flux, 'photon_flux': photon_flux, \
730 730 'energy_flux': energy_flux, \
731 731 'sensitivity': sensitivity}
732   -
  732 +
733 733 # Return result
734 734 return result
735 735  
... ... @@ -756,7 +756,6 @@ if __name__ == &#39;__main__&#39;:
756 756 """
757 757 # Create instance of application
758 758 app = cssens(sys.argv)
759   -
  759 +
760 760 # Run application
761 761 app.execute()
762   -
763 762 \ No newline at end of file
... ...
cscripts/csspec.py
  View file @5c53c14
... ... @@ -37,14 +37,14 @@ class csspec(ctools.cscript):
37 37 """
38 38 Constructor.
39 39 """
40   -
  40 +
41 41 # Set name
42 42 self.name = "csspec"
43 43 self.version = "1.0.0"
44   -
  44 +
45 45 # Initialise some members
46 46 self.obs = None
47   -
  47 +
48 48 # Initialise some members
49 49 if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations):
50 50 self.obs = argv[0]
... ... @@ -53,7 +53,7 @@ class csspec(ctools.cscript):
53 53 self.obs = gammalib.GObservations()
54 54 self.obs.clear()
55 55 self.m_outfile = ""
56   -
  56 +
57 57 # Make sure that parfile exists
58 58 file = self.parfile()
59 59  
... ... @@ -71,7 +71,7 @@ class csspec(ctools.cscript):
71 71  
72 72 # Return
73 73 return
74   -
  74 +
75 75 def __del__(self):
76 76 """
77 77 Destructor.
... ... @@ -79,7 +79,7 @@ class csspec(ctools.cscript):
79 79 # Write separator into logger
80 80 if self.logTerse():
81 81 self.log("\n")
82   -
  82 +
83 83 # Return
84 84 return
85 85  
... ... @@ -90,13 +90,13 @@ class csspec(ctools.cscript):
90 90 """
91 91 # Set parfile name
92 92 parfile = self.name+".par"
93   -
  93 +
94 94 try:
95 95 pars = gammalib.GApplicationPars(parfile)
96 96 except:
97 97 # Signal if parfile was not found
98 98 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
99   -
  99 +
100 100 # Create default parfile
101 101 pars = gammalib.GApplicationPars()
102 102 pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file"))
... ... @@ -128,10 +128,10 @@ class csspec(ctools.cscript):
128 128 pars.append_standard()
129 129 pars.append(gammalib.GApplicationPar("logfile","f","h","csspec.log","","","Log filename"))
130 130 pars.save(parfile)
131   -
  131 +
132 132 # Return
133 133 return
134   -
  134 +
135 135 def get_parameters(self):
136 136 """
137 137 Get parameters from parfile and setup the observation.
... ... @@ -144,10 +144,10 @@ class csspec(ctools.cscript):
144 144 # Set models if we have none
145 145 if self.obs.models().size() == 0:
146 146 self.obs.models(self["inmodel"].filename())
147   -
  147 +
148 148 # Get source name
149 149 self.m_srcname = self["srcname"].string()
150   -
  150 +
151 151 # Get ebounds
152 152 self.m_ebounds = self.create_ebounds()
153 153  
... ... @@ -165,32 +165,32 @@ class csspec(ctools.cscript):
165 165  
166 166 # Read other parameters
167 167 self.m_outfile = self["outfile"].filename()
168   -
  168 +
169 169 # Get other parameeters
170 170 self.m_calc_ulimit = self["calc_ulim"].boolean()
171 171 self.m_calc_ts = self["calc_ts"].boolean()
172 172 self.m_fix_bkg = self["fix_bkg"].boolean()
173 173 self.m_fix_srcs = self["fix_srcs"].boolean()
174   -
  174 +
175 175 # Set some fixed parameters
176 176 self.m_log = False # Logging in client tools
177 177 self.m_chatter = self["chatter"].integer()
178 178 self.m_clobber = self["clobber"].boolean()
179 179 self.m_debug = self["debug"].boolean()
180   -
  180 +
181 181 # Return
182 182 return
183   -
  183 +
184 184 def models(self, models):
185 185 """
186 186 Set model.
187 187 """
188 188 # Copy models
189 189 self.obs.models(models.clone())
190   -
  190 +
191 191 # Return
192 192 return
193   -
  193 +
194 194 def execute(self):
195 195 """
196 196 Execute the script.
... ... @@ -200,7 +200,7 @@ class csspec(ctools.cscript):
200 200  
201 201 # Save residual map
202 202 self.fits.saveto(self.m_outfile, self.m_clobber)
203   -
  203 +
204 204 # Return
205 205 return
206 206  
... ... @@ -214,12 +214,12 @@ class csspec(ctools.cscript):
214 214  
215 215 # Get parameters
216 216 self.get_parameters()
217   -
  217 +
218 218 # Write input parameters into logger
219 219 if self.logTerse():
220 220 self.log_parameters()
221 221 self.log("\n")
222   -
  222 +
223 223 # Write observation into logger
224 224 if self.logTerse():
225 225 self.log("\n")
... ... @@ -234,15 +234,15 @@ class csspec(ctools.cscript):
234 234  
235 235 # Adjust model parameters dependent on input user parameters
236 236 for model in self.obs.models():
237   -
  237 +
238 238 # Set TS flag for all models to false.
239 239 # Source of interest will be set to true later
240 240 model.tscalc(False)
241   -
  241 +
242 242 # Log model name
243 243 if self.logExplicit():
244 244 self.log.header3(model.name())
245   -
  245 +
246 246 # Deal with the source of interest
247 247 if model.name() == self.m_srcname:
248 248 for par in model:
... ... @@ -254,32 +254,32 @@ class csspec(ctools.cscript):
254 254 model.spectral()[0].free()
255 255 if self.m_calc_ts:
256 256 model.tscalc(True)
257   -
  257 +
258 258 elif self.m_fix_bkg and not model.classname() == "GModelSky":
259 259 for par in model:
260 260 if par.is_free() and self.logExplicit():
261 261 self.log(" Fixing \""+par.name()+"\"\n")
262 262 par.fix()
263   -
  263 +
264 264 elif self.m_fix_srcs and model.classname() == "GModelSky":
265 265 for par in model:
266 266 if par.is_free() and self.logExplicit():
267 267 self.log(" Fixing \""+par.name()+"\"\n")
268 268 par.fix()
269   -
  269 +
270 270 # Write header
271 271 if self.logTerse():
272 272 self.log("\n")
273 273 self.log.header1("Generate spectrum")
274   -
  274 +
275 275 # Initialise FITS Table with extension "SPECTRUM"
276 276 table = gammalib.GFitsBinTable(self.m_ebounds.size())
277 277 table.extname("SPECTRUM")
278   -
  278 +
279 279 # Add Header for compatibility with gammalib.GMWLSpectrum
280 280 table.card("INSTRUME", "CTA", "Name of Instrument")
281 281 table.card("TELESCOP", "CTA", "Name of Telescope")
282   -
  282 +
283 283 # Create FITS table columns
284 284 energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size())
285 285 energy.unit("TeV")
... ... @@ -295,7 +295,7 @@ class csspec(ctools.cscript):
295 295 ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size())
296 296 ulim_values.unit("erg/cm2/s")
297 297 Npred_values = gammalib.GFitsTableDoubleCol("Npred", self.m_ebounds.size())
298   -
  298 +
299 299  
300 300 # Loop over energy bins
301 301 for i in range(self.m_ebounds.size()):
... ... @@ -310,18 +310,18 @@ class csspec(ctools.cscript):
310 310 emax = self.m_ebounds.emax(i)
311 311 elogmean = self.m_ebounds.elogmean(i)
312 312 elogmean2 = elogmean.MeV() * elogmean.MeV()
313   -
  313 +
314 314 # Store energy as TeV
315 315 energy[i] = elogmean.TeV()
316   -
  316 +
317 317 # Store energy errors
318 318 energy_low[i] = (elogmean - emin).TeV()
319 319 energy_high[i] = (emax - elogmean).TeV()
320   -
  320 +
321 321 # Log information
322 322 if self.logExplicit():
323   - self.log.header3("Selecting events")
324   -
  323 + self.log.header3("Selecting events")
  324 +
325 325 # Select events
326 326 select = ctools.ctselect(self.obs)
327 327 select["emin"] = emin.TeV()
... ... @@ -342,7 +342,7 @@ class csspec(ctools.cscript):
342 342 # Header
343 343 if self.logTerse():
344 344 self.log.header3("Binning events")
345   -
  345 +
346 346 # Bin events
347 347 bin = ctools.ctbin(select.obs())
348 348 bin["usepnt"] = False
... ... @@ -358,11 +358,11 @@ class csspec(ctools.cscript):
358 358 bin["coordsys"] = self.m_coordsys
359 359 bin["proj"] = self.m_proj
360 360 bin.run()
361   -
  361 +
362 362 # Header
363 363 if self.logTerse():
364 364 self.log.header3("Creating exposure cube")
365   -
  365 +
366 366 # Create exposure cube
367 367 expcube = ctools.ctexpcube(select.obs())
368 368 expcube["incube"] = "NONE"
... ... @@ -379,11 +379,11 @@ class csspec(ctools.cscript):
379 379 expcube["coordsys"] = self.m_coordsys
380 380 expcube["proj"] = self.m_proj
381 381 expcube.run()
382   -
  382 +
383 383 # Header
384 384 if self.logTerse():
385 385 self.log.header3("Creating PSF cube")
386   -
  386 +
387 387 # Create psf cube
388 388 psfcube = ctools.ctpsfcube(select.obs())
389 389 psfcube["incube"] = "NONE"
... ... @@ -400,11 +400,11 @@ class csspec(ctools.cscript):
400 400 psfcube["coordsys"] = self.m_coordsys
401 401 psfcube["proj"] = self.m_proj
402 402 psfcube.run()
403   -
  403 +
404 404 # Header
405 405 if self.logTerse():
406 406 self.log.header3("Creating background cube")
407   -
  407 +
408 408 # Create background cube
409 409 bkgcube = ctools.ctbkgcube(select.obs())
410 410 bkgcube["incube"] = "NONE"
... ... @@ -421,37 +421,37 @@ class csspec(ctools.cscript):
421 421 bkgcube["coordsys"] = self.m_coordsys
422 422 bkgcube["proj"] = self.m_proj
423 423 bkgcube.run()
424   -
  424 +
425 425 # Set new binned observation
426 426 obs = bin.obs()
427   -
  427 +
428 428 # Set precomputed binned response
429 429 obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube())
430 430  
431 431 # Get new models
432 432 models = bkgcube.models()
433   -
  433 +
434 434 # Fix background models if required
435 435 if self.m_fix_bkg:
436 436 for model in models:
437 437 if not model.classname() == "GModelSky":
438 438 for par in model:
439 439 par.fix()
440   -
  440 +
441 441 # Set new models to binned observation
442 442 obs.models(models)
443   -
  443 +
444 444 # Header
445 445 if self.logTerse():
446 446 self.log.header3("Performing fit")
447   -
  447 +
448 448 # Likelihood
449 449 like = ctools.ctlike(obs)
450 450 like.run()
451   -
  451 +
452 452 # Skip bin if no event was present
453 453 if like.obs().logL() == 0.0:
454   -
  454 +
455 455 # Log information
456 456 if self.logTerse():
457 457 self.log("No event in this bin. Bin is skipped\n")
... ... @@ -463,7 +463,7 @@ class csspec(ctools.cscript):
463 463 ulim_values[i] = 0.0
464 464 Npred_values[i] = 0.0
465 465 continue
466   -
  466 +
467 467 # Get results
468 468 fitted_models = like.obs().models()
469 469 source = fitted_models[self.m_srcname]
... ... @@ -471,16 +471,16 @@ class csspec(ctools.cscript):
471 471 # Calculate Upper Limit
472 472 ulimit_value = -1.0
473 473 if self.m_calc_ulimit:
474   -
  474 +
475 475 # Logging information
476 476 if self.logTerse():
477 477 self.log.header3("Computing upper limit")
478   -
  478 +
479 479 # Create upper limit object
480 480 ulimit = ctools.ctulimit(like.obs())
481 481 ulimit["srcname"] = self.m_srcname
482 482 ulimit["eref"] = elogmean.TeV()
483   -
  483 +
484 484 # Try to run upper limit and catch exceptions
485 485 try:
486 486 ulimit.run()
... ... @@ -489,37 +489,37 @@ class csspec(ctools.cscript):
489 489 if self.logTerse():
490 490 self.log("Upper limit calculation failed\n")
491 491 ulimit_value = -1.0
492   -
  492 +
493 493 # Get TS value
494 494 TS = -1.0
495 495 if self.m_calc_ts:
496 496 TS = source.ts()
497   -
  497 +
498 498 # Compute Npred value
499 499 Npred = 0.0
500 500 for observation in like.obs():
501 501 Npred += observation.npred(source)
502   -
  502 +
503 503 # Get differential flux
504 504 fitted_flux = source.spectral().eval(elogmean,gammalib.GTime())
505   -
  505 +
506 506 # Compute flux error
507 507 parvalue = source.spectral()[0].value()
508 508 rel_error = source.spectral()[0].error()/parvalue
509 509 e_flux = fitted_flux*rel_error
510   -
  510 +
511 511 # Set values for storage
512 512 TSvalues[i] = TS
513   -
  513 +
514 514 # Set npred values
515 515 Npred_values[i] = Npred
516   -
  516 +
517 517 # Convert fluxes to nuFnu
518 518 flux[i] = fitted_flux * elogmean2 * gammalib.MeV2erg
519 519 flux_err[i] = e_flux * elogmean2 * gammalib.MeV2erg
520 520 if ulimit_value > 0.0:
521 521 ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg
522   -
  522 +
523 523 # Log information
524 524 if self.logExplicit():
525 525 self.log("Bin "+str(i)+" ["+str(emin.TeV())+"-"+str(emax.TeV())+"] TeV: ")
... ... @@ -540,14 +540,14 @@ class csspec(ctools.cscript):
540 540 table.append(TSvalues)
541 541 table.append(ulim_values)
542 542 table.append(Npred_values)
543   -
  543 +
544 544 # Create the FITS file now
545 545 self.fits = gammalib.GFits()
546 546 self.fits.append(table)
547   -
  547 +
548 548 # Return
549 549 return
550   -
  550 +
551 551  
552 552 # ======================== #
553 553 # Main routine entry point #
... ... @@ -558,10 +558,9 @@ if __name__ == &#39;__main__&#39;:
558 558 """
559 559 # Create instance of application
560 560 app = csspec(sys.argv)
561   -
  561 +
562 562 # Open logfile
563 563 app.logFileOpen()
564   -
  564 +
565 565 # Execute application
566 566 app.execute()
567   -
568 567 \ No newline at end of file
... ...
cscripts/cstsdist.py
  View file @5c53c14
... ... @@ -43,12 +43,12 @@ class cstsdist(ctools.cscript):
43 43 # Set name
44 44 self.name = "cstsdist"
45 45 self.version = "1.0.0"
46   -
  46 +
47 47 # Initialise some members
48 48 self.obs = None
49 49 self.bkg_model = None
50 50 self.full_model = None
51   -
  51 +
52 52 # Make sure that parfile exists
53 53 file = self.parfile()
54 54  
... ... @@ -66,7 +66,7 @@ class cstsdist(ctools.cscript):
66 66  
67 67 # Return
68 68 return
69   -
  69 +
70 70 def __del__(self):
71 71 """
72 72 Destructor.
... ... @@ -74,7 +74,7 @@ class cstsdist(ctools.cscript):
74 74 # Write separator into logger
75 75 if self.logTerse():
76 76 self.log("\n")
77   -
  77 +
78 78 # Return
79 79 return
80 80  
... ... @@ -85,13 +85,13 @@ class cstsdist(ctools.cscript):
85 85 """
86 86 # Set parfile name
87 87 parfile = self.name+".par"
88   -
  88 +
89 89 try:
90 90 pars = gammalib.GApplicationPars(parfile)
91 91 except:
92 92 # Signal if parfile was not found
93 93 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
94   -
  94 +
95 95 # Create default parfile
96 96 pars = gammalib.GApplicationPars()
97 97 pars.append(gammalib.GApplicationPar("inobs","f","h","NONE","","","Event list, counts cube, or observation definition file"))
... ... @@ -121,15 +121,15 @@ class cstsdist(ctools.cscript):
121 121 pars.append_standard()
122 122 pars.append(gammalib.GApplicationPar("logfile","f","h","cstsdist.log","","","Log filename"))
123 123 pars.save(parfile)
124   -
  124 +
125 125 # Return
126 126 return
127   -
  127 +
128 128 def get_parameters(self):
129 129 """
130 130 Get parameters from parfile and setup the observation.
131 131 """
132   -
  132 +
133 133 # Set observation if not done before
134 134 if self.obs == None or self.obs.size() == 0:
135 135 self.obs = self.get_observations()
... ... @@ -142,10 +142,10 @@ class cstsdist(ctools.cscript):
142 142  
143 143 # Get source name
144 144 self.m_srcname = self["srcname"].string()
145   -
  145 +
146 146 # Get number of energy bins
147 147 self.m_enumbins = self["enumbins"].integer()
148   -
  148 +
149 149 # Read parameters for binned if requested
150 150 if not self.m_enumbins == 0:
151 151 self.m_npix = self["npix"].integer()
... ... @@ -158,7 +158,7 @@ class cstsdist(ctools.cscript):
158 158 # Set models if we have none
159 159 if self.obs.models().size() == 0:
160 160 self.obs.models(self["inmodel"].filename())
161   -
  161 +
162 162 # Get other parameters
163 163 self.m_outfile = self["outfile"].filename()
164 164 self.m_ntrials = self["ntrials"].integer()
... ... @@ -171,24 +171,24 @@ class cstsdist(ctools.cscript):
171 171  
172 172 # Return
173 173 return
174   -
  174 +
175 175 def models(self, models):
176 176 """
177 177 Set model.
178 178 """
179 179 # Copy models
180 180 self.obs.models(models.copy())
181   -
  181 +
182 182 # Return
183 183 return
184   -
  184 +
185 185 def execute(self):
186 186 """
187 187 Execute the script.
188 188 """
189 189 # Run the script
190 190 self.run()
191   -
  191 +
192 192 # Return
193 193 return
194 194  
... ... @@ -202,7 +202,7 @@ class cstsdist(ctools.cscript):
202 202  
203 203 # Get parameters
204 204 self.get_parameters()
205   -
  205 +
206 206 # Write input parameters into logger
207 207 if self.logTerse():
208 208 self.log_parameters()
... ... @@ -221,7 +221,7 @@ class cstsdist(ctools.cscript):
221 221 self.log.header2("Full model")
222 222 self.log(str(full_model))
223 223 self.log("\n")
224   -
  224 +
225 225 # Write observation into logger
226 226 if self.logTerse():
227 227 self.log("\n")
... ... @@ -243,10 +243,10 @@ class cstsdist(ctools.cscript):
243 243  
244 244 # Loop over trials
245 245 for seed in range(self.m_ntrials):
246   -
  246 +
247 247 # Make a trial
248 248 result = self.trial(seed, full_model, bkg_model)
249   -
  249 +
250 250 # Write out result immediately
251 251 if seed == 0:
252 252 file = open(self.m_outfile, 'w')
... ... @@ -257,10 +257,10 @@ class cstsdist(ctools.cscript):
257 257 writer = csv.DictWriter(file, result['colnames'])
258 258 writer.writerow(result['values'])
259 259 file.close()
260   -
  260 +
261 261 # Return
262 262 return
263   -
  263 +
264 264 def set_models(self, fitpos=False, fitspec=False):
265 265 """
266 266 Set full and background model.
... ... @@ -270,7 +270,7 @@ class cstsdist(ctools.cscript):
270 270  
271 271 # Get source model
272 272 model = full_model[self.m_srcname]
273   -
  273 +
274 274 # Check that model has a Prefactor
275 275 if not model.has_par("Prefactor"):
276 276 msg = "Model \""+self.m_srcname+"\" has no parameter \"Prefactor\"."+ \
... ... @@ -323,16 +323,16 @@ class cstsdist(ctools.cscript):
323 323 def set_obs(self):
324 324 """
325 325 Returns an observation container with a set of CTA observations.
326   -
  326 +
327 327 Keywords:
328 328 """
329   -
  329 +
330 330 # Setup observation definition list
331 331 obsdeflist = obsutils.set_obs_patterns(self.m_pattern,
332 332 ra=self["ra"].real(),
333 333 dec=self["dec"].real(),
334 334 offset=self["offset"].real())
335   -
  335 +
336 336 # Create list of observations
337 337 obs = obsutils.set_obs_list(obsdeflist,
338 338 tstart=self["tmin"].real(),
... ... @@ -342,14 +342,14 @@ class cstsdist(ctools.cscript):
342 342 emax=self["emax"].real(),
343 343 rad=self["rad"].real(),
344 344 irf=self["irf"].string(), caldb=self["caldb"].string())
345   -
  345 +
346 346 # Return observation container
347 347 return obs
348 348  
349 349 def trial(self, seed, full_model, bkg_model):
350 350 """
351 351 Create the TS for a single trial.
352   -
  352 +
353 353 Parameters:
354 354 seed - Random number generator seed
355 355 """
... ... @@ -364,7 +364,7 @@ class cstsdist(ctools.cscript):
364 364 binsz=self.m_binsz,
365 365 npix=self.m_npix,
366 366 log=self.m_log, debug=self.m_debug)
367   -
  367 +
368 368 # Determine number of events in simulation
369 369 nevents = 0.0
370 370 for run in sim:
... ... @@ -376,7 +376,7 @@ class cstsdist(ctools.cscript):
376 376 self.log.parformat("Number of simulated events")
377 377 self.log(nevents)
378 378 self.log("\n")
379   -
  379 +
380 380 # Fit background only
381 381 sim.models(bkg_model)
382 382 like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
... ... @@ -407,7 +407,7 @@ class cstsdist(ctools.cscript):
407 407 LogL_all = like_all.opt().value()
408 408 npred_all = like_all.obs().npred()
409 409 ts = 2.0*(LogL_bgm-LogL_all)
410   -
  410 +
411 411 # Write background and test source fit results
412 412 if self.logExplicit():
413 413 self.log.header3("Background and test source model fit")
... ... @@ -426,7 +426,7 @@ class cstsdist(ctools.cscript):
426 426 self.log("\n")
427 427 for par in model:
428 428 self.log(str(par)+"\n")
429   -
  429 +
430 430 # Write result
431 431 elif self.logTerse():
432 432 self.log.parformat("Trial "+str(seed))
... ... @@ -437,11 +437,11 @@ class cstsdist(ctools.cscript):
437 437 self.log("+/-")
438 438 self.log(result_all[self.m_srcname]["Prefactor"].error())
439 439 self.log("\n")
440   -
  440 +
441 441 # Initialise results
442 442 colnames = []
443 443 values = {}
444   -
  444 +
445 445 # Set TS value
446 446 colnames.append("TS")
447 447 values["TS"] = ts
... ... @@ -465,7 +465,7 @@ class cstsdist(ctools.cscript):
465 465 # Set Npred for full fit
466 466 colnames.append("Npred_all")
467 467 values["Npred_all"] = npred_all
468   -
  468 +
469 469 # Gather free full fit parameters
470 470 for i in range(result_all.size()):
471 471 model = result_all[i]
... ... @@ -473,22 +473,22 @@ class cstsdist(ctools.cscript):
473 473 for k in range(model.size()):
474 474 par = model[k]
475 475 if par.is_free():
476   -
  476 +
477 477 # Set parameter name
478 478 name = model_name+"_"+par.name()
479   -
  479 +
480 480 # Append value
481 481 colnames.append(name)
482 482 values[name] = par.value()
483   -
  483 +
484 484 # Append error
485 485 name = "Unc_"+name
486 486 colnames.append(name)
487 487 values[name] = par.error()
488   -
  488 +
489 489 # Bundle together results
490 490 result = {'colnames': colnames, 'values': values}
491   -
  491 +
492 492 # Return
493 493 return result
494 494  
... ... @@ -502,10 +502,9 @@ if __name__ == &#39;__main__&#39;:
502 502 """
503 503 # Create instance of application
504 504 app = cstsdist(sys.argv)
505   -
  505 +
506 506 # Open logfile
507 507 app.logFileOpen()
508   -
  508 +
509 509 # Execute application
510 510 app.execute()
511   -
512 511 \ No newline at end of file
... ...
cscripts/cta_kb_bkgrates.py
  View file @5c53c14
... ... @@ -32,7 +32,7 @@ import sys
32 32 def make_file_function(irfname, filename):
33 33 """
34 34 Generate file function from model.
35   -
  35 +
36 36 Parameters:
37 37 irfname - IRF filename
38 38 filename - File function filename
... ... @@ -40,16 +40,16 @@ def make_file_function(irfname, filename):
40 40 # Open performance file and file function
41 41 irf = open(irfname, "r")
42 42 file = open(filename, "w")
43   -
  43 +
44 44 # Read data
45 45 for row in irf:
46   -
  46 +
47 47 # Skip header and trailer
48 48 if row.find("log(E)") != -1:
49 49 continue
50 50 if row.find("--------") != -1:
51 51 break
52   -
  52 +
53 53 # Extract data
54 54 fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",")
55 55 logE = float(fields[0])
... ... @@ -59,16 +59,16 @@ def make_file_function(irfname, filename):
59 59 eres = float(fields[4])
60 60 bgd = float(fields[5])
61 61 sens = float(fields[6])
62   -
  62 +
63 63 # Compute energy (in MeV)
64 64 energy = math.pow(10.0, logE)*1.0e6
65 65 emin = math.pow(10.0, logE-0.1)*1.0e6
66 66 emax = math.pow(10.0, logE+0.1)*1.0e6
67 67 ewidth = emax-emin
68   -
  68 +
69 69 # Compute solid angle of r80 region in sr
70 70 omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80)))
71   -
  71 +
72 72 # Compute background rate per steradian and MeV
73 73 if omega > 0.0:
74 74 bkg_rate = bgd/omega/ewidth
... ... @@ -77,11 +77,11 @@ def make_file_function(irfname, filename):
77 77  
78 78 # Write
79 79 file.write(str(energy)+" "+str(bkg_rate)+"\n")
80   -
  80 +
81 81 # Close files
82 82 irf.close()
83 83 file.close()
84   -
  84 +
85 85 # Return
86 86 return
87 87  
... ... @@ -111,6 +111,6 @@ if __name__ == &#39;__main__&#39;:
111 111 irf = tail.strip(".dat")
112 112 filename = "bkg_"+irf+".txt"
113 113 print(filename)
114   -
  114 +
115 115 # Make background file
116 116 make_file_function(irfname, filename)
... ...
cscripts/cta_kb_root2caldb.py
  View file @5c53c14
... ... @@ -50,7 +50,7 @@ class caldb():
50 50 """
51 51 # Store root path
52 52 self.path = path
53   -
  53 +
54 54 # Initialise some members
55 55 self.cif = None
56 56 self.irf = None
... ... @@ -62,7 +62,7 @@ class caldb():
62 62  
63 63 # Store UTC date of job execution
64 64 self.utc = datetime.utcnow().isoformat()[:19]
65   -
  65 +
66 66 # Set CALDB information
67 67 self.cal_tel = "CTA"
68 68 self.cal_inst = inst.upper()
... ... @@ -99,13 +99,13 @@ class caldb():
99 99 self.bgd_doc = "???"
100 100 self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
101 101 self.bgd_desc = "CTA background"
102   -
  102 +
103 103 # Create directory structure
104 104 self.make_dirs()
105   -
  105 +
106 106 # Return
107 107 return
108   -
  108 +
109 109 def __del__(self):
110 110 """
111 111 Destructor. Close all FITS files.
... ... @@ -115,18 +115,18 @@ class caldb():
115 115  
116 116 # Return
117 117 return
118   -
  118 +
119 119 def make_dirs(self):
120 120 """
121 121 Generate CALDB directory structure for one observation identifier.
122 122 The structure is given by
123   -
  123 +
124 124 data/<tel>/<inst>/bcf/
125   -
  125 +
126 126 where <tel> is "cta" and <inst> is the instrument specified in the
127 127 CALDB constructor (the instrument may be used for different array
128 128 configurations).
129   -
  129 +
130 130 Parameters:
131 131 None
132 132 Keywords:
... ... @@ -139,7 +139,7 @@ class caldb():
139 139  
140 140 # Set directory names
141 141 self.bcf_dir = self.base_dir+"/bcf"
142   -
  142 +
143 143 # Optionally prefix path
144 144 if len(self.path) > 0:
145 145 self.base_path = self.path+"/"+self.base_dir
... ... @@ -151,17 +151,17 @@ class caldb():
151 151 # Create directory structure
152 152 if not os.path.isdir(self.bcf_path):
153 153 os.makedirs(self.bcf_path)
154   -
  154 +
155 155 # Return
156 156 return
157   -
  157 +
158 158 def add(self, rspname, split=False, clobber=True):
159 159 """
160 160 Add new calibration. The actual version will put
161 161 all calibrations in the same file, although each part of the response
162 162 function will have its own logical name. We can thus easily modify
163 163 the script to put each calibration information in a separate file.
164   -
  164 +
165 165 Parameters:
166 166 rspname - Response name
167 167 Keywords:
... ... @@ -179,21 +179,21 @@ class caldb():
179 179 self.psf_file = "irf_"+rspname+".dat"
180 180 self.edisp_file = "irf_"+rspname+".dat"
181 181 self.bgd_file = "irf_"+rspname+".dat"
182   -
  182 +
183 183 # Open calibration database index
184 184 if self.cif == None:
185 185 self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True)
186   -
  186 +
187 187 # If file has no CIF extension than create it now
188 188 try:
189 189 self.hdu_cif = self.cif.table("CIF")
190 190 except:
191 191 self.create_cif()
192 192 self.hdu_cif = self.cif.table("CIF")
193   -
  193 +
194 194 # Set response name
195 195 self.cal_name = "NAME("+rspname+")"
196   -
  196 +
197 197 # Return
198 198 return
199 199  
... ... @@ -210,7 +210,7 @@ class caldb():
210 210 # necessarily have all the information at hand (in particular
211 211 # about the boundaries)
212 212 #self.add_cif_info()
213   -
  213 +
214 214 # Close CIF
215 215 if self.cif != None:
216 216 self.cif.save(True)
... ... @@ -222,7 +222,7 @@ class caldb():
222 222  
223 223 # Return
224 224 return
225   -
  225 +
226 226 def create_cif(self):
227 227 """
228 228 Create Calibration Index File (CIF) extension in CIF FITS file.
... ... @@ -234,9 +234,9 @@ class caldb():
234 234 """
235 235 # Create binary table
236 236 table = gammalib.GFitsBinTable()
237   -
  237 +
238 238 # Set boundary
239   -
  239 +
240 240 # Attach columns. Reference: CAL/GEN/92-008
241 241 table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10))
242 242 table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10))
... ... @@ -256,20 +256,20 @@ class caldb():
256 256 table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0))
257 257 table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8))
258 258 table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70))
259   -
  259 +
260 260 # Set keywords. Reference: CAL/GEN/92-008
261 261 table.extname("CIF")
262 262 table.card("CIFVERSN", "1992a", "Version of CIF format")
263   -
  263 +
264 264 # Attach table to FITS file. Note that at this moment the
265 265 # FITS table gets connected to the FITS file. Yet since
266 266 # nothing was yet written to the FITS file, we cannot read
267 267 # anything from it.
268 268 self.cif.append(table)
269   -
  269 +
270 270 # Return
271 271 return
272   -
  272 +
273 273 def add_cif_info(self):
274 274 """
275 275 Add information to CIF extension.
... ... @@ -287,13 +287,13 @@ class caldb():
287 287  
288 288 # Append 4 rows to CIF extension
289 289 self.hdu_cif.append_rows(4)
290   -
  290 +
291 291 # Add generic information for these 4 rows
292 292 for i in range(4):
293   -
  293 +
294 294 # Set row number
295 295 row = i+self.hdu_cif.nrows()-4
296   -
  296 +
297 297 # Set element
298 298 self.hdu_cif["TELESCOP"][row] = self.cal_tel
299 299 self.hdu_cif["INSTRUME"][row] = self.cal_inst
... ... @@ -307,7 +307,7 @@ class caldb():
307 307 self.hdu_cif["REF_TIME"][row] = self.ref_time
308 308 self.hdu_cif["CAL_QUAL"][row] = self.cal_qual
309 309 self.hdu_cif["CAL_DATE"][row] = self.cal_date
310   -
  310 +
311 311 # Add effective area information
312 312 row = self.hdu_cif.nrows()-4
313 313 self.hdu_cif["CAL_DIR"][row] = self.bcf_dir
... ... @@ -359,14 +359,14 @@ class caldb():
359 359 self.hdu_cif["CAL_CBD"][row, n] = "NONE"
360 360 else:
361 361 self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n]
362   -
  362 +
363 363 # Return
364 364 return
365   -
  365 +
366 366 def root2caldb(self, filename):
367 367 """
368 368 Generate files.
369   -
  369 +
370 370 Parameters:
371 371 filename - Performance file name
372 372 """
... ... @@ -383,16 +383,16 @@ class caldb():
383 383 if self.ea_file == self.psf_file and \
384 384 self.ea_file == self.edisp_file and \
385 385 self.ea_file == self.bgd_file :
386   - split = False
387   - fprf = open(ea_filename, "w")
388   - files = [fprf]
  386 + split = False
  387 + fprf = open(ea_filename, "w")
  388 + files = [fprf]
389 389 else:
390   - split = True
391   - fea = open(ea_filename, "w")
392   - fpsf = open(psf_filename, "w")
393   - fedisp = open(edisp_filename, "w")
394   - fbgd = open(bgd_filename, "w")
395   - files = [fea, fpsf, fedisp, fbgd]
  390 + split = True
  391 + fea = open(ea_filename, "w")
  392 + fpsf = open(psf_filename, "w")
  393 + fedisp = open(edisp_filename, "w")
  394 + fbgd = open(bgd_filename, "w")
  395 + files = [fea, fpsf, fedisp, fbgd]
396 396  
397 397 # Write header
398 398 for file in files:
... ... @@ -420,12 +420,12 @@ class caldb():
420 420 rootfile.GetObject("AngRes", angres)
421 421 rootfile.GetObject("AngRes80", angres80)
422 422 rootfile.GetObject("ERes", eres)
423   -
  423 +
424 424 # Setup vectors
425 425 energies = bgrate.GetXaxis()
426 426 nbins_eng = energies.GetNbins()
427 427 for i in range(nbins_eng):
428   -
  428 +
429 429 # Get logE and background rate per squared degree
430 430 logE = energies.GetBinCenter(i+1)
431 431 bgd = bgrate.GetBinContent(i+1,1)
... ... @@ -436,7 +436,7 @@ class caldb():
436 436 r80 = angres80.GetBinContent(i+1,1)
437 437 dE = eres.GetBinContent(i+1,1)
438 438 diffs = sens.GetBinContent(i+1,1)
439   -
  439 +
440 440 # Skip any NaNs
441 441 if math.isnan(area80):
442 442 continue
... ... @@ -446,10 +446,10 @@ class caldb():
446 446 emin = math.pow(10.0, logE-0.1)*1.0e6
447 447 emax = math.pow(10.0, logE+0.1)*1.0e6
448 448 ewidth = emax-emin
449   -
  449 +
450 450 # Convert into background rate per steradian and MeV
451 451 bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth
452   -
  452 +
453 453 # Compute control background rate (this only works for KB files
454 454 # as they are for 80% containment radius). But we don't use this
455 455 # anyways, it's just for control on display ...
... ... @@ -458,7 +458,7 @@ class caldb():
458 458 bkg_rate2 = bgd / omega / ewidth
459 459 else:
460 460 bkg_rate2 = 0.0
461   -
  461 +
462 462 # Compute full effective area by dividing area80/0.80
463 463 area = area80/0.80
464 464  
... ... @@ -494,7 +494,7 @@ class caldb():
494 494 # necessarily have all the information at hand (in particular
495 495 # about the boundaries)
496 496 self.add_cif_info()
497   -
  497 +
498 498 # Return
499 499 return
500 500  
... ... @@ -521,7 +521,7 @@ if __name__ == &#39;__main__&#39;:
521 521  
522 522 # Loop over all files
523 523 for file in files:
524   -
  524 +
525 525 # Get filename
526 526 head, tail = os.path.split(file)
527 527  
... ... @@ -540,6 +540,6 @@ if __name__ == &#39;__main__&#39;:
540 540  
541 541 # Open calibration file
542 542 irf.add(name)
543   -
  543 +
544 544 # Translate ROOT to CALDB information
545 545 irf.root2caldb(file)
... ...
cscripts/cta_make_nodes.py
  View file @5c53c14
... ... @@ -56,7 +56,7 @@ def write_node(file, energy, scale):
56 56  
57 57 # Return
58 58 return
59   -
  59 +
60 60  
61 61 # ============ #
62 62 # Create nodes #
... ... @@ -72,29 +72,29 @@ def create_nodes(emin, emax, enumbins):
72 72 elogmin = math.log10(float(emin))
73 73 elogmax = math.log10(float(emax))
74 74 elogbin = (elogmax - elogmin)/(float(enumbins)-1.0)
75   -
  75 +
76 76 # Fill arrays
77 77 for i in range(int(enumbins)):
78   -
  78 +
79 79 # Compute energy
80 80 energy = math.pow(10.0, i*elogbin+elogmin)
81   -
  81 +
82 82 # Compute scale (HESS Crab spectrum)
83 83 scale = 3.45e-17 * math.pow(energy, -2.63)
84   -
  84 +
85 85 # Write node
86 86 write_node(file, energy, scale)
87   -
  87 +
88 88 # Debug
89 89 #sys.stdout.write(energy)
90 90 #sys.stdout.write(math.pow(10.0, int(math.log10(scale))-1.0)ยฑ"\n")
91   -
  91 +
92 92 # Close file
93 93 file.close()
94   -
  94 +
95 95 # Return
96 96 return
97   -
  97 +
98 98  
99 99 # ======================= #
100 100 # Main script entry point #
... ... @@ -108,7 +108,7 @@ if __name__ == &#39;__main__&#39;:
108 108 if len(sys.argv) < 3:
109 109 sys.stdout.write(usage+"\n")
110 110 sys.exit()
111   -
  111 +
112 112 # Extract parameters
113 113 emin = sys.argv[1]
114 114 emax = sys.argv[2]
... ... @@ -116,4 +116,3 @@ if __name__ == &#39;__main__&#39;:
116 116  
117 117 # Create nodes
118 118 create_nodes(emin, emax, enumbins)
119   -
... ...
cscripts/cta_root2caldb.py
  View file @5c53c14
... ... @@ -52,7 +52,7 @@ class caldb():
52 52 """
53 53 # Store root path
54 54 self.path = path
55   -
  55 +
56 56 # Initialise some members
57 57 self.cif = None
58 58 self.irf = None
... ... @@ -64,7 +64,7 @@ class caldb():
64 64  
65 65 # Store UTC date of job execution
66 66 self.utc = datetime.utcnow().isoformat()[:19]
67   -
  67 +
68 68 # Set CALDB information
69 69 self.cal_tel = "CTA"
70 70 self.cal_inst = inst.upper()
... ... @@ -102,13 +102,13 @@ class caldb():
102 102 self.bgd_doc = "???"
103 103 self.bgd_bounds = [self.cal_name, self.cal_version, self.cal_cut, self.cal_analysis]
104 104 self.bgd_desc = "CTA background"
105   -
  105 +
106 106 # Create directory structure
107 107 self.make_dirs()
108   -
  108 +
109 109 # Return
110 110 return
111   -
  111 +
112 112 def __del__(self):
113 113 """
114 114 Destructor. Close all FITS files.
... ... @@ -118,18 +118,18 @@ class caldb():
118 118  
119 119 # Return
120 120 return
121   -
  121 +
122 122 def make_dirs(self):
123 123 """
124 124 Generate CALDB directory structure for one observation identifier.
125 125 The structure is given by
126   -
  126 +
127 127 data/<tel>/<inst>/bcf/<obsid>
128   -
  128 +
129 129 where <tel> is "cta" and <inst> is the instrument specified in the
130 130 CALDB constructor (the instrument may be used for different array
131 131 configurations).
132   -
  132 +
133 133 Parameters:
134 134 None
135 135 Keywords:
... ... @@ -147,7 +147,7 @@ class caldb():
147 147 self.psf_dir = self.obs_dir
148 148 self.edisp_dir = self.obs_dir
149 149 self.bgd_dir = self.obs_dir
150   -
  150 +
151 151 # Optionally prefix path
152 152 if len(self.path) > 0:
153 153 self.base_path = self.path+"/"+self.base_dir
... ... @@ -175,17 +175,17 @@ class caldb():
175 175 os.makedirs(self.edisp_path)
176 176 if not os.path.isdir(self.bgd_path):
177 177 os.makedirs(self.bgd_path)
178   -
  178 +
179 179 # Return
180 180 return
181   -
  181 +
182 182 def open(self, version, split=False, clobber=True):
183 183 """
184 184 Open existing or create new calibration. The actual version will put
185 185 all calibrations in the same file, although each part of the response
186 186 function will have its own logical name. We can thus easily modify
187 187 the script to put each calibration information in a separate file.
188   -
  188 +
189 189 Parameters:
190 190 version - Filename version
191 191 Keywords:
... ... @@ -203,17 +203,17 @@ class caldb():
203 203 self.psf_file = "irf_"+version+".fits"
204 204 self.edisp_file = "irf_"+version+".fits"
205 205 self.bgd_file = "irf_"+version+".fits"
206   -
  206 +
207 207 # Open calibration database index
208 208 self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True)
209   -
  209 +
210 210 # If file has no CIF extension than create it now
211 211 try:
212 212 self.hdu_cif = self.cif.table("CIF")
213 213 except:
214 214 self.create_cif()
215 215 self.hdu_cif = self.cif.table("CIF")
216   -
  216 +
217 217 # Set filenames
218 218 ea_filename = self.ea_path+"/"+self.ea_file
219 219 psf_filename = self.psf_path+"/"+self.psf_file
... ... @@ -232,13 +232,13 @@ class caldb():
232 232 self.psf = self.irf
233 233 self.edisp = self.irf
234 234 self.bgd = self.irf
235   -
  235 +
236 236 # Open HDUs
237 237 self.open_ea()
238 238 self.open_psf()
239 239 self.open_edisp()
240 240 self.open_bgd()
241   -
  241 +
242 242 # Return
243 243 return
244 244  
... ... @@ -255,7 +255,7 @@ class caldb():
255 255 # necessarily have all the information at hand (in particular
256 256 # about the boundaries)
257 257 self.add_cif_info()
258   -
  258 +
259 259 # Close CIF
260 260 if self.cif != None:
261 261 self.cif.save(True)
... ... @@ -272,11 +272,11 @@ class caldb():
272 272 self.psf = None
273 273 self.edisp = None
274 274 self.bgd = None
275   -
  275 +
276 276 # ... otherwise we have split files, so we have to close them
277 277 # all
278 278 else:
279   -
  279 +
280 280 # Close effective area file
281 281 if self.ea != None:
282 282 self.ea.save(True)
... ... @@ -300,10 +300,10 @@ class caldb():
300 300 self.bgd.save(True)
301 301 self.bgd.close()
302 302 self.bgd = None
303   -
  303 +
304 304 # Return
305 305 return
306   -
  306 +
307 307 def create_cif(self):
308 308 """
309 309 Create Calibration Index File (CIF) extension in CIF FITS file.
... ... @@ -315,9 +315,9 @@ class caldb():
315 315 """
316 316 # Create binary table
317 317 table = gammalib.GFitsBinTable()
318   -
  318 +
319 319 # Set boundary
320   -
  320 +
321 321 # Attach columns. Reference: CAL/GEN/92-008
322 322 table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10))
323 323 table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10))
... ... @@ -337,20 +337,20 @@ class caldb():
337 337 table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0))
338 338 table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8))
339 339 table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70))
340   -
  340 +
341 341 # Set keywords. Reference: CAL/GEN/92-008
342 342 table.extname("CIF")
343 343 table.card("CIFVERSN", "1992a", "Version of CIF format")
344   -
  344 +
345 345 # Attach table to FITS file. Note that at this moment the
346 346 # FITS table gets connected to the FITS file. Yet since
347 347 # nothing was yet written to the FITS file, we cannot read
348 348 # anything from it.
349 349 self.cif.append(table)
350   -
  350 +
351 351 # Return
352 352 return
353   -
  353 +
354 354 def add_cif_info(self):
355 355 """
356 356 Add information to CIF extension.
... ... @@ -362,13 +362,13 @@ class caldb():
362 362 """
363 363 # Append 4 rows to CIF extension
364 364 self.hdu_cif.append_rows(4)
365   -
  365 +
366 366 # Add generic information for these 4 rows
367 367 for i in range(4):
368   -
  368 +
369 369 # Set row number
370 370 row = i+self.hdu_cif.nrows()-4
371   -
  371 +
372 372 # Set element
373 373 self.hdu_cif["TELESCOP"][row] = self.cal_tel
374 374 self.hdu_cif["INSTRUME"][row] = self.cal_inst
... ... @@ -382,7 +382,7 @@ class caldb():
382 382 self.hdu_cif["REF_TIME"][row] = self.ref_time
383 383 self.hdu_cif["CAL_QUAL"][row] = self.cal_qual
384 384 self.hdu_cif["CAL_DATE"][row] = self.cal_date
385   -
  385 +
386 386 # Add effective area information
387 387 row = self.hdu_cif.nrows()-4
388 388 self.hdu_cif["CAL_DIR"][row] = self.ea_dir
... ... @@ -434,7 +434,7 @@ class caldb():
434 434 self.hdu_cif["CAL_CBD"][row, n] = "NONE"
435 435 else:
436 436 self.hdu_cif["CAL_CBD"][row, n] = self.bgd_bounds[n]
437   -
  437 +
438 438 # Return
439 439 return
440 440  
... ... @@ -455,20 +455,20 @@ class caldb():
455 455 except:
456 456 # Create binary table
457 457 table = gammalib.GFitsBinTable()
458   -
  458 +
459 459 # Set extension name
460 460 table.extname("EFFECTIVE AREA")
461   -
  461 +
462 462 # Set OGIP keywords
463 463 self.set_ogip_keywords(table, self.ea_doc, \
464 464 ["RESPONSE", self.ea_name])
465   -
  465 +
466 466 # Append table
467 467 self.ea.append(table)
468   -
  468 +
469 469 # Get extension
470 470 self.hdu_ea = self.ea.table("EFFECTIVE AREA")
471   -
  471 +
472 472 # Return
473 473 return
474 474  
... ... @@ -489,20 +489,20 @@ class caldb():
489 489 except:
490 490 # Create binary table
491 491 table = gammalib.GFitsBinTable()
492   -
  492 +
493 493 # Set extension name
494 494 table.extname("POINT SPREAD FUNCTION")
495   -
  495 +
496 496 # Set OGIP keywords
497 497 self.set_ogip_keywords(table, self.psf_doc, \
498 498 ["RESPONSE", self.psf_name])
499   -
  499 +
500 500 # Append table
501 501 self.psf.append(table)
502 502  
503 503 # Get extension
504 504 self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION")
505   -
  505 +
506 506 # Return
507 507 return
508 508  
... ... @@ -523,10 +523,10 @@ class caldb():
523 523 except:
524 524 # Create binary table
525 525 table = gammalib.GFitsBinTable()
526   -
  526 +
527 527 # Set extension name
528 528 table.extname("ENERGY DISPERSION")
529   -
  529 +
530 530 # Set OGIP keywords
531 531 self.set_ogip_keywords(table, self.edisp_doc, \
532 532 ["RESPONSE", self.edisp_name])
... ... @@ -536,7 +536,7 @@ class caldb():
536 536  
537 537 # Get extension
538 538 self.hdu_edisp = self.edisp.table("ENERGY DISPERSION")
539   -
  539 +
540 540 # Return
541 541 return
542 542  
... ... @@ -557,10 +557,10 @@ class caldb():
557 557 except:
558 558 # Create binary table
559 559 table = gammalib.GFitsBinTable()
560   -
  560 +
561 561 # Set extension name
562 562 table.extname("BACKGROUND")
563   -
  563 +
564 564 # Set OGIP keywords
565 565 self.set_ogip_keywords(table, self.bgd_doc, \
566 566 ["RESPONSE", self.bgd_name])
... ... @@ -570,7 +570,7 @@ class caldb():
570 570  
571 571 # Get extension
572 572 self.hdu_bgd = self.bgd.table("BACKGROUND")
573   -
  573 +
574 574 # Return
575 575 return
576 576  
... ... @@ -597,7 +597,7 @@ class caldb():
597 597 key = "HDUCLAS%d" % (i+1)
598 598 hdu.card(key, item, "")
599 599 hdu.card("HDUVERS", "1.0.0", "")
600   -
  600 +
601 601 # Return
602 602 return
603 603  
... ... @@ -608,7 +608,7 @@ class caldb():
608 608 offset angle columns, this method will simply add another data
609 609 column.
610 610 If name==None, the method will not append any data column.
611   -
  611 +
612 612 Parameters:
613 613 array - ROOT 2D histogram.
614 614 hdu - FITS HDU.
... ... @@ -690,7 +690,7 @@ class caldb():
690 690 offset angle columns, this method will simply add another data
691 691 column.
692 692 If name==None, the method will not append any data column.
693   -
  693 +
694 694 Parameters:
695 695 array - ROOT 2D histogram.
696 696 hdu - FITS HDU.
... ... @@ -701,7 +701,7 @@ class caldb():
701 701 """
702 702 # Set constants
703 703 deg2sr = 0.01745329*0.01745329
704   -
  704 +
705 705 # Extract energy and offset angle vectors
706 706 energies = array.GetXaxis()
707 707 neng = energies.GetNbins()
... ... @@ -825,7 +825,7 @@ class caldb():
825 825 histogram. If the HDU has already the columns, this method will
826 826 simply add another data column.
827 827 If name==None, the method will not append any data column.
828   -
  828 +
829 829 Parameters:
830 830 array - ROOT 3D histogram.
831 831 hdu - FITS HDU.
... ... @@ -925,11 +925,11 @@ class caldb():
925 925  
926 926 # Return boundary information
927 927 return bounds
928   -
  928 +
929 929 def root2caldb(self, filename, rebin=False, psftype="Gauss", eascale=1.0, bgdscale=1.0):
930 930 """
931 931 Translate ROOT to CALDB information.
932   -
  932 +
933 933 Parameters:
934 934 filename - ROOT 2D performance filename.
935 935 Keywords:
... ... @@ -943,7 +943,7 @@ class caldb():
943 943  
944 944 # Create effective area
945 945 self.root2ea(file, rebin=rebin, eascale=eascale)
946   -
  946 +
947 947 # Create point spread function
948 948 self.root2psf(file, psftype)
949 949  
... ... @@ -953,7 +953,7 @@ class caldb():
953 953 # Create background
954 954 #self.root2bgd(file)
955 955 self.root2bgd3D(file, bgdscale=bgdscale)
956   -
  956 +
957 957 # Return
958 958 return
959 959  
... ... @@ -1076,7 +1076,7 @@ class caldb():
1076 1076 for ieng in range(neng):
1077 1077 sigma = r68.GetBinContent(ieng+1,ioff+1) * r68_to_sigma
1078 1078 r68.SetBinContent(ieng+1,ioff+1,sigma)
1079   -
  1079 +
1080 1080 # Compute scale histogram
1081 1081 scale = r68.Clone()
1082 1082 for ioff in range(noffset):
... ... @@ -1087,7 +1087,7 @@ class caldb():
1087 1087 else:
1088 1088 value = 0.0
1089 1089 scale.SetBinContent(ieng+1,ioff+1,value)
1090   -
  1090 +
1091 1091 # Set zero histogram
1092 1092 zero = r68.Clone()
1093 1093 for ioff in range(noffset):
... ... @@ -1172,12 +1172,12 @@ class caldb():
1172 1172  
1173 1173 # Continue only if both radii are positive
1174 1174 if r_68 > 0 and r_80 > 0:
1175   -
  1175 +
1176 1176 # Derive constants for equation to solve
1177 1177 a = 1.0 - 0.68
1178 1178 b = 1.0 - 0.80
1179 1179 c = r_68*r_68/(r_80*r_80)
1180   -
  1180 +
1181 1181 # Solve equation (a^x-1)/(b^x-1)=c for x using secant
1182 1182 # method. Stop when we are better than 1e-6.
1183 1183 x1 = -0.5
... ... @@ -1202,7 +1202,7 @@ class caldb():
1202 1202 gamma = 1.0 - 1.0/x
1203 1203 else:
1204 1204 gamma = 1.0
1205   -
  1205 +
1206 1206 # Compute sigma
1207 1207 denominator = 2.0 * gamma * (math.pow(a, x) - 1.0)
1208 1208 if denominator > 0.0:
... ... @@ -1228,7 +1228,7 @@ class caldb():
1228 1228  
1229 1229 # Show results on console
1230 1230 #print(ioff,ieng,gamma,sigma,x,f,iter)
1231   -
  1231 +
1232 1232 # Set bin contents
1233 1233 gamma2D.SetBinContent(ieng+1,ioff+1,gamma)
1234 1234 sigma2D.SetBinContent(ieng+1,ioff+1,sigma)
... ... @@ -1352,7 +1352,7 @@ class caldb():
1352 1352  
1353 1353 # BGD_RECO
1354 1354 self.make_2D(array, self.hdu_bgd, "BGD_RECO", "")
1355   -
  1355 +
1356 1356 # Return
1357 1357 return
1358 1358  
... ... @@ -1384,7 +1384,7 @@ class caldb():
1384 1384  
1385 1385 # BGD (reconstructed energy)
1386 1386 self.make_3D(array, self.hdu_bgd, "BGD", "1/s/MeV/sr", scale=bgdscale)
1387   -
  1387 +
1388 1388 # Return
1389 1389 return
1390 1390  
... ... @@ -1397,7 +1397,7 @@ class caldb():
1397 1397 hdu.card("CCLS0001", self.cal_class, "Calibration class")
1398 1398 hdu.card("CDTP0001", self.cal_type, "Calibration type")
1399 1399 hdu.card("CCNM0001", name, "Calibration name")
1400   -
  1400 +
1401 1401 # Set boundary keywords
1402 1402 for n in range(9):
1403 1403 keyname = "CBD%d0001" % (n+1)
... ... @@ -1406,12 +1406,12 @@ class caldb():
1406 1406 else:
1407 1407 value = bounds[n]
1408 1408 hdu.card(keyname, value, "Calibration boundary")
1409   -
  1409 +
1410 1410 # Set validity keywords
1411 1411 hdu.card("CVSD0001", self.val_date, "Calibration validity start date (UTC)")
1412 1412 hdu.card("CVST0001", self.val_time, "Calibration validity start time (UTC)")
1413 1413 hdu.card("CDEC0001", desc, "Calibration description")
1414   -
  1414 +
1415 1415 # Return
1416 1416 return
1417 1417  
... ... @@ -1551,7 +1551,7 @@ if __name__ == &#39;__main__&#39;:
1551 1551 This script translates a CTA ROOT 2D performance files into a CALDB
1552 1552 compliant CTA response. The script will create the response for 3
1553 1553 observations to illustrate the run-wise file structure.
1554   -
  1554 +
1555 1555 The CTA ROOT 2D performance files contain the following histograms:
1556 1556 - DiffSens_offaxis (2D) - Differential sensitivity
1557 1557 - EffectiveArea_offaxis (2D) - Effective area (full containment?, reco energy?)
... ... @@ -1592,9 +1592,9 @@ if __name__ == &#39;__main__&#39;:
1592 1592  
1593 1593 # Allocate caldb entry
1594 1594 irf = caldb(inst, id)
1595   -
  1595 +
1596 1596 # Open calibration file
1597 1597 irf.open("file")
1598   -
  1598 +
1599 1599 # Translate ROOT to CALDB information
1600 1600 irf.root2caldb(filename, rebin=rebin, psftype=psftype, eascale=eascale, bgdscale=bgdscale)
... ...
cscripts/cta_root2irf.py
  View file @5c53c14
... ... @@ -40,7 +40,7 @@ import os
40 40 def sigma_lima(Non, Noff, alpha=0.2):
41 41 """
42 42 Compute Eq. (17) of Li & Ma (1983).
43   -
  43 +
44 44 Parameters:
45 45 Non - Number of on counts
46 46 Noff - Number of off counts
... ... @@ -55,7 +55,7 @@ def sigma_lima(Non, Noff, alpha=0.2):
55 55 term1 = Non * math.log((alpha1/alpha)*arg1)
56 56 term2 = Noff * math.log(alpha1*arg2)
57 57 sigma = math.sqrt(2.0 * (term1 + term2))
58   -
  58 +
59 59 # Return sensitivity
60 60 return sigma
61 61  
... ... @@ -66,7 +66,7 @@ def sigma_lima(Non, Noff, alpha=0.2):
66 66 def Non_lima(sigma, Noff, alpha=0.2):
67 67 """
68 68 Solve Eq. (17) of Li & Ma (1983) for Non
69   -
  69 +
70 70 Parameters:
71 71 sigma - Required significance
72 72 Noff - Number of off counts
... ... @@ -76,18 +76,18 @@ def Non_lima(sigma, Noff, alpha=0.2):
76 76 # Set initial guess for Non
77 77 S = 25.0
78 78 Non = S + alpha*Noff
79   -
  79 +
80 80 # Compute initial sigma estimate
81 81 s = sigma_lima(Non, Noff, alpha=alpha)
82 82 delta = s - sigma
83   -
  83 +
84 84 # Iterate
85 85 while abs(delta) > 1.0e-6:
86 86 S *= sigma/s
87 87 Non = S + alpha*Noff
88 88 s = sigma_lima(Non, Noff, alpha=alpha)
89 89 delta = s - sigma
90   -
  90 +
91 91 # Return
92 92 return Non
93 93  
... ... @@ -101,7 +101,7 @@ def get_aeff(bgrate, emin, emax, \
101 101 Computes effective area assuming a Crab spectrum. This does not
102 102 include any hadron cut efficiencies, and thus is not very useful in a first
103 103 place ...
104   -
  104 +
105 105 Parameters:
106 106 bgrate - Background rate (TeV)
107 107 emin - Minimum energy (TeV)
... ... @@ -118,14 +118,14 @@ def get_aeff(bgrate, emin, emax, \
118 118 Noff = B/alpha
119 119 Non = Non_lima(sigma, Noff)
120 120 S = (Non - B)
121   -
  121 +
122 122 # Compute Crab flux in specified energy band
123 123 inx = index+1.0
124 124 cu = norm/inx * (math.pow(emax, inx) - math.pow(emin, inx))
125   -
  125 +
126 126 # Compute effective area
127 127 area = S / duration / cu
128   -
  128 +
129 129 # Return effective area
130 130 return area
131 131  
... ... @@ -136,7 +136,7 @@ def get_aeff(bgrate, emin, emax, \
136 136 def make_files(filename, prfname, bkgname):
137 137 """
138 138 Generate files.
139   -
  139 +
140 140 Parameters:
141 141 filename - Performance file name
142 142 prfname - Performance file name
... ... @@ -146,7 +146,7 @@ def make_files(filename, prfname, bkgname):
146 146 file = TFile(filename)
147 147 fprf = open(prfname, "w")
148 148 fbkg = open(bkgname, "w")
149   -
  149 +
150 150 # Write header
151 151 fprf.write("log(E) Area r68 r80 ERes. BG Rate Diff Sens\n")
152 152 #sys.stdout.write("log(E) Area Area80 r68 r80 ERes. BG Rate Diff Sens BGRateSqDeg BGInt (BGIntControl)\n")
... ... @@ -173,12 +173,12 @@ def make_files(filename, prfname, bkgname):
173 173 file.GetObject("AngRes", angres)
174 174 file.GetObject("AngRes80", angres80)
175 175 file.GetObject("ERes", eres)
176   -
  176 +
177 177 # Setup vectors
178 178 energies = bgrate.GetXaxis()
179 179 nbins_eng = energies.GetNbins()
180 180 for i in range(nbins_eng):
181   -
  181 +
182 182 # Get logE and background rate per squared degree
183 183 logE = energies.GetBinCenter(i+1)
184 184 bgd = bgrate.GetBinContent(i+1,1)
... ... @@ -189,7 +189,7 @@ def make_files(filename, prfname, bkgname):
189 189 r80 = angres80.GetBinContent(i+1,1)
190 190 dE = eres.GetBinContent(i+1,1)
191 191 diffs = sens.GetBinContent(i+1,1)
192   -
  192 +
193 193 # Skip any NaNs
194 194 if math.isnan(area80):
195 195 continue
... ... @@ -199,10 +199,10 @@ def make_files(filename, prfname, bkgname):
199 199 emin = math.pow(10.0, logE-0.1)*1.0e6
200 200 emax = math.pow(10.0, logE+0.1)*1.0e6
201 201 ewidth = emax-emin
202   -
  202 +
203 203 # Convert into background rate per steradian and MeV
204 204 bkg_rate = bgdsq / (0.01745329*0.01745329) / ewidth
205   -
  205 +
206 206 # Compute control background rate (this only works for KB files
207 207 # as they are for 80% containment radius). But we don't use this
208 208 # anyways, it's just for control on display ...
... ... @@ -211,7 +211,7 @@ def make_files(filename, prfname, bkgname):
211 211 bkg_rate2 = bgd / omega / ewidth
212 212 else:
213 213 bkg_rate2 = 0.0
214   -
  214 +
215 215 # Compute full effective area by dividing area80/0.80
216 216 area = area80/0.80
217 217  
... ... @@ -241,7 +241,7 @@ def make_files(filename, prfname, bkgname):
241 241 # Close files
242 242 fprf.close()
243 243 fbkg.close()
244   -
  244 +
245 245 # Return
246 246 return
247 247  
... ... @@ -263,13 +263,13 @@ if __name__ == &#39;__main__&#39;:
263 263 # Get list of all performance files
264 264 files = glob.glob(caldb+"/*.root")
265 265 files.sort()
266   -
  266 +
267 267 # Loop over all files
268 268 for file in files:
269   -
  269 +
270 270 # Get filename
271 271 head, tail = os.path.split(file)
272   -
  272 +
273 273 # Build name
274 274 name = tail.strip(".root")
275 275 if name.find("kb_") != -1:
... ... @@ -281,7 +281,7 @@ if __name__ == &#39;__main__&#39;:
281 281 i = name.find("hours")
282 282 if i != -1:
283 283 name = name[0:i+1]
284   -
  284 +
285 285 # Build output filenames
286 286 prfname = name + ".dat"
287 287 bkgname = "bkg_"+name+".txt"
... ...
cscripts/lsmodel.py
  View file @5c53c14
... ... @@ -39,7 +39,7 @@ class lsmodel(gammalib.GApplication):
39 39 # Set name
40 40 self.name = "lsmodel"
41 41 self.version = "0.1.0"
42   -
  42 +
43 43 # Make sure that parfile exists
44 44 file = self.parfile()
45 45  
... ... @@ -57,7 +57,7 @@ class lsmodel(gammalib.GApplication):
57 57  
58 58 # Return
59 59 return
60   -
  60 +
61 61 def __del__(self):
62 62 """
63 63 Destructor.
... ... @@ -65,7 +65,7 @@ class lsmodel(gammalib.GApplication):
65 65 # Write separator into logger
66 66 if self.logTerse():
67 67 self.log("\n")
68   -
  68 +
69 69 # Return
70 70 return
71 71  
... ... @@ -76,13 +76,13 @@ class lsmodel(gammalib.GApplication):
76 76 """
77 77 # Set parfile name
78 78 parfile = self.name+".par"
79   -
  79 +
80 80 try:
81 81 pars = gammalib.GApplicationPars(parfile)
82 82 except:
83 83 # Signal if parfile was not found
84 84 sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
85   -
  85 +
86 86 # Create default parfile
87 87 pars = gammalib.GApplicationPars()
88 88 pars.append(gammalib.GApplicationPar("cntmap","f","a","cntmap.fits","","","Counts map"))
... ... @@ -94,10 +94,10 @@ class lsmodel(gammalib.GApplication):
94 94 pars.append(gammalib.GApplicationPar("outfile","f","a","model_cube.fits","","","Output file name"))
95 95 pars.append_standard()
96 96 pars.save(parfile)
97   -
  97 +
98 98 # Return
99 99 return
100   -
  100 +
101 101 def get_parameters(self):
102 102 """
103 103 Get parameters from parfile and setup the observation.
... ... @@ -114,17 +114,17 @@ class lsmodel(gammalib.GApplication):
114 114 # Set some fixed parameters
115 115 self.m_log = False # Logging in client tools
116 116 self.m_debug = False # Debugging in client tools
117   -
  117 +
118 118 # Return
119 119 return
120   -
  120 +
121 121 def execute(self):
122 122 """
123 123 Execute the script.
124 124 """
125 125 # Run the script
126 126 self.run()
127   -
  127 +
128 128 # Return
129 129 return
130 130  
... ... @@ -138,12 +138,12 @@ class lsmodel(gammalib.GApplication):
138 138  
139 139 # Get parameters
140 140 self.get_parameters()
141   -
  141 +
142 142 # Write input parameters into logger
143 143 if self.logTerse():
144 144 self.log_parameters()
145 145 self.log("\n")
146   -
  146 +
147 147 # Write header
148 148 if self.logTerse():
149 149 self.log("\n")
... ... @@ -194,10 +194,9 @@ if __name__ == &#39;__main__&#39;:
194 194 """
195 195 # Create instance of application
196 196 app = lsmodel(sys.argv)
197   -
  197 +
198 198 # Open logfile
199 199 app.logFileOpen()
200   -
  200 +
201 201 # Execute application
202 202 app.execute()
203   -
204 203 \ No newline at end of file
... ...
cscripts/obsutils.py
  View file @5c53c14
... ... @@ -30,7 +30,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
30 30 binsz=0.05, npix=200, proj="TAN", coord="GAL"):
31 31 """
32 32 Simulate events for all observations in the container.
33   -
  33 +
34 34 Parameters:
35 35 obs - Observation container
36 36 Keywords:
... ... @@ -42,12 +42,12 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
42 42 binsz - Pixel size for binned simulation (deg/pixel)
43 43 npix - Number of pixels in X and Y for binned simulation
44 44 """
45   -
  45 +
46 46 # Allocate ctobssim application and set parameters
47 47 sim = ctools.ctobssim(obs)
48 48 sim["seed"] = seed
49 49 sim["edisp"] = edisp
50   -
  50 +
51 51 # Optionally open the log file
52 52 if log:
53 53 sim.logFileOpen()
... ... @@ -58,16 +58,16 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
58 58  
59 59 # Set chatter level
60 60 sim["chatter"] = chatter
61   -
  61 +
62 62 # Run ctobssim application. This will loop over all observations in the
63 63 # container and simulation the events for each observation. Note that
64 64 # events are not added together, they still apply to each observation
65 65 # separately.
66 66 sim.run()
67   -
  67 +
68 68 # Binned option?
69 69 if nbins > 0:
70   -
  70 +
71 71 # Determine common energy boundaries for observations
72 72 emin = None
73 73 emax = None
... ... @@ -82,7 +82,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
82 82 emax = run_emax
83 83 elif run_emax > emax:
84 84 emax = run_emax
85   -
  85 +
86 86 # Allocate ctbin application and set parameters
87 87 bin = ctools.ctbin(sim.obs())
88 88 bin["ebinalg"] = "LOG"
... ... @@ -95,7 +95,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
95 95 bin["binsz"] = binsz
96 96 bin["coordsys"] = coord
97 97 bin["proj"] = proj
98   -
  98 +
99 99 # Optionally open the log file
100 100 if log:
101 101 bin.logFileOpen()
... ... @@ -110,22 +110,22 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
110 110 # Run ctbin application. This will loop over all observations in
111 111 # the container and bin the events in counts maps
112 112 bin.run()
113   -
  113 +
114 114 # Make a deep copy of the observation that will be returned
115 115 # (the ctbin object will go out of scope one the function is
116 116 # left)
117 117 obs = bin.obs().copy()
118 118  
119 119 else:
120   -
  120 +
121 121 # Make a deep copy of the observation that will be returned
122 122 # (the ctobssim object will go out of scope one the function is
123 123 # left)
124 124 obs = sim.obs().copy()
125   -
  125 +
126 126 # Delete the simulation
127 127 del sim
128   -
  128 +
129 129 # Return observation container
130 130 return obs
131 131  
... ... @@ -136,7 +136,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
136 136 def fit(obs, log=False, debug=False, chatter=2, edisp=False):
137 137 """
138 138 Perform maximum likelihood fitting of observations in the container.
139   -
  139 +
140 140 Parameters:
141 141 obs - Observation container
142 142 Keywords:
... ... @@ -147,24 +147,24 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False):
147 147 """
148 148 # Allocate ctlike application
149 149 like = ctools.ctlike(obs)
150   -
  150 +
151 151 # Optionally open the log file
152 152 if log:
153 153 like.logFileOpen()
154   -
  154 +
155 155 # Optionally switch-on debugging model
156 156 if debug:
157 157 like["debug"] = True
158 158  
159 159 # Set chatter level
160 160 like["chatter"] = chatter
161   -
  161 +
162 162 # Optionally apply energy dispersion
163 163 like["edisp"] = edisp
164 164  
165 165 # Run ctlike application.
166 166 like.run()
167   -
  167 +
168 168 # Return observations
169 169 return like
170 170  
... ... @@ -175,7 +175,7 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False):
175 175 def cterror(obs, srcname, log=False, debug=False, chatter=2):
176 176 """
177 177 Perform maximum likelihood fitting of observations in the container.
178   -
  178 +
179 179 Parameters:
180 180 obs - Observation container
181 181 srcname - Source name
... ... @@ -190,11 +190,11 @@ def cterror(obs, srcname, log=False, debug=False, chatter=2):
190 190  
191 191 # Set cterror parameters
192 192 error["srcname"] = srcname
193   -
  193 +
194 194 # Optionally open the log file
195 195 if log:
196 196 error.logFileOpen()
197   -
  197 +
198 198 # Optionally switch-on debugging model
199 199 if debug:
200 200 error["debug"] = True
... ... @@ -204,7 +204,7 @@ def cterror(obs, srcname, log=False, debug=False, chatter=2):
204 204  
205 205 # Run cterror application.
206 206 error.run()
207   -
  207 +
208 208 # Return observations
209 209 return error
210 210  
... ... @@ -218,7 +218,7 @@ def cntmap(obs, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
218 218 """
219 219 Creates a counts map by combining the events of all observations.
220 220 The counts map will be a summed map over all energies.
221   -
  221 +
222 222 Parameters:
223 223 obs - Observation container
224 224 Keywords:
... ... @@ -236,25 +236,25 @@ def cntmap(obs, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
236 236  
237 237 # Set maximum pixel number
238 238 maxpixel = nxpix * nypix
239   -
  239 +
240 240 # Fill all observations
241 241 for run in obs:
242   -
  242 +
243 243 # Loop over all events
244 244 for event in run.events():
245   -
  245 +
246 246 # Determine sky pixel
247 247 skydir = event.dir().dir()
248 248 pixel = map.dir2inx(skydir)
249   -
  249 +
250 250 # Set pixel
251 251 if pixel < maxpixel:
252 252 map[pixel] += event.counts()
253   -
  253 +
254 254 # Save sky map. The clobber flag is set to True, so any existing FITS
255 255 # file will be overwritten.
256 256 map.save(outname, True)
257   -
  257 +
258 258 # Return counts map
259 259 return map
260 260  
... ... @@ -269,7 +269,7 @@ def modmap(obs, eref=0.1, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
269 269 Make model map for a given reference energy by combining all observations.
270 270 The model map will be evaluated for a given reference energy 'eref' and will
271 271 be given in units of [counts/(sr MeV s)].
272   -
  272 +
273 273 Parameters:
274 274 obs - Observation container
275 275 Keywords:
... ... @@ -285,7 +285,7 @@ def modmap(obs, eref=0.1, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
285 285 """
286 286 # Allocate model map
287 287 map = gammalib.GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1)
288   -
  288 +
289 289 # Set reference energy, time and direction. The time is not initialised and is
290 290 # in fact not used (as the IRF is assumed to be time independent for now).
291 291 # The sky direction is set later using the pixel values.
... ... @@ -293,33 +293,33 @@ def modmap(obs, eref=0.1, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
293 293 time = gammalib.GTime()
294 294 instdir = gammalib.GCTAInstDir()
295 295 energy.TeV(eref)
296   -
  296 +
297 297 # Loop over all map pixels
298 298 for pixel in range(map.npix()):
299   -
  299 +
300 300 # Get sky direction
301 301 skydir = map.inx2dir(pixel)
302 302 instdir.dir(skydir)
303   -
  303 +
304 304 # Create event atom for map pixel
305 305 atom = gammalib.GCTAEventAtom()
306 306 atom.dir(instdir)
307 307 atom.energy(energy)
308 308 atom.time(time)
309   -
  309 +
310 310 # Initialise model value
311 311 value = 0.0
312   -
  312 +
313 313 # Loop over all observations
314 314 for run in obs:
315 315 value += obs.models().eval(atom, run)
316   -
  316 +
317 317 # Set map value
318 318 map[pixel] = value
319   -
  319 +
320 320 # Save sky map
321 321 map.save(outname, True)
322   -
  322 +
323 323 # Return model map
324 324 return map
325 325  
... ... @@ -335,15 +335,15 @@ def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
335 335 """
336 336 # Print warning
337 337 print("Warning: obsutils.set is obsolete, use obsutils.set_obs instead.")
338   -
  338 +
339 339 # Call new function
340 340 obs_cta = set_obs(pntdir, tstart=tstart, duration=duration, deadc=deadc, \
341 341 emin=emin, emax=emax, rad=rad, \
342 342 irf=irf, caldb=caldb)
343   -
  343 +
344 344 # Return CTA observation
345 345 return obs_cta
346   -
  346 +
347 347  
348 348 # ======================= #
349 349 # Set one CTA observation #
... ... @@ -355,7 +355,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
355 355 Returns a single CTA observation containing an empty CTA event list.
356 356 By looping over this function you can add CTA observations to the
357 357 observation container.
358   -
  358 +
359 359 Parameters:
360 360 pntdir - Pointing direction [GSkyDir]
361 361 Keywords:
... ... @@ -379,23 +379,23 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
379 379 db.rootdir(caldb)
380 380 else:
381 381 db.open("cta", caldb)
382   -
  382 +
383 383 # Set pointing direction
384 384 pnt = gammalib.GCTAPointing()
385 385 pnt.dir(pntdir)
386 386 obs_cta.pointing(pnt)
387   -
  387 +
388 388 # Set ROI
389 389 roi = gammalib.GCTARoi()
390 390 instdir = gammalib.GCTAInstDir()
391 391 instdir.dir(pntdir)
392 392 roi.centre(instdir)
393 393 roi.radius(rad)
394   -
  394 +
395 395 # Set GTI
396 396 gti = gammalib.GGti()
397 397 gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration))
398   -
  398 +
399 399 # Set energy boundaries
400 400 ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \
401 401 gammalib.GEnergy(emax, "TeV"))
... ... @@ -406,16 +406,16 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
406 406 events.gti(gti)
407 407 events.ebounds(ebounds)
408 408 obs_cta.events(events)
409   -
  409 +
410 410 # Set instrument response
411 411 obs_cta.response(irf, db)
412   -
  412 +
413 413 # Set ontime, livetime, and deadtime correction factor
414 414 obs_cta.ontime(duration)
415 415 obs_cta.livetime(duration*deadc)
416 416 obs_cta.deadc(deadc)
417 417 obs_cta.id(id)
418   -
  418 +
419 419 # Return CTA observation
420 420 return obs_cta
421 421  
... ... @@ -434,7 +434,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
434 434 over other observation proposerties, such as duration, deadtime correction,
435 435 energy range, etc. If an optional keyword is not specified, the function
436 436 keyword is used instead.
437   -
  437 +
438 438 Parameters:
439 439 obsdeflist - Observation definition list [{'ra': x.xx, 'dec': x.xx}]
440 440 The directory can take the following optional keywords:
... ... @@ -470,7 +470,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
470 470 # Set pointing direction
471 471 pntdir = gammalib.GSkyDir()
472 472 pntdir.radec_deg(obsdef['ra'], obsdef['dec'])
473   -
  473 +
474 474 # Set duration (use default if not found in definition list)
475 475 if 'duration' in obsdef:
476 476 obs_duration = obsdef['duration']
... ... @@ -512,7 +512,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
512 512 obs_irf = obsdef['irf']
513 513 else:
514 514 obs_irf = irf
515   -
  515 +
516 516 # Generate identifier string
517 517 id = "%6.6d" % obs_id
518 518  
... ... @@ -520,14 +520,14 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
520 520 obs_cta = set_obs(pntdir, tstart=obs_start, duration=obs_duration, \
521 521 deadc=obs_deadc, emin=obs_emin, emax=obs_emax, \
522 522 rad=obs_rad, irf=obs_irf, caldb=obs_caldb, id=id)
523   -
  523 +
524 524 # Append to container
525 525 obs.append(obs_cta)
526 526  
527 527 # Update start time and identifier
528 528 obs_start += obs_duration
529 529 obs_id += 1
530   -
  530 +
531 531 # Return observation container
532 532 return obs
533 533  
... ... @@ -543,7 +543,7 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
543 543 pattern - Observation pattern. Possible options are:
544 544 - "single": single pointing
545 545 - "four": four pointings 'offset' around pattern centre
546   -
  546 +
547 547 Keywords:
548 548 ra - Right Ascension of pattern centre [deg] (default: 83.6331)
549 549 dec - Declination of pattern centre [deg] (default: 22.0145)
... ... @@ -551,7 +551,7 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
551 551 """
552 552 # Initialise observation definition list
553 553 obsdeflist = []
554   -
  554 +
555 555 # Add patterns
556 556 if pattern == "single":
557 557 obsdef = {'ra': ra, 'dec': dec}
... ... @@ -560,7 +560,7 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
560 560 # Set pattern centre
561 561 centre = gammalib.GSkyDir()
562 562 centre.radec_deg(ra, dec)
563   -
  563 +
564 564 # Append pointings
565 565 for phi in [0.0, 90.0, 180.0, 270.0]:
566 566 pntdir = centre.copy()
... ... @@ -568,6 +568,6 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
568 568 obsdeflist.append({'ra': pntdir.ra_deg(), 'dec': pntdir.dec_deg()})
569 569 else:
570 570 print("Warning: Observation pattern '"+str(pattern)+"' not recognized.")
571   -
  571 +
572 572 # Return observation definition list
573 573 return obsdeflist
... ...
cscripts/tsmerge.py
  View file @5c53c14
... ... @@ -25,27 +25,27 @@ def print_help():
25 25 # Maps class holding the different sky maps
26 26 # created by cttsmap
27 27 class maps:
28   -
  28 +
29 29 # Constructor
30 30 def __init__(self,fitsfile):
31   -
  31 +
32 32 self.filename = fitsfile
33 33 if not os.path.isfile(fitsfile):
34 34 self.error("Maps instance cannot be constructed: File \""+self.filename+"\" does not exist")
35   -
  35 +
36 36 # open FITS file
37 37 fits = gammalib.GFits(fitsfile)
38   -
  38 +
39 39 # get TS map
40 40 self.tsmap = gammalib.GSkymap()
41 41 self.tsmap.read(fits[0])
42   -
  42 +
43 43 # Get status map
44 44 if not fits.contains("STATUS MAP"):
45 45 self.error("No extension \"STATUS MAP\" found in \""+self.filename+"\"")
46 46 self.statusmap = gammalib.GSkymap()
47 47 self.statusmap.read(fits["STATUS MAP"])
48   -
  48 +
49 49 # Get other maps
50 50 self.maps = []
51 51 self.mapnames = []
... ... @@ -72,30 +72,30 @@ class maps:
72 72 check = False
73 73 break
74 74 return check
75   -
  75 +
76 76 # Save maps to outfile
77 77 def save(self,outfile):
78   -
  78 +
79 79 fits = gammalib.GFits()
80 80 self.tsmap.write(fits)
81   -
  81 +
82 82 for i in range(len(self.maps)):
83 83 self.maps[i].write(fits)
84   -
  84 +
85 85 for i in range(len(self.mapnames)):
86 86 fits[i+1].extname(self.mapnames[i])
87   -
  87 +
88 88 if not self.check_status():
89 89 self.statusmap.write(fits)
90 90 fits[fits.size()-1].extname("STATUS MAP")
91 91 fits.saveto(outfile,True)
92   -
  92 +
93 93 # Add map instances to each other
94 94 # in this way the TS maps get merged
95 95 def add(self,addmaps):
96 96 if not len(addmaps) == len(self.maps):
97 97 self.error("Maps cannot be added: Maps from \""+self.filename+"\" do not match from \""+addmaps.filename+"\"")
98   -
  98 +
99 99 # Loop over bins
100 100 for i in range(self.tsmap.npix()):
101 101 if addmaps.statusmap[i] > 0.5:
... ... @@ -107,8 +107,8 @@ class maps:
107 107 def error(self,message):
108 108 print("*** ERROR ***: "+message)
109 109 sys.exit(-1)
110   -
111   -
  110 +
  111 +
112 112 def pars_from_argv(argv):
113 113  
114 114 if argv is None:
... ... @@ -139,10 +139,10 @@ def get_kwargs(argv):
139 139 if s.count('=') == 1:
140 140 key, value = s.split('=', 1)
141 141 kwargs[key] = value
142   -
  142 +
143 143 else:
144 144 print_help()
145   -
  145 +
146 146 return kwargs
147 147  
148 148 if __name__ == "__main__":
... ... @@ -158,10 +158,10 @@ if __name__ == &quot;__main__&quot;:
158 158 allfiles = open(files).read().splitlines()
159 159 else:
160 160 allfiles = glob.glob(files)
161   -
  161 +
162 162 # Initialise files to be merged
163 163 mergefiles = []
164   -
  164 +
165 165 # Test files for the entry status map
166 166 # and add them to the file list
167 167 for fitsfile in allfiles:
... ... @@ -169,14 +169,13 @@ if __name__ == &quot;__main__&quot;:
169 169 if fits.contains("STATUS MAP"):
170 170 mergefiles.append(fitsfile)
171 171 fits.close()
172   -
  172 +
173 173 # create a maps instance from the first file
174 174 finalmaps = maps(mergefiles[0])
175   -
  175 +
176 176 # add other files to the maps instance
177 177 for fitsfile in mergefiles[1:]:
178 178 finalmaps.add(maps(fitsfile))
179   -
  179 +
180 180 # save the filled maps to file
181 181 finalmaps.save(outfile=outfile)
182   -
... ...
examples/benchmark_cta_analysis.py
  View file @5c53c14
... ... @@ -88,7 +88,7 @@ def unbinned_pipeline(duration):
88 88 tstop = time.clock()
89 89 telapsed = tstop - tstart
90 90 tctlike = tstop - tctlike
91   -
  91 +
92 92 # Return
93 93 return telapsed, tctlike
94 94  
... ... @@ -161,7 +161,7 @@ def binned_pipeline(duration):
161 161 tstop = time.clock()
162 162 telapsed = tstop - tstart
163 163 tctlike = tstop - tctlike
164   -
  164 +
165 165 # Return
166 166 return telapsed, tctlike
167 167  
... ... @@ -293,7 +293,7 @@ def stacked_pipeline(duration):
293 293 tstop = time.clock()
294 294 telapsed = tstop - tstart
295 295 tctlike = tstop - tctlike
296   -
  296 +
297 297 # Return
298 298 return telapsed, tctlike
299 299  
... ... @@ -344,7 +344,7 @@ if __name__ == &#39;__main__&#39;:
344 344  
345 345 # Increment duration
346 346 duration *= 2.0
347   -
  347 +
348 348 # Optionally plot results
349 349 if has_matplotlib:
350 350 plt.figure(1)
... ... @@ -359,4 +359,3 @@ if __name__ == &#39;__main__&#39;:
359 359 plt.ylabel("CPU time (seconds)")
360 360 plt.legend(loc="lower right")
361 361 plt.show()
362   -
363 362 \ No newline at end of file
... ...
examples/make_pointings.py
  View file @5c53c14
... ... @@ -330,7 +330,7 @@ def write_obsdef(filename, obsdef):
330 330  
331 331 # Return
332 332 return
333   -
  333 +
334 334  
335 335 # ======================== #
336 336 # Main routine entry point #
... ... @@ -341,7 +341,7 @@ if __name__ == &#39;__main__&#39;:
341 341 """
342 342 # Initialise flags
343 343 need_help = False
344   -
  344 +
345 345 # Test for command line arguments
346 346 if (len(sys.argv) > 1):
347 347 if sys.argv[1] == "-h":
... ...
examples/make_pull_at_sensitivity_limit.py
  View file @5c53c14
... ... @@ -33,10 +33,10 @@ from cscripts import cspull
33 33 # Try importing processing module
34 34 has_processing = False
35 35 try:
36   - import processing
37   - has_processing = True
  36 + import processing
  37 + has_processing = True
38 38 except:
39   - pass
  39 + pass
40 40  
41 41  
42 42 # ======================== #
... ... @@ -44,170 +44,170 @@ except:
44 44 # ======================== #
45 45 def create_pull(loge, emin, emax, models, ntrials=100, duration=180000.0,
46 46 enumbins=0, log=False):
47   - """
48   - Create pull distribution
49   -
50   - Parameters:
51   - emin - Minimum energy (TeV)
52   - emax - Maximum energy (TeV)
53   - Keywords:
54   - log - Create log file(s)
55   - """
56   - # Generate output filename
57   - outfile = "pull_"+loge+".dat"
58   -
59   - # Setup cspull tool
60   - pull = cspull()
61   - pull.models(models)
62   - pull["outfile"] = outfile
63   - pull["ntrials"] = ntrials
64   - pull["caldb"] = "prod2"
65   - pull["irf"] = "South_50h"
66   - pull["ra"] = 83.6331
67   - pull["dec"] = 22.0145
68   - pull["emin"] = emin
69   - pull["emax"] = emax
70   - pull["enumbins"] = enumbins
71   - pull["duration"] = duration
72   -
73   - # Optionally open the log file
74   - if log:
75   - pull.logFileOpen()
76   -
77   - # Run tool
78   - pull.run()
79   -
80   - # Return
81   - return
  47 + """
  48 + Create pull distribution
  49 +
  50 + Parameters:
  51 + emin - Minimum energy (TeV)
  52 + emax - Maximum energy (TeV)
  53 + Keywords:
  54 + log - Create log file(s)
  55 + """
  56 + # Generate output filename
  57 + outfile = "pull_"+loge+".dat"
  58 +
  59 + # Setup cspull tool
  60 + pull = cspull()
  61 + pull.models(models)
  62 + pull["outfile"] = outfile
  63 + pull["ntrials"] = ntrials
  64 + pull["caldb"] = "prod2"
  65 + pull["irf"] = "South_50h"
  66 + pull["ra"] = 83.6331
  67 + pull["dec"] = 22.0145
  68 + pull["emin"] = emin
  69 + pull["emax"] = emax
  70 + pull["enumbins"] = enumbins
  71 + pull["duration"] = duration
  72 +
  73 + # Optionally open the log file
  74 + if log:
  75 + pull.logFileOpen()
  76 +
  77 + # Run tool
  78 + pull.run()
  79 +
  80 + # Return
  81 + return
82 82  
83 83  
84 84 # ============ #
85 85 # Create model #
86 86 # ============ #
87 87 def create_model(flux, index=-2.48, fitidx=False):
88   - """
89   - Add standard CTA background model to observations container.
90   - We use a simple power law here, scaled to Konrad's E configuration
91   - performance table.
92   -
93   - Parameters:
94   - flux - Flux in Crab units
95   - """
96   - # Define background model
97   - bgd_radial = gammalib.GCTAModelRadialGauss(3.0)
98   - bgd_spectrum = gammalib.GModelSpectralPlaw(61.8, -1.85)
99   - bgd_spectrum["Prefactor"].scale(1.0e-6)
100   - bgd_spectrum["PivotEnergy"].value(1.0)
101   - bgd_spectrum["PivotEnergy"].scale(1.0e6)
102   - if fitidx:
103   - bgd_spectrum["Index"].free()
104   - else:
105   - bgd_spectrum["Index"].fix()
106   - bgd_model = gammalib.GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum)
107   - bgd_model.name("Background")
108   - bgd_model.instruments("CTA")
109   -
110   - # Define source spectrum
111   - location = gammalib.GSkyDir()
112   - location.radec_deg(83.6331, 22.0145)
113   - src_spatial = gammalib.GModelSpatialPtsrc(location)
114   - src_spectrum = gammalib.GModelSpectralPlaw(flux, index)
115   - src_spectrum["Prefactor"].scale(5.7e-16)
116   - src_spectrum["PivotEnergy"].value(0.3)
117   - src_spectrum["PivotEnergy"].scale(1.0e6)
118   - if fitidx:
119   - src_spectrum["Index"].free()
120   - else:
121   - src_spectrum["Index"].fix()
122   - src_model = gammalib.GModelPointSource(src_spatial, src_spectrum)
123   - src_model.name("Test")
124   -
125   - # Add models to container
126   - models = gammalib.GModels()
127   - models.append(bgd_model)
128   - models.append(src_model)
129   -
130   - # Return model container
131   - return models
  88 + """
  89 + Add standard CTA background model to observations container.
  90 + We use a simple power law here, scaled to Konrad's E configuration
  91 + performance table.
  92 +
  93 + Parameters:
  94 + flux - Flux in Crab units
  95 + """
  96 + # Define background model
  97 + bgd_radial = gammalib.GCTAModelRadialGauss(3.0)
  98 + bgd_spectrum = gammalib.GModelSpectralPlaw(61.8, -1.85)
  99 + bgd_spectrum["Prefactor"].scale(1.0e-6)
  100 + bgd_spectrum["PivotEnergy"].value(1.0)
  101 + bgd_spectrum["PivotEnergy"].scale(1.0e6)
  102 + if fitidx:
  103 + bgd_spectrum["Index"].free()
  104 + else:
  105 + bgd_spectrum["Index"].fix()
  106 + bgd_model = gammalib.GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum)
  107 + bgd_model.name("Background")
  108 + bgd_model.instruments("CTA")
  109 +
  110 + # Define source spectrum
  111 + location = gammalib.GSkyDir()
  112 + location.radec_deg(83.6331, 22.0145)
  113 + src_spatial = gammalib.GModelSpatialPtsrc(location)
  114 + src_spectrum = gammalib.GModelSpectralPlaw(flux, index)
  115 + src_spectrum["Prefactor"].scale(5.7e-16)
  116 + src_spectrum["PivotEnergy"].value(0.3)
  117 + src_spectrum["PivotEnergy"].scale(1.0e6)
  118 + if fitidx:
  119 + src_spectrum["Index"].free()
  120 + else:
  121 + src_spectrum["Index"].fix()
  122 + src_model = gammalib.GModelPointSource(src_spatial, src_spectrum)
  123 + src_model.name("Test")
  124 +
  125 + # Add models to container
  126 + models = gammalib.GModels()
  127 + models.append(bgd_model)
  128 + models.append(src_model)
  129 +
  130 + # Return model container
  131 + return models
132 132  
133 133  
134 134 # ======================== #
135 135 # Main routine entry point #
136 136 # ======================== #
137 137 if __name__ == '__main__':
138   - """
139   - Create pull distributions for the energy bands and fluxes given in the
140   - file specified in the argument. The file needs to be created using the
141   - cssens script.
142   - """
143   - # Get input arguments
144   - usage = """
  138 + """
  139 + Create pull distributions for the energy bands and fluxes given in the
  140 + file specified in the argument. The file needs to be created using the
  141 +cssens script.
  142 + """
  143 + # Get input arguments
  144 + usage = """
145 145 make_pull_at_sensitivity_limit filename [max_threads]
146 146  
147 147 Run cssens to create the input file.
148 148 """
149   - if len(sys.argv) < 2 or len(sys.argv) > 3:
150   - print(usage)
151   - sys.exit()
152   -
153   - # Get sensitivity file
154   - filename = sys.argv[1]
155   -
156   - # Set maximum number of threads
157   - if len(sys.argv) > 2:
158   - max_threads = int(sys.argv[2])
159   - else:
160   - max_threads = 1
161   -
162   - # Load sensitivity file
163   - file = open(filename, 'r')
164   - reader = csv.DictReader(file)
165   -
166   - # Loop over rows
167   - for row in reader:
168   -
169   - # Get parameters
170   - loge = float(row['loge'])
171   - if loge < 0:
172   - loge = "m"+str(abs(loge))
173   - else:
174   - loge = "p"+str(abs(loge))
175   - emin = float(row['emin'])
176   - emax = float(row['emax'])
177   - flux = float(row['crab'])
178   -
179   - # Generate model
180   - models = create_model(flux)
181   -
182   - # Processing support?
183   - if has_processing:
184   -
185   - # Wait until one thread has finished
186   - while len(processing.activeChildren()) >= max_threads:
187   - time.sleep(60)
188   -
189   - # Set arguments
190   - args = (loge, emin, emax, models)
191   - kwargs = {}
192   -
193   - # Generate pull distribution
194   - p = processing.Process(target=create_pull, args=args, kwargs=kwargs)
195   - p.start()
196   - print("Process emin=%.4f emax=%.4f started." % (emin, emax))
197   -
198   - # Wait some time
199   - time.sleep(1)
200   -
201   - # ... no
202   - else:
203   - create_pull(loge, emin, emax, models)
204   -
205   - # Processing support
206   - if has_processing:
207   -
208   - # Wait until all threads finished
209   - while len(processing.activeChildren()) > 0:
210   - time.sleep(60)
211   -
212   - # Close file
213   - file.close()
  149 + if len(sys.argv) < 2 or len(sys.argv) > 3:
  150 + print(usage)
  151 + sys.exit()
  152 +
  153 + # Get sensitivity file
  154 + filename = sys.argv[1]
  155 +
  156 + # Set maximum number of threads
  157 + if len(sys.argv) > 2:
  158 + max_threads = int(sys.argv[2])
  159 + else:
  160 + max_threads = 1
  161 +
  162 + # Load sensitivity file
  163 + file = open(filename, 'r')
  164 + reader = csv.DictReader(file)
  165 +
  166 + # Loop over rows
  167 + for row in reader:
  168 +
  169 + # Get parameters
  170 + loge = float(row['loge'])
  171 + if loge < 0:
  172 + loge = "m"+str(abs(loge))
  173 + else:
  174 + loge = "p"+str(abs(loge))
  175 + emin = float(row['emin'])
  176 + emax = float(row['emax'])
  177 + flux = float(row['crab'])
  178 +
  179 + # Generate model
  180 + models = create_model(flux)
  181 +
  182 + # Processing support?
  183 + if has_processing:
  184 +
  185 + # Wait until one thread has finished
  186 + while len(processing.activeChildren()) >= max_threads:
  187 + time.sleep(60)
  188 +
  189 + # Set arguments
  190 + args = (loge, emin, emax, models)
  191 + kwargs = {}
  192 +
  193 + # Generate pull distribution
  194 + p = processing.Process(target=create_pull, args=args, kwargs=kwargs)
  195 + p.start()
  196 + print("Process emin=%.4f emax=%.4f started." % (emin, emax))
  197 +
  198 + # Wait some time
  199 + time.sleep(1)
  200 +
  201 + # ... no
  202 + else:
  203 + create_pull(loge, emin, emax, models)
  204 +
  205 + # Processing support
  206 + if has_processing:
  207 +
  208 + # Wait until all threads finished
  209 + while len(processing.activeChildren()) > 0:
  210 + time.sleep(60)
  211 +
  212 + # Close file
  213 + file.close()
... ...
examples/make_spectrum.py
  View file @5c53c14
... ... @@ -22,10 +22,10 @@ import gammalib
22 22 import ctools
23 23 from cscripts import obsutils
24 24 try:
25   - import matplotlib.pyplot as plt
26   - has_matplotlib = True
  25 + import matplotlib.pyplot as plt
  26 + has_matplotlib = True
27 27 except:
28   - has_matplotlib = False
  28 + has_matplotlib = False
29 29  
30 30  
31 31 # =============== #
... ... @@ -67,7 +67,7 @@ def make_spectrum():
67 67  
68 68 # Generate spectral points
69 69 spectrum = obsutils.spectrum(sim.obs(), "Crab", ebounds)
70   -
  70 +
71 71 # Return spectrum
72 72 return spectrum
73 73  
... ... @@ -100,7 +100,7 @@ def plot_spectrum(spectrum):
100 100  
101 101 # Return
102 102 return
103   -
  103 +
104 104  
105 105 #==========================#
106 106 # Main routine entry point #
... ...
examples/make_survey.py
  View file @5c53c14
... ... @@ -42,33 +42,33 @@ def plot_counts(obs):
42 42 try:
43 43 # Import matplotlib
44 44 import matplotlib.pyplot as plt
45   -
  45 +
46 46 # Set legend fontsize
47 47 params = {'legend.fontsize': 10}
48 48 plt.rcParams.update(params)
49   -
  49 +
50 50 # Set plot styles
51 51 styles = ['b-', 'g-', 'y-', 'n-']
52   -
  52 +
53 53 # Dump header
54 54 print("")
55 55 print("Make plots (using matplotlib):")
56 56 print("==============================")
57   -
  57 +
58 58 # Create figure 1
59 59 plt.figure(1,figsize=(12,6))
60 60 plt.subplots_adjust(hspace=.7)
61   -
  61 +
62 62 # Create subplot 1
63 63 plt.subplot(121)
64 64 plt.title("Spectrum (summed over all pixels)")
65   -
  65 +
66 66 # Loop over observations
67 67 for run in obs:
68   -
  68 +
69 69 # Get event list
70 70 list = run.events()
71   -
  71 +
72 72 # Create energy axis
73 73 ebounds = gammalib.GEbounds()
74 74 emin = gammalib.GEnergy(0.1, "TeV")
... ... @@ -77,7 +77,7 @@ def plot_counts(obs):
77 77 #emax.TeV(100.0)
78 78 ebounds.setlog(emin, emax, 10)
79 79 energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
80   -
  80 +
81 81 # Create spectrum
82 82 print("Extract data:")
83 83 counts = [0.0 for i in range(ebounds.size())]
... ... @@ -85,14 +85,14 @@ def plot_counts(obs):
85 85 for atom in list:
86 86 index = ebounds.index(atom.energy())
87 87 counts[index] = counts[index] + 1.0
88   -
  88 +
89 89 # Create error bars
90 90 error = [math.sqrt(c) for c in counts]
91   -
  91 +
92 92 # Plot spectrum
93 93 plt.loglog(energy, counts, 'ro', label='data')
94 94 #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
95   -
  95 +
96 96 # Extract models
97 97 print("Extract models:")
98 98 sum_model = [0.0 for i in range(ebounds.size())]
... ... @@ -107,26 +107,26 @@ def plot_counts(obs):
107 107 sum_model[i] = sum_model[i] + model[i]
108 108 #plt.loglog(energy, model, styles[k], label=m.name())
109 109 #plt.loglog(energy, sum_model, 'r-', label='total')
110   -
  110 +
111 111 # Put labels
112 112 plt.xlabel("Energy (TeV)")
113 113 plt.ylabel("Counts")
114 114 plt.legend(loc="lower left")
115   -
  115 +
116 116 # Create subplot 2
117 117 plt.subplot(122)
118 118 plt.title("Offset (summed over all energies)")
119   -
  119 +
120 120 # Set Crab direction
121 121 crab = gammalib.GSkyDir()
122 122 crab.radec_deg(83.63, 22.01)
123   -
  123 +
124 124 # Loop over observations
125 125 for run in obs:
126   -
  126 +
127 127 # Get event list
128 128 list = run.events()
129   -
  129 +
130 130 # Create offset histogram
131 131 print("Extract data:")
132 132 nx = 30
... ... @@ -139,13 +139,13 @@ def plot_counts(obs):
139 139 index = int(off2/doffset2)
140 140 if index < nx:
141 141 counts[index] = counts[index] + 1.0
142   -
  142 +
143 143 # Create error bars
144 144 error = [math.sqrt(c) for c in counts]
145   -
  145 +
146 146 # Plot distribution
147 147 #plt.semilogy(offset2, counts, 'ro', label='data')
148   -
  148 +
149 149 # Extract models
150 150 print("Extract models:")
151 151 sum_model = [0.0 for i in range(nx)]
... ... @@ -164,15 +164,15 @@ def plot_counts(obs):
164 164 #plt.plot(offset2, model, styles[k], label=m.name())
165 165 #plt.plot(offset2, sum_model, 'r-', label='total')
166 166 #plt.ylim(ymin=0.1)
167   -
  167 +
168 168 # Put labels
169 169 plt.xlabel("Offset (deg^2)")
170 170 plt.ylabel("Counts")
171 171 #plt.legend(loc="upper right")
172   -
  172 +
173 173 # Show counts spectra
174 174 plt.show()
175   -
  175 +
176 176 except ImportError:
177 177 print("Matplotlib is not (correctly) installed on your system.")
178 178  
... ... @@ -186,26 +186,26 @@ def plot_counts(obs):
186 186 def add_background_model(obs):
187 187 """
188 188 Add standard CTA background model to observations container.
189   -
  189 +
190 190 We use a simple power law here, scaled to Konrad's E configuration
191 191 performance table. The model needs still to be validated.
192 192 """
193 193 # Recover models from observation
194 194 models = obs.models()
195   -
  195 +
196 196 # Define background model
197 197 bgd_radial = gammalib.GCTAModelRadialGauss(3.0)
198 198 bgd_spectrum = gammalib.GModelSpectralPlaw(61.8e-6, -1.85, gammalib.GEnergy(1.0, "TeV"))
199 199 bgd_model = gammalib.GCTAModelRadialAcceptance(bgd_radial, bgd_spectrum)
200 200 bgd_model.name("Background")
201 201 bgd_model.instruments("CTA")
202   -
  202 +
203 203 # Add background model to container
204 204 models.append(bgd_model)
205   -
  205 +
206 206 # Put container back in observation container
207 207 obs.models(models)
208   -
  208 +
209 209 # Return observation container
210 210 return obs
211 211  
... ... @@ -220,7 +220,7 @@ def crab_spec():
220 220 """
221 221 # Set parameters
222 222 spectrum = gammalib.GModelSpectralPlaw(5.7e-16, -2.48, gammalib.GEnergy(0.3, "TeV"))
223   -
  223 +
224 224 # Return spectrum
225 225 return spectrum
226 226  
... ... @@ -234,13 +234,13 @@ def survey_single():
234 234 """
235 235 # Allocate observation container
236 236 obs = gammalib.GObservations()
237   -
  237 +
238 238 # Set single pointing at galactic centre
239 239 pntdir = gammalib.GSkyDir()
240 240 pntdir.lb_deg(0.0, 0.0)
241 241 run = obsutils.set_obs(pntdir)
242 242 obs.append(run)
243   -
  243 +
244 244 # Define single point source with Crab flux at galactic centre
245 245 center = gammalib.GSkyDir()
246 246 center.lb_deg(0.0, 0.0)
... ... @@ -248,12 +248,12 @@ def survey_single():
248 248 point_spectrum = crab_spec()
249 249 point = gammalib.GModelSky(point_spatial, point_spectrum)
250 250 point.name('GC source')
251   -
  251 +
252 252 # Create model container
253 253 models = gammalib.GModels()
254 254 models.append(point)
255 255 obs.models(models)
256   -
  256 +
257 257 # Return observation container
258 258 return obs
259 259  
... ... @@ -264,24 +264,24 @@ def survey_single():
264 264 def survey_gplane(lrange=10, lstep=2):
265 265 """
266 266 Creates a single observation survey for test purposes.
267   -
  267 +
268 268 Keywords:
269 269 lrange - Longitude range (integer deg)
270 270 lstep - Longitude step size (integer deg)
271 271 """
272 272 # Allocate observation container
273 273 obs = gammalib.GObservations()
274   -
  274 +
275 275 # Loop over longitudes
276 276 for l in range(-lrange,lrange+lstep,lstep):
277   -
  277 +
278 278 # Set pointing
279 279 pntdir = gammalib.GSkyDir()
280 280 pntdir.lb_deg(l, 0.0)
281 281 run = obsutils.set_obs(pntdir)
282 282 run.id(str(l))
283 283 obs.append(run)
284   -
  284 +
285 285 # Define single point source with Crab flux at galactic centre
286 286 center = gammalib.GSkyDir()
287 287 center.lb_deg(0.0, 0.0)
... ... @@ -289,12 +289,12 @@ def survey_gplane(lrange=10, lstep=2):
289 289 point_spectrum = crab_spec()
290 290 point = gammalib.GModelSky(point_spatial, point_spectrum)
291 291 point.name('GC source')
292   -
  292 +
293 293 # Create model container
294 294 models = gammalib.GModels()
295 295 models.append(point)
296 296 obs.models(models)
297   -
  297 +
298 298 # Return observation container
299 299 return obs
300 300  
... ... @@ -308,7 +308,7 @@ if __name__ == &#39;__main__&#39;:
308 308 """
309 309 # Initialise flags
310 310 need_help = False
311   -
  311 +
312 312 # Test for command line arguments
313 313 print(sys.argv[0])
314 314 if (len(sys.argv) > 1):
... ... @@ -322,42 +322,42 @@ if __name__ == &#39;__main__&#39;:
322 322 print("Usage: example_survey.py [OPTIONS]")
323 323 print(" -h Display this usage message")
324 324 sys.exit()
325   -
  325 +
326 326 # Dump header
327 327 print("***********************")
328 328 print("* Simulate CTA survey *")
329 329 print("***********************")
330   -
  330 +
331 331 # Remove any existing result files
332 332 list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
333 333 for files in list:
334 334 for file in files:
335 335 os.remove(file)
336   -
  336 +
337 337 # Setup single observation survey
338 338 #obs = survey_single()
339 339 obs = survey_gplane()
340   -
  340 +
341 341 # Add background model
342 342 obs = add_background_model(obs)
343   -
  343 +
344 344 # Simulate events
345 345 print("Simulate events")
346 346 obs = obsutils.sim(obs)
347   -
  347 +
348 348 # Make counts map
349 349 print("Make counts map")
350 350 cntmap = obsutils.cntmap(obs)
351   -
  351 +
352 352 # Fit observations
353 353 print("Fit observations")
354 354 like = obsutils.fit(obs)
355 355 #print like.opt()
356 356 #print like.obs().models()
357   -
  357 +
358 358 # Make model map
359 359 print("Make model map (this step will take some time)")
360 360 modmap = obsutils.modmap(obs)
361   -
  361 +
362 362 # Show fit results
363 363 #plot_counts(like.obs())
... ...
examples/make_ts_distributions.py
  View file @5c53c14
... ... @@ -35,155 +35,155 @@ from cscripts import cstsdist
35 35 # Try importing processing module
36 36 has_processing = False
37 37 try:
38   - import processing
39   - has_processing = True
  38 + import processing
  39 + has_processing = True
40 40 except:
41   - pass
  41 + pass
42 42  
43 43  
44 44 # ====================== #
45 45 # Create TS distribution #
46 46 # ====================== #
47 47 def create_ts(loge, emin, emax, ntrials=100, duration=180000.0,
48   - enumbins=0, log=False):
49   - """
50   - Create TS distribution.
51   -
52   - Parameters:
53   - loge - Logarithm of mean energy in TeV
54   - emin - Minimum energy (TeV)
55   - emax - Maximum energy (TeV)
56   - Keywords:
57   - log - Create log file(s)
58   - """
59   - # Generate output filename
60   - outfile = "ts_"+str(loge)+".dat"
61   -
62   - # Setup cstsdist tool
63   - tsdist = cstsdist()
64   - tsdist["inmodel"] = "$CTOOLS/share/models/crab.xml"
65   - tsdist["srcname"] = "Crab"
66   - tsdist["outfile"] = outfile
67   - tsdist["ntrials"] = ntrials
68   - tsdist["caldb"] = "prod2"
69   - tsdist["irf"] = "South_50h"
70   - tsdist["ra"] = 83.63
71   - tsdist["dec"] = 22.01
72   - tsdist["emin"] = emin
73   - tsdist["emax"] = emax
74   - tsdist["enumbins"] = enumbins
75   - tsdist["tmin"] = 0
76   - tsdist["tmax"] = duration
77   - tsdist["rad"] = 5.0
78   - tsdist["npix"] = 200
79   - tsdist["binsz"] = 0.05
80   - #tsdist["debug"].boolean(True)
81   -
82   - # Optionally open the log file
83   - if log:
84   - tsdist.logFileOpen()
85   -
86   - # Run tool
87   - tsdist.run()
88   -
89   - # Return
90   - return
  48 + enumbins=0, log=False):
  49 + """
  50 + Create TS distribution.
  51 +
  52 + Parameters:
  53 + loge - Logarithm of mean energy in TeV
  54 + emin - Minimum energy (TeV)
  55 + emax - Maximum energy (TeV)
  56 + Keywords:
  57 + log - Create log file(s)
  58 + """
  59 + # Generate output filename
  60 + outfile = "ts_"+str(loge)+".dat"
  61 +
  62 + # Setup cstsdist tool
  63 + tsdist = cstsdist()
  64 + tsdist["inmodel"] = "$CTOOLS/share/models/crab.xml"
  65 + tsdist["srcname"] = "Crab"
  66 + tsdist["outfile"] = outfile
  67 + tsdist["ntrials"] = ntrials
  68 + tsdist["caldb"] = "prod2"
  69 + tsdist["irf"] = "South_50h"
  70 + tsdist["ra"] = 83.63
  71 + tsdist["dec"] = 22.01
  72 + tsdist["emin"] = emin
  73 + tsdist["emax"] = emax
  74 + tsdist["enumbins"] = enumbins
  75 + tsdist["tmin"] = 0
  76 + tsdist["tmax"] = duration
  77 + tsdist["rad"] = 5.0
  78 + tsdist["npix"] = 200
  79 + tsdist["binsz"] = 0.05
  80 + #tsdist["debug"].boolean(True)
  81 +
  82 + # Optionally open the log file
  83 + if log:
  84 + tsdist.logFileOpen()
  85 +
  86 + # Run tool
  87 + tsdist.run()
  88 +
  89 + # Return
  90 + return
91 91  
92 92  
93 93 # ======================== #
94 94 # Main routine entry point #
95 95 # ======================== #
96 96 if __name__ == '__main__':
97   - """
98   - Create TS distribution in a number of energy bands.
99   - """
100   - # Get input arguments
101   - usage = "make_ts_distributions [-n ntrials] [-e enumbins] [-m max_threads]"
102   - if len(sys.argv) < 1:
103   - print(usage)
104   - sys.exit()
105   -
106   - # Set default parameters
107   - ntrials = 100
108   - enumbins = 0
109   - duration = 180000.0
110   - max_threads = 1
111   -
112   - # Parameter dictionnary
113   - pars = [{'option': '-n', 'value': ntrials},
114   - {'option': '-e', 'value': enumbins},
115   - {'option': '-d', 'value': duration},
116   - {'option': '-m', 'value': max_threads}]
117   -
118   - # Gather parameters from command line
119   - i = 1
120   - while i < len(sys.argv):
121   -
122   - # Search for options
123   - for par in pars:
124   - if sys.argv[i] == par['option']:
125   - if len(sys.argv) > i+1:
126   - i += 1
127   - try:
128   - par['value'] = int(sys.argv[i])
129   - except:
130   - print(usage)
131   - sys.exit()
132   - else:
133   - print(usage)
134   - sys.exit()
135   -
136   - # Next item
137   - i += 1
138   -
139   - # Recover parameters
140   - ntrials = pars[0]['value']
141   - enumbins = pars[1]['value']
142   - duration = pars[2]['value']
143   - max_threads = pars[3]['value']
144   -
145   - # Loop over energy bands. The energy bands are those that are also
146   - # used for sensitivity computation.
147   - for ieng in range(20):
148   -
149   - # Set energies
150   - loge = -1.7 + ieng * 0.2
151   - emean = pow(10.0, loge)
152   - emin = pow(10.0, loge-0.1)
153   - emax = pow(10.0, loge+0.1)
154   - if loge < 0:
155   - loge = "m"+str(abs(loge))
156   - else:
157   - loge = "p"+str(abs(loge))
158   -
159   - # Processing support?
160   - if has_processing:
161   -
162   - # Wait until one thread has finished
163   - while len(processing.activeChildren()) >= max_threads:
164   - time.sleep(10)
165   -
166   - # Set arguments
167   - args = (loge, emin, emax)
168   - kwargs = {'ntrials': ntrials, 'enumbins': enumbins,
169   - 'duration': duration}
170   -
171   - # Generate pull distribution
172   - p = processing.Process(target=create_ts, args=args, kwargs=kwargs)
173   - p.start()
174   - print("Process emin=%.4f emax=%.4f started." % (emin, emax))
175   -
176   - # Wait a short time to allow process to start
177   - time.sleep(1)
178   -
179   - # ... no
180   - else:
181   - create_ts(loge, emin, emax, ntrials=ntrials, enumbins=enumbins,
182   - duration=duration)
183   -
184   - # Processing support
185   - if has_processing:
186   -
187   - # Wait until all threads finished
188   - while len(processing.activeChildren()) > 0:
189   - time.sleep(10)
  97 + """
  98 + Create TS distribution in a number of energy bands.
  99 + """
  100 + # Get input arguments
  101 + usage = "make_ts_distributions [-n ntrials] [-e enumbins] [-m max_threads]"
  102 + if len(sys.argv) < 1:
  103 + print(usage)
  104 + sys.exit()
  105 +
  106 + # Set default parameters
  107 + ntrials = 100
  108 + enumbins = 0
  109 + duration = 180000.0
  110 + max_threads = 1
  111 +
  112 + # Parameter dictionnary
  113 + pars = [{'option': '-n', 'value': ntrials},
  114 + {'option': '-e', 'value': enumbins},
  115 + {'option': '-d', 'value': duration},
  116 + {'option': '-m', 'value': max_threads}]
  117 +
  118 + # Gather parameters from command line
  119 + i = 1
  120 + while i < len(sys.argv):
  121 +
  122 + # Search for options
  123 + for par in pars:
  124 + if sys.argv[i] == par['option']:
  125 + if len(sys.argv) > i+1:
  126 + i += 1
  127 + try:
  128 + par['value'] = int(sys.argv[i])
  129 + except:
  130 + print(usage)
  131 + sys.exit()
  132 + else:
  133 + print(usage)
  134 + sys.exit()
  135 +
  136 + # Next item
  137 + i += 1
  138 +
  139 + # Recover parameters
  140 + ntrials = pars[0]['value']
  141 + enumbins = pars[1]['value']
  142 + duration = pars[2]['value']
  143 + max_threads = pars[3]['value']
  144 +
  145 + # Loop over energy bands. The energy bands are those that are also
  146 + # used for sensitivity computation.
  147 + for ieng in range(20):
  148 +
  149 + # Set energies
  150 + loge = -1.7 + ieng * 0.2
  151 + emean = pow(10.0, loge)
  152 + emin = pow(10.0, loge-0.1)
  153 + emax = pow(10.0, loge+0.1)
  154 + if loge < 0:
  155 + loge = "m"+str(abs(loge))
  156 + else:
  157 + loge = "p"+str(abs(loge))
  158 +
  159 + # Processing support?
  160 + if has_processing:
  161 +
  162 + # Wait until one thread has finished
  163 + while len(processing.activeChildren()) >= max_threads:
  164 + time.sleep(10)
  165 +
  166 + # Set arguments
  167 + args = (loge, emin, emax)
  168 + kwargs = {'ntrials': ntrials, 'enumbins': enumbins,
  169 + 'duration': duration}
  170 +
  171 + # Generate pull distribution
  172 + p = processing.Process(target=create_ts, args=args, kwargs=kwargs)
  173 + p.start()
  174 + print("Process emin=%.4f emax=%.4f started." % (emin, emax))
  175 +
  176 + # Wait a short time to allow process to start
  177 + time.sleep(1)
  178 +
  179 + # ... no
  180 + else:
  181 + create_ts(loge, emin, emax, ntrials=ntrials, enumbins=enumbins,
  182 + duration=duration)
  183 +
  184 + # Processing support
  185 + if has_processing:
  186 +
  187 + # Wait until all threads finished
  188 + while len(processing.activeChildren()) > 0:
  189 + time.sleep(10)
... ...
examples/pipeline_binned_disk.py
  View file @5c53c14
... ... @@ -39,7 +39,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
39 39 caldb="prod2", irf="South_50h"):
40 40 """
41 41 Returns an observation container.
42   -
  42 +
43 43 Keywords:
44 44 pattern - Pointing pattern, either "single" or "four"
45 45 ra - RA of pattern centre [deg] (default: 83.6331)
... ... @@ -58,7 +58,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
58 58 ra=ra,
59 59 dec=dec,
60 60 offset=offset)
61   -
  61 +
62 62 # Get observation container
63 63 obs = obsutils.set_obs_list(obs_def_list,
64 64 duration=duration,
... ... @@ -78,13 +78,13 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
78 78 def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
79 79 """
80 80 Setup model for analysis.
81   -
  81 +
82 82 Keywords:
83 83 model - Model Xml file
84 84 """
85 85 # Append model from file to observation container
86 86 obs.models(gammalib.GModels(model))
87   -
  87 +
88 88 # Return observation container
89 89 return obs
90 90  
... ... @@ -170,7 +170,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0,
170 170 like["irf"] = irf
171 171 like["debug"] = True # Switch this always on for results in console
172 172 like.execute()
173   -
  173 +
174 174 # Return
175 175 return
176 176  
... ...
examples/pipeline_binned_mem.py
  View file @5c53c14
... ... @@ -39,7 +39,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
39 39 caldb="prod2", irf="South_50h"):
40 40 """
41 41 Returns an observation container.
42   -
  42 +
43 43 Keywords:
44 44 pattern - Pointing pattern, either "single" or "four"
45 45 ra - RA of pattern centre [deg] (default: 83.6331)
... ... @@ -58,7 +58,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
58 58 ra=ra,
59 59 dec=dec,
60 60 offset=offset)
61   -
  61 +
62 62 # Get observation container
63 63 obs = obsutils.set_obs_list(obs_def_list,
64 64 duration=duration,
... ... @@ -78,13 +78,13 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
78 78 def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
79 79 """
80 80 Setup model for analysis.
81   -
  81 +
82 82 Keywords:
83 83 model - Model Xml file
84 84 """
85 85 # Append model from file to observation container
86 86 obs.models(gammalib.GModels(model))
87   -
  87 +
88 88 # Return observation container
89 89 return obs
90 90  
... ... @@ -144,7 +144,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0,
144 144 like = ctools.ctlike(obs)
145 145 like["debug"] = True # Switch this always on for results in console
146 146 like.run()
147   -
  147 +
148 148 # Return
149 149 return
150 150  
... ...
examples/pipeline_stacked_disk.py
  View file @5c53c14
... ... @@ -39,7 +39,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
39 39 caldb="prod2", irf="South_50h"):
40 40 """
41 41 Returns an observation container.
42   -
  42 +
43 43 Keywords:
44 44 pattern - Pointing pattern, either "single" or "four"
45 45 ra - RA of pattern centre [deg] (default: 83.6331)
... ... @@ -58,7 +58,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
58 58 ra=ra,
59 59 dec=dec,
60 60 offset=offset)
61   -
  61 +
62 62 # Get observation container
63 63 obs = obsutils.set_obs_list(obs_def_list,
64 64 duration=duration,
... ... @@ -78,13 +78,13 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
78 78 def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
79 79 """
80 80 Setup model for analysis.
81   -
  81 +
82 82 Keywords:
83 83 model - Model Xml file
84 84 """
85 85 # Append model from file to observation container
86 86 obs.models(gammalib.GModels(model))
87   -
  87 +
88 88 # Return observation container
89 89 return obs
90 90  
... ... @@ -215,7 +215,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0,
215 215 like["irf"] = irf
216 216 like["debug"] = True # Switch this always on for results in console
217 217 like.execute()
218   -
  218 +
219 219 # Return
220 220 return
221 221  
... ...
examples/pipeline_stacked_mem.py
  View file @5c53c14
... ... @@ -39,7 +39,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
39 39 caldb="prod2", irf="South_50h"):
40 40 """
41 41 Returns an observation container.
42   -
  42 +
43 43 Keywords:
44 44 pattern - Pointing pattern, either "single" or "four"
45 45 ra - RA of pattern centre [deg] (default: 83.6331)
... ... @@ -58,7 +58,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
58 58 ra=ra,
59 59 dec=dec,
60 60 offset=offset)
61   -
  61 +
62 62 # Get observation container
63 63 obs = obsutils.set_obs_list(obs_def_list,
64 64 duration=duration,
... ... @@ -78,13 +78,13 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
78 78 def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
79 79 """
80 80 Setup model for analysis.
81   -
  81 +
82 82 Keywords:
83 83 model - Model Xml file
84 84 """
85 85 # Append model from file to observation container
86 86 obs.models(gammalib.GModels(model))
87   -
  87 +
88 88 # Return observation container
89 89 return obs
90 90  
... ... @@ -194,7 +194,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0,
194 194 like = ctools.ctlike(bin.obs())
195 195 like["debug"] = True # Switch this always on for results in console
196 196 like.run()
197   -
  197 +
198 198 # Return
199 199 return
200 200  
... ...
examples/pipeline_unbinned_disk.py
  View file @5c53c14
... ... @@ -39,7 +39,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
39 39 caldb="prod2", irf="South_50h"):
40 40 """
41 41 Returns an observation container.
42   -
  42 +
43 43 Keywords:
44 44 pattern - Pointing pattern, either "single" or "four"
45 45 ra - RA of pattern centre [deg] (default: 83.6331)
... ... @@ -58,7 +58,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
58 58 ra=ra,
59 59 dec=dec,
60 60 offset=offset)
61   -
  61 +
62 62 # Get observation container
63 63 obs = obsutils.set_obs_list(obs_def_list,
64 64 duration=duration,
... ... @@ -78,13 +78,13 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
78 78 def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
79 79 """
80 80 Setup model for analysis.
81   -
  81 +
82 82 Keywords:
83 83 model - Model Xml file
84 84 """
85 85 # Append model from file to observation container
86 86 obs.models(gammalib.GModels(model))
87   -
  87 +
88 88 # Return observation container
89 89 return obs
90 90  
... ... @@ -145,7 +145,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0,
145 145 like["irf"] = irf
146 146 like["debug"] = True # Switch this always on for results in console
147 147 like.execute()
148   -
  148 +
149 149 # Return
150 150 return
151 151  
... ...
examples/pipeline_unbinned_mem.py
  View file @5c53c14
... ... @@ -39,7 +39,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
39 39 caldb="prod2", irf="South_50h"):
40 40 """
41 41 Returns an observation container.
42   -
  42 +
43 43 Keywords:
44 44 pattern - Pointing pattern, either "single" or "four"
45 45 ra - RA of pattern centre [deg] (default: 83.6331)
... ... @@ -58,7 +58,7 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
58 58 ra=ra,
59 59 dec=dec,
60 60 offset=offset)
61   -
  61 +
62 62 # Get observation container
63 63 obs = obsutils.set_obs_list(obs_def_list,
64 64 duration=duration,
... ... @@ -78,13 +78,13 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
78 78 def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
79 79 """
80 80 Setup model for analysis.
81   -
  81 +
82 82 Keywords:
83 83 model - Model Xml file
84 84 """
85 85 # Append model from file to observation container
86 86 obs.models(gammalib.GModels(model))
87   -
  87 +
88 88 # Return observation container
89 89 return obs
90 90  
... ... @@ -130,7 +130,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0,
130 130 like = ctools.ctlike(select.obs())
131 131 like["debug"] = True # Switch this always on for results in console
132 132 like.run()
133   -
  133 +
134 134 # Return
135 135 return
136 136  
... ...
examples/show_butterfly.py
  View file @5c53c14
... ... @@ -25,46 +25,46 @@ try:
25 25 import matplotlib.pyplot as plt
26 26 except:
27 27 sys.exit("This script needs matplotlib")
28   -
  28 +
29 29 # =============== #
30 30 # Script entry #
31 31 # =============== #
32 32 if __name__ == "__main__":
33   -
  33 +
34 34 # check for given butterfly file
35 35 if not len(sys.argv) == 2:
36 36 sys.exit("Usage: show_butterfly.py butterfly.txt")
37   -
  37 +
38 38 # Read given butterfly file
39 39 filename = sys.argv[1]
40 40 csv = gammalib.GCsv(filename)
41   -
  41 +
42 42 # initialise arrays to be filled
43 43 butterfly_x = []
44 44 butterfly_y = []
45   -
  45 +
46 46 line_x = []
47 47 line_y = []
48   -
  48 +
49 49 # loop over rows of the file
50 50 nrows = csv.nrows()
51 51 for row in range(nrows):
52 52 # Compute upper edge of confidence band
53 53 butterfly_x.append(csv.real(row,0))
54 54 butterfly_y.append(csv.real(row,1)+csv.real(row,2))
55   -
  55 +
56 56 # Set line values
57 57 line_x.append(csv.real(row,0))
58 58 line_y.append(csv.real(row,1))
59   -
  59 +
60 60 # Loop over the rows backwards to compute the lower edge
61 61 # of the confidence band
62 62 for row in range(nrows):
63 63 index = nrows - 1 - row
64   -
  64 +
65 65 butterfly_x.append(csv.real(index,0))
66 66 butterfly_y.append(csv.real(index,1)-csv.real(index,2))
67   -
  67 +
68 68 # plot the butterfly and spectral line
69 69 plt.figure()
70 70 plt.loglog()
... ... @@ -73,4 +73,4 @@ if __name__ == &quot;__main__&quot;:
73 73 plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5)
74 74 plt.xlabel("Energy [MeV]")
75 75 plt.ylabel(r"dN/dE [MeV$^{-1}$ s$^{-1}$ cm$^{-2}$]")
76   - plt.show()
77 76 \ No newline at end of file
  77 + plt.show()
... ...
examples/show_model.py
  View file @5c53c14
... ... @@ -35,7 +35,7 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100):
35 35 """
36 36 # Get spectral component
37 37 spectrum = model.spectral()
38   -
  38 +
39 39 # Setup energy axis
40 40 e_min = gammalib.GEnergy(emin, "TeV")
41 41 e_max = gammalib.GEnergy(emax, "TeV")
... ... @@ -61,7 +61,7 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100):
61 61 if min < 0.0001*max:
62 62 min = 0.0001*max
63 63 min = math.pow(10.0, math.log10(min)-1)
64   -
  64 +
65 65  
66 66 # Show spectrum
67 67 plt.figure()
... ... @@ -73,11 +73,11 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100):
73 73 plt.ylabel(r"dN/dE [ph s$^{-1}$ cm$^{-2}$ MeV$^{-1}$]")
74 74 plt.ylim([min,max])
75 75 plt.show()
76   -
  76 +
77 77 # Return
78 78 return
79   -
80   -
  79 +
  80 +
81 81 # ============= #
82 82 # Script entry #
83 83 # ============= #
... ... @@ -100,6 +100,6 @@ if __name__ == &quot;__main__&quot;:
100 100  
101 101 # Extract relevant model
102 102 model = models[index]
103   -
  103 +
104 104 # Plot spectrum
105 105 plot_spectrum(model)
... ...
examples/show_pha.py
  View file @5c53c14
... ... @@ -53,7 +53,7 @@ def show_pha(pha):
53 53 for i in range(pha.size()):
54 54 counts[i] = pha[i]
55 55 error = [math.sqrt(c) for c in counts]
56   -
  56 +
57 57 # Plot data
58 58 if channels:
59 59 plt.semilogy(energy, counts, 'ro')
... ... @@ -70,11 +70,11 @@ def show_pha(pha):
70 70  
71 71 # Show plot
72 72 plt.show()
73   -
  73 +
74 74 # Return
75 75 #return
76   -
77   -
  76 +
  77 +
78 78 # ======================== #
79 79 # Main routine entry point #
80 80 # ======================== #
... ... @@ -107,4 +107,3 @@ if __name__ == &#39;__main__&#39;:
107 107 # Show PHA spectrum
108 108 if has_matplotlib:
109 109 show_pha(pha)
110   -
111 110 \ No newline at end of file
... ...
examples/show_pull_histogram.py
  View file @5c53c14
... ... @@ -39,85 +39,84 @@ import csv
39 39 # Read pull distribution #
40 40 # ====================== #
41 41 def read_pull(filename, parname):
42   - """
43   - Read pull distribution.
44   - """
45   - # Initialise list
46   - values = []
47   -
48   - # Open reader
49   - reader = csv.reader(open(filename, 'r'), delimiter=',')
50   -
51   - # Read rows
52   - first = True
53   - index = -1
54   - for row in reader:
55   -
56   - # Get column index if first row
57   - if first:
58   - try:
59   - index = row.index(parname)
60   - except:
61   - sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n')
62   - for p in row:
63   - sys.stdout.write(' "'+p+'"\n')
64   - raise NameError(parname)
65   -
66   - # Handle data rows
67   - else:
68   - values.append(float(row[index]))
69   -
70   - # Flag that first row has been passed
71   - first = False
72   -
73   - # Create numpy array
74   - a = np.array(values)
75   -
76   - # Return array
77   - return a
  42 + """
  43 + Read pull distribution.
  44 + """
  45 + # Initialise list
  46 + values = []
  47 +
  48 + # Open reader
  49 + reader = csv.reader(open(filename, 'r'), delimiter=',')
  50 +
  51 + # Read rows
  52 + first = True
  53 + index = -1
  54 + for row in reader:
  55 +
  56 + # Get column index if first row
  57 + if first:
  58 + try:
  59 + index = row.index(parname)
  60 + except:
  61 + sys.stdout.write('ERROR: Parameter "'+parname+'" not found in list:\n')
  62 + for p in row:
  63 + sys.stdout.write(' "'+p+'"\n')
  64 + raise NameError(parname)
  65 +
  66 + # Handle data rows
  67 + else:
  68 + values.append(float(row[index]))
  69 +
  70 + # Flag that first row has been passed
  71 + first = False
  72 +
  73 + # Create numpy array
  74 + a = np.array(values)
  75 +
  76 + # Return array
  77 + return a
78 78  
79 79  
80 80 # ======================== #
81 81 # Main routine entry point #
82 82 # ======================== #
83 83 if __name__ == '__main__':
84   - """
85   - Display pull histogram generated by cspull using matplotlib.
86   - """
87   - # Print usage information
88   - usage = "Usage: show_pull_histogram filename parname [bins]"
89   - if len(sys.argv) < 3:
90   - sys.stdout.write(usage+"\n")
91   - sys.exit()
92   -
93   - # Extract parameters
94   - filename = sys.argv[1]
95   - parname = sys.argv[2]
96   - if len(sys.argv) > 3:
97   - nbins = int(sys.argv[3])
98   - else:
99   - nbins = 50
100   -
101   - # Read values from CSV file
102   - values = read_pull(filename, parname)
103   -
104   - # Create histogram
105   - n, bins, patches = plt.hist(values, nbins, range=[-4.0,4.0],
106   - normed=True, facecolor='green')
107   -
108   - # Create expected distribution
109   - y = mlab.normpdf(bins, 0.0, 1.0)
110   - l = plt.plot(bins, y, 'r-', linewidth=2)
111   -
112   - # Set parname for plotting
113   - name = parname
114   -
115   - # Set plot
116   - plt.xlabel('Pull ('+name+')')
117   - plt.ylabel('Arbitrary units')
118   - plt.title(name)
119   - plt.grid(True)
120   -
121   - # Show histogram
122   - plt.show()
123   -
124 84 \ No newline at end of file
  85 + """
  86 + Display pull histogram generated by cspull using matplotlib.
  87 + """
  88 + # Print usage information
  89 + usage = "Usage: show_pull_histogram filename parname [bins]"
  90 + if len(sys.argv) < 3:
  91 + sys.stdout.write(usage+"\n")
  92 + sys.exit()
  93 +
  94 + # Extract parameters
  95 + filename = sys.argv[1]
  96 + parname = sys.argv[2]
  97 + if len(sys.argv) > 3:
  98 + nbins = int(sys.argv[3])
  99 + else:
  100 + nbins = 50
  101 +
  102 + # Read values from CSV file
  103 + values = read_pull(filename, parname)
  104 +
  105 + # Create histogram
  106 + n, bins, patches = plt.hist(values, nbins, range=[-4.0,4.0],
  107 + normed=True, facecolor='green')
  108 +
  109 + # Create expected distribution
  110 + y = mlab.normpdf(bins, 0.0, 1.0)
  111 + l = plt.plot(bins, y, 'r-', linewidth=2)
  112 +
  113 + # Set parname for plotting
  114 + name = parname
  115 +
  116 + # Set plot
  117 + plt.xlabel('Pull ('+name+')')
  118 + plt.ylabel('Arbitrary units')
  119 + plt.title(name)
  120 + plt.grid(True)
  121 +
  122 + # Show histogram
  123 + plt.show()
... ...
examples/show_response.py
  View file @5c53c14
... ... @@ -32,58 +32,58 @@ import matplotlib.pyplot as plt
32 32 # Compute Eq. (17) of Li & Ma (1983) #
33 33 # ================================== #
34 34 def sigma_lima(Non, Noff, alpha=0.2):
35   - """
36   - Compute Eq. (17) of Li & Ma (1983).
37   -
38   - Parameters:
39   - Non - Number of on counts
40   - Noff - Number of off counts
41   - Keywords:
42   - alpha - Ratio of on-to-off exposure
43   - """
44   - # Compute sensitivity
45   - alpha1 = alpha + 1.0
46   - sum = Non + Noff
47   - arg1 = Non / sum
48   - arg2 = Noff / sum
49   - term1 = Non * math.log((alpha1/alpha)*arg1)
50   - term2 = Noff * math.log(alpha1*arg2)
51   - sigma = math.sqrt(2.0 * (term1 + term2))
52   -
53   - # Return sensitivity
54   - return sigma
  35 + """
  36 + Compute Eq. (17) of Li & Ma (1983).
  37 +
  38 + Parameters:
  39 + Non - Number of on counts
  40 + Noff - Number of off counts
  41 + Keywords:
  42 + alpha - Ratio of on-to-off exposure
  43 + """
  44 + # Compute sensitivity
  45 + alpha1 = alpha + 1.0
  46 + sum = Non + Noff
  47 + arg1 = Non / sum
  48 + arg2 = Noff / sum
  49 + term1 = Non * math.log((alpha1/alpha)*arg1)
  50 + term2 = Noff * math.log(alpha1*arg2)
  51 + sigma = math.sqrt(2.0 * (term1 + term2))
  52 +
  53 + # Return sensitivity
  54 + return sigma
55 55  
56 56  
57 57 # ======================================== #
58 58 # Solve Eq. (17) of Li & Ma (1983) for Non #
59 59 # ======================================== #
60 60 def Non_lima(sigma, Noff, alpha=0.2):
61   - """
62   - Solve Eq. (17) of Li & Ma (1983) for Non
63   -
64   - Parameters:
65   - sigma - Required significance
66   - Noff - Number of off counts
67   - Keywords:
68   - alpha - Ratio of on-to-off exposure
69   - """
70   - # Set initial guess for Non
71   - S = 25.0
72   - Non = S + alpha*Noff
73   -
74   - # Compute initial sigma estimate
75   - s = sigma_lima(Non, Noff, alpha=alpha)
76   - delta = s - sigma
77   -
78   - # Iterate
79   - while abs(delta) > 1.0e-6:
80   - S *= sigma/s
81   - Non = S + alpha*Noff
82   - s = sigma_lima(Non, Noff, alpha=alpha)
83   - delta = s - sigma
84   -
85   - # Return
86   - return Non
  61 + """
  62 + Solve Eq. (17) of Li & Ma (1983) for Non
  63 +
  64 + Parameters:
  65 + sigma - Required significance
  66 + Noff - Number of off counts
  67 + Keywords:
  68 + alpha - Ratio of on-to-off exposure
  69 + """
  70 + # Set initial guess for Non
  71 + S = 25.0
  72 + Non = S + alpha*Noff
  73 +
  74 + # Compute initial sigma estimate
  75 + s = sigma_lima(Non, Noff, alpha=alpha)
  76 + delta = s - sigma
  77 +
  78 + # Iterate
  79 + while abs(delta) > 1.0e-6:
  80 + S *= sigma/s
  81 + Non = S + alpha*Noff
  82 + s = sigma_lima(Non, Noff, alpha=alpha)
  83 + delta = s - sigma
  84 +
  85 + # Return
  86 + return Non
87 87  
88 88  
89 89 # ======================= #
... ... @@ -95,7 +95,7 @@ def show_one_effective_area(rsp, name, color=&quot;r&quot;):
95 95 """
96 96 # Set to figure 1
97 97 plt.figure(1)
98   -
  98 +
99 99 # Generate logE(TeV) vector
100 100 nbins = 20
101 101 logE = [-1.7+0.2*i for i in range(nbins)]
... ... @@ -130,7 +130,7 @@ def show_one_background_rate(rsp, name, color=&quot;r&quot;):
130 130 """
131 131 # Set to figure 2
132 132 plt.figure(2)
133   -
  133 +
134 134 # Generate logE(TeV) vector
135 135 nbins = 20
136 136 logE = [-1.7+0.2*i for i in range(nbins)]
... ... @@ -142,7 +142,7 @@ def show_one_background_rate(rsp, name, color=&quot;r&quot;):
142 142 bgrate[i] = rsp.background()(logE[i], 0.0, 0.0)
143 143 if bgrate[i] < 1.0e-10:
144 144 bgrate[i] = 0.0
145   -
  145 +
146 146 # Plot data
147 147 plt.loglog(E, bgrate, color+'-', label=name)
148 148 plt.loglog(E, bgrate, color+'o')
... ... @@ -171,7 +171,7 @@ def show_one_sensitivity(rsp, name, color=&quot;r&quot;, duration=180000.0, alpha=0.2, sig
171 171  
172 172 # Set to figure 3
173 173 plt.figure(3)
174   -
  174 +
175 175 # Generate logE(TeV) vector
176 176 nbins = 20
177 177 logE = [-1.7+0.2*i for i in range(nbins)]
... ... @@ -200,7 +200,7 @@ def show_one_sensitivity(rsp, name, color=&quot;r&quot;, duration=180000.0, alpha=0.2, sig
200 200 plaw = gammalib.GModelSpectralPlaw(1.0e-6, -2.6, epivot)
201 201 conv = TeV2erg*E[i]*E[i]/plaw.flux(emin, emax)
202 202 flux[i] = conv * src_counts / (duration * aeff*0.68)
203   -
  203 +
204 204 # Plot data
205 205 plt.loglog(E, flux, color+'-', label=name)
206 206 plt.loglog(E, flux, color+'o')
... ... @@ -243,11 +243,11 @@ def show_one_response(rspname, dbname, name, rootdir=None, color=&quot;r&quot;):
243 243  
244 244 # Show sensitivity
245 245 show_one_sensitivity(rsp, name, color=color)
246   -
  246 +
247 247 # Return
248 248 return
249   -
250   -
  249 +
  250 +
251 251 # ======================== #
252 252 # Main routine entry point #
253 253 # ======================== #
... ... @@ -291,7 +291,7 @@ if __name__ == &#39;__main__&#39;:
291 291 rootdir = rsp['rootdir']
292 292 else:
293 293 rootdir = None
294   -
  294 +
295 295 # Show response
296 296 show_one_response(rspname, dbname, name, rootdir=rootdir, color=color)
297 297  
... ...
examples/show_spectrum.py
  View file @5c53c14
... ... @@ -26,18 +26,18 @@ try:
26 26 except:
27 27 sys.exit("This script needs matplotlib")
28 28  
29   -
  29 +
30 30 # ============ #
31 31 # Script entry #
32 32 # ============ #
33 33 if __name__ == "__main__":
34 34 """
35 35 """
36   -
  36 +
37 37 # Check for given spectrum file
38 38 if not len(sys.argv) == 2:
39 39 sys.exit("Usage: show_spectrum.py spectrum.fits")
40   -
  40 +
41 41 # Read spectrum file
42 42 filename = sys.argv[1]
43 43 fits = gammalib.GFits(filename)
... ... @@ -49,7 +49,7 @@ if __name__ == &quot;__main__&quot;:
49 49 c_eflux = table["e_Flux"]
50 50 c_ts = table["TS"]
51 51 c_upper = table["UpperLimit"]
52   -
  52 +
53 53 # Initialise arrays to be filled
54 54 energies = []
55 55 flux = []
... ... @@ -60,7 +60,7 @@ if __name__ == &quot;__main__&quot;:
60 60 ul_ed_engs = []
61 61 ul_eu_engs = []
62 62 ul_flux = []
63   -
  63 +
64 64 # Loop over rows of the file
65 65 nrows = table.nrows()
66 66 for row in range(nrows):
... ... @@ -69,7 +69,7 @@ if __name__ == &quot;__main__&quot;:
69 69 ts = c_ts.real(row)
70 70 flx = c_flux.real(row)
71 71 e_flx = c_eflux.real(row)
72   -
  72 +
73 73 # Switch
74 74 if ts > 9.0 and e_flx < flx:
75 75  
... ... @@ -88,7 +88,7 @@ if __name__ == &quot;__main__&quot;:
88 88 ul_flux.append(c_upper.real(row))
89 89 ul_ed_engs.append(c_ed.real(row))
90 90 ul_eu_engs.append(c_eu.real(row))
91   -
  91 +
92 92 # Plot the spectrum
93 93 plt.figure()
94 94 plt.loglog()
... ... @@ -98,4 +98,3 @@ if __name__ == &quot;__main__&quot;:
98 98 plt.xlabel("Energy [TeV]")
99 99 plt.ylabel(r"dN/dE [erg cm$^{-2}$ s$^{-1}$]")
100 100 plt.show()
101   -
102 101 \ No newline at end of file
... ...
examples/show_ts_distribution.py
  View file @5c53c14
... ... @@ -34,249 +34,249 @@ import math
34 34 # Read TS distribution #
35 35 # ==================== #
36 36 def read_ts(filename, tsname="TS"):
37   - """
38   - Read TS distribution values from CSV file.
39   - """
40   - # Initialise list
41   - values = []
42   -
43   - # Open reader
44   - reader = csv.reader(open(filename, 'r'), delimiter=',')
45   -
46   - # Read rows
47   - first = True
48   - index = -1
49   - for row in reader:
50   -
51   - # Get column index if first row
52   - if first:
53   - try:
54   - index = row.index(tsname)
55   - except:
56   - sys.stdout.write('ERROR: Column "TS" not found in file\n')
57   - sys.stdout.write(row+"\n")
58   - raise NameError("TS")
59   -
60   - # Handle data rows
61   - else:
62   - values.append(float(row[index]))
63   -
64   - # Flag that first row has been passed
65   - first = False
66   -
67   - # Create numpy array
68   - a = np.array(values)
69   -
70   - # Return array
71   - return a
  37 + """
  38 + Read TS distribution values from CSV file.
  39 + """
  40 + # Initialise list
  41 + values = []
  42 +
  43 + # Open reader
  44 + reader = csv.reader(open(filename, 'r'), delimiter=',')
  45 +
  46 + # Read rows
  47 + first = True
  48 + index = -1
  49 + for row in reader:
  50 +
  51 + # Get column index if first row
  52 + if first:
  53 + try:
  54 + index = row.index(tsname)
  55 + except:
  56 + sys.stdout.write('ERROR: Column "TS" not found in file\n')
  57 + sys.stdout.write(row+"\n")
  58 + raise NameError("TS")
  59 +
  60 + # Handle data rows
  61 + else:
  62 + values.append(float(row[index]))
  63 +
  64 + # Flag that first row has been passed
  65 + first = False
  66 +
  67 + # Create numpy array
  68 + a = np.array(values)
  69 +
  70 + # Return array
  71 + return a
72 72  
73 73  
74 74 # =========== #
75 75 # Compute erf #
76 76 # =========== #
77 77 def erf(x):
78   - """
79   - Compute error function.
80   - """
81   - # save the sign of x
82   - sign = 1
83   - if x < 0:
84   - sign = -1
85   - x = abs(x)
86   -
87   - # constants
88   - a1 = 0.254829592
89   - a2 = -0.284496736
90   - a3 = 1.421413741
91   - a4 = -1.453152027
92   - a5 = 1.061405429
93   - p = 0.3275911
94   -
95   - # A&S formula 7.1.26
96   - t = 1.0/(1.0 + p*x)
97   - y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)
98   -
99   - # Return result
100   - return sign*y # erf(-x) = -erf(x)
  78 + """
  79 + Compute error function.
  80 + """
  81 + # save the sign of x
  82 + sign = 1
  83 + if x < 0:
  84 + sign = -1
  85 + x = abs(x)
  86 +
  87 + # constants
  88 + a1 = 0.254829592
  89 + a2 = -0.284496736
  90 + a3 = 1.421413741
  91 + a4 = -1.453152027
  92 + a5 = 1.061405429
  93 + p = 0.3275911
  94 +
  95 + # A&S formula 7.1.26
  96 + t = 1.0/(1.0 + p*x)
  97 + y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)
  98 +
  99 + # Return result
  100 + return sign*y # erf(-x) = -erf(x)
101 101  
102 102  
103 103 # =============================== #
104 104 # Display cumulative distribution #
105 105 # =============================== #
106 106 def dist_cdf(values, nbins, title):
107   - """
108   - Displays cumulative TS distribution and overplots expectation.
109   - """
110   - # Set range and adapt number of bins. We make sure that we have a
111   - # histogram bin centred on 0 that will capture the TS=0 part.
112   - max_value = max(values)
113   - binsize = max_value/nbins
114   - min_value = 0.5*binsize
115   - while min_value > min(values):
116   - min_value -= binsize
117   - nbins += 1
118   -
119   - # Create histogram
120   - hist_x = []
121   - hist_y = []
122   - hist_expect = []
123   - norm = 0.5*len(values)
124   - for i in range(nbins):
125   - x_min = min_value+i*binsize
126   - x_max = min_value+(i+1)*binsize
127   - x_val = 0.5*(x_min+x_max)
128   - y_val = 0.0
129   - for ts in values:
130   - if ts > x_min:
131   - y_val += 1.0
132   - hist_x.append(x_min)
133   - hist_y.append(y_val)
134   - hist_x.append(x_max)
135   - hist_y.append(y_val)
136   - if x_min > 0.0:
137   - value = norm*erf(math.sqrt(0.5*x_min)) + norm
138   - hist_expect.append(len(values)-value)
139   - hist_expect.append(len(values)-value)
140   - else:
141   - hist_expect.append(norm)
142   - hist_expect.append(norm)
143   -
144   - # Show expected distribution (half)
145   - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label="simulations")
146   - l = plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label="expectation")
147   - l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2)
148   -
149   - # Set plot
150   - plt.xlabel('TS')
151   - plt.ylabel('Number of trials > TS')
152   - plt.title(title)
153   - plt.grid(True)
154   - plt.legend(loc="upper right")
155   -
156   - # Show histogram
157   - plt.show()
158   -
159   - # Return
160   - return
  107 + """
  108 + Displays cumulative TS distribution and overplots expectation.
  109 + """
  110 + # Set range and adapt number of bins. We make sure that we have a
  111 + # histogram bin centred on 0 that will capture the TS=0 part.
  112 + max_value = max(values)
  113 + binsize = max_value/nbins
  114 + min_value = 0.5*binsize
  115 + while min_value > min(values):
  116 + min_value -= binsize
  117 + nbins += 1
  118 +
  119 + # Create histogram
  120 + hist_x = []
  121 + hist_y = []
  122 + hist_expect = []
  123 + norm = 0.5*len(values)
  124 + for i in range(nbins):
  125 + x_min = min_value+i*binsize
  126 + x_max = min_value+(i+1)*binsize
  127 + x_val = 0.5*(x_min+x_max)
  128 + y_val = 0.0
  129 + for ts in values:
  130 + if ts > x_min:
  131 + y_val += 1.0
  132 + hist_x.append(x_min)
  133 + hist_y.append(y_val)
  134 + hist_x.append(x_max)
  135 + hist_y.append(y_val)
  136 + if x_min > 0.0:
  137 + value = norm*erf(math.sqrt(0.5*x_min)) + norm
  138 + hist_expect.append(len(values)-value)
  139 + hist_expect.append(len(values)-value)
  140 + else:
  141 + hist_expect.append(norm)
  142 + hist_expect.append(norm)
  143 +
  144 + # Show expected distribution (half)
  145 + l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2, label="simulations")
  146 + l = plt.semilogy(hist_x, hist_expect, 'k-', linewidth=2, label="expectation")
  147 + l = plt.semilogy(hist_x, hist_y, 'r-', linewidth=2)
  148 +
  149 + # Set plot
  150 + plt.xlabel('TS')
  151 + plt.ylabel('Number of trials > TS')
  152 + plt.title(title)
  153 + plt.grid(True)
  154 + plt.legend(loc="upper right")
  155 +
  156 + # Show histogram
  157 + plt.show()
  158 +
  159 + # Return
  160 + return
161 161  
162 162  
163 163 # ======================================== #
164 164 # Display probability density distribution #
165 165 # ======================================== #
166 166 def dist_pdf(values, nbins, title):
167   - """
168   - Display probability density TS distribution and overplots
169   - expectation.
170   - """
171   - # Set range and adapt number of bins. We make sure that we have a
172   - # histogram bin centred on 0 that will capture the TS=0 part.
173   - max_value = max(values)
174   - binsize = max_value/nbins
175   - min_value = 0.5*binsize
176   - while min_value > min(values):
177   - min_value -= binsize
178   - nbins += 1
179   -
180   - # Create histogram
181   - n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value],
182   - align='mid', facecolor='green',
183   - log=True, label="simulations")
184   -
185   - # Create expected distribution (onyl for positive TS). We compute the
186   - # full and half of the distribution, as for positively constrained
187   - # quantities half of the values have TS=0
188   - x = []
189   - y_full = []
190   - y_half = []
191   - width = bins[1]-bins[0]
192   - norm = len(values)/(math.sqrt(2.0)*math.sqrt(math.pi))*width
193   - x.append(0.0)
194   - y_full.append(len(values))
195   - y_half.append(0.5*len(values))
196   - for bin in bins:
197   - ts = bin + 0.5*width
198   - if ts > 0:
199   - y = norm*math.pow(ts,-0.5)*math.exp(-0.5*ts)
200   - x.append(ts)
201   - y_full.append(y)
202   - y_half.append(0.5*y)
203   -
204   - # Show expected distribution (half)
205   - l = plt.semilogy(x, y_half, 'ro', linewidth=2, label="expectation")
206   -
207   - # Set plot
208   - plt.xlabel('TS')
209   - plt.ylabel('Number of trials')
210   - plt.title(title)
211   - plt.grid(True)
212   - plt.legend(loc="upper right")
213   -
214   - # Show histogram
215   - plt.show()
216   -
217   -
  167 + """
  168 + Display probability density TS distribution and overplots
  169 + expectation.
  170 + """
  171 + # Set range and adapt number of bins. We make sure that we have a
  172 + # histogram bin centred on 0 that will capture the TS=0 part.
  173 + max_value = max(values)
  174 + binsize = max_value/nbins
  175 + min_value = 0.5*binsize
  176 + while min_value > min(values):
  177 + min_value -= binsize
  178 + nbins += 1
  179 +
  180 + # Create histogram
  181 + n, bins, patches = plt.hist(values, nbins, range=[min_value, max_value],
  182 + align='mid', facecolor='green',
  183 + log=True, label="simulations")
  184 +
  185 + # Create expected distribution (onyl for positive TS). We compute the
  186 + # full and half of the distribution, as for positively constrained
  187 + # quantities half of the values have TS=0
  188 + x = []
  189 + y_full = []
  190 + y_half = []
  191 + width = bins[1]-bins[0]
  192 + norm = len(values)/(math.sqrt(2.0)*math.sqrt(math.pi))*width
  193 + x.append(0.0)
  194 + y_full.append(len(values))
  195 + y_half.append(0.5*len(values))
  196 + for bin in bins:
  197 + ts = bin + 0.5*width
  198 + if ts > 0:
  199 + y = norm*math.pow(ts,-0.5)*math.exp(-0.5*ts)
  200 + x.append(ts)
  201 + y_full.append(y)
  202 + y_half.append(0.5*y)
  203 +
  204 + # Show expected distribution (half)
  205 + l = plt.semilogy(x, y_half, 'ro', linewidth=2, label="expectation")
  206 +
  207 + # Set plot
  208 + plt.xlabel('TS')
  209 + plt.ylabel('Number of trials')
  210 + plt.title(title)
  211 + plt.grid(True)
  212 + plt.legend(loc="upper right")
  213 +
  214 + # Show histogram
  215 + plt.show()
  216 +
  217 +
218 218 # ======================== #
219 219 # Main routine entry point #
220 220 # ======================== #
221 221 if __name__ == '__main__':
222   - """
223   - Display TS distribution generated using cstsdist using matplotlib.
224   - """
225   - # Print usage information
226   - usage = "Usage: show_ts_distribution filename" \
227   - " [-n bins] [-c column] [-t title] [-p plot]"
228   - if len(sys.argv) < 2:
229   - sys.stdout.write(usage+"\n")
230   - sys.exit()
231   -
232   - # Extract parameters
233   - filename = sys.argv[1]
234   -
235   - # Set default parameters
236   - nbins = 30
237   - tsname = "TS"
238   - title = "TS distribution"
239   - plot = "pdf"
240   -
241   - # Parameter dictionnary
242   - pars = [{'option': '-n', 'value': nbins},
243   - {'option': '-c', 'value': tsname},
244   - {'option': '-t', 'value': title},
245   - {'option': '-p', 'value': plot}]
246   -
247   - # Gather parameters from command line
248   - i = 2
249   - while i < len(sys.argv):
250   -
251   - # Search for options
252   - for par in pars:
253   - if sys.argv[i] == par['option']:
254   - if len(sys.argv) > i+1:
255   - i += 1
256   - try:
257   - par['value'] = sys.argv[i]
258   - except:
259   - sys.stdout.write(usage+"\n")
260   - sys.exit()
261   - else:
262   - sys.stdout.write(usage+"\n")
263   - sys.exit()
264   -
265   - # Next item
266   - i += 1
267   -
268   - # Recover parameters
269   - nbins = int(pars[0]['value'])
270   - tsname = pars[1]['value']
271   - title = pars[2]['value']
272   - plot = pars[3]['value']
273   -
274   - # Read values from CSV file
275   - values = read_ts(filename, tsname=tsname)
276   - sys.stdout.write(len(values), "values read.\n")
277   -
278   - # Show histogram
279   - if plot == "pdf":
280   - dist_pdf(values, nbins, title)
281   - else:
282   - dist_cdf(values, nbins, title)
  222 + """
  223 + Display TS distribution generated using cstsdist using matplotlib.
  224 + """
  225 + # Print usage information
  226 + usage = "Usage: show_ts_distribution filename" \
  227 + " [-n bins] [-c column] [-t title] [-p plot]"
  228 + if len(sys.argv) < 2:
  229 + sys.stdout.write(usage+"\n")
  230 + sys.exit()
  231 +
  232 + # Extract parameters
  233 + filename = sys.argv[1]
  234 +
  235 + # Set default parameters
  236 + nbins = 30
  237 + tsname = "TS"
  238 + title = "TS distribution"
  239 + plot = "pdf"
  240 +
  241 + # Parameter dictionnary
  242 + pars = [{'option': '-n', 'value': nbins},
  243 + {'option': '-c', 'value': tsname},
  244 + {'option': '-t', 'value': title},
  245 + {'option': '-p', 'value': plot}]
  246 +
  247 + # Gather parameters from command line
  248 + i = 2
  249 + while i < len(sys.argv):
  250 +
  251 + # Search for options
  252 + for par in pars:
  253 + if sys.argv[i] == par['option']:
  254 + if len(sys.argv) > i+1:
  255 + i += 1
  256 + try:
  257 + par['value'] = sys.argv[i]
  258 + except:
  259 + sys.stdout.write(usage+"\n")
  260 + sys.exit()
  261 + else:
  262 + sys.stdout.write(usage+"\n")
  263 + sys.exit()
  264 +
  265 + # Next item
  266 + i += 1
  267 +
  268 + # Recover parameters
  269 + nbins = int(pars[0]['value'])
  270 + tsname = pars[1]['value']
  271 + title = pars[2]['value']
  272 + plot = pars[3]['value']
  273 +
  274 + # Read values from CSV file
  275 + values = read_ts(filename, tsname=tsname)
  276 + sys.stdout.write(len(values), "values read.\n")
  277 +
  278 + # Show histogram
  279 + if plot == "pdf":
  280 + dist_pdf(values, nbins, title)
  281 + else:
  282 + dist_cdf(values, nbins, title)
... ...
examples/simulate_events.py
  View file @5c53c14
... ... @@ -86,10 +86,10 @@ def createobs(ra=86.171648, dec=-1.4774586, rad=5.0,
86 86 events.gti(gti)
87 87 events.ebounds(ebounds)
88 88 obs.events(events)
89   -
  89 +
90 90 # Set instrument response
91 91 obs.response(irf, db)
92   -
  92 +
93 93 # Set ontime, livetime, and deadtime correction factor
94 94 obs.ontime(duration)
95 95 obs.livetime(duration*deadc)
... ... @@ -155,7 +155,7 @@ def simulate_ctobssim(obs, xmlname, seed=0):
155 155  
156 156 # Run simulator
157 157 sim.run()
158   -
  158 +
159 159 # Retrieve events
160 160 events = sim.obs()[0].events().copy()
161 161  
... ... @@ -169,7 +169,7 @@ def simulate_ctobssim(obs, xmlname, seed=0):
169 169  
170 170 # Return events
171 171 return events
172   -
  172 +
173 173  
174 174 # =========== #
175 175 # Show events #
... ...
test/test_ctbin.py
  View file @5c53c14
... ... @@ -77,7 +77,7 @@ class Test(gammalib.GPythonTestSuite):
77 77 bin["proj"] = "CAR"
78 78 bin["xref"] = 83.63
79 79 bin["yref"] = 22.01
80   -
  80 +
81 81 # Run tool
82 82 self.test_try("Run ctbin")
83 83 try:
... ...
test/test_ctbkgcube.py
  View file @5c53c14
... ... @@ -85,7 +85,7 @@ class Test(gammalib.GPythonTestSuite):
85 85 bkgcube["proj"] = "CAR"
86 86 bkgcube["xref"] = 83.63
87 87 bkgcube["yref"] = 22.01
88   -
  88 +
89 89 # Run tool
90 90 self.test_try("Run ctbkgcube")
91 91 try:
... ...
test/test_ctbutterfly.py
  View file @5c53c14
... ... @@ -75,7 +75,7 @@ class Test(gammalib.GPythonTestSuite):
75 75 butterfly["emin"] = 0.1
76 76 butterfly["emax"] = 100
77 77 butterfly["outfile"] = "butterfly.txt"
78   -
  78 +
79 79 # Run tool
80 80 self.test_try("Run ctbutterfly")
81 81 try:
... ... @@ -91,6 +91,6 @@ class Test(gammalib.GPythonTestSuite):
91 91 self.test_try_success()
92 92 except:
93 93 self.test_try_failure("Exception occured in saving results.")
94   -
  94 +
95 95 # Return
96 96 return
... ...
test/test_ctcubemask.py
  View file @5c53c14
... ... @@ -73,7 +73,7 @@ class Test(gammalib.GPythonTestSuite):
73 73 mask["rad"] = 2.0
74 74 mask["emin"] = 0.1
75 75 mask["emax"] = 100
76   -
  76 +
77 77 # Run tool
78 78 self.test_try("Run ctcubemask")
79 79 try:
... ...
test/test_cterror.py
  View file @5c53c14
... ... @@ -74,7 +74,7 @@ class Test(gammalib.GPythonTestSuite):
74 74 error["srcname"] = "Crab"
75 75 error["caldb"] = self.caldb
76 76 error["irf"] = self.irf
77   -
  77 +
78 78 # Run tool
79 79 self.test_try("Run cterror")
80 80 try:
... ... @@ -90,6 +90,6 @@ class Test(gammalib.GPythonTestSuite):
90 90 self.test_try_success()
91 91 except:
92 92 self.test_try_failure("Exception occured in saving results.")
93   -
  93 +
94 94 # Return
95 95 return
... ...
test/test_ctexpcube.py
  View file @5c53c14
... ... @@ -82,7 +82,7 @@ class Test(gammalib.GPythonTestSuite):
82 82 expcube["proj"] = "CAR"
83 83 expcube["xref"] = 83.63
84 84 expcube["yref"] = 22.01
85   -
  85 +
86 86 # Run tool
87 87 self.test_try("Run ctexpcube")
88 88 try:
... ...
test/test_ctlike.py
  View file @5c53c14
... ... @@ -88,6 +88,6 @@ class Test(gammalib.GPythonTestSuite):
88 88 self.test_try_success()
89 89 except:
90 90 self.test_try_failure("Exception occured in saving results.")
91   -
  91 +
92 92 # Return
93 93 return
... ...
test/test_ctmodel.py
  View file @5c53c14
... ... @@ -90,7 +90,7 @@ class Test(gammalib.GPythonTestSuite):
90 90 model["proj"] = "CAR"
91 91 model["xref"] = 83.63
92 92 model["yref"] = 22.01
93   -
  93 +
94 94 # Run tool
95 95 self.test_try("Run ctmodel")
96 96 try:
... ...
test/test_ctobssim.py
  View file @5c53c14
... ... @@ -78,7 +78,7 @@ class Test(gammalib.GPythonTestSuite):
78 78 sim["tmax"] = 1800.0
79 79 sim["emin"] = 0.1
80 80 sim["emax"] = 100.0
81   -
  81 +
82 82 # Run tool
83 83 self.test_try("Run ctobssim")
84 84 try:
... ... @@ -98,7 +98,7 @@ class Test(gammalib.GPythonTestSuite):
98 98 self.test_value(pnt.dir().ra_deg(), 83.63, 1.0e-6, "ROI Right Ascension is not 83.63 deg")
99 99 self.test_value(pnt.dir().dec_deg(), 22.01, 1.0e-6, "ROI Declination is not 22.01 deg")
100 100 self.test_value(evt.size(), 4134, "Number of events is not 4134")
101   -
  101 +
102 102 # Save events
103 103 self.test_try("Save events")
104 104 try:
... ... @@ -114,7 +114,7 @@ class Test(gammalib.GPythonTestSuite):
114 114 """
115 115 # Set-up observation container
116 116 obs = self.set_obs(4)
117   -
  117 +
118 118 # Set-up ctobssim
119 119 sim = ctools.ctobssim(obs)
120 120 sim["outevents"] = "sim_events.xml"
... ... @@ -185,18 +185,18 @@ class Test(gammalib.GPythonTestSuite):
185 185 pnt = gammalib.GCTAPointing()
186 186 pnt.dir(pntdir)
187 187 obs_cta.pointing(pnt)
188   -
  188 +
189 189 # Set ROI
190 190 roi = gammalib.GCTARoi()
191 191 instdir = gammalib.GCTAInstDir()
192 192 instdir.dir(pntdir)
193 193 roi.centre(instdir)
194 194 roi.radius(rad)
195   -
  195 +
196 196 # Set GTI
197 197 gti = gammalib.GGti()
198 198 gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration))
199   -
  199 +
200 200 # Set energy boundaries
201 201 ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), \
202 202 gammalib.GEnergy(emax, "TeV"))
... ... @@ -207,16 +207,15 @@ class Test(gammalib.GPythonTestSuite):
207 207 events.gti(gti)
208 208 events.ebounds(ebounds)
209 209 obs_cta.events(events)
210   -
  210 +
211 211 # Set instrument response
212 212 obs_cta.response(self.irf, db)
213   -
  213 +
214 214 # Set ontime, livetime, and deadtime correction factor
215 215 obs_cta.ontime(duration)
216 216 obs_cta.livetime(duration*deadc)
217 217 obs_cta.deadc(deadc)
218 218 obs_cta.id(id)
219   -
  219 +
220 220 # Return CTA observation
221 221 return obs_cta
222   -
223 222 \ No newline at end of file
... ...
test/test_ctpsfcube.py
  View file @5c53c14
... ... @@ -84,7 +84,7 @@ class Test(gammalib.GPythonTestSuite):
84 84 psfcube["yref"] = 22.01
85 85 psfcube["amax"] = 0.3
86 86 psfcube["anumbins"] = 10
87   -
  87 +
88 88 # Run tool
89 89 self.test_try("Run ctpsfcube")
90 90 try:
... ...
test/test_ctselect.py
  View file @5c53c14
... ... @@ -89,6 +89,6 @@ class Test(gammalib.GPythonTestSuite):
89 89 self.test_try_success()
90 90 except:
91 91 self.test_try_failure("Exception occured in saving events.")
92   -
  92 +
93 93 # Return
94 94 return
... ...
test/test_ctskymap.py
  View file @5c53c14
... ... @@ -75,7 +75,7 @@ class Test(gammalib.GPythonTestSuite):
75 75 skymap["proj"] = "CAR"
76 76 skymap["xref"] = 83.63
77 77 skymap["yref"] = 22.01
78   -
  78 +
79 79 # Run tool
80 80 self.test_try("Run ctskymap")
81 81 try:
... ...
test/test_cttsmap.py
  View file @5c53c14
... ... @@ -96,6 +96,6 @@ class Test(gammalib.GPythonTestSuite):
96 96 self.test_try_success()
97 97 except:
98 98 self.test_try_failure("Exception occured in saving results.")
99   -
  99 +
100 100 # Return
101 101 return
... ...
test/test_ctulimit.py
  View file @5c53c14
... ... @@ -72,7 +72,7 @@ class Test(gammalib.GPythonTestSuite):
72 72 ulimit["srcname"] = "Crab"
73 73 ulimit["caldb"] = self.caldb
74 74 ulimit["irf"] = self.irf
75   -
  75 +
76 76 # Run tool
77 77 self.test_try("Run ctulimit")
78 78 try:
... ... @@ -88,6 +88,6 @@ class Test(gammalib.GPythonTestSuite):
88 88 self.test_try_success()
89 89 except:
90 90 self.test_try_failure("Exception occured in saving results.")
91   -
  91 +
92 92 # Return
93 93 return
... ...
test/test_pipelines.py
  View file @5c53c14
... ... @@ -75,7 +75,7 @@ class Test(gammalib.GPythonTestSuite):
75 75 emin = 0.1
76 76 emax = 100.0
77 77 rad_select = 3.0
78   -
  78 +
79 79 # Simulate events
80 80 sim = ctools.ctobssim()
81 81 sim["inmodel"] = model_name
... ... @@ -95,7 +95,7 @@ class Test(gammalib.GPythonTestSuite):
95 95 self.test_try_success()
96 96 except:
97 97 self.test_try_failure("Exception occured in ctobssim.")
98   -
  98 +
99 99 # Select events
100 100 select = ctools.ctselect()
101 101 select["inobs"] = events_name
... ... @@ -113,7 +113,7 @@ class Test(gammalib.GPythonTestSuite):
113 113 self.test_try_success()
114 114 except:
115 115 self.test_try_failure("Exception occured in ctselect.")
116   -
  116 +
117 117 # Perform maximum likelihood fitting
118 118 like = ctools.ctlike()
119 119 like["inobs"] = selected_events_name
... ... @@ -146,7 +146,7 @@ class Test(gammalib.GPythonTestSuite):
146 146 emin = 0.1
147 147 emax = 100.0
148 148 rad_select = 3.0
149   -
  149 +
150 150 # Simulate events
151 151 sim = ctools.ctobssim()
152 152 sim["inmodel"] = model_name
... ... @@ -165,7 +165,7 @@ class Test(gammalib.GPythonTestSuite):
165 165 self.test_try_success()
166 166 except:
167 167 self.test_try_failure("Exception occured in ctobssim.")
168   -
  168 +
169 169 # Select events
170 170 select = ctools.ctselect(sim.obs())
171 171 select["ra"] = ra
... ... @@ -181,7 +181,7 @@ class Test(gammalib.GPythonTestSuite):
181 181 self.test_try_success()
182 182 except:
183 183 self.test_try_failure("Exception occured in ctselect.")
184   -
  184 +
185 185 # Perform maximum likelihood fitting
186 186 like = ctools.ctlike(select.obs())
187 187 self.test_try("Run ctlike")
... ...
test/test_python.py
  View file @5c53c14
... ... @@ -110,7 +110,7 @@ if __name__ == &#39;__main__&#39;:
110 110  
111 111 # Run test suite
112 112 success = suites.run()
113   -
  113 +
114 114  
115 115 # Save test results
116 116 suites.save("reports/ctools.xml")
... ...