Luigi Tibaldo / ctools

Commit 230a31bb2a963669c414c59f336808e6085a3db7

Authored by Michael Mayer
1 parent d124c4ae

Changed cttsmap.par and code to support TS map creation from given xml models

Inline Side-by-side
Showing 3 changed files with 94 additions and 182 deletions
src/cttsmap/cttsmap.cpp
  View file @230a31b
... ... @@ -33,6 +33,7 @@
33 33 #include "GTools.hpp"
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
  36 +#define G_GET_PARAMETERS "cttsmap::get_parameters()"
36 37  
37 38 /* __ Debug definitions __________________________________________________ */
38 39  
... ... @@ -272,9 +273,6 @@ void cttsmap::run(void)
272 273 // Initialise maps
273 274 init_maps();
274 275  
275   - // Initialise point source model
276   - init_testsource();
277   -
278 276 // Get bins to be computed
279 277 // To split the computation on several jobs,
280 278 // the user is able to provide a subset of bins
... ... @@ -282,18 +280,28 @@ void cttsmap::run(void)
282 280 int binmin = (m_binmin == -1) ? 0 : m_binmin;
283 281 int binmax = (m_binmax == -1) ? m_tsmap.npix() : m_binmax;
284 282  
  283 + // Initialise optimizer
285 284 GOptimizerLM* opt = (logTerse()) ? new GOptimizerLM(log)
286 285 : new GOptimizerLM();
287 286  
288   - // Get likelihood of null hypothesis if not given before
289   - if (m_logL0 == 0.0) {
290   - m_obs.optimize(*opt);
291   - m_logL0 = -(opt->value());
292   - }
293   -
  287 + // Store initial models
294 288 GModels models_orig = m_obs.models();
  289 +
  290 + // Get model instance for further computations
295 291 GModels models = m_obs.models();
296 292  
  293 + // Remove test source
  294 + models.remove(m_srcname);
  295 +
  296 + // Get likelihood of null hypothesis if not given before
  297 + if (m_logL0 == 0.0) {
  298 +
  299 + // Compute likelihood without the test source
  300 + m_obs.models(models);
  301 + m_obs.optimize(*opt);
  302 + m_logL0 = -(opt->value());
  303 + }
  304 +
297 305 // Loop over all observations in the container
298 306 for (int i = binmin; i < binmax; ++i) {
299 307  
... ... @@ -309,10 +317,10 @@ void cttsmap::run(void)
309 317 m_testsource["RA"].value(bincentre.ra_deg());
310 318 m_testsource["DEC"].value(bincentre.dec_deg());
311 319  
312   - // append test source to model container
  320 + // Add test source to model container
313 321 models.append(m_testsource);
314 322  
315   - // Assign models to observation
  323 + // Assign models to observations
316 324 m_obs.models(models);
317 325  
318 326 // Optimize observation container
... ... @@ -326,22 +334,29 @@ void cttsmap::run(void)
326 334  
327 335 // Get test source model instance
328 336 GModels best_fit_model = m_obs.models();
329   - GModel* testsource = best_fit_model["TestSource"];
  337 + GModel* testsource = best_fit_model[m_srcname];
330 338  
331 339 // Assign values to the maps
332 340 m_tsmap(i) = ts;
333   - m_fluxmap(i) = (*testsource)["Integral"].value();
334   - m_indexmap(i) = (*testsource)["Index"].value();
335   - m_statusmap(i) = 1.0;
336 341  
337   - // Remove test source from container
338   - models.remove("TestSource");
  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();
  347 + }
  348 +
  349 + // Set status of bin to true
  350 + m_statusmap(i) = 1.0;
339 351  
340   - // Bring models to initial state
341   - m_obs.models(models_orig);
  352 + // Remove model from container
  353 + models.remove(m_srcname);
342 354  
343 355 } // endfor: looped over bins
344 356  
  357 + // Bring models to initial state
  358 + m_obs.models(models_orig);
  359 +
