Konstancja Satalecka / ctools

Commit d975bd1a3c9f1e6eec60bd1495a1c0135ea1ae1c

Authored by Jürgen Knödlseder
1 parent 2ed89b08

ctselect now respects any pre-existing energy or ROI cut. It automatically adjus… 

…ts the intervals, leading in the worst case to an empty event list if no overlap exists.
Inline Side-by-side
Showing 1 changed file with 254 additions and 154 deletions
src/ctselect/ctselect.cpp
  View file @d975bd1
... ... @@ -373,6 +373,92 @@ void ctselect::save(void)
373 373 }
374 374  
375 375  
  376 +/*==========================================================================
  377 + = =
  378 + = Private methods =
  379 + = =
  380 + ==========================================================================*/
  381 +
  382 +/***********************************************************************//**
  383 + * @brief Initialise class members
  384 + ***************************************************************************/
  385 +void ctselect::init_members(void)
  386 +{
  387 + // Initialise parameters
  388 + m_infile.clear();
  389 + m_outfile.clear();
  390 + m_prefix.clear();
  391 + m_usepnt = false;
  392 + m_ra = -1.0;
  393 + m_dec = -1.0;
  394 + m_rad = -1.0;
  395 + m_tmin = 0.0;
  396 + m_tmax = 0.0;
  397 + m_emin = 0.0;
  398 + m_emax = 0.0;
  399 + m_expr.clear();
  400 +
  401 + // Initialise protected members
  402 + m_obs.clear();
  403 + m_infiles.clear();
  404 + m_timemin.clear();
  405 + m_timemax.clear();
  406 + m_use_xml = false;
  407 +
  408 + // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD,
  409 + // which is defined in GCTALib.hpp. This is somehow a kluge. We need
  410 + // a better mechanism to implement the CTA reference MJD.
  411 + m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL");
  412 +
  413 + // Return
  414 + return;
  415 +}
  416 +
  417 +
  418 +/***********************************************************************//**
  419 + * @brief Copy class members
  420 + *
  421 + * @param[in] app Application.
  422 + ***************************************************************************/
  423 +void ctselect::copy_members(const ctselect& app)
  424 +{
  425 + // Copy parameters
  426 + m_infile = app.m_infile;
  427 + m_outfile = app.m_outfile;
  428 + m_prefix = app.m_prefix;
  429 + m_usepnt = app.m_usepnt;
  430 + m_ra = app.m_ra;
  431 + m_dec = app.m_dec;
  432 + m_rad = app.m_rad;
  433 + m_tmin = app.m_tmin;
  434 + m_tmax = app.m_tmax;
  435 + m_emin = app.m_emin;
  436 + m_emax = app.m_emax;
  437 + m_expr = app.m_expr;
  438 +
  439 + // Copy protected members
  440 + m_obs = app.m_obs;
  441 + m_infiles = app.m_infiles;
  442 + m_timemin = app.m_timemin;
  443 + m_timemax = app.m_timemax;
  444 + m_use_xml = app.m_use_xml;
  445 + m_cta_ref = app.m_cta_ref;
  446 +
  447 + // Return
  448 + return;
  449 +}
  450 +
  451 +
  452 +/***********************************************************************//**
  453 + * @brief Delete class members
  454 + ***************************************************************************/
  455 +void ctselect::free_members(void)
  456 +{
  457 + // Return
  458 + return;
  459 +}
  460 +
  461 +
