GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMEventBin.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMEventBin.cpp - COMPTEL event bin class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2018 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCOMEventBin.cpp
23  * @brief COMPTEL event bin class implementation.
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <string>
32 #include <cmath>
33 #include "GCOMEventBin.hpp"
34 #include "GTools.hpp"
35 
36 /* __ Method name definitions ____________________________________________ */
37 #define G_DIR_GET "GCOMEventBin::dir()"
38 #define G_ENERGY_GET "GCOMEventBin::energy()"
39 #define G_TIME_GET "GCOMEventBin::time()"
40 #define G_COUNTS_GET "GCOMEventBin::counts()"
41 #define G_SOLIDANGLE_GET "GCOMEventBin::solidangle()"
42 #define G_EWIDTH_GET "GCOMEventBin::ewidth()"
43 #define G_ONTIME_GET "GCOMEventBin::ontime()"
44 #define G_DIR_SET "GCOMEventBin::dir(GCOMInstDir&)"
45 #define G_ENERGY_SET "GCOMEventBin::energy(GEnergy&)"
46 #define G_TIME_SET "GCOMEventBin::time(GTime&)"
47 #define G_COUNTS_SET "GCOMEventBin::counts(double&)"
48 #define G_SOLIDANGLE_SET "GCOMEventBin::solidangle(double&)"
49 #define G_EWIDTH_SET "GCOMEventBin::ewidth(GEnergy&)"
50 #define G_ONTIME_SET "GCOMEventBin::ontime(double&)"
51 
52 /* __ Macros _____________________________________________________________ */
53 
54 /* __ Coding definitions _________________________________________________ */
55 
56 /* __ Debug definitions __________________________________________________ */
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void constructor
67  ***************************************************************************/
69 {
70  // Initialise class members for clean destruction
71  init_members();
72 
73  // Return
74  return;
75 }
76 
77 
78 /***********************************************************************//**
79  * @brief Copy constructor
80  *
81  * @param[in] bin Event bin.
82  ***************************************************************************/
84 {
85  // Initialise class members for clean destruction
86  init_members();
87 
88  // Copy members
89  copy_members(bin);
90 
91  // Return
92  return;
93 }
94 
95 
96 /***********************************************************************//**
97  * @brief Destructor
98  ***************************************************************************/
100 {
101  // Free members
102  free_members();
103 
104  // Return
105  return;
106 }
107 
108 
109 /*==========================================================================
110  = =
111  = Operators =
112  = =
113  ==========================================================================*/
114 
115 /***********************************************************************//**
116  * @brief Assignment operator
117  *
118  * @param[in] bin Event bin.
119  * @return Event bin.
120  ***************************************************************************/
122 {
123  // Execute only if object is not identical
124  if (this != &bin) {
125 
126  // Copy base class members
127  this->GEventBin::operator=(bin);
128 
129  // Free members
130  free_members();
131 
132  // Initialise private members for clean destruction
133  init_members();
134 
135  // Copy members
136  copy_members(bin);
137 
138  } // endif: object was not identical
139 
140  // Return this object
141  return *this;
142 }
143 
144 
145 /*==========================================================================
146  = =
147  = Public methods =
148  = =
149  ==========================================================================*/
150 
151 /***********************************************************************//**
152  * @brief Clear instance
153  *
154  * This method properly resets the instance to an initial state.
155  ***************************************************************************/
157 {
158  // Free class members (base and derived classes, derived class first)
159  free_members();
160  this->GEventBin::free_members();
161  this->GEvent::free_members();
162 
163  // Initialise members
164  this->GEvent::init_members();
165  this->GEventBin::init_members();
166  init_members();
167 
168  // Return
169  return;
170 }
171 
172 
173 /***********************************************************************//**
174  * @brief Clone instance
175  *
176  * @return Pointer to deep copy of event bin.
177  ***************************************************************************/
179 {
180  return new GCOMEventBin(*this);
181 }
182 
183 
184 /***********************************************************************//**
185  * @brief Return size of event bin
186  *
187  * @return Size of event bin (sr MeV s)
188  *
189  * The size of the event bin (units: sr MeV s) is given by
190  * \f[size = \Omega \times \Delta E \times \Delta T\f]
191  * where
192  * \f$\Omega\f$ is the size of the spatial bin in sr,
193  * \f$\Delta E\f$ is the size of the energy bin in MeV, and
194  * \f$\Delta T\f$ is the ontime of the observation in seconds.
195  ***************************************************************************/
196 double GCOMEventBin::size(void) const
197 {
198  // Compute bin size
199  double size = solidangle() * ewidth().MeV() * ontime();
200 
201  // Return bin size
202  return size;
203 }
204 
205 
206 /***********************************************************************//**
207  * @brief Return instrument direction of event bin
208  *
209  * @return Instrument direction of event bin
210  *
211  * @exception GException::invalid_value
212  * Invalid instrument direction pointer.
213  *
214  * Returns reference to the instrument direction of the event bin.
215  ***************************************************************************/
216 const GCOMInstDir& GCOMEventBin::dir(void) const
217 {
218  // Throw an exception if instrument direction pointer is not valid
219  if (m_dir == NULL) {
220  std::string msg = "No valid instrument direction found in event bin";
222  }
223 
224  // Return instrument direction
225  return *m_dir;
226 }
227 
228 
229 /***********************************************************************//**
230  * @brief Return energy of event bin
231  *
232  * @return Energy of event bin
233  *
234  * @exception GException::invalid_value
235  * Invalid energy pointer.
236  *
237  * Returns reference to the energy of the event bin.
238  ***************************************************************************/
239 const GEnergy& GCOMEventBin::energy(void) const
240 {
241  // Throw an exception if energy pointer is not valid
242  if (m_energy == NULL) {
243  std::string msg = "No valid energy found in event bin";
245  }
246 
247  // Return energy
248  return *m_energy;
249 }
250 
251 
252 /***********************************************************************//**
253  * @brief Return time of event bin
254  *
255  * @return Time of event bin
256  *
257  * @exception GException::invalid_value
258  * Invalid time pointer.
259  *
260  * Returns reference to the time of the event bin.
261  ***************************************************************************/
262 const GTime& GCOMEventBin::time(void) const
263 {
264  // Throw an exception if time pointer is not valid
265  if (m_energy == NULL) {
266  std::string msg = "No valid time found in event bin";
268  }
269 
270  // Return time
271  return *m_time;
272 }
273 
274 
275 /***********************************************************************//**
276  * @brief Return number of counts in event bin
277  *
278  * @return Number of counts in event bin
279  *
280  * @exception GCTAException::invalid_value
281  * Invalid counts pointer.
282  *
283  * Returns reference to the number of counts in the event bin.
284  ***************************************************************************/
285 double GCOMEventBin::counts(void) const
286 {
287  // Throw an exception if counts pointer is not valid
288  if (m_counts == NULL) {
289  std::string msg = "No valid counts found in event bin";
291  }
292 
293  // Return counts
294  return *m_counts;
295 }
296 
297 
298 /***********************************************************************//**
299  * @brief Return error in number of counts
300  *
301  * @return Error in number of counts in event bin.
302  *
303  * Returns \f$\sqrt(counts+delta)\f$ as the uncertainty in the number of
304  * counts in the bin. Adding delta avoids uncertainties of 0 which will
305  * lead in the optimisation step to the exlusion of the corresponding bin.
306  * In the actual implementation delta=1e-50.
307  *
308  * @todo The choice of delta has been made somewhat arbitrary, mainly
309  * because the optimizer routines filter error^2 below 1e-100.
310  ***************************************************************************/
311 double GCOMEventBin::error(void) const
312 {
313  // Compute uncertainty
314  double error = std::sqrt(counts()+1.0e-50);
315 
316  // Return error
317  return error;
318 }
319 
320 
321 /***********************************************************************//**
322  * @brief Return solid angle of event bin
323  *
324  * @return Solid angle of event bin
325  *
326  * @exception GException::invalid_value
327  * Invalid solid angle pointer.
328  *
329  * Returns reference to the solid angle of the event bin.
330  ***************************************************************************/
331 const double& GCOMEventBin::solidangle(void) const
332 {
333  // Throw an exception if counts pointer is not valid
334  if (m_solidangle == NULL) {
335  std::string msg = "No valid solid angle found in event bin";
337  }
338 
339  // Return solid angle
340  return *m_solidangle;
341 }
342 
343 
344 /***********************************************************************//**
345  * @brief Return energy width of event bin
346  *
347  * @return Energy width of event bin
348  *
349  * @exception GException::invalid_value
350  * Invalid energy width pointer.
351  *
352  * Returns reference to the energy width of the event bin.
353  ***************************************************************************/
354 const GEnergy& GCOMEventBin::ewidth(void) const
355 {
356  // Throw an exception if energy width pointer is not valid
357  if (m_ewidth == NULL) {
358  std::string msg = "No valid energy width found in event bin";
360  }
361 
362  // Return energy width
363  return *m_ewidth;
364 }
365 
366 
367 /***********************************************************************//**
368  * @brief Return ontime of event bin
369  *
370  * @return Ontime of event bin
371  *
372  * @exception GException::invalid_value
373  * Invalid ontime pointer.
374  *
375  * Returns reference to the ontime of the event bin.
376  ***************************************************************************/
377 const double& GCOMEventBin::ontime(void) const
378 {
379  // Throw an exception if energy width pointer is not valid
380  if (m_ontime == NULL) {
381  std::string msg = "No valid ontime found in event bin";
383  }
384 
385  // Return ontime
386  return *m_ontime;
387 }
388 
389 
390 /***********************************************************************//**
391  * @brief Set instrument direction of event bin
392  *
393  * @param[in] dir Instrument direction of event bin
394  *
395  * @exception GException::invalid_value
396  * No memory available to hold instrument direction.
397  *
398  * Sets the instrument direction of the event bin.
399  ***************************************************************************/
401 {
402  // Throw an exception if no memory has been allocated
403  if (m_dir == NULL) {
404  std::string msg = "No memory available to hold instrument direction.";
406  }
407 
408  // Set instrument direction
409  *m_dir = dir;
410 
411  // Return
412  return;
413 }
414 
415 
416 /***********************************************************************//**
417  * @brief Set energy of event bin
418  *
419  * @param[in] energy Energy of event bin
420  *
421  * @exception GException::invalid_value
422  * No memory available to hold energy.
423  *
424  * Sets the energy of the event bin.
425  ***************************************************************************/
426 void GCOMEventBin::energy(const GEnergy& energy)
427 {
428  // Throw an exception if no memory has been allocated
429  if (m_energy == NULL) {
430  std::string msg = "No memory available to hold energy.";
432  }
433 
434  // Set energy
435  *m_energy = energy;
436 
437  // Return
438  return;
439 }
440 
441 
442 /***********************************************************************//**
443  * @brief Set time of event bin
444  *
445  * @param[in] time Time of event bin
446  *
447  * @exception GException::invalid_value
448  * No memory available to hold instrument direction.
449  *
450  * Sets the time of the event bin.
451  ***************************************************************************/
452 void GCOMEventBin::time(const GTime& time)
453 {
454  // Throw an exception if no memory has been allocated
455  if (m_time == NULL) {
456  std::string msg = "No memory available to hold time.";
458  }
459 
460  // Set time
461  *m_time = time;
462 
463  // Return
464  return;
465 }
466 
467 
468 /***********************************************************************//**
469  * @brief Set number of counts in event bin
470  *
471  * @param[in] counts Number of counts.
472  *
473  * @exception GException::invalid_value
474  * No memory available to hold counts.
475  *
476  * Set the number of counts in the event bin.
477  ***************************************************************************/
478 void GCOMEventBin::counts(const double& counts)
479 {
480  // Throw an exception if counts pointer is not valid
481  if (m_counts == NULL) {
482  std::string msg = "No memory available to hold counts.";
484  }
485 
486  // Set number of counts in event bin
487  *m_counts = counts;
488 
489  // Return
490  return;
491 }
492 
493 
494 /***********************************************************************//**
495  * @brief Set solid angle of event bin
496  *
497  * @param[in] solidangle Solid angle of event bin
498  *
499  * @exception GException::invalid_value
500  * No memory available to hold solid angle.
501  *
502  * Sets the solid angle of the event bin.
503  ***************************************************************************/
504 void GCOMEventBin::solidangle(const double& solidangle)
505 {
506  // Throw an exception if no memory has been allocated
507  if (m_solidangle == NULL) {
508  std::string msg = "No memory available to hold solid angle.";
510  }
511 
512  // Set solid angle
514 
515  // Return
516  return;
517 }
518 
519 
520 /***********************************************************************//**
521  * @brief Set energy width of event bin
522  *
523  * @param[in] ewidth Energy width of event bin
524  *
525  * @exception GException::invalid_value
526  * No memory available to hold energy width.
527  *
528  * Sets the energy width of the event bin.
529  ***************************************************************************/
530 void GCOMEventBin::ewidth(const GEnergy& ewidth)
531 {
532  // Throw an exception if no memory has been allocated
533  if (m_ewidth == NULL) {
534  std::string msg = "No memory available to hold energy width.";
536  }
537 
538  // Set energy width
539  *m_ewidth = ewidth;
540 
541  // Return
542  return;
543 }
544 
545 
546 /***********************************************************************//**
547  * @brief Set ontime of event bin
548  *
549  * @param[in] ontime Ontime of event bin (sec).
550  *
551  * @exception GException::invalid_value
552  * No memory available to hold ontime.
553  *
554  * Sets the ontime of the event bin.
555  ***************************************************************************/
556 void GCOMEventBin::ontime(const double& ontime)
557 {
558  // Throw an exception if no memory has been allocated
559  if (m_ontime == NULL) {
560  std::string msg = "No memory available to hold ontime.";
562  }
563 
564  // Set solid angle
565  *m_ontime = ontime;
566 
567  // Return
568  return;
569 }
570 
571 
572 /***********************************************************************//**
573  * @brief Print event information
574  *
575  * @param[in] chatter Chattiness (defaults to NORMAL).
576  * @return String containing event information.
577  ***************************************************************************/
578 std::string GCOMEventBin::print(const GChatter& chatter) const
579 {
580  // Initialise result string
581  std::string result;
582 
583  // Continue only if chatter is not silent
584  if (chatter != SILENT) {
585 
586  // Append number of counts
587  result.append(gammalib::str(counts()));
588 
589  } // endif: chatter was not silent
590 
591  // Return result
592  return result;
593 }
594 
595 
596 /*==========================================================================
597  = =
598  = Private methods =
599  = =
600  ==========================================================================*/
601 
602 /***********************************************************************//**
603  * @brief Initialise class members
604  *
605  * This method allocates memory for all event bin attributes and intialises
606  * the attributes to well defined initial values.
607  *
608  * The method assumes that on entry no memory is hold by the member pointers.
609  ***************************************************************************/
611 {
612  // Allocate members
613  m_alloc = true;
614  m_index = -1; // Signals that event bin does not correspond to cube
615  m_dir = new GCOMInstDir;
616  m_time = new GTime;
617  m_energy = new GEnergy;
618  m_ewidth = new GEnergy;
619  m_counts = new double;
620  m_solidangle = new double;
621  m_ontime = new double;
622 
623  // Initialise members
624  m_dir->clear();
625  m_time->clear();
626  m_energy->clear();
627  m_ewidth->clear();
628  *m_counts = 0.0;
629  *m_solidangle = 0.0;
630  *m_ontime = 0.0;
631 
632  // Return
633  return;
634 }
635 
636 
637 /***********************************************************************//**
638  * @brief Copy class members
639  *
640  * @param[in] bin Event bin.
641  ***************************************************************************/
643 {
644  // First de-allocate existing memory if needed
645  free_members();
646 
647  // Copy members by cloning
648  m_dir = new GCOMInstDir(*bin.m_dir);
649  m_time = new GTime(*bin.m_time);
650  m_energy = new GEnergy(*bin.m_energy);
651  m_ewidth = new GEnergy(*bin.m_ewidth);
652  m_counts = new double(*bin.m_counts);
653  m_solidangle = new double(*bin.m_solidangle);
654  m_ontime = new double(*bin.m_ontime);
655 
656  // Copy non-pointer members
657  m_index = bin.m_index;
658 
659  // Signal memory allocation
660  m_alloc = true;
661 
662  // Return
663  return;
664 }
665 
666 
667 /***********************************************************************//**
668  * @brief Delete class members
669  *
670  * This method frees all memory of the class attributes and sets the member
671  * pointers to NULL. This method should only be called if new memory is
672  * allocated immediately afterwards (for example by cloning another event
673  * bin), or upon destruction of the object.
674  *
675  * Note that some logic has been implemented that frees only memory that also
676  * has indeed been allocated by the class. Thus if the class only serves as
677  * container to hold memory pointer allocated by someone else (for example
678  * the GCOMEventCube class), no memory is freed.
679  ***************************************************************************/
681 {
682  // If memory was allocated then free members now
683  if (m_alloc) {
684  if (m_dir != NULL) delete m_dir;
685  if (m_time != NULL) delete m_time;
686  if (m_energy != NULL) delete m_energy;
687  if (m_ewidth != NULL) delete m_ewidth;
688  if (m_counts != NULL) delete m_counts;
689  if (m_solidangle != NULL) delete m_solidangle;
690  if (m_ontime != NULL) delete m_ontime;
691  }
692 
693  // Signal member pointers as free
694  m_dir = NULL;
695  m_time = NULL;
696  m_energy = NULL;
697  m_ewidth = NULL;
698  m_counts = NULL;
699  m_solidangle = NULL;
700  m_ontime = NULL;
701 
702  // Signal memory de-allocation
703  m_alloc = false;
704 
705  // Return
706  return;
707 }
const double & ontime(void) const
Return ontime of event bin.
#define G_TIME_SET
#define G_EWIDTH_GET
virtual std::string print(const GChatter &chatter=NORMAL) const
Print event information.
#define G_COUNTS_SET
int m_index
Dataspace index.
#define G_SOLIDANGLE_GET
GCOMInstDir * m_dir
Pointer to bin direction.
#define G_ONTIME_GET
Abstract interface for the event bin class.
Definition: GEventBin.hpp:64
virtual double error(void) const
Return error in number of counts.
#define G_ENERGY_SET
GTime * m_time
Pointer to bin time.
void copy_members(const GCOMEventBin &bin)
Copy class members.
void clear(void)
Clear time.
Definition: GTime.cpp:251
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
#define G_COUNTS_GET
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
virtual const GTime & time(void) const
Return time of event bin.
COMPTEL event bin class.
void free_members(void)
Delete class members.
virtual void clear(void)
Clear instance.
COMPTEL event bin class interface definition.
virtual GCOMEventBin * clone(void) const
Clone instance.
double * m_ontime
Pointer to ontime of bin (seconds)
virtual ~GCOMEventBin(void)
Destructor.
const GEnergy & ewidth(void) const
Return energy width of event bin.
bool m_alloc
Signals proper memory allocation.
#define G_ENERGY_GET
#define G_SOLIDANGLE_SET
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
#define G_ONTIME_SET
GEnergy * m_ewidth
Pointer to energy width of bin.
#define G_EWIDTH_SET
const double & solidangle(void) const
Return solid angle of event bin.
virtual GCOMEventBin & operator=(const GCOMEventBin &bin)
Assignment operator.
GChatter
Definition: GTypemaps.hpp:33
void init_members(void)
Initialise class members.
Definition: GEvent.cpp:141
void free_members(void)
Delete class members.
Definition: GEventBin.cpp:182
#define G_TIME_GET
virtual void clear(void)
Clear instance.
virtual double counts(void) const
Return number of counts in event bin.
void init_members(void)
Initialise class members.
virtual const GCOMInstDir & dir(void) const
Return instrument direction of event bin.
GEnergy * m_energy
Pointer to bin energy.
void free_members(void)
Delete class members.
Definition: GEvent.cpp:163
#define G_DIR_GET
double * m_counts
Pointer to number of counts.
GCOMEventBin(void)
Void constructor.
virtual GEventBin & operator=(const GEventBin &bin)
Assignment operator.
Definition: GEventBin.cpp:105
virtual double size(void) const
Return size of event bin.
void init_members(void)
Initialise class members.
Definition: GEventBin.cpp:160
#define G_DIR_SET
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
Interface for the COMPTEL instrument direction class.
Definition: GCOMInstDir.hpp:41
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
double * m_solidangle
Pointer to solid angle of pixel (sr)
virtual const GEnergy & energy(void) const
Return energy of event bin.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413