Maria Haupt / ctools

Commit a448964cd5bb31bd5ef515cadb62ca5f0d9df91d

Authored by Jürgen Knödlseder
1 parent 08e0327b

Add logL and event statistics to pull distribution output file.

Inline Side-by-side
Showing 1 changed file with 311 additions and 305 deletions
scripts/cspull.py
  View file @a448964
@@ -30,332 +30,338 @@ import math @@ -30,332 +30,338 @@ import math
30 # cspull class # 30 # cspull class #
31 # ============ # 31 # ============ #
32 class cspull(GApplication): 32 class cspull(GApplication):
33 - """  
34 - This class implements the pull distribution generation script. It derives  
35 - from the GammaLib::GApplication class which provides support for parameter  
36 - files, command line arguments, and logging. In that way the Python  
37 - script behaves just as a regular ctool.  
38 - """  
39 - def __init__(self, *argv):  
40 - """  
41 - Constructor.  
42 - """  
43 - # Set name  
44 - self.name = "cspull"  
45 - self.version = "0.2.0"  
46 -  
47 - # Initialise some members  
48 - self.obs = None  
49 - self.model = None  
50 - self.m_srcmdl = None  
51 -  
52 - # Make sure that parfile exists  
53 - file = self.parfile() 33 + """
  34 + This class implements the pull distribution generation script. It derives
  35 + from the GammaLib::GApplication class which provides support for parameter
  36 + files, command line arguments, and logging. In that way the Python
  37 + script behaves just as a regular ctool.
  38 + """
  39 + def __init__(self, *argv):
  40 + """
  41 + Constructor.
  42 + """
  43 + # Set name
  44 + self.name = "cspull"
  45 + self.version = "0.2.0"
  46 +
  47 + # Initialise some members
  48 + self.obs = None
  49 + self.model = None
  50 + self.m_srcmdl = None
  51 +
  52 + # Make sure that parfile exists
  53 + file = self.parfile()
54 54
55 - # Initialise application  
56 - if len(argv) == 0:  
57 - GApplication.__init__(self, self.name, self.version)  
58 - elif len(argv) ==1:  
59 - GApplication.__init__(self, self.name, self.version, *argv)  
60 - else:  
61 - raise TypeError("Invalid number of arguments given.") 55 + # Initialise application
  56 + if len(argv) == 0:
  57 + GApplication.__init__(self, self.name, self.version)
  58 + elif len(argv) ==1:
  59 + GApplication.__init__(self, self.name, self.version, *argv)
  60 + else:
  61 + raise TypeError("Invalid number of arguments given.")
62 62
63 - # Set logger properties  
64 - self.log_header()  
65 - self.log.date(True) 63 + # Set logger properties
  64 + self.log_header()
  65 + self.log.date(True)
66 66
67 - # Return  
68 - return  
69 -  
70 - def __del__(self):  
71 - """  
72 - Destructor.  
73 - """  
74 - # Write separator into logger  
75 - if self.logTerse():  
76 - self.log("\n")  
77 -  
78 - # Return  
79 - return 67 + # Return
  68 + return
  69 +
  70 + def __del__(self):
  71 + """
  72 + Destructor.
  73 + """
  74 + # Write separator into logger
  75 + if self.logTerse():
  76 + self.log("\n")
  77 +
  78 + # Return
  79 + return