376 462 /***********************************************************************//**
377 463 * @brief Get application parameters
378 464 *
... ... @@ -458,44 +544,55 @@ void ctselect::get_parameters(void)
458 544 * @param[in] obs CTA observation.
459 545 * @param[in] filename File name.
460 546 *
  547 + * @exception GException::invalid_value
  548 + * No EVENTS extension found in FITS file.
  549 + *
461 550 * Select events from a FITS file by making use of the selection possibility
462 551 * of the cfitsio library on loading a file. A selection string is created
463 552 * from the specified criteria that is appended to the filename so that
464 553 * cfitsio will automatically filter the event data. This selection string
465 554 * is then applied when opening the FITS file. The opened FITS file is then
466 555 * saved into a temporary file which is the loaded into the actual CTA
467   - * observation, overwriting the old CTA observation. The ROI, GTI and EBounds
  556 + * observation, overwriting the old CTA observation. The ROI, GTI and Ebounds
468 557 * of the CTA event list are then set accordingly to the specified selection.
469 558 * Finally, the temporary file created during this process is removed.
470 559 *
471 560 * Good Time Intervals of the observation will be limited to the time
472 561 * interval [m_tmin, m_tmax]. If m_tmin=m_tmax=0, no time selection is
473 562 * performed.
474   - *
475   - * @todo Use INDEF instead of 0.0 for pointing as RA/DEC selection
476 563 ***************************************************************************/
477 564 void ctselect::select_events(GCTAObservation* obs, const std::string& filename)
478 565 {
479   - // Allocate selection string
  566 + // Initialise selection string
480 567 std::string selection;
481   - char cmin[80];
482   - char cmax[80];
483   - char cra[80];
484   - char cdec[80];
485   - char crad[80];
486   -
487   - // Set requested selections
488   - bool select_time = (m_tmin != 0.0 || m_tmax != 0.0);
489   - bool select_roi = ((m_dec != -1.0 && m_ra != -1.0 && m_rad != -1.0) || m_usepnt);
  568 +
  569 + // Initialise selection flags
  570 + bool remove_all = false;
  571 +
  572 + // Analyse task parameters to see what selections are required
  573 + bool select_time = (m_tmin != 0.0 || m_tmax != 0.0);
  574 + bool select_energy = true;
  575 + bool select_roi = ((m_dec != -1.0 && m_ra != -1.0 && m_rad != -1.0) || m_usepnt);
  576 + bool select_expr = (gammalib::strip_whitespace(m_expr).length() > 0);
  577 +
  578 + // Set energy selection
  579 + double emin = m_emin;
  580 + double emax = m_emax;
490 581  
491   - // Set RA/DEC selection
  582 + // Set RA/DEC selection. If the "usepnt" parameter is set to true then
  583 + // use the pointing direction as the ROI centre.
492 584 double ra = m_ra;
493 585 double dec = m_dec;
  586 + double rad = m_rad;
494 587 if (m_usepnt) {
495 588 const GCTAPointing &pnt = obs->pointing();
496 589 ra = pnt.dir().ra_deg();
497 590 dec = pnt.dir().dec_deg();
498 591 }
  592 +
  593 + // Get CTA event list pointer
  594 + GCTAEventList* list =
  595 + static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
499 596  
500 597 // Set time selection interval. We make sure here that the time selection
501 598 // interval cannot be wider than the GTIs covering the data. This is done
... ... @@ -505,12 +602,12 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
505 602 // Reduce GTIs to specified time interval. The complicated cast is
506 603 // necessary here because the gti() method is declared const, so
507 604 // we're not officially allowed to modify the GTIs.
508   - ((GGti*)(&obs->events()->gti()))->reduce(m_timemin, m_timemax);
  605 + ((GGti*)(&list->gti()))->reduce(m_timemin, m_timemax);
509 606  
510 607 } // endif: time selection was required
511 608  
512 609 // Save GTI for later usage
513   - GGti gti = obs->events()->gti();
  610 + GGti gti = list->gti();
