GammaLib  2.1.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-2021 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file 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  *
377  * @todo Make sure that spatial model value is taken into account.
378  ***************************************************************************/
380  const GTime& time,
381  GRan& ran) const
382 {
383  // Fetch cube
384  fetch_cube();
385 
386  // Determine number of skymap pixels
387  int npix = cube().npix();
388 
389  // Throw an exception if there are no sky map pixels
390  if (npix <= 0) {
391  std::string msg = "No map cube defined. Please specify a valid map cube.";
392  throw GException::invalid_value(G_MC, msg);
393  }
394 
395  // Throw an exception if no energy boundaries are defined
396  if (m_ebounds.size() < 1) {
397  std::string msg = "The energy boundaries of the maps in the cube have "
398  "not been defined. Maybe the map cube file is missing "
399  "the \"ENERGIES\" extension which defines the energy "
400  "of each map in the cube. Please provide the energy "
401  "information.";
402  throw GException::invalid_value(G_MC, msg);
403  }
404 
405  // Set energy for interpolation
406  m_logE.set_value(energy.log10MeV());
407 
408  // Compute maximum map value within simulation cone by log-log interpolation
409  // of the maximum values that are stored for each map
410  double max = 0.0;
411  double max_left = m_mc_max[m_logE.inx_left()];
412  double max_right = m_mc_max[m_logE.inx_right()];
413  if (max_left > 0.0 && max_right > 0.0) {
414  double max_log = m_logE.wgt_left() * std::log(max_left) +
415  m_logE.wgt_right() * std::log(max_right);
416  max = std::exp(max_log);
417  }
418  else if (max_left > 0.0) {
419  max = max_left;
420  }
421  else if (max_right > 0.0) {
422  max = max_right;
423  }
424 
425  // Throw an exception if the maximum MC intensity is not positive. This
426  // can happen if the simulation cone has not been defined or if there is
427  // no overlap with the sky map or if the sky map is empty
428  if (max <= 0.0) {
429  std::string msg;
430  if (m_mc_max.size() == 0) {
431  msg = "Simulation cone has not been defined. Please specify a "
432  "valid simulation cone before calling the method.";
433  }
434  else {
435  msg = "The map cube is empty at "+energy.print()+" within the "
436  "simulation cone. Please specify a valid map cube or "
437  "restrict the simulated energies to values where the map "
438  "cube is non-zero within the simulation cone.";
439  }
441  }
442 
443  // Allocate sky direction
444  GSkyDir dir;
445 
446  // Debug option: initialise counter
447  #if defined(G_DEBUG_MC)
448  int num_iterations = 0;
449  #endif
450 
451  // Get sky direction
452  while (true) {
453 
454  // Debug option: increment counter
455  #if defined(G_DEBUG_MC)
456  num_iterations++;
457  #endif
458 
459  // Simulate random sky direction within Monte Carlo simulation cone
460  double theta = std::acos(1.0 - ran.uniform() * m_mc_one_minus_cosrad) *
462  double phi = 360.0 * ran.uniform();
463  dir = m_mc_cone.centre();
464  dir.rotate_deg(phi, theta);
465 
466  // Get map value at simulated sky direction
467  double value = 0.0;
468  double value_left = m_cube(dir, m_logE.inx_left());
469  double value_right = m_cube(dir, m_logE.inx_right());
470  if (value_left > 0.0 && value_right > 0.0) {
471  double value_log = m_logE.wgt_left() * std::log(value_left) +
472  m_logE.wgt_right() * std::log(value_right);
473  value = std::exp(value_log);
474  }
475  else if (value_left > 0.0) {
476  value = value_left;
477  }
478  else if (value_right > 0.0) {
479  value = value_right;
480  }
481 
482  // If map value is non-positive then simulate a new sky direction
483  if (value <= 0.0) {
484  continue;
485  }
486 
487  // Get uniform random number
488  double uniform = ran.uniform() * max;
489 
490  // Exit loop if the random number is not larger than the map cube value
491  if (uniform <= value) {
492  break;
493  }
494 
495  } // endwhile: loop until sky direction was accepted
496 
497  // Debug option: log counter
498  #if defined(G_DEBUG_MC)
499  std::cout << num_iterations << " ";
500  #endif
501 
502  // Return sky direction
503  return dir;
504 }
505 
506 
507 /***********************************************************************//**
508  * @brief Signals whether model contains sky direction
509  *
510  * @param[in] dir Sky direction.
511  * @param[in] margin Margin to be added to sky direction (degrees)
512  * @return True if the model contains the sky direction.
513  *
514  * Signals whether a sky direction falls within the bounding circle of
515  * the diffuse cube.
516  ***************************************************************************/
518  const double& margin) const
519 {
520  // Make sure that cube in online
521  fetch_cube();
522 
523  // Initialise containment flag
524  bool contains = false;
525 
526  // Continue only if radius is positive
527  if (m_region.radius() > 0.0) {
528 
529  // Compute distance to centre
530  double distance = m_region.centre().dist_deg(dir);
531 
532  // If distance is smaller than radius plus margin we consider
533  // the position to be contained within the bounding circle
534  if (distance < m_region.radius() + margin) {
535  contains = true;
536  }
537 
538  }
539 
540  // Return containment
541  return (contains);
542 }
543 
544 
545 /***********************************************************************//**
546  * @brief Read model from XML element
547  *
548  * @param[in] xml XML element.
549  *
550  * @exception GException::invalid_value
551  * Model parameters not found in XML element.
552  *
553  * Read the map cube information from an XML element. The XML element should
554  * have the format
555  *
556  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
557  * <parameter name="Normalization" ../>
558  * </spatialModel>
559  ***************************************************************************/
561 {
562  // Verify number of model parameters
564 
565  // If "Normalization" parameter exists then read parameter from this
566  // XML element
567  if (gammalib::xml_has_par(xml, "Normalization")) {
568  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Normalization");
569  m_value.read(*par);
570  }
571 
572  // ... otherwise try reading parameter from "Value" parameter
573  #if defined(G_LEGACY_XML_FORMAT)
574  else if (gammalib::xml_has_par(xml, "Value")) {
575  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Value");
576  m_value.read(*par);
577  }
578  #endif
579 
580  // ... otherwise throw an exception
581  else {
582  #if defined(G_LEGACY_XML_FORMAT)
583  std::string msg = "Diffuse map cube model requires either "
584  "\"Normalization\" or \"Value\" parameter.";
585  #else
586  std::string msg = "Diffuse map cube model requires \"Normalization\" "
587  "parameter.";
588  #endif
589  throw GException::invalid_value(G_READ, msg);
590  }
591 
592  // Save filename
593  m_filename = gammalib::xml_file_expand(xml, xml.attribute("file"));
594 
595  // Return
596  return;
597 }
598 
599 
600 /***********************************************************************//**
601  * @brief Write model into XML element
602  *
603  * @param[in] xml XML element.
604  *
605  * Write the map cube information into an XML element. The XML element will
606  * have either the format
607  *
608  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
609  * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
610  * </spatialModel>
611  *
612  * or alternatively
613  *
614  * <spatialModel type="DiffuseMapCube" file="test_file.fits">
615  * <parameter name="Value" scale="1" value="1" min="0.1" max="10" free="0"/>
616  * </spatialModel>
617  *
618  * The latter format is the default for newly written XML elements.
619  ***************************************************************************/
621 {
622  // Verify model type
624 
625  // Get or create parameter
627 
628  // Write parameter
629  m_value.write(*value);
630 
631  // Set filename
633 
634  // Return
635  return;
636 }
637 
638 
639 /***********************************************************************//**
640  * @brief Set map cube
641  *
642  * @param[in] cube Sky map.
643  *
644  * Set the map cube of the spatial map cube model.
645  ***************************************************************************/
647 {
648  // Clear filename and signal map as unloaded
649  m_filename.clear();
650  m_loaded = false;
651 
652  // Assign map
653  m_cube = cube;
654 
655  // Update Monte Carlo cache
656  update_mc_cache();
657 
658  // Return
659  return;
660 }
661 
662 
663 /***********************************************************************//**
664  * @brief Set energies for map cube
665  *
666  * @param[in] energies Sky map energies.
667  *
668  * @exception GException::invalid_argument
669  * Specified sky map energies incompatible with map cube.
670  *
671  * Sets the energies for the map cube.
672  ***************************************************************************/
674 {
675  // Initialise energies
676  m_logE.clear();
677 
678  // Fetch cube
679  fetch_cube();
680 
681  // Extract number of energies in vector
682  int num = energies.size();
683 
684  // Check if energy binning is consistent with number of maps in the cube
685  if (num != m_cube.nmaps() ) {
686  std::string msg = "Number of specified energies ("+gammalib::str(num)+")"
687  " does not match the number of maps ("
688  ""+gammalib::str(m_cube.nmaps())+" in the map cube.\n"
689  "The energies argument shall provide a vector of length"
690  " "+gammalib::str(m_cube.nmaps())+".";
692  }
693 
694  // Set log10(energy) nodes, where energy is in units of MeV
695  for (int i = 0; i < num; ++i) {
696  m_logE.append(energies[i].log10MeV());
697  }
698 
699  // Set energy boundaries
701 
702  // Update MC cache
703  update_mc_cache();
704 
705  // Return
706  return;
707 }
708 
709 
710 /***********************************************************************//**
711  * @brief Returns map cube energies
712  *
713  * @return Map cube energies.
714  *
715  * Returns the energies for the map cube in a vector.
716  ***************************************************************************/
718 {
719  // Initialise energies container
721 
722  // Fetch cube
723  fetch_cube();
724 
725  // Get number of map energies
726  int num = m_logE.size();
727 
728  // Continue only if there are maps in the cube
729  if (num > 0) {
730 
731  // Reserve space for all energies
732  energies.reserve(num);
733 
734  // Set log10(energy) nodes, where energy is in units of MeV
735  for (int i = 0; i < num; ++i) {
736  GEnergy energy;
737  energy.log10MeV(m_logE[i]);
738  energies.append(energy);
739  }
740 
741  } // endif: there were maps in the cube
742 
743  // Return energies
744  return energies;
745 }
746 
747 
748 /***********************************************************************//**
749  * @brief Set Monte Carlo simulation cone
750  *
751  * @param[in] cone Monte Carlo simulation cone.
752  *
753  * Sets the simulation cone that defines the directions that will be
754  * simulated using the mc() method and pre-computes the maximum intensity
755  * and the spatially integrated flux of each map within the simulation cone
756  * region.
757  ***************************************************************************/
759 {
760  // Continue only if the simulation cone has changed
761  if (cone != m_mc_cone) {
762 
763  // Save simulation cone definition
764  m_mc_cone = cone;
765 
766  // Pre-compute 1 - cosine of radius
768 
769  // Initialise Monte Carlo cache
770  m_mc_max.clear();
772 
773  // Fetch cube
774  fetch_cube();
775 
776  // Determine number of cube pixels and maps
777  int npix = cube().npix();
778  int nmaps = cube().nmaps();
779 
780  // Continue only if there are pixels and maps
781  if (npix > 0 && nmaps > 0) {
782 
783  // Initial vectors for total flux and maximum intensity
784  std::vector<double> total_flux(nmaps, 0.0);
785  std::vector<double> max_intensity(nmaps, 0.0);
786 
787  // Loop over all map pixels
788  for (int k = 0; k < npix; ++k) {
789 
790  // Continue only if pixel is within MC cone
791  double distance = m_mc_cone.centre().dist_deg(m_cube.pix2dir(k));
792  if (distance <= m_mc_cone.radius()) {
793 
794  // Loop over all maps
795  for (int i = 0; i < nmaps; ++i) {
796 
797  // Update total flux
798  double flux = m_cube.flux(k,i);
799  if (flux > 0.0) {
800  total_flux[i] += flux;
801  }
802 
803  // Update maximum intensity
804  double value = m_cube(k,i);
805  if (value > max_intensity[i]) {
806  max_intensity[i] = value;
807  }
808 
809  } // endfor: looped over all maps
810 
811  } // endif: pixel within MC cone
812 
813  } // endfor: looped over all pixels
814 
815  // Reserve space in cache
816  m_mc_max.reserve(nmaps);
817  m_mc_spectrum.reserve(nmaps);
818 
819  // Set maximum intensity and spectrum vectors
820  for (int i = 0; i < nmaps; ++i) {
821 
822  // Store maximum map intensity
823  m_mc_max.push_back(max_intensity[i]);
824 
825  // Store flux as spectral node
826  if (m_logE.size() == nmaps) {
827 
828  // Set map energy
829  GEnergy energy;
830  energy.log10MeV(m_logE[i]);
831 
832  // Only append node if the integrated flux is positive
833  // (as we do a log-log interpolation)
834  if (total_flux[i] > 0.0) {
835  m_mc_spectrum.append(energy, total_flux[i]);
836  }
837 
838  }
839 
840  } // endfor: looped over all maps
841 
842  // Log maximum intensity and total flux for debugging
843  #if defined(G_DEBUG_MC_CACHE)
844  std::cout << "GModelSpatialDiffuseCube::mc_cone:" << std::endl;
845  std::cout << " Maximum map intensity:" << std::endl;
846  for (int i = 0; i < m_mc_max.size(); ++i) {
847  GEnergy energy;
848  energy.log10MeV(m_logE[i]);
849  std::cout << " " << i;
850  std::cout << " " << energy;
851  std::cout << " " << m_mc_max[i] << " ph/cm2/s/sr/MeV";
852  std::cout << std::endl;
853  }
854  std::cout << " Spatially integrated flux:" << std::endl;
855  for (int i = 0; i < m_mc_spectrum.nodes(); ++i) {
856  std::cout << " " << i;
857  std::cout << " " << m_mc_spectrum.energy(i);
858  std::cout << " " << m_mc_spectrum.intensity(i);
859  std::cout << " ph/cm2/s/MeV" << std::endl;
860  }
861  #endif
862 
863  } // endif: there were cube pixels and maps
864 
865  } // endif: simulation cone has changed
866 
867  // Return
868  return;
869 }
870 
871 
872 /***********************************************************************//**
873  * @brief Load cube into the model class
874  *
875  * @param[in] filename cube file.
876  *
877  * Loads cube into the model class.
878  ***************************************************************************/
880 {
881  // Load cube
882  load_cube(filename);
883 
884  // Update Monte Carlo cache
885  update_mc_cache();
886 
887  // Return
888  return;
889 }
890 
891 
892 /***********************************************************************//**
893  * @brief Save cube into FITS file
894  *
895  * @param[in] filename Cube FITS file name.
896  * @param[in] clobber Overwrite existing FITS file (default: false).
897  *
898  * Saves spatial cube model into a FITS file.
899  ***************************************************************************/
901  const bool& clobber) const
902 {
903  // Create empty FITS file
904  GFits fits;
905 
906  // Write event cube
907  write(fits);
908 
909  // Save FITS file
910  fits.saveto(filename, clobber);
911 
912  // Return
913  return;
914 }
915 
916 
917 /***********************************************************************//**
918  * @brief Read diffuse cube from FITS file
919  *
920  * @param[in] fits FITS file.
921  *
922  * @exception GException::invalid_value
923  * ENERGIES extension incompatible with map cube.
924  *
925  * Reads the diffuse cube sky map from the first extension in the @p fits
926  * file and the energy extension from the "ENERGIES" extension.
927  ***************************************************************************/
929 {
930  // Initialise skymap
931  m_cube.clear();
932  m_logE.clear();
933 
934  // Read cube from the first valid WCS image
935  for (int extno = 0; extno < fits.size(); ++extno) {
936  const GFitsHDU& hdu = *fits.at(extno);
937  if ((hdu.exttype() == GFitsHDU::HT_IMAGE) &&
938  (hdu.has_card("NAXIS") && (hdu.integer("NAXIS") >= 2))) {
939  m_cube.read(hdu);
940  break;
941  }
942  }
943 
944  // Read energies
946  energies.read(*(fits.table(gammalib::extname_energies)));
947 
948  // Extract number of energy bins
949  int num = energies.size();
950 
951  // Check if energy binning is consistent with primary image hdu
952  if (num != m_cube.nmaps() ) {
953  std::string msg = "Number of energies in \"ENERGIES\" extension ("+
954  gammalib::str(num)+") does not match the number of "
955  "maps ("+gammalib::str(m_cube.nmaps())+" in the "
956  "map cube. The \"ENERGIES\" extension table shall "
957  "provide one enegy value for each map in the cube.";
959  }
960 
961  // Set log10(energy) nodes, where energy is in units of MeV
962  for (int i = 0; i < num; ++i) {
963  m_logE.append(energies[i].log10MeV());
964  }
965 
966  // Signal that cube has been loaded
967  m_loaded = true;
968 
969  // Set energy boundaries
971 
972  // Return
973  return;
974 }
975 
976 
977 /***********************************************************************//**
978  * @brief Write diffuse cube into FITS file
979  *
980  * @param[in] fits FITS file.
981  *
982  * Writes the diffuse cube sky map and the energy intervals into @p fits
983  * file.
984  ***************************************************************************/
986 {
987  // Make sure that cube is fetched
988  fetch_cube();
989 
990  // Write cube
991  m_cube.write(fits);
992 
993  // Create energy intervals
995  for (int i = 0; i < m_logE.size(); ++i) {
996  double energy = std::pow(10.0, m_logE[i]);
997  energies.append(GEnergy(energy, "MeV"));
998  }
999 
1000  // Write energy intervals
1001  energies.write(fits);
1002 
1003  // Return
1004  return;
1005 }
1006 
1007 
1008 /***********************************************************************//**
1009  * @brief Returns diffuse cube flux integrated in sky region
1010  *
1011  * @param[in] region Sky region.
1012  * @param[in] srcEng Energy.
1013  * @param[in] srcTime Time.
1014  * @return Flux (adimensional or ph/cm2/s).
1015  *
1016  * Returns diffuse cube flux within a sky region. The flux is computed by
1017  * log-log interpolation of the flux of the sky maps that encompass the
1018  * specified @p srcEng energy. Specifically
1019  *
1020  * \f[
1021  * F = {\tt value} \times \left \{
1022  * \begin{array}{l l}
1023  * 0 & \mbox{if $F_l=0$ and $F_r=0$} \\
1024  * F_l & \mbox{if $F_l\ne0$ and $F_r=0$} \\
1025  * F_r & \mbox{if $F_l=0$ and $F_r\ne0$} \\
1026  * \exp \left( w_l \log F_l + w_r \log F_r \right) & \mbox{else}
1027  * \end{array}
1028  * \right .
1029  * \f]
1030  *
1031  * where
1032  * - \f$F\f$ is the diffuse cube flux
1033  * - \f${\tt value}\f$ is the cube normalisation that is returned by the
1034  * value() method,
1035  * - \f$w_l\f$ and \f$w_r\f$ are the weighting factors from the logarithmic
1036  * interpolation in @p srcEng, with \f$w_l+w_r=1\f$, and
1037  * - \f$F_l\f$ and \f$F_r\f$ are the fluxes in the sky maps that encompass
1038  * the @p srcEng energy.
1039  *
1040  * This method fetches the diffuse cube if it was not yet fetched before.
1041  ***************************************************************************/
1043  const GEnergy& srcEng,
1044  const GTime& srcTime) const
1045 {
1046  // Fetch cube
1047  fetch_cube();
1048 
1049  // Initialise flux value
1050  double flux = 0.0;
1051 
1052  // Continue only if there is energy information for the map cube
1053  if (m_logE.size() > 0) {
1054 
1055  // Set energy for interpolation in log-energy
1056  m_logE.set_value(srcEng.log10MeV());
1057 
1058  // Compute map cube flux for the left and right map
1059  double left_flux = m_cube.flux(region, m_logE.inx_left());
1060  double right_flux = m_cube.flux(region, m_logE.inx_right());
1061 
1062  // Perform log-log interpolation
1063  if (left_flux > 0.0 && right_flux > 0.0) {
1064  double log_flux = m_logE.wgt_left() * std::log(left_flux) +
1065  m_logE.wgt_right() * std::log(right_flux);
1066  flux = std::exp(log_flux);
1067  }
1068  else if (left_flux > 0.0) {
1069  flux = left_flux;
1070  }
1071  else if (right_flux > 0.0) {
1072  flux = right_flux;
1073  }
1074 
1075  // Multiply by cube normalisation
1076  flux *= m_value.value();
1077 
1078  } // endif: energy information was available
1079 
1080  // Return flux
1081  return flux;
1082 }
1083 
1084 
1085 /***********************************************************************//**
1086  * @brief Print map cube information
1087  *
1088  * @param[in] chatter Chattiness.
1089  * @return String containing model information.
1090  ***************************************************************************/
1091 std::string GModelSpatialDiffuseCube::print(const GChatter& chatter) const
1092 {
1093  // Initialise result string
1094  std::string result;
1095 
1096  // Continue only if chatter is not silent
1097  if (chatter != SILENT) {
1098 
1099  // Append header
1100  result.append("=== GModelSpatialDiffuseCube ===");
1101 
1102  // Append parameters
1103  result.append("\n"+gammalib::parformat("Map cube file"));
1104  result.append(m_filename);
1105  if (m_loaded) {
1106  result.append(" [loaded]");
1107  }
1108  result.append("\n"+gammalib::parformat("Number of parameters"));
1109  result.append(gammalib::str(size()));
1110  for (int i = 0; i < size(); ++i) {
1111  result.append("\n"+m_pars[i]->print(chatter));
1112  }
1113 
1114  // Append detailed information only if a map cube exists
1115  if (m_cube.npix() > 0) {
1116 
1117  // NORMAL: Append sky map
1118  if (chatter >= NORMAL) {
1119  result.append("\n"+m_cube.print(chatter));
1120  }
1121 
1122  // EXPLICIT: Append energy nodes
1123  if (chatter >= EXPLICIT && m_logE.size() > 0) {
1124  result.append("\n"+gammalib::parformat("Map energy values"));
1125  if (m_logE.size() > 0) {
1126  for (int i = 0; i < m_logE.size(); ++i) {
1127  result.append("\n"+gammalib::parformat(" Map "+gammalib::str(i+1)));
1128  result.append(gammalib::str(std::pow(10.0, m_logE[i])));
1129  result.append(" MeV (log10E=");
1130  result.append(gammalib::str(m_logE[i]));
1131  result.append(")");
1132  if (m_ebounds.size() == m_logE.size()) {
1133  result.append(" [");
1134  result.append(m_ebounds.emin(i).print());
1135  result.append(", ");
1136  result.append(m_ebounds.emax(i).print());
1137  result.append("]");
1138  }
1139  }
1140  }
1141  else {
1142  result.append("not specified");
1143  }
1144  }
1145 
1146  // VERBOSE: Append MC cache
1147  if (chatter >= VERBOSE) {
1148  result.append("\n"+gammalib::parformat("Map flux"));
1149  if (m_mc_spectrum.nodes() > 0) {
1150  for (int i = 0; i < m_mc_spectrum.nodes(); ++i) {
1151  result.append("\n"+gammalib::parformat(" Map "+gammalib::str(i+1)));
1152  result.append(gammalib::str(m_mc_spectrum.intensity(i)));
1153  }
1154  }
1155  else {
1156  result.append("not specified");
1157  }
1158  }
1159 
1160  } // endif: map cube exists
1161 
1162  } // endif: chatter was not silent
1163 
1164  // Return result
1165  return result;
1166 }
1167 
1168 
1169 /*==========================================================================
1170  = =
1171  = Private methods =
1172  = =
1173  ==========================================================================*/
1174 
1175 /***********************************************************************//**
1176  * @brief Initialise class members
1177  ***************************************************************************/
1179 {
1180  // Initialise model type
1181  m_type = "DiffuseMapCube";
1182 
1183  // Initialise Value
1184  m_value.clear();
1185  m_value.name("Normalization");
1186  m_value.value(1.0);
1187  m_value.scale(1.0);
1188  m_value.range(0.001, 1000.0);
1189  m_value.gradient(0.0);
1190  m_value.fix();
1191  m_value.has_grad(true);
1192 
1193  // Set parameter pointer(s)
1194  m_pars.clear();
1195  m_pars.push_back(&m_value);
1196 
1197  // Initialise other members
1198  m_filename.clear();
1199  m_cube.clear();
1200  m_logE.clear();
1201  m_ebounds.clear();
1202  m_loaded = false;
1203 
1204  // Initialise MC cache
1205  m_mc_cone.clear();
1206  m_mc_one_minus_cosrad = 1.0;
1207  m_mc_max.clear();
1208  m_mc_spectrum.clear();
1209 
1210  // Return
1211  return;
1212 }
1213 
1214 
1215 /***********************************************************************//**
1216  * @brief Copy class members
1217  *
1218  * @param[in] model Spatial map cube model.
1219  ***************************************************************************/
1221 {
1222  // Copy members
1223  m_type = model.m_type; // Needed to conserve model type
1224  m_value = model.m_value;
1225  m_filename = model.m_filename;
1226  m_cube = model.m_cube;
1227  m_logE = model.m_logE;
1228  m_ebounds = model.m_ebounds;
1229  m_loaded = model.m_loaded;
1230 
1231  // Copy MC cache
1232  m_mc_cone = model.m_mc_cone;
1234  m_mc_max = model.m_mc_max;
1235  m_mc_spectrum = model.m_mc_spectrum;
1236 
1237  // Set parameter pointer(s)
1238  m_pars.clear();
1239  m_pars.push_back(&m_value);
1240 
1241  // Return
1242  return;
1243 }
1244 
1245 
1246 /***********************************************************************//**
1247  * @brief Delete class members
1248  ***************************************************************************/
1250 {
1251  // Return
1252  return;
1253 }
1254 
1255 
1256 /***********************************************************************//**
1257  * @brief Fetch cube
1258  *
1259  * Load diffuse cube if it is not yet loaded. The loading is thread save.
1260  ***************************************************************************/
1262 {
1263  // Load cube if it is not yet loaded
1264  if (!m_loaded && !m_filename.is_empty()) {
1265 
1266  // Put in a OMP critical zone
1267  #pragma omp critical(GModelSpatialDiffuseCube_fetch_cube)
1268  {
1269  const_cast<GModelSpatialDiffuseCube*>(this)->load_cube(m_filename);
1270  } // end of pragma
1271 
1272  } // endif: file not
1273 
1274  // Return
1275  return;
1276 }
1277 
1278 
1279 /***********************************************************************//**
1280  * @brief Load map cube
1281  *
1282  * @param[in] filename Map cube file.
1283  *
1284  * @exception GException::invalid_value
1285  * Number of maps in cube mismatches number of energy bins.
1286  *
1287  * Load diffuse map cube.
1288  ***************************************************************************/
1290 {
1291  // Initialise sky map cube, region and log(energy) nodes
1292  m_cube.clear();
1293  m_region.clear();
1294  m_logE.clear();
1295 
1296  // Store filename of cube (for XML writing). Note that we do
1297  // not expand any environment variable at this level, so that
1298  // if we write back the XML element we write the filepath with
1299  // the environment variables
1300  m_filename = filename;
1301 
1302  // Load cube
1303  m_cube.load(filename);
1304 
1305  // Load energies
1306  GEnergies energies(filename);
1307 
1308  // Extract number of energy bins
1309  int num = energies.size();
1310 
1311  // Check if energy binning is consistent with primary image hdu
1312  if (num != m_cube.nmaps() ) {
1313  std::string msg = "Number of energies in \"ENERGIES\" extension ("+
1314  gammalib::str(num)+") does not match the number of "
1315  "maps ("+gammalib::str(m_cube.nmaps())+" in the "
1316  "map cube. The \"ENERGIES\" extension table shall "
1317  "provide one enegy value for each map in the cube.";
1319  }
1320 
1321  // Set log10(energy) nodes, where energy is in units of MeV
1322  for (int i = 0; i < num; ++i) {
1323  m_logE.append(energies[i].log10MeV());
1324  }
1325 
1326  // Signal that cube has been loaded
1327  m_loaded = true;
1328 
1329  // Set energy boundaries
1331 
1332  // Get region circle
1334 
1335  // Set simulation cone
1336  mc_cone(m_region);
1337 
1338  // Return
1339  return;
1340 }
1341 
1342 
1343 /***********************************************************************//**
1344  * @brief Set energy boundaries
1345  *
1346  * Computes energy boundaries from the energy nodes. The boundaries will be
1347  * computed as the centre in log10(energy) between the nodes. If only a
1348  * single map exists, no boundaries will be computed.
1349  ***************************************************************************/
1351 {
1352  // Initialise energy boundaries
1353  m_ebounds.clear();
1354 
1355  // Determine number of energy bins
1356  int num = m_logE.size();
1357 
1358  // Continue only if there are at least two energy nodes
1359  if (num > 1) {
1360 
1361  // Loop over all nodes
1362  for (int i = 0; i < num; ++i) {
1363 
1364  // Calculate minimum energy
1365  double e_min = (i == 0) ? m_logE[i] - 0.5 * (m_logE[i+1] - m_logE[i])
1366  : m_logE[i] - 0.5 * (m_logE[i] - m_logE[i-1]);
1367 
1368  // Calculate maximum energy
1369  double e_max = (i < num-1) ? m_logE[i] + 0.5 * (m_logE[i+1] - m_logE[i])
1370  : m_logE[i] + 0.5 * (m_logE[i] - m_logE[i-1]);
1371 
1372  // Set energy boundaries
1373  GEnergy emin;
1374  GEnergy emax;
1375  emin.log10MeV(e_min);
1376  emax.log10MeV(e_max);
1377 
1378  // Append energy bin to energy boundary arra
1379  m_ebounds.append(emin, emax);
1380 
1381  } // endfor: looped over energy nodes
1382 
1383  } // endif: there were at least two energy nodes
1384 
1385  // Return
1386  return;
1387 }
1388 
1389 
1390 /***********************************************************************//**
1391  * @brief Update Monte Carlo cache
1392  *
1393  * Initialise the cache for Monte Carlo sampling of the map cube. See the
1394  * mc_cone() for more information.
1395  ***************************************************************************/
1397 {
1398  // Initialise sky region to all sky (dummy)
1399  GSkyRegionCircle cone(0.0, 0.0, 360.0);
1400 
1401  // Set cone
1402  mc_cone(cone);
1403 
1404  // Return
1405  return;
1406 }
1407 
1408 
1409 /***********************************************************************//**
1410  * @brief Compute cube intensity by log-log interpolation
1411  *
1412  * @param[in] photon Incident photon.
1413  * @return Cube intensity value.
1414  ***************************************************************************/
1416 {
1417  // Fetch cube
1418  fetch_cube();
1419 
1420  // Initialise intensity value
1421  double intensity = 0.0;
1422 
1423  // Continue only if there is energy information for the map cube
1424  if (m_logE.size() > 0) {
1425 
1426  // Set energy for interpolation in log-energy
1427  m_logE.set_value(photon.energy().log10MeV());
1428 
1429  // Compute map cube intensity for the left and right map
1430  double left_intensity = m_cube(photon.dir(), m_logE.inx_left());
1431  double right_intensity = m_cube(photon.dir(), m_logE.inx_right());
1432 
1433  // Perform log-log interpolation
1434  if (left_intensity > 0.0 && right_intensity > 0.0) {
1435  double log_intensity = m_logE.wgt_left() * std::log(left_intensity) +
1436  m_logE.wgt_right() * std::log(right_intensity);
1437  intensity = std::exp(log_intensity);
1438  }
1439  else if (left_intensity > 0.0) {
1440  intensity = left_intensity;
1441  }
1442  else if (right_intensity > 0.0) {
1443  intensity = right_intensity;
1444  }
1445 
1446  } // endif: energy information was available
1447 
1448  // Return intensity
1449  return intensity;
1450 }
1451 
1452 
1453 /***********************************************************************//**
1454  * @brief Set boundary sky region
1455  ***************************************************************************/
1457 {
1458  // Make sure that cube in online
1459  fetch_cube();
1460 
1461  // Return
1462  return;
1463 }
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:286
double flux(const int &index, const int &map=0) const
Returns flux in pixel.
Definition: GSkyMap.cpp:1551
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:615
virtual void set_region(void) const
Set boundary sky region.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
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:659
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.
virtual double flux(const GSkyRegion &region, const GEnergy &srcEng=GEnergy(), const GTime &srcTime=GTime()) const
Returns diffuse cube flux integrated in sky region.
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:163
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:301
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
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:48
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.
Time class.
Definition: GTime.hpp:55
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:587
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:1946
Interface definition for the spatial model registry class.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2664
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:1293
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2787
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:1360
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void init_members(void)
Initialise class members.
Abstract interface for the sky region class.
Definition: GSkyRegion.hpp:57
const 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:391
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:2697
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:424
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:187
void fix(void)
Fix a parameter.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:170
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1637
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
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.
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:1100
#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:268
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:915
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
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:1596
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:2401
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.
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:1422
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:1232
void init_members(void)
Initialise class members.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:199
Abstract diffuse spatial model base class.
const double rad2deg
Definition: GMath.hpp:44
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void 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:2511
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:347
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:1689
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:1889
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
FITS table abstract base class interface definition.
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819