ctools  2.0.0
 All Classes Namespaces Files Functions Variables Macros Pages
ctpsfcube.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * ctpsfcube - PSF cube generation tool *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2014-2022 by Chia-Chun Lu *
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 ctpsfcube.cpp
23  * @brief PSF cube generation tool implementation
24  * @author Chia-Chun Lu
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cstdio>
32 #include "ctpsfcube.hpp"
33 #include "GTools.hpp"
34 #include "GWcs.hpp"
35 
36 /* __ Method name definitions ____________________________________________ */
37 
38 /* __ Debug definitions __________________________________________________ */
39 
40 /* __ Coding definitions _________________________________________________ */
41 
42 
43 /*==========================================================================
44  = =
45  = Constructors/destructors =
46  = =
47  ==========================================================================*/
48 
49 /***********************************************************************//**
50  * @brief Void constructor
51  ***************************************************************************/
53 {
54  // Initialise members
55  init_members();
56 
57  // Return
58  return;
59 }
60 
61 
62 /***********************************************************************//**
63  * @brief Observations constructor
64  *
65  * @param[in] obs Observation container.
66  *
67  * This method creates an instance of the class by copying an existing
68  * observations container.
69  ***************************************************************************/
70 ctpsfcube::ctpsfcube(const GObservations& obs) :
71  ctobservation(CTPSFCUBE_NAME, VERSION, obs)
72 {
73  // Initialise members
74  init_members();
75 
76  // Return
77  return;
78 }
79 
80 
81 
82 /***********************************************************************//**
83  * @brief Command line constructor
84  *
85  * @param[in] argc Number of arguments in command line.
86  * @param[in] argv Array of command line arguments.
87  ***************************************************************************/
88 ctpsfcube::ctpsfcube(int argc, char *argv[]) :
89  ctobservation(CTPSFCUBE_NAME, VERSION, argc, argv)
90 {
91  // Initialise members
92  init_members();
93 
94  // Return
95  return;
96 }
97 
98 
99 /***********************************************************************//**
100  * @brief Copy constructor
101  *
102  * @param[in] app Application.
103  ***************************************************************************/
105 {
106  // Initialise members
107  init_members();
108 
109  // Copy members
110  copy_members(app);
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief Destructor
119  ***************************************************************************/
121 {
122  // Free members
123  free_members();
124 
125  // Return
126  return;
127 }
128 
129 
130 /*==========================================================================
131  = =
132  = Operators =
133  = =
134  ==========================================================================*/
135 
136 /***********************************************************************//**
137  * @brief Assignment operator
138  *
139  * @param[in] app Application.
140  * @return Application.
141  ***************************************************************************/
143 {
144  // Execute only if object is not identical
145  if (this != &app) {
146 
147  // Copy base class members
148  this->ctobservation::operator=(app);
149 
150  // Free members
151  free_members();
152 
153  // Initialise members
154  init_members();
155 
156  // Copy members
157  copy_members(app);
158 
159  } // endif: object was not identical
160 
161  // Return this object
162  return *this;
163 }
164 
165 
166 /*==========================================================================
167  = =
168  = Public methods =
169  = =
170  ==========================================================================*/
171 
172 /***********************************************************************//**
173  * @brief Clear ctpsfcube tool
174  *
175  * Clears ctpsfcube tool.
176  ***************************************************************************/
178 {
179  // Free members
180  free_members();
182  this->ctool::free_members();
183 
184  // Clear base class (needed to conserve tool name and version)
185  this->GApplication::clear();
186 
187  // Initialise members
188  this->ctool::init_members();
190  init_members();
191 
192  // Write header into logger
193  log_header();
194 
195  // Return
196  return;
197 }
198 
199 
200 /***********************************************************************//**
201  * @brief Generate the model map(s)
202  *
203  * This method reads the task parameters from the parfile, sets up the
204  * observation container, loops over all CTA observations in the container
205  * and generates a PSF cube from the CTA observations.
206  ***************************************************************************/
208 {
209  // Get task parameters
210  get_parameters();
211 
212  // Write input observation container into logger
213  log_observations(NORMAL, m_obs, "Input observation");
214 
215  // Initialise PSF cube
216  init_cube();
217 
218  // Write header into logger
219  log_header1(TERSE, "Generate PSF cube");
220 
221  // Set pointer to logger dependent on chattiness
222  GLog* logger = (logNormal()) ? &log : NULL;
223 
224  // Fill PSF cube
225  m_psfcube.fill(m_obs, logger);
226 
227  // Write PSF cube into logger
228  log_string(EXPLICIT, m_psfcube.print(m_chatter));
229 
230  // Return
231  return;
232 }
233 
234 
235 /***********************************************************************//**
236  * @brief Save PSF cube
237  *
238  * Saves the PSF cube into a FITS file. A file is only created if the
239  * "outcube" parameter is not empty and if a PSF cube has been computed.
240  ***************************************************************************/
241 void ctpsfcube::save(void)
242 {
243  // Write header into logger
244  log_header1(TERSE, "Save PSF cube");
245 
246  // Save PSF cube if filename is valid and the PSF cube is not empty
247  if ((*this)["outcube"].is_valid() && !m_psfcube.cube().is_empty()) {
248 
249  // Get output filename
250  m_outcube = (*this)["outcube"].filename();
251 
252  // Save PSF cube
253  m_psfcube.save(m_outcube, clobber());
254 
255  // Stamp PSF cube
256  stamp(m_outcube);
257 
258  }
259 
260  // Write into logger what has been done
261  std::string fname = (m_outcube.is_empty()) ? "NONE" : m_outcube.url();
262  if (m_psfcube.cube().is_empty()) {
263  fname.append(" (cube is empty, no file created)");
264  }
265  log_value(NORMAL, "PSF cube file", fname);
266 
267  // Return
268  return;
269 }
270 
271 
272 /*==========================================================================
273  = =
274  = Private methods =
275  = =
276  ==========================================================================*/
277 
278 /***********************************************************************//**
279  * @brief Initialise class members
280  ***************************************************************************/
282 {
283  // Initialise user parameters
284  m_outcube.clear();
285  m_addbounds = true;
286  m_chatter = static_cast<GChatter>(2);
287 
288  // Initialise protected members
289  m_psfcube.clear();
290 
291  // Return
292  return;
293 }
294 
295 
296 /***********************************************************************//**
297  * @brief Copy class members
298  *
299  * @param[in] app Application.
300  ***************************************************************************/
302 {
303  // Copy user parameters
304  m_outcube = app.m_outcube;
305  m_addbounds = app.m_addbounds;
306  m_chatter = app.m_chatter;
307 
308  // Copy protected members
309  m_psfcube = app.m_psfcube;
310 
311  // Return
312  return;
313 }
314 
315 
316 /***********************************************************************//**
317  * @brief Delete class members
318  ***************************************************************************/
320 {
321  // Return
322  return;
323 }
324 
325 
326 /***********************************************************************//**
327  * @brief Get application parameters
328  *
329  * Get all task parameters from parameter file or (if required) by querying
330  * the user. The parameters are read in the correct order.
331  ***************************************************************************/
333 {
334  // Setup observations from "inobs" parameter. Do not accept counts cubes.
335  setup_observations(m_obs, true, true, false);
336 
337  // Get additional binning parameters
338  double amax = (*this)["amax"].real();
339  int anumbins = (*this)["anumbins"].integer();
340 
341  // If the "incube" file name is valid then setup the PSF cube from the
342  // counts cube. Otherwise create a counts cube from the user parameters
343  GCTAEventCube cube = (*this)["incube"].is_valid()
344  ? GCTAEventCube((*this)["incube"].filename())
345  : create_cube(m_obs);
346 
347  // Define PSF cube
348  m_psfcube = GCTACubePsf(cube, amax, anumbins);
349 
350  // Get remaining parameters
351  m_addbounds = (*this)["addbounds"].boolean();
352  m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
353 
354  // If needed later, query output filename now
355  if (read_ahead()) {
356  (*this)["outcube"].query();
357  }
358 
359  // Write parameters into logger
360  log_parameters(TERSE);
361 
362  // Return
363  return;
364 }
365 
366 
367 /***********************************************************************//**
368  * @brief Initialise PSF cube
369  *
370  * Initialise the PSF cube.
371  ***************************************************************************/
373 {
374  // Write header into logger
375  log_header1(TERSE, "Initialise PSF cube");
376 
377  // Extract PSF cube definition
378  const GWcs* proj = static_cast<const GWcs*>(m_psfcube.cube().projection());
379  std::string wcs = m_psfcube.cube().projection()->code();
380  std::string coords = m_psfcube.cube().projection()->coordsys();
381  double x = proj->crval(0);
382  double y = proj->crval(1);
383  double dx = proj->cdelt(0);
384  double dy = proj->cdelt(1);
385  int nx = m_psfcube.cube().nx();
386  int ny = m_psfcube.cube().ny();
387  double dmax = (*this)["amax"].real();
388  int ndbins = (*this)["anumbins"].integer();
389 
390  // Extract energies
391  GEnergies energies = m_psfcube.energies();
392 
393  // If requested, insert energies at all event list energy boundaries
394  if (m_addbounds) {
395 
396  // Loop over all observations
397  for (int i = 0; i < m_obs.size(); ++i) {
398 
399  // Get observation and continue only if it is a CTA observation
400  const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>
401  (m_obs[i]);
402 
403  // Skip observation if it's not a CTA observation
404  if (cta == NULL) {
405  continue;
406  }
407 
408  // Skip observation if it does not contain an event list
409  if (cta->eventtype() != "EventList") {
410  continue;
411  }
412 
413  // Insert energy boundaries
414  energies = insert_energy_boundaries(energies, *cta);
415 
416  } // endfor: looped over all observations
417 
418  } // endif: energy bin insertion requested
419 
420  // Setup PSF cube
421  m_psfcube = GCTACubePsf(wcs, coords, x, y, dx, dy, nx, ny, energies,
422  dmax, ndbins);
423 
424  // Write PSF cube in logger
425  log_string(NORMAL, m_psfcube.print(m_chatter));
426 
427  // Return
428  return;
429 }
ctpsfcube(void)
Void constructor.
Definition: ctpsfcube.cpp:52
void setup_observations(GObservations &obs, const bool &response=true, const bool &list=true, const bool &cube=true)
Setup observation container.
Definition: ctool.cpp:431
GChatter m_chatter
Chattiness.
Definition: ctpsfcube.hpp:74
void init_members(void)
Initialise class members.
Definition: ctool.cpp:321
void process(void)
Generate the model map(s)
Definition: ctpsfcube.cpp:207
void get_parameters(void)
Get application parameters.
Definition: ctpsfcube.cpp:332
void free_members(void)
Delete class members.
Definition: ctpsfcube.cpp:319
void save(void)
Save PSF cube.
Definition: ctpsfcube.cpp:241
void init_members(void)
Initialise class members.
Definition: ctpsfcube.cpp:281
virtual ~ctpsfcube(void)
Destructor.
Definition: ctpsfcube.cpp:120
void log_parameters(const GChatter &chatter)
Log application parameters.
Definition: ctool.cpp:1208
Base class for observation tools.
void free_members(void)
Delete class members.
PSF cube generation tool.
Definition: ctpsfcube.hpp:44
void free_members(void)
Delete class members.
Definition: ctool.cpp:357
void init_members(void)
Initialise class members.
GCTACubePsf m_psfcube
PSF cube.
Definition: ctpsfcube.hpp:77
PSF cube generation tool definition.
void clear(void)
Clear ctpsfcube tool.
Definition: ctpsfcube.cpp:177
GCTAEventCube create_cube(const GObservations &obs)
Create a CTA event cube from user parameters.
Definition: ctool.cpp:1006
ctobservation & operator=(const ctobservation &app)
Assignment operator.
void copy_members(const ctpsfcube &app)
Copy class members.
Definition: ctpsfcube.cpp:301
void init_cube(void)
Initialise PSF cube.
Definition: ctpsfcube.cpp:372
GFilename m_outcube
Output PSF cube file name.
Definition: ctpsfcube.hpp:72
const bool & read_ahead(void) const
Signal whether parameters should be read ahead.
Definition: ctool.hpp:177
bool m_addbounds
Add energies at boundaries?
Definition: ctpsfcube.hpp:73
#define CTPSFCUBE_NAME
Definition: ctpsfcube.hpp:36
ctpsfcube & operator=(const ctpsfcube &app)
Assignment operator.
Definition: ctpsfcube.cpp:142
void log_observations(const GChatter &chatter, const GObservations &obs, const std::string &what="Observation")
Log observation container.
Definition: ctool.cpp:1251
GObservations m_obs
Observation container.
GEnergies insert_energy_boundaries(const GEnergies &energies, const GCTAObservation &obs)
Insert observation energy boundaries into list of energies.
Definition: ctool.cpp:2119