GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelSpatialLookup.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelSpatialLookup.cpp - Spatial lookup table model *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2019 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 GCTAModelSpatialLookup.cpp
23  * @brief Spatial lookup table model interface implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GEnergy.hpp"
35 #include "GTime.hpp"
36 #include "GEbounds.hpp"
37 #include "GFits.hpp"
38 #include "GFitsTable.hpp"
39 #include "GFitsBinTable.hpp"
40 #include "GObservations.hpp"
41 #include "GCTAObservation.hpp"
42 #include "GCTAInstDir.hpp"
45 
46 /* __ Constants __________________________________________________________ */
47 
48 /* __ Globals ____________________________________________________________ */
50 const GCTAModelSpatialRegistry g_cta_spatial_lookup_registry(&g_cta_spatial_lookup_seed);
51 
52 /* __ Method name definitions ____________________________________________ */
53 #define G_CONSTRUCTOR "GCTAModelSpatialLookup(double&, double&, GEbounds&)"
54 #define G_READ_XML "GCTAModelSpatialLookup::read(GXmlElement&)"
55 #define G_WRITE_XML "GCTAModelSpatialLookup::write(GXmlElement&)"
56 #define G_FILL_1 "GCTAModelSpatialLookup::fill(GCTAObservation&)"
57 #define G_FILL_2 "GCTAModelSpatialLookup::fill(GObservations&)"
58 #define G_READ_TABLE "GCTAModelSpatialLookup::read(GFitsTable&)"
59 #define G_LOAD "GCTAModelSpatialLookup::load(GFilename&)"
60 #define G_FILL_BUFFER "GCTAModelSpatialLookup::fill_buffer(GCTAObservation&,"\
61  " std::vector<double>&)"
62 
63 /* __ Macros _____________________________________________________________ */
64 
65 /* __ Coding definitions _________________________________________________ */
66 
67 /* __ Debug definitions __________________________________________________ */
68 
69 
70 /*==========================================================================
71  = =
72  = Constructors/destructors =
73  = =
74  ==========================================================================*/
75 
76 /***********************************************************************//**
77  * @brief Void constructor
78  ***************************************************************************/
80 {
81  // Initialise members
82  init_members();
83 
84  // Return
85  return;
86 }
87 
88 
89 /***********************************************************************//**
90  * @brief Lookup table file constructor
91  *
92  * @param[in] filename Lookup table file name.
93  *
94  * Creates instance of a spatial lookup table model from a lookup table FITS
95  * file.
96  ***************************************************************************/
99 {
100  // Initialise members
101  init_members();
102 
103  // Load lookup table
104  load(filename);
105 
106  // Return
107  return;
108 }
109 
110 
111 /***********************************************************************//**
112  * @brief Response table constructor
113  *
114  * @param[in] table Response table.
115  *
116  * Creates instance of a spatial lookup table model from a response table.
117  ***************************************************************************/
120 {
121  // Initialise members
122  init_members();
123 
124  // Set lookup table
125  this->table(table);
126 
127  // Return
128  return;
129 }
130 
131 
132 /***********************************************************************//**
133  * @brief XML constructor
134  *
135  * @param[in] xml XML element.
136  *
137  * Creates instance of a spatial lookup table model from an XML element. See
138  * the read() method for information about the expected structure of the XML
139  * element.
140  ***************************************************************************/
143 {
144  // Initialise members
145  init_members();
146 
147  // Read information from XML element
148  read(xml);
149 
150  // Return
151  return;
152 }
153 
154 
155 /***********************************************************************//**
156  * @brief Lookup table constructor
157  *
158  * @param[in] maxrad Maximum radius (deg).
159  * @param[in] radbin Radial binsize (deg).
160  * @param[in] ebds Energy boundaries.
161  *
162  * @exception GException::invalid_argument
163  * Non positive radial binsize specified.
164  *
165  * Creates empty instance of a spatial lookup table model using the maximum
166  * lookup radius @p maxrad, the radial lookup binsize @p radbin, and some
167  * energy boundaries @p ebds. Following construction, the instance may be
168  * filled using the fill() methods.
169  *
170  * The method throws an exception if a non-positive radial bin size is
171  * specified.
172  ***************************************************************************/
174  const double& radbin,
175  const GEbounds& ebds) :
177 {
178  // Throw an exception if the radian binsize is not positive
179  if (radbin <= 0.0) {
180  std::string msg = "Non-positive radial bin size "+
181  gammalib::str(radbin)+" specified.";
183  }
184 
185  // Initialise members
186  init_members();
187 
188  // Setup energy axis vectors
189  std::vector<double> eng_lo(ebds.size());
190  std::vector<double> eng_hi(ebds.size());
191  for (int i = 0; i < ebds.size(); ++i) {
192  eng_lo[i] = ebds.emin(i).TeV();
193  eng_hi[i] = ebds.emax(i).TeV();
194  }
195 
196  // Setup radial axis vector
197  int nrad = int(maxrad/radbin+0.5);
198  std::vector<double> rad_lo(nrad);
199  std::vector<double> rad_hi(nrad);
200  for (int i = 0; i < nrad; ++i) {
201  rad_lo[i] = double(i)*radbin;
202  rad_hi[i] = double(i+1)*radbin;
203  }
204 
205  // Setup response table
206  m_lookup.append_axis(eng_lo, eng_hi, "ENERG", "TeV");
207  m_lookup.append_axis(rad_lo, rad_hi, "THETA", "deg");
208  m_lookup.append_table("BKG", "1/sr");
209 
210  // Prepare table
211  prepare_table();
212 
213  // Return
214  return;
215 }
216 
217 
218 /***********************************************************************//**
219  * @brief Copy constructor
220  *
221  * @param[in] model Lookup table spatial model.
222  ***************************************************************************/
224  GCTAModelSpatial(model)
225 {
226  // Initialise members
227  init_members();
228 
229  // Copy members
230  copy_members(model);
231 
232  // Return
233  return;
234 }
235 
236 
237 /***********************************************************************//**
238  * @brief Destructor
239  ***************************************************************************/
241 {
242  // Free members
243  free_members();
244 
245  // Return
246  return;
247 }
248 
249 
250 /*==========================================================================
251  = =
252  = Operators =
253  = =
254  ==========================================================================*/
255 
256 /***********************************************************************//**
257  * @brief Assignment operator
258  *
259  * @param[in] model Lookup table spatial model.
260  ***************************************************************************/
262 {
263  // Execute only if object is not identical
264  if (this != &model) {
265 
266  // Copy base class members
267  this->GCTAModelSpatial::operator=(model);
268 
269  // Free members
270  free_members();
271 
272  // Initialise members
273  init_members();
274 
275  // Copy members
276  copy_members(model);
277 
278  } // endif: object was not identical
279 
280  // Return
281  return *this;
282 }
283 
284 
285 /*==========================================================================
286  = =
287  = Public methods =
288  = =
289  ==========================================================================*/
290 
291 /***********************************************************************//**
292  * @brief Clear instance
293  ***************************************************************************/
295 {
296  // Free class members (base and derived classes, derived class first)
297  free_members();
299 
300  // Initialise members
302  init_members();
303 
304  // Return
305  return;
306 }
307 
308 
309 /***********************************************************************//**
310  * @brief Clone instance
311  ***************************************************************************/
313 {
314  return new GCTAModelSpatialLookup(*this);
315 }
316 
317 
318 /***********************************************************************//**
319  * @brief Evaluate function
320  *
321  * @param[in] dir Event direction.
322  * @param[in] energy Event energy.
323  * @param[in] time Event time.
324  * @param[in] gradients Compute gradients?
325  * @return Function value
326  *
327  * Evaluates the lookup table model for a given event direction and energy.
328  * The evaluation is done using bi-linear interpolation in the logarithm
329  * of the energy and the offset angle radius. The lookup table value is
330  * multiplied with a normalisation parameter. The method always returns a
331  * non-negative result.
332  *
333  * If the @p gradients flag is true the method will also compute the partial
334  * derivatives of the normalisation parameter.
335  ***************************************************************************/
337  const GEnergy& energy,
338  const GTime& time,
339  const bool& gradients) const
340 {
341  // Compute offset angle in radians
342  double offset = dir.theta();
343 
344  // Determine function value by bilinear interpolation
345  double value = ((m_lookup.elements() > 0) &&
346  (m_lookup.tables() > 0))
347  ? m_lookup(0, energy.log10TeV(), offset) : 0.0;
348 
349  // Make sure that function does not become negative
350  if (value < 0.0) {
351  value = 0.0;
352  }
353 
354  // Optionally compute partial derivatives
355  if (gradients) {
356 
357  // Compute partial derivative of the normalisation parameter
358  double g_norm = (m_norm.is_free()) ? value * m_norm.scale() : 0.0;
359 
360  // Set gradient
361  m_norm.factor_gradient(g_norm);
362 
363  } // endif: computed partial derivatives
364 
365  // Return intensity times normalization factor
366  return (value * m_norm.value());
367 }
368 
369 
370 /***********************************************************************//**
371  * @brief Read model from XML element
372  *
373  * @param[in] xml XML element.
374  *
375  * Read the lookup table model information from an XML element. The XML
376  * element needs to be of the following format:
377  *
378  * <spatialModel type="LookupTable" file="lookuptable.fits">
379  * <parameter name="Normalization" ../>
380  * </spatialModel>
381  *
382  * The @p file attribute provides the filename of the lookup table FITS file.
383  * The filename may be either an absolute filename (starting with '/') or a
384  * relative filename. If no access path is given, the file is expected to
385  * reside in the same location as the XML file.
386  ***************************************************************************/
388 {
389  // Clear model
390  clear();
391 
392  // Load lookup table
393  load(gammalib::xml_file_expand(xml, xml.attribute("file")));
394 
395  // Get parameter pointers
397 
398  // Read parameters
399  m_norm.read(*norm);
400 
401  // Return
402  return;
403 }
404 
405 
406 /***********************************************************************//**
407  * @brief Write model into XML element
408  *
409  * @param[in] xml XML element.
410  *
411  * @exception GException::invalid_value
412  * Spatial model is not of valid type.
413  *
414  * Write the lookup table model information into an XML element. The XML
415  * element will be of the following format:
416  *
417  * <spatialModel type="LookupTable" file="lookuptable.fits">
418  * <parameter name="Normalization" ../>
419  * </spatialModel>
420  *
421  * The @p file attribute provides the filename of the lookup table FITS file.
422  * The filename may be either an absolute filename (starting with '/') or a
423  * relative filename. If no access path is given, the file is expected to
424  * reside in the same location as the XML file.
425  ***************************************************************************/
427 {
428  // Set model type
429  if (xml.attribute("type") == "") {
430  xml.attribute("type", type());
431  }
432 
433  // Verify model type
434  if (xml.attribute("type") != type()) {
435  std::string msg = "Spatial model \""+xml.attribute("type")+
436  "\" is not of type \""+type()+"\".";
438  }
439 
440  // Set lookup table file name
442 
443  // Get XML parameters
445 
446  // Write parameters
447  m_norm.write(*norm);
448 
449  // Return
450  return;
451 }
452 
453 
454 /***********************************************************************//**
455  * @brief Fill lookup table with events from one CTA observation
456  *
457  * @param[in] obs CTA observation.
458  *
459  * @exception GException::invalid_value
460  * No lookup table has been defined.
461  *
462  * Fills lookup table with events from one CTA observation. The method needs
463  * that a lookup table is defined before. In case that no lookup table is
464  * defined the method throws an exception.
465  ***************************************************************************/
467 {
468  // Throw an exception if response table is not initialised
469  if (m_lookup.elements() == 0) {
470  std::string msg = "Attempting to fill a yet undefined lookup table. "
471  "Please allocate a valid lookup table before "
472  "calling this method.";
474  }
475  if (m_lookup.tables() == 0) {
476  std::string msg = "Attempting to fill non-existing lookup table. "
477  "While the lookup table axes were defined, no table "
478  "exists. Please allocate a valid lookup table before "
479  "calling this method.";
481  }
482 
483  // Allocate buffer for events
484  std::vector<double> counts(m_lookup.elements());
485 
486  // Fill counts into buffer
487  fill_buffer(obs, counts);
488 
489  // Set lookup table from buffer
490  set_from_buffer(counts);
491 
492  // Return
493  return;
494 }
495 
496 
497 /***********************************************************************//**
498  * @brief Fill lookup table with CTA events in observation container
499  *
500  * @param[in] obs Observation container.
501  *
502  * @exception GException::invalid_value
503  * No lookup table has been defined.
504  *
505  * Fills lookup table with CTA events in observation container. The method
506  * needs that a lookup table is defined before. In case that no lookup table
507  * is defined the method throws an exception.
508  ***************************************************************************/
510 {
511  // Throw an exception if response table is not initialised
512  if (m_lookup.elements() == 0) {
513  std::string msg = "Attempting to fill a yet undefined lookup table. "
514  "Please allocate a valid lookup table before "
515  "calling this method.";
517  }
518  if (m_lookup.tables() == 0) {
519  std::string msg = "Attempting to fill non-existing lookup table. "
520  "While the lookup table axes were defined, no table "
521  "exists. Please allocate a valid lookup table before "
522  "calling this method.";
524  }
525 
526  // Allocate buffer for events
527  std::vector<double> counts(m_lookup.elements());
528 
529  // Loop over all observations in container
530  for (int i = 0; i < obs.size(); ++i) {
531 
532  // Get observation and continue only if it is a CTA observation
533  const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
534 
535  // Skip observation if it's not CTA
536  if (cta == NULL) {
537  continue;
538  }
539 
540  // Skip observation if we have a binned observation
541  if (cta->eventtype() == "CountsCube") {
542  continue;
543  }
544 
545  // Fill counts into buffer
546  fill_buffer(*cta, counts);
547 
548  }
549 
550  // Set lookup table from buffer
551  set_from_buffer(counts);
552 
553  // Return
554  return;
555 }
556 
557 
558 /***********************************************************************//**
559  * @brief Assign lookup table
560  *
561  * @param[in] table Lookup table.
562  *
563  * Assigns lookup table.
564  ***************************************************************************/
566 {
567  // Assign lookup table
568  m_lookup = table;
569 
570  // Prepare table
571  prepare_table();
572 
573  // Return
574  return;
575 }
576 
577 
578 /***********************************************************************//**
579  * @brief Read lookup table from FITS table
580  *
581  * @param[in] table FITS table.
582  *
583  * @exception GException::invalid_value
584  * Response table is not two-dimensional.
585  *
586  * Reads the lookup table form the FITS @p table. The following column names
587  * are mandatory:
588  *
589  * ENERG_LO - Energy lower bin boundaries
590  * ENERG_HI - Energy upper bin boundaries
591  * THETA_LO - Offset angle lower bin boundaries
592  * THETA_HI - Offset angle upper bin boundaries
593  *
594  * After reading, the method assures that the lookup table is properly
595  * normalised, i.e. that for each energy vector the radial profile has a
596  * maximum value of one.
597  ***************************************************************************/
599 {
600  // Clear lookup table
601  m_lookup.clear();
602 
603  // Read lookup table
604  m_lookup.read(table);
605 
606  // Throw an exception if the table is not two-dimensional
607  if (m_lookup.axes() != 2) {
608  std::string msg = "Expected two-dimensional lookup table but found "+
609  gammalib::str(m_lookup.axes())+" dimensions. Please "
610  "specify a two-dimensional lookup table.";
612  }
613 
614  // Prepare table
615  prepare_table();
616 
617  // Return
618  return;
619 }
620 
621 
622 /***********************************************************************//**
623  * @brief Write lookup table into FITS binary table
624  *
625  * @param[in] table FITS binary table.
626  *
627  * Writes lookup table into a FITS binary @p table.
628  ***************************************************************************/
630 {
631  // Write lookup table
632  m_lookup.write(table);
633 
634  // Return
635  return;
636 }
637 
638 
639 /***********************************************************************//**
640  * @brief Load lookup table
641  *
642  * @param[in] filename Load lookup table file name.
643  *
644  * Loads lookup table from FITS file.
645  ***************************************************************************/
647 {
648  // Clear sigma spectrum
649  clear();
650 
651  // Continue only if filename is not empty
652  if (!filename.is_empty()) {
653 
654  // Open FITS file
655  GFits fits(filename);
656 
657  // Open FITS table
658  const GFitsTable* table = NULL;
659  if (filename.has_extname()) {
660  table = fits.table(filename.extname());
661  }
662  else if (filename.has_extno()) {
663  table = fits.table(filename.extno());
664  }
665  else {
666  table = fits.table(1);
667  }
668 
669  // If no FITS table is found then throw an exception
670  if (table == NULL) {
671  std::string msg = "No lookup table found in file \""+
672  filename.url()+"\".";
673  throw GException::file_error(G_LOAD, msg);
674  }
675 
676  // Read lookup table
677  read(*table);
678 
679  // Close FITS file
680  fits.close();
681 
682  // Store filename
683  m_filename = filename;
684 
685  } // endif: filename was not empty
686 
687  // Return
688  return;
689 }
690 
691 
692 /***********************************************************************//**
693  * @brief Save lookup table into FITS file
694  *
695  * @param[in] filename FITS file name.
696  * @param[in] clobber Overwrite existing file?
697  *
698  * Saves lookup table into a FITS file. If a file with the given @p filename
699  * does not yet exist it will be created, otherwise the method opens the
700  * existing file. The method will create a (or replace an existing)
701  * effective area extension. The extension name can be specified as part
702  * of the @p filename, or if no extension name is given, is assumed to be
703  * `RADIAL BACKGROUND LOOKUP`.
704  *
705  * An existing file will only be modified if the @p clobber flag is set to
706  * true.
707  ***************************************************************************/
708 void GCTAModelSpatialLookup::save(const GFilename& filename, const bool& clobber) const
709 {
710  // Get extension name
711  std::string extname = filename.extname(gammalib::extname_cta_spatial_lookup);
712 
713  // Open or create FITS file (without extension name since the requested
714  // extension may not yet exist in the file)
715  GFits fits(filename.url(), true);
716 
717  // Remove extension if it exists already
718  if (fits.contains(extname)) {
719  fits.remove(extname);
720  }
721 
722  // Create binary table
724 
725  // Write the lookup table
726  write(table);
727 
728  // Set binary table extension name
729  table.extname(extname);
730 
731  // Append table to FITS file
732  fits.append(table);
733 
734  // Save to file
735  fits.save(clobber);
736 
737  // Return
738  return;
739 }
740 
741 
742 /***********************************************************************//**
743  * @brief Print lookup table information
744  *
745  * @param[in] chatter Chattiness.
746  * @return String containing lookup table information.
747  ***************************************************************************/
748 std::string GCTAModelSpatialLookup::print(const GChatter& chatter) const
749 {
750  // Initialise result string
751  std::string result;
752 
753  // Continue only if chatter is not silent
754  if (chatter != SILENT) {
755 
756  // Initialise information
757  int nebins = 0;
758  int nthetabins = 0;
759  double emin = 0.0;
760  double emax = 0.0;
761  double omin = 0.0;
762  double omax = 0.0;
763 
764  // Extract information if there are axes in the response table
765  if (m_lookup.axes() > 0) {
766  nebins = m_lookup.axis_bins(m_inx_energy);
767  nthetabins = m_lookup.axis_bins(m_inx_theta);
768  emin = m_lookup.axis_lo(m_inx_energy,0);
769  emax = m_lookup.axis_hi(m_inx_energy,nebins-1);
770  omin = m_lookup.axis_lo(m_inx_theta,0);
771  omax = m_lookup.axis_hi(m_inx_theta,nthetabins-1);
772  }
773 
774  // Append header
775  result.append("=== GCTAModelSpatialLookup ===");
776 
777  // Append information
778  result.append("\n"+gammalib::parformat("Filename")+m_filename);
779  result.append("\n"+gammalib::parformat("Number of energy bins") +
780  gammalib::str(nebins));
781  result.append("\n"+gammalib::parformat("Number of offset bins") +
782  gammalib::str(nthetabins));
783  result.append("\n"+gammalib::parformat("Energy range"));
784  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
785  result.append("\n"+gammalib::parformat("Offset angle range"));
786  result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
787 
788  } // endif: chatter was not silent
789 
790  // Return result
791  return result;
792 }
793 
794 
795 /*==========================================================================
796  = =
797  = Private methods =
798  = =
799  ==========================================================================*/
800 
801 /***********************************************************************//**
802  * @brief Initialise class members
803  ***************************************************************************/
805 {
806  // Initialise members
807  m_filename.clear();
808  m_lookup.clear();
809  m_inx_energy = 0;
810  m_inx_theta = 1;
811 
812  // Initialise lookup table normalisation
813  m_norm.clear();
814  m_norm.name("Normalization");
815  m_norm.value(1.0);
816  m_norm.scale(1.0);
817  m_norm.range(0.001, 1000.0);
818  m_norm.gradient(0.0);
819  m_norm.fix();
820  m_norm.has_grad(true);
821 
822  // Set parameter pointer(s)
823  m_pars.clear();
824  m_pars.push_back(&m_norm);
825 
826  // Return
827  return;
828 }
829 
830 
831 /***********************************************************************//**
832  * @brief Copy class members
833  *
834  * @param[in] model Lookup table model.
835  ***************************************************************************/
837 {
838  // Copy members
839  m_filename = model.m_filename;
840  m_lookup = model.m_lookup;
841  m_norm = model.m_norm;
842  m_inx_energy = model.m_inx_energy;
843  m_inx_theta = model.m_inx_theta;
844 
845  // Set parameter pointer(s)
846  m_pars.clear();
847  m_pars.push_back(&m_norm);
848 
849  // Return
850  return;
851 }
852 
853 
854 /***********************************************************************//**
855  * @brief Delete class members
856  ***************************************************************************/
858 {
859  // Return
860  return;
861 }
862 
863 
864 /***********************************************************************//**
865  * @brief Prepare lookup table indices
866  *
867  * Prepares a lookup table for usage.
868  *
869  * The method sets the data members m_inx_energy and m_inx_theta which
870  * determine the axes order.
871  *
872  * The method furthermore sets the energy axis to logarithmic scale and the
873  * offset angles to radians.
874  *
875  * Finally, the method assures that the maximum values for each radial
876  * profile are one.
877  ***************************************************************************/
879 {
880  // Get mandatory indices (throw exception if not found)
881  m_inx_energy = m_lookup.axis("ENERG");
882  m_inx_theta = m_lookup.axis("THETA");
883 
884  // Set energy axis to logarithmic scale
886 
887  // Set offset angle axis to radians
889 
890  // Normalise radial profiles
891  normalise_table();
892 
893  // Return
894  return;
895 }
896 
897 
898 /***********************************************************************//**
899  * @brief Normalise lookup table
900  *
901  * Normalises lookup table so that each radial profile has a maximum value
902  * of one.
903  ***************************************************************************/
905 {
906  // Get axes dimensions
907  int energy_size = m_lookup.axis_bins(m_inx_energy);
908  int theta_size = m_lookup.axis_bins(m_inx_theta);
909 
910  // Loop over energies
911  for (int i_energy = 0; i_energy < energy_size; ++i_energy) {
912 
913  // Initialise maximum radial profile value
914  double rad_max = 0.0;
915 
916  // Determine maximum radial profile value
917  for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
918  int i = table_index(i_energy, i_theta);
919  if (m_lookup(0,i) > rad_max) {
920  rad_max = m_lookup(0,i);
921  }
922  }
923 
924  // If maximum value is positive then normalise profile
925  if (rad_max > 0.0) {
926  for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
927  int i = table_index(i_energy, i_theta);
928  m_lookup(0,i) /= rad_max;
929  }
930  }
931 
932  } // endfor: looped over all energies
933 
934  // Return
935  return;
936 }
937 
938 
939 /***********************************************************************//**
940  * @brief Return index of lookup table element
941  *
942  * @param[in] ienergy Energy index.
943  * @param[in] itheta Offset index.
944  * @return Index of lookup table element.
945  ***************************************************************************/
946 int GCTAModelSpatialLookup::table_index(const int& ienergy,
947  const int& itheta) const
948 {
949  // Set index vector
950  int inx[2];
951  inx[m_inx_energy] = ienergy;
952  inx[m_inx_theta] = itheta;
953 
954  // Compute index
955  int index = inx[0] + inx[1] * m_lookup.axis_bins(0);
956 
957  // Return index
958  return index;
959 }
960 
961 
962 /***********************************************************************//**
963  * @brief Fill buffer from events in CTA observation
964  *
965  * @param[in] obs CTA observation.
966  * @param[in,out] buffer Counts buffer.
967  ***************************************************************************/
969  std::vector<double>& buffer)
970 {
971  // Make sure that the observation holds a CTA event list. If this is
972  // not the case then throw an exception.
973  const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
974  if (events == NULL) {
975  std::string msg = "CTA Observation does not contain an event "
976  "list. An event list is needed to fill the "
977  "lookup table.";
979  }
980 
981  // Get axes dimensions
982  int energy_size = m_lookup.axis_bins(m_inx_energy);
983  int theta_size = m_lookup.axis_bins(m_inx_theta);
984 
985  // Fill buffer
986  for (int i = 0; i < events->size(); ++i) {
987 
988  // Get event
989  const GCTAEventAtom* event = (*events)[i];
990 
991  // Determine offset angle in degrees
992  GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
993  double theta = inst->theta() * gammalib::rad2deg;
994 
995  // Determine offset angle bin
996  int i_theta = 0;
997  for (; i_theta < theta_size; ++i_theta) {
998  if ((theta >= m_lookup.axis_lo(m_inx_theta,i_theta)) &&
999  (theta <= m_lookup.axis_hi(m_inx_theta,i_theta))) {
1000  break;
1001  }
1002  }
1003 
1004  // Skip event if no valid offset angle bin was found
1005  if (i_theta >= theta_size) {
1006  continue;
1007  }
1008 
1009  // Determine event energy in TeV
1010  double energy = event->energy().TeV();
1011 
1012  // Determine energy bin
1013  int i_energy = 0;
1014  for (; i_energy < energy_size; ++i_energy) {
1015  if ((energy >= m_lookup.axis_lo(m_inx_energy,i_energy)) &&
1016  (energy <= m_lookup.axis_hi(m_inx_energy,i_energy))) {
1017  break;
1018  }
1019  }
1020 
1021  // Skip event if no valid offset angle bin was found
1022  if (i_energy >= theta_size) {
1023  continue;
1024  }
1025 
1026  // Fill event into bugger
1027  int inx = table_index(i_energy, i_theta);
1028  buffer[inx] += 1.0;
1029 
1030  } // endfor: looped over all events
1031 
1032  // Return
1033  return;
1034 }
1035 
1036 
1037 /***********************************************************************//**
1038  * @brief Set lookup table from buffer
1039  *
1040  * @param[in] buffer Counts buffer.
1041  ***************************************************************************/
1042 void GCTAModelSpatialLookup::set_from_buffer(const std::vector<double>& buffer)
1043 {
1044  // Get axes dimensions
1045  int energy_size = m_lookup.axis_bins(m_inx_energy);
1046  int theta_size = m_lookup.axis_bins(m_inx_theta);
1047 
1048  // Loop over offset angle
1049  for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
1050 
1051  // Compute solid angle for offset angle bin
1052  double theta_min = m_lookup.axis_lo(m_inx_theta,i_theta) * gammalib::deg2rad;
1053  double theta_max = m_lookup.axis_hi(m_inx_theta,i_theta) * gammalib::deg2rad;
1054  double area = gammalib::pi * (theta_max*theta_max - theta_min*theta_min);
1055 
1056  // Fill energy vector into lookup table
1057  for (int i_energy = 0; i_energy < energy_size; ++i_energy) {
1058  int inx = table_index(i_energy, i_theta);
1059  m_lookup(0,inx) = buffer[inx] / area;
1060  }
1061 
1062  } // endfor: looped over offset angles
1063 
1064  // Normalise radial profiles
1065  normalise_table();
1066 
1067  // Return
1068  return;
1069 }
Abstract spatial model class.
double norm(void) const
Get lookup table model normalisation.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print lookup table information.
#define G_FILL_2
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
void append_table(const std::string &name, const std::string &unit)
Append table to response table.
const double & factor_gradient(void) const
Return parameter gradient factor.
Energy value class definition.
const std::string & name(void) const
Return parameter name.
const double pi
Definition: GMath.hpp:35
double gradient(void) const
Return parameter gradient.
virtual std::string type(void) const
Return model type.
bool has_extno(void) const
Signal if filename has an extension number.
Definition: GFilename.hpp:267
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
#define G_FILL_1
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
void prepare_table(void)
Prepare lookup table indices.
CTA event list class.
virtual GCTAModelSpatial & operator=(const GCTAModelSpatial &model)
Assignment operator.
XML element node class.
Definition: GXmlElement.hpp:47
bool has_extname(void) const
Signal if filename has an extension name.
Definition: GFilename.hpp:255
#define G_CONSTRUCTOR
void fill_buffer(const GCTAObservation &obs, std::vector< double > &buffer)
Fill buffer from events in CTA observation.
double TeV(void) const
Return energy in TeV.
Definition: GEnergy.cpp:348
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
const int & tables(void) const
Return number of tables.
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
void axis_radians(const int &axis)
Set nodes for a radians axis.
FITS file class interface definition.
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1713
void read(const GFitsTable &table)
Read response table from FITS table HDU.
#define G_LOAD
bool is_free(void) const
Signal if parameter is free.
Interface definition for the spatial model registry class.
const GCTAResponseTable & table(void) const
Return lookup table.
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
const double & scale(void) const
Return parameter scale.
#define G_READ_TABLE
void write(GFitsTable &table) const
Write response table into FITS table HDU.
int axis(const std::string &name) const
Determine index of an axis.
virtual void read(const GXmlElement &xml)
Read model from XML element.
const double deg2rad
Definition: GMath.hpp:43
void set_from_buffer(const std::vector< double > &buffer)
Set lookup table from buffer.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void init_members(void)
Initialise class members.
Energy boundaries container class.
Definition: GEbounds.hpp:60
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
const GXmlAttribute * attribute(const int &index) const
Return attribute.
GCTAModelSpatialLookup(void)
Void constructor.
double theta(void) const
Return offset angle (in radians)
void save(const GFilename &filename, const bool &clobber=false) const
Save lookup table into FITS file.
void clear(void)
Clear response table.
void load(const GFilename &filename)
Load lookup table.
Filename class.
Definition: GFilename.hpp:62
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
void free_members(void)
Delete class members.
void fix(void)
Fix a parameter.
GModelPar m_norm
Normalization factor.
GFilename m_filename
Name of lookup table.
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:1513
virtual GCTAModelSpatialLookup & operator=(const GCTAModelSpatialLookup &model)
Assignment operator.
int table_index(const int &ienergy, const int &itheta) const
Return index of lookup table element.
void init_members(void)
Initialise class members.
void clear(void)
Clear parameter.
CTA instrument direction class interface definition.
virtual void clear(void)
Clear instance.
const int & elements(void) const
Return number of elements per table.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GChatter
Definition: GTypemaps.hpp:33
const std::string extname_cta_spatial_lookup
void normalise_table(void)
Normalise lookup table.
int size(void) const
Return number of observations in container.
virtual void write(GXmlElement &xml) const
Write model into XML element.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void axis_log10(const int &axis)
Set nodes for a logarithmic (base 10) axis.
CTA observation class interface definition.
Observation container class.
#define G_FILL_BUFFER
virtual GCTAModelSpatialLookup * clone(void) const
Clone instance.
Spatial lookup table model class.
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
int axis_bins(const int &axis) const
Return number bins in an axis.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void append_axis(const std::vector< double > &axis_lo, const std::vector< double > &axis_hi, const std::string &name, const std::string &unit)
Append an axis to the response table.
#define G_READ_XML
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.
FITS binary table class.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
CTA event atom class.
const GCTAModelSpatialLookup g_cta_spatial_lookup_seed
void fill(const GCTAObservation &obs)
Fill lookup table with events from one CTA observation.
int extno(const int &defaultno=-1) const
Return extension number.
Definition: GFilename.cpp:410
virtual GEvents * events(void)
Return events.
Energy boundaries class interface definition.
std::string eventtype(void) const
Return event type string.
virtual int size(void) const
Return number of events in list.
GCTAResponseTable m_lookup
Lookup table.
FITS binary table class definition.
CTA response table class.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
const double rad2deg
Definition: GMath.hpp:44
CTA observation class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
std::vector< GModelPar * > m_pars
Parameter pointers.
virtual double eval(const GCTAInstDir &dir, const GEnergy &energy, const GTime &time, const bool &gradients=false) const
Evaluate function.
Observations container class interface definition.
Spatial lookup table model interface definition.
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
const int & axes(void) const
Return number of axes of the tables.
Time class interface definition.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1562
void free_members(void)
Delete class members.
Spatial model registry class definition.
virtual ~GCTAModelSpatialLookup(void)
Destructor.
#define G_WRITE_XML
void copy_members(const GCTAModelSpatialLookup &model)
Copy class members.
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1683
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
FITS table abstract base class interface definition.