GammaLib  2.0.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-2020 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_cone.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  mc_cone(GSkyRegionCircle(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_region.radius() > 0.0) {
479 
480  // Compute distance to centre
481  double distance = m_region.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_region.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] cone Monte Carlo simulation cone.
660  *
661  * Sets the simulation cone that defines the directions that will be
662  * simulated using the mc() method and pre-computes the maximum intensity and
663  * the spatially integrated flux of the map within the simulation cone
664  * region.
665  ***************************************************************************/
667 {
668  // Continue only if the simulation cone has changed
669  if (cone != m_mc_cone) {
670 
671  // Save simulation cone definition
672  m_mc_cone = cone;
673 
674  // Pre-compute 1 - cosine of radius
676 
677  // Initialise map maximum and normalisation
678  m_mc_max = 0.0;
679  m_mc_norm = 0.0;
680 
681  // Determine number of skymap pixels
682  int npix = m_map.npix();
683 
684  // Continue only if there are pixels
685  if (npix > 0) {
686 
687  // Compute flux and maximum map intensity within the simulation cone
688  double sum = 0.0;
689  double sum_map = 0.0;
690  for (int i = 0; i < npix; ++i) {
691 
692  // Get map flux, intensity and distance from MC cone centre
693  double flux = m_map.flux(i);
694  double intensity = m_map(i);
695  double distance = m_mc_cone.centre().dist_deg(m_map.pix2dir(i));
696 
697  // Add flux if positive
698  if (flux > 0.0) {
699  if (distance <= m_mc_cone.radius()) {
700  sum += flux; // flux within simulation cone
701  }
702  sum_map += flux; // total flux
703  }
704 
705  // Update maximum intensity
706  if (distance <= m_mc_cone.radius()) {
707  if (intensity > m_mc_max) {
708  m_mc_max = intensity;
709  }
710  }
711 
712  } // endfor: looped over all map pixels
713 
714  // Set the normalization factor for the MC simulations. In case
715  // that the map is normalised, this is the fraction of the flux
716  // that is comprised within the simulation cone. For non-normalised
717  // maps, this is simply the flux comprised within the simulation
718  // cone.
719  if (sum_map > 0.0) {
720  if (m_normalize) {
721  m_mc_norm = sum / sum_map;
722  }
723  else {
724  m_mc_norm = sum;
725  }
726  }
727 
728  // Log maximum intensity and total flux for debugging
729  #if defined(G_DEBUG_MC_CACHE)
730  std::cout << "GModelSpatialDiffuseMap::mc_cone:" << std::endl;
731  std::cout << " Maximum map intensity:";
732  std::cout << m_mc_max << " ph/cm2/s/sr" << std::endl;
733  std::cout << " Spatially integrated flux:" << std::endl;
734  std::cout << sum << " ph/cm2/s" << std::endl;
735  std::cout << " Map normalisation:" << std::endl;
736  std::cout << m_mc_norm << " ph/cm2/s" << std::endl;
737  #endif
738 
739  } // endif: there were map pixels
740 
741  } // endif: simulation cone has changed
742 
743  // Return
744  return;
745 }
746 
747 
748 /***********************************************************************//**
749  * @brief Print map information
750  *
751  * @param[in] chatter Chattiness.
752  * @return String with diffuse map model information.
753  ***************************************************************************/
754 std::string GModelSpatialDiffuseMap::print(const GChatter& chatter) const
755 {
756  // Initialise result string
757  std::string result;
758 
759  // Continue only if chatter is not silent
760  if (chatter != SILENT) {
761 
762  // Append header
763  result.append("=== GModelSpatialDiffuseMap ===");
764 
765  // Append parameters
766  result.append("\n"+gammalib::parformat("Sky map file"));
767  result.append(m_filename);
768  result.append("\n"+gammalib::parformat("Map normalization"));
769  result.append(gammalib::str(m_mc_norm)+" ph/cm2/s");
770  if (normalize()) {
771  result.append(" [normalized]");
772  }
773  result.append("\n"+gammalib::parformat("Map centre"));
774  result.append(m_region.centre().print());
775  result.append("\n"+gammalib::parformat("Map radius"));
776  result.append(gammalib::str(m_region.radius())+" deg");
777  result.append("\n"+gammalib::parformat("Number of parameters"));
778  result.append(gammalib::str(size()));
779  for (int i = 0; i < size(); ++i) {
780  result.append("\n"+m_pars[i]->print(chatter));
781  }
782 
783  } // endif: chatter was not silent
784 
785  // Return result
786  return result;
787 }
788 
789 
790 /***********************************************************************//**
791  * @brief Load skymap into the model class
792  *
793  * @param[in] filename Sky map file.
794  *
795  * Loads skymap into the model class. The method calls the protected method
796  * prepare_map() that prepares the map for usage by the class.
797  *
798  * If the @p filename is empty, no map will be loaded.
799  ***************************************************************************/
801 {
802  // Initialise skymap
803  m_map.clear();
804 
805  // Continue only if filename is not empty
806  if (!filename.is_empty()) {
807 
808  // Store filename of skymap (for XML writing). Note that we do not
809  // expand any environment variable at this level, so that if we
810  // write back the XML element we write the filepath with the
811  // environment variables
813 
814  // Load sky map
816 
817  // Prepare sky map
818  prepare_map();
819 
820  } // endif: filename was not empty
821 
822  // Return
823  return;
824 }
825 
826 
827 /*==========================================================================
828  = =
829  = Private methods =
830  = =
831  ==========================================================================*/
832 
833 /***********************************************************************//**
834  * @brief Initialise class members
835  ***************************************************************************/
837 {
838  // Initialise model type
839  m_type = "DiffuseMap";
840 
841  // Initialise Value
842  m_value.clear();
843  m_value.name("Normalization");
844  m_value.value(1.0);
845  m_value.scale(1.0);
846  m_value.range(0.001, 1000.0);
847  m_value.gradient(0.0);
848  m_value.fix();
849  m_value.has_grad(true);
850 
851  // Set parameter pointer(s)
852  m_pars.clear();
853  m_pars.push_back(&m_value);
854 
855  // Initialise other members
856  m_map.clear();
857  m_filename.clear();
858  m_normalize = true;
859  m_has_normalize = false;
860 
861  // Initialise MC cache
862  m_mc_cone.clear();
863  m_mc_one_minus_cosrad = 1.0;
864  m_mc_norm = 0.0;
865  m_mc_max = 0.0;
866 
867  // Return
868  return;
869 }
870 
871 
872 /***********************************************************************//**
873  * @brief Copy class members
874  *
875  * @param[in] model Spatial map cube model.
876  ***************************************************************************/
878 {
879  // Copy members
880  m_type = model.m_type; // Needed to conserve model type
881  m_value = model.m_value;
882  m_map = model.m_map;
883  m_filename = model.m_filename;
884  m_normalize = model.m_normalize;
886 
887  // Copy MC cache
888  m_mc_cone = model.m_mc_cone;
890  m_mc_norm = model.m_mc_norm;
891  m_mc_max = model.m_mc_max;
892 
893  // Set parameter pointer(s)
894  m_pars.clear();
895  m_pars.push_back(&m_value);
896 
897  // Return
898  return;
899 }
900 
901 
902 /***********************************************************************//**
903  * @brief Delete class members
904  ***************************************************************************/
906 {
907  // Return
908  return;
909 }
910 
911 
912 /***********************************************************************//**
913  * @brief Prepare sky map after loading
914  *
915  * Prepares a sky map after loading. Negative, infinite or undefined skymap
916  * pixels are set to zero intensity. The method also determine the centre
917  * and radius of a circle enclosing the map, and the Monte Carlo simulation
918  * cone is set to this circle.
919  *
920  * If normalize() is true, the map is furthermore normalised so that the
921  * total flux in the map amounts to 1 ph/cm2/s. Negative skymap pixels are
922  * set to zero intensity.
923  ***************************************************************************/
925 {
926  // Initialise region
927  m_region.clear();
928 
929  // Determine number of skymap pixels
930  int npix = m_map.npix();
931 
932  // Continue only if there are skymap pixels
933  if (npix > 0) {
934 
935  // Initialise flux sum
936  double sum = 0.0;
937 
938  // Compute flux sum and set negative or invalid pixels to zero
939  // intensity.
940  for (int i = 0; i < npix; ++i) {
941  double flux = m_map.flux(i);
942  if (flux < 0.0 ||
943  gammalib::is_notanumber(flux) ||
944  gammalib::is_infinite(flux)) {
945  m_map(i) = 0.0;
946  flux = 0.0;
947  }
948  sum += flux;
949  }
950 
951  // Optionally normalize the sky map.
952  if (sum > 0.0) {
953  if (normalize()) {
954  for (int i = 0; i < npix; ++i) {
955  m_map(i) /= sum;
956  }
957  }
958  }
959 
960  // Get region circle
962 
963  // Set simulation cone
964  mc_cone(m_region);
965 
966  } // endif: there were skymap pixels
967 
968  // Return
969  return;
970 }
971 
972 
973 /***********************************************************************//**
974  * @brief Set boundary sky region
975  ***************************************************************************/
977 {
978  // Return
979  return;
980 }
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:284
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 double & factor_gradient(void) const
Return parameter factor gradient.
virtual void set_region(void) const
Set boundary sky region.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:828
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.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1107
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:908
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
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.
Interface for the circular sky region class.
std::string print(const GChatter &chatter=NORMAL) const
Print sky direction information.
Definition: GSkyDir.cpp:1159
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1733
Interface definition for the spatial model registry class.
HealPix projection class definition.
bool is_free(void) const
Signal if parameter is free.
GSkyRegionCircle m_region
Bounding circle.
virtual void write(GXmlElement &xml) const
Write model into XML element.
GModelSpatialDiffuseMap(void)
Void constructor.
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)
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.
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.
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:389
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
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.
std::string type(void) const
Return model type.
#define G_READ
const GSkyRegionCircle & mc_cone(void) const
Get Monte Carlo simulation cone.
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1486
virtual ~GModelSpatialDiffuseMap(void)
Destructor.
GSkyRegionCircle region_circle(void) const
Return sky region circle that encloses the sky map.
Definition: GSkyMap.cpp:2287
void free_members(void)
Delete class members.
void free_members(void)
Delete class members.
const GSkyDir & centre(void) const
Return circular region centre.
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.
std::string m_type
Spatial model type.
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
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:845
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:1033
void load(const GFilename &filename)
Load skymap from FITS file.
Definition: GSkyMap.cpp:2397
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:335
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:1579
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.
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:1703
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415
GSkyRegionCircle m_mc_cone
MC cone.