GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpatialDiffuseCube.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpatialDiffuseCube.cpp - Spatial map cube model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2018 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 GModelSpatialDiffuseCube.cpp
23  * @brief Spatial map cube 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"
35 #include "GEnergies.hpp"
36 #include "GFits.hpp"
37 #include "GFitsTable.hpp"
38 #include "GFitsTableCol.hpp"
39 #include "GXmlElement.hpp"
42 
43 /* __ Constants __________________________________________________________ */
44 
45 /* __ Globals ____________________________________________________________ */
47 const GModelSpatialRegistry g_spatial_cube_registry(&g_spatial_cube_seed);
48 #if defined(G_LEGACY_XML_FORMAT)
49 const GModelSpatialDiffuseCube g_spatial_cube_legacy_seed(true, "MapCubeFunction");
50 const GModelSpatialRegistry g_spatial_cube_legacy_registry(&g_spatial_cube_legacy_seed);
51 #endif
52 
53 /* __ Method name definitions ____________________________________________ */
54 #define G_MC "GModelSpatialDiffuseCube::mc(GEnergy&, GTime&, GRan&)"
55 #define G_READ "GModelSpatialDiffuseCube::read(GXmlElement&)"
56 #define G_WRITE "GModelSpatialDiffuseCube::write(GXmlElement&)"
57 #define G_ENERGIES "GModelSpatialDiffuseCube::energies(GEnergies&)"
58 #define G_LOAD_CUBE "GModelSpatialDiffuseCube::load_cube(GFilename&)"
59 #define G_READ_FITS "GModelSpatialDiffuseCube::read(GFits&)"
60 
61 /* __ Macros _____________________________________________________________ */
62 
63 /* __ Coding definitions _________________________________________________ */
64 
65 /* __ Debug definitions __________________________________________________ */
66 //#define G_DEBUG_MC //!< Debug MC method
67 //#define G_DEBUG_MC_CACHE //!< Debug MC cache
68 
69 
70 /*==========================================================================
71  = =
72  = Constructors/destructors =
73  = =
74  ==========================================================================*/
75 
76 /***********************************************************************//**
77  * @brief Void constructor
78  *
79  * Constructs empty map cube model.
80  ***************************************************************************/
83 {
84  // Initialise members
85  init_members();
86 
87  // Return
88  return;
89 }
90 
91 
92 /***********************************************************************//**
93  * @brief Model type constructor
94  *
95  * @param[in] dummy Dummy flag.
96  * @param[in] type Model type.
97  *
98  * Constructs empty map cube model by specifying a model @p type.
99  ***************************************************************************/
101  const std::string& type) :
103 {
104  // Initialise members
105  init_members();
106 
107  // Set model type
108  m_type = type;
109 
110  // Return
111  return;
112 }
113 
114 
115 /***********************************************************************//**
116  * @brief XML constructor
117  *
118  * @param[in] xml XML element.
119  *
120  * Constructs map cube model by extracting information from an XML element.
121  * See the read() method for more information about the expected structure
122  * of the XML element.
123  ***************************************************************************/
126 {
127  // Initialise members
128  init_members();
129 
130  // Read information from XML element
131  read(xml);
132 
133  // Return
134  return;
135 }
136 
137 
138 /***********************************************************************//**
139  * @brief Filename constructor
140  *
141  * @param[in] filename File name.
142  * @param[in] value Normalization factor (defaults to 1).
143  *
144  * Constructs map cube model by loading a map cube from @p filename and by
145  * assigning the normalization @p value.
146  ***************************************************************************/
148  const double& value) :
150 {
151  // Initialise members
152  init_members();
153 
154  // Set parameter
155  m_value.value(value);
156 
157  // Perform autoscaling of parameter
158  autoscale();
159 
160  // Store filename
162 
163  // Return
164  return;
165 }
166 
167 
168 /***********************************************************************//**
169  * @brief Sky map constructor
170  *
171  * @param[in] cube Sky map cube.
172  * @param[in] energies Sky map energies.
173  * @param[in] value Normalization factor (defaults to 1).
174  *
175  * Constructs map cube model by extracting a @p cube from a sky map. The
176  * constructor also assigns the energy values for all maps and sets the
177  * scaling @p value. The filename will remain blank.
178  ***************************************************************************/
180  const GEnergies& energies,
181  const double& value) :
183 {
184  // Initialise members
185  init_members();
186 
187  // Set parameter
188  m_value.value(value);
189 
190  // Perform autoscaling of parameter
191  autoscale();
192 
193  // Set map cube
194  this->cube(cube);
195 
196  // Set energies
197  this->energies(energies);
198 
199  // Return
200  return;
201 }
202 
203 
204 /***********************************************************************//**
205  * @brief Copy constructor
206  *
207  * @param[in] model Map cube model.
208  ***************************************************************************/
210  GModelSpatialDiffuse(model)
211 {
212  // Initialise members
213  init_members();
214 
215  // Copy members
216  copy_members(model);
217 
218  // Return
219  return;
220 }
221 
222 
223 /***********************************************************************//**
224  * @brief Destructor
225  ***************************************************************************/
227 {
228  // Free members
229  free_members();
230 
231  // Return
232  return;
233 }
234 
235 
236 /*==========================================================================
237  = =
238  = Operators =
239  = =
240  ==========================================================================*/
241 
242 /***********************************************************************//**
243  * @brief Assignment operator
244  *
245  * @param[in] model Map cube model.
246  * @return Map cube model.
247  ***************************************************************************/
249 {
250  // Execute only if object is not identical
251  if (this != &model) {
252 
253  // Copy base class members
254  this->GModelSpatialDiffuse::operator=(model);
255 
256  // Free members
257  free_members();
258 
259  // Initialise members
260  init_members();
261 
262  // Copy members
263  copy_members(model);
264 
265  } // endif: object was not identical
266 
267  // Return
268  return *this;
269 }
270 
271 
272 /*==========================================================================
273  = =
274  = Public methods =
275  = =
276  ==========================================================================*/
277 
278 /***********************************************************************//**
279  * @brief Clear map cube model
280  ***************************************************************************/
282 {
283  // Free class members (base and derived classes, derived class first)
284  free_members();
287 
288  // Initialise members
291  init_members();
292 
293  // Return
294  return;
295 }
296 
297 
298 /***********************************************************************//**
299  * @brief Clone map cube model
300  *
301  * @return Pointer to deep copy of map cube model.
302  ***************************************************************************/
304 {
305  // Clone map cube model
306  return new GModelSpatialDiffuseCube(*this);
307 }
308 
309 
310 /***********************************************************************//**
311  * @brief Evaluate function
312  *
313  * @param[in] photon Incident photon.
314  * @param[in] gradients Compute gradients?
315  * @return Sky map intensity (\f$\mbox{ph cm}^{-2}\mbox{sr}^{-1}\mbox{s}^{-1}\f$)
316  *
317  * Computes the spatial diffuse model as function of photon parameters.
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 log-log interpolated cube intensity
326  double intensity = cube_intensity(photon);
327 
328  // Set the intensity times the scaling factor as model value
329  double value = intensity * m_value.value();
330 
331  // Optionally compute partial derivatives
332  if (gradients) {
333 
334  // Compute partial derivatives of the parameter value. In case that
335  // the value is negative set the gradient to zero.
336  double g_value = (m_value.is_free()) ? intensity * m_value.scale() : 0.0;
337  if (value < 0.0) {
338  g_value = 0.0;
339  }
340 
341  // Set gradient
342  m_value.factor_gradient(g_value);
343 
344  } // endif: computed partial derivatives
345 
346  // Make sure that value is not negative
347  if (value < 0.0) {
348  value = 0.0;
349  }
350 
351  // Return value
352  return value;
353 }
354 
355 
356 /***********************************************************************//**
357  * @brief Returns MC sky direction
358  *
359  * @param[in] energy Photon energy.
360  * @param[in] time Photon arrival time.
361  * @param[in,out] ran Random number generator.
362  * @return Sky direction.
363  *
364  * @exception GException::invalid_value
365  * No map cube defined.
366  * No energy boundaries defined.
367  * @exception GException::invalid_return_value
368  * Simulation cone not defined, does not overlap with map cube or
369  * map cube is empty for the specified energy.
370  *
371  * Returns a random sky direction according to the intensity distribution of
372  * the model sky map and the specified energy. The method uses a rejection
373  * method to determine the sky direction. If no sky direction could be
374  * determined, the method throws an GException::invalid_return_value
375  * exception.
376  ***************************************************************************/
378  const GTime& time,
379  GRan& ran) const
380 {
381  // Fetch cube
382  fetch_cube();
383 
384  // Determine number of skymap pixels
385  int npix = pixels();
386 
387  // Throw an exception if there are no sky map pixels
388  if (npix <= 0) {
389  std::string msg = "No map cube defined. Please specify a valid map cube.";
390  throw GException::invalid_value(G_MC, msg);
391  }
392 
393  // Throw an exception if no energy boundaries are defined
394  if (m_ebounds.size() < 1) {
395  std::string msg = "The energy boundaries of the maps in the cube have "
396  "not been defined. Maybe the map cube file is missing "
397  "the \"ENERGIES\" extension which defines the energy "
398  "of each map in the cube. Please provide the energy "
399  "information.";
400  throw GException::invalid_value(G_MC, msg);
401  }
402 
403  // Set energy for interpolation
404  m_logE.set_value(energy.log10MeV());
405 
406  // Compute maximum map value within simulation cone by log-log interpolation
407  // of the maximum values that are stored for each map
408  double max = 0.0;
409  double max_left = m_mc_max[m_logE.inx_left()];
410  double max_right = m_mc_max[m_logE.inx_right()];
411  if (max_left > 0.0 && max_right > 0.0) {
412  double max_log = m_logE.wgt_left() * std::log(max_left) +
413  m_logE.wgt_right() * std::log(max_right);
414  max = std::exp(max_log);
415  }
416  else if (max_left > 0.0) {
417  max = max_left;
418  }
419  else if (max_right > 0.0) {
420  max = max_right;
421  }
422 
423  // Throw an exception if the maximum MC intensity is not positive. This
424  // can happen if the simulation cone has not been defined or if there is
425  // no overlap with the sky map or if the sky map is empty
426  if (max <= 0.0) {
427  std::string msg;
428  if (m_mc_max.size() == 0) {
429  msg = "Simulation cone has not been defined. Please specify a "
430  "valid simulation cone before calling the method.";
431  }
432  else {
433  msg = "The map cube is empty at "+energy.print()+" within the "
434  "simulation cone. Please specify a valid map cube or "
435  "restrict the simulated energies to values where the map "
436  "cube is non-zero within the simulation cone.";
437  }
439  }
440 
441  // Allocate sky direction
442  GSkyDir dir;
443 
444  // Debug option: initialise counter
445  #if defined(G_DEBUG_MC)
446  int num_iterations = 0;
447  #endif
448 
449  // Get sky direction
450  while (true) {
451 
452  // Debug option: increment counter
453  #if defined(G_DEBUG_MC)
454  num_iterations++;
455  #endif
456 
457  // Simulate random sky direction within Monte Carlo simulation cone
458  double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
460  double phi = 360.0 * ran.uniform();
461  dir = m_mc_centre;
462  dir.rotate_deg(phi, theta);
463 
464  // Get map value at simulated sky direction
465  double value = 0.0;
466  double value_left = m_cube(dir, m_logE.inx_left());
467  double value_right = m_cube(dir, m_logE.inx_right());
468  if (value_left > 0.0 && value_right > 0.0) {
469  double value_log = m_logE.wgt_left() * std::log(value_left) +
470  m_logE.wgt_right() * std::log(value_right);
471  value = std::exp(value_log);
472  }
473  else if (value_left > 0.0) {
474  value = value_left;
475  }
476  else if (value_right > 0.0) {
477  value = value_right;
478  }
479 
480  // If map value is non-positive then simulate a new sky direction
481  if (value <= 0.0) {
482  continue;
483  }
484 
485  // Get uniform random number
486  double uniform = ran.uniform() * max;
487 
488  // Exit loop if the random number is not larger than the map cube value
489  if (uniform <= value) {
490  break;
491  }
492 
493  } // endwhile: loop until sky direction was accepted
494 
495  // Debug option: log counter
496  #if defined(G_DEBUG_MC)
497  std::cout << num_iterations << " ";
498  #endif
499 
500  // Return sky direction
501  return dir;
502 }
503 
504 
505 /***********************************************************************//**
506  * @brief Signals whether model contains sky direction
507  *
508  * @param[in] dir Sky direction.
509  * @param[in] margin Margin to be added to sky direction (degrees)
510  * @return True if the model contains the sky direction.
511  *
512  * Signals whether a sky direction falls within the bounding circle of
513  * the diffuse cube.
514  *
515  * @todo To be implemented.
516  ***************************************************************************/
518  const double& margin) const
519 {
520  return (true);
521 }
522 
523 
524 /***********************************************************************//**
525  * @brief Read model from XML element
526  *
527  * @param[in] xml XML element.
528  *
529  * @exception GException::invalid_value
530  * Model parameters not found in XML element.
531  *
532  * Read the map cube information from an XML element. The XML element should
533  * have the format
534  *
535  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
536  * <parameter name="Normalization" ../>
537  * </spatialModel>
538  ***************************************************************************/
540 {
541  // If "Normalization" parameter exists then read parameter from this
542  // XML element
543  if (gammalib::xml_has_par(xml, "Normalization")) {
544  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Normalization");
545  m_value.read(*par);
546  }
547 
548  // ... otherwise try reading parameter from "Value" parameter
549  #if defined(G_LEGACY_XML_FORMAT)
550  else if (gammalib::xml_has_par(xml, "Value")) {
551  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Value");
552  m_value.read(*par);
553  }
554  #endif
555 
556  // ... otherwise throw an exception
557  else {
558  #if defined(G_LEGACY_XML_FORMAT)
559  std::string msg = "Diffuse map cube model requires either "
560  "\"Normalization\" or \"Value\" parameter.";
561  #else
562  std::string msg = "Diffuse map cube model requires \"Normalization\" "
563  "parameter.";
564  #endif
565  throw GException::invalid_value(G_READ, msg);
566  }
567 
568  // Save filename
569  m_filename = gammalib::xml_file_expand(xml, xml.attribute("file"));
570 
571  // Return
572  return;
573 }
574 
575 
576 /***********************************************************************//**
577  * @brief Write model into XML element
578  *
579  * @param[in] xml XML element.
580  *
581  * @exception GException::model_invalid_spatial
582  * Existing XML element is not of type "MapCubeFunction"
583  * @exception GException::model_invalid_parnum
584  * Invalid number of model parameters found in XML element.
585  * @exception GException::model_invalid_parnames
586  * Invalid model parameter names found in XML element.
587  *
588  * Write the map cube information into an XML element. The XML element will
589  * have either the format
590  *
591  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
592  * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
593  * </spatialModel>
594  *
595  * or alternatively
596  *
597  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
598  * <parameter name="Value" scale="1" value="1" min="0.1" max="10" free="0"/>
599  * </spatialModel>
600  *
601  * The latter format is the default for newly written XML elements.
602  ***************************************************************************/
604 {
605  // Set model type
606  if (xml.attribute("type") == "") {
607  xml.attribute("type", type());
608  }
609 
610  // Verify model type
611  if (xml.attribute("type") != type()) {
613  "Spatial model is not of type \""+type()+"\".");
614  }
615 
616  // If XML element has 0 nodes then append parameter node. The name
617  // of the node is "Normalization" as this is the Fermi-LAT standard.
618  // We thus assure that the XML files will be compatible with
619  // Fermi-LAT.
620  if (xml.elements() == 0) {
621  xml.append(GXmlElement("parameter name=\"Normalization\""));
622  }
623 
624  // Verify that XML element has exactly 1 parameter
625  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
627  "Map cube spatial model requires exactly 1 parameter.");
628  }
629 
630  // Get pointers on model parameter
631  GXmlElement* par = xml.element("parameter", 0);
632 
633  // Set or update parameter
634  if (par->attribute("name") == "Normalization" ||
635  par->attribute("name") == "Value") {
636  m_value.write(*par);
637  }
638  else {
640  "Map cube spatial model requires either \"Value\" or"
641  " \"Normalization\" parameter.");
642  }
643 
644  // Set filename
646 
647  // Return
648  return;
649 }
650 
651 
652 /***********************************************************************//**
653  * @brief Set map cube
654  *
655  * @param[in] cube Sky map.
656  *
657  * Set the map cube of the spatial map cube model.
658  ***************************************************************************/
660 {
661  // Clear filename and signal map as unloaded
662  m_filename.clear();
663  m_loaded = false;
664 
665  // Assign map
666  m_cube = cube;
667 
668  // Update Monte Carlo cache
669  update_mc_cache();
670 
671  // Return
672  return;
673 }
674 
675 
676 /***********************************************************************//**
677  * @brief Set energies for map cube
678  *
679  * @param[in] energies Sky map energies.
680  *
681  * @exception GException::invalid_argument
682  * Specified sky map energies incompatible with map cube.
683  *
684  * Sets the energies for the map cube.
685  ***************************************************************************/
687 {
688  // Initialise energies
689  m_logE.clear();
690 
691  // Fetch cube
692  fetch_cube();
693 
694  // Extract number of energies in vector
695  int num = energies.size();
696 
697  // Check if energy binning is consistent with number of maps in the cube
698  if (num != m_cube.nmaps() ) {
699  std::string msg = "Number of specified energies ("+gammalib::str(num)+")"
700  " does not match the number of maps ("
701  ""+gammalib::str(m_cube.nmaps())+" in the map cube.\n"
702  "The energies argument shall provide a vector of length"
703  " "+gammalib::str(m_cube.nmaps())+".";
705  }
706 
707  // Set log10(energy) nodes, where energy is in units of MeV
708  for (int i = 0; i < num; ++i) {
709  m_logE.append(energies[i].log10MeV());
710  }
711 
712  // Set energy boundaries
714 
715  // Update MC cache
716  update_mc_cache();
717 
718  // Return
719  return;
720 }
721 
722 
723 /***********************************************************************//**
724  * @brief Returns map cube energies
725  *
726  * @return Map cube energies.
727  *
728  * Returns the energies for the map cube in a vector.
729  ***************************************************************************/
731 {
732  // Initialise energies container
734 
735  // Fetch cube
736  fetch_cube();
737 
738  // Get number of map energies
739  int num = m_logE.size();
740 
741  // Continue only if there are maps in the cube
742  if (num > 0) {
743 
744  // Reserve space for all energies
745  energies.reserve(num);
746 
747  // Set log10(energy) nodes, where energy is in units of MeV
748  for (int i = 0; i < num; ++i) {
749  GEnergy energy;
750  energy.log10MeV(m_logE[i]);
751  energies.append(energy);
752  }
753 
754  } // endif: there were maps in the cube
755 
756  // Return energies
757  return energies;
758 }
759 
760 
761 /***********************************************************************//**
762  * @brief Set Monte Carlo simulation cone
763  *
764  * @param[in] centre Simulation cone centre.
765  * @param[in] radius Simulation cone radius (degrees).
766  *
767  * Sets the simulation cone centre and radius that defines the directions
768  * that will be simulated using the mc() method and pre-computes the maximum
769  * intensity and the spatially integrated flux of each map within the
770  * simulation cone region.
771  ***************************************************************************/
773  const double& radius) const
774 {
775  // Continue only if the simulation cone has changed
776  if ((centre != m_mc_centre) || (radius != m_mc_radius)) {
777 
778  // Save simulation cone definition
780  m_mc_radius = radius;
781 
782  // Pre-compute 1 - cosine of radius
784 
785  // Initialise Monte Carlo cache
786  m_mc_max.clear();
788 
789  // Fetch cube
790  fetch_cube();
791 
792  // Determine number of cube pixels and maps
793  int npix = pixels();
794  int nmaps = maps();
795 
796  // Continue only if there are pixels and maps
797  if (npix > 0 && nmaps > 0) {
798 
799  // Initial vectors for total flux and maximum intensity
800  std::vector<double> total_flux(nmaps, 0.0);
801  std::vector<double> max_intensity(nmaps, 0.0);
802 
803  // Loop over all map pixels
804  for (int k = 0; k < npix; ++k) {
805 
806  // Continue only if pixel is within MC cone
807  double distance = centre.dist_deg(m_cube.pix2dir(k));
808  if (distance <= radius) {
809 
810  // Loop over all maps
811  for (int i = 0; i < nmaps; ++i) {
812 
813  // Update total flux
814  double flux = m_cube.flux(k,i);
815  if (flux > 0.0) {
816  total_flux[i] += flux;
817  }
818 
819  // Update maximum intensity
820  double value = m_cube(k,i);
821  if (value > max_intensity[i]) {
822  max_intensity[i] = value;
823  }
824 
825  } // endfor: looped over all maps
826 
827  } // endif: pixel within MC cone
828 
829  } // endfor: looped over all pixels
830 
831  // Reserve space in cache
832  m_mc_max.reserve(nmaps);
833  m_mc_spectrum.reserve(nmaps);
834 
835  // Set maximum intensity and spectrum vectors
836  for (int i = 0; i < nmaps; ++i) {
837 
838  // Store maximum map intensity
839  m_mc_max.push_back(max_intensity[i]);
840 
841  // Store flux as spectral node
842  if (m_logE.size() == nmaps) {
843 
844  // Set map energy
845  GEnergy energy;
846  energy.log10MeV(m_logE[i]);
847 
848  // Only append node if the integrated flux is positive
849  // (as we do a log-log interpolation)
850  if (total_flux[i] > 0.0) {
851  m_mc_spectrum.append(energy, total_flux[i]);
852  }
853 
854  }
855 
856  } // endfor: looped over all maps
857 
858  // Log maximum intensity and total flux for debugging
859  #if defined(G_DEBUG_MC_CACHE)
860  std::cout << "GModelSpatialDiffuseCube::set_mc_cone:" << std::endl;
861  std::cout << " Maximum map intensity:" << std::endl;
862  for (int i = 0; i < m_mc_max.size(); ++i) {
863  GEnergy energy;
864  energy.log10MeV(m_logE[i]);
865  std::cout << " " << i;
866  std::cout << " " << energy;
867  std::cout << " " << m_mc_max[i] << " ph/cm2/s/sr/MeV";
868  std::cout << std::endl;
869  }
870  std::cout << " Spatially integrated flux:" << std::endl;
871  for (int i = 0; i < m_mc_spectrum.nodes(); ++i) {
872  std::cout << " " << i;
873  std::cout << " " << m_mc_spectrum.energy(i);
874  std::cout << " " << m_mc_spectrum.intensity(i);
875  std::cout << " ph/cm2/s/MeV" << std::endl;
876  }
877  #endif
878 
879  } // endif: there were cube pixels and maps
880 
881  } // endif: simulation cone has changed
882 
883  // Return
884  return;
885 }
886 
887 
888 /***********************************************************************//**
889  * @brief Load cube into the model class
890  *
891  * @param[in] filename cube file.
892  *
893  * Loads cube into the model class.
894  ***************************************************************************/
896 {
897  // Load cube
898  load_cube(filename);
899 
900  // Update Monte Carlo cache
901  update_mc_cache();
902 
903  // Return
904  return;
905 }
906 
907 
908 /***********************************************************************//**
909  * @brief Save cube into FITS file
910  *
911  * @param[in] filename Cube FITS file name.
912  * @param[in] clobber Overwrite existing FITS file (default: false).
913  *
914  * Saves spatial cube model into a FITS file.
915  ***************************************************************************/
917  const bool& clobber) const
918 {
919  // Create empty FITS file
920  GFits fits;
921 
922  // Write event cube
923  write(fits);
924 
925  // Save FITS file
926  fits.saveto(filename, clobber);
927 
928  // Return
929  return;
930 }
931 
932 
933 /***********************************************************************//**
934  * @brief Read diffuse cube from FITS file
935  *
936  * @param[in] fits FITS file.
937  *
938  * @exception GException::invalid_value
939  * ENERGIES extension incompatible with map cube.
940  *
941  * Reads the diffuse cube sky map from the first extension in the @p fits
942  * file and the energy extension from the "ENERGIES" extension.
943  ***************************************************************************/
945 {
946  // Initialise skymap
947  m_cube.clear();
948  m_logE.clear();
949 
950  // Read cube from the first valid WCS image
951  for (int extno = 0; extno < fits.size(); ++extno) {
952  const GFitsHDU& hdu = *fits.at(extno);
953  if ((hdu.exttype() == GFitsHDU::HT_IMAGE) &&
954  (hdu.has_card("NAXIS") && (hdu.integer("NAXIS") >= 2))) {
955  m_cube.read(hdu);
956  break;
957  }
958  }
959 
960  // Read energies
962  energies.read(*(fits.table(gammalib::extname_energies)));
963 
964  // Extract number of energy bins
965  int num = energies.size();
966 
967  // Check if energy binning is consistent with primary image hdu
968  if (num != m_cube.nmaps() ) {
969  std::string msg = "Number of energies in \"ENERGIES\" extension ("+
970  gammalib::str(num)+") does not match the number of "
971  "maps ("+gammalib::str(m_cube.nmaps())+" in the "
972  "map cube. The \"ENERGIES\" extension table shall "
973  "provide one enegy value for each map in the cube.";
975  }
976 
977  // Set log10(energy) nodes, where energy is in units of MeV
978  for (int i = 0; i < num; ++i) {
979  m_logE.append(energies[i].log10MeV());
980  }
981 
982  // Signal that cube has been loaded
983  m_loaded = true;
984 
985  // Set energy boundaries
987 
988  // Return
989  return;
990 }
991 
992 
993 /***********************************************************************//**
994  * @brief Write diffuse cube into FITS file
995  *
996  * @param[in] fits FITS file.
997  *
998  * Writes the diffuse cube sky map and the energy intervals into @p fits
999  * file.
1000  ***************************************************************************/
1002 {
1003  // Make sure that cube is fetched
1004  fetch_cube();
1005 
1006  // Write cube
1007  m_cube.write(fits);
1008 
1009  // Create energy intervals
1011  for (int i = 0; i < m_logE.size(); ++i) {
1012  double energy = std::pow(10.0, m_logE[i]);
1013  energies.append(GEnergy(energy, "MeV"));
1014  }
1015 
1016  // Write energy intervals
1017  energies.write(fits);
1018 
1019  // Return
1020  return;
1021 }
1022 
1023 
1024 /***********************************************************************//**
1025  * @brief Print map cube information
1026  *
1027  * @param[in] chatter Chattiness.
1028  * @return String containing model information.
1029  ***************************************************************************/
1030 std::string GModelSpatialDiffuseCube::print(const GChatter& chatter) const
1031 {
1032  // Initialise result string
1033  std::string result;
1034 
1035  // Continue only if chatter is not silent
1036  if (chatter != SILENT) {
1037 
1038  // Append header
1039  result.append("=== GModelSpatialDiffuseCube ===");
1040 
1041  // Append parameters
1042  result.append("\n"+gammalib::parformat("Map cube file"));
1043  result.append(m_filename);
1044  if (m_loaded) {
1045  result.append(" [loaded]");
1046  }
1047  result.append("\n"+gammalib::parformat("Number of parameters"));
1048  result.append(gammalib::str(size()));
1049  for (int i = 0; i < size(); ++i) {
1050  result.append("\n"+m_pars[i]->print(chatter));
1051  }
1052 
1053  // Append detailed information only if a map cube exists
1054  if (m_cube.npix() > 0) {
1055 
1056  // NORMAL: Append sky map
1057  if (chatter >= NORMAL) {
1058  result.append("\n"+m_cube.print(chatter));
1059  }
1060 
1061  // EXPLICIT: Append energy nodes
1062  if (chatter >= EXPLICIT && m_logE.size() > 0) {
1063  result.append("\n"+gammalib::parformat("Map energy values"));
1064  if (m_logE.size() > 0) {
1065  for (int i = 0; i < m_logE.size(); ++i) {
1066  result.append("\n"+gammalib::parformat(" Map "+gammalib::str(i+1)));
1067  result.append(gammalib::str(std::pow(10.0, m_logE[i])));
1068  result.append(" MeV (log10E=");
1069  result.append(gammalib::str(m_logE[i]));
1070  result.append(")");
1071  if (m_ebounds.size() == m_logE.size()) {
1072  result.append(" [");
1073  result.append(m_ebounds.emin(i).print());
1074  result.append(", ");
1075  result.append(m_ebounds.emax(i).print());
1076  result.append("]");
1077  }
1078  }
1079  }
1080  else {
1081  result.append("not specified");
1082  }
1083  }
1084 
1085  // VERBOSE: Append MC cache
1086  if (chatter >= VERBOSE) {
1087  result.append("\n"+gammalib::parformat("Map flux"));
1088  if (m_mc_spectrum.nodes() > 0) {
1089  for (int i = 0; i < m_mc_spectrum.nodes(); ++i) {
1090  result.append("\n"+gammalib::parformat(" Map "+gammalib::str(i+1)));
1091  result.append(gammalib::str(m_mc_spectrum.intensity(i)));
1092  }
1093  }
1094  else {
1095  result.append("not specified");
1096  }
1097  }
1098 
1099  } // endif: map cube exists
1100 
1101  } // endif: chatter was not silent
1102 
1103  // Return result
1104  return result;
1105 }
1106 
1107 
1108 /*==========================================================================
1109  = =
1110  = Private methods =
1111  = =
1112  ==========================================================================*/
1113 
1114 /***********************************************************************//**
1115  * @brief Initialise class members
1116  ***************************************************************************/
1118 {
1119  // Initialise model type
1120  m_type = "DiffuseMapCube";
1121 
1122  // Initialise Value
1123  m_value.clear();
1124  m_value.name("Normalization");
1125  m_value.value(1.0);
1126  m_value.scale(1.0);
1127  m_value.range(0.001, 1000.0);
1128  m_value.gradient(0.0);
1129  m_value.fix();
1130  m_value.has_grad(true);
1131 
1132  // Set parameter pointer(s)
1133  m_pars.clear();
1134  m_pars.push_back(&m_value);
1135 
1136  // Initialise other members
1137  m_filename.clear();
1138  m_cube.clear();
1139  m_logE.clear();
1140  m_ebounds.clear();
1141  m_loaded = false;
1142  m_region.clear();
1143 
1144  // Initialise MC cache
1145  m_mc_centre.clear();
1146  m_mc_radius = -1.0; //!< Signal that initialisation is needed
1147  m_mc_one_minus_cosrad = 1.0;
1148  m_mc_max.clear();
1149  m_mc_spectrum.clear();
1150 
1151  // Return
1152  return;
1153 }
1154 
1155 
1156 /***********************************************************************//**
1157  * @brief Copy class members
1158  *
1159  * @param[in] model Spatial map cube model.
1160  ***************************************************************************/
1162 {
1163  // Copy members
1164  m_type = model.m_type;
1165  m_value = model.m_value;
1166  m_filename = model.m_filename;
1167  m_cube = model.m_cube;
1168  m_logE = model.m_logE;
1169  m_ebounds = model.m_ebounds;
1170  m_loaded = model.m_loaded;
1171  m_region = model.m_region;
1172 
1173  // Copy MC cache
1174  m_mc_centre = model.m_mc_centre;
1175  m_mc_radius = model.m_mc_radius;
1177  m_mc_max = model.m_mc_max;
1178  m_mc_spectrum = model.m_mc_spectrum;
1179 
1180  // Set parameter pointer(s)
1181  m_pars.clear();
1182  m_pars.push_back(&m_value);
1183 
1184  // Return
1185  return;
1186 }
1187 
1188 
1189 /***********************************************************************//**
1190  * @brief Delete class members
1191  ***************************************************************************/
1193 {
1194  // Return
1195  return;
1196 }
1197 
1198 
1199 /***********************************************************************//**
1200  * @brief Fetch cube
1201  *
1202  * Load diffuse cube if it is not yet loaded. The loading is thread save.
1203  ***************************************************************************/
1205 {
1206  // Load cube if it is not yet loaded
1207  if (!m_loaded && !m_filename.is_empty()) {
1208 
1209  // Put in a OMP critical zone
1210  #pragma omp critical(GModelSpatialDiffuseCube_fetch_cube)
1211  {
1212  const_cast<GModelSpatialDiffuseCube*>(this)->load_cube(m_filename);
1213  } // end of pragma
1214 
1215  } // endif: file not
1216 
1217  // Return
1218  return;
1219 }
1220 
1221 
1222 /***********************************************************************//**
1223  * @brief Load map cube
1224  *
1225  * @param[in] filename Map cube file.
1226  *
1227  * @exception GException::invalid_value
1228  * Number of maps in cube mismatches number of energy bins.
1229  *
1230  * Load diffuse map cube.
1231  ***************************************************************************/
1233 {
1234  // Initialise skymap
1235  m_cube.clear();
1236  m_logE.clear();
1237 
1238  // Store filename of cube (for XML writing). Note that we do
1239  // not expand any environment variable at this level, so that
1240  // if we write back the XML element we write the filepath with
1241  // the environment variables
1242  m_filename = filename;
1243 
1244  // Load cube
1245  m_cube.load(filename);
1246 
1247  // Load energies
1248  GEnergies energies(filename);
1249 
1250  // Extract number of energy bins
1251  int num = energies.size();
1252 
1253  // Check if energy binning is consistent with primary image hdu
1254  if (num != m_cube.nmaps() ) {
1255  std::string msg = "Number of energies in \"ENERGIES\" extension ("+
1256  gammalib::str(num)+") does not match the number of "
1257  "maps ("+gammalib::str(m_cube.nmaps())+" in the "
1258  "map cube. The \"ENERGIES\" extension table shall "
1259  "provide one enegy value for each map in the cube.";
1261  }
1262 
1263  // Set log10(energy) nodes, where energy is in units of MeV
1264  for (int i = 0; i < num; ++i) {
1265  m_logE.append(energies[i].log10MeV());
1266  }
1267 
1268  // Signal that cube has been loaded
1269  m_loaded = true;
1270 
1271  // Set energy boundaries
1273 
1274  // Return
1275  return;
1276 }
1277 
1278 
1279 /***********************************************************************//**
1280  * @brief Set energy boundaries
1281  *
1282  * Computes energy boundaries from the energy nodes. The boundaries will be
1283  * computed as the centre in log10(energy) between the nodes. If only a
1284  * single map exists, no boundaries will be computed.
1285  ***************************************************************************/
1287 {
1288  // Initialise energy boundaries
1289  m_ebounds.clear();
1290 
1291  // Determine number of energy bins
1292  int num = m_logE.size();
1293 
1294  // Continue only if there are at least two energy nodes
1295  if (num > 1) {
1296 
1297  // Loop over all nodes
1298  for (int i = 0; i < num; ++i) {
1299 
1300  // Calculate minimum energy
1301  double e_min = (i == 0) ? m_logE[i] - 0.5 * (m_logE[i+1] - m_logE[i])
1302  : m_logE[i] - 0.5 * (m_logE[i] - m_logE[i-1]);
1303 
1304  // Calculate maximum energy
1305  double e_max = (i < num-1) ? m_logE[i] + 0.5 * (m_logE[i+1] - m_logE[i])
1306  : m_logE[i] + 0.5 * (m_logE[i] - m_logE[i-1]);
1307 
1308  // Set energy boundaries
1309  GEnergy emin;
1310  GEnergy emax;
1311  emin.log10MeV(e_min);
1312  emax.log10MeV(e_max);
1313 
1314  // Append energy bin to energy boundary arra
1315  m_ebounds.append(emin, emax);
1316 
1317  } // endfor: looped over energy nodes
1318 
1319  } // endif: there were at least two energy nodes
1320 
1321  // Return
1322  return;
1323 }
1324 
1325 
1326 /***********************************************************************//**
1327  * @brief Update Monte Carlo cache
1328  *
1329  * Initialise the cache for Monte Carlo sampling of the map cube. The Monte
1330  * Carlo cache consists of a linear array that maps a value between 0 and 1
1331  * into the skymap pixel for all maps in the cube.
1332  ***************************************************************************/
1334 {
1335  // Set centre and radius to all sky
1336  GSkyDir centre;
1337  double radius = 360.0;
1338 
1339  // Compute cache
1340  set_mc_cone(centre, radius);
1341 
1342  // Return
1343  return;
1344 }
1345 
1346 
1347 /***********************************************************************//**
1348  * @brief Compute cube intensity by log-log interpolation
1349  *
1350  * @param[in] photon Incident photon.
1351  * @return Cube intensity value.
1352  ***************************************************************************/
1354 {
1355  // Fetch cube
1356  fetch_cube();
1357 
1358  // Initialise intensity value
1359  double intensity = 0.0;
1360 
1361  // Continue only if there is energy information for the map cube
1362  if (m_logE.size() > 0) {
1363 
1364  // Set energy for interpolation in log-energy
1365  m_logE.set_value(photon.energy().log10MeV());
1366 
1367  // Compute map cube intensity for the left and right map
1368  double left_intensity = m_cube(photon.dir(), m_logE.inx_left());
1369  double right_intensity = m_cube(photon.dir(), m_logE.inx_right());
1370 
1371  // Perform log-log interpolation
1372  if (left_intensity > 0.0 && right_intensity > 0.0) {
1373  double log_intensity = m_logE.wgt_left() * std::log(left_intensity) +
1374  m_logE.wgt_right() * std::log(right_intensity);
1375  intensity = std::exp(log_intensity);
1376  }
1377  else if (left_intensity > 0.0) {
1378  intensity = left_intensity;
1379  }
1380  else if (right_intensity > 0.0) {
1381  intensity = right_intensity;
1382  }
1383 
1384  } // endif: energy information was available
1385 
1386  // Return intensity
1387  return intensity;
1388 }
1389 
1390 
1391 /***********************************************************************//**
1392  * @brief Set boundary sky region
1393  *
1394  * @todo Implement determination of the cube boundary circle
1395  ***************************************************************************/
1397 {
1398  // Set sky region centre to (0,0)
1399  m_region.centre(0.0, 0.0);
1400 
1401  // Set sky region radius to 180 degrees (all points included)
1402  m_region.radius(180.0);
1403 
1404  // Return
1405  return;
1406 }
GFitsHDU * at(const int &extno)
Get pointer to HDU.
Definition: GFits.cpp:233
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
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 m_loaded
Signals if map cube has been loaded.
#define G_WRITE
void read(const GFitsTable &table)
Read energies from FITS table.
Definition: GEnergies.cpp:681
void set_region(void) const
Set boundary sky region.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
void free_members(void)
Delete class members.
const double & factor_gradient(void) const
Return parameter gradient factor.
void fetch_cube(void) const
Fetch cube.
#define G_LOAD_CUBE
void set_mc_cone(const GSkyDir &centre, const double &radius) const
Set Monte Carlo simulation cone.
const std::string & name(void) const
Return parameter name.
void write(GFits &file, const std::string &extname=gammalib::extname_energies) const
Write energies into FITS object.
Definition: GEnergies.cpp:725
XML element node class interface definition.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
const GFilename & filename(void) const
Get file name.
void init_members(void)
Initialise class members.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
GFilename m_filename
Name of map cube.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:302
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
void reserve(const int &num)
Reserves space for energies in container.
Definition: GEnergies.hpp:190
int maps(void) const
Return number of maps in cube.
XML element node class.
Definition: GXmlElement.hpp:47
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
#define G_ENERGIES
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
virtual HDUType exttype(void) const =0
virtual std::string print(const GChatter &chatter=NORMAL) const
Print map cube information.
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
int pixels(void) const
Return number of pixels in cube.
Gammalib tools definition.
virtual ~GModelSpatialDiffuseCube(void)
Destructor.
FITS file class.
Definition: GFits.hpp:63
void reserve(const int &num)
Reserve space for nodes.
FITS file class interface definition.
virtual void read(const GXmlElement &xml)
Read model from XML element.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
double m_mc_one_minus_cosrad
1-cosine of radius
virtual std::string type(void) const
Return spatial model type.
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.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2434
HealPix projection class definition.
bool is_free(void) const
Signal if parameter is free.
FITS table column abstract base class definition.
GNodeArray m_logE
Log10(energy) values of the maps.
Energy container class.
Definition: GEnergies.hpp:60
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
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
virtual GModelSpatialDiffuseCube & operator=(const GModelSpatialDiffuseCube &model)
Assignment operator.
const double & radius(void) const
Return circular region radius (in degrees)
void load(const GFilename &filename)
Load cube into the model class.
void copy_members(const GModelSpatialDiffuseCube &model)
Copy class members.
Class that handles photons.
Definition: GPhoton.hpp:47
const double & scale(void) const
Return parameter scale.
std::vector< double > m_mc_max
Maximum values for MC.
void append(const GEnergy &energy, const double &intensity)
Append node.
double value(void) const
Get model value.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1267
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2557
const double deg2rad
Definition: GMath.hpp:43
virtual void write(GXmlElement &xml) const
Write model into XML element.
GSkyDir pix2dir(const GSkyPixel &pixel) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1350
GSkyDir m_mc_centre
Centre of MC cone.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void init_members(void)
Initialise class members.
GSkyRegionCircle m_region
Bounding circle.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
void save(const GFilename &filename, const bool &clobber=false) const
Save cube into FITS file.
double m_mc_radius
Radius of MC cone (degrees)
const GXmlAttribute * attribute(const int &index) const
Return attribute.
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:378
void load_cube(const GFilename &filename)
Load map cube.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
double cube_intensity(const GPhoton &photon) const
Compute cube intensity by log-log interpolation.
void update_mc_cache(void)
Update Monte Carlo cache.
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2467
Filename class.
Definition: GFilename.hpp:62
std::vector< GModelPar * > m_pars
Parameter pointers.
Energy container class definition.
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:321
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
void fix(void)
Fix a parameter.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:162
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
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.
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
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
GEnergy energy(const int &index) const
Return node energy.
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
#define G_READ
virtual GModelSpatialDiffuse & operator=(const GModelSpatialDiffuse &model)
Assignment operator.
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:269
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:872
const GModelSpatialDiffuseCube g_spatial_cube_seed
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
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
const std::string extname_energies
Definition: GEnergies.hpp:44
double intensity(const int &index) const
Return node intensity.
void free_members(void)
Delete class members.
GModelSpatialDiffuseCube(void)
Void constructor.
virtual GSkyDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC sky direction.
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
virtual void clear(void)
Clear spectral nodes model.
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const
Signals whether model contains sky direction.
void set_energy_boundaries(void)
Set energy boundaries.
GEnergies energies(void)
Returns map cube energies.
Spatial map cube model class interface definition.
int nodes(void) const
Return number of nodes.
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.
int size(void) const
Return number of HDUs in FITS file.
Definition: GFits.hpp:237
GModelSpectralNodes m_mc_spectrum
Map cube spectrum.
virtual double eval(const GPhoton &photon, const bool &gradients=false) const
Evaluate function.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
GEnergy & append(const GEnergy &energy)
Append energy to container.
Definition: GEnergies.cpp:303
Exception handler interface definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
void init_members(void)
Initialise class members.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
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 append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
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
Sky direction class.
Definition: GSkyDir.hpp:62
#define G_MC
const GSkyMap & cube(void) const
Get map cube.
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
void autoscale(void)
Autoscale parameters.
virtual void clear(void)
Clear map cube model.
GEbounds m_ebounds
Energy bounds of the maps.
virtual GModelSpatialDiffuseCube * clone(void) const
Clone map cube model.
Mathematical function definitions.
void clear(void)
Clear instance.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
#define G_READ_FITS
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
FITS table abstract base class interface definition.