345 360 // Return
346 361 return;
347 362 }
... ... @@ -373,18 +388,19 @@ void cttsmap::save(void)
373 388  
374 389 // Write the sky maps to the FITS file
375 390 m_tsmap.write(fitsfile);
376   - m_fluxmap.write(fitsfile);
377   - m_indexmap.write(fitsfile);
378 391  
379   - // Change names of extensions
380   - fitsfile[0]->extname("TS MAP");
381   - fitsfile[1]->extname("FLUX MAP");
382   - fitsfile[2]->extname("SPECTRAL INDEX MAP");
  392 + for (int i = 0; i < m_mapnames.size(); i++) {
  393 + m_maps[i].write(fitsfile);
  394 + }
  395 +
  396 + for (int i = 0; i < m_mapnames.size(); i++) {
  397 + fitsfile[i+1]->extname(m_mapnames[i]);
  398 + }
383 399  
384 400 // Add computation log if not all bins are computed
385 401 if (m_binmin != -1 || m_binmax != -1) {
386 402 m_statusmap.write(fitsfile);
387   - fitsfile[3]->extname("STATUS MAP");
  403 + fitsfile[m_mapnames.size()+1]->extname("STATUS MAP");
388 404 }
389 405  
390 406 // Save FITS file
... ... @@ -424,19 +440,10 @@ void cttsmap::init_members(void)
424 440 m_binmax = -1;
425 441 m_logL0 = 0.0;
426 442 m_tsmap.clear();
427   - m_fluxmap.clear();
428   - m_indexmap.clear();
  443 + m_maps.clear();
  444 + m_mapnames.clear();
429 445 m_statusmap.clear();
430 446 m_testsource.clear();
431   - m_source_integral = 1e-12;
432   - m_source_integral_min = 1e-16;
433   - m_source_integral_max = 1e-8;
434   - m_source_index = -2.5;
435   - m_source_index_min = - 5.0;
436   - m_source_index_max = -1.0;
437   - m_source_emin = GEnergy(1,"TeV");
438   - m_source_emax = GEnergy(100, "TeV");
439   - m_free_index = false;
440 447  
441 448 // Set logger properties
442 449 log.date(true);
... ... @@ -471,20 +478,10 @@ void cttsmap::copy_members(const cttsmap&amp; app)
471 478 m_obs = app.m_obs;
472 479 m_read_ahead = app.m_read_ahead;
473 480 m_tsmap = app.m_tsmap;
474   - m_fluxmap = app.m_fluxmap;
475   - m_indexmap = app.m_indexmap;
  481 + m_mapnames = app.m_mapnames;
  482 + m_maps = app.m_maps;
476 483 m_statusmap = app.m_statusmap;
477   - m_testsource = app.m_testsource;
478   - m_source_integral = app.m_source_integral;
479   - m_source_integral_min = app.m_source_integral_min;
480   - m_source_integral_max = app.m_source_integral_max;
481   - m_source_index = app.m_source_integral_max;
482   - m_source_index_min = app.m_source_index_min;
483   - m_source_index_max = app.m_source_index_max;
484   - m_source_emin = app.m_source_emin;
485   - m_source_emax = app.m_source_emax;
486   - m_free_index = app.m_free_index;
487   -
  484 + m_testsource = app.m_testsource;
488 485 // Return
489 486 return;
490 487 }
... ... @@ -549,6 +546,27 @@ void cttsmap::get_parameters(void)
549 546  
550 547 } // endif: there was no observation in the container
551 548  
  549 + // Get name of test source and check container for this name
  550 + m_srcname = (*this)["srcname"].string();
  551 + if (!m_obs.models().contains(m_srcname)) {
  552 + throw GException::model_not_found(G_GET_PARAMETERS, m_srcname);
  553 + }
  554 +
  555 + // 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! ");
  564 + }
  565 +
  566 + // Fix the spatial parameters
  567 + m_testsource["RA"].fix();
  568 + m_testsource["DEC"].fix();
  569 +
