GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialDiffuseMap.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialDiffuseMap.cpp - Spatial map model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2016 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GModelSpatialDiffuseMap.cpp
23  * @brief Spatial map model class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GException.hpp"
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GHealpix.hpp"
37 
38 /* __ Constants __________________________________________________________ */
39 
40 /* __ Globals ____________________________________________________________ */
42 const GModelSpatialRegistry g_spatial_map_registry(&g_spatial_map_seed);
43 #if defined(G_LEGACY_XML_FORMAT)
44 const GModelSpatialDiffuseMap g_spatial_map_legacy_seed(true, "SpatialMap");
45 const GModelSpatialRegistry g_spatial_map_legacy_registry(&g_spatial_map_legacy_seed);
46 #endif
47 
48 /* __ Method name definitions ____________________________________________ */
49 #define G_MC "GModelSpatialDiffuseMap::mc(GEnergy&, GTime&, GRan&)"
50 #define G_READ "GModelSpatialDiffuseMap::read(GXmlElement&)"
51 #define G_WRITE "GModelSpatialDiffuseMap::write(GXmlElement&)"
52 
53 /* __ Macros _____________________________________________________________ */
54 
55 /* __ Coding definitions _________________________________________________ */
56 
57 /* __ Debug definitions __________________________________________________ */
58 //#define G_DEBUG_MC //!< Debug MC method
59 //#define G_DEBUG_MC_CACHE //!< Debug MC cache
60 
61 
62 /*==========================================================================
63  = =
64  = Constructors/destructors =
65  = =
66  ==========================================================================*/
67 
68 /***********************************************************************//**
69  * @brief Void constructor
70  *
71  * Constructs empty spatial map model.
72  ***************************************************************************/
75 {
76  // Initialise members
77  init_members();
78 
79  // Return
80  return;
81 }
82 
83 
84 /***********************************************************************//**
85  * @brief Model type constructor
86  *
87  * @param[in] dummy Dummy flag.
88  * @param[in] type Model type.
89  *
90  * Constructs empty spatial map model by specifying a model @p type.
91  ***************************************************************************/
93  const std::string& type) :
95 {
96  // Initialise members
97  init_members();
98 
99  // Set model type
100  m_type = type;
101 
102  // Return
103  return;
104 }
105 
106 
107 /***********************************************************************//**
108  * @brief XML constructor
109  *
110  * @param[in] xml XML element.
111  *
112  * Constructs spatial map model by extracting information from an XML
113  * element. See the read() method for more information about the expected
114  * structure of the XML element.
115  ***************************************************************************/
118 {
119  // Initialise members
120  init_members();
121 
122  // Read information from XML element
123  read(xml);
124 
125  // Return
126  return;
127 }
128 
129 
130 /***********************************************************************//**
131  * @brief File name constructor
132  *
133  * @param[in] filename File name.
134  * @param[in] value Normalization factor.
135  * @param[in] normalize Normalize map?
136  *
137  * Constructs spatial map model by loading a sky map from the file specified
138  * by @p filename and by setting the @p value by which the map will be
139  * multiplied (or normalized).
140  ***************************************************************************/
142  const double& value,
143  const bool& normalize) :
145 {
146  // Initialise members
147  init_members();
148 
149  // Set normalisation parameter
150  m_value.value(value);
151 
152  // Set normalization flag
154 
155  // Load sky map
156  load(filename);
157 
158  // Return
159  return;
160 }
161 
162 
163 /***********************************************************************//**
164  * @brief Sky map constructor
165  *
166  * @param[in] map Sky map.
167  * @param[in] value Normalization factor.
168  * @param[in] normalize Normalize map.
169  *
170  * Constructs spatial map model by setting the sky @p map and by setting the
171  * @p value by which the map will be multiplied (or normalized).
172  ***************************************************************************/
174  const double& value,
175  const bool& normalize) :
177 {
178  // Initialise members
179  init_members();
180 
181  // Set normalization parameter
182  m_value.value(value);
183 
184  // Set normalization flag
186 
187  // Set sky map
188  m_map = map;
189 
190  // Prepare sky map
191  prepare_map();
192 
193  // Return
194  return;
195 }
196 
197 
198 /***********************************************************************//**
199  * @brief Copy constructor
200  *
201  * @param[in] model Spatial map model.
202  *
203  * Construct a spatial map model by copying another spatial map model.
204  ***************************************************************************/
206  GModelSpatialDiffuse(model)
207 {
208  // Initialise members
209  init_members();
210 
211  // Copy members
212  copy_members(model);
213 
214  // Return
215  return;
216 }
217 
218 
219 /***********************************************************************//**
220  * @brief Destructor
221  ***************************************************************************/
223 {
224  // Free members
225  free_members();
226 
227  // Return
228  return;
229 }
230 
231 
232 /*==========================================================================
233  = =
234  = Operators =
235  = =
236  ==========================================================================*/
237 
238 /***********************************************************************//**
239  * @brief Assignment operator
240  *
241  * @param[in] model Spatial map model.
242  * @return Spatial map model.
243  *
244  * Assigns a spatial map model to another spatial map model.
245  ***************************************************************************/
247 {
248  // Execute only if object is not identical
249  if (this != &model) {
250 
251  // Copy base class members
252  this->GModelSpatialDiffuse::operator=(model);
253 
254  // Free members
255  free_members();
256 
257  // Initialise members
258  init_members();
259 
260  // Copy members
261  copy_members(model);
262 
263  } // endif: object was not identical
264 
265  // Return
266  return *this;
267 }
268 
269 
270 /*==========================================================================
271  = =
272  = Public methods =
273  = =
274  ==========================================================================*/
275 
276 /***********************************************************************//**
277  * @brief Clear spatial map model
278  ***************************************************************************/
280 {
281  // Free class members (base and derived classes, derived class first)
282  free_members();
285 
286  // Initialise members
289  init_members();
290 
291  // Return
292  return;
293 }
294 
295 
296 /***********************************************************************//**
297  * @brief Clone spatial map model
298  *
299  * @return Pointer to deep copy of spatial map model.
300  ***************************************************************************/
302 {
303  // Clone diffuse map
304  return new GModelSpatialDiffuseMap(*this);
305 }
306 
307 
308 /***********************************************************************//**
309  * @brief Return intensity of sky map
310  *
311  * @param[in] photon Incident photon.
312  * @param[in] gradients Compute gradients?
313  * @return Sky map intensity (\f$\mbox{ph cm}^{-2}\mbox{sr}^{-1}\mbox{s}^{-1}\f$)
314  *
315  * Returns the intensity of the sky map at the specified sky direction
316  * multiplied by the normalization factor. If the sky direction falls outside
317  * the sky map, an intensity of 0 is returned.
318  *
319  * If the @p gradients flag is true the method will also evaluate the partial
320  * derivatives of the model.
321  ***************************************************************************/
323  const bool& gradients) const
324 {
325  // Get skymap intensity
326  double intensity = m_map(photon.dir());
327 
328  // Optionally compute partial derivatives
329  if (gradients) {
330 
331  // Compute partial derivatives of the parameter values
332  double g_value = (m_value.is_free()) ? intensity * m_value.scale() : 0.0;
333 
334  // Set gradient
335  m_value.factor_gradient(g_value);
336 
337  } // endif: computed partial derivatives
338 
339  // Return intensity times normalization factor
340  return (intensity * m_value.value());
341 }
342 
343 
344 /***********************************************************************//**
345  * @brief Returns MC sky direction
346  *
347  * @param[in] energy Photon energy (ignored).
348  * @param[in] time Photon arrival time (ignored).
349  * @param[in,out] ran Random number generator.
350  * @return Sky direction.
351  *
352  * @exception GException::invalid_value
353  * No sky map defined.
354  * @exception GException::invalid_return_value
355  * Simulation cone not defined, does not overlap with map or map
356  * is empty.
357  *
358  * Returns a random sky direction according to the intensity distribution of
359  * the model sky map. The method uses a rejection method to determine the sky
360  * direction. If no sky direction could be determined, the method throws an
361  * GException::invalid_return_value exception.
362  ***************************************************************************/
364  const GTime& time,
365  GRan& ran) const
366 {
367  // Throw an exception if the maximum MC intensity is not positive. This
368  // can be the case because the simulation cone has not been defined or
369  // because it does not overlap with the sky map
370  if (m_mc_max <= 0.0) {
371  std::string msg = "Simulation cone has not been defined or does not "
372  "overlap with the sky map. Please specify a valid "
373  "simulation cone.";
375  }
376 
377  // Determine number of skymap pixels
378  int npix = m_map.npix();
379 
380  // Throw an exception if there are no sky map pixels
381  if (npix <= 0) {
382  std::string msg = "No sky map defined. Please specify a valid sky map.";
383  throw GException::invalid_value(G_MC, msg);
384  }
385 
386  // Allocate sky direction
387  GSkyDir dir;
388 
389  // Debug option: initialise counter
390  #if defined(G_DEBUG_MC)
391  int num_iterations = 0;
392  #endif
393 
394  // Get sky direction
395  while (true) {
396 
397  // Debug option: increment counter
398  #if defined(G_DEBUG_MC)
399  num_iterations++;
400  #endif
401 
402  // Simulate random sky direction within Monte Carlo simulation cone
403  double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
405  double phi = 360.0 * ran.uniform();
406  dir = m_mc_centre;
407  dir.rotate_deg(phi, theta);
408 
409  // Get map value at simulated sky direction. If the map value is non-
410  // positive then simulate a new sky direction.
411  double value = m_map(dir);
412  if (value <= 0.0) {
413  continue;
414  }
415 
416  // Get uniform random number
417  double uniform = ran.uniform() * m_mc_max;
418 
419  // Exit loop if the random number is not larger than the sky map value
420  if (uniform <= value) {
421  break;
422  }
423 
424  } // endwhile: loop until sky direction was accepted
425 
426  // Debug option: log counter
427  #if defined(G_DEBUG_MC)
428  std::cout << num_iterations << " ";
429  #endif
430 
431  // Return sky direction
432  return dir;
433 }
434 
435 
436 /***********************************************************************//**
437  * @brief Return normalization of diffuse map for Monte Carlo simulations
438  *
439  * @param[in] dir Centre of simulation cone.
440  * @param[in] radius Radius of simulation cone (degrees).
441  * @return Normalization.
442  *
443  * Returns the normalization of a diffuse map. The normalization is given
444  * by the model value times the integrated flux in the sky map over a
445  * circular region.
446  ***************************************************************************/
448  const double& radius) const
449 {
450  // Set the MC cone
451  set_mc_cone(dir, radius);
452 
453  // Retrieve normalization
454  double norm = m_mc_norm * value();
455 
456  // Return normalization
457  return norm;
458 }
459 
460 
461 /***********************************************************************//**
462  * @brief Signals whether model contains sky direction
463  *
464  * @param[in] dir Sky direction.
465  * @param[in] margin Margin to be added to sky direction (degrees)
466  * @return True if the model contains the sky direction.
467  *
468  * Signals whether a sky direction falls within the bounding circle of
469  * the diffuse map.
470  ***************************************************************************/
472  const double& margin) const
473 {
474  // Initialise containment flag
475  bool contains = false;
476 
477  // Continue only if radius is positive
478  if (m_radius > 0.0) {
479 
480  // Compute distance to centre
481  double distance = m_centre.dist_deg(dir);
482 
483  // If distance is smaller than radius plus margin we consider
484  // the position to be contained within the bounding circle
485  if (distance < m_radius + margin) {
486  contains = true;
487  }
488 
489  }
490 
491  // Return containment
492  return (contains);
493 }
494 
495 
496 /***********************************************************************//**
497  * @brief Read model from XML element
498  *
499  * @param[in] xml XML element.
500  *
501  * @exception GException::invalid_value
502  * Model parameters not found in XML element.
503  *
504  * Reads the spatial information for a diffuse map from an XML element. The
505  * expected format of the XML element is
506  *
507  * <spatialModel type="DiffuseMap" file="myfile.fits" normalize="1">
508  * <parameter name="Normalization" ../>
509  * </spatialModel>
510  *
511  * The @p file attribute provides the filename of the diffuse map FITS file.
512  * The filename may be either an absolute filename (starting with '/') or a
513  * relative filename. If no access path is given, the file is expected to
514  * reside in the same location as the XML file.
515  *
516  * The @p normalize attribute specifies whether the sky map should be
517  * normalized to unity flux or not. If the attribute is not given, the map
518  * will be automatically normalized. To prevent normalization,
519  * @p normalize="0" needs to be specified.
520  ***************************************************************************/
522 {
523  // If "Normalization" parameter exists then read parameter from this
524  // XML element
525  if (gammalib::xml_has_par(xml, "Normalization")) {
526  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Normalization");
527  m_value.read(*par);
528  }
529 
530  // ... otherwise try reading parameter from "Prefactor" parameter
531  #if defined(G_LEGACY_XML_FORMAT)
532  else if (gammalib::xml_has_par(xml, "Prefactor")) {
533  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Prefactor");
534  m_value.read(*par);
535  }
536  #endif
537 
538  // ... otherwise throw an exception
539  else {
540  #if defined(G_LEGACY_XML_FORMAT)
541  std::string msg = "Diffuse map model requires either "
542  "\"Normalization\" or \"Prefactor\" parameter.";
543  #else
544  std::string msg = "Diffuse map model requires \"Normalization\" "
545  "parameter.";
546  #endif
547  throw GException::invalid_value(G_READ, msg);
548  }
549 
550  // Get optional normalization attribute
551  m_normalize = true;
552  m_has_normalize = false;
553  if (xml.has_attribute("normalize")) {
554  m_has_normalize = true;
555  std::string arg = xml.attribute("normalize");
556  if (arg == "0" || gammalib::tolower(arg) == "false") {
557  m_normalize = false;
558  }
559  }
560 
561  // Load sky map.
562  load(gammalib::xml_file_expand(xml, xml.attribute("file")));
563 
564  // Return
565  return;
566 }
567 
568 
569 /***********************************************************************//**
570  * @brief Write model into XML element
571  *
572  * @param[in] xml XML element.
573  *
574  * @exception GException::model_invalid_spatial
575  * Existing XML element is not of type "SpatialMap"
576  * @exception GException::model_invalid_parnum
577  * Invalid number of model parameters found in XML element.
578  * @exception GException::model_invalid_parnames
579  * Invalid model parameter names found in XML element.
580  *
581  * Writes the spatial information for a diffuse map into an XML element. The
582  * format of the XML element is
583  *
584  * <spatialModel type="DiffuseMap" file="myfile.fits" normalize="1">
585  * <parameter name="Prefactor" value="1" min="0.1" max="10" scale="1" free="0"/>
586  * </spatialModel>
587  *
588  * The @p file attribute provides the filename of the diffuse map FITS file.
589  * The filename may be either an absolute filename (starting with '/') or a
590  * relative filename. If no access path is given, the file is expected to
591  * reside in the same location as the XML file.
592  *
593  * The @p normalize attribute specifies whether the sky map should be
594  * normalized to unity flux or not. The attribute will only be written if the
595  * normalization is disabled.
596  ***************************************************************************/
598 {
599  // Set model type
600  if (xml.attribute("type") == "") {
601  xml.attribute("type", type());
602  }
603 
604  // Verify model type
605  if (xml.attribute("type") != type()) {
607  "Spatial model is not of type \""+type()+"\".");
608  }
609 
610  // Set sky map file name
612 
613  // If XML element has 0 nodes then append parameter node. The name
614  // of the node is "Prefactor" as this is the Fermi/LAT standard.
615  // We thus assure that the XML files will be compatible with
616  // Fermi/LAT.
617  if (xml.elements() == 0) {
618  xml.append(GXmlElement("parameter name=\"Prefactor\""));
619  }
620 
621  // Verify that XML element has exactly 1 parameter
622  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
624  "Spatial map model requires exactly 1 parameter.");
625  }
626 
627  // Get pointer on model parameter
628  GXmlElement* par = xml.element("parameter", 0);
629 
630  // Set or update parameter
631  if (par->attribute("name") == "Normalization" ||
632  par->attribute("name") == "Prefactor") {
633  m_value.write(*par);
634  }
635  else {
637  "Spatial map model requires either \"Prefactor\" or"
638  " \"Normalization\" parameter.");
639  }
640 
641  // Set optional normalization attribute
642  if (m_has_normalize || !m_normalize) {
643  if (m_normalize) {
644  xml.attribute("normalize", "1");
645  }
646  else {
647  xml.attribute("normalize", "0");
648  }
649  }
650 
651  // Return
652  return;
653 }
654 
655 
656 /***********************************************************************//**
657  * @brief Set Monte Carlo simulation cone
658  *
659  * @param[in] centre Simulation cone centre.
660  * @param[in] radius Simulation cone radius (degrees).
661  *
662  * Sets the simulation cone centre and radius that defines the directions
663  * that will be simulated using the mc() method and pre-computes the maximum
664  * intensity and the spatially integrated flux of the map within the
665  * simulation cone region.
666  ***************************************************************************/
668  const double& radius) const
669 {
670  // Continue only if the simulation cone has changed
671  if ((centre != m_mc_centre) || (radius != m_mc_radius)) {
672 
673  // Save simulation cone definition
675  m_mc_radius = radius;
676 
677  // Pre-compute 1 - cosine of radius
679 
680  // Initialise map maximum and normalisation
681  m_mc_max = 0.0;
682  m_mc_norm = 0.0;
683 
684  // Determine number of skymap pixels
685  int npix = m_map.npix();
686 
687  // Continue only if there are pixels
688  if (npix > 0) {
689 
690  // Compute flux and maximum map intensity within the simulation cone
691  double sum = 0.0;
692  double sum_map = 0.0;
693  for (int i = 0; i < npix; ++i) {
694 
695  // Get map flux, intensity and distance from MC cone centre
696  double flux = m_map.flux(i);
697  double intensity = m_map(i);
698  double distance = centre.dist_deg(m_map.pix2dir(i));
699 
700  // Add flux if positive
701  if (flux > 0.0) {
702  if (distance <= radius) {
703  sum += flux; // flux within simulation cone
704  }
705  sum_map += flux; // total flux
706  }
707 
708  // Update maximum intensity
709  if (distance <= radius) {
710  if (intensity > m_mc_max) {
711  m_mc_max = intensity;
712  }
713  }
714 
715  } // endfor: looped over all map pixels
716 
717  // Set the normalization factor for the MC simulations. In case
718  // that the map is normalised, this is the fraction of the flux
719  // that is comprised within the simulation cone. For non-normalised
720  // maps, this is simply the flux comprised within the simulation
721  // cone.
722  if (sum_map > 0.0) {
723  if (m_normalize) {
724  m_mc_norm = sum / sum_map;
725  }
726  else {
727  m_mc_norm = sum;
728  }
729  }
730 
731  // Log maximum intensity and total flux for debugging
732  #if defined(G_DEBUG_MC_CACHE)
733  std::cout << "GModelSpatialDiffuseMap::set_mc_cone:" << std::endl;
734  std::cout << " Maximum map intensity:";
735  std::cout << m_mc_max << " ph/cm2/s/sr" << std::endl;
736  std::cout << " Spatially integrated flux:" << std::endl;
737  std::cout << sum << " ph/cm2/s" << std::endl;
738  std::cout << " Map normalisation:" << std::endl;
739  std::cout << m_mc_norm << " ph/cm2/s" << std::endl;
740  #endif
741 
742  } // endif: there were map pixels
743 
744  } // endif: simulation cone has changed
745 
746  // Return
747  return;
748 }
749 
750 
751 /***********************************************************************//**
752  * @brief Print map information
753  *
754  * @param[in] chatter Chattiness.
755  * @return String with diffuse map model information.
756  ***************************************************************************/
757 std::string GModelSpatialDiffuseMap::print(const GChatter& chatter) const
758 {
759  // Initialise result string
760  std::string result;
761 
762  // Continue only if chatter is not silent
763  if (chatter != SILENT) {
764 
765  // Append header
766  result.append("=== GModelSpatialDiffuseMap ===");
767 
768  // Append parameters
769  result.append("\n"+gammalib::parformat("Sky map file"));
770  result.append(m_filename);
771  result.append("\n"+gammalib::parformat("Map normalization"));
772  result.append(gammalib::str(m_mc_norm)+" ph/cm2/s");
773  if (normalize()) {
774  result.append(" [normalized]");
775  }
776  result.append("\n"+gammalib::parformat("Map centre"));
777  result.append(m_centre.print());
778  result.append("\n"+gammalib::parformat("Map radius"));
779  result.append(gammalib::str(m_radius)+" deg");
780  result.append("\n"+gammalib::parformat("Number of parameters"));
781  result.append(gammalib::str(size()));
782  for (int i = 0; i < size(); ++i) {
783  result.append("\n"+m_pars[i]->print(chatter));
784  }
785 
786  } // endif: chatter was not silent
787 
788  // Return result
789  return result;
790 }
791 
792 
793 /***********************************************************************//**
794  * @brief Load skymap into the model class
795  *
796  * @param[in] filename Sky map file.
797  *
798  * Loads skymap into the model class. The method calls the protected method
799  * prepare_map() that prepares the map for usage by the class.
800  *
801  * If the @p filename is empty, no map will be loaded.
802  ***************************************************************************/
804 {
805  // Initialise skymap
806  m_map.clear();
807 
808  // Continue only if filename is not empty
809  if (!filename.is_empty()) {
810 
811  // Store filename of skymap (for XML writing). Note that we do not
812  // expand any environment variable at this level, so that if we
813  // write back the XML element we write the filepath with the
814  // environment variables
816 
817  // Load sky map
819 
820  // Prepare sky map
821  prepare_map();
822 
823  } // endif: filename was not empty
824 
825  // Return
826  return;
827 }
828 
829 
830 /*==========================================================================
831  = =
832  = Private methods =
833  = =
834  ==========================================================================*/
835 
836 /***********************************************************************//**
837  * @brief Initialise class members
838  ***************************************************************************/
840 {
841  // Initialise model type
842  m_type = "DiffuseMap";
843 
844  // Initialise Value
845  m_value.clear();
846  m_value.name("Normalization");
847  m_value.value(1.0);
848  m_value.scale(1.0);
849  m_value.range(0.001, 1000.0);
850  m_value.gradient(0.0);
851  m_value.fix();
852  m_value.has_grad(true);
853 
854  // Set parameter pointer(s)
855  m_pars.clear();
856  m_pars.push_back(&m_value);
857 
858  // Initialise other members
859  m_map.clear();
860  m_filename.clear();
861  m_normalize = true;
862  m_has_normalize = false;
863  m_centre.clear();
864  m_radius = 0.0;
865  m_region.clear();
866 
867  // Initialise MC cache
868  m_mc_centre.clear();
869  m_mc_radius = -1.0; //!< Signal initialisation
870  m_mc_one_minus_cosrad = 1.0;
871  m_mc_norm = 0.0;
872  m_mc_max = 0.0;
873 
874  // Return
875  return;
876 }
877 
878 
879 /***********************************************************************//**
880  * @brief Copy class members
881  *
882  * @param[in] model Spatial map cube model.
883  ***************************************************************************/
885 {
886  // Copy members
887  m_type = model.m_type;
888  m_value = model.m_value;
889  m_map = model.m_map;
890  m_filename = model.m_filename;
891  m_normalize = model.m_normalize;
893  m_centre = model.m_centre;
894  m_radius = model.m_radius;
895  m_region = model.m_region;
896 
897  // Copy MC cache
898  m_mc_centre = model.m_mc_centre;
899  m_mc_radius = model.m_mc_radius;
901  m_mc_norm = model.m_mc_norm;
902  m_mc_max = model.m_mc_max;
903 
904  // Set parameter pointer(s)
905  m_pars.clear();
906  m_pars.push_back(&m_value);
907 
908  // Return
909  return;
910 }
911 
912 
913 /***********************************************************************//**
914  * @brief Delete class members
915  ***************************************************************************/
917 {
918  // Return
919  return;
920 }
921 
922 
923 /***********************************************************************//**
924  * @brief Prepare sky map after loading
925  *
926  * Prepares a sky map after loading. Negative, infinite or undefined skymap
927  * pixels are set to zero intensity. The method also determine the centre
928  * and radius of a circle enclosing the map, and the Monte Carlo simulation
929  * cone is set to this circle.
930  *
931  * If normalize() is true, the map is furthermore normalised so that the
932  * total flux in the map amounts to 1 ph/cm2/s. Negative skymap pixels are
933  * set to zero intensity.
934  ***************************************************************************/
936 {
937  // Initialise centre and radius
938  m_centre.clear();
939  m_radius = 0.0;
940 
941  // Determine number of skymap pixels
942  int npix = m_map.npix();
943 
944  // Continue only if there are skymap pixels
945  if (npix > 0) {
946 
947  // Initialise flux sum
948  double sum = 0.0;
949 
950  // Compute flux sum and set negative or invalid pixels to zero
951  // intensity.
952  for (int i = 0; i < npix; ++i) {
953  double flux = m_map.flux(i);
954  if (flux < 0.0 ||
955  gammalib::is_notanumber(flux) ||
956  gammalib::is_infinite(flux)) {
957  m_map(i) = 0.0;
958  flux = 0.0;
959  }
960  sum += flux;
961  }
962 
963  // Optionally normalize the sky map.
964  if (sum > 0.0) {
965  if (normalize()) {
966  for (int i = 0; i < npix; ++i) {
967  m_map(i) /= sum;
968  }
969  }
970  }
971 
972  // If we have a HealPix map then set radius to 180 deg
973  if (m_map.projection()->code() == "HPX") {
974  m_radius = 180.0;
975  }
976 
977  // ... otherwise compute map centre and radius
978  else {
979 
980  // Get map centre
981  GSkyPixel centre(m_map.nx()/2.0, m_map.ny()/2.0);
982  m_centre = m_map.pix2dir(centre);
983 
984  // Determine map radius
985  for (int i = 0; i < npix; ++i) {
986  double radius = m_map.inx2dir(i).dist_deg(m_centre);
987  if (radius > m_radius) {
988  m_radius = radius;
989  }
990  }
991 
992  } // endelse: computed map centre and radius
993 
994  // Set simulation cone
996 
997  } // endif: there were skymap pixels
998 
999  // Return
1000  return;
1001 }
1002 
1003 
1004 /***********************************************************************//**
1005  * @brief Set boundary sky region
1006  ***************************************************************************/
1008 {
1009  // Set sky region centre to bounding circle centre
1011 
1012  // Set sky region radius to bounding circle radius
1014 
1015  // Return
1016  return;
1017 }
virtual double eval(const GPhoton &photon, const bool &gradients=false) const
Return intensity of sky map.
bool m_normalize
Normalize map (default: true)
Sky map class.
Definition: GSkyMap.hpp:89
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:280
double flux(const int &index, const int &map=0) const
Returns flux in pixel.
Definition: GSkyMap.cpp:1541
bool normalize(void) const
Return normalization flag.
double m_mc_one_minus_cosrad
1-cosine of radius
virtual double mc_norm(const GSkyDir &dir, const double &radius) const
Return normalization of diffuse map for Monte Carlo simulations.
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition: GSkyMap.hpp:366
std::string m_type
Model type.
const double & factor_gradient(void) const
Return parameter gradient factor.
void set_region(void) const
Set boundary sky region.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:414
const std::string & name(void) const
Return parameter name.
void load(const GFilename &filename)
Load skymap into the model class.
double value(void) const
Get model value.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
GSkyRegionCircle m_region
Bounding circle.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
const GSkyMap & map(void) const
Get map.
virtual GModelSpatialDiffuseMap * clone(void) const
Clone spatial map model.
void prepare_map(void)
Prepare sky map after loading.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
virtual void clear(void)
Clear spatial map model.
XML element node class.
Definition: GXmlElement.hpp:47
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
Random number generator class.
Definition: GRan.hpp:44
double m_mc_radius
Radius of MC cone.
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
double m_mc_max
Map maximum for MC.
std::string print(const GChatter &chatter=NORMAL) const
Print sky direction information.
Definition: GSkyDir.cpp:1108
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1713
Interface definition for the spatial model registry class.
HealPix projection class definition.
bool is_free(void) const
Signal if parameter is free.
std::string centre(const std::string &s, const int &n, const char &c= ' ')
Centre string to achieve a length of n characters.
Definition: GTools.cpp:997
virtual void write(GXmlElement &xml) const
Write model into XML element.
GModelSpatialDiffuseMap(void)
Void constructor.
void set_mc_cone(const GSkyDir &centre, const double &radius) const
Set Monte Carlo simulation cone.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
const double & radius(void) const
Return circular region radius (in degrees)
virtual std::string code(void) const =0
Class that handles photons.
Definition: GPhoton.hpp:47
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
double m_mc_norm
Map normalization.
void init_members(void)
Initialise class members.
Sky map pixel class.
Definition: GSkyPixel.hpp:74
const double deg2rad
Definition: GMath.hpp:43
const GModelSpatialDiffuseMap g_spatial_map_seed
GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1350
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void init_members(void)
Initialise class members.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
GSkyDir inx2dir(const int &index) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1321
Filename class.
Definition: GFilename.hpp:62
std::vector< GModelPar * > m_pars
Parameter pointers.
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:321
void fix(void)
Fix a parameter.
bool has_attribute(const std::string &name) const
Check if element has a given attribute.
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC sky direction.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
virtual std::string type(void) const
Return spatial model type.
void clear(void)
Clear parameter.
Spatial model registry class definition.
const GFilename & filename(void) const
Get file name.
double flux(const GSkyRegion *reg, const GEnergy &srcEng=GEnergy(), const GTime &srcTime=GTime()) const
Returns model flux integrated in circular sky region.
GChatter
Definition: GTypemaps.hpp:33
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
virtual GModelSpatialDiffuse & operator=(const GModelSpatialDiffuse &model)
Assignment operator.
#define G_READ
double m_radius
Radius of bounding circle.
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1475
virtual ~GModelSpatialDiffuseMap(void)
Destructor.
void free_members(void)
Delete class members.
GSkyDir m_centre
Centre of bounding circle.
void free_members(void)
Delete class members.
const GSkyDir & centre(void) const
Return circular region centre.
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:142
void free_members(void)
Delete class members.
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
bool m_has_normalize
XML has normalize attribute.
void range(const double &min, const double &max)
Set minimum and maximum parameter boundaries.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
Exception handler interface definition.
#define G_WRITE
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition: GSkyMap.hpp:350
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:834
void init_members(void)
Initialise class members.
#define G_MC
Abstract diffuse spatial model base class.
const double rad2deg
Definition: GMath.hpp:44
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void load(const GFilename &filename)
Load skymap from FITS file.
Definition: GSkyMap.cpp:2281
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:334
void copy_members(const GModelSpatialDiffuseMap &model)
Copy class members.
GFilename m_filename
Name of skymap.
Sky direction class.
Definition: GSkyDir.hpp:62
virtual GModelSpatialDiffuseMap & operator=(const GModelSpatialDiffuseMap &model)
Assignment operator.
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1562
virtual void read(const GXmlElement &xml)
Read model from XML element.
Spatial map model class interface definition.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print map information.
GSkyDir m_mc_centre
Centre of MC cone.
Mathematical function definitions.
void clear(void)
Clear instance.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const
Signals whether model contains sky direction.
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1683
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413