GammaLib  2.0.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-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 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 = cube().npix();
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_cone.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  ***************************************************************************/
516  const double& margin) const
517 {
518  // Make sure that cube in online
519  fetch_cube();
520 
521  // Initialise containment flag
522  bool contains = false;
523 
524  // Continue only if radius is positive
525  if (m_region.radius() > 0.0) {
526 
527  // Compute distance to centre
528  double distance = m_region.centre().dist_deg(dir);
529 
530  // If distance is smaller than radius plus margin we consider
531  // the position to be contained within the bounding circle
532  if (distance < m_region.radius() + margin) {
533  contains = true;
534  }
535 
536  }
537 
538  // Return containment
539  return (contains);
540 }
541 
542 
543 /***********************************************************************//**
544  * @brief Read model from XML element
545  *
546  * @param[in] xml XML element.
547  *
548  * @exception GException::invalid_value
549  * Model parameters not found in XML element.
550  *
551  * Read the map cube information from an XML element. The XML element should
552  * have the format
553  *
554  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
555  * <parameter name="Normalization" ../>
556  * </spatialModel>
557  ***************************************************************************/
559 {
560  // If "Normalization" parameter exists then read parameter from this
561  // XML element
562  if (gammalib::xml_has_par(xml, "Normalization")) {
563  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Normalization");
564  m_value.read(*par);
565  }
566 
567  // ... otherwise try reading parameter from "Value" parameter
568  #if defined(G_LEGACY_XML_FORMAT)
569  else if (gammalib::xml_has_par(xml, "Value")) {
570  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Value");
571  m_value.read(*par);
572  }
573  #endif
574 
575  // ... otherwise throw an exception
576  else {
577  #if defined(G_LEGACY_XML_FORMAT)
578  std::string msg = "Diffuse map cube model requires either "
579  "\"Normalization\" or \"Value\" parameter.";
580  #else
581  std::string msg = "Diffuse map cube model requires \"Normalization\" "
582  "parameter.";
583  #endif
584  throw GException::invalid_value(G_READ, msg);
585  }
586 
587  // Save filename
588  m_filename = gammalib::xml_file_expand(xml, xml.attribute("file"));
589 
590  // Return
591  return;
592 }
593 
594 
595 /***********************************************************************//**
596  * @brief Write model into XML element
597  *
598  * @param[in] xml XML element.
599  *
600  * @exception GException::model_invalid_spatial
601  * Existing XML element is not of type "MapCubeFunction"
602  * @exception GException::model_invalid_parnum
603  * Invalid number of model parameters found in XML element.
604  * @exception GException::model_invalid_parnames
605  * Invalid model parameter names found in XML element.
606  *
607  * Write the map cube information into an XML element. The XML element will
608  * have either the format
609  *
610  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
611  * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
612  * </spatialModel>
613  *
614  * or alternatively
615  *
616  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
617  * <parameter name="Value" scale="1" value="1" min="0.1" max="10" free="0"/>
618  * </spatialModel>
619  *
620  * The latter format is the default for newly written XML elements.
621  ***************************************************************************/
623 {
624  // Set model type
625  if (xml.attribute("type") == "") {
626  xml.attribute("type", type());
627  }
628 
629  // Verify model type
630  if (xml.attribute("type") != type()) {
632  "Spatial model is not of type \""+type()+"\".");
633  }
634 
635  // If XML element has 0 nodes then append parameter node. The name
636  // of the node is "Normalization" as this is the Fermi-LAT standard.
637  // We thus assure that the XML files will be compatible with
638  // Fermi-LAT.
639  if (xml.elements() == 0) {
640  xml.append(GXmlElement("parameter name=\"Normalization\""));
641  }
642 
643  // Verify that XML element has exactly 1 parameter
644  if (xml.elements() != 1 || xml.elements("parameter") != 1) {
646  "Map cube spatial model requires exactly 1 parameter.");
647  }
648 
649  // Get pointers on model parameter
650  GXmlElement* par = xml.element("parameter", 0);
651 
652  // Set or update parameter
653  if (par->attribute("name") == "Normalization" ||
654  par->attribute("name") == "Value") {
655  m_value.write(*par);
656  }
657  else {
659  "Map cube spatial model requires either \"Value\" or"
660  " \"Normalization\" parameter.");
661  }
662 
663  // Set filename
665 
666  // Return
667  return;
668 }
669 
670 
671 /***********************************************************************//**
672  * @brief Set map cube
673  *
674  * @param[in] cube Sky map.
675  *
676  * Set the map cube of the spatial map cube model.
677  ***************************************************************************/
679 {
680  // Clear filename and signal map as unloaded
681  m_filename.clear();
682  m_loaded = false;
683 
684  // Assign map
685  m_cube = cube;
686 
687  // Update Monte Carlo cache
688  update_mc_cache();
689 
690  // Return
691  return;
692 }
693 
694 
695 /***********************************************************************//**
696  * @brief Set energies for map cube
697  *
698  * @param[in] energies Sky map energies.
699  *
700  * @exception GException::invalid_argument
701  * Specified sky map energies incompatible with map cube.
702  *
703  * Sets the energies for the map cube.
704  ***************************************************************************/
706 {
707  // Initialise energies
708  m_logE.clear();
709 
710  // Fetch cube
711  fetch_cube();
712 
713  // Extract number of energies in vector
714  int num = energies.size();
715 
716  // Check if energy binning is consistent with number of maps in the cube
717  if (num != m_cube.nmaps() ) {
718  std::string msg = "Number of specified energies ("+gammalib::str(num)+")"
719  " does not match the number of maps ("
720  ""+gammalib::str(m_cube.nmaps())+" in the map cube.\n"
721  "The energies argument shall provide a vector of length"
722  " "+gammalib::str(m_cube.nmaps())+".";
724  }
725 
726  // Set log10(energy) nodes, where energy is in units of MeV
727  for (int i = 0; i < num; ++i) {
728  m_logE.append(energies[i].log10MeV());
729  }
730 
731  // Set energy boundaries
733 
734  // Update MC cache
735  update_mc_cache();
736 
737  // Return
738  return;
739 }
740 
741 
742 /***********************************************************************//**
743  * @brief Returns map cube energies
744  *
745  * @return Map cube energies.
746  *
747  * Returns the energies for the map cube in a vector.
748  ***************************************************************************/
750 {
751  // Initialise energies container
753 
754  // Fetch cube
755  fetch_cube();
756 
757  // Get number of map energies
758  int num = m_logE.size();
759 
760  // Continue only if there are maps in the cube
761  if (num > 0) {
762 
763  // Reserve space for all energies
764  energies.reserve(num);
765 
766  // Set log10(energy) nodes, where energy is in units of MeV
767  for (int i = 0; i < num; ++i) {
768  GEnergy energy;
769  energy.log10MeV(m_logE[i]);
770  energies.append(energy);
771  }
772 
773  } // endif: there were maps in the cube
774 
775  // Return energies
776  return energies;
777 }
778 
779 
780 /***********************************************************************//**
781  * @brief Set Monte Carlo simulation cone
782  *
783  * @param[in] cone Monte Carlo simulation cone.
784  *
785  * Sets the simulation cone that defines the directions that will be
786  * simulated using the mc() method and pre-computes the maximum intensity
787  * and the spatially integrated flux of each map within the simulation cone
788  * region.
789  ***************************************************************************/
791 {
792  // Continue only if the simulation cone has changed
793  if (cone != m_mc_cone) {
794 
795  // Save simulation cone definition
796  m_mc_cone = cone;
797 
798  // Pre-compute 1 - cosine of radius
800 
801  // Initialise Monte Carlo cache
802  m_mc_max.clear();
804 
805  // Fetch cube
806  fetch_cube();
807 
808  // Determine number of cube pixels and maps
809  int npix = cube().npix();
810  int nmaps = cube().nmaps();
811 
812  // Continue only if there are pixels and maps
813  if (npix > 0 && nmaps > 0) {
814 
815  // Initial vectors for total flux and maximum intensity
816  std::vector<double> total_flux(nmaps, 0.0);
817  std::vector<double> max_intensity(nmaps, 0.0);
818 
819  // Loop over all map pixels
820  for (int k = 0; k < npix; ++k) {
821 
822  // Continue only if pixel is within MC cone
823  double distance = m_mc_cone.centre().dist_deg(m_cube.pix2dir(k));
824  if (distance <= m_mc_cone.radius()) {
825 
826  // Loop over all maps
827  for (int i = 0; i < nmaps; ++i) {
828 
829  // Update total flux
830  double flux = m_cube.flux(k,i);
831  if (flux > 0.0) {
832  total_flux[i] += flux;
833  }
834 
835  // Update maximum intensity
836  double value = m_cube(k,i);
837  if (value > max_intensity[i]) {
838  max_intensity[i] = value;
839  }
840 
841  } // endfor: looped over all maps
842 
843  } // endif: pixel within MC cone
844 
845  } // endfor: looped over all pixels
846 
847  // Reserve space in cache
848  m_mc_max.reserve(nmaps);
849  m_mc_spectrum.reserve(nmaps);
850 
851  // Set maximum intensity and spectrum vectors
852  for (int i = 0; i < nmaps; ++i) {
853 
854  // Store maximum map intensity
855  m_mc_max.push_back(max_intensity[i]);
856 
857  // Store flux as spectral node
858  if (m_logE.size() == nmaps) {
859 
860  // Set map energy
861  GEnergy energy;
862  energy.log10MeV(m_logE[i]);
863 
864  // Only append node if the integrated flux is positive
865  // (as we do a log-log interpolation)
866  if (total_flux[i] > 0.0) {
867  m_mc_spectrum.append(energy, total_flux[i]);
868  }
869 
870  }
871 
872  } // endfor: looped over all maps
873 
874  // Log maximum intensity and total flux for debugging
875  #if defined(G_DEBUG_MC_CACHE)
876  std::cout << "GModelSpatialDiffuseCube::mc_cone:" << std::endl;
877  std::cout << " Maximum map intensity:" << std::endl;
878  for (int i = 0; i < m_mc_max.size(); ++i) {
879  GEnergy energy;
880  energy.log10MeV(m_logE[i]);
881  std::cout << " " << i;
882  std::cout << " " << energy;
883  std::cout << " " << m_mc_max[i] << " ph/cm2/s/sr/MeV";
884  std::cout << std::endl;
885  }
886  std::cout << " Spatially integrated flux:" << std::endl;
887  for (int i = 0; i < m_mc_spectrum.nodes(); ++i) {
888  std::cout << " " << i;
889  std::cout << " " << m_mc_spectrum.energy(i);
890  std::cout << " " << m_mc_spectrum.intensity(i);
891  std::cout << " ph/cm2/s/MeV" << std::endl;
892  }
893  #endif
894 
895  } // endif: there were cube pixels and maps
896 
897  } // endif: simulation cone has changed
898 
899  // Return
900  return;
901 }
902 
903 
904 /***********************************************************************//**
905  * @brief Load cube into the model class
906  *
907  * @param[in] filename cube file.
908  *
909  * Loads cube into the model class.
910  ***************************************************************************/
912 {
913  // Load cube
914  load_cube(filename);
915 
916  // Update Monte Carlo cache
917  update_mc_cache();
918 
919  // Return
920  return;
921 }
922 
923 
924 /***********************************************************************//**
925  * @brief Save cube into FITS file
926  *
927  * @param[in] filename Cube FITS file name.
928  * @param[in] clobber Overwrite existing FITS file (default: false).
929  *
930  * Saves spatial cube model into a FITS file.
931  ***************************************************************************/
933  const bool& clobber) const
934 {
935  // Create empty FITS file
936  GFits fits;
937 
938  // Write event cube
939  write(fits);
940 
941  // Save FITS file
942  fits.saveto(filename, clobber);
943 
944  // Return
945  return;
946 }
947 
948 
949 /***********************************************************************//**
950  * @brief Read diffuse cube from FITS file
951  *
952  * @param[in] fits FITS file.
953  *
954  * @exception GException::invalid_value
955  * ENERGIES extension incompatible with map cube.
956  *
957  * Reads the diffuse cube sky map from the first extension in the @p fits
958  * file and the energy extension from the "ENERGIES" extension.
959  ***************************************************************************/
961 {
962  // Initialise skymap
963  m_cube.clear();
964  m_logE.clear();
965 
966  // Read cube from the first valid WCS image
967  for (int extno = 0; extno < fits.size(); ++extno) {
968  const GFitsHDU& hdu = *fits.at(extno);
969  if ((hdu.exttype() == GFitsHDU::HT_IMAGE) &&
970  (hdu.has_card("NAXIS") && (hdu.integer("NAXIS") >= 2))) {
971  m_cube.read(hdu);
972  break;
973  }
974  }
975 
976  // Read energies
978  energies.read(*(fits.table(gammalib::extname_energies)));
979 
980  // Extract number of energy bins
981  int num = energies.size();
982 
983  // Check if energy binning is consistent with primary image hdu
984  if (num != m_cube.nmaps() ) {
985  std::string msg = "Number of energies in \"ENERGIES\" extension ("+
986  gammalib::str(num)+") does not match the number of "
987  "maps ("+gammalib::str(m_cube.nmaps())+" in the "
988  "map cube. The \"ENERGIES\" extension table shall "
989  "provide one enegy value for each map in the cube.";
991  }
992 
993  // Set log10(energy) nodes, where energy is in units of MeV
994  for (int i = 0; i < num; ++i) {
995  m_logE.append(energies[i].log10MeV());
996  }
997 
998  // Signal that cube has been loaded
999  m_loaded = true;
1000 
1001  // Set energy boundaries
1003 
1004  // Return
1005  return;
1006 }
1007 
1008 
1009 /***********************************************************************//**
1010  * @brief Write diffuse cube into FITS file
1011  *
1012  * @param[in] fits FITS file.
1013  *
1014  * Writes the diffuse cube sky map and the energy intervals into @p fits
1015  * file.
1016  ***************************************************************************/
1018 {
1019  // Make sure that cube is fetched
1020  fetch_cube();
1021 
1022  // Write cube
1023  m_cube.write(fits);
1024 
1025  // Create energy intervals
1027  for (int i = 0; i < m_logE.size(); ++i) {
1028  double energy = std::pow(10.0, m_logE[i]);
1029  energies.append(GEnergy(energy, "MeV"));
1030  }
1031 
1032  // Write energy intervals
1033  energies.write(fits);
1034 
1035  // Return
1036  return;
1037 }
1038 
1039 
1040 /***********************************************************************//**
1041  * @brief Print map cube information
1042  *
1043  * @param[in] chatter Chattiness.
1044  * @return String containing model information.
1045  ***************************************************************************/
1046 std::string GModelSpatialDiffuseCube::print(const GChatter& chatter) const
1047 {
1048  // Initialise result string
1049  std::string result;
1050 
1051  // Continue only if chatter is not silent
1052  if (chatter != SILENT) {
1053 
1054  // Append header
1055  result.append("=== GModelSpatialDiffuseCube ===");
1056 
1057  // Append parameters
1058  result.append("\n"+gammalib::parformat("Map cube file"));
1059  result.append(m_filename);
1060  if (m_loaded) {
1061  result.append(" [loaded]");
1062  }
1063  result.append("\n"+gammalib::parformat("Number of parameters"));
1064  result.append(gammalib::str(size()));
1065  for (int i = 0; i < size(); ++i) {
1066  result.append("\n"+m_pars[i]->print(chatter));
1067  }
1068 
1069  // Append detailed information only if a map cube exists
1070  if (m_cube.npix() > 0) {
1071 
1072  // NORMAL: Append sky map
1073  if (chatter >= NORMAL) {
1074  result.append("\n"+m_cube.print(chatter));
1075  }
1076 
1077  // EXPLICIT: Append energy nodes
1078  if (chatter >= EXPLICIT && m_logE.size() > 0) {
1079  result.append("\n"+gammalib::parformat("Map energy values"));
1080  if (m_logE.size() > 0) {
1081  for (int i = 0; i < m_logE.size(); ++i) {
1082  result.append("\n"+gammalib::parformat(" Map "+gammalib::str(i+1)));
1083  result.append(gammalib::str(std::pow(10.0, m_logE[i])));
1084  result.append(" MeV (log10E=");
1085  result.append(gammalib::str(m_logE[i]));
1086  result.append(")");
1087  if (m_ebounds.size() == m_logE.size()) {
1088  result.append(" [");
1089  result.append(m_ebounds.emin(i).print());
1090  result.append(", ");
1091  result.append(m_ebounds.emax(i).print());
1092  result.append("]");
1093  }
1094  }
1095  }
1096  else {
1097  result.append("not specified");
1098  }
1099  }
1100 
1101  // VERBOSE: Append MC cache
1102  if (chatter >= VERBOSE) {
1103  result.append("\n"+gammalib::parformat("Map flux"));
1104  if (m_mc_spectrum.nodes() > 0) {
1105  for (int i = 0; i < m_mc_spectrum.nodes(); ++i) {
1106  result.append("\n"+gammalib::parformat(" Map "+gammalib::str(i+1)));
1107  result.append(gammalib::str(m_mc_spectrum.intensity(i)));
1108  }
1109  }
1110  else {
1111  result.append("not specified");
1112  }
1113  }
1114 
1115  } // endif: map cube exists
1116 
1117  } // endif: chatter was not silent
1118 
1119  // Return result
1120  return result;
1121 }
1122 
1123 
1124 /*==========================================================================
1125  = =
1126  = Private methods =
1127  = =
1128  ==========================================================================*/
1129 
1130 /***********************************************************************//**
1131  * @brief Initialise class members
1132  ***************************************************************************/
1134 {
1135  // Initialise model type
1136  m_type = "DiffuseMapCube";
1137 
1138  // Initialise Value
1139  m_value.clear();
1140  m_value.name("Normalization");
1141  m_value.value(1.0);
1142  m_value.scale(1.0);
1143  m_value.range(0.001, 1000.0);
1144  m_value.gradient(0.0);
1145  m_value.fix();
1146  m_value.has_grad(true);
1147 
1148  // Set parameter pointer(s)
1149  m_pars.clear();
1150  m_pars.push_back(&m_value);
1151 
1152  // Initialise other members
1153  m_filename.clear();
1154  m_cube.clear();
1155  m_logE.clear();
1156  m_ebounds.clear();
1157  m_loaded = false;
1158 
1159  // Initialise MC cache
1160  m_mc_cone.clear();
1161  m_mc_one_minus_cosrad = 1.0;
1162  m_mc_max.clear();
1163  m_mc_spectrum.clear();
1164 
1165  // Return
1166  return;
1167 }
1168 
1169 
1170 /***********************************************************************//**
1171  * @brief Copy class members
1172  *
1173  * @param[in] model Spatial map cube model.
1174  ***************************************************************************/
1176 {
1177  // Copy members
1178  m_type = model.m_type; // Needed to conserve model type
1179  m_value = model.m_value;
1180  m_filename = model.m_filename;
1181  m_cube = model.m_cube;
1182  m_logE = model.m_logE;
1183  m_ebounds = model.m_ebounds;
1184  m_loaded = model.m_loaded;
1185 
1186  // Copy MC cache
1187  m_mc_cone = model.m_mc_cone;
1189  m_mc_max = model.m_mc_max;
1190  m_mc_spectrum = model.m_mc_spectrum;
1191 
1192  // Set parameter pointer(s)
1193  m_pars.clear();
1194  m_pars.push_back(&m_value);
1195 
1196  // Return
1197  return;
1198 }
1199 
1200 
1201 /***********************************************************************//**
1202  * @brief Delete class members
1203  ***************************************************************************/
1205 {
1206  // Return
1207  return;
1208 }
1209 
1210 
1211 /***********************************************************************//**
1212  * @brief Fetch cube
1213  *
1214  * Load diffuse cube if it is not yet loaded. The loading is thread save.
1215  ***************************************************************************/
1217 {
1218  // Load cube if it is not yet loaded
1219  if (!m_loaded && !m_filename.is_empty()) {
1220 
1221  // Put in a OMP critical zone
1222  #pragma omp critical(GModelSpatialDiffuseCube_fetch_cube)
1223  {
1224  const_cast<GModelSpatialDiffuseCube*>(this)->load_cube(m_filename);
1225  } // end of pragma
1226 
1227  } // endif: file not
1228 
1229  // Return
1230  return;
1231 }
1232 
1233 
1234 /***********************************************************************//**
1235  * @brief Load map cube
1236  *
1237  * @param[in] filename Map cube file.
1238  *
1239  * @exception GException::invalid_value
1240  * Number of maps in cube mismatches number of energy bins.
1241  *
1242  * Load diffuse map cube.
1243  ***************************************************************************/
1245 {
1246  // Initialise sky map cube, region and log(energy) nodes
1247  m_cube.clear();
1248  m_region.clear();
1249  m_logE.clear();
1250 
1251  // Store filename of cube (for XML writing). Note that we do
1252  // not expand any environment variable at this level, so that
1253  // if we write back the XML element we write the filepath with
1254  // the environment variables
1255  m_filename = filename;
1256 
1257  // Load cube
1258  m_cube.load(filename);
1259 
1260  // Load energies
1261  GEnergies energies(filename);
1262 
1263  // Extract number of energy bins
1264  int num = energies.size();
1265 
1266  // Check if energy binning is consistent with primary image hdu
1267  if (num != m_cube.nmaps() ) {
1268  std::string msg = "Number of energies in \"ENERGIES\" extension ("+
1269  gammalib::str(num)+") does not match the number of "
1270  "maps ("+gammalib::str(m_cube.nmaps())+" in the "
1271  "map cube. The \"ENERGIES\" extension table shall "
1272  "provide one enegy value for each map in the cube.";
1274  }
1275 
1276  // Set log10(energy) nodes, where energy is in units of MeV
1277  for (int i = 0; i < num; ++i) {
1278  m_logE.append(energies[i].log10MeV());
1279  }
1280 
1281  // Signal that cube has been loaded
1282  m_loaded = true;
1283 
1284  // Set energy boundaries
1286 
1287  // Get region circle
1289 
1290  // Set simulation cone
1291  mc_cone(m_region);
1292 
1293  // Return
1294  return;
1295 }
1296 
1297 
1298 /***********************************************************************//**
1299  * @brief Set energy boundaries
1300  *
1301  * Computes energy boundaries from the energy nodes. The boundaries will be
1302  * computed as the centre in log10(energy) between the nodes. If only a
1303  * single map exists, no boundaries will be computed.
1304  ***************************************************************************/
1306 {
1307  // Initialise energy boundaries
1308  m_ebounds.clear();
1309 
1310  // Determine number of energy bins
1311  int num = m_logE.size();
1312 
1313  // Continue only if there are at least two energy nodes
1314  if (num > 1) {
1315 
1316  // Loop over all nodes
1317  for (int i = 0; i < num; ++i) {
1318 
1319  // Calculate minimum energy
1320  double e_min = (i == 0) ? m_logE[i] - 0.5 * (m_logE[i+1] - m_logE[i])
1321  : m_logE[i] - 0.5 * (m_logE[i] - m_logE[i-1]);
1322 
1323  // Calculate maximum energy
1324  double e_max = (i < num-1) ? m_logE[i] + 0.5 * (m_logE[i+1] - m_logE[i])
1325  : m_logE[i] + 0.5 * (m_logE[i] - m_logE[i-1]);
1326 
1327  // Set energy boundaries
1328  GEnergy emin;
1329  GEnergy emax;
1330  emin.log10MeV(e_min);
1331  emax.log10MeV(e_max);
1332 
1333  // Append energy bin to energy boundary arra
1334  m_ebounds.append(emin, emax);
1335 
1336  } // endfor: looped over energy nodes
1337 
1338  } // endif: there were at least two energy nodes
1339 
1340  // Return
1341  return;
1342 }
1343 
1344 
1345 /***********************************************************************//**
1346  * @brief Update Monte Carlo cache
1347  *
1348  * Initialise the cache for Monte Carlo sampling of the map cube. See the
1349  * mc_cone() for more information.
1350  ***************************************************************************/
1352 {
1353  // Initialise sky region to all sky (dummy)
1354  GSkyRegionCircle cone(0.0, 0.0, 360.0);
1355 
1356  // Set cone
1357  mc_cone(cone);
1358 
1359  // Return
1360  return;
1361 }
1362 
1363 
1364 /***********************************************************************//**
1365  * @brief Compute cube intensity by log-log interpolation
1366  *
1367  * @param[in] photon Incident photon.
1368  * @return Cube intensity value.
1369  ***************************************************************************/
1371 {
1372  // Fetch cube
1373  fetch_cube();
1374 
1375  // Initialise intensity value
1376  double intensity = 0.0;
1377 
1378  // Continue only if there is energy information for the map cube
1379  if (m_logE.size() > 0) {
1380 
1381  // Set energy for interpolation in log-energy
1382  m_logE.set_value(photon.energy().log10MeV());
1383 
1384  // Compute map cube intensity for the left and right map
1385  double left_intensity = m_cube(photon.dir(), m_logE.inx_left());
1386  double right_intensity = m_cube(photon.dir(), m_logE.inx_right());
1387 
1388  // Perform log-log interpolation
1389  if (left_intensity > 0.0 && right_intensity > 0.0) {
1390  double log_intensity = m_logE.wgt_left() * std::log(left_intensity) +
1391  m_logE.wgt_right() * std::log(right_intensity);
1392  intensity = std::exp(log_intensity);
1393  }
1394  else if (left_intensity > 0.0) {
1395  intensity = left_intensity;
1396  }
1397  else if (right_intensity > 0.0) {
1398  intensity = right_intensity;
1399  }
1400 
1401  } // endif: energy information was available
1402 
1403  // Return intensity
1404  return intensity;
1405 }
1406 
1407 
1408 /***********************************************************************//**
1409  * @brief Set boundary sky region
1410  ***************************************************************************/
1412 {
1413  // Make sure that cube in online
1414  fetch_cube();
1415 
1416  // Return
1417  return;
1418 }
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:192
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 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:611
virtual 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 factor gradient.
void fetch_cube(void) const
Fetch cube.
#define G_LOAD_CUBE
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:655
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:351
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:162
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:300
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1107
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:198
XML element node class.
Definition: GXmlElement.hpp:47
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:264
#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
GSkyRegionCircle m_mc_cone
MC cone.
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
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.
Interface for the circular sky region class.
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
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.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2550
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
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
GSkyRegionCircle m_region
Bounding circle.
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:2673
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
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void init_members(void)
Initialise class members.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:279
void save(const GFilename &filename, const bool &clobber=false) const
Save cube into FITS file.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:379
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:2583
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:389
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:186
void fix(void)
Fix a parameter.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:170
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1191
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
const GSkyRegionCircle & mc_cone(void) const
Get Monte Carlo simulation cone.
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:235
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.
std::string type(void) const
Return model type.
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:267
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:879
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:1486
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:249
const std::string extname_energies
Definition: GEnergies.hpp:44
double intensity(const int &index) const
Return node intensity.
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.
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 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.
std::string m_type
Spatial model type.
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:1339
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
GEnergy & append(const GEnergy &energy)
Append energy to container.
Definition: GEnergies.cpp:305
Exception handler interface definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1149
void init_members(void)
Initialise class members.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:198
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 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:2397
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:335
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:1579
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:1703
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415
FITS table abstract base class interface definition.