552 570 // Get map parameters
553 571 m_xref = (*this)["xref"].real();
554 572 m_yref = (*this)["yref"].real();
... ... @@ -558,17 +576,6 @@ void cttsmap::get_parameters(void)
558 576 m_nxpix = (*this)["nxpix"].integer();
559 577 m_nypix = (*this)["nypix"].integer();
560 578  
561   - // Get test source parameters
562   - m_source_integral = (*this)["source_integral"].real();
563   - m_source_integral_min = (*this)["source_integral_min"].real();
564   - m_source_integral_max = (*this)["source_integral_max"].real();
565   - m_source_index = (*this)["source_index"].real();
566   - m_source_index_min = (*this)["source_index_min"].real();
567   - m_source_index_max = (*this)["source_index_max"].real();
568   - m_source_emin = GEnergy((*this)["source_emin"].real(),"TeV");
569   - m_source_emax = GEnergy((*this)["source_emax"].real(),"TeV");
570   - m_free_index = (*this)["free_index"].boolean();
571   -
572 579 // Get optional splitting parameters
573 580 m_binmin = (*this)["binmin"].integer();
574 581 m_binmax = (*this)["binmax"].integer();
... ... @@ -601,65 +608,27 @@ void cttsmap::init_maps(void)
601 608 m_nxpix, m_nypix, 1);
602 609  
603 610 // Initialise map information
604   - m_fluxmap.clear();
605   -
606   - // Create skymap
607   - m_fluxmap = GSkymap(m_proj, m_coordsys,
608   - m_xref, m_yref, -m_binsz, m_binsz,
609   - m_nxpix, m_nypix, 1);
610   -
611   - // Initialise map information
612   - m_indexmap.clear();
613   -
614   - // Create skymap
615   - m_indexmap = GSkymap(m_proj, m_coordsys,
616   - m_xref, m_yref, -m_binsz, m_binsz,
617   - m_nxpix, m_nypix, 1);
618   -
619   - // Initialise map information
620   - m_statusmap.clear();
621   -
622   - // Create status map
623   - m_statusmap = GSkymap(m_proj, m_coordsys,
624   - m_xref, m_yref, -m_binsz, m_binsz,
625   - m_nxpix, m_nypix, 1);
626   -
627   - // Return
628   - return;
629   -}
630   -
631   -
632   -/***********************************************************************//**
633   - * @brief Initialise the test source
634   - *
635   - * Initialises the test point source model.
636   - ***************************************************************************/
637   -void cttsmap::init_testsource(void)
638   -{
639   - // Create spatial model
640   - GModelSpatialPointSource spat = GModelSpatialPointSource(m_tsmap.inx2dir(0));
641   - spat["RA"].fix();
642   - spat["DEC"].fix();
643   -
644   - // Create spectral model
645   - GModelSpectralPlaw2 plaw2 = GModelSpectralPlaw2(m_source_integral,m_source_index,m_source_emin,m_source_emax);
646   -
647   - // Create sky model
648   - m_testsource = GModelSky(spat,plaw2);
649   - m_testsource.name("TestSource");
650   -
651   - // Set source parameter boundaries
652   - m_testsource["Integral"].min(m_source_integral_min);
653   - m_testsource["Integral"].max(m_source_integral_max);
654   - m_testsource["Index"].min(m_source_index_min);
655   - m_testsource["Index"].max(m_source_index_max);
656   -
657   - // Free spectral index if required
658   - if (m_free_index) {
659   - m_testsource["Index"].free();
660   - }
661   - else {
662   - m_testsource["Index"].fix();
  611 + m_statusmap.clear();
  612 +
  613 + // Create status map
  614 + 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());
663 632 }
664 633  
665 634 // Return
... ...
src/cttsmap/cttsmap.hpp
  View file @230a31b
... ... @@ -72,9 +72,6 @@ public:
72 72 void save(void);
73 73 const GObservations& obs(void) const;
74 74 const GSkymap& tsmap(void) const;
75   - const GSkymap& fluxmap(void) const;
76   - const GSkymap& indexmap(void) const;
77   - const GSkymap& statusmap(void) const;
78 75  
79 76 protected:
80 77 // Protected methods
... ... @@ -83,10 +80,10 @@ protected:
83 80 void free_members(void);
84 81 void get_parameters(void);
85 82 void init_maps(void);
86   - void init_testsource(void);
87 83  
88 84 // User parameters
89 85 std::string m_evfile; //!< Input event list or XML file
  86 + std::string m_srcname; //!< Name of source which is moved around
