Commit 429914889e59fce3cb1b3d7587cbc1a2f56a2bee

Authored by Jürgen Knödlseder
1 parent db071294

Increase sampling of area computation in ctobssim

Showing 1 changed file with 14 additions and 6 deletions
src/ctobssim/ctobssim.cpp
@@ -1229,7 +1229,10 @@ GEbounds ctobssim::get_ebounds(const GEbounds& ebounds) const @@ -1229,7 +1229,10 @@ GEbounds ctobssim::get_ebounds(const GEbounds& ebounds) const
1229 * @param[in] emax Maximum energy. 1229 * @param[in] emax Maximum energy.
1230 * @return Simulation area (cm^2). 1230 * @return Simulation area (cm^2).
1231 * 1231 *
1232 - * Get the simulation area for an energy interval in units of cm^2. 1232 + * Get the simulation area for an energy interval in units of cm^2. This is
  1233 + * done by extracting the maximum effective area value within the energy
  1234 + * range [emin,emax] and by multiplying this value by 2 for security. The
  1235 + * effective area is sampled at 10 energy values within the energy interval.
1233 ***************************************************************************/ 1236 ***************************************************************************/
1234 double ctobssim::get_area(GCTAObservation* obs, 1237 double ctobssim::get_area(GCTAObservation* obs,
1235 const GEnergy& emin, 1238 const GEnergy& emin,
@@ -1247,11 +1250,16 @@ double ctobssim::get_area(GCTAObservation* obs, @@ -1247,11 +1250,16 @@ double ctobssim::get_area(GCTAObservation* obs,
1247 } 1250 }
1248 1251
1249 // Compute effective area at minimum and maximum energy 1252 // Compute effective area at minimum and maximum energy
1250 - double area_emin = rsp->aeff()->max(emin.log10TeV(), 0.0, 0.0);  
1251 - double area_emax = rsp->aeff()->max(emax.log10TeV(), 0.0, 0.0);  
1252 -  
1253 - // Use maximum of both  
1254 - double area = (area_emin > area_emax) ? area_emin : area_emax; 1253 + const int nbins = 10;
  1254 + double logE = emin.log10TeV();
  1255 + double logEbin = (emax.log10TeV() - logE)/double(nbins-1);
  1256 + double area = 0.0;
  1257 + for (int i = 0; i < nbins; ++i, logE += logEbin) {
  1258 + double aeff = rsp->aeff()->max(logE, 0.0, 0.0);
  1259 + if (aeff > area) {
  1260 + area = aeff;
  1261 + }
  1262 + }
1255 1263
1256 // Multiply by security factor 1264 // Multiply by security factor
1257 area *= 2.0; 1265 area *= 2.0;