Commit 46feb7955ce9a4d68e3abdf0ca2d663d4782a9af

Authored by knodlseder
1 parent 2cc87dd0

Supports now also reading of observations from an observation definition XML file.

Binned and unbinned analysis mode is now automatically detected from the input file type. The method parameter has been removed, and a infile parameter has been introduced instead of the old evfile and cntfile parameters.
Showing 1 changed file with 68 additions and 104 deletions
src/ctlike/ctlike.cpp
1 /*************************************************************************** 1 /***************************************************************************
2 * ctlike - CTA maximum likelihood tool * 2 * ctlike - CTA maximum likelihood tool *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 - * copyright (C) 2010-2011 by Jurgen Knodlseder * 4 + * copyright (C) 2010-2012 by Jurgen Knodlseder *
5 * ----------------------------------------------------------------------- * 5 * ----------------------------------------------------------------------- *
6 * * 6 * *
7 * This program is free software: you can redistribute it and/or modify * 7 * This program is free software: you can redistribute it and/or modify *
@@ -229,8 +229,9 @@ void ctlike::execute(void) @@ -229,8 +229,9 @@ void ctlike::execute(void)
229 void ctlike::run(void) 229 void ctlike::run(void)
230 { 230 {
231 // Switch screen logging on in debug mode 231 // Switch screen logging on in debug mode
232 - if (logDebug()) 232 + if (logDebug()) {
233 log.cout(true); 233 log.cout(true);
  234 + }
234 235
235 // Get parameters 236 // Get parameters
236 get_parameters(); 237 get_parameters();
@@ -244,10 +245,12 @@ void ctlike::run(void) @@ -244,10 +245,12 @@ void ctlike::run(void)
244 // Write observation(s) into logger 245 // Write observation(s) into logger
245 if (logTerse()) { 246 if (logTerse()) {
246 log << std::endl; 247 log << std::endl;
247 - if (m_obs.size() > 1) 248 + if (m_obs.size() > 1) {
248 log.header1("Observations"); 249 log.header1("Observations");
249 - else 250 + }
  251 + else {
250 log.header1("Observation"); 252 log.header1("Observation");
  253 + }
251 log << m_obs << std::endl; 254 log << m_obs << std::endl;
252 } 255 }
253 256
@@ -283,8 +286,9 @@ void ctlike::save(void) @@ -283,8 +286,9 @@ void ctlike::save(void)
283 m_outmdl = (*this)["outmdl"].value(); 286 m_outmdl = (*this)["outmdl"].value();
284 287
285 // Write results out as XML model 288 // Write results out as XML model
286 - if (toupper(m_outmdl) != "NONE") 289 + if (toupper(m_outmdl) != "NONE") {
287 m_obs.models().save(m_outmdl); 290 m_obs.models().save(m_outmdl);
  291 + }
288 292
289 // Return 293 // Return
290 return; 294 return;
@@ -294,76 +298,88 @@ void ctlike::save(void) @@ -294,76 +298,88 @@ void ctlike::save(void)
294 /***********************************************************************//** 298 /***********************************************************************//**
295 * @brief Get application parameters 299 * @brief Get application parameters
296 * 300 *
297 - * GException::par_error  
298 - * Invalid analysis method specified (neither BINNED nor UNBINNED).  
299 - *  
300 * Get all required task parameters from the parameter file or (if specified) 301 * Get all required task parameters from the parameter file or (if specified)
301 * by querying the user. Observation dependent parameters will only be read 302 * by querying the user. Observation dependent parameters will only be read
302 * if the observation container is actually empty. Observation dependent 303 * if the observation container is actually empty. Observation dependent
303 * parameters are: 304 * parameters are:
304 - * "method" (analysis method),  
305 * "stat" (statistics to be used for observation), 305 * "stat" (statistics to be used for observation),
306 * "caldb" (calibration database), 306 * "caldb" (calibration database),
307 - * "irf" (instrument response function),  
308 - * "evfile" (event file name), and  
309 - * "cntmap" (counts map file name). 307 + * "irf" (instrument response function), and
  308 + * "infile" (input file name).
310 * The model will only be loaded if no model components exist in the 309 * The model will only be loaded if no model components exist in the
311 * observation container. 310 * observation container.
  311 + *
  312 + * This method handles both loading of FITS files and of handling XML
  313 + * observation definition files.
312 ***************************************************************************/ 314 ***************************************************************************/
313 void ctlike::get_parameters(void) 315 void ctlike::get_parameters(void)
314 { 316 {
315 - // If there is are no models associated with the observations then  
316 - // load now the model definition  
317 - if (m_obs.models().size() == 0) { 317 + // If there are no observations in container then add a single CTA
  318 + // observation using the parameters from the parameter file
  319 + if (m_obs.size() == 0) {
318 320
319 - // Get models XML filename  
320 - std::string filename = (*this)["srcmdl"].filename(); 321 + // Allocate CTA observation
  322 + GCTAObservation obs;
321 323
322 - // Setup models for optimizing.  
323 - m_obs.models(GModels(filename)); 324 + // Get event file name
  325 + std::string filename = (*this)["infile"].filename();
324 326
325 - } // endif: no models were associated with observations 327 + // Try first to open as unbinned FITS file
  328 + try {
326 329
327 - // If there are no observations in container then add a single CTA  
328 - // observation using the parameters from the parameter file  
329 - if (m_obs.size() == 0) { 330 + // Determine whether FITS file has EVENTS extension
  331 + GFits file(filename);
  332 + bool is_unbinned = file.hashdu("EVENTS");
  333 + file.close();
330 334
331 - // Get observation parameters  
332 - m_method = toupper((*this)["method"].string());  
333 - m_stat = toupper((*this)["stat"].string());  
334 - m_caldb = (*this)["caldb"].string();  
335 - m_irf = (*this)["irf"].string(); 335 + // If FITS file has an events header then load as unbinned
  336 + if (is_unbinned) {
  337 + obs.load_unbinned(filename);
  338 + }
336 339
337 - // Case A: set-up unbinned CTA observation  
338 - if (m_method == "UNBINNED") { 340 + // ... otherwise load as binned
  341 + else {
  342 + obs.load_binned(filename);
  343 + }
339 344
340 - // Get event file name  
341 - std::string filename = (*this)["evfile"].filename(); 345 + // Get task parameters for event file loading
  346 + m_stat = toupper((*this)["stat"].string());
  347 + m_caldb = (*this)["caldb"].string();
  348 + m_irf = (*this)["irf"].string();
342 349
343 - // Load and append unbinned CTA observation  
344 - load_unbinned(filename); 350 + // Set statistics
  351 + obs.statistics(m_stat);
345 352
346 - } // endif: unbinned analysis mode 353 + // Set reponse
  354 + obs.response(m_irf, m_caldb);
347 355
348 - // Case B: set-up binned CTA observation  
349 - else if (m_method == "BINNED") { 356 + // Append observation to container
  357 + m_obs.append(obs);
350 358
351 - // Get counts map file name  
352 - std::string filename = (*this)["cntmap"].filename(); 359 + }
353 360
354 - // Load and append binned CTA observation  
355 - load_binned(filename); 361 + // ... otherwise try to open as XML file
  362 + catch (GException::fits_open_error) {
356 363
357 - } // endif: binned analysis mode 364 + // Load observations from XML file
  365 + m_obs.load(filename);
358 366
359 - // ... otherwise signal an invalid analysis method  
360 - else {  
361 - throw GException::par_error(G_GET_PARAMETERS, "method",  
362 - "only \"BINNED\" or \"UNBINNED\" supported.");  
363 } 367 }
364 368
365 } // endif: there was no observation in the container 369 } // endif: there was no observation in the container
366 370
  371 + // If there is are no models associated with the observations then
  372 + // load now the model definition
  373 + if (m_obs.models().size() == 0) {
  374 +
  375 + // Get models XML filename
  376 + std::string filename = (*this)["srcmdl"].filename();
  377 +
  378 + // Setup models for optimizing.
  379 + m_obs.models(GModels(filename));
  380 +
  381 + } // endif: no models were associated with observations
  382 +
367 // Get standard parameters 383 // Get standard parameters
368 m_refit = (*this)["refit"].boolean(); 384 m_refit = (*this)["refit"].boolean();
369 385
@@ -404,8 +420,9 @@ void ctlike::optimize_lm(void) @@ -404,8 +420,9 @@ void ctlike::optimize_lm(void)
404 m_obs.optimize(*opt); 420 m_obs.optimize(*opt);
405 421
406 // Optionally refit 422 // Optionally refit
407 - if (m_refit) 423 + if (m_refit) {
408 m_obs.optimize(*opt); 424 m_obs.optimize(*opt);
  425 + }
409 426
410 // Store maximum log likelihood value 427 // Store maximum log likelihood value
411 m_logL = -(opt->value()); 428 m_logL = -(opt->value());
@@ -423,60 +440,6 @@ void ctlike::optimize_lm(void) @@ -423,60 +440,6 @@ void ctlike::optimize_lm(void)
423 } 440 }
424 441
425 442
426 -/***********************************************************************//**  
427 - * @brief Load unbinned observation and append it to observations container  
428 - *  
429 - * @param[in] evfile Events file name.  
430 - ***************************************************************************/  
431 -void ctlike::load_unbinned(const std::string& evfile)  
432 -{  
433 - // Declare CTA observation  
434 - GCTAObservation obs;  
435 -  
436 - // Set statistics  
437 - obs.statistics(m_stat);  
438 -  
439 - // Load data  
440 - obs.load_unbinned(evfile);  
441 -  
442 - // Set reponse  
443 - obs.response(m_irf, m_caldb);  
444 -  
445 - // Append observation to contained  
446 - m_obs.append(obs);  
447 -  
448 - // Return  
449 - return;  
450 -}  
451 -  
452 -  
453 -/***********************************************************************//**  
454 - * @brief Load binned observation and append it to observations container  
455 - *  
456 - * @param[in] cntmap Counts map file name.  
457 - ***************************************************************************/  
458 -void ctlike::load_binned(const std::string& cntmap)  
459 -{  
460 - // Declare CTA observation  
461 - GCTAObservation obs;  
462 -  
463 - // Set statistics  
464 - obs.statistics(m_stat);  
465 -  
466 - // Load data  
467 - obs.load_binned(cntmap);  
468 -  
469 - // Set reponse  
470 - obs.response(m_irf, m_caldb);  
471 -  
472 - // Append observation to contained  
473 - m_obs.append(obs);  
474 -  
475 - // Return  
476 - return;  
477 -}  
478 -  
479 -  
480 /*========================================================================== 443 /*==========================================================================
481 = = 444 = =
482 = Private methods = 445 = Private methods =
@@ -489,7 +452,7 @@ void ctlike::load_binned(const std::string&amp; cntmap) @@ -489,7 +452,7 @@ void ctlike::load_binned(const std::string&amp; cntmap)
489 void ctlike::init_members(void) 452 void ctlike::init_members(void)
490 { 453 {
491 // Initialise members 454 // Initialise members
492 - m_method.clear(); 455 + //m_method.clear();
493 m_stat.clear(); 456 m_stat.clear();
494 m_caldb.clear(); 457 m_caldb.clear();
495 m_irf.clear(); 458 m_irf.clear();
@@ -517,7 +480,7 @@ void ctlike::init_members(void) @@ -517,7 +480,7 @@ void ctlike::init_members(void)
517 void ctlike::copy_members(const ctlike& app) 480 void ctlike::copy_members(const ctlike& app)
518 { 481 {
519 // Copy attributes 482 // Copy attributes
520 - m_method = app.m_method; 483 + //m_method = app.m_method;
521 m_stat = app.m_stat; 484 m_stat = app.m_stat;
522 m_refit = app.m_refit; 485 m_refit = app.m_refit;
523 m_caldb = app.m_caldb; 486 m_caldb = app.m_caldb;
@@ -546,8 +509,9 @@ void ctlike::free_members(void) @@ -546,8 +509,9 @@ void ctlike::free_members(void)
546 m_opt = NULL; 509 m_opt = NULL;
547 510
548 // Write separator into logger 511 // Write separator into logger
549 - if (logTerse()) 512 + if (logTerse()) {
550 log << std::endl; 513 log << std::endl;
  514 + }
551 515
552 // Return 516 // Return
553 return; 517 return;