GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAEventBin.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAEventBin.cpp - CTA event bin class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2016 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 GCTAEventBin.cpp
23  * @brief CTA 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 "GCTAEventBin.hpp"
34 #include "GCTAException.hpp"
35 #include "GTools.hpp"
36 
37 /* __ Method name definitions ____________________________________________ */
38 #define G_DIR_GET "GCTAEventBin::dir()"
39 #define G_ENERGY "GCTAEventBin::energy()"
40 #define G_TIME "GCTAEventBin::time()"
41 #define G_COUNTS_GET "GCTAEventBin::counts()"
42 #define G_SOLIDANGLE "GCTAEventBin::solidangle()"
43 #define G_EWIDTH "GCTAEventBin::ewidth()"
44 #define G_ONTIME "GCTAEventBin::ontime()"
45 #define G_WEIGHT "GCTAEventBin::weight()"
46 #define G_DIR_SET "GCTAEventBin::dir(GCTAInstDir&)"
47 #define G_ENERGY_SET "GCTAEventBin::energy(GEnergy&)"
48 #define G_TIME_SET "GCTAEventBin::time(GTime&)"
49 #define G_COUNTS_SET "GCTAEventBin::counts(double&)"
50 #define G_SOLIDANGLE_SET "GCTAEventBin::solidangle(double&)"
51 #define G_EWIDTH_SET "GCTAEventBin::ewidth(GEnergy&)"
52 #define G_ONTIME_SET "GCTAEventBin::ontime(double&)"
53 #define G_WEIGHT_SET "GCTAEventBin::weight(double&)"
54 
55 /* __ Macros _____________________________________________________________ */
56 
57 /* __ Coding definitions _________________________________________________ */
58 
59 /* __ Debug definitions __________________________________________________ */
60 
61 
62 /*==========================================================================
63  = =
64  = Constructors/destructors =
65  = =
66  ==========================================================================*/
67 
68 /***********************************************************************//**
69  * @brief Void constructor
70  ***************************************************************************/
72 {
73  // Initialise class members for clean destruction
74  init_members();
75 
76  // Return
77  return;
78 }
79 
80 
81 /***********************************************************************//**
82  * @brief Copy constructor
83  *
84  * @param[in] bin Event bin.
85  ***************************************************************************/
87 {
88  // Initialise class members for clean destruction
89  init_members();
90 
91  // Copy members
92  copy_members(bin);
93 
94  // Return
95  return;
96 }
97 
98 
99 /***********************************************************************//**
100  * @brief Destructor
101  ***************************************************************************/
103 {
104  // Free members
105  free_members();
106 
107  // Return
108  return;
109 }
110 
111 
112 /*==========================================================================
113  = =
114  = Operators =
115  = =
116  ==========================================================================*/
117 
118 /***********************************************************************//**
119  * @brief Assignment operator
120  *
121  * @param[in] bin Event bin.
122  * @return Event bin.
123  ***************************************************************************/
125 {
126  // Execute only if object is not identical
127  if (this != &bin) {
128 
129  // Copy base class members
130  this->GEventBin::operator=(bin);
131 
132  // Free members
133  free_members();
134 
135  // Initialise private members for clean destruction
136  init_members();
137 
138  // Copy members
139  copy_members(bin);
140 
141  } // endif: object was not identical
142 
143  // Return this object
144  return *this;
145 }
146 
147 
148 /*==========================================================================
149  = =
150  = Public methods =
151  = =
152  ==========================================================================*/
153 
154 /***********************************************************************//**
155  * @brief Clear eventbin
156  ***************************************************************************/
158 {
159  // Free class members (base and derived classes, derived class first)
160  free_members();
161  this->GEventBin::free_members();
162  this->GEvent::free_members();
163 
164  // Initialise members
165  this->GEvent::init_members();
166  this->GEventBin::init_members();
167  init_members();
168 
169  // Return
170  return;
171 }
172 
173 
174 /***********************************************************************//**
175  * @brief Clone event bin
176  *
177  * @return Pointer to deep copy of event bin.
178  ***************************************************************************/
180 {
181  return new GCTAEventBin(*this);
182 }
183 
184 
185 /***********************************************************************//**
186  * @brief Return size of event bin
187  *
188  * The size of the event bin (units sr MeV s) is given by
189  * \f[size = \Omega \times \Delta E \times \Delta T \times W\f]
190  * where
191  * \f$\Omega\f$ is the size of the spatial bin in sr,
192  * \f$\Delta E\f$ is the size of the energy bin in MeV,
193  * \f$\Delta T\f$ is the ontime of the observation in seconds, and
194  * \f$W\f$ is the weight of the bin.
195  ***************************************************************************/
196 double GCTAEventBin::size(void) const
197 {
198  // Compute bin size
199  double size = solidangle() * ewidth().MeV() * ontime() * weight();
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 GCTAInstDir& GCTAEventBin::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& GCTAEventBin::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& GCTAEventBin::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";
267  throw GException::invalid_value(G_TIME, msg);
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 GCTAEventBin::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 GCTAEventBin::error(void) const
312 {
313  // Compute uncertainty
314  double error = 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& GCTAEventBin::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& GCTAEventBin::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& GCTAEventBin::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 Return weight of event bin
392  *
393  * @return Weight of event bin
394  *
395  * @exception GException::invalid_value
396  * Invalid weight pointer.
397  *
398  * Returns reference to the weight of the event bin.
399  ***************************************************************************/
400 const double& GCTAEventBin::weight(void) const
401 {
402  // Throw an exception if weight pointer is not valid
403  if (m_weight == NULL) {
404  std::string msg = "No valid weight found in event bin";
406  }
407 
408  // Return weight
409  return *m_weight;
410 }
411 
412 
413 /***********************************************************************//**
414  * @brief Set instrument direction of event bin
415  *
416  * @param[in] dir Instrument direction of event bin
417  *
418  * @exception GException::invalid_value
419  * No memory available to hold instrument direction.
420  *
421  * Sets the instrument direction of the event bin.
422  ***************************************************************************/
424 {
425  // Throw an exception if no memory has been allocated
426  if (m_dir == NULL) {
427  std::string msg = "No memory available to hold instrument direction.";
429  }
430 
431  // Set instrument direction
432  *m_dir = dir;
433 
434  // Return
435  return;
436 }
437 
438 
439 /***********************************************************************//**
440  * @brief Set energy of event bin
441  *
442  * @param[in] energy Energy of event bin
443  *
444  * @exception GException::invalid_value
445  * No memory available to hold energy.
446  *
447  * Sets the energy of the event bin.
448  ***************************************************************************/
449 void GCTAEventBin::energy(const GEnergy& energy)
450 {
451  // Throw an exception if no memory has been allocated
452  if (m_energy == NULL) {
453  std::string msg = "No memory available to hold energy.";
455  }
456 
457  // Set energy
458  *m_energy = energy;
459 
460  // Return
461  return;
462 }
463 
464 
465 /***********************************************************************//**
466  * @brief Set time of event bin
467  *
468  * @param[in] time Time of event bin
469  *
470  * @exception GException::invalid_value
471  * No memory available to hold instrument direction.
472  *
473  * Sets the time of the event bin.
474  ***************************************************************************/
475 void GCTAEventBin::time(const GTime& time)
476 {
477  // Throw an exception if no memory has been allocated
478  if (m_time == NULL) {
479  std::string msg = "No memory available to hold time.";
481  }
482 
483  // Set time
484  *m_time = time;
485 
486  // Return
487  return;
488 }
489 
490 
491 /***********************************************************************//**
492  * @brief Set number of counts in event bin
493  *
494  * @param[in] counts Number of counts.
495  *
496  * @exception GException::invalid_value
497  * No memory available to hold counts.
498  *
499  * Set the number of counts in the event bin.
500  ***************************************************************************/
501 void GCTAEventBin::counts(const double& counts)
502 {
503  // Throw an exception if counts pointer is not valid
504  if (m_counts == NULL) {
505  std::string msg = "No memory available to hold counts.";
507  }
508 
509  // Set number of counts in event bin
510  *m_counts = counts;
511 
512  // Return
513  return;
514 }
515 
516 
517 /***********************************************************************//**
518  * @brief Set solid angle of event bin
519  *
520  * @param[in] solidangle Solid angle of event bin
521  *
522  * @exception GException::invalid_value
523  * No memory available to hold solid angle.
524  *
525  * Sets the solid angle of the event bin.
526  ***************************************************************************/
527 void GCTAEventBin::solidangle(const double& solidangle)
528 {
529  // Throw an exception if no memory has been allocated
530  if (m_solidangle == NULL) {
531  std::string msg = "No memory available to hold solid angle.";
533  }
534 
535  // Set solid angle
537 
538  // Return
539  return;
540 }
541 
542 
543 /***********************************************************************//**
544  * @brief Set energy width of event bin
545  *
546  * @param[in] ewidth Energy width of event bin
547  *
548  * @exception GException::invalid_value
549  * No memory available to hold energy width.
550  *
551  * Sets the energy width of the event bin.
552  ***************************************************************************/
553 void GCTAEventBin::ewidth(const GEnergy& ewidth)
554 {
555  // Throw an exception if no memory has been allocated
556  if (m_ewidth == NULL) {
557  std::string msg = "No memory available to hold energy width.";
559  }
560 
561  // Set energy width
562  *m_ewidth = ewidth;
563 
564  // Return
565  return;
566 }
567 
568 
569 /***********************************************************************//**
570  * @brief Set ontime of event bin
571  *
572  * @param[in] ontime Ontime of event bin (sec).
573  *
574  * @exception GException::invalid_value
575  * No memory available to hold ontime.
576  *
577  * Sets the ontime of the event bin.
578  ***************************************************************************/
579 void GCTAEventBin::ontime(const double& ontime)
580 {
581  // Throw an exception if no memory has been allocated
582  if (m_ontime == NULL) {
583  std::string msg = "No memory available to hold ontime.";
585  }
586 
587  // Set solid angle
588  *m_ontime = ontime;
589 
590  // Return
591  return;
592 }
593 
594 
595 /***********************************************************************//**
596  * @brief Set weight of event bin
597  *
598  * @param[in] weight Weight angle of event bin
599  *
600  * @exception GException::invalid_value
601  * No memory available to hold weight.
602  *
603  * Sets the weight of the event bin.
604  ***************************************************************************/
605 void GCTAEventBin::weight(const double& weight)
606 {
607  // Throw an exception if no memory has been allocated
608  if (m_weight == NULL) {
609  std::string msg = "No memory available to hold weight.";
611  }
612 
613  // Set weight
614  *m_weight = weight;
615 
616  // Return
617  return;
618 }
619 
620 
621 /***********************************************************************//**
622  * @brief Print event information
623  *
624  * @param[in] chatter Chattiness (defaults to NORMAL).
625  * @return String containing event information.
626  ***************************************************************************/
627 std::string GCTAEventBin::print(const GChatter& chatter) const
628 {
629  // Initialise result string
630  std::string result;
631 
632  // Continue only if chatter is not silent
633  if (chatter != SILENT) {
634 
635  // Append number of counts
636  result.append(gammalib::str(counts()));
637 
638  } // endif: chatter was not silent
639 
640  // Return result
641  return result;
642 }
643 
644 
645 /*==========================================================================
646  = =
647  = Private methods =
648  = =
649  ==========================================================================*/
650 
651 /***********************************************************************//**
652  * @brief Initialise class members
653  *
654  * This method allocates memory for all event bin attributes and intialises
655  * the attributes to well defined initial values.
656  *
657  * The method assumes that on entry no memory is hold by the member pointers.
658  ***************************************************************************/
660 {
661  // Initialise members
662  m_alloc = true;
663  m_ipix = -1; //!< Not part of an event cube
664  m_ieng = -1; //!< Not part of an event cube
665  m_dir = new GCTAInstDir;
666  m_time = new GTime;
667  m_energy = new GEnergy;
668  m_ewidth = new GEnergy;
669  m_counts = new double;
670  m_solidangle = new double;
671  m_ontime = new double;
672  m_weight = new double;
673 
674  // Initialise members
675  m_dir->clear();
676  m_time->clear();
677  m_energy->clear();
678  m_ewidth->clear();
679  *m_counts = 0.0;
680  *m_solidangle = 0.0;
681  *m_ontime = 0.0;
682  *m_weight = 0.0;
683 
684  // Return
685  return;
686 }
687 
688 
689 /***********************************************************************//**
690  * @brief Copy class members
691  *
692  * @param[in] bin Event bin.
693  ***************************************************************************/
695 {
696  // First de-allocate existing memory if needed
697  free_members();
698 
699  // Copy members by cloning
700  m_dir = new GCTAInstDir(*bin.m_dir);
701  m_time = new GTime(*bin.m_time);
702  m_energy = new GEnergy(*bin.m_energy);
703  m_ewidth = new GEnergy(*bin.m_ewidth);
704  m_counts = new double(*bin.m_counts);
705  m_solidangle = new double(*bin.m_solidangle);
706  m_ontime = new double(*bin.m_ontime);
707  m_weight = new double(*bin.m_weight);
708 
709  // Copy non-pointer members
710  m_ipix = bin.m_ipix;
711  m_ieng = bin.m_ieng;
712 
713  // Signal memory allocation
714  m_alloc = true;
715 
716  // Return
717  return;
718 }
719 
720 
721 /***********************************************************************//**
722  * @brief Delete class members
723  *
724  * This method frees all memory of the class attributes and sets the member
725  * pointers to NULL. This method should only be called if new memory is
726  * allocated immediately afterwards (for example by cloning another event
727  * bin), or upon destruction of the object.
728  *
729  * Note that some logic has been implemented that frees only memory that also
730  * has indeed been allocated by the class. Thus if the class only serves as
731  * container to hold memory pointer allocated by someone else (for example
732  * the GCTAEventCube class), no memory is freed.
733  ***************************************************************************/
735 {
736  // If memory was allocated then free members now
737  if (m_alloc) {
738  if (m_dir != NULL) delete m_dir;
739  if (m_time != NULL) delete m_time;
740  if (m_energy != NULL) delete m_energy;
741  if (m_ewidth != NULL) delete m_ewidth;
742  if (m_counts != NULL) delete m_counts;
743  if (m_solidangle != NULL) delete m_solidangle;
744  if (m_ontime != NULL) delete m_ontime;
745  if (m_weight != NULL) delete m_weight;
746  }
747 
748  // Signal member pointers as free
749  m_dir = NULL;
750  m_time = NULL;
751  m_energy = NULL;
752  m_ewidth = NULL;
753  m_counts = NULL;
754  m_solidangle = NULL;
755  m_ontime = NULL;
756  m_weight = NULL;
757 
758  // Signal memory de-allocation
759  m_alloc = false;
760 
761  // Return
762  return;
763 }
double * m_ontime
Pointer to ontime of bin (seconds)
virtual ~GCTAEventBin(void)
Destructor.
double * m_counts
Pointer to number of counts.
virtual const GCTAInstDir & dir(void) const
Return instrument direction of event bin.
const double & weight(void) const
Return weight of event bin.
#define G_COUNTS_SET
GCTAEventBin class interface definition.
Abstract interface for the event bin class.
Definition: GEventBin.hpp:64
int m_ipix
Index in spatial map.
CTA event bin class interface definition.
#define G_TIME_SET
GCTAEventBin(void)
Void constructor.
void clear(void)
Clear time.
Definition: GTime.cpp:251
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
virtual const GEnergy & energy(void) const
Return energy of event bin.
int m_ieng
Index of energy layer.
#define G_SOLIDANGLE
GCTAInstDir * m_dir
Pointer to bin direction.
#define G_ENERGY_SET
#define G_EWIDTH
#define G_DIR_GET
virtual double size(void) const
Return size of event bin.
#define G_ONTIME_SET
virtual GCTAEventBin * clone(void) const
Clone event bin.
GTime * m_time
Pointer to bin time.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
#define G_ENERGY
const GEnergy & ewidth(void) const
Return energy width of event bin.
#define G_SOLIDANGLE_SET
#define G_COUNTS_GET
const double & ontime(void) const
Return ontime of event bin.
virtual const GTime & time(void) const
Return time of event bin.
#define G_WEIGHT_SET
GChatter
Definition: GTypemaps.hpp:33
void free_members(void)
Delete class members.
void init_members(void)
Initialise class members.
Definition: GEvent.cpp:141
void init_members(void)
Initialise class members.
CTA exception handler interface definition.
bool m_alloc
Signals proper memory allocation.
void free_members(void)
Delete class members.
Definition: GEventBin.cpp:182
virtual void clear(void)
Clear eventbin.
double * m_weight
Pointer to weight of bin.
#define G_ONTIME
virtual void clear(void)
Clear CTA instrument direction.
#define G_DIR_SET
virtual GCTAEventBin & operator=(const GCTAEventBin &bin)
Assignment operator.
void free_members(void)
Delete class members.
Definition: GEvent.cpp:163
double * m_solidangle
Pointer to solid angle of pixel (sr)
GEnergy * m_ewidth
Pointer to energy width of bin.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
virtual GEventBin & operator=(const GEventBin &bin)
Assignment operator.
Definition: GEventBin.cpp:105
GEnergy * m_energy
Pointer to bin energy.
const double & solidangle(void) const
Return solid angle of event bin.
#define G_EWIDTH_SET
void copy_members(const GCTAEventBin &bin)
Copy class members.
#define G_WEIGHT
void init_members(void)
Initialise class members.
Definition: GEventBin.cpp:160
#define G_TIME
virtual double counts(void) const
Return number of counts in event bin.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
virtual double error(void) const
Return error in number of counts.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
virtual std::string print(const GChatter &chatter=NORMAL) const
Print event information.