43#define G_MC "GCTABackground2D::mc(GEnergy&, GTime&, GRan&)"
44#define G_SET_MEMBERS "GCTABackground2D::set_members()"
45#define G_INIT_MC_CACHE "GCTABackground2D::init_mc_cache()"
196 const double& dety)
const
205 double theta = std::sqrt(detx * detx + dety * dety);
217 double rate_emin = this->
rate(energy_nodes.
inx_left(), theta);
218 double rate_emax = this->
rate(energy_nodes.
inx_right(), theta);
222 if (rate_emin > 0.0 && rate_emax > 0.0) {
225 rate = std::exp(energy_nodes.
wgt_left() * std::log(rate_emin) +
226 energy_nodes.
wgt_right() * std::log(rate_emax));
312 for (
int i = 0; i <
table.ncols(); ++i) {
313 std::string colname =
table[i]->name();
314 if ((colname !=
"THETA_LO") && (colname !=
"THETA_HI") &&
315 (colname !=
"ENERG_LO") && (colname !=
"ENERG_HI") &&
316 (colname !=
"BKG")) {
376 if (extname.empty()) {
434 table.extname(extname);
473 if (energy < emin || energy > emax) {
474 std::string msg =
"Event energy "+energy.
print()+
" is outside the "
475 "energy range ["+emin.
print()+
", "+emax.
print()+
"] "
476 "covered by the background response table. Please "
477 "restrict the energy range of the simulation to the "
478 "validity range of the background response table.";
498 int inx_left = energy_nodes.
inx_left();
499 int inx_right = energy_nodes.
inx_right();
503 double max_rate =
m_mc_max[inx_left];
504 if (
m_mc_max[inx_right] > max_rate) {
522 double value = this->
operator()(logE, detx, dety);
529 if (value > max_rate) {
530 std::string msg =
"Background rate "+
gammalib::str(value)+
" "
531 "for "+energy.
print()+
" at "
535 "is larger than the maximum expected rate "+
537 "wrong with the code.";
542 double uniform = ran.
uniform() * max_rate;
545 if (uniform <= value) {
588 double theta = std::sqrt(dir.
detx() * dir.
detx() + dir.
dety() * dir.
dety());
594 double x1 = emin.
MeV();
616 f2 = this->
rate(i, theta);
620 if (f1 > 0.0 && f2 > 0.0) {
637 if (x1 < emax.
MeV()) {
640 if (f1 > 0.0 && f2 > 0.0) {
669 result.append(
"=== GCTABackground2D ===");
803 std::string msg =
"Expected two-dimensional background "
804 "response table but found "+
806 " dimensions. Please specify a two-dimensional "
870 const int& iebin)
const
896 const int nsubbins = 10;
945 for (
int i = 0; i < energies.
size(); ++i) {
948 double logE = energies[i].log10TeV();
953 double total_flux = 0.0;
954 double dx = 0.5 * detx_bin;
955 double dy = 0.5 * dety_bin;
957 for (
int ix = 0; ix < 2*
m_num_theta; ++ix, detx += detx_bin) {
959 for (
int iy = 0; iy < 2*
m_num_theta; ++iy, dety += dety_bin) {
962 double i0 = (*this)(logE, detx, dety);
963 double i1 = (*this)(logE, detx-dx, dety-dy);
964 double i2 = (*this)(logE, detx, dety-dy);
965 double i3 = (*this)(logE, detx+dx, dety-dy);
966 double i4 = (*this)(logE, detx+dx, dety);
967 double i5 = (*this)(logE, detx+dx, dety+dy);
968 double i6 = (*this)(logE, detx, dety+dy);
969 double i7 = (*this)(logE, detx-dx, dety+dy);
970 double i8 = (*this)(logE, detx-dx, dety);
973 double s1 =
solid_angle(detx, dety, detx-dx, dety-dy, detx, dety-dy);
974 double s2 =
solid_angle(detx, dety, detx, dety-dy, detx+dx, dety-dy);
975 double s3 =
solid_angle(detx, dety, detx+dx, dety-dy, detx+dx, dety);
976 double s4 =
solid_angle(detx, dety, detx+dx, dety, detx+dx, dety+dy);
977 double s5 =
solid_angle(detx, dety, detx+dx, dety+dy, detx, dety+dy);
978 double s6 =
solid_angle(detx, dety, detx, dety+dy, detx-dx, dety+dy);
979 double s7 =
solid_angle(detx, dety, detx-dx, dety+dy, detx-dx, dety);
980 double s8 =
solid_angle(detx, dety, detx-dx, dety, detx-dx, dety-dy);
983 double flux1 = s1 * (i1 + i2 + i0);
984 double flux2 = s2 * (i2 + i3 + i0);
985 double flux3 = s3 * (i3 + i4 + i0);
986 double flux4 = s4 * (i4 + i5 + i0);
987 double flux5 = s5 * (i5 + i6 + i0);
988 double flux6 = s6 * (i6 + i7 + i0);
989 double flux7 = s7 * (i7 + i8 + i0);
990 double flux8 = s8 * (i8 + i1 + i0);
991 double flux = (flux1 + flux2 + flux3 + flux4 +
992 flux5 + flux6 + flux7 + flux8) / 3.0;
1004 if (total_flux > 0.0) {
1009 #if defined(G_DEBUG_MC_INIT)
1010 std::cout <<
"Energy=" << energies[i];
1011 std::cout <<
" Rate_total=" << total_flux;
1012 std::cout <<
" events/(s MeV)" << std::endl;
1020 std::string msg =
"Background response table is empty. Please provide "
1021 "a valid three-dimensional background template.";
1046 double max_rate = 0.0;
1049 for (
int itheta = 0; itheta <
m_num_theta; ++itheta) {
1052 int inx =
index(itheta, iebin);
1059 if (
rate > max_rate) {
1107 const double& detx2,
const double& dety2,
1108 const double& detx3,
const double& dety3)
const
1111 double theta1 = std::sqrt(detx1 * detx1 + dety1 * dety1);
1112 double phi1 = std::atan2(dety1, detx1);
1113 double theta2 = std::sqrt(detx2 * detx2 + dety2 * dety2);
1114 double phi2 = std::atan2(dety2, detx2);
1115 double theta3 = std::sqrt(detx3 * detx3 + dety3 * dety3);
1116 double phi3 = std::atan2(dety3, detx3);
1125 double a12 = dir1.
dist(dir2);
1126 double a13 = dir1.
dist(dir3);
1127 double a23 = dir2.
dist(dir3);
1130 double s = 0.5 * (a12 + a23 + a13);
1131 double solidangle = 4.0 * std::atan(std::sqrt(std::tan(0.5*s) *
1132 std::tan(0.5*(s-a12)) *
1133 std::tan(0.5*(s-a23)) *
1134 std::tan(0.5*(s-a13))));
CTA 2D background class definition.
CTA instrument direction class interface definition.
Definition of support function used by CTA classes.
Exception handler interface definition.
Filename class interface definition.
FITS binary table class definition.
FITS file class interface definition.
Mathematical function definitions.
Random number generator class definition.
double rate(const int &iebin, const double &theta) const
Return background rate for a given energy bin and offset angle value (events/s/MeV/sr)
void save(const GFilename &filename, const bool &clobber=false) const
Save background into FITS file.
std::string print(const GChatter &chatter=NORMAL) const
Print background information.
const GCTAResponseTable & table(void) const
Return response table.
int index(const int &itheta, const int &iebin) const
Return background rate bin index.
virtual ~GCTABackground2D(void)
Destructor.
GFilename filename(void) const
Return filename.
int m_inx_theta
THETA axis index.
bool is_valid(void) const
Return validity of background model.
GCTABackground2D(void)
Void constructor.
GFilename m_filename
Name of background response file.
void init_mc_max_rate(void) const
Initialise array of maximum background rate.
void init_members(void)
Initialise class members.
void clear(void)
Clear background.
GCTAResponseTable m_background
Background response table.
void read(const GFitsTable &table)
Read background from FITS table.
double m_theta_min
THETA minimum (radians)
void free_members(void)
Delete class members.
void write(GFitsBinTable &table) const
Write background into FITS table.
void copy_members(const GCTABackground2D &bgd)
Copy class members.
GEnergies m_energy
Vector of energies.
int m_num_energy
Number of energy bins.
double m_logE_min
Log10(E/TeV) minimum.
void init_mc_cache(void) const
Initialise Monte Carlo cache.
double m_theta_max
THETA maximum (radians)
GCTABackground2D & operator=(const GCTABackground2D &bgd)
Assignment operator.
GCTAInstDir mc(const GEnergy &energy, const GTime &time, GRan &ran) const
Returns MC instrument direction.
int m_num[2]
Array of number of bins.
double solid_angle(const double &detx1, const double &dety1, const double &detx2, const double &dety2, const double &detx3, const double &dety3) const
Compute solid angle of pixel wedge.
GCTABackground2D * clone(void) const
Clone background.
int m_inx_bgd
Background index.
void set_members(void)
Set members from background table.
int m_num_theta
Number of THETA bins.
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 .
int m_inx_energy
Energy axis index.
double m_logE_max
Log10(E/TeV) maximum.
std::vector< double > m_mc_max
Maximum background rate.
GModelSpectralNodes m_mc_spectrum
Response cube spectrum.
virtual double operator()(const double &logE, const double &detx, const double &dety) const
Return background rate in units of events MeV s sr .
void load(const GFilename &filename)
Load background from FITS file.
Abstract base class for the CTA background model.
void init_members(void)
Initialise class members.
GCTABackground & operator=(const GCTABackground &bgd)
Assignment operator.
void free_members(void)
Delete class members.
CTA instrument direction class.
void dety(const double &y)
Set DETY coordinate (in radians)
void detx(const double &x)
Set DETX coordinate (in radians)
GCTAResponseTable * clone(void) const
Clone response table.
void read(const GFitsTable &table)
Read response table from FITS table HDU.
const std::string & axis_lo_unit(const int &axis) const
Return lower bin boundary unit for axis.
int table(const std::string &name) const
Determine index of table.
const std::string & axis_hi_unit(const int &axis) const
Return upper bin boundary unit for axis.
void axis_log10(const int &axis)
Set nodes for a logarithmic (base 10) axis.
void axis_radians(const int &axis)
Set nodes for a radians axis.
int axis(const std::string &name) const
Determine index of an axis.
const int & axes(void) const
Return number of axes of the tables.
void write(GFitsTable &table) const
Write response table into FITS table HDU.
int axis_bins(const int &axis) const
Return number bins in an axis.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
const GNodeArray & axis_nodes(const int &axis) const
Return axis nodes.
void clear(void)
Clear response table.
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
void clear(void)
Clear energy container.
int size(void) const
Return number of energies in container.
void extend(const GEnergies &energies)
Append energy container.
GEnergy & append(const GEnergy &energy)
Append energy to container.
void remove(const int &index)
Remove energy from container.
Class that handles energies in a unit independent way.
double log10TeV(void) const
Return log10 of energy in TeV.
double MeV(void) const
Return energy in MeV.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
double TeV(void) const
Return energy in TeV.
std::string url(void) const
Return Uniform Resource Locator (URL)
std::string extname(const std::string &defaultname="") const
Return extension name.
void clear(void)
Clear file name.
Abstract interface for FITS table.
void remove(const int &colnum)
Remove column from the table.
bool contains(const int &extno) const
Check if HDU exists in FITS file.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
void close(void)
Close FITS file.
void remove(const int &extno)
Remove HDU from FITS file.
void save(const bool &clobber=false)
Saves FITS file.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
virtual void clear(void)
Clear spectral nodes model.
int nodes(void) const
Return number of nodes.
void append(const GEnergy &energy, const double &intensity)
Append node.
GEnergy energy(const int &index) const
Return node energy.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
const int & inx_right(void) const
Returns right node index.
const int & inx_left(void) const
Returns left node index.
const double & wgt_right(void) const
Returns right node weight.
const double & wgt_left(void) const
Returns left node weight.
Random number generator class.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
void radec(const double &ra, const double &dec)
Set equatorial sky direction (radians)
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
double plaw_integral(const double &x1, const double &f1, const double &x2, const double &f2)
Returns the integral of a power law.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
const std::string extname_cta_background2d
std::string gadf_hduclas4(const GFits &fits, const std::string &hduclas4)
Return extension name for GADF response table of given HDU class 4.
void warning(const std::string &origin, const std::string &message)
Emits warning.