Commit 07a98b27f51824dc08acd63cc6e5dff62f01fc7c

Authored by Jürgen Knödlseder
1 parent 230a31bb

Integrate code:

- change evfile parameter to infile
- added irf and caldb parameters
- use pointer to GModel for test source
- clean exception handling
- minor code adjustments
pyext/cttsmap.i
1 /*************************************************************************** 1 /***************************************************************************
2 - * cttsmap - CTA data binning tool * 2 + * cttsmap - TS map calculation tool *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 - * copyright (C) 2010-2014 by Juergen Knoedlseder * 4 + * copyright (C) 2014 by Michael Mayer *
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 *
@@ -20,8 +20,8 @@ @@ -20,8 +20,8 @@
20 ***************************************************************************/ 20 ***************************************************************************/
21 /** 21 /**
22 * @file cttsmap.i 22 * @file cttsmap.i
23 - * @brief CTA data binning tool Python interface definition  
24 - * @author Juergen Knoedlseder 23 + * @brief TS map calculation tool interface definition
  24 + * @author Michael Mayer
25 */ 25 */
26 %{ 26 %{
27 /* Put headers and other declarations here that are needed for compilation */ 27 /* Put headers and other declarations here that are needed for compilation */
@@ -30,12 +30,11 @@ @@ -30,12 +30,11 @@
30 %} 30 %}
31 31
32 32
33 -/***********************************************************************  
34 -//**  
35 -* @class cttsmap  
36 -*  
37 -* @brief CTA data binning tool Python interface  
38 -***************************************************************************/ 33 +/***********************************************************************//**
  34 + * @class cttsmap
  35 + *
  36 + * @brief TS map calculation
  37 + ***************************************************************************/
39 class cttsmap : public GApplication { 38 class cttsmap : public GApplication {
40 public: 39 public:
41 // Constructors and destructors 40 // Constructors and destructors
@@ -51,19 +50,13 @@ public: @@ -51,19 +50,13 @@ public:
51 void run(void); 50 void run(void);
52 void save(void); 51 void save(void);
53 const GObservations& obs(void) const; 52 const GObservations& obs(void) const;
54 - const GSkymap& tsmap(void) const;  
55 - const GSkymap& fluxmap(void) const;  
56 - const GSkymap& indexmap(void) const;  
57 - const GSkymap& statusmap(void) const;  
58 -  
59 - 53 + const GSkymap& tsmap(void) const;
60 }; 54 };
61 55
62 56
63 -/***********************************************************************  
64 - //**  
65 -* @brief CTA data binning tool Python extension  
66 -***************************************************************************/ 57 +/***********************************************************************//**
  58 + * @brief TS map calculation tool Python extensions
  59 + ***************************************************************************/
