GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMSelection.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMSelection.cpp - COMPTEL selection set class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017 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 GCOMSelection.cpp
23  * @brief COMPTEL selection set class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GCOMSelection.hpp"
32 #include "GCOMEventAtom.hpp"
33 
34 /* __ Method name definitions ____________________________________________ */
35 
36 /* __ Macros _____________________________________________________________ */
37 
38 /* __ Coding definitions _________________________________________________ */
39 
40 /* __ Debug definitions __________________________________________________ */
41 
42 
43 
44 /*==========================================================================
45  = =
46  = Constructors/destructors =
47  = =
48  ==========================================================================*/
49 
50 /***********************************************************************//**
51  * @brief Void constructor
52  ***************************************************************************/
54 {
55  // Initialise class members
56  init_members();
57 
58  // Return
59  return;
60 }
61 
62 
63 /***********************************************************************//**
64  * @brief Copy constructor
65  *
66  * @param[in] select COMPTEL selection set.
67  ***************************************************************************/
69 {
70  // Initialise class members
71  init_members();
72 
73  // Copy members
74  copy_members(select);
75 
76  // Return
77  return;
78 }
79 
80 
81 /***********************************************************************//**
82  * @brief Destructor
83  ***************************************************************************/
85 {
86  // Free members
87  free_members();
88 
89  // Return
90  return;
91 }
92 
93 
94 /*==========================================================================
95  = =
96  = Operators =
97  = =
98  ==========================================================================*/
99 
100 /***********************************************************************//**
101  * @brief Assignment operator
102  *
103  * @param[in] select COMPTEL selection set.
104  * @return COMPTEL selection set.
105  ***************************************************************************/
107 {
108  // Execute only if object is not identical
109  if (this != &select) {
110 
111  // Free members
112  free_members();
113 
114  // Initialise private members
115  init_members();
116 
117  // Copy members
118  copy_members(select);
119 
120  } // endif: object was not identical
121 
122  // Return this object
123  return *this;
124 }
125 
126 
127 /*==========================================================================
128  = =
129  = Public methods =
130  = =
131  ==========================================================================*/
132 
133 /***********************************************************************//**
134  * @brief Clear COMPTEL selection set
135  ***************************************************************************/
137 {
138  // Free members
139  free_members();
140 
141  // Initialise private members
142  init_members();
143 
144  // Return
145  return;
146 }
147 
148 
149 /***********************************************************************//**
150  * @brief Clone COMPTEL selection set
151  *
152  * @return Pointer to deep copy of COMPTEL selection set.
153  ***************************************************************************/
155 {
156  return new GCOMSelection(*this);
157 }
158 
159 
160 /***********************************************************************//**
161  * @brief Initialise selection statistics
162  ***************************************************************************/
164 {
165  // Initialise statistics members
167  m_num_events_used = 0;
169  m_num_e1_min = 0;
170  m_num_e1_max = 0;
171  m_num_e2_min = 0;
172  m_num_e2_max = 0;
173  m_num_tof_min = 0;
174  m_num_tof_max = 0;
175  m_num_psd_min = 0;
176  m_num_psd_max = 0;
177  m_num_zeta_min = 0;
178  m_num_zeta_max = 0;
179  m_num_reflag_min = 0;
180  m_num_reflag_max = 0;
181  m_num_vetoflag_min = 0;
182  m_num_vetoflag_max = 0;
183  m_num_no_scatter = 0;
185 
186  // Return
187  return;
188 }
189 
190 
191 /***********************************************************************//**
192  * @brief Check if event should be used
193  *
194  * @param[in] event Event.
195  * @return True if event should be used, false otherwise.
196  ***************************************************************************/
197 bool GCOMSelection::use_event(const GCOMEventAtom& event) const
198 {
199  // Initialise usage flag
200  bool use = true;
201 
202  // Increment number of checked events
204 
205  // Compute zeta for event
206  double zeta = event.eha() - event.phibar();
207 
208  // Check for bad minitelescopes
209  if (event.modcom() < 1 && event.modcom() > 98) {
211  use = false;
212  }
213 
214  // Check whether the event has a scatter angle determined. This is
215  // signaled by a scatter angle -10e20 in radians.
216  else if (event.theta() < -1.0e3) {
218  use = false;
219  }
220 
221  // Apply event selection
222  else if (event.e1() < m_e1_min) {
223  m_num_e1_min++;
224  use = false;
225  }
226  else if (event.e1() > m_e1_max) {
227  m_num_e1_max++;
228  use = false;
229  }
230  else if (event.e2() < m_e2_min) {
231  m_num_e2_min++;
232  use = false;
233  }
234  else if (event.e2() > m_e2_max) {
235  m_num_e2_max++;
236  use = false;
237  }
238  else if (event.tof() < m_tof_min) {
239  m_num_tof_min++;
240  use = false;
241  }
242  else if (event.tof() > m_tof_max) {
243  m_num_tof_max++;
244  use = false;
245  }
246  else if (event.psd() < m_psd_min) {
247  m_num_psd_min++;
248  use = false;
249  }
250  else if (event.psd() > m_psd_max) {
251  m_num_psd_max++;
252  use = false;
253  }
254  else if (zeta < m_zeta_min) {
255  m_num_zeta_min++;
256  use = false;
257  }
258  else if (zeta > m_zeta_max) {
259  m_num_zeta_max++;
260  use = false;
261  }
262  else if (event.reflag() < m_reflag_min) {
264  use = false;
265  }
266  else if (event.reflag() > m_reflag_max) {
268  use = false;
269  }
270  else if (event.veto() < m_vetoflag_min) {
272  use = false;
273  }
274  else if (event.veto() > m_vetoflag_max) {
276  use = false;
277  }
278 
279  // Update acceptance and rejection statistics
280  if (use) {
282  }
283  else {
285  }
286 
287  // Return usage flag
288  return use;
289 }
290 
291 
292 /***********************************************************************//**
293  * @brief Print COMPTEL selection set
294  *
295  * @param[in] chatter Chattiness.
296  * @return String containing COMPTEL selection set information.
297  ***************************************************************************/
298 std::string GCOMSelection::print(const GChatter& chatter) const
299 {
300  // Initialise result string
301  std::string result;
302 
303  // Continue only if chatter is not silent
304  if (chatter != SILENT) {
305 
306  // Append header
307  result.append("=== GCOMSelection ===");
308 
309  // Append number of checked events
310  result.append("\n"+gammalib::parformat("Checked events"));
311  result.append(gammalib::str(m_num_events_checked));
312  result.append("\n"+gammalib::parformat("Accepted events"));
313  result.append(gammalib::str(m_num_events_used));
314  result.append("\n"+gammalib::parformat("Rejected events"));
315  result.append(gammalib::str(m_num_events_rejected));
316 
317  // Append E1 selection statistics
318  result.append("\n"+gammalib::parformat("E1 selection"));
319  result.append(gammalib::str(m_num_e1_min)+" < [");
320  result.append(gammalib::str(m_e1_min)+" - ");
321  result.append(gammalib::str(m_e1_max)+" MeV] < ");
322  result.append(gammalib::str(m_num_e1_max));
323 
324  // Append E2 selection statistics
325  result.append("\n"+gammalib::parformat("E2 selection"));
326  result.append(gammalib::str(m_num_e2_min)+" < [");
327  result.append(gammalib::str(m_e2_min)+" - ");
328  result.append(gammalib::str(m_e2_max)+" MeV] < ");
329  result.append(gammalib::str(m_num_e2_max));
330 
331  // Append TOF selection statistics
332  result.append("\n"+gammalib::parformat("TOF selection"));
333  result.append(gammalib::str(m_num_tof_min)+" < [");
334  result.append(gammalib::str(m_tof_min)+" - ");
335  result.append(gammalib::str(m_tof_max)+"] < ");
336  result.append(gammalib::str(m_num_tof_max));
337 
338  // Append PSD selection statistics
339  result.append("\n"+gammalib::parformat("PSD selection"));
340  result.append(gammalib::str(m_num_psd_min)+" < [");
341  result.append(gammalib::str(m_psd_min)+" - ");
342  result.append(gammalib::str(m_psd_max)+"] < ");
343  result.append(gammalib::str(m_num_psd_max));
344 
345  // Append zeta angle selection statistics
346  result.append("\n"+gammalib::parformat("Zeta angle selection"));
347  result.append(gammalib::str(m_num_zeta_min)+" < [");
348  result.append(gammalib::str(m_zeta_min)+" - ");
349  result.append(gammalib::str(m_zeta_max)+" deg] < ");
350  result.append(gammalib::str(m_num_zeta_max));
351 
352  // Append rejection flag selection statistics
353  result.append("\n"+gammalib::parformat("Rejection flag selection"));
354  result.append(gammalib::str(m_num_reflag_min)+" < [");
355  result.append(gammalib::str(m_reflag_min)+" - ");
356  result.append(gammalib::str(m_reflag_max)+"] < ");
357  result.append(gammalib::str(m_num_reflag_max));
358 
359  // Append veto flag selection statistics
360  result.append("\n"+gammalib::parformat("Veto flag selection"));
361  result.append(gammalib::str(m_num_vetoflag_min)+" < [");
362  result.append(gammalib::str(m_vetoflag_min)+" - ");
363  result.append(gammalib::str(m_vetoflag_max)+"] < ");
364  result.append(gammalib::str(m_num_vetoflag_max));
365 
366  // Append other statistics
367  result.append("\n"+gammalib::parformat("No scatter angle"));
368  result.append(gammalib::str(m_num_no_scatter));
369  result.append("\n"+gammalib::parformat("Invalid minitelescope"));
370  result.append(gammalib::str(m_num_invalid_modcom));
371 
372  } // endif: chatter was not silent
373 
374  // Return result
375  return result;
376 }
377 
378 
379 /*==========================================================================
380  = =
381  = Private methods =
382  = =
383  ==========================================================================*/
384 
385 /***********************************************************************//**
386  * @brief Initialise class members
387  ***************************************************************************/
389 {
390  // Initialise members
391  m_e1_min = 0.070; //!< Minimum D1 energy deposit (MeV)
392  m_e1_max = 20.0; //!< Maximum D1 energy deposit (MeV)
393  m_e2_min = 0.650; //!< Minimum D2 energy deposit (MeV)
394  m_e2_max = 30.0; //!< Maximum D2 energy deposit (MeV)
395  m_tof_min = 115.0; //!< Minimum TOF window
396  m_tof_max = 130.0; //!< Maximum TOF window
397  m_psd_min = 0.0; //!< Minimum PSD window
398  m_psd_max = 110.0; //!< Maximum PSD window
399  m_zeta_min = 0.0; //!< Minimum Earth horizon angle - Phibar window
400  m_zeta_max = 180.0; //!< Maximum Earth horizon angle - Phibar window
401  m_reflag_min = 1; //!< Minimum rejection flag
402  m_reflag_max = 1000; //!< Maximum rejection flag
403  m_vetoflag_min = 0; //!< Minimum veto flag
404  m_vetoflag_max = 0; //!< Maximum veto flag
405 
406  // Initialise statistics
407  init_statistics();
408 
409  // Return
410  return;
411 }
412 
413 
414 /***********************************************************************//**
415  * @brief Copy class members
416  *
417  * @param[in] select COMPTEL selection set.
418  ***************************************************************************/
420 {
421  // Copy members
422  m_e1_min = select.m_e1_min;
423  m_e1_max = select.m_e1_max;
424  m_e2_min = select.m_e2_min;
425  m_e2_max = select.m_e2_max;
426  m_tof_min = select.m_tof_min;
427  m_tof_max = select.m_tof_max;
428  m_psd_min = select.m_psd_min;
429  m_psd_max = select.m_psd_max;
430  m_zeta_min = select.m_zeta_min;
431  m_zeta_max = select.m_zeta_max;
432  m_reflag_min = select.m_reflag_min;
433  m_reflag_max = select.m_reflag_max;
436 
437  // Copy statistics
441  m_num_e1_min = select.m_num_e1_min;
442  m_num_e1_max = select.m_num_e1_max;
443  m_num_e2_min = select.m_num_e2_min;
444  m_num_e2_max = select.m_num_e2_max;
445  m_num_tof_min = select.m_num_tof_min;
446  m_num_tof_max = select.m_num_tof_max;
447  m_num_psd_min = select.m_num_psd_min;
448  m_num_psd_max = select.m_num_psd_max;
457 
458  // Return
459  return;
460 }
461 
462 
463 /***********************************************************************//**
464  * @brief Delete class members
465  ***************************************************************************/
467 {
468  // Return
469  return;
470 }
void modcom(const int &modcom)
Set mini telescope.
void copy_members(const GCOMSelection &select)
Copy class members.
int m_vetoflag_max
Maximum veto flag.
int m_num_tof_max
Number of events above TOF threshold.
int m_reflag_min
Minimum rejection flag.
COMPTEL event atom class definition.
void reflag(const int &reflag)
Set rejection flag.
void theta(const float &theta)
Set scatter zenith angle.
void tof(const float &tof)
Set TOF value.
double m_psd_min
Minimum PSD window.
double m_e2_max
Maximum D2 energy deposit (MeV)
double m_e2_min
Minimum D2 energy deposit (MeV)
void e2(const float &e2)
Set D2 module energy deposit.
COMPTEL selection set class.
int m_reflag_max
Maximum rejection flag.
double m_psd_max
Maximum PSD window.
int m_vetoflag_min
Minimum veto flag.
COMPTEL selection set class definition.
double m_zeta_max
Maximum Earth horizon angle - Phibar window.
void free_members(void)
Delete class members.
int m_num_zeta_max
Number of events above Zeta threshold.
double m_zeta_min
Minimum Earth horizon angle - Phibar window.
void veto(const int &veto)
Set veto flag.
int m_num_reflag_max
Number of events above rejection flag threshold.
GCOMSelection(void)
Void constructor.
int m_num_e1_min
Number of events below E1 threshold.
void init_members(void)
Initialise class members.
void init_statistics(void) const
Initialise selection statistics.
GCOMSelection & operator=(const GCOMSelection &select)
Assignment operator.
virtual ~GCOMSelection(void)
Destructor.
GChatter
Definition: GTypemaps.hpp:33
virtual std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL selection set.
void psd(const float &psd)
Set PSD value.
double m_e1_max
Maximum D1 energy deposit (MeV)
int m_num_e2_max
Number of events above E2 threshold.
double m_tof_max
Maximum TOF window.
int m_num_events_rejected
Number of rejected events.
virtual GCOMSelection * clone(void) const
Clone COMPTEL selection set.
int m_num_vetoflag_min
Number of events below rejection flag threshold.
int m_num_reflag_min
Number of events below rejection flag threshold.
bool use_event(const GCOMEventAtom &event) const
Check if event should be used.
void e1(const float &e1)
Set D1 module energy deposit.
int m_num_invalid_modcom
Number of events with invalid minitelescopes.
int m_num_events_checked
Number of checked events.
double m_tof_min
Minimum TOF window.
int m_num_tof_min
Number of events below TOF threshold.
virtual void clear(void)
Clear COMPTEL selection set.
int m_num_psd_max
Number of events above PSD threshold.
int m_num_no_scatter
Number of events without scatter angle.
int m_num_psd_min
Number of events below PSD threshold.
double m_e1_min
Minimum D1 energy deposit (MeV)
int m_num_zeta_min
Number of events below Zeta threshold.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
int m_num_e2_min
Number of events below E2 threshold.
int m_num_events_used
Number of used events.
int m_num_vetoflag_max
Number of events above rejection flag threshold.
int m_num_e1_max
Number of events above E1 threshold.
COMPTEL event atom class.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413