Commit 5fba3818bb7c930f5c28735bd82d0902ca9c587d

Authored by Michael Mayer
1 parent 83984a06

added two constructors to GCTAMeanPsf and GCTAExposure

inst/cta/include/GCTAExposure.hpp
... ... @@ -34,6 +34,7 @@
34 34 #include "GSkymap.hpp"
35 35 #include "GObservations.hpp"
36 36 #include "GCTAObservation.hpp"
  37 +#include "GCTAEventCube.hpp"
37 38  
38 39  
39 40 /***********************************************************************//**
... ... @@ -53,6 +54,7 @@ public:
53 54 GCTAExposure(void);
54 55 GCTAExposure(const GCTAExposure& cube);
55 56 explicit GCTAExposure(const std::string& filename);
  57 + GCTAExposure(const GCTAEventCube& cube);
56 58 GCTAExposure(const std::string& wcs,
57 59 const std::string& coords,
58 60 const double& x,
... ...
inst/cta/include/GCTAMeanPsf.hpp
... ... @@ -35,7 +35,7 @@
35 35 #include "GObservations.hpp"
36 36 #include "GCTAObservation.hpp"
37 37 #include "GNodeArray.hpp"
38   -
  38 +#include "GCTAEventCube.hpp"
39 39  
40 40 /***********************************************************************//**
41 41 * @class GCTAMeanPsf
... ... @@ -55,6 +55,7 @@ public:
55 55 GCTAMeanPsf(void);
56 56 GCTAMeanPsf(const GCTAMeanPsf& cube);
57 57 explicit GCTAMeanPsf(const std::string& filename);
  58 + GCTAMeanPsf(const GCTAEventCube& cube, const double& dmax, const int& ndbins);
58 59 GCTAMeanPsf(const std::string& wcs,
59 60 const std::string& coords,
60 61 const double& x,
... ...
inst/cta/pyext/GCTAExposure.i
... ... @@ -46,6 +46,7 @@ public:
46 46 GCTAExposure(void);
47 47 GCTAExposure(const GCTAExposure& cube);
48 48 explicit GCTAExposure(const std::string& filename);
  49 + GCTAExposure(const GCTAEventCube& cube);
49 50 GCTAExposure(const std::string& wcs,
50 51 const std::string& coords,
51 52 const double& x,
... ...
inst/cta/pyext/GCTAMeanPsf.i
... ... @@ -49,6 +49,7 @@ public:
49 49 GCTAMeanPsf(void);
50 50 GCTAMeanPsf(const GCTAMeanPsf& cube);
51 51 explicit GCTAMeanPsf(const std::string& filename);
  52 + GCTAMeanPsf(const GCTAEventCube& cube, const double& dmax, const int& ndbins);
52 53 GCTAMeanPsf(const std::string& wcs,
53 54 const std::string& coords,
54 55 const double& x,
... ...
inst/cta/src/GCTAExposure.cpp
... ... @@ -104,6 +104,25 @@ GCTAExposure::GCTAExposure(const std::string& filename)
104 104 return;
105 105 }
106 106  
  107 +GCTAExposure::GCTAExposure(const GCTAEventCube& cube)
  108 +{
  109 + // Initialise class members
  110 + init_members();
  111 +
  112 + // Store energy boundaries
  113 + m_ebounds = cube.ebounds();
  114 +
  115 + // Set GNodeArray used for interpolation
  116 + set_eng_axis();
  117 +
  118 + // Get sky map
  119 + m_cube = cube.map();
  120 +
  121 + // Return
  122 + return;
  123 +
  124 +}
  125 +
107 126  
108 127 /***********************************************************************//**
109 128 * @brief Exposure cube constructor
... ...
inst/cta/src/GCTAMeanPsf.cpp
... ... @@ -35,6 +35,8 @@
35 35 #include "GCTAExposure.hpp"
36 36 #include "GMath.hpp"
37 37 #include "GTools.hpp"
  38 +#include "GWcsRegistry.hpp"
  39 +#include "GWcs.hpp"
38 40  
39 41 /* __ Method name definitions ____________________________________________ */
40 42 #define G_SET "GCTAMeanPsf::set(GCTAObservation&)"
... ... @@ -106,6 +108,52 @@ GCTAMeanPsf::GCTAMeanPsf(const std::string& filename)
106 108 return;
107 109 }
108 110  
  111 +GCTAMeanPsf::GCTAMeanPsf(const GCTAEventCube& cube, const double& dmax, const int& ndbins)
  112 +{
  113 + // Initialise class members
  114 + init_members();
  115 +
  116 + // Store energy boundaries
  117 + m_ebounds = cube.ebounds();
  118 +
  119 + // Set GNodeArray used for interpolation
  120 + set_eng_axis();
  121 +
  122 + // Set delta node array
  123 +
  124 + m_deltas.clear();
  125 + for (int i = 0; i < ndbins; ++i) {
  126 + #if defined(G_QUADRATIC_BINNING)
  127 + double binsize = std::sqrt(dmax) / double(ndbins);
  128 + double delta = binsize * (double(i) + 0.5); // avoid central singularity
  129 + delta *= delta;
  130 + #else
  131 + double binsize = dmax / double(ndbins);
  132 + double delta = binsize * (double(i) + 0.5); // avoid central singularity
  133 + #endif
  134 + m_deltas.append(delta);
  135 + }
  136 +
  137 + // Compute number of sky maps
  138 + int nmaps = m_ebounds.size() * m_deltas.size();
  139 +
  140 +
  141 + // Create sky map
  142 + const GWcs* wcs = dynamic_cast<const GWcs*>(cube.map().projection());// reg.alloc(cube.map().projection()->code());
  143 +
  144 + std::string coordsys = wcs->coordsys();
  145 + std::string proj = wcs->code();
  146 + double xref = wcs->crval(0);
  147 + double yref = wcs->crval(1);
  148 + double xbinsz = wcs->cdelt(0);
  149 + double ybinsz = wcs->cdelt(1);
  150 + m_cube = GSkymap(proj, coordsys, xref, yref, xbinsz, ybinsz, cube.map().nx(), cube.map().ny(), nmaps);
  151 +
  152 + // Return
  153 + return;
  154 +
  155 +}
  156 +
109 157  
110 158 /***********************************************************************//**
111 159 * @brief Mean PSF cube constructor
... ... @@ -431,7 +479,7 @@ void GCTAMeanPsf::fill(const GObservations&amp; obs)
431 479  
432 480 // Set map index
433 481 int imap = offset(idelta, iebin);
434   -
  482 +
435 483 // Add on PSF cube
436 484 m_cube(pixel, imap) +=
437 485 rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
... ...