Josh Cardenzana / gammalib

Commit 3f9c60f8c57874aed836795fe5841c4bfd0e5d6c

Authored by Luigi Tibaldo
1 parent e1ae67f5

special case for point sources

Inline Side-by-side
Showing 2 changed files with 76 additions and 49 deletions
src/model/GModelSpatial.cpp
  View file @3f9c60f
... ... @@ -30,6 +30,7 @@
30 30 #endif
31 31 #include "GException.hpp"
32 32 #include "GModelSpatial.hpp"
  33 +#include "GModelSpatialPointSource.hpp"
33 34  
34 35 /* __ Method name definitions ____________________________________________ */
35 36 #define G_ACCESS "GModelSpatial::operator[](std::string&)"
... ... @@ -347,63 +348,81 @@ double GModelSpatial::flux(const GSkyRegion* reg,
347 348 const GSkyRegion* model_reg = this->region();
348 349 if (model_reg->overlaps(*reg)) {
349 350  
350   - // Thrown an exception if region is not a sky circle
351   - const GSkyRegionCircle* reg_circle =
352   - dynamic_cast<const GSkyRegionCircle*>(reg);
353   - if (reg_circle == NULL) {
  351 + // Check if the model is a point source
  352 + const GModelSpatialPointSource* ps_model =
  353 + dynamic_cast<const GModelSpatialPointSource*>(this);
  354 +
  355 + // If the model is not a point source integrate over circular regions overlap
  356 + if (ps_model == NULL) {
  357 +
  358 + // Thrown an exception if region is not a sky circle
  359 + const GSkyRegionCircle* reg_circle =
  360 + dynamic_cast<const GSkyRegionCircle*>(reg);
  361 + if (reg_circle == NULL) {
354 362 std::string msg = "Flux can only be computed for a circular "
355   - "region.";
  363 + "region.";
356 364 throw GException::feature_not_implemented(G_FLUX,msg);
357   - }
  365 + }
358 366  
359   - // Thrown an exception if model region is not a sky circle
360   - const GSkyRegionCircle* model_circle =
361   - dynamic_cast<const GSkyRegionCircle*>(model_reg);
362   - if (model_circle == NULL) {
  367 + // Thrown an exception if model region is not a sky circle
  368 + const GSkyRegionCircle* model_circle =
  369 + dynamic_cast<const GSkyRegionCircle*>(model_reg);
  370 + if (model_circle == NULL) {
363 371 std::string msg = "Flux can only be computed for spatial model "
364   - "with region defined as circle.";
  372 + "with region defined as circle.";
365 373 throw GException::feature_not_implemented(G_FLUX,msg);
366   - }
  374 + }
367 375  
368   - // Model centre and radius in radians
369   - GSkyDir model_centre = model_circle->centre();
370   - double model_radius = model_circle->radius() * gammalib::deg2rad;
  376 + // Model centre and radius in radians
  377 + GSkyDir model_centre = model_circle->centre();
  378 + double model_radius = model_circle->radius() * gammalib::deg2rad;
371 379  
372   - // Distance between region and model centres in radians
373   - double distance = model_centre.dist(reg_circle->centre());
  380 + // Distance between region and model centres in radians
  381 + double distance = model_centre.dist(reg_circle->centre());
374 382  
375   - // Take the distance between centers minus radius of model region as
376   - // the minimum radial integration boundary (in radians)
377   - double rho_min = distance - model_radius;
  383 + // Take the distance between centers minus radius of model region as
  384 + // the minimum radial integration boundary (in radians)
  385 + double rho_min = distance - model_radius;
378 386  
379   - // Make sure that rho_min does not become negative
380   - if (rho_min < 0.0) {
  387 + // Make sure that rho_min does not become negative
  388 + if (rho_min < 0.0) {
381 389 rho_min = 0.0;
382   - }
383   -
384   - // Take the region radius as the maximum radial integration boundary
385   - // (in radians)
386   - double rho_max = reg_circle->radius() * gammalib::deg2rad;
387   -
388   - // Pre-compute some quantities for arclength
389   - double cosdist = std::cos(distance);
390   - double sindist = std::sin(distance);
391   - double cosmodrad = std::cos(model_radius);
392   -
393   - // Setup integration kernel
394   - GModelSpatial::circle_int_kern_rho integrand(this,
395   - reg_circle,
396   - srcEng,
397   - srcTime,
398   - distance,
399   - cosdist,
400   - sindist,
401   - model_radius,
402   - cosmodrad);
403   - GIntegral integral(&integrand);
404   -
405   - // Perform integration
406   - flux = integral.romberg(rho_min, rho_max);
  390 + }
  391 +
  392 + // Take the region radius as the maximum radial integration boundary
  393 + // (in radians)
  394 + double rho_max = reg_circle->radius() * gammalib::deg2rad;
  395 +
  396 + // Pre-compute some quantities for arclength
  397 + double cosdist = std::cos(distance);
  398 + double sindist = std::sin(distance);
  399 + double cosmodrad = std::cos(model_radius);
  400 +
  401 + // Setup integration kernel
  402 + GModelSpatial::circle_int_kern_rho integrand(this,
  403 + reg_circle,
  404 + srcEng,
  405 + srcTime,
  406 + distance,
  407 + cosdist,
  408 + sindist,
  409 + model_radius,
  410 + cosmodrad);
  411 + GIntegral integral(&integrand);
  412 +
  413 + // Perform integration
  414 + flux = integral.romberg(rho_min, rho_max);
  415 +
  416 + }
  417 +
  418 + // Otherwise, if model is point source since it overlaps
  419 + // with region set flux to 1
  420 + else {
  421 +
  422 + flux = 1.0;
  423 +
  424 + } // endif, model is point source
  425 +
407 426  
408 427 } // endif: model overlaps with region
409 428  
... ...
test/test_GModel.cpp
  View file @3f9c60f
... ... @@ -3472,14 +3472,14 @@ void TestGModel::test_flux(void)
3472 3472 centre1.radec_deg(0.0, 0.5);
3473 3473 GModelSpatialRadialDisk model = GModelSpatialRadialDisk(centre1, 0.25);
3474 3474  
3475   - // Define entre of test ROI
  3475 + // Define centre of test ROI
3476 3476 GSkyDir centre2;
3477 3477 centre2.radec_deg(0.0, 0.0);
3478 3478  
3479 3479 // Test non overlapping region
3480 3480 GSkyRegionCircle circle1 = GSkyRegionCircle(centre2, 0.1);
3481 3481 GSkyRegion* roi1 = &circle1;
3482   - test_value(model.flux(roi1), 0.0, 1.e-4, "Check model not overlapping with region");
  3482 + test_value(model.flux(roi1), 0.0, 1.e-7, "Check model not overlapping with region");
3483 3483  
3484 3484 // Test fully contained model
3485 3485 GSkyRegionCircle circle2 = GSkyRegionCircle(centre2, 1.0);
... ... @@ -3491,6 +3491,14 @@ void TestGModel::test_flux(void)
3491 3491 GSkyRegion* roi3 = &circle3;
3492 3492 test_value(model.flux(roi3), 0.44654, 1.e-4, "Check model partially contained in region");
3493 3493  
  3494 + // Test point source in ROI
  3495 + GModelSpatialPointSource ps = GModelSpatialPointSource(centre1);
  3496 + test_value(ps.flux(roi2), 1.0, 1.e-7, "Check point source contained in region");
  3497 +
  3498 + // Test point source outside ROI
  3499 + test_value(ps.flux(roi1), 0.0, 1.e-7, "Check point source not contained in region");
  3500 +
  3501 +
3494 3502 // Exit test
3495 3503 return;
3496 3504 }
... ...