514 611  
515 612 // Make time selection
516 613 if (select_time) {
... ... @@ -521,6 +618,8 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
521 618 double tmax = gti.tstop().convert(gti.reference());
522 619  
523 620 // Format time with sufficient accuracy and add to selection string
  621 + char cmin[80];
  622 + char cmax[80];
524 623 sprintf(cmin, "%.8f", tmin);
525 624 sprintf(cmax, "%.8f", tmax);
526 625 selection = "TIME >= "+std::string(cmin)+" && TIME <= "+std::string(cmax);
... ... @@ -528,56 +627,147 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
528 627 log << gammalib::parformat("Time range");
529 628 log << tmin << " - " << tmax << " s" << std::endl;
530 629 }
  630 +
  631 + } // endif: made time selection
  632 +
  633 + // Make energy selection
  634 + if (select_energy) {
  635 +
  636 + // If we have aready a selection than add an "&&" operator
531 637 if (selection.length() > 0) {
532 638 selection += " && ";
533 639 }
534   - }
535 640  
536   - // Make energy selection
537   - sprintf(cmin, "%.8f", m_emin);
538   - sprintf(cmax, "%.8f", m_emax);
539   - selection += "ENERGY >= "+std::string(cmin)+" && ENERGY <= "+std::string(cmax);
540   - if (logTerse()) {
541   - log << gammalib::parformat("Energy range");
542   - log << m_emin << " - " << m_emax << " TeV" << std::endl;
543   - }
  641 + // Log the requested energy range
  642 + if (logTerse()) {
  643 + log << gammalib::parformat("Requested energy range");
  644 + log << m_emin << " - " << m_emax << " TeV" << std::endl;
  645 + }
  646 +
  647 + // Set the energy boundaries. Make sure that they overlap with any
  648 + // existing energy boundaries
  649 + const GEbounds& ebounds = list->ebounds();
  650 + if (ebounds.size() >= 1) {
  651 + if (logTerse()) {
  652 + log << gammalib::parformat("Energy range of data");
  653 + log << ebounds.emin(0).TeV() << " - ";
  654 + log << ebounds.emax(0).TeV() << " TeV" << std::endl;
  655 + }
  656 + if (emin < ebounds.emin(0).TeV()) {
  657 + emin = ebounds.emin(0).TeV();
  658 + }
  659 + if (emax > ebounds.emax(0).TeV()) {
  660 + emax = ebounds.emax(0).TeV();
  661 + }
  662 + }
  663 +
  664 + // Log the selected energy range
  665 + if (logTerse()) {
  666 + log << gammalib::parformat("Selected energy range");
  667 + if (emin >= emax) {
  668 + log << "None. There is no overlap between existing ";
  669 + log << "and requested energy range." << std::endl;
  670 + }
  671 + else {
  672 + log << emin << " - " << emax << " TeV" << std::endl;
  673 + }
  674 + }
  675 +
  676 + // Format energy selection
  677 + char cmin[80];
  678 + char cmax[80];
  679 + sprintf(cmin, "%.8f", emin);
  680 + sprintf(cmax, "%.8f", emax);
  681 + selection += "ENERGY >= "+std::string(cmin)+" && ENERGY <= "+std::string(cmax);
  682 +
  683 + } // endif: made energy selection