90 87 std::string m_outfile; //!< Output counts map or XML file
91 88 std::string m_proj; //!< WCS projection
92 89 std::string m_coordsys; //!< Coordinate system
... ... @@ -104,20 +101,11 @@ protected:
104 101 // Protected members
105 102 GObservations m_obs; //!< Observation container
106 103 bool m_read_ahead; //!< Read ahead parameters
107   - GSkymap m_tsmap; //!< ts map
108   - GSkymap m_fluxmap; //!< flux map
109   - GSkymap m_indexmap; //!< spectral index map
  104 + GSkymap m_tsmap; //!< ts map
110 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
111 108 GModelSky m_testsource; //!< test source for TS computation
112   - double m_source_integral; //!< Integral flux of test source
113   - double m_source_integral_min; //!< minimum boundary of integral flux
114   - double m_source_integral_max; //!< maximum boundary of integral flux
115   - double m_source_index; //!< assumed spectral index of test source
116   - double m_source_index_min; //!< minimum boundary of spectral index of test source
117   - double m_source_index_max; //!< maximum boundary of spectral index of test source
118   - GEnergy m_source_emin; //!< Minimum energy for source integral
119   - GEnergy m_source_emax; //!< Maximum energy for source integral
120   - bool m_free_index; //!< flag if spectral index of test source is left free or fixed
121 109  
122 110 };
123 111  
... ... @@ -143,38 +131,5 @@ const GSkymap&amp; cttsmap::tsmap(void) const
143 131 return m_tsmap;
144 132 }
145 133  
146   -/***********************************************************************//**
147   - * @brief Return flux skymap
148   - *
149   - * @return Reference to flux skymap
150   - ***************************************************************************/
151   -inline
152   -const GSkymap& cttsmap::fluxmap(void) const
153   -{
154   - return m_fluxmap;
155   -}
156   -
157   -/***********************************************************************//**
158   - * @brief Return spectral index skymap
159   - *
160   - * @return Reference to spectral index skymap
161   - ***************************************************************************/
162   -inline
163   -const GSkymap& cttsmap::indexmap(void) const
164   -{
165   - return m_indexmap;
166   -}
167   -
168   -/***********************************************************************//**
169   - * @brief Return status skymap
170   - *
171   - * @return Reference to status skymap
172   - ***************************************************************************/
173   -inline
174   -const GSkymap& cttsmap::statusmap(void) const
175   -{
176   - return m_statusmap;
177   -}
178   -
179 134  
180 135 #endif /* CTTSMAP_HPP */
... ...
src/cttsmap/cttsmap.par
  View file @230a31b
... ... @@ -19,6 +19,7 @@
19 19 #===================
20 20 evfile, f, a, "events.fits",,,"Input event list or observation definition file"
21 21 srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  22 +srcname, s,a,"TestSource",,,"Test source"
22 23 outfile, f, a, "tsmap.fits",,, "Output counts cube"
23 24  
24 25 #
... ... @@ -33,19 +34,6 @@ yref, r, a, 22.01,,, &quot;Second coordinate of image center in degrees (DEC or g
33 34 proj, s, a, "CAR", AIT|AZP|CAR|MER|STG|TAN,, "Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"
34 35  
35 36 #
36   -# Test source parameters
37   -#===========================
38   -source_integral, r, h, 1e-10,,, "Integral start parameter [cm-2 s-1]"
39   -source_integral_min, r, h, 1e-16,,, "Integral minimum boundary [cm-2 s-1]"
40   -source_integral_max, r, h, 1e-3,,, "Integral maximum boundary [cm-2 s-1]"
41   -source_index, r, h, -2.5,,, "Index start parameter"
42   -source_index_min, r, h, -5.0,,, "Index minimum boundary"
43   -source_index_max, r, h, -1.0,,, "Index maximum boundary"
44   -source_emin, r, h, 1.0,,, "Minimum energy for flux computation"
45   -source_emax, r, h, 100.0,,, "Maximum energy for flux computation"
46   -free_index, b, h, no,,, "leave source index free?"
47   -
48   -#
49 37 # Parameters for splitting and speed purpose
50 38 #===========================
51 39 binmin, i, h, -1,,, "First bin to compute"
... ...