80 80
81 - def parfile(self):  
82 - """  
83 - Check if parfile exists. If parfile does not exist then create a  
84 - default parfile. This kluge avoids shipping the cscript with a parfile.  
85 - """  
86 - # Set parfile name  
87 - parfile = self.name+".par"  
88 -  
89 - try:  
90 - pars = GApplicationPars(parfile)  
91 - except:  
92 - # Signal if parfile was not found  
93 - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")  
94 -  
95 - # Create default parfile  
96 - pars = GApplicationPars()  
97 - pars.append(GApplicationPar("srcmdl","f","a","$GAMMALIB/share/models/crab.xml","","","Source model"))  
98 - pars.append(GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))  
99 - pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))  
100 - pars.append(GApplicationPar("caldb","s","h","","","","Calibration database"))  
101 - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))  
102 - pars.append(GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)"))  
103 - pars.append(GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))  
104 - pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))  
105 - pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))  
106 - pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))  
107 - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))  
108 - pars.append(GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))  
109 - pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))  
110 - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))  
111 - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))  
112 - pars.append_standard()  
113 - pars.save(parfile)  
114 -  
115 - # Return  
116 - return  
117 -  
118 - def get_parameters(self):  
119 - """  
120 - Get parameters from parfile and setup the observation.  
121 - """  
122 - # Get parameters  
123 - if self.model == None:  
124 - self.m_srcmdl = self["srcmdl"].filename()  
125 - self.m_outfile = self["outfile"].filename()  
126 - self.m_ntrials = self["ntrials"].integer()  
127 - self.m_caldb = self["caldb"].string()  
128 - self.m_irf = self["irf"].string()  
129 - self.m_ra = self["ra"].real()  
130 - self.m_dec = self["dec"].real()  
131 - self.m_emin = self["emin"].real()  
132 - self.m_emax = self["emax"].real()  
133 - self.m_enumbins = self["enumbins"].integer()  
134 - self.m_duration = self["duration"].real()  
135 - self.m_deadc = self["deadc"].real()  
136 - self.m_rad = self["rad"].real()  
137 - self.m_npix = self["npix"].integer()  
138 - self.m_binsz = self["binsz"].real() 81 + def parfile(self):
  82 + """
  83 + Check if parfile exists. If parfile does not exist then create a
  84 + default parfile. This kluge avoids shipping the cscript with a parfile.
  85 + """
  86 + # Set parfile name
  87 + parfile = self.name+".par"
  88 +
  89 + try:
  90 + pars = GApplicationPars(parfile)
  91 + except:
  92 + # Signal if parfile was not found
  93 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  94 +
  95 + # Create default parfile
  96 + pars = GApplicationPars()
  97 + pars.append(GApplicationPar("srcmdl","f","a","$GAMMALIB/share/models/crab.xml","","","Source model"))
  98 + pars.append(GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))
  99 + pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))
  100 + pars.append(GApplicationPar("caldb","s","h","","","","Calibration database"))
  101 + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  102 + pars.append(GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)"))
  103 + pars.append(GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))
  104 + pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
  105 + pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  106 + pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
  107 + pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  108 + pars.append(GApplicationPar("deadc","r","h","0.95","","","Deadtime correction factor"))
  109 + pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  110 + pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned"))
  111 + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  112 + pars.append_standard()
  113 + pars.save(parfile)
  114 +
  115 + # Return
  116 + return
  117 +
  118 + def get_parameters(self):
  119 + """
  120 + Get parameters from parfile and setup the observation.
  121 + """
  122 + # Get parameters
  123 + if self.model == None:
  124 + self.m_srcmdl = self["srcmdl"].filename()
  125 + self.m_outfile = self["outfile"].filename()
  126 + self.m_ntrials = self["ntrials"].integer()
  127 + self.m_caldb = self["caldb"].string()
  128 + self.m_irf = self["irf"].string()
  129 + self.m_ra = self["ra"].real()
  130 + self.m_dec = self["dec"].real()
  131 + self.m_emin = self["emin"].real()
  132 + self.m_emax = self["emax"].real()
  133 + self.m_enumbins = self["enumbins"].integer()
  134 + self.m_duration = self["duration"].real()
  135 + self.m_deadc = self["deadc"].real()
  136 + self.m_rad = self["rad"].real()
  137 + self.m_npix = self["npix"].integer()
  138 + self.m_binsz = self["binsz"].real()
139 139
140 - # Set some fixed parameters  
141 - self.m_log = False # Logging in client tools  
142 - self.m_debug = False # Debugging in client tools  
143 -  
144 - # Setup observation  
145 - self.obs = self.set_obs(ra=self.m_ra, dec=self.m_dec, \  
146 - emin=self.m_emin, emax=self.m_emax)  
147 -  
148 - # Load source model  
149 - if self.m_srcmdl != None:  
150 - self.model = GModels(self.m_srcmdl)  
151 -  
152 - # Append source model to observation  
153 - self.obs.models(self.model) 140 + # Set some fixed parameters
  141 + self.m_log = False # Logging in client tools
  142 + self.m_debug = False # Debugging in client tools
  143 +
  144 + # Setup observation
  145 + self.obs = self.set_obs(ra=self.m_ra, dec=self.m_dec, \
  146 + emin=self.m_emin, emax=self.m_emax)
  147 +
  148 + # Load source model
  149 + if self.m_srcmdl != None:
  150 + self.model = GModels(self.m_srcmdl)
  151 +
  152 + # Append source model to observation
  153 + self.obs.models(self.model)