67 %extend cttsmap { 60 %extend cttsmap {
68 cttsmap copy() { 61 cttsmap copy() {
69 return (*self); 62 return (*self);
src/cttsmap/cttsmap.cpp
1 /*************************************************************************** 1 /***************************************************************************
2 - * cttsmap - TS map computation tool * 2 + * cttsmap - TS map calculation tool *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 - * copyright (C) 2010-2014 by Juergen Knoedlseder * 4 + * copyright (C) 2014 by Michael Mayer *
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 *
@@ -19,8 +19,8 @@ @@ -19,8 +19,8 @@
19 * * 19 * *
20 ***************************************************************************/ 20 ***************************************************************************/
21 /** 21 /**
22 - * @file cttsmap.cpp  
23 - * @brief TS map computation tool 22 + * @file cttsmap.hpp
  23 + * @brief TS map calculation tool interface implementation
24 * @author Michael Mayer 24 * @author Michael Mayer
25 */ 25 */
26 26
@@ -33,7 +33,7 @@ @@ -33,7 +33,7 @@
33 #include "GTools.hpp" 33 #include "GTools.hpp"
34 34
35 /* __ Method name definitions ____________________________________________ */ 35 /* __ Method name definitions ____________________________________________ */
36 -#define G_GET_PARAMETERS "cttsmap::get_parameters()" 36 +#define G_GET_PARAMETERS "cttsmap::get_parameters()"
37 37
38 /* __ Debug definitions __________________________________________________ */ 38 /* __ Debug definitions __________________________________________________ */
39 39
@@ -71,7 +71,7 @@ cttsmap::cttsmap(void) : GApplication(CTTSMAP_NAME, CTTSMAP_VERSION) @@ -71,7 +71,7 @@ cttsmap::cttsmap(void) : GApplication(CTTSMAP_NAME, CTTSMAP_VERSION)
71 * observations container. 71 * observations container.
72 ***************************************************************************/ 72 ***************************************************************************/
73 cttsmap::cttsmap(const GObservations& obs) : 73 cttsmap::cttsmap(const GObservations& obs) :
74 - GApplication(CTTSMAP_NAME, CTTSMAP_VERSION) 74 + GApplication(CTTSMAP_NAME, CTTSMAP_VERSION)
75 { 75 {
76 // Initialise members 76 // Initialise members
77 init_members(); 77 init_members();
@@ -95,7 +95,7 @@ cttsmap::cttsmap(const GObservations& obs) : @@ -95,7 +95,7 @@ cttsmap::cttsmap(const GObservations& obs) :
95 * @param[in] argv Array of command line arguments. 95 * @param[in] argv Array of command line arguments.
96 ***************************************************************************/ 96 ***************************************************************************/
97 cttsmap::cttsmap(int argc, char *argv[]) : 97 cttsmap::cttsmap(int argc, char *argv[]) :
98 - GApplication(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv) 98 + GApplication(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
99 { 99 {
100 // Initialise members 100 // Initialise members
101 init_members(); 101 init_members();
@@ -204,8 +204,8 @@ void cttsmap::clear(void) @@ -204,8 +204,8 @@ void cttsmap::clear(void)
204 * This is the main execution method of the cttsmap class. It is invoked when 204 * This is the main execution method of the cttsmap class. It is invoked when
205 * the executable is called from command line. 205 * the executable is called from command line.
206 * 206 *
207 - * The method reads the task parameters, computed the TS values  
208 - * and stores it into map(s), and writes the results into FITS files on disk. 207 + * The method reads the task parameters, computes the TS values, stores the
  208 + * results into map(s), and writes the maps into a FITS file on disk.
209 ***************************************************************************/ 209 ***************************************************************************/
210 void cttsmap::execute(void) 210 void cttsmap::execute(void)
211 { 211 {
@@ -226,10 +226,9 @@ void cttsmap::execute(void) @@ -226,10 +226,9 @@ void cttsmap::execute(void)
226 /***********************************************************************//** 226 /***********************************************************************//**
227 * @brief Computes the TS maps 227 * @brief Computes the TS maps
228 * 228 *
229 - * This method moves a point-like source along a grid of coordinates  
230 - * and refits the given models including this additional source.  
231 - * The fit parameters of the source and its TS values are filled into  
232 - * respective position in the map. 229 + * This method moves a source along a grid of coordinates and refits the
  230 + * given models including this additional source. The fit parameters of the
  231 + * source and its TS values are filled into respective position in the map.
233 ***************************************************************************/ 232 ***************************************************************************/
234 void cttsmap::run(void) 233 void cttsmap::run(void)
235 { 234 {
@@ -262,12 +261,7 @@ void cttsmap::run(void) @@ -262,12 +261,7 @@ void cttsmap::run(void)
262 // Write header 261 // Write header
263 if (logTerse()) { 262 if (logTerse()) {
264 log << std::endl; 263 log << std::endl;
265 - if (m_obs.size() > 1) {  
266 - log.header1("Bin observations");  
267 - }  
268 - else {  
269 - log.header1("Bin observation");  
270 - } 264 + log.header1("Initialise TS map");
271 } 265 }
272 266
273 // Initialise maps 267 // Initialise maps
@@ -281,8 +275,8 @@ void cttsmap::run(void) @@ -281,8 +275,8 @@ void cttsmap::run(void)
281 int binmax = (m_binmax == -1) ? m_tsmap.npix() : m_binmax; 275 int binmax = (m_binmax == -1) ? m_tsmap.npix() : m_binmax;
282 276
283 // Initialise optimizer 277 // Initialise optimizer
284 - GOptimizerLM* opt = (logTerse()) ? new GOptimizerLM(log)  
285 - : new GOptimizerLM(); 278 + GOptimizerLM* opt = (logExplicit()) ? new GOptimizerLM(log)
  279 + : new GOptimizerLM();
286 280
287 // Store initial models 281 // Store initial models
288 GModels models_orig = m_obs.models(); 282 GModels models_orig = m_obs.models();
@@ -296,29 +290,45 @@ void cttsmap::run(void) @@ -296,29 +290,45 @@ void cttsmap::run(void)
296 // Get likelihood of null hypothesis if not given before 290 // Get likelihood of null hypothesis if not given before
297 if (m_logL0 == 0.0) { 291 if (m_logL0 == 0.0) {
298 292
  293 + // Write header
  294 + if (logTerse()) {
  295 + log << std::endl;
  296 + log.header1("Compute NULL Hypothesis for TS computation");
  297 + }
  298 +
299 // Compute likelihood without the test source 299 // Compute likelihood without the test source
300 m_obs.models(models); 300 m_obs.models(models);
301 m_obs.optimize(*opt); 301 m_obs.optimize(*opt);
302 m_logL0 = -(opt->value()); 302 m_logL0 = -(opt->value());
  303 +
303 } 304 }
304 305
305 - // Loop over all observations in the container 306 + // Write header
  307 + if (logTerse()) {
  308 + log << std::endl;
  309 + log.header1("Generate TS map");
  310 + }
  311 +
  312 + // Loop over grid positions
306 for (int i = binmin; i < binmax; ++i) { 313 for (int i = binmin; i < binmax; ++i) {
307 314
308 // Get the coordinate of current bin 315 // Get the coordinate of current bin
309 GSkyDir bincentre = m_tsmap.inx2dir(i); 316 GSkyDir bincentre = m_tsmap.inx2dir(i);
310 317
311 - // Verbosity option  
312 - if ((*this)["chatter"].integer()>2) {  
313 - std::cout<<"Computing TS for bin number "<<i<<" at "<<bincentre.print()<<std::endl; 318 + // Header for verbose logging
  319 + if (logExplicit()) {
  320 + std::string msg = "Computing TS for bin number "+gammalib::str(i)+
  321 + " at "+bincentre.print();
  322 + log << std::endl;
  323 + log.header2(msg);
314 } 324 }
315 325
316 - // Set the point source at current bin position  
317 - m_testsource["RA"].value(bincentre.ra_deg());  
318 - m_testsource["DEC"].value(bincentre.dec_deg());  
319 -  
320 - // Add test source to model container  
321 - models.append(m_testsource); 326 + // Add test source at current bin position
  327 + if (m_testsource != NULL) {
  328 + (*m_testsource)["RA"].value(bincentre.ra_deg());
  329 + (*m_testsource)["DEC"].value(bincentre.dec_deg());
  330 + models.append(*m_testsource);
  331 + }
322 332
323 // Assign models to observations 333 // Assign models to observations
324 m_obs.models(models); 334 m_obs.models(models);
@@ -330,20 +340,34 @@ void cttsmap::run(void) @@ -330,20 +340,34 @@ void cttsmap::run(void)
330 double logL1 = -(opt->value()); 340 double logL1 = -(opt->value());
331 341
332 // Compute TS value 342 // Compute TS value
333 - double ts = 2.0 * (logL1-m_logL0); 343 + double ts = 2.0 * (logL1 - m_logL0);
  344 +
  345 + // Log information
  346 + if (logExplicit()) {
  347 + log << " TS value ....: ";
  348 + log << ts << std::endl;
  349 + }
  350 + else if (logTerse()) {
  351 + log << "TS for bin number ";
  352 + log << i;
  353 + log << " at ";
  354 + log << bincentre.print();
  355 + log << ": ";
  356 + log << ts << std::endl;
  357 + }
334 358
335 // Get test source model instance 359 // Get test source model instance
336 GModels best_fit_model = m_obs.models(); 360 GModels best_fit_model = m_obs.models();
337 - GModel* testsource = best_fit_model[m_srcname]; 361 + GModel* testsource = best_fit_model[m_srcname];
338 362
339 // Assign values to the maps 363 // Assign values to the maps
340 m_tsmap(i) = ts; 364 m_tsmap(i) = ts;
341 365
342 - // Get fit values  
343 - for (int j = 0; j < m_mapnames.size(); j++) {  
344 -  
345 - // get best fit values and assign to maps  
346 - m_maps[j](i) = (*testsource)[m_mapnames[j]].value(); 366 + // Extract fitted test source parameters
  367 + if (m_testsource != NULL) {
  368 + for (int j = 0; j < m_mapnames.size(); ++j) {
  369 + m_maps[j](i) = (*testsource)[m_mapnames[j]].value();
  370 + }
347 } 371 }
348 372
349 // Set status of bin to true 373 // Set status of bin to true
@@ -352,7 +376,7 @@ void cttsmap::run(void) @@ -352,7 +376,7 @@ void cttsmap::run(void)
352 // Remove model from container 376 // Remove model from container
353 models.remove(m_srcname); 377 models.remove(m_srcname);
354 378
355 - } // endfor: looped over bins 379 + } // endfor: looped over grid positions
356 380
357 // Bring models to initial state 381 // Bring models to initial state
358 m_obs.models(models_orig); 382 m_obs.models(models_orig);
@@ -372,12 +396,7 @@ void cttsmap::save(void) @@ -372,12 +396,7 @@ void cttsmap::save(void)
372 // Write header 396 // Write header
373 if (logTerse()) { 397 if (logTerse()) {
374 log << std::endl; 398 log << std::endl;
375 - if (m_obs.size() > 1) {  
376 - log.header1("Save observations");  
377 - }  
378 - else {  
379 - log.header1("Save observation");  
380 - } 399 + log.header1("Save TS map");
381 } 400 }
382 401
383 // Get output filename 402 // Get output filename
@@ -388,11 +407,11 @@ void cttsmap::save(void) @@ -388,11 +407,11 @@ void cttsmap::save(void)
388 407
389 // Write the sky maps to the FITS file 408 // Write the sky maps to the FITS file
390 m_tsmap.write(fitsfile); 409 m_tsmap.write(fitsfile);
391 -  
392 for (int i = 0; i < m_mapnames.size(); i++) { 410 for (int i = 0; i < m_mapnames.size(); i++) {
393 m_maps[i].write(fitsfile); 411 m_maps[i].write(fitsfile);
394 } 412 }
395 413
  414 + // Set extension name for all maps
396 for (int i = 0; i < m_mapnames.size(); i++) { 415 for (int i = 0; i < m_mapnames.size(); i++) {
397 fitsfile[i+1]->extname(m_mapnames[i]); 416 fitsfile[i+1]->extname(m_mapnames[i]);
398 } 417 }
@@ -423,7 +442,9 @@ void cttsmap::save(void) @@ -423,7 +442,9 @@ void cttsmap::save(void)
423 void cttsmap::init_members(void) 442 void cttsmap::init_members(void)
424 { 443 {
425 // Initialise members 444 // Initialise members
426 - m_evfile.clear(); 445 + m_infile.clear();
  446 + m_caldb.clear();
  447 + m_irf.clear();
427 m_outfile.clear(); 448 m_outfile.clear();
428 m_proj.clear(); 449 m_proj.clear();
429 m_coordsys.clear(); 450 m_coordsys.clear();
@@ -436,14 +457,14 @@ void cttsmap::init_members(void) @@ -436,14 +457,14 @@ void cttsmap::init_members(void)
436 // Initialise protected members 457 // Initialise protected members
437 m_obs.clear(); 458 m_obs.clear();
438 m_read_ahead = false; 459 m_read_ahead = false;
439 - m_binmin = -1;  
440 - m_binmax = -1;  
441 - m_logL0 = 0.0; 460 + m_binmin = -1;
  461 + m_binmax = -1;
  462 + m_logL0 = 0.0;
442 m_tsmap.clear(); 463 m_tsmap.clear();
443 - m_maps.clear();  
444 - m_mapnames.clear();  
445 m_statusmap.clear(); 464 m_statusmap.clear();
446 - m_testsource.clear(); 465 + m_mapnames.clear();
  466 + m_maps.clear();
  467 + m_testsource = NULL;
447 468
448 // Set logger properties 469 // Set logger properties
449 log.date(true); 470 log.date(true);
@@ -461,7 +482,9 @@ void cttsmap::init_members(void) @@ -461,7 +482,9 @@ void cttsmap::init_members(void)
461 void cttsmap::copy_members(const cttsmap& app) 482 void cttsmap::copy_members(const cttsmap& app)
462 { 483 {
463 // Copy attributes 484 // Copy attributes
464 - m_evfile = app.m_evfile; 485 + m_infile = app.m_infile;
  486 + m_caldb = app.m_caldb;
  487 + m_irf = app.m_irf;
465 m_outfile = app.m_outfile; 488 m_outfile = app.m_outfile;
466 m_proj = app.m_proj; 489 m_proj = app.m_proj;
467 m_coordsys = app.m_coordsys; 490 m_coordsys = app.m_coordsys;
@@ -472,16 +495,19 @@ void cttsmap::copy_members(const cttsmap&amp; app) @@ -472,16 +495,19 @@ void cttsmap::copy_members(const cttsmap&amp; app)
472 m_nypix = app.m_nypix; 495 m_nypix = app.m_nypix;
473 496
474 // Copy protected members 497 // Copy protected members
475 - m_binmin = app.m_binmin;  
476 - m_binmax = app.m_binmax;  
477 - m_logL0 = app.m_logL0; 498 + m_binmin = app.m_binmin;
  499 + m_binmax = app.m_binmax;
  500 + m_logL0 = app.m_logL0;
478 m_obs = app.m_obs; 501 m_obs = app.m_obs;
479 m_read_ahead = app.m_read_ahead; 502 m_read_ahead = app.m_read_ahead;
480 - m_tsmap = app.m_tsmap;  
481 - m_mapnames = app.m_mapnames; 503 + m_tsmap = app.m_tsmap;
  504 + m_mapnames = app.m_mapnames;
482 m_maps = app.m_maps; 505 m_maps = app.m_maps;
483 - m_statusmap = app.m_statusmap;  
484 - m_testsource = app.m_testsource; 506 + m_statusmap = app.m_statusmap;
  507 +
  508 + // Clone protected members
  509 + m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL;
  510 +
485 // Return 511 // Return
486 return; 512 return;
487 } 513 }
@@ -505,6 +531,10 @@ void cttsmap::free_members(void) @@ -505,6 +531,10 @@ void cttsmap::free_members(void)
505 /***********************************************************************//** 531 /***********************************************************************//**
506 * @brief Get application parameters 532 * @brief Get application parameters
507 * 533 *
  534 + * @exception GException::invalid_value
  535 + * Test source not found or no RA/DEC parameters found for test
  536 + * source.
  537 + *
508 * Get all task parameters from parameter file or (if required) by querying 538 * Get all task parameters from parameter file or (if required) by querying
509 * the user. Most parameters are only required if no observation exists so 539 * the user. Most parameters are only required if no observation exists so
510 * far in the observation container. In this case, a single CTA observation 540 * far in the observation container. In this case, a single CTA observation
@@ -517,19 +547,38 @@ void cttsmap::get_parameters(void) @@ -517,19 +547,38 @@ void cttsmap::get_parameters(void)
517 // observation using the parameters from the parameter file 547 // observation using the parameters from the parameter file
518 if (m_obs.size() == 0) { 548 if (m_obs.size() == 0) {
519 549
520 - // Get name of CTA events file  
521 - m_evfile = (*this)["evfile"].filename();  
522 -  
523 // Allocate CTA observation 550 // Allocate CTA observation
524 GCTAObservation obs; 551 GCTAObservation obs;
525 552
  553 + // Get event file name
  554 + std::string filename = (*this)["infile"].filename();
  555 +
526 // Try first to open as FITS file 556 // Try first to open as FITS file
527 try { 557 try {
528 558
529 - // Load event list in CTA observation  
530 - obs.load(m_evfile);  
531 -  
532 - // Append CTA observation to container 559 + // Load data
  560 + obs.load(filename);
  561 +
  562 + // Get other task parameters
  563 + m_caldb = (*this)["caldb"].string();
  564 + m_irf = (*this)["irf"].string();
  565 +
  566 + // Set calibration database. If specified parameter is a
  567 + // directory then use this as the pathname to the calibration
  568 + // database. Otherwise interpret this as the instrument name,
  569 + // the mission being "cta"
  570 + GCaldb caldb;
  571 + if (gammalib::dir_exists(m_caldb)) {
  572 + caldb.rootdir(m_caldb);
  573 + }
  574 + else {
  575 + caldb.open("cta", m_caldb);
  576 + }
  577 +
  578 + // Set reponse
  579 + obs.response(m_irf, caldb);
  580 +
  581 + // Append observation to container
533 m_obs.append(obs); 582 m_obs.append(obs);
534 583
535 } 584 }
@@ -538,38 +587,95 @@ void cttsmap::get_parameters(void) @@ -538,38 +587,95 @@ void cttsmap::get_parameters(void)
538 catch (GException::fits_open_error &e) { 587 catch (GException::fits_open_error &e) {
539 588
540 // Load observations from XML file 589 // Load observations from XML file
541 - m_obs.load(m_evfile); 590 + m_obs.load(filename);
542 591
543 - } 592 + // Check if all observations have response information. If
  593 + // not, get the calibration database parameters and set
  594 + // the response properly
  595 + bool asked_for_response = false;
  596 + for (int i = 0; i < m_obs.size(); ++i) {
  597 +
  598 + // Get CTA observation
  599 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
544 600
545 - m_obs.models((*this)["srcmdl"].filename()); 601 + // Continue only if observation is a CTA observation
  602 + if (obs != NULL) {
  603 +
  604 + // If response is not valid then set response from
  605 + // task parameters
  606 + if (!obs->has_response()) {
  607 +
  608 + // Get calibration parameters
  609 + if (!asked_for_response) {
  610 + m_caldb = (*this)["caldb"].string();
  611 + m_irf = (*this)["irf"].string();
  612 + asked_for_response = true;
  613 + }
  614 +
  615 + // Set response
  616 + GCaldb caldb;
  617 + if (gammalib::dir_exists(m_caldb)) {
  618 + caldb.rootdir(m_caldb);
  619 + }
  620 + else {
  621 + caldb.open("cta", m_caldb);
  622 + }
  623 + obs->response(m_irf, caldb);
  624 +
  625 + } // endif: there was no response
  626 + } // endif: observation was a CTA observation
  627 +
  628 + } // endfor: looped over observations
  629 +
  630 + } // endcatch: file was an XML file
546 631
547 } // endif: there was no observation in the container 632 } // endif: there was no observation in the container
548 633
  634 + // If there is are no models associated with the observations then
  635 + // load now the model definition
  636 + if (m_obs.models().size() == 0) {
  637 +
  638 + // Get models XML filename
  639 + std::string filename = (*this)["srcmdl"].filename();
  640 +
  641 + // Setup models for optimizing.
  642 + m_obs.models(GModels(filename));
  643 +
  644 + } // endif: no models were associated with observations
  645 +
549 // Get name of test source and check container for this name 646 // Get name of test source and check container for this name
550 m_srcname = (*this)["srcname"].string(); 647 m_srcname = (*this)["srcname"].string();
551 if (!m_obs.models().contains(m_srcname)) { 648 if (!m_obs.models().contains(m_srcname)) {
552 - throw GException::model_not_found(G_GET_PARAMETERS, m_srcname); 649 + std::string msg = "Source \""+m_srcname+"\" not found in model "
  650 + "container. Please add a source with that name "
  651 + "or check for a possible typos.";
  652 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
553 } 653 }
554 654
555 // Get Model and store it as protected member 655 // Get Model and store it as protected member
556 - GModel* model = m_obs.models()[m_srcname]->clone();  
557 - m_testsource = *dynamic_cast<GModelSky*>(model);  
558 -  
559 - // Check if model has RA and DEC parameters which are necessary to continue  
560 - if (!m_testsource.has_par("RA") || !m_testsource.has_par("DEC")) {  
561 -  
562 - // Is this the correct exception? Do we need a new one?  
563 - throw GException::model_invalid(G_GET_PARAMETERS, m_srcname, "RA and DEC parameters required! "); 656 + m_testsource = m_obs.models()[m_srcname]->clone();
  657 +
  658 + // Check if model has RA and DEC parameters which are necessary to
  659 + // continue
  660 + if (m_testsource != NULL) {
  661 + if (!m_testsource->has_par("RA") || !m_testsource->has_par("DEC")) {
  662 + std::string msg = "Source \""+m_srcname+"\" has no \"RA\" and "
  663 + "\"DEC\" parameters. Only sources with \"RA\" "
  664 + " and \"DEC\" parameters can be used as test "
  665 + "sources.";
  666 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  667 + }
564 } 668 }
565 669
566 - // Fix the spatial parameters  
567 - m_testsource["RA"].fix();  
568 - m_testsource["DEC"].fix(); 670 + // Fix the spatial parameters of the test source
  671 + if (m_testsource != NULL) {
  672 + (*m_testsource)["RA"].fix();
  673 + (*m_testsource)["DEC"].fix();
  674 + }
569 675
570 // Get map parameters 676 // Get map parameters
571 - m_xref = (*this)["xref"].real();  
572 - m_yref = (*this)["yref"].real(); 677 + m_xref = (*this)["xref"].real();
  678 + m_yref = (*this)["yref"].real();
573 m_proj = (*this)["proj"].string(); 679 m_proj = (*this)["proj"].string();
574 m_coordsys = (*this)["coordsys"].string(); 680 m_coordsys = (*this)["coordsys"].string();
575 m_binsz = (*this)["binsz"].real(); 681 m_binsz = (*this)["binsz"].real();
@@ -579,7 +685,7 @@ void cttsmap::get_parameters(void) @@ -579,7 +685,7 @@ void cttsmap::get_parameters(void)
579 // Get optional splitting parameters 685 // Get optional splitting parameters
580 m_binmin = (*this)["binmin"].integer(); 686 m_binmin = (*this)["binmin"].integer();
581 m_binmax = (*this)["binmax"].integer(); 687 m_binmax = (*this)["binmax"].integer();
582 - m_logL0 = (*this)["logL0"].real(); 688 + m_logL0 = (*this)["logL0"].real();
583 689
584 // Optionally read ahead parameters so that they get correctly 690 // Optionally read ahead parameters so that they get correctly
585 // dumped into the log file 691 // dumped into the log file
@@ -593,9 +699,9 @@ void cttsmap::get_parameters(void) @@ -593,9 +699,9 @@ void cttsmap::get_parameters(void)
593 699
594 700
595 /***********************************************************************//** 701 /***********************************************************************//**
596 - * @brief Initialise map information 702 + * @brief Initialise skymaps
597 * 703 *
598 - * Initialises the skymaps. 704 + * Initialises skymaps that will contain map information.
599 ***************************************************************************/ 705 ***************************************************************************/
600 void cttsmap::init_maps(void) 706 void cttsmap::init_maps(void)
601 { 707 {
@@ -604,34 +710,39 @@ void cttsmap::init_maps(void) @@ -604,34 +710,39 @@ void cttsmap::init_maps(void)
604 710
605 // Create skymap 711 // Create skymap
606 m_tsmap = GSkymap(m_proj, m_coordsys, 712 m_tsmap = GSkymap(m_proj, m_coordsys,
607 - m_xref, m_yref, -m_binsz, m_binsz,  
608 - m_nxpix, m_nypix, 1); 713 + m_xref, m_yref, -m_binsz, m_binsz,
  714 + m_nxpix, m_nypix, 1);
609 715
610 // Initialise map information 716 // Initialise map information
611 m_statusmap.clear(); 717 m_statusmap.clear();
612 718
613 // Create status map 719 // Create status map
614 m_statusmap = GSkymap(m_proj, m_coordsys, 720 m_statusmap = GSkymap(m_proj, m_coordsys,
615 - m_xref, m_yref, -m_binsz, m_binsz,  
616 - m_nxpix, m_nypix, 1);  
617 -  
618 - // Initialise maps of other free parameters  
619 - for (int i = 0; i < m_testsource.size(); i++) {  
620 - if (m_testsource[i].is_fixed()) {  
621 - continue;  
622 - }  
623 -  
624 - // create sky map  
625 - GSkymap map = GSkymap(m_proj, m_coordsys,  
626 - m_xref, m_yref, -m_binsz, m_binsz,  
627 - m_nxpix, m_nypix, 1);  
628 - m_maps.push_back(map);  
629 -  
630 - // Store name of parameter  
631 - m_mapnames.push_back(m_testsource[i].name());  
632 - } 721 + m_xref, m_yref, -m_binsz, m_binsz,
  722 + m_nxpix, m_nypix, 1);
  723 +
  724 + // Initialise maps of free parameters
  725 + if (m_testsource != NULL) {
  726 +
  727 + // Loop over all model parameters
  728 + for (int i = 0; i < m_testsource->size(); i++) {
  729 +
  730 + // Skip fixed parameters
  731 + if ((*m_testsource)[i].is_fixed()) {
  732 + continue;
  733 + }
  734 +
  735 + // Add sky map for free parameters. Note that push_back will
  736 + // create a copy.
  737 + m_maps.push_back(m_tsmap);
  738 +
  739 + // Store parameter name
  740 + m_mapnames.push_back((*m_testsource)[i].name());
  741 +
  742 + } // endfor: looped over all model parameters
  743 +
  744 + } // endif: test source was valid
633 745
634 // Return 746 // Return
635 return; 747 return;
636 } 748 }
637 -  
src/cttsmap/cttsmap.hpp
1 /*************************************************************************** 1 /***************************************************************************
2 - * cttsmap - TS map tool * 2 + * cttsmap - TS map calculation tool *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 - * copyright (C) 2010-2014 by Juergen Knoedlseder * 4 + * copyright (C) 2014 by Michael Mayer *
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 *
@@ -20,7 +20,7 @@ @@ -20,7 +20,7 @@
20 ***************************************************************************/ 20 ***************************************************************************/
21 /** 21 /**
22 * @file cttsmap.hpp 22 * @file cttsmap.hpp
23 - * @brief TS map calculation tool 23 + * @brief TS map calculation tool interface definition
24 * @author Michael Mayer 24 * @author Michael Mayer
25 */ 25 */
26 26
@@ -39,19 +39,20 @@ @@ -39,19 +39,20 @@
39 /***********************************************************************//** 39 /***********************************************************************//**
40 * @class cttsmap 40 * @class cttsmap
41 * 41 *
42 - * @brief TS map calculation interface 42 + * @brief TS map calculation
43 * 43 *
44 * This class computes a set of maps from any kind of observations: 44 * This class computes a set of maps from any kind of observations:
45 - * - TS map  
46 - * - flux map  
47 - * - map of spectral index  
48 - * The class operates on predefined observation containers, an individual  
49 - * event list or an observation definition XML file.  
50 * 45 *
51 - * During the computation a putative point-like source is moved along a  
52 - * grid of coordinates. The best fit results (TS, flux, index) are stored in maps  
53 - * which are saved in the output FITS files. 46 + * - TS map
  47 + * - Flux map
  48 + * - Map of spectral index
54 * 49 *
  50 + * The class operates on predefined observation containers, an individual
  51 + * event list or an observation definition XML file.
  52 + *
  53 + * During the computation a putative point-like source is moved along a
  54 + * grid of coordinates. The best fit results (TS, flux, index) are stored in
  55 + * maps which are saved in the output FITS files.
55 ***************************************************************************/ 56 ***************************************************************************/
56 class cttsmap : public GApplication { 57 class cttsmap : public GApplication {
57 public: 58 public:
@@ -71,7 +72,7 @@ public: @@ -71,7 +72,7 @@ public:
71 void run(void); 72 void run(void);
72 void save(void); 73 void save(void);
73 const GObservations& obs(void) const; 74 const GObservations& obs(void) const;
74 - const GSkymap& tsmap(void) const; 75 + const GSkymap& tsmap(void) const;
75 76
76 protected: 77 protected:
77 // Protected methods 78 // Protected methods
@@ -82,8 +83,10 @@ protected: @@ -82,8 +83,10 @@ protected:
82 void init_maps(void); 83 void init_maps(void);
83 84
84 // User parameters 85 // User parameters
85 - std::string m_evfile; //!< Input event list or XML file  
86 - std::string m_srcname; //!< Name of source which is moved around 86 + std::string m_infile; //!< Input file
  87 + std::string m_srcname; //!< Name of source which is moved around
  88 + std::string m_caldb; //!< Calibration database
  89 + std::string m_irf; //!< Instrument response functions
87 std::string m_outfile; //!< Output counts map or XML file 90 std::string m_outfile; //!< Output counts map or XML file
88 std::string m_proj; //!< WCS projection 91 std::string m_proj; //!< WCS projection
89 std::string m_coordsys; //!< Coordinate system 92 std::string m_coordsys; //!< Coordinate system
@@ -96,19 +99,19 @@ protected: @@ -96,19 +99,19 @@ protected:
96 // Parameters to control speed and job splitting 99 // Parameters to control speed and job splitting
97 int m_binmin; //!< Map bin number from which computation should start 100 int m_binmin; //!< Map bin number from which computation should start
98 int m_binmax; //!< Map bin number where map computation should end 101 int m_binmax; //!< Map bin number where map computation should end
99 - double m_logL0; //!< Likelihood value of null hypothesis 102 + double m_logL0; //!< Likelihood value of null hypothesis
100 103
101 // Protected members 104 // Protected members
102 GObservations m_obs; //!< Observation container 105 GObservations m_obs; //!< Observation container
103 bool m_read_ahead; //!< Read ahead parameters 106 bool m_read_ahead; //!< Read ahead parameters
104 - GSkymap m_tsmap; //!< ts map  
105 - GSkymap m_statusmap; //!< map of computed bins  
106 - std::vector<std::string> m_mapnames; // names of free parameters  
107 - std::vector<GSkymap> m_maps; // sky maps for each free parameter  
108 - GModelSky m_testsource; //!< test source for TS computation  
109 - 107 + GSkymap m_tsmap; //!< TS map
  108 + GSkymap m_statusmap; //!< Map of computed bins
  109 + std::vector<std::string> m_mapnames; //!< Names of free parameters
  110 + std::vector<GSkymap> m_maps; //!< Sky maps for each free parameter
  111 + GModel* m_testsource; //!< Pointer to test source for TS computation
110 }; 112 };
111 113
  114 +
112 /***********************************************************************//** 115 /***********************************************************************//**
113 * @brief Return observation container 116 * @brief Return observation container
114 * 117 *
@@ -120,6 +123,7 @@ const GObservations&amp; cttsmap::obs(void) const @@ -120,6 +123,7 @@ const GObservations&amp; cttsmap::obs(void) const
120 return m_obs; 123 return m_obs;
121 } 124 }
122 125
  126 +
123 /***********************************************************************//** 127 /***********************************************************************//**
124 * @brief Return TS skymap 128 * @brief Return TS skymap
125 * 129 *
@@ -131,5 +135,4 @@ const GSkymap&amp; cttsmap::tsmap(void) const @@ -131,5 +135,4 @@ const GSkymap&amp; cttsmap::tsmap(void) const
131 return m_tsmap; 135 return m_tsmap;
132 } 136 }
133 137
134 -  
135 #endif /* CTTSMAP_HPP */ 138 #endif /* CTTSMAP_HPP */
src/cttsmap/cttsmap.par
@@ -17,10 +17,12 @@ @@ -17,10 +17,12 @@
17 # 17 #
18 # General parameters 18 # General parameters
19 #=================== 19 #===================
20 -evfile, f, a, "events.fits",,,"Input event list or observation definition file"  
21 -srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"  
22 -srcname, s,a,"TestSource",,,"Test source"  
23 -outfile, f, a, "tsmap.fits",,, "Output counts cube" 20 +infile, f, a, "events.fits",,,"Input event list, cube or observation definition file"
  21 +srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  22 +srcname, s, a, "TestSource",,,"Test source"
  23 +caldb, s, a, "dummy",,, "Calibration database"
  24 +irf, s, a, "cta_dummy_irf",,, "Instrument response function"
  25 +outfile, f, a, "tsmap.fits",,, "Output counts cube"
24 26
25 # 27 #
26 # Spatial binning parameters 28 # Spatial binning parameters
@@ -36,7 +38,7 @@ proj, s, a, &quot;CAR&quot;, AIT|AZP|CAR|MER|STG|TAN,, &quot;Projection method e.g. AIT|AZP @@ -36,7 +38,7 @@ proj, s, a, &quot;CAR&quot;, AIT|AZP|CAR|MER|STG|TAN,, &quot;Projection method e.g. AIT|AZP
36 # 38 #
37 # Parameters for splitting and speed purpose 39 # Parameters for splitting and speed purpose
38 #=========================== 40 #===========================
39 -binmin, i, h, -1,,, "First bin to compute" 41 +binmin, i, h, -1,,, "First bin to compute"
40 binmax, i, h, -1,,, "Last bin to compute" 42 binmax, i, h, -1,,, "Last bin to compute"
41 logL0, r, h, 0.0,,, "LogLikelihood value of null hypothesis" 43 logL0, r, h, 0.0,,, "LogLikelihood value of null hypothesis"
42 44
@@ -48,4 +50,4 @@ chatter, i, h, 2,0,4, &quot;Chattiness of output&quot; @@ -48,4 +50,4 @@ chatter, i, h, 2,0,4, &quot;Chattiness of output&quot;
48 clobber, b, h, yes,,, "Overwrite existing output files with new output files?" 50 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
49 debug, b, h, no,,, "Debugging mode activated" 51 debug, b, h, no,,, "Debugging mode activated"
50 mode, s, h, "ql",,, "Mode of automatic parameters" 52 mode, s, h, "ql",,, "Mode of automatic parameters"
51 -logfile, f, h, "ctbin.log",,, "Log filename" 53 +logfile, f, h, "cttsmap.log",,, "Log filename"
src/cttsmap/main.cpp
1 /*************************************************************************** 1 /***************************************************************************
2 - * ctbin - CTA data binning main code * 2 + * cttsmap - TS map calculation tool *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 - * copyright (C) 2010-2011 by Jurgen Knodlseder * 4 + * copyright (C) 2014 by Michael Mayer *
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 *
@@ -19,9 +19,9 @@ @@ -19,9 +19,9 @@
19 * * 19 * *
20 ***************************************************************************/ 20 ***************************************************************************/
21 /** 21 /**
22 - * @file ctbin/main.cpp  
23 - * @brief CTA data binning tool main code  
24 - * @author J. Knodlseder 22 + * @file cttsmap/main.cpp
  23 + * @brief CTA TS map generation
  24 + * @author Michael Mayer
25 */ 25 */
26 26
27 /* __ Includes ___________________________________________________________ */ 27 /* __ Includes ___________________________________________________________ */