GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTABackgroundPerfTable.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTABackgroundPerfTable.cpp - CTA performance table background class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2014-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 GCTABackgroundPerfTable.cpp
23  * @brief CTA performance table background class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cstdio> // std::fopen, std::fgets, and std::fclose
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GIntegral.hpp"
35 #include "GRan.hpp"
36 #include "GCTAInstDir.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_LOAD "GCTABackgroundPerfTable::load(GFilename&)"
41 #define G_MC "GCTABackgroundPerfTable::mc(GEnergy&, GTime&, GRan&)"
42 
43 /* __ Macros _____________________________________________________________ */
44 
45 /* __ Coding definitions _________________________________________________ */
46 
47 /* __ Debug definitions __________________________________________________ */
48 //#define G_DEBUG_MC //!< Debug Monte Carlo method
49 #define G_LOG_INTERPOLATION //!< Energy interpolate log(background rate)
50 
51 /* __ Constants __________________________________________________________ */
52 
53 
54 /*==========================================================================
55  = =
56  = Constructors/destructors =
57  = =
58  ==========================================================================*/
59 
60 /***********************************************************************//**
61  * @brief Void constructor
62  ***************************************************************************/
64 {
65  // Initialise class members
66  init_members();
67 
68  // Return
69  return;
70 }
71 
72 
73 /***********************************************************************//**
74  * @brief File constructor
75  *
76  * @param[in] filename Performance table file name.
77  *
78  * Construct instance by loading the background information from a
79  * performance table.
80  ***************************************************************************/
83 {
84  // Initialise class members
85  init_members();
86 
87  // Load background from file
88  load(filename);
89 
90  // Return
91  return;
92 }
93 
94 
95 /***********************************************************************//**
96  * @brief Copy constructor
97  *
98  * @param[in] bgd Background.
99  ***************************************************************************/
101  GCTABackground(bgd)
102 {
103  // Initialise class members
104  init_members();
105 
106  // Copy members
107  copy_members(bgd);
108 
109  // Return
110  return;
111 }
112 
113 
114 /***********************************************************************//**
115  * @brief Destructor
116  ***************************************************************************/
118 {
119  // Free members
120  free_members();
121 
122  // Return
123  return;
124 }
125 
126 
127 /*==========================================================================
128  = =
129  = Operators =
130  = =
131  ==========================================================================*/
132 
133 /***********************************************************************//**
134  * @brief Assignment operator
135  *
136  * @param[in] bgd Background.
137  * @return Background.
138  ***************************************************************************/
140 {
141  // Execute only if object is not identical
142  if (this != &bgd) {
143 
144  // Copy base class members
145  this->GCTABackground::operator=(bgd);
146 
147  // Free members
148  free_members();
149 
150  // Initialise private members
151  init_members();
152 
153  // Copy members
154  copy_members(bgd);
155 
156  } // endif: object was not identical
157 
158  // Return this object
159  return *this;
160 }
161 
162 
163 /***********************************************************************//**
164  * @brief Return background rate in units of events MeV\f$^{-1}\f$
165  * s\f$^{-1}\f$ sr\f$^{-1}\f$
166  *
167  * @param[in] logE Log10 of the true photon energy (TeV).
168  * @param[in] detx Tangential coordinate in nominal system (rad).
169  * @param[in] dety Tangential coordinate in nominal system (rad).
170  * @return Background rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
171  *
172  * Returns the background rate for a given energy and detector coordinates.
173  * The rate is given per ontime. The method assures that the background rate
174  * never becomes negative.
175  *
176  * If the performance table contains less than 2 nodes an empty value is
177  * returned.
178  ***************************************************************************/
179 double GCTABackgroundPerfTable::operator()(const double& logE,
180  const double& detx,
181  const double& dety) const
182 {
183  // Initialise rate
184  double rate = this->rate(logE);
185 
186  // Optionally add in Gaussian offset angle dependence
187  if (m_sigma != 0.0) {
188  double theta = std::sqrt(detx*detx + dety*dety) * gammalib::rad2deg;
189  double arg = theta * theta / m_sigma;
190  double scale = std::exp(-0.5 * arg * arg);
191  rate *= scale;
192  }
193 
194  // Return background rate
195  return rate;
196 }
197 
198 
199 /*==========================================================================
200  = =
201  = Public methods =
202  = =
203  ==========================================================================*/
204 
205 /***********************************************************************//**
206  * @brief Clear background
207  *
208  * This method properly resets the background to an initial state.
209  ***************************************************************************/
211 {
212  // Free class members (base and derived classes, derived class first)
213  free_members();
215 
216  // Initialise members
218  init_members();
219 
220  // Return
221  return;
222 }
223 
224 
225 /***********************************************************************//**
226  * @brief Clone background
227  *
228  * @return Pointer to deep copy of background.
229  ***************************************************************************/
231 {
232  return new GCTABackgroundPerfTable(*this);
233 }
234 
235 
236 /***********************************************************************//**
237  * @brief Load background from performance table
238  *
239  * @param[in] filename Performance table file name.
240  *
241  * @exception GException::file_open_error
242  * File could not be opened for read access.
243  *
244  * Loads the background information from a performance table.
245  ***************************************************************************/
247 {
248  // Clear arrays
249  m_energy.clear();
251  m_background.clear();
252  m_log_background.clear();
253 
254  // Allocate line buffer
255  const int n = 1000;
256  char line[n];
257 
258  // Open performance table readonly
259  FILE* fptr = std::fopen(filename.url().c_str(), "r");
260  if (fptr == NULL) {
261  throw GException::file_open_error(G_LOAD, filename.url());
262  }
263 
264  // Read lines
265  while (std::fgets(line, n, fptr) != NULL) {
266 
267  // Split line in elements. Strip empty elements from vector.
268  std::vector<std::string> elements = gammalib::split(line, " ");
269  for (int i = elements.size()-1; i >= 0; i--) {
270  if (gammalib::strip_whitespace(elements[i]).length() == 0) {
271  elements.erase(elements.begin()+i);
272  }
273  }
274 
275  // Skip header
276  if (elements[0].find("log(E)") != std::string::npos) {
277  continue;
278  }
279 
280  // Break loop if end of data table has been reached
281  if (elements[0].find("----------") != std::string::npos) {
282  break;
283  }
284 
285  // Determine on-axis background rate (counts/s/MeV/sr)
286  double logE = gammalib::todouble(elements[0]);
287  double r80 = gammalib::todouble(elements[3]) * gammalib::deg2rad;
288  double bgrate = gammalib::todouble(elements[5]); // in Hz
289  double emin = std::pow(10.0, logE-0.1) * 1.0e6;
290  double emax = std::pow(10.0, logE+0.1) * 1.0e6;
291  double ewidth = emax - emin; // in MeV
292  double solidangle = gammalib::twopi * (1.0 - std::cos(r80)); // in sr
293  if (solidangle > 0.0) {
294  bgrate /= (solidangle * ewidth); // counts/s/MeV/sr
295  }
296  else {
297  bgrate = 0.0;
298  }
299 
300  // Set energy
301  GEnergy energy;
302  energy.log10TeV(logE);
303 
304  // Push elements in node array and vectors
305  m_energy.append(energy);
306  m_log10_energy.append(logE);
307  m_background.push_back(bgrate);
308  m_log_background.push_back(std::log(bgrate));
309 
310  } // endwhile: looped over lines
311 
312  // Close file
313  std::fclose(fptr);
314 
315  // Store filename
317 
318  // Return
319  return;
320 }
321 
322 
323 /***********************************************************************//**
324  * @brief Returns MC instrument direction
325  *
326  * @param[in] energy Photon energy.
327  * @param[in] time Photon arrival time.
328  * @param[in,out] ran Random number generator.
329  * @return CTA instrument direction.
330  ***************************************************************************/
332  const GTime& time,
333  GRan& ran) const
334 {
335  // Simulate theta angle
336  #if defined(G_DEBUG_MC)
337  int n_samples = 0;
338  #endif
339  double sigma_max = 4.0 * std::sqrt(sigma());
340  double u_max = std::sin(sigma_max * gammalib::deg2rad);
341  double value = 0.0;
342  double u = 1.0;
343  double theta = 0.0;
344  do {
345  theta = ran.uniform() * sigma_max;
346  double arg = theta * theta / sigma();
347  double arg2 = arg * arg;
348  value = std::sin(theta * gammalib::deg2rad) * exp(-0.5 * arg2);
349  u = ran.uniform() * u_max;
350  #if defined(G_DEBUG_MC)
351  n_samples++;
352  #endif
353  } while (u > value);
354  theta *= gammalib::deg2rad;
355  #if defined(G_DEBUG_MC)
356  std::cout << "#=" << n_samples << " ";
357  #endif
358 
359  // Simulate azimuth angle
360  double phi = gammalib::twopi * ran.uniform();
361 
362  // Compute detx and dety
363  double detx(0.0);
364  double dety(0.0);
365  if (theta > 0.0 ) {
366  detx = theta * std::cos(phi);
367  dety = theta * std::sin(phi);
368  }
369 
370  // Set instrument direction (in radians)
371  GCTAInstDir dir;
372  dir.detx(detx);
373  dir.dety(dety);
374 
375  // Return instrument direction
376  return dir;
377 }
378 
379 
380 /***********************************************************************//**
381  * @brief Returns background rate integrated over energy interval in units
382  * of events s\f$^{-1}\f$ sr\f$^{-1}\f$
383  *
384  * @param[in] dir Instrument direction.
385  * @param[in] emin Minimum energy of energy interval.
386  * @param[in] emax Maximum energy of energy interval.
387  * @return Integrated background count rate (events s\f$^{-1}\f$ sr\f$^{-1}\f$).
388  *
389  * Returns the background count rate for a given instrument direction that
390  * is integrated over a specified energy interval. The rate is given per
391  * ontime.
392  *
393  * If the energy interval is not positive, a zero background rate is
394  * returned.
395  ***************************************************************************/
397  const GEnergy& emin,
398  const GEnergy& emax) const
399 {
400  // Initialise rate
401  double rate = 0.0;
402 
403  // Continue only if energy interval is positive
404  if (emax > emin) {
405 
406  // Initialise first and second node
407  double x1 = emin.MeV();
408  double x2 = 0.0;
409  double f1 = this->rate(emin.log10TeV());
410  double f2 = 0.0;
411 
412  // Loop over all nodes
413  for (int i = 0; i < size(); ++i) {
414 
415  // If node energy is below x1 then skip node
416  if (m_energy[i].MeV() <= x1) {
417  continue;
418  }
419 
420  // If node energy is above emax then use emax as energy
421  if (m_energy[i] > emax) {
422  x2 = emax.MeV();
423  f2 = this->rate(emax.log10TeV());
424  }
425 
426  // ... otherwise use node energy
427  else {
428  x2 = m_energy[i].MeV();
429  f2 = m_background[i];
430  }
431 
432  // Compute integral
433  rate += gammalib::plaw_integral(x1, f1, x2, f2);
434 
435  // Set second node as first node
436  x1 = x2;
437  f1 = f2;
438 
439  // If node energy is above emax then break now
440  if (m_energy[i] > emax) {
441  break;
442  }
443 
444  } // endfor: looped over all nodes
445 
446  // If last node energy is below emax then compute last part of
447  // integral up to emax
448  if (x1 < emax.MeV()) {
449  x2 = emax.MeV();
450  f2 = this->rate(emax.log10TeV());
451  rate += gammalib::plaw_integral(x1, f1, x2, f2);
452  }
453 
454  // Optionally add in Gaussian offset angle dependence
455  if (m_sigma != 0.0) {
456  double theta = dir.theta() * gammalib::rad2deg;
457  double arg = theta * theta / m_sigma;
458  double scale = std::exp(-0.5 * arg * arg);
459  rate *= scale;
460  }
461 
462  } // endif: energy interval was positive
463 
464  // Return background rate
465  return rate;
466 }
467 
468 
469 /***********************************************************************//**
470  * @brief Print background information
471  *
472  * @param[in] chatter Chattiness.
473  * @return String containing background information.
474  ***************************************************************************/
475 std::string GCTABackgroundPerfTable::print(const GChatter& chatter) const
476 {
477  // Initialise result string
478  std::string result;
479 
480  // Continue only if chatter is not silent
481  if (chatter != SILENT) {
482 
483  // Compute energy boundaries in TeV
484  double emin = (size() > 0) ? m_energy[0].TeV() : 0.0;
485  double emax = (size() > 1) ? m_energy[size()-1].TeV() : 0.0;
486 
487  // Append header
488  result.append("=== GCTABackgroundPerfTable ===");
489 
490  // Append information
491  result.append("\n"+gammalib::parformat("Filename")+m_filename);
492  result.append("\n"+gammalib::parformat("Number of energy bins") +
493  gammalib::str(size()));
494  result.append("\n"+gammalib::parformat("Energy range"));
495  result.append(gammalib::str(emin) +
496  " - " +
497  gammalib::str(emax) +
498  " TeV");
499 
500  // Append offset angle dependence
501  if (m_sigma == 0) {
502  result.append("\n"+gammalib::parformat("Offset angle dependence") +
503  "none");
504  }
505  else {
506  std::string txt = "Fixed sigma=" + gammalib::str(m_sigma);
507  result.append("\n"+gammalib::parformat("Offset angle dependence") +
508  txt);
509  }
510 
511  } // endif: chatter was not silent
512 
513  // Return result
514  return result;
515 }
516 
517 
518 /*==========================================================================
519  = =
520  = Private methods =
521  = =
522  ==========================================================================*/
523 
524 /***********************************************************************//**
525  * @brief Initialise class members
526  ***************************************************************************/
528 {
529  // Initialise members
530  m_filename.clear();
531  m_energy.clear();
533  m_background.clear();
534  m_log_background.clear();
535  m_sigma = 3.0;
536 
537  // Initialise MC cache
539 
540  // Return
541  return;
542 }
543 
544 
545 /***********************************************************************//**
546  * @brief Copy class members
547  *
548  * @param[in] bgd Background.
549  ***************************************************************************/
551 {
552  // Copy members
553  m_filename = bgd.m_filename;
554  m_energy = bgd.m_energy;
558  m_sigma = bgd.m_sigma;
559 
560  // Copy MC cache
562 
563  // Return
564  return;
565 }
566 
567 
568 /***********************************************************************//**
569  * @brief Delete class members
570  ***************************************************************************/
572 {
573  // Return
574  return;
575 }
576 
577 
578 /***********************************************************************//**
579  * @brief Returns integral over radial model (in steradians)
580  *
581  * Computes
582  * \f[\Omega = 2 \pi \int_0^{\pi} \sin \theta f(\theta) d\theta\f]
583  * where
584  * \f[f(\theta) = \exp \left(-\frac{1}{2}
585  * \left( \frac{\theta^2}{\sigma} \right)^2 \right)\f]
586  * \f$\theta\f$ is the offset angle (in degrees), and
587  * \f$\sigma\f$ is the Gaussian width (in degrees\f$^2\f$).
588  *
589  * The integration is performed numerically, and the upper integration bound
590  * \f$\pi\f$
591  * is set to
592  * \f$\sqrt(10 \sigma)\f$
593  * to reduce inaccuracies in the numerical integration.
594  ***************************************************************************/
596 {
597  // Allocate integrand
601 
602  // Allocate intergal
603  GIntegral integral(&integrand);
604 
605  // Set upper integration boundary
606  double offset_max = std::sqrt(10.0*sigma()) * gammalib::deg2rad;
607  if (offset_max > gammalib::pi) {
608  offset_max = gammalib::pi;
609  }
610 
611  // Perform numerical integration
612  double solidangle = integral.romberg(0.0, offset_max) * gammalib::twopi;
613 
614  // Return integral
615  return solidangle;
616 }
617 
618 
619 /***********************************************************************//**
620  * @brief Initialise Monte Carlo cache
621  *
622  * @todo Verify assumption made about the solid angles of the response table
623  * elements.
624  * @todo Add optional sampling on a finer spatial grid.
625  ***************************************************************************/
627 {
628  // Initialise cache
630 
631  // Compute solid angle of model
632  double solidangle = this->solidangle();
633 
634  // Loop over nodes
635  for (int i = 0; i < size(); ++i) {
636 
637  // Compute total rate
638  double total_rate = m_background[i] * solidangle;
639 
640  // Set node
641  if (total_rate > 0.0) {
642  m_mc_spectrum.append(m_energy[i], total_rate);
643  }
644 
645  }
646 
647  // Return
648  return;
649 }
650 
651 
652 /***********************************************************************//**
653  * @brief Return background rate for a given energy (events/s/MeV/sr)
654  *
655  * @param[in] logE Log10 of the true photon energy (TeV).
656  * @return Background rate (events/s/MeV/sr)
657  ***************************************************************************/
658 double GCTABackgroundPerfTable::rate(const double& logE) const
659 {
660  // Initialise rate
661  double rate = 0.0;
662 
663  // Continue only if there are at least two nodes
664  if (size() > 1) {
665 
666  // Get background rate
667  #if defined(G_LOG_INTERPOLATION)
669  #else
671  #endif
672 
673  // Make sure that background rate is not negative
674  if (rate < 0.0) {
675  rate = 0.0;
676  }
677 
678  } // endif: there were enough rates
679 
680  // Return
681  return rate;
682 }
GCTABackgroundPerfTable(void)
Void constructor.
GCTABackgroundPerfTable & operator=(const GCTABackgroundPerfTable &bgd)
Assignment operator.
std::vector< double > m_log_background
log(background rate)
void clear(void)
Clear energy container.
Definition: GEnergies.cpp:225
void init_members(void)
Initialise class members.
double solidangle(void) const
Returns integral over radial model (in steradians)
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
void detx(const double &x)
Set DETX coordinate (in radians)
Random number generator class definition.
const double pi
Definition: GMath.hpp:35
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
GCTABackground & operator=(const GCTABackground &bgd)
Assignment operator.
virtual ~GCTABackgroundPerfTable(void)
Destructor.
CTA performance table background class definition.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
void free_members(void)
Delete class members.
Random number generator class.
Definition: GRan.hpp:44
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:862
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
void load(const GFilename &filename)
Load background from performance table.
Time class.
Definition: GTime.hpp:54
GFilename filename(void) const
Return filename.
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
void copy_members(const GCTABackgroundPerfTable &bgd)
Copy class members.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
GCTABackgroundPerfTable * clone(void) const
Clone background.
void append(const GEnergy &energy, const double &intensity)
Append node.
GNodeArray m_log10_energy
log10(E) nodes for background interpolation
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
GEnergies m_energy
Vector of energies.
const double & sigma(void) const
Return sigma for offset angle dependence.
double rate(const double &logE) const
Return background rate for a given energy (events/s/MeV/sr)
double theta(void) const
Return offset angle (in radians)
Filename class.
Definition: GFilename.hpp:62
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
void init_mc_cache(void) const
Initialise Monte Carlo cache.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
CTA instrument direction class interface definition.
GChatter
Definition: GTypemaps.hpp:33
int size(void) const
Return number of node energies in response.
virtual double operator()(const double &logE, const double &detx, const double &dety) const
Return background rate in units of events MeV s sr .
std::string print(const GChatter &chatter=NORMAL) const
Print background information.
GModelSpectralNodes m_mc_spectrum
Background spectrum.
void init_members(void)
Initialise class members.
GCTAInstDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC instrument direction.
double plaw_integral(const double &x1, const double &f1, const double &x2, const double &f2)
Returns the integral of a power law.
Definition: GMath.cpp:566
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
virtual void clear(void)
Clear spectral nodes model.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
void clear(void)
Clear background.
GEnergy & append(const GEnergy &energy)
Append energy to container.
Definition: GEnergies.cpp:303
Abstract base class for the CTA background model.
CTA performance table background class.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
double rate_ebin(const GCTAInstDir &dir, const GEnergy &emin, const GEnergy &emax) const
Returns background rate integrated over energy interval in units of events s sr .
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
void dety(const double &y)
Set DETY coordinate (in radians)
const double twopi
Definition: GMath.hpp:36
const double rad2deg
Definition: GMath.hpp:44
std::vector< double > m_background
Background rate.
double m_sigma
Sigma for offset angle computation (0=none)
GFilename m_filename
Name of background response file.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
Integration class interface definition.
Mathematical function definitions.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:805
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
#define G_LOAD
void free_members(void)
Delete class members.