154 154
155 - # Return  
156 - return  
157 -  
158 - def models(self, models):  
159 - """  
160 - Set model.  
161 - """  
162 - # Copy models  
163 - self.model = models.copy()  
164 -  
165 - # Return  
166 - return  
167 -  
168 - def execute(self):  
169 - """  
170 - Execute the script.  
171 - """  
172 - # Run the script  
173 - self.run()  
174 -  
175 - # Return  
176 - return 155 + # Return
  156 + return
  157 +
  158 + def models(self, models):
  159 + """
  160 + Set model.
  161 + """
  162 + # Copy models
  163 + self.model = models.copy()
  164 +
  165 + # Return
  166 + return
  167 +
  168 + def execute(self):
  169 + """
  170 + Execute the script.
  171 + """
  172 + # Run the script
  173 + self.run()
  174 +
  175 + # Return
  176 + return
177 177
178 - def run(self):  
179 - """  
180 - Run the script.  
181 - """  
182 - # Switch screen logging on in debug mode  
183 - if self.logDebug():  
184 - self.log.cout(True) 178 + def run(self):
  179 + """
  180 + Run the script.
  181 + """
  182 + # Switch screen logging on in debug mode
  183 + if self.logDebug():
  184 + self.log.cout(True)
185 185
186 - # Get parameters  
187 - self.get_parameters()  
188 -  
189 - # Write input parameters into logger  
190 - if self.logTerse():  
191 - self.log_parameters()  
192 - self.log("\n")  
193 -  
194 - # Write observation into logger  
195 - if self.logTerse():  
196 - self.log("\n")  
197 - self.log.header1("Observation")  
198 - self.log(str(self.obs))  
199 - self.log("\n") 186 + # Get parameters
  187 + self.get_parameters()
  188 +
  189 + # Write input parameters into logger
  190 + if self.logTerse():
  191 + self.log_parameters()
  192 + self.log("\n")
  193 +
  194 + # Write observation into logger
  195 + if self.logTerse():
  196 + self.log("\n")
  197 + self.log.header1("Observation")
  198 + self.log(str(self.obs))
  199 + self.log("\n")
200 200
201 - # Write header  
202 - if self.logTerse():  
203 - self.log("\n")  
204 - self.log.header1("Generate pull distribution") 201 + # Write header
  202 + if self.logTerse():
  203 + self.log("\n")
  204 + self.log.header1("Generate pull distribution")
205 205
206 - # Loop over trials  
207 - for seed in range(self.m_ntrials):  
208 -  
209 - # Make a trial  
210 - result = self.trial(seed)  
211 -  
212 - # Write out result immediately  
213 - if seed == 0:  
214 - file = open(self.m_outfile, 'w')  
215 - writer = csv.DictWriter(file, result['colnames'])  
216 - writer.writerow(dict((_,_) for _ in result['colnames']))  
217 - else:  
218 - file = open(self.m_outfile, 'a')  
219 - writer = csv.DictWriter(file, result['colnames'])  
220 - writer.writerow(result['values'])  
221 - file.close()  
222 -  
223 - # Return  
224 - return  
225 -  
226 - def set_obs(self, ra=83.6331, dec=22.0145, emin=0.1, emax=100.0):  
227 - """  
228 - Returns an observation container with a single CTA observation.  
229 -  
230 - Keywords:  
231 - ra - Right Ascension of pointing [deg] (default: 83.6331)  
232 - dec - Declination of pointing [deg] (default: 22.0145)  
233 - emin - Minimum energy [TeV] (default: 0.1)  
234 - emax - Maximum energy [TeV] (default: 100.0)  
235 - """  
236 - # Allocate observation container  
237 - obs = GObservations()  
238 -  
239 - # Set single pointing  
240 - pntdir = GSkyDir()  
241 - pntdir.radec_deg(ra, dec)  
242 -  
243 - # Create CTA observation  
244 - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \  
245 - duration=self.m_duration, deadc=self.m_deadc, \ 206 + # Loop over trials
  207 + for seed in range(self.m_ntrials):
  208 +
  209 + # Make a trial
  210 + result = self.trial(seed)
  211 +
  212 + # Write out result immediately
  213 + if seed == 0:
  214 + file = open(self.m_outfile, 'w')
  215 + writer = csv.DictWriter(file, result['colnames'])
  216 + writer.writerow(dict((_,_) for _ in result['colnames']))
  217 + else:
  218 + file = open(self.m_outfile, 'a')
  219 + writer = csv.DictWriter(file, result['colnames'])
  220 + writer.writerow(result['values'])
  221 + file.close()
  222 +
  223 + # Return
  224 + return
  225 +
  226 + def set_obs(self, ra=83.6331, dec=22.0145, emin=0.1, emax=100.0):
  227 + """
  228 + Returns an observation container with a single CTA observation.
  229 +
  230 + Keywords:
  231 + ra - Right Ascension of pointing [deg] (default: 83.6331)
  232 + dec - Declination of pointing [deg] (default: 22.0145)
  233 + emin - Minimum energy [TeV] (default: 0.1)
  234 + emax - Maximum energy [TeV] (default: 100.0)
  235 + """
  236 + # Allocate observation container
  237 + obs = GObservations()
  238 +
  239 + # Set single pointing
  240 + pntdir = GSkyDir()
  241 + pntdir.radec_deg(ra, dec)
  242 +
  243 + # Create CTA observation
  244 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  245 + duration=self.m_duration, deadc=self.m_deadc, \