544 684  
545 685 // Make ROI selection
546   - if (select_roi){
  686 + if (select_roi) {
  687 +
  688 + // If we have aready a selection than add an "&&" operator
547 689 if (selection.length() > 0) {
548 690 selection += " && ";
549 691 }
  692 +
  693 + // Log the requested ROI
  694 + if (logTerse()) {
  695 + log << gammalib::parformat("Requested ROI");
  696 + log << "Centre(RA,DEC)=(" << ra << ", " << dec << ") deg, ";
  697 + log << "Radius=" << m_rad << " deg" << std::endl;
  698 + }
  699 +
  700 + // If we have already an ROI then make sure that the selected
  701 + // ROI overlaps with the existing ROI
  702 + double roi_radius = list->roi().radius();
  703 + if (roi_radius > 0.0) {
  704 + GSkyDir roi_centre = list->roi().centre().dir();
  705 + if (logTerse()) {
  706 + log << gammalib::parformat("ROI of data");
  707 + log << "Centre(RA,DEC)=(" << roi_centre.ra_deg() << ", ";
  708 + log << roi_centre.dec_deg() << ") deg, ";
  709 + log << "Radius=" << roi_radius << " deg" << std::endl;
  710 + }
  711 + GSkyDir centre;
  712 + centre.radec_deg(ra, dec);
  713 + double distance = centre.dist_deg(roi_centre);
  714 + if (distance + rad > roi_radius) {
  715 + rad = roi_radius - distance;
  716 + }
  717 + }
  718 +
  719 + // Log the selected ROI
  720 + if (logTerse()) {
  721 + log << gammalib::parformat("Selected ROI");
  722 + if (rad <= 0.0) {
  723 + log << "None. There is no overlap between existing ";
  724 + log << "and requested ROI." << std::endl;
  725 + }
  726 + else {
  727 + log << "Centre(RA,DEC)=(" << ra << ", " << dec << ") deg, ";
  728 + log << "Radius=" << rad << " deg" << std::endl;
  729 + }
  730 + }
550 731  
  732 + // Format ROI selection
  733 + char cra[80];
  734 + char cdec[80];
  735 + char crad[80];
551 736 sprintf(cra, "%.6f", ra);
552 737 sprintf(cdec, "%.6f", dec);
553   - sprintf(crad, "%.6f", m_rad);
  738 + sprintf(crad, "%.6f", rad);
554 739 selection += "ANGSEP("+std::string(cra)+"," +
555 740 std::string(cdec)+",RA,DEC) <= " +
556 741 std::string(crad);
557   - if (logTerse()) {
558   - log << gammalib::parformat("Acceptance cone centre");
559   - log << "RA=" << ra << ", DEC=" << dec << " deg" << std::endl;
560   - log << gammalib::parformat("Acceptance cone radius");
561   - log << m_rad << " deg" << std::endl;
562   - }
563   - if (logTerse()) {
564   - log << gammalib::parformat("cfitsio selection");
565   - log << selection << std::endl;
566   - }
567   - }
568 742  
569   - // Add additional expression
570   - if (gammalib::strip_whitespace(m_expr).length() > 0) {
  743 + } // endif: made ROI selection
  744 +
  745 + // Make an expression selection
  746 + if (select_expr) {
  747 +
  748 + // If we have aready a selection than add an "&&" operator
571 749 if (selection.length() > 0) {
572 750 selection += " && ";
573 751 }
  752 +
  753 + // Append the expression
574 754 selection += "("+gammalib::strip_whitespace(m_expr)+")";
  755 +
  756 + } // endif: made expression selection
  757 +
  758 + // Dump cfitsio selection string
  759 + if (logTerse()) {
  760 + log << gammalib::parformat("cfitsio selection");
  761 + log << selection << std::endl;
575 762 }
576 763  
577 764 // Build input filename including selection expression
578 765 std::string expression = filename;
579   - if (selection.length() > 0)
  766 + if (selection.length() > 0) {
580 767 expression += "[EVENTS]["+selection+"]";
  768 + }
  769 +
  770 + // Dump FITS filename including selection expression
581 771 if (logTerse()) {
582 772 log << gammalib::parformat("FITS filename");
583 773 log << expression << std::endl;
... ... @@ -595,17 +785,17 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
595 785  
596 786 // Check if we have an EVENTS HDU
597 787 if (!file.contains("EVENTS")) {
598   - std::string message = "No \"EVENTS\" extension found in FITS file "+
599   - expression+".";
600   - throw GException::app_error(G_SELECT_EVENTS, message);
  788 + std::string msg = "No \"EVENTS\" extension found in FITS file \""+
  789 + expression+"\".";
  790 + throw GException::invalid_value(G_SELECT_EVENTS, msg);
601 791 }
602 792  
603 793 // Determine number of events in EVENTS HDU
604 794 int nevents = file.table("EVENTS")->nrows();
605 795  
606   - // If the selected event list is empty then append an empty event list
607   - // to the observation. Otherwise load the data from the temporary file.
608   - if (nevents < 1) {
  796 + // If the selected event list is empty or if removal of all events
  797 + // has been requesten then append an empty event list to the observation.
  798 + if ((nevents < 1) || (remove_all)) {
609 799  
610 800 // Create empty event list
611 801 GCTAEventList eventlist;
... ... @@ -614,127 +804,37 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
614 804 obs->events(eventlist);
615 805  
616 806 }
617   - else {
618 807  
619   - // Re-read observation from FITS object
  808 + // ... otherwise load the data from the temporary file
  809 + else {
620 810 obs->read(file);
621   -
622 811 }
623 812  
624 813 // Get CTA event list pointer
625   - GCTAEventList* list =
626   - static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
627   -
628   - if (select_roi){
629   - // Set ROI
  814 + list = static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
  815 +
  816 + // If ROI selection has been applied then set the event list ROI
  817 + if (select_roi) {
630 818 GCTAInstDir instdir;
631 819 instdir.dir().radec_deg(ra, dec);
632   - list->roi(GCTARoi(instdir, m_rad));
  820 + list->roi(GCTARoi(instdir, rad));
633 821 }
634 822  
635   - // Set GTI
  823 + // Set event list GTI (in any case as any event list has a GTI)
636 824 list->gti(gti);
637 825  
638   - // Set energy boundaries
639   - GEbounds ebounds;
640   - GEnergy emin;
641   - GEnergy emax;
642   - emin.TeV(m_emin);
643   - emax.TeV(m_emax);
644   - ebounds.append(emin, emax);
645   - list->ebounds(ebounds);
  826 + // If an energy selection has been applied then set the energy boundaries
  827 + if (select_energy) {
  828 + GEbounds ebounds;
  829 + ebounds.append(GEnergy(emin, "TeV"), GEnergy(emax, "TeV"));
  830 + list->ebounds(ebounds);
  831 + }
646 832  
647 833 // Recompute ontime and livetime.
648   - GTime meantime = 0.5 * (gti.tstart() + gti.tstop());
649   - obs->ontime(gti.ontime());
650   - obs->livetime(gti.ontime() * obs->deadc(meantime));
651   -
652   - // Return
653   - return;
654   -}
655   -
  834 + GTime meantime = 0.5 * (list->gti().tstart() + list->gti().tstop());
  835 + obs->ontime(list->gti().ontime());
  836 + obs->livetime(list->gti().ontime() * obs->deadc(meantime));
656 837  
657   -/*==========================================================================
658   - = =
659   - = Private methods =
660   - = =
661   - ==========================================================================*/
662   -
663   -/***********************************************************************//**
664   - * @brief Initialise class members
665   - ***************************************************************************/
666   -void ctselect::init_members(void)
667   -{
668   - // Initialise parameters
669   - m_infile.clear();
670   - m_outfile.clear();
671   - m_prefix.clear();
672   - m_usepnt = false;
673   - m_ra = -1.0;
674   - m_dec = -1.0;
675   - m_rad = -1.0;
676   - m_tmin = 0.0;
677   - m_tmax = 0.0;
678   - m_emin = 0.0;
679   - m_emax = 0.0;
680   - m_expr.clear();
681   -
682   - // Initialise protected members
683   - m_obs.clear();
684   - m_infiles.clear();
685   - m_timemin.clear();
686   - m_timemax.clear();
687   - m_use_xml = false;
688   -
689   - // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD,
690   - // which is defined in GCTALib.hpp. This is somehow a kluge. We need
691   - // a better mechanism to implement the CTA reference MJD.
692   - m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL");
693   -
694   - // Return
695   - return;
696   -}
697   -
698   -
699   -/***********************************************************************//**
700   - * @brief Copy class members
701   - *
702   - * @param[in] app Application.
703   - ***************************************************************************/
704   -void ctselect::copy_members(const ctselect& app)
705   -{
706   - // Copy parameters
707   - m_infile = app.m_infile;
708   - m_outfile = app.m_outfile;
709   - m_prefix = app.m_prefix;
710   - m_usepnt = app.m_usepnt;
711   - m_ra = app.m_ra;
712   - m_dec = app.m_dec;
713   - m_rad = app.m_rad;
714   - m_tmin = app.m_tmin;
715   - m_tmax = app.m_tmax;
716   - m_emin = app.m_emin;
717   - m_emax = app.m_emax;
718   - m_expr = app.m_expr;
719   -
720   - // Copy protected members
721   - m_obs = app.m_obs;
722   - m_infiles = app.m_infiles;
723   - m_timemin = app.m_timemin;
724   - m_timemax = app.m_timemax;
725   - m_use_xml = app.m_use_xml;
726   - m_cta_ref = app.m_cta_ref;
727   -
728   - // Return
729   - return;
730   -}
731   -
732   -
733   -/***********************************************************************//**
734   - * @brief Delete class members
735   - ***************************************************************************/
736   -void ctselect::free_members(void)
737   -{
738 838 // Return
739 839 return;
740 840 }
... ...