GammaLib  2.0.0
 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-2021 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  // Verify number of model parameters
525 
526  // If "Normalization" parameter exists then read parameter from this
527  // XML element
528  if (gammalib::xml_has_par(xml, "Normalization")) {
529  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Normalization");
530  m_value.read(*par);
531  }
532 
533  // ... otherwise try reading parameter from "Prefactor" parameter
534  #if defined(G_LEGACY_XML_FORMAT)
535  else if (gammalib::xml_has_par(xml, "Prefactor")) {
536  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Prefactor");
537  m_value.read(*par);
538  }
539  #endif
540 
541  // ... otherwise throw an exception
542  else {
543  #if defined(G_LEGACY_XML_FORMAT)
544  std::string msg = "Diffuse map model requires either "
545  "\"Normalization\" or \"Prefactor\" parameter.";
546  #else
547  std::string msg = "Diffuse map model requires \"Normalization\" "
548  "parameter.";
549  #endif
550  throw GException::invalid_value(G_READ, msg);
551  }
552 
553  // Get optional normalization attribute
554  m_normalize = true;
555  m_has_normalize = false;
556  if (xml.has_attribute("normalize")) {
557  m_has_normalize = true;
558  std::string arg = xml.attribute("normalize");
559  if (arg == "0" || gammalib::tolower(arg) == "false") {
560  m_normalize = false;
561  }
562  }
563 
564  // Load sky map.
565  load(gammalib::xml_file_expand(xml, xml.attribute("file")));
566 
567  // Return
568  return;
569 }
570 
571 
572 /***********************************************************************//**
573  * @brief Write model into XML element
574  *
575  * @param[in] xml XML element.
576  *
577  * @exception GException::model_invalid_spatial
578  * Existing XML element is not of type "SpatialMap"
579  * @exception GException::model_invalid_parnum
580  * Invalid number of model parameters found in XML element.
581  * @exception GException::model_invalid_parnames
582  * Invalid model parameter names found in XML element.
583  *
584  * Writes the spatial information for a diffuse map into an XML element. The
585  * format of the XML element is
586  *
587  * <spatialModel type="DiffuseMap" file="myfile.fits" normalize="1">
588  * <parameter name="Prefactor" value="1" min="0.1" max="10" scale="1" free="0"/>
589  * </spatialModel>
590  *
591  * The @p file attribute provides the filename of the diffuse map FITS file.
592  * The filename may be either an absolute filename (starting with '/') or a
593  * relative filename. If no access path is given, the file is expected to
594  * reside in the same location as the XML file.
595  *
596  * The @p normalize attribute specifies whether the sky map should be
597  * normalized to unity flux or not. The attribute will only be written if the
598  * normalization is disabled.
599  ***************************************************************************/
601 {
602  // Verify model type
604 
605  // Get or create parameter
607 
608  // Write parameters
609  m_value.write(*value);
610 
611  // Set sky map file name
613 
614  // Set optional normalization attribute
615  if (m_has_normalize || !m_normalize) {
616  if (m_normalize) {
617  xml.attribute("normalize", "1");
618  }
619  else {
620  xml.attribute("normalize", "0");
621  }
622  }
623 
624  // Return
625  return;
626 }
627 
628 
629 /***********************************************************************//**
630  * @brief Set Monte Carlo simulation cone
631  *
632  * @param[in] cone Monte Carlo simulation cone.
633  *
634  * Sets the simulation cone that defines the directions that will be
635  * simulated using the mc() method and pre-computes the maximum intensity and
636  * the spatially integrated flux of the map within the simulation cone
637  * region.
638  ***************************************************************************/
640 {
641  // Continue only if the simulation cone has changed
642  if (cone != m_mc_cone) {
643 
644  // Save simulation cone definition
645  m_mc_cone = cone;
646 
647  // Pre-compute 1 - cosine of radius
649 
650  // Initialise map maximum and normalisation
651  m_mc_max = 0.0;
652  m_mc_norm = 0.0;
653 
654  // Determine number of skymap pixels
655  int npix = m_map.npix();
656 
657  // Continue only if there are pixels
658  if (npix > 0) {
659 
660  // Compute flux and maximum map intensity within the simulation cone
661  double sum = 0.0;
662  double sum_map = 0.0;
663  for (int i = 0; i < npix; ++i) {
664 
665  // Get map flux, intensity and distance from MC cone centre
666  double flux = m_map.flux(i);
667  double intensity = m_map(i);
668  double distance = m_mc_cone.centre().dist_deg(m_map.pix2dir(i));
669 
670  // Add flux if positive
671  if (flux > 0.0) {
672  if (distance <= m_mc_cone.radius()) {
673  sum += flux; // flux within simulation cone
674  }
675  sum_map += flux; // total flux
676  }
677 
678  // Update maximum intensity
679  if (distance <= m_mc_cone.radius()) {
680  if (intensity > m_mc_max) {
681  m_mc_max = intensity;
682  }
683  }
684 
685  } // endfor: looped over all map pixels
686 
687  // Set the normalization factor for the MC simulations. In case
688  // that the map is normalised, this is the fraction of the flux
689  // that is comprised within the simulation cone. For non-normalised
690  // maps, this is simply the flux comprised within the simulation
691  // cone.
692  if (sum_map > 0.0) {
693  if (m_normalize) {
694  m_mc_norm = sum / sum_map;
695  }
696  else {
697  m_mc_norm = sum;
698  }
699  }
700 
701  // Log maximum intensity and total flux for debugging
702  #if defined(G_DEBUG_MC_CACHE)
703  std::cout << "GModelSpatialDiffuseMap::mc_cone:" << std::endl;
704  std::cout << " Maximum map intensity:";
705  std::cout << m_mc_max << " ph/cm2/s/sr" << std::endl;
706  std::cout << " Spatially integrated flux:" << std::endl;
707  std::cout << sum << " ph/cm2/s" << std::endl;
708  std::cout << " Map normalisation:" << std::endl;
709  std::cout << m_mc_norm << " ph/cm2/s" << std::endl;
710  #endif
711 
712  } // endif: there were map pixels
713 
714  } // endif: simulation cone has changed
715 
716  // Return
717  return;
718 }
719 
720 
721 /***********************************************************************//**
722  * @brief Returns diffuse map flux integrated in sky region
723  *
724  * @param[in] region Sky region.
725  * @param[in] srcEng Energy.
726  * @param[in] srcTime Time.
727  * @return Flux (adimensional or ph/cm2/s).
728  *
729  * Returns diffuse map flux within a sky region. The flux is defined as the
730  * sum of the flux in all diffuse map pixels that are contained within the
731  * sky region multiplied by the diffuse map normalisation factor.
732  ***************************************************************************/
734  const GEnergy& srcEng,
735  const GTime& srcTime) const
736 {
737  // Compute flux in sky region
738  double flux = m_map.flux(region) * m_value.value();
739 
740  // Return flux
741  return flux;
742 }
743 
744 
745 /***********************************************************************//**
746  * @brief Print map information
747  *
748  * @param[in] chatter Chattiness.
749  * @return String with diffuse map model information.
750  ***************************************************************************/
751 std::string GModelSpatialDiffuseMap::print(const GChatter& chatter) const
752 {
753  // Initialise result string
754  std::string result;
755 
756  // Continue only if chatter is not silent
757  if (chatter != SILENT) {
758 
759  // Append header
760  result.append("=== GModelSpatialDiffuseMap ===");
761 
762  // Append parameters
763  result.append("\n"+gammalib::parformat("Sky map file"));
764  result.append(m_filename);
765  result.append("\n"+gammalib::parformat("Map normalization"));
766  result.append(gammalib::str(m_mc_norm)+" ph/cm2/s");
767  if (normalize()) {
768  result.append(" [normalized]");
769  }
770  result.append("\n"+gammalib::parformat("Map centre"));
771  result.append(m_region.centre().print());
772  result.append("\n"+gammalib::parformat("Map radius"));
773  result.append(gammalib::str(m_region.radius())+" deg");
774  result.append("\n"+gammalib::parformat("Number of parameters"));
775  result.append(gammalib::str(size()));
776  for (int i = 0; i < size(); ++i) {
777  result.append("\n"+m_pars[i]->print(chatter));
778  }
779 
780  } // endif: chatter was not silent
781 
782  // Return result
783  return result;
784 }
785 
786 
787 /***********************************************************************//**
788  * @brief Load skymap into the model class
789  *
790  * @param[in] filename Sky map file.
791  *
792  * Loads skymap into the model class. The method calls the protected method
793  * prepare_map() that prepares the map for usage by the class.
794  *
795  * If the @p filename is empty, no map will be loaded.
796  ***************************************************************************/
798 {
799  // Initialise skymap
800  m_map.clear();
801 
802  // Continue only if filename is not empty
803  if (!filename.is_empty()) {
804 
805  // Store filename of skymap (for XML writing). Note that we do not
806  // expand any environment variable at this level, so that if we
807  // write back the XML element we write the filepath with the
808  // environment variables
810 
811  // Load sky map
813 
814  // Prepare sky map
815  prepare_map();
816 
817  } // endif: filename was not empty
818 
819  // Return
820  return;
821 }
822 
823 
824 /*==========================================================================
825  = =
826  = Private methods =
827  = =
828  ==========================================================================*/
829 
830 /***********************************************************************//**
831  * @brief Initialise class members
832  ***************************************************************************/
834 {
835  // Initialise model type
836  m_type = "DiffuseMap";
837 
838  // Initialise Value
839  m_value.clear();
840  m_value.name("Normalization");
841  m_value.value(1.0);
842  m_value.scale(1.0);
843  m_value.range(0.001, 1000.0);
844  m_value.gradient(0.0);
845  m_value.fix();
846  m_value.has_grad(true);
847 
848  // Set parameter pointer(s)
849  m_pars.clear();
850  m_pars.push_back(&m_value);
851 
852  // Initialise other members
853  m_map.clear();
854  m_filename.clear();
855  m_normalize = true;
856  m_has_normalize = false;
857 
858  // Initialise MC cache
859  m_mc_cone.clear();
860  m_mc_one_minus_cosrad = 1.0;
861  m_mc_norm = 0.0;
862  m_mc_max = 0.0;
863 
864  // Return
865  return;
866 }
867 
868 
869 /***********************************************************************//**
870  * @brief Copy class members
871  *
872  * @param[in] model Spatial map cube model.
873  ***************************************************************************/
875 {
876  // Copy members
877  m_type = model.m_type; // Needed to conserve model type
878  m_value = model.m_value;
879  m_map = model.m_map;
880  m_filename = model.m_filename;
881  m_normalize = model.m_normalize;
883 
884  // Copy MC cache
885  m_mc_cone = model.m_mc_cone;
887  m_mc_norm = model.m_mc_norm;
888  m_mc_max = model.m_mc_max;
889 
890  // Set parameter pointer(s)
891  m_pars.clear();
892  m_pars.push_back(&m_value);
893 
894  // Return
895  return;
896 }
897 
898 
899 /***********************************************************************//**
900  * @brief Delete class members
901  ***************************************************************************/
903 {
904  // Return
905  return;
906 }
907 
908 
909 /***********************************************************************//**
910  * @brief Prepare sky map after loading
911  *
912  * Prepares a sky map after loading. Negative, infinite or undefined skymap
913  * pixels are set to zero intensity. The method also determine the centre
914  * and radius of a circle enclosing the map, and the Monte Carlo simulation
915  * cone is set to this circle.
916  *
917  * If normalize() is true, the map is furthermore normalised so that the
918  * total flux in the map amounts to 1 ph/cm2/s. Negative skymap pixels are
919  * set to zero intensity.
920  ***************************************************************************/
922 {
923  // Initialise region
924  m_region.clear();
925 
926  // Determine number of skymap pixels
927  int npix = m_map.npix();
928 
929  // Continue only if there are skymap pixels
930  if (npix > 0) {
931 
932  // Initialise flux sum
933  double sum = 0.0;
934 
935  // Compute flux sum and set negative or invalid pixels to zero
936  // intensity.
937  for (int i = 0; i < npix; ++i) {
938  double flux = m_map.flux(i);
939  if (flux < 0.0 ||
940  gammalib::is_notanumber(flux) ||
941  gammalib::is_infinite(flux)) {
942  m_map(i) = 0.0;
943  flux = 0.0;
944  }
945  sum += flux;
946  }
947 
948  // Optionally normalize the sky map.
949  if (sum > 0.0) {
950  if (normalize()) {
951  for (int i = 0; i < npix; ++i) {
952  m_map(i) /= sum;
953  }
954  }
955  }
956 
957  // Get region circle
959 
960  // Set simulation cone
961  mc_cone(m_region);
962 
963  } // endif: there were skymap pixels
964 
965  // Return
966  return;
967 }
968 
969 
970 /***********************************************************************//**
971  * @brief Set boundary sky region
972  ***************************************************************************/
974 {
975  // Return
976  return;
977 }
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:286
double flux(const int &index, const int &map=0) const
Returns flux in pixel.
Definition: GSkyMap.cpp:1551
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:864
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:1190
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:48
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:944
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
Time class.
Definition: GTime.hpp:55
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:1227
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1946
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:201
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:184
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:1360
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void init_members(void)
Initialise class members.
Abstract interface for the sky region class.
Definition: GSkyRegion.hpp:57
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:424
void fix(void)
Fix a parameter.
bool has_attribute(const std::string &name) const
Check if element has a given attribute.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1637
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.
GChatter
Definition: GTypemaps.hpp:33
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1100
virtual GModelSpatialDiffuse & operator=(const GModelSpatialDiffuse &model)
Assignment operator.
std::string type(void) const
Return model type.
#define G_READ
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
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:1596
virtual ~GModelSpatialDiffuseMap(void)
Destructor.
GSkyRegionCircle region_circle(void) const
Return sky region circle that encloses the sky map.
Definition: GSkyMap.cpp:2401
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.
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.
virtual double flux(const GSkyRegion &region, const GEnergy &srcEng=GEnergy(), const GTime &srcTime=GTime()) const
Returns diffuse map flux integrated in sky region.
#define G_WRITE
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:955
void init_members(void)
Initialise class members.
#define G_MC
Abstract diffuse spatial model base class.
const double rad2deg
Definition: GMath.hpp:44
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void load(const GFilename &filename)
Load skymap from FITS file.
Definition: GSkyMap.cpp:2511
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:345
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:1689
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:1889
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
GSkyRegionCircle m_mc_cone
MC cone.
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819