246 emin=emin, emax=emax, rad=self.m_rad) 246 emin=emin, emax=emax, rad=self.m_rad)
247 -  
248 - # Append observation to container  
249 - obs.append(run)  
250 -  
251 - # Return observation container  
252 - return obs  
253 -  
254 - def trial(self, seed):  
255 - """  
256 - Create the pull for a single trial.  
257 -  
258 - Parameters:  
259 - seed - Random number generator seed  
260 - """  
261 - # Write header  
262 - if self.logExplicit():  
263 - self.log.header2("Trial "+str(seed+1)) 247 +
  248 + # Append observation to container
  249 + obs.append(run)
  250 +
  251 + # Return observation container
  252 + return obs
  253 +
  254 + def trial(self, seed):
  255 + """
  256 + Create the pull for a single trial.
  257 +
  258 + Parameters:
  259 + seed - Random number generator seed
  260 + """
  261 + # Write header
  262 + if self.logExplicit():
  263 + self.log.header2("Trial "+str(seed+1))
264 264
265 - # Simulate events  
266 - obs = obsutils.sim(self.obs, \  
267 - nbins=self.m_enumbins, \  
268 - seed=seed, \  
269 - binsz=self.m_binsz, \  
270 - npix=self.m_npix, \  
271 - log=self.m_log, debug=self.m_debug) 265 + # Simulate events
  266 + obs = obsutils.sim(self.obs, \
  267 + nbins=self.m_enumbins, \
  268 + seed=seed, \
  269 + binsz=self.m_binsz, \
  270 + npix=self.m_npix, \
  271 + log=self.m_log, debug=self.m_debug)
272 272
273 - # Determine number of events in simulation  
274 - nevents = 0.0  
275 - for run in obs:  
276 - nevents += run.events().number() 273 + # Determine number of events in simulation
  274 + nevents = 0.0
  275 + for run in obs:
  276 + nevents += run.events().number()
277 277
278 - # Write simulation results  
279 - if self.logExplicit():  
280 - self.log.header3("Simulation")  
281 - self.log.parformat("Number of simulated events")  
282 - self.log(nevents)  
283 - self.log("\n") 278 + # Write simulation results
  279 + if self.logExplicit():
  280 + self.log.header3("Simulation")
  281 + self.log.parformat("Number of simulated events")
  282 + self.log(nevents)
  283 + self.log("\n")
284 284
285 - # Fit model  
286 - like = obsutils.fit(obs, log=self.m_log, debug=self.m_debug) 285 + # Fit model
  286 + like = obsutils.fit(obs, log=self.m_log, debug=self.m_debug)
287 287
288 - # Store results  
289 - logL = like.opt().value()  
290 - npred = like.obs().npred()  
291 - models = like.obs().models() 288 + # Store results
  289 + logL = like.opt().value()
  290 + npred = like.obs().npred()
  291 + models = like.obs().models()
