GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMObservation.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMObservation.cpp - COMPTEL Observation class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-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 GCOMObservation.cpp
23  * @brief COMPTEL observation class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <typeinfo>
32 #include "GObservationRegistry.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GException.hpp"
36 #include "GFits.hpp"
37 #include "GFitsHDU.hpp"
38 #include "GCaldb.hpp"
39 #include "GSource.hpp"
40 #include "GModels.hpp"
41 #include "GModelSky.hpp"
42 #include "GModelSpectralConst.hpp"
43 #include "GXmlElement.hpp"
44 #include "GCOMObservation.hpp"
45 #include "GCOMEventBin.hpp"
46 #include "GCOMStatus.hpp"
47 #include "GCOMSupport.hpp"
48 
49 /* __ Globals ____________________________________________________________ */
51 const GObservationRegistry g_obs_com_registry(&g_obs_com_seed);
52 
53 /* __ Method name definitions ____________________________________________ */
54 #define G_RESPONSE "GCOMObservation::response(GResponse&)"
55 #define G_READ "GCOMObservation::read(GXmlElement&)"
56 #define G_WRITE "GCOMObservation::write(GXmlElement&)"
57 #define G_LOAD_DRB "GCOMObservation::load_drb(GFilename&)"
58 #define G_LOAD_DRG "GCOMObservation::load_drg(GFilename&)"
59 #define G_DRM "GCOMObservation::drm(GModels&)"
60 #define G_ADD_DRM "GCOMObservation::add_drm(GSource&)"
61 
62 /* __ Macros _____________________________________________________________ */
63 
64 /* __ Coding definitions _________________________________________________ */
65 
66 /* __ Debug definitions __________________________________________________ */
67 
68 
69 /*==========================================================================
70  = =
71  = Constructors/destructors =
72  = =
73  ==========================================================================*/
74 
75 /***********************************************************************//**
76  * @brief Void constructor
77  *
78  * Creates an empty COMPTEL observation.
79  ***************************************************************************/
81 {
82  // Initialise members
83  init_members();
84 
85  // Return
86  return;
87 }
88 
89 
90 /***********************************************************************//**
91  * @brief XML constructor
92  *
93  * @param[in] xml XML element.
94  *
95  * Constructs a COMPTEL observation from the information that is found in an
96  * XML element.
97  ***************************************************************************/
99 {
100  // Initialise members
101  init_members();
102 
103  // Read XML
104  read(xml);
105 
106  // Return
107  return;
108 }
109 
110 
111 /***********************************************************************//**
112  * @brief Binned observation constructor
113  *
114  * @param[in] drename Event cube name.
115  * @param[in] drbname Background cube name.
116  * @param[in] drgname Geometry cube name.
117  * @param[in] drxname Exposure map name.
118  *
119  * Creates COMPTEL observation by loading the following FITS files:
120  *
121  * DRE - Events cube
122  * DRB - Background model cube
123  * DRG - Geometry factors cube
124  * DRX - Exposure map
125  *
126  * Each of the four files is mandatory.
127  ***************************************************************************/
129  const GFilename& drbname,
130  const GFilename& drgname,
131  const GFilename& drxname) : GObservation()
132 {
133  // Initialise members
134  init_members();
135 
136  // Load observation
137  load(drename, drbname, drgname, drxname);
138 
139  // Return
140  return;
141 }
142 
143 
144 /***********************************************************************//**
145  * @brief Unbinned observation constructor
146  *
147  * @param[in] evpname Event list FITS file name.
148  * @param[in] timname Good Time Intervals FITS file name.
149  * @param[in] oadnames List of Orbit Aspect Data FITS file names.
150  *
151  * Creates a COMPTEL unbinned observation by loading the event list, Good
152  * Time Interval and Orbit Aspect Data from FITS files. All files are
153  * mandatory.
154  ***************************************************************************/
156  const GFilename& timname,
157  const std::vector<GFilename>& oadnames)
158 {
159  // Initialise members
160  init_members();
161 
162  // Load observation
163  load(evpname, timname, oadnames);
164 
165  // Return
166  return;
167 }
168 
169 
170 /***********************************************************************//**
171  * @brief Copy constructor
172  *
173  * @param[in] obs COMPTEL observation.
174  *
175  * Creates COMPTEL observation by copying an existing COMPTEL observation.
176  ***************************************************************************/
178 {
179  // Initialise members
180  init_members();
181 
182  // Copy members
183  copy_members(obs);
184 
185  // Return
186  return;
187 }
188 
189 
190 /***********************************************************************//**
191  * @brief Destructor
192  ***************************************************************************/
194 {
195  // Free members
196  free_members();
197 
198  // Return
199  return;
200 }
201 
202 
203 /*==========================================================================
204  = =
205  = Operators =
206  = =
207  ==========================================================================*/
208 
209 /***********************************************************************//**
210  * @brief Assignment operator
211  *
212  * @param[in] obs COMPTEL observation.
213  * @return COMPTEL observation.
214  *
215  * Assign COMPTEL observation to this object.
216  ***************************************************************************/
218 {
219  // Execute only if object is not identical
220  if (this != &obs) {
221 
222  // Copy base class members
223  this->GObservation::operator=(obs);
224 
225  // Free members
226  free_members();
227 
228  // Initialise members
229  init_members();
230 
231  // Copy members
232  copy_members(obs);
233 
234  } // endif: object was not identical
235 
236  // Return this object
237  return *this;
238 }
239 
240 
241 /*==========================================================================
242  = =
243  = Public methods =
244  = =
245  ==========================================================================*/
246 
247 /***********************************************************************//**
248  * @brief Clear COMPTEL observation
249  ***************************************************************************/
251 {
252  // Free members
253  free_members();
255 
256  // Initialise members
258  init_members();
259 
260  // Return
261  return;
262 }
263 
264 
265 /***********************************************************************//**
266  * @brief Clone COMPTEL observation
267  *
268  * @return Pointer to deep copy of COMPTEL observation.
269  ***************************************************************************/
271 {
272  return new GCOMObservation(*this);
273 }
274 
275 
276 /***********************************************************************//**
277  * @brief Set response function
278  *
279  * @param[in] rsp Response function.
280  *
281  * @exception GException::invalid_argument
282  * Specified response is not a COMPTEL response.
283  *
284  * Sets the response function for the observation.
285  ***************************************************************************/
287 {
288  // Get pointer on COM response
289  const GCOMResponse* comrsp = dynamic_cast<const GCOMResponse*>(&rsp);
290  if (comrsp == NULL) {
291  std::string cls = std::string(typeid(&rsp).name());
292  std::string msg = "Response of type \""+cls+"\" is "
293  "not a COMPTEL response. Please specify a COMPTEL "
294  "response as argument.";
296  }
297 
298  // Clone response function
299  m_response = *comrsp;
300 
301  // Return
302  return;
303 }
304 
305 
306 /***********************************************************************//**
307  * @brief Set response function
308  *
309  * @param[in] caldb Calibration database.
310  * @param[in] rspname Name of COMPTEL response.
311  *
312  * Sets the response function by loading the response information from the
313  * calibration database.
314  ***************************************************************************/
315 void GCOMObservation::response(const GCaldb& caldb, const std::string& rspname)
316 {
317  // Clear COM response function
318  m_response.clear();
319 
320  // Set calibration database
321  m_response.caldb(caldb);
322 
323  // Load instrument response function
324  m_response.load(rspname);
325 
326  // Return
327  return;
328 }
329 
330 
331 /***********************************************************************//**
332  * @brief Read observation from XML element
333  *
334  * @param[in] xml XML element.
335  *
336  * Reads information for a COMPTEL observation from an XML element. The
337  * method supports both an unbinned and a binned observation.
338  *
339  * For an unbinned observation the XML format is
340  *
341  * <observation name="Crab" id="000001" instrument="COM">
342  * <parameter name="EVP" file="m16992_evp.fits"/>
343  * <parameter name="TIM" file="m10695_tim.fits"/>
344  * <parameter name="OAD" file="m20039_oad.fits"/>
345  * <parameter name="OAD" file="m20041_oad.fits"/>
346  * ...
347  * </observation>
348  *
349  * where the observation can contain an arbitrary number of OAD file
350  * parameters. The @p file attribute provide either absolute or relative
351  * file name. If a file name includes no access path it is assumed that
352  * the file resides in the same location as the XML file.
353  *
354  * For a binned observation the XML format is
355  *
356  * <observation name="Crab" id="000001" instrument="COM">
357  * <parameter name="DRE" file="m50438_dre.fits"/>
358  * <parameter name="DRB" file="m34997_drg.fits"/>
359  * <parameter name="DRG" file="m34997_drg.fits"/>
360  * <parameter name="DRX" file="m32171_drx.fits"/>
361  * <parameter name="IAQ" value="UNH(1.0-3.0)MeV"/>
362  * </observation>
363  ***************************************************************************/
365 {
366  // Clear observation
367  clear();
368 
369  // Extract instrument name
370  m_instrument = xml.attribute("instrument");
371 
372  // If the XML elements has a "EVP" attribute we have an unbinned
373  // observation
374  if (gammalib::xml_has_par(xml, "EVP")) {
375 
376  // Get EVP and TIM file names
377  std::string evpname = gammalib::xml_get_attr(G_READ, xml, "EVP", "file");
378  std::string timname = gammalib::xml_get_attr(G_READ, xml, "TIM", "file");
379 
380  // Expand EVP and TIM file names
381  evpname = gammalib::xml_file_expand(xml, evpname);
382  timname = gammalib::xml_file_expand(xml, timname);
383 
384  // Get OAD file names
385  std::vector<GFilename> oadnames;
386  for (int i = 0; i < xml.elements("parameter"); ++i) {
387  const GXmlElement* element = xml.element("parameter", i);
388  if (element->attribute("name") == "OAD") {
389  std::string oadname = element->attribute("file");
390  oadname = gammalib::xml_file_expand(xml, oadname);
391  oadnames.push_back(GFilename(oadname));
392  }
393  }
394 
395  // Load observation
396  load(evpname, timname, oadnames);
397 
398  } // endif: unbinned observation
399 
400  // ... otherwise we have a binned observation
401  else {
402 
403  // Get parameters
404  std::string drename = gammalib::xml_get_attr(G_READ, xml, "DRE", "file");
405  std::string drbname = gammalib::xml_get_attr(G_READ, xml, "DRB", "file");
406  std::string drgname = gammalib::xml_get_attr(G_READ, xml, "DRG", "file");
407  std::string drxname = gammalib::xml_get_attr(G_READ, xml, "DRX", "file");
408  std::string iaqname = gammalib::xml_get_attr(G_READ, xml, "IAQ", "value");
409 
410  // Expand file names
411  drename = gammalib::xml_file_expand(xml, drename);
412  drbname = gammalib::xml_file_expand(xml, drbname);
413  drgname = gammalib::xml_file_expand(xml, drgname);
414  drxname = gammalib::xml_file_expand(xml, drxname);
415 
416  // Load observation
417  load(drename, drbname, drgname, drxname);
418 
419  // Load IAQ
420  response(GCaldb("cgro", "comptel"), iaqname);
421 
422  } // endelse: binned observation
423 
424  // Return
425  return;
426 }
427 
428 
429 /***********************************************************************//**
430  * @brief Write observation into XML element
431  *
432  * @param[in] xml XML element.
433  *
434  * Writes information for a COMPTEL observation into an XML element. The
435  * method supports both an unbinned and a binned observation.
436  *
437  * For an unbinned observation the XML format is
438  *
439  * <observation name="Crab" id="000001" instrument="COM">
440  * <parameter name="EVP" file="m16992_evp.fits"/>
441  * <parameter name="TIM" file="m10695_tim.fits"/>
442  * <parameter name="OAD" file="m20039_oad.fits"/>
443  * <parameter name="OAD" file="m20041_oad.fits"/>
444  * ...
445  * </observation>
446  *
447  * where the observation can contain an arbitrary number of OAD file
448  * parameters. The @p file attribute provide either absolute or relative
449  * file name. If a file name includes no access path it is assumed that
450  * the file resides in the same location as the XML file.
451  *
452  * For a binned observation the XML format is
453  *
454  * <observation name="Crab" id="000001" instrument="COM">
455  * <parameter name="DRE" file="m50438_dre.fits"/>
456  * <parameter name="DRB" file="m34997_drg.fits"/>
457  * <parameter name="DRG" file="m34997_drg.fits"/>
458  * <parameter name="DRX" file="m32171_drx.fits"/>
459  * <parameter name="IAQ" value="UNH(1.0-3.0)MeV"/>
460  * </observation>
461  ***************************************************************************/
463 {
464  // Allocate XML element pointer
465  GXmlElement* par;
466 
467  // Handle unbinned observation
468  if (is_unbinned()) {
469 
470  // Set EVP parameter
471  par = gammalib::xml_need_par(G_WRITE, xml, "EVP");
472  par->attribute("file", gammalib::xml_file_reduce(xml, m_evpname));
473 
474  // Set TIM parameter
475  par = gammalib::xml_need_par(G_WRITE, xml, "TIM");
476  par->attribute("file", gammalib::xml_file_reduce(xml, m_timname));
477 
478  // Remove all existing OAD parameters
479  for (int i = 0; i < xml.elements("parameter"); ++i) {
480  GXmlElement* element = xml.element("parameter", i);
481  if (element->attribute("name") == "OAD") {
482  xml.remove(i);
483  }
484  }
485 
486  // Set OAD parameters
487  for (int i = 0; i < m_oadnames.size(); ++i) {
488  par = static_cast<GXmlElement*>(xml.append(GXmlElement("parameter name=\"OAD\"")));
489  par->attribute("file", gammalib::xml_file_reduce(xml, m_oadnames[i]));
490  }
491 
492  } // endif: observation was unbinned
493 
494  // ... otherwise handle a binned observation
495  else if (is_binned()) {
496 
497  // Set DRE parameter
498  par = gammalib::xml_need_par(G_WRITE, xml, "DRE");
499  par->attribute("file", gammalib::xml_file_reduce(xml, m_drename));
500 
501  // Set DRB parameter
502  par = gammalib::xml_need_par(G_WRITE, xml, "DRB");
503  par->attribute("file", gammalib::xml_file_reduce(xml, m_drbname));
504 
505  // Set DRG parameter
506  par = gammalib::xml_need_par(G_WRITE, xml, "DRG");
507  par->attribute("file", gammalib::xml_file_reduce(xml, m_drgname));
508 
509  // Set DRX parameter
510  par = gammalib::xml_need_par(G_WRITE, xml, "DRX");
511  par->attribute("file", gammalib::xml_file_reduce(xml, m_drxname));
512 
513  // Set IAQ parameter
514  par = gammalib::xml_need_par(G_WRITE, xml, "IAQ");
515  par->attribute("value", m_response.rspname());
516 
517  } // endif: observation was binned
518 
519  // Return
520  return;
521 }
522 
523 
524 /***********************************************************************//**
525  * @brief Load data for a binned observation
526  *
527  * @param[in] drename Event cube name.
528  * @param[in] drbname Background cube name.
529  * @param[in] drgname Geometry cube name.
530  * @param[in] drxname Exposure map name.
531  *
532  * Load event cube from DRE file, background model from DRB file, geometry
533  * factors from DRG file and the exposure map from the DRX file. All files
534  * are mandatory.
535  ***************************************************************************/
536 void GCOMObservation::load(const GFilename& drename,
537  const GFilename& drbname,
538  const GFilename& drgname,
539  const GFilename& drxname)
540 {
541  // Load DRE
542  load_dre(drename);
543 
544  // Load DRB
545  load_drb(drbname);
546 
547  // Load DRG
548  load_drg(drgname);
549 
550  // Load DRX
551  load_drx(drxname);
552 
553  // Return
554  return;
555 }
556 
557 
558 /***********************************************************************//**
559  * @brief Load data for an unbinned observation
560  *
561  * @param[in] evpname Event list FITS file name.
562  * @param[in] timname Good Time Intervals FITS file name.
563  * @param[in] oadnames List of Orbit Aspect Data FITS file names.
564  *
565  * Loads the event list, Good Time Interval and Orbit Aspect Data for an
566  * unbinned observation. All files are mandatory.
567  ***************************************************************************/
568 void GCOMObservation::load(const GFilename& evpname,
569  const GFilename& timname,
570  const std::vector<GFilename>& oadnames)
571 {
572  // Clear object
573  clear();
574 
575  // Extract observation information from event list
576  GFits fits(evpname);
577  read_attributes(fits[1]);
578  fits.close();
579 
580  // Allocate event list
581  GCOMEventList *evp = new GCOMEventList(evpname);
582  m_events = evp;
583 
584  // Load TIM data
585  m_tim.load(timname);
586 
587  // Extract ontime from TIM and compute livetime assuming 15% deadtime
588  m_ontime = m_tim.gti().ontime();
589  m_deadc = 0.85;
591 
592  // Initialise intermediate vector for OADs
593  std::vector<GCOMOads> oads;
594 
595  // Load OAD data
596  for (int i = 0; i < oadnames.size(); ++i) {
597 
598  // Load OAD file
599  GCOMOads oad(oadnames[i]);
600 
601  // Skip file if it is empty
602  if (oad.size() < 1) {
603  continue;
604  }
605 
606  // Find index of where to insert OAD file
607  int index = 0;
608  for (; index < oads.size(); ++index) {
609  if (oad[0].tstop() < oads[index][oads[index].size()-1].tstart()) {
610  break;
611  }
612  }
613 
614  // Inserts Orbit Aspect Data
615  if (index < oads.size()) {
616  oads.insert(oads.begin()+index, oad);
617  }
618  else {
619  oads.push_back(oad);
620  }
621 
622  }
623 
624  // Now put all OAD into a single container
625  for (int i = 0; i < oads.size(); ++i) {
626  m_oads.extend(oads[i]);
627  }
628 
629  // Store filenames
630  m_evpname = evpname;
631  m_timname = timname;
632  m_oadnames = oadnames;
633 
634  // Return
635  return;
636 }
637 
638 
639 /***********************************************************************//**
640  * @brief Return source model DRM cube
641  *
642  * @param[in] source Source.
643  *
644  * @return Source model DRM cube.
645  ***************************************************************************/
646 const GCOMDri& GCOMObservation::drm(const GSource& source) const
647 {
648  // Search source in DRM cache
649  int index = 0;
650  for (; index < m_drms.size(); ++index) {
651  if (source.name() == m_drms[index].name()) {
652  break;
653  }
654  }
655 
656  // If no source was found then add a cube to the DRM cache
657  if (index >= m_drms.size()) {
658  const_cast<GCOMObservation*>(this)->add_drm(source);
659  index = m_drms.size() - 1;
660  }
661 
662  // Return reference to DRM cube
663  return (m_drms[index]);
664 }
665 
666 
667 /***********************************************************************//**
668  * @brief Compute DRM cube
669  *
670  * @param[in] models Model container.
671  *
672  * @exception GException::invalid_value
673  * Observation does not contain an event cube.
674  *
675  * Computes a COMPTEL DRM cube from the information provided in a model
676  * container.
677  ***************************************************************************/
679 {
680  // Get pointer on COMPTEL event cube
681  const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(m_events);
682  if (cube == NULL) {
683  std::string cls = std::string(typeid(m_events).name());
684  std::string msg = "Events of type \""+cls+"\" is not a COMPTEL event "
685  "cube. Please specify a COMPTEL event cube when "
686  "using this method.";
687  throw GException::invalid_value(G_DRM, msg);
688  }
689 
690  // Create DRM cube as copy of DRE cube
691  GCOMEventCube drm_cube = *cube;
692 
693  // Loop over all cube bins
694  for (int i = 0; i < drm_cube.size(); ++i) {
695 
696  // Get pointer to cube bin
697  GCOMEventBin* bin = drm_cube[i];
698 
699  // Compute model value for cube bin
700  double model = models.eval(*bin, *this) * bin->size();
701 
702  // Store model value in cube bin
703  bin->counts(model);
704 
705  } // endfor: looped over DRM cube bins
706 
707  // Get a copy of the DRM
708  GCOMDri drm = drm_cube.dre();
709 
710  // Return DRM
711  return drm;
712 }
713 
714 
715 /***********************************************************************//**
716  * @brief Remove response cache for model
717  *
718  * @param[in] name Model name.
719  *
720  * Remove response cache for model @p name from response cache.
721  ***************************************************************************/
722 void GCOMObservation::remove_response_cache(const std::string& name)
723 {
724  // Search source in DRM cache
725  int index = 0;
726  for (; index < m_drms.size(); ++index) {
727  if (name == m_drms[index].name()) {
728  break;
729  }
730  }
731 
732  // If source was found then remove cache entry for source
733  if (index < m_drms.size()) {
734  m_drms.erase(m_drms.begin()+index);
735  }
736 
737  // Return
738  return;
739 }
740 
741 
742 /***********************************************************************//**
743  * @brief Print observation information
744  *
745  * @param[in] chatter Chattiness.
746  * @return String containing observation information.
747  ***************************************************************************/
748 std::string GCOMObservation::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  // Append header
757  result.append("=== GCOMObservation ===");
758 
759  // Append information
760  result.append("\n"+gammalib::parformat("Name")+name());
761  result.append("\n"+gammalib::parformat("Identifier")+id());
762  result.append("\n"+gammalib::parformat("Instrument")+instrument());
763  result.append("\n"+gammalib::parformat("Statistic")+statistic());
764  result.append("\n"+gammalib::parformat("Ontime"));
765  result.append(gammalib::str(ontime())+" sec");
766  result.append("\n"+gammalib::parformat("Livetime"));
767  result.append(gammalib::str(livetime())+" sec");
768  result.append("\n"+gammalib::parformat("Deadtime correction"));
769  result.append(gammalib::str(m_deadc));
770 
771  // Append response (if available)
772  if (response()->rspname().length() > 0) {
773  result.append("\n"+response()->print(gammalib::reduce(chatter)));
774  }
775 
776  // Append events
777  if (m_events != NULL) {
778  result.append("\n"+m_events->print(gammalib::reduce(chatter)));
779  }
780  else {
781  result.append("\n"+gammalib::parformat("Events")+"undefined");
782  }
783 
784  // Append TIM (if available)
785  if (m_tim.gti().size() > 0) {
786  result.append("\n"+m_tim.print(gammalib::reduce(chatter)));
787  }
788  else {
789  result.append("\n"+gammalib::parformat("TIM")+"undefined");
790  }
791 
792  // Append OADs (if available)
793  if (!m_oads.is_empty()) {
794  result.append("\n"+m_oads.print(gammalib::reduce(chatter)));
795  }
796  else {
797  result.append("\n"+gammalib::parformat("OADs")+"undefined");
798  }
799 
800  // Append DRB, DRG and DRX
801  //result.append("\n"+m_drb.print(chatter));
802  //result.append("\n"+m_drg.print(chatter));
803  //result.append("\n"+m_drx.print(chatter));
804 
805  } // endif: chatter was not silent
806 
807  // Return result
808  return result;
809 }
810 
811 
812 /*==========================================================================
813  = =
814  = Private methods =
815  = =
816  ==========================================================================*/
817 
818 /***********************************************************************//**
819  * @brief Initialise class members
820  ***************************************************************************/
822 {
823  // Initialise members
824  m_instrument = "COM";
825  m_response.clear();
826  m_obs_id = 0;
827  m_ontime = 0.0;
828  m_livetime = 0.0;
829  m_deadc = 0.0;
830 
831  // Initialise members for binned observation
832  m_drename.clear();
833  m_drbname.clear();
834  m_drgname.clear();
835  m_drxname.clear();
836  m_drb.clear();
837  m_drg.clear();
838  m_drx.clear();
839  m_ewidth = 0.0;
840 
841  // Initialise members for unbinned observation
842  m_evpname.clear();
843  m_timname.clear();
844  m_oadnames.clear();
845  m_tim.clear();
846  m_oads.clear();
847 
848  // Initialise members for response cache
849  m_drms.clear();
850 
851  // Return
852  return;
853 }
854 
855 
856 /***********************************************************************//**
857  * @brief Copy class members
858  *
859  * @param[in] obs COMPTEL observation.
860  ***************************************************************************/
862 {
863  // Copy members
865  m_response = obs.m_response;
866  m_obs_id = obs.m_obs_id;
867  m_ontime = obs.m_ontime;
868  m_livetime = obs.m_livetime;
869  m_deadc = obs.m_deadc;
870 
871  // Copy members for binned observation
872  m_drename = obs.m_drename;
873  m_drbname = obs.m_drbname;
874  m_drgname = obs.m_drgname;
875  m_drxname = obs.m_drxname;
876  m_drb = obs.m_drb;
877  m_drg = obs.m_drg;
878  m_drx = obs.m_drx;
879  m_ewidth = obs.m_ewidth;
880 
881  // Copy members for unbinned observation
882  m_evpname = obs.m_evpname;
883  m_timname = obs.m_timname;
884  m_oadnames = obs.m_oadnames;
885  m_tim = obs.m_tim;
886  m_oads = obs.m_oads;
887 
888  // Copy members for response cache
889  m_drms = obs.m_drms;
890 
891  // Return
892  return;
893 }
894 
895 
896 /***********************************************************************//**
897  * @brief Delete class members
898  ***************************************************************************/
900 {
901  // Return
902  return;
903 }
904 
905 
906 /***********************************************************************//**
907  * @brief Load event cube data from DRE file
908  *
909  * @param[in] drename DRE filename.
910  *
911  * Loads the event cube from a DRE file.
912  *
913  * The ontime is extracted from the Good Time Intervals. The deadtime
914  * fraction is fixed to 15%, hence DEADC=0.85. The livetime is then computed
915  * by multiplying the deadtime correction by the ontime, i.e.
916  * LIVETIME = ONTIME * DEADC.
917  ***************************************************************************/
919 {
920  // Delete any existing event container (do not call clear() as we do not
921  // want to delete the response function)
922  if (m_events != NULL) delete m_events;
923  m_events = NULL;
924 
925  // Allocate event cube
926  m_events = new GCOMEventCube;
927 
928  // Open FITS file
929  GFits fits(drename);
930 
931  // Read event cube
932  m_events->read(fits);
933 
934  // Extract ontime
935  m_ontime = m_events->gti().ontime();
936 
937  // Set fixed deadtime fraction
938  m_deadc = 0.85;
939 
940  // Compute livetime
942 
943  // Compute energy width
944  m_ewidth = m_events->emax().MeV() - m_events->emin().MeV();
945 
946  // Read additional observation attributes from primary extension
947  read_attributes(fits[0]);
948 
949  // Close FITS file
950  fits.close();
951 
952  // Store event filename
953  m_drename = drename;
954 
955  // Return
956  return;
957 }
958 
959 
960 /***********************************************************************//**
961  * @brief Load background model from DRB file
962  *
963  * @param[in] drbname DRB filename.
964  *
965  * @exception GException::invalid_value
966  * DRB data space incompatible with DRE data space.
967  *
968  * Load the background model from the primary image of the specified FITS
969  * file.
970  ***************************************************************************/
972 {
973  // Open FITS file
974  GFits fits(drbname);
975 
976  // Get image
977  const GFitsImage& image = *fits.image("Primary");
978 
979  // Load background model as sky map
980  m_drb.read(image);
981 
982  // Correct WCS projection (HEASARC data format kluge)
984 
985  // Close FITS file
986  fits.close();
987 
988  // Check map dimensions
989  if (!check_map(m_drb)) {
990  std::string msg = "DRB data cube \""+drbname+"\" incompatible with "
991  "DRE data cube \""+m_drename+"\".";
993  }
994 
995  // Store DRB filename
996  m_drbname = drbname;
997 
998  // Return
999  return;
1000 }
1001 
1002 
1003 /***********************************************************************//**
1004  * @brief Load geometry factors from DRG file
1005  *
1006  * @param[in] drgname DRG filename.
1007  *
1008  * @exception GException::invalid_value
1009  * DRG data space incompatible with DRE data space.
1010  *
1011  * Load the geometry factors from the primary image of the specified FITS
1012  * file.
1013  ***************************************************************************/
1015 {
1016  // Open FITS file
1017  GFits fits(drgname);
1018 
1019  // Get image
1020  const GFitsImage& image = *fits.image("Primary");
1021 
1022  // Load geometry factors as sky map
1023  m_drg.read(image);
1024 
1025  // Correct WCS projection (HEASARC data format kluge)
1027 
1028  // Close FITS file
1029  fits.close();
1030 
1031  // Check map dimensions
1032  if (!check_map(m_drg)) {
1033  std::string msg = "DRG data cube \""+drgname+"\" incompatible with "
1034  "DRE data cube \""+m_drename+"\".";
1036  }
1037 
1038  // Store DRG filename
1039  m_drgname = drgname;
1040 
1041  // Return
1042  return;
1043 }
1044 
1045 
1046 /***********************************************************************//**
1047  * @brief Load exposure from DRX file
1048  *
1049  * @param[in] drxname DRX filename.
1050  *
1051  * Load the exposure map from the primary image of the specified FITS file.
1052  ***************************************************************************/
1054 {
1055  // Open FITS file
1056  GFits fits(drxname);
1057 
1058  // Get HDU
1059  const GFitsImage& image = *fits.image("Primary");
1060 
1061  // Load exposure map as sky map
1062  m_drx.read(image);
1063 
1064  // Correct WCS projection (HEASARC data format kluge)
1066 
1067  // Close FITS file
1068  fits.close();
1069 
1070  // Store DRX filename
1071  m_drxname = drxname;
1072 
1073  // Return
1074  return;
1075 }
1076 
1077 
1078 /***********************************************************************//**
1079  * @brief Add DRM cube to observation response cache
1080  *
1081  * @param[in] source Source.
1082  *
1083  * @exception GException::invalid_value
1084  * Observation does not contain an event cube.
1085  *
1086  * Adds a DRM cube based on the given source to the observation response
1087  * cache.
1088  ***************************************************************************/
1090 {
1091  // Get pointer on COMPTEL event cube
1092  const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(m_events);
1093  if (cube == NULL) {
1094  std::string cls = std::string(typeid(m_events).name());
1095  std::string msg = "Events of type \""+cls+"\" is not a COMPTEL event "
1096  "cube. Please specify a COMPTEL event cube when "
1097  "using this method.";
1099  }
1100 
1101  // Initialise DRM cube based on DRE cube
1102  GCOMDri drm = cube->dre();
1103 
1104  // Setup sky model for DRM computation
1105  GModelSky model(*(source.model()), GModelSpectralConst());
1106 
1107  // Compute DRM
1108  drm.compute_drm((*this), model);
1109 
1110  // Set model name
1111  drm.name(source.name());
1112 
1113  // Push model on stack
1114  m_drms.push_back(drm);
1115 
1116  // Return
1117  return;
1118 }
1119 
1120 
1121 /***********************************************************************//**
1122  * @brief Check if sky map is compatible with event cube
1123  *
1124  * @param[in] map Sky map.
1125  * @return True if map is compatible, false otherwise.
1126  *
1127  * Compares the dimension and the WCS definition of a sky map to that of the
1128  * event cube. If both are identical, true is returned, false otherwise.
1129  ***************************************************************************/
1130 bool GCOMObservation::check_map(const GSkyMap& map) const
1131 {
1132  // Get reference to event cube map
1133  const GSkyMap& ref = dynamic_cast<GCOMEventCube*>(m_events)->dre().map();
1134 
1135  // Compare dimensions
1136  bool same_dimension = ((map.nx() == ref.nx()) &&
1137  (map.ny() == ref.ny()) &&
1138  (map.nmaps() == ref.nmaps()));
1139 
1140  // Compare projections
1141  bool same_projection = (*(map.projection()) == *(ref.projection()));
1142 
1143  // Return
1144  return (same_dimension && same_projection);
1145 }
1146 
1147 
1148 /***********************************************************************//**
1149  * @brief Read observation attributes
1150  *
1151  * @param[in] hdu FITS HDU pointer
1152  *
1153  * Reads optional attributes are
1154  *
1155  * OBS_ID - Observation identifier
1156  * OBJECT - Object
1157  *
1158  * Nothing is done if the HDU pointer is NULL.
1159  ***************************************************************************/
1161 {
1162  // Continue only if HDU is valid
1163  if (hdu != NULL) {
1164 
1165  // Read observation information
1166  m_obs_id = (hdu->has_card("OBS_ID")) ? hdu->real("OBS_ID") : 0;
1167  m_name = (hdu->has_card("OBJECT")) ? hdu->string("OBJECT") : "unknown";
1168 
1169  } // endif: HDU was valid
1170 
1171  // Return
1172  return;
1173 }
1174 
1175 
1176 /***********************************************************************//**
1177  * @brief Write observation attributes
1178  *
1179  * @param[in] hdu FITS HDU pointer
1180  *
1181  * Nothing is done if the HDU pointer is NULL.
1182  *
1183  * @todo Implement method.
1184  ***************************************************************************/
1186 {
1187  // Continue only if HDU is valid
1188  if (hdu != NULL) {
1189 
1190  } // endif: HDU was valid
1191 
1192  // Return
1193  return;
1194 }
#define G_LOAD_DRB
virtual void clear(void)
Clear COMPTEL observation.
std::vector< GCOMDri > m_drms
Convolved model cubes.
const std::string & statistic(void) const
Return optimizer statistic.
double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const
Evaluate sum of all models.
Definition: GModels.cpp:871
const std::string & rspname(void) const
Return response name.
COMPTEL instrument status class definition.
Sky map class.
Definition: GSkyMap.hpp:89
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
virtual void clear(void)
Clear instance.
const GCOMOads & oads(void) const
Return Orbit Aspect Data.
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition: GSkyMap.hpp:366
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:414
std::string m_name
Observation name.
double m_obs_id
Observation ID.
#define G_RESPONSE
GSkyMap m_drb
Background model.
XML element node class interface definition.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL Good Time Intervals.
Definition: GCOMTim.cpp:404
void copy_members(const GCOMObservation &obs)
Copy class members.
#define G_ADD_DRM
Model container class definition.
GFilename m_drbname
DRB filename.
void init_members(void)
Initialise class members.
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
double m_deadc
Deadtime correction.
GFilename m_timname
TIM filename.
std::vector< GFilename > m_oadnames
OAD filenames.
GEvents * m_events
Pointer to event container.
virtual int size(void) const
Return number of bins in event cube.
virtual double ontime(void) const
Return ontime.
GFilename m_evpname
EVP filename.
#define G_DRM
void caldb(const GCaldb &caldb)
Set calibration database.
COMPTEL observation class interface definition.
XML element node class.
Definition: GXmlElement.hpp:47
void gti(const GGti &gti)
Set Good Time Intervals.
Definition: GEvents.cpp:154
bool is_unbinned(void) const
Check whether observation is unbinned.
bool is_empty(void) const
Signals if there are no Orbit Aspect Data in container.
Definition: GCOMOads.hpp:156
#define G_WRITE
Interface for the COMPTEL instrument response function.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
void load_dre(const GFilename &drename)
Load event cube data from DRE file.
std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL Orbit Aspect Data container.
Definition: GCOMOads.cpp:492
Gammalib tools definition.
void free_members(void)
Delete class members.
int size(void) const
Return number of Orbit Aspect Data in container.
Definition: GCOMOads.hpp:141
FITS file class.
Definition: GFits.hpp:63
void clear(void)
Clear COMPTEL Orbit Aspect Data container.
Definition: GCOMOads.cpp:166
FITS file class interface definition.
const std::string & name(void) const
Return model name.
Definition: GSource.hpp:114
COMPTEL event bin class.
GCOMTim m_tim
COMPTEL Good Time Intervals.
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 GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2434
Source class definition.
Implementation of support function used by COMPTEL classes.
virtual void remove_response_cache(const std::string &name)
Remove response cache for model.
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:153
const GCOMObservation g_obs_com_seed
virtual std::string print(const GChatter &chatter=NORMAL) const
Print observation information.
COMPTEL event bin class interface definition.
void free_members(void)
Delete class members.
void load(const GFilename &drename, const GFilename &drbname, const GFilename &drgname, const GFilename &drxname)
Load data for a binned observation.
virtual const GCOMResponse * response(void) const
Return response function.
Abstract FITS extension base class definition.
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
Model container class.
Definition: GModels.hpp:150
virtual double model(const GModels &models, const GEvent &event, GVector *gradient=NULL) const
Return model value and (optionally) gradient.
Calibration database class.
Definition: GCaldb.hpp:66
virtual void read(const GFits &file)=0
const GXmlAttribute * attribute(const int &index) const
Return attribute.
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:378
virtual void write(GXmlElement &xml) const
Write observation into XML element.
virtual void read(const GXmlElement &xml)
Read observation from XML element.
Filename class.
Definition: GFilename.hpp:62
void compute_drm(const GCOMObservation &obs, const GModel &model)
Compute DRM model.
Definition: GCOMDri.cpp:741
void com_wcs_mer2car(GSkyMap &map)
Changes Mercator&#39;s projection to cartesian projection.
Definition: GCOMSupport.cpp:58
void load_drb(const GFilename &drbname)
Load background model from DRB file.
double m_ontime
Ontime (sec)
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
bool is_binned(void) const
Check whether observation is binned.
Constant spectral model class interface definition.
Interface definition for the observation registry class.
const std::string & name(void) const
Return observation name.
#define G_READ
void read_attributes(const GFitsHDU *hdu)
Read observation attributes.
Observation registry class definition.
bool check_map(const GSkyMap &map) const
Check if sky map is compatible with event cube.
virtual void clear(void)
Clear COMPTEL good time intervals.
Definition: GCOMTim.cpp:184
virtual std::string instrument(void) const
Return instrument.
GChatter
Definition: GTypemaps.hpp:33
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
double m_livetime
Livetime (sec)
const GModelSpatial * model(void) const
Return spatial model component.
Definition: GSource.hpp:128
void init_members(void)
Initialise class members.
Abstract observation base class.
virtual GCOMObservation * clone(void) const
Clone COMPTEL observation.
void add_drm(const GSource &source)
Add DRM cube to observation response cache.
COMPTEL event bin container class.
void extend(const GCOMOads &oads)
Append Orbit Aspect Data container.
Definition: GCOMOads.cpp:328
Calibration database class interface definition.
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1475
GSkyMap m_drg
Geometry factors.
const GCOMDri & dre(void) const
Return reference to DRE data.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:360
double m_ewidth
Energy width (MeV)
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
Sky model class interface definition.
virtual double counts(void) const
Return number of counts in event bin.
Sky model class.
Definition: GModelSky.hpp:120
const GEnergy & emin(void) const
Return minimum energy.
Definition: GEvents.hpp:170
GCOMOads m_oads
Orbit Aspect Data.
void load_drx(const GFilename &drxname)
Load exposure from DRX file.
Interface class for COMPTEL observations.
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
void load(const GFilename &filename, const std::string &usage="YES", const std::string &mode="NORMAL")
Load COMPTEL Good Time Intervals from FITS file.
Definition: GCOMTim.cpp:217
virtual GObservation & operator=(const GObservation &obs)
Assignment operator.
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
Exception handler interface definition.
COMPTEL event list class.
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition: GSkyMap.hpp:350
COMPTEL Data Space class.
Definition: GCOMDri.hpp:60
virtual double livetime(void) const
Return livetime.
COMPTEL Orbit Aspect Data container class.
Definition: GCOMOads.hpp:51
std::string xml_get_attr(const std::string &origin, const GXmlElement &xml, const std::string &name, const std::string &attribute)
Return attribute value for a given parameter in XML element.
Definition: GTools.cpp:1608
void write_attributes(GFitsHDU *hdu) const
Write observation attributes.
GSkyMap m_drx
Exposure map.
const std::string & name(void) const
Return DRI cube name.
Definition: GCOMDri.hpp:257
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
Abstract instrument response base class.
Definition: GResponse.hpp:67
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
virtual double size(void) const
Return size of event bin.
virtual GCOMObservation & operator=(const GCOMObservation &obs)
Assignment operator.
Class that handles gamma-ray sources.
Definition: GSource.hpp:53
GCOMObservation(void)
Void constructor.
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
const GCOMDri & drm(const GSource &source) const
Return source model DRM cube.
const GEnergy & emax(void) const
Return maximum energy.
Definition: GEvents.hpp:182
GFilename m_drxname
DRX filename.
virtual ~GCOMObservation(void)
Destructor.
const GGti & gti(void) const
Return Good Time Intervals.
Definition: GCOMTim.hpp:127
void load_drg(const GFilename &drgname)
Load geometry factors from DRG file.
const double & ontime(void) const
Returns ontime.
Definition: GGti.hpp:239
Constant spectral model class.
GCOMResponse m_response
Response functions.
GFilename m_drgname
DRG filename.
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
std::string m_instrument
Instrument name.
Mathematical function definitions.
GFilename m_drename
DRE filename.
void load(const std::string &rspname)
Load COMPTEL response.
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1683
virtual void remove(const int &index)
Remove XML child node.
Definition: GXmlNode.cpp:476
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413