292 292
293 - # Write result header  
294 - if self.logExplicit():  
295 - self.log.header3("Pulls")  
296 -  
297 - # Gather results  
298 - colnames = []  
299 - values = {}  
300 - for i in range(models.size()):  
301 - model = models[i]  
302 - model_name = model.name()  
303 - for k in range(model.size()):  
304 - par = model[k]  
305 - if par.is_free():  
306 -  
307 - # Set parameter name  
308 - name = model_name+"_"+par.name()  
309 -  
310 - # Append parameter, Pull_parameter and Unc_parameter  
311 - colnames.append(name)  
312 - colnames.append("Pull_"+name)  
313 - colnames.append("Unc_"+name)  
314 -  
315 - # Compute pull  
316 - fitted_value = par.value()  
317 - real_value = self.model[i][k].value()  
318 - error = par.error()  
319 - if error != 0.0:  
320 - pull = (fitted_value - real_value) / error  
321 - else:  
322 - pull = 99.0  
323 -  
324 - # Store results  
325 - values[name] = fitted_value  
326 - values["Pull_"+name] = pull  
327 - values["Unc_"+name] = error 293 + # Write result header
  294 + if self.logExplicit():
  295 + self.log.header3("Pulls")
  296 +
  297 + # Gather results
  298 + colnames = []
  299 + values = {}
  300 + colnames.append("LogL")
  301 + colnames.append("Sim_Events")
  302 + colnames.append("Npred_Events")
  303 + values["LogL"] = logL
  304 + values["Sim_Events"] = nevents
  305 + values["Npred_Events"] = npred
  306 + for i in range(models.size()):
  307 + model = models[i]
  308 + model_name = model.name()
  309 + for k in range(model.size()):
  310 + par = model[k]
  311 + if par.is_free():
  312 +
  313 + # Set parameter name
  314 + name = model_name+"_"+par.name()
  315 +
  316 + # Append parameter, Pull_parameter and Unc_parameter
  317 + colnames.append(name)
  318 + colnames.append("Pull_"+name)
  319 + colnames.append("Unc_"+name)
  320 +
  321 + # Compute pull
  322 + fitted_value = par.value()
  323 + real_value = self.model[i][k].value()
  324 + error = par.error()
  325 + if error != 0.0:
  326 + pull = (fitted_value - real_value) / error
  327 + else:
  328 + pull = 99.0
  329 +
  330 + # Store results
  331 + values[name] = fitted_value
  332 + values["Pull_"+name] = pull
  333 + values["Unc_"+name] = error
328 334
329 - # Write result  
330 - if self.logExplicit():  
331 - self.log.parformat(name)  
332 - self.log(pull)  
333 - self.log(" (")  
334 - self.log(fitted_value)  
335 - self.log(" +/- ")  
336 - self.log(error)  
337 - self.log(")\n")  
338 -  
339 - # Bundle together results  
340 - result = {'colnames': colnames, 'values': values}  
341 -  
342 - # Return  
343 - return result 335 + # Write result
  336 + if self.logExplicit():
  337 + self.log.parformat(name)
  338 + self.log(pull)
  339 + self.log(" (")
  340 + self.log(fitted_value)
  341 + self.log(" +/- ")
  342 + self.log(error)
  343 + self.log(")\n")
  344 +
  345 + # Bundle together results
  346 + result = {'colnames': colnames, 'values': values}
  347 +
  348 + # Return
  349 + return result
344 350
345 351
346 # ======================== # 352 # ======================== #
347 # Main routine entry point # 353 # Main routine entry point #
348 # ======================== # 354 # ======================== #
349 if __name__ == '__main__': 355 if __name__ == '__main__':
350 - """  
351 - Generates pull distribution for a source model.  
352 - """  
353 - # Create instance of application  
354 - app = cspull(sys.argv)  
355 -  
356 - # Open logfile  
357 - app.logFileOpen()  
358 -  
359 - # Execute application  
360 - app.execute()  
361 -  
362 \ No newline at end of file 356 \ No newline at end of file
  357 + """
  358 + Generates pull distribution for a source model.
  359 + """
  360 + # Create instance of application
  361 + app = cspull(sys.argv)
  362 +
  363 + # Open logfile
  364 + app.logFileOpen()
  365 +
  366 + # Execute application
  367 + app.execute()
  368 +
363 \ No newline at end of file 369 \ No newline at end of file