Francois Brun / ctools

Commit 28b601b25f583537d41dc0b9ce178fbd5d1a1d31

Authored by Jürgen Knödlseder
1 parent 63ac93b6

Remove prefix parameters and make sure that bins outside RoI are not filled. 

We do not want to fill cube bins outside the RoI to be consistent with the response and modelling part.
Inline Side-by-side
Showing 3 changed files with 142 additions and 546 deletions
src/ctbin/ctbin.cpp
  View file @28b601b
... ... @@ -33,15 +33,13 @@
33 33 #include "GTools.hpp"
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
36   -#define G_BIN_EVENTS "ctbin::bin_events(GCTAObservation*)"
37   -#define G_GET_EBOUNDS "ctbin::get_ebounds()"
38 36 #define G_INIT_CUBE "ctbin::init_cube()"
39 37 #define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)"
  38 +#define G_GET_EBOUNDS "ctbin::get_ebounds()"
40 39  
41 40 /* __ Debug definitions __________________________________________________ */
42 41  
43 42 /* __ Coding definitions _________________________________________________ */
44   -#define G_MERGE_EVENTS
45 43  
46 44  
47 45 /*==========================================================================
... ... @@ -154,7 +152,7 @@ ctbin::~ctbin(void)
154 152 *
155 153 * @param[in] app Application.
156 154 ***************************************************************************/
157   -ctbin& ctbin::operator= (const ctbin& app)
  155 +ctbin& ctbin::operator=(const ctbin& app)
158 156 {
159 157 // Execute only if object is not identical
160 158 if (this != &app) {
... ... @@ -274,20 +272,12 @@ void ctbin::run(void)
274 272 }
275 273 }
276 274  
277   - // Initialise observation counter
278   - int n_observations = 0;
279   -
280   - // Compile option: if cube merging is requested, initialise the cube
281   - #if defined(G_MERGE_EVENTS)
  275 + // Initialise counts cube
282 276 init_cube();
283   - #endif
284 277  
285 278 // Loop over all observations in the container
286 279 for (int i = 0; i < m_obs.size(); ++i) {
287 280  
288   - // Initialise event input and output filenames
289   - m_infiles.push_back("");
290   -
291 281 // Get CTA observation
292 282 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
293 283  
... ... @@ -304,36 +294,15 @@ void ctbin::run(void)
304 294 }
305 295 }
306 296  
307   - // Increment number of observations
308   - n_observations++;
309   -
310   - // Save event file name (for possible saving)
311   - m_infiles[i] = obs->eventfile();
312   -
313   - // Compile option: if cube merging is requested, fill the cube;
314   - // otherwise just bin the events
315   - #if defined(G_MERGE_EVENTS)
  297 + // Fill the cube
316 298 fill_cube(obs);
317   - #else
318   - bin_events(obs);
319   - #endif
320 299  
321 300 } // endif: CTA observation found
322 301  
323 302 } // endfor: looped over observations
324 303  
325   - // Compile option: if cube merging is requested, set a single cube in
326   - // the observation container
327   - #if defined(G_MERGE_EVENTS)
  304 + // Set a single cube in the observation container
328 305 obs_cube();
329   - n_observations = 1;
330   - #endif
331   -
332   - // If more than a single observation has been handled then make sure
333   - // that an XML file will be used for storage
334   - if (n_observations > 1) {
335   - m_use_xml = true;
336   - }
337 306  
338 307 // Write observation(s) into logger
339 308 if (logTerse()) {
... ... @@ -353,21 +322,9 @@ void ctbin::run(void)
353 322  
354 323  
355 324 /***********************************************************************//**
356   - * @brief Save counts map(s)
357   - *
358   - * This method saves the counts map(s) into FITS file(s). There are two
359   - * modes, depending on the m_use_xml flag.
360   - *
361   - * If m_use_xml is true, all counts map(s) will be saved into FITS files,
362   - * where the output filenames are constructued from the input filenames by
363   - * prepending the m_prefix string to name. Any path information will be
364   - * stripped form the input name, hence event files will be written into the
365   - * local working directory (unless some path information is present
366   - * in the prefix). In addition, an XML file will be created that gathers
367   - * the filename information for the counts map(s). If an XML file was present
368   - * on input, all metadata information will be copied from this input file.
  325 + * @brief Save counts cube
369 326 *
370   - * If m_use_xml is false, the counts map will be saved into a FITS file.
  327 + * This method saves the counts cube into a FITS file.
371 328 ***************************************************************************/
372 329 void ctbin::save(void)
373 330 {
... ... @@ -382,26 +339,114 @@ void ctbin::save(void)
382 339 }
383 340 }
384 341  
385   - // Compile option: if merging was requested, always save into a single
386   - // FITS file
387   - #if defined(G_MERGE_EVENTS)
388   - save_fits();
389   - #else
  342 + // Get output filename
  343 + m_outfile = (*this)["outfile"].filename();
390 344  
391   - /*
392   - // Case A: Save counts map(s) and XML metadata information
393   - if (m_use_xml) {
394   - save_xml();
  345 + // Get CTA observation from observation container
  346 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
  347 +
  348 + // Save only if observation is valid
  349 + if (obs != NULL) {
  350 + obs->save(m_outfile, clobber());
395 351 }
396   - // Case B: Save counts map as FITS file
397   - else {
398   - save_fits();
  352 +
  353 + // Return
  354 + return;
  355 +}
  356 +
  357 +
  358 +/*==========================================================================
  359 + = =
  360 + = Private methods =
  361 + = =
  362 + ==========================================================================*/
  363 +
  364 +/***********************************************************************//**
  365 + * @brief Initialise class members
  366 + ***************************************************************************/
  367 +void ctbin::init_members(void)
  368 +{
  369 + // Initialise members
  370 + m_evfile.clear();
  371 + m_outfile.clear();
  372 + m_proj.clear();
  373 + m_coordsys.clear();
  374 + m_ebinalg.clear();
  375 + m_ebinfile.clear();
  376 + m_usepnt = false;
  377 + m_emin = 0.0;
  378 + m_emax = 0.0;
  379 + m_enumbins = 0;
  380 + m_xref = 0.0;
  381 + m_yref = 0.0;
  382 + m_binsz = 0.0;
  383 + m_nxpix = 0;
  384 + m_nypix = 0;
  385 +
  386 + // Initialise protected members
  387 + m_obs.clear();
  388 + m_read_ahead = false;
  389 + m_cube.clear();
  390 + m_ebounds.clear();
  391 + m_gti.clear();
  392 + m_ontime = 0.0;
  393 + m_livetime = 0.0;
  394 +
  395 + // Set logger properties
  396 + log.date(true);
  397 +
  398 + // Return
  399 + return;
  400 +}
  401 +
  402 +
  403 +/***********************************************************************//**
  404 + * @brief Copy class members
  405 + *
  406 + * @param[in] app Application.
  407 + ***************************************************************************/
  408 +void ctbin::copy_members(const ctbin& app)
  409 +{
  410 + // Copy attributes
  411 + m_evfile = app.m_evfile;
  412 + m_outfile = app.m_outfile;
  413 + m_proj = app.m_proj;
  414 + m_coordsys = app.m_coordsys;
  415 + m_ebinalg = app.m_ebinalg;
  416 + m_ebinfile = app.m_ebinfile;
  417 + m_usepnt = app.m_usepnt;
  418 + m_emin = app.m_emin;
  419 + m_emax = app.m_emax;
  420 + m_enumbins = app.m_enumbins;
  421 + m_xref = app.m_xref;
  422 + m_yref = app.m_yref;
  423 + m_binsz = app.m_binsz;
  424 + m_nxpix = app.m_nxpix;
  425 + m_nypix = app.m_nypix;
  426 +
  427 + // Copy protected members
  428 + m_obs = app.m_obs;
  429 + m_read_ahead = app.m_read_ahead;
  430 + m_cube = app.m_cube;
  431 + m_ebounds = app.m_ebounds;
  432 + m_gti = app.m_gti;
  433 + m_ontime = app.m_ontime;
  434 + m_livetime = app.m_livetime;
  435 +
  436 + // Return
  437 + return;
  438 +}
  439 +
  440 +
  441 +/***********************************************************************//**
  442 + * @brief Delete class members
  443 + ***************************************************************************/
  444 +void ctbin::free_members(void)
  445 +{
  446 + // Write separator into logger
  447 + if (logTerse()) {
  448 + log << std::endl;
399 449 }
400   - */
401   -
402   - // Always save FITS and an XML summary
403   - save_xml();
404   - #endif
405 450  
406 451 // Return
407 452 return;
... ... @@ -438,9 +483,6 @@ void ctbin::get_parameters(void)
438 483 // Append CTA observation to container
439 484 m_obs.append(obs);
440 485  
441   - // Signal that no XML file should be used for storage
442   - m_use_xml = false;
443   -
444 486 }
445 487  
446 488 // ... otherwise try to open as XML file
... ... @@ -449,9 +491,6 @@ void ctbin::get_parameters(void)
449 491 // Load observations from XML file
450 492 m_obs.load(m_evfile);
451 493  
452   - // Signal that XML file should be used for storage
453   - m_use_xml = true;
454   -
455 494 }
456 495  
457 496 } // endif: there was no observation in the container
... ... @@ -485,7 +524,6 @@ void ctbin::get_parameters(void)
485 524 // dumped into the log file
486 525 if (m_read_ahead) {
487 526 m_outfile = (*this)["outfile"].filename();
488   - m_prefix = (*this)["prefix"].string();
489 527 }
490 528  
491 529 // Return
... ... @@ -494,203 +532,6 @@ void ctbin::get_parameters(void)
494 532  
495 533  
496 534 /***********************************************************************//**
497   - * @brief Bin events into a counts map
498   - *
499   - * @param[in] obs CTA observation.
500   - *
501   - * @exception GException::no_list
502   - * No event list found in observation.
503   - *
504   - * This method bins the events found in a CTA events list into a counts map
505   - * and replaces the event list by the counts map in the observation. The
506   - * energy boundaries of the counts map are also stored in the observation's
507   - * energy boundary member.
508   - *
509   - * If the reference values for the map centre (m_xref, m_yref) are 9999.0,
510   - * the pointing direction of the observation is taken as the map centre.
511   - * Otherwise, the specified reference value is used.
512   - ***************************************************************************/
513   -void ctbin::bin_events(GCTAObservation* obs)
514   -{
515   - // Continue only if observation pointer is valid
516   - if (obs != NULL) {
517   -
518   - // Make sure that the observation holds a CTA event list. If this
519   - // is not the case then throw an exception.
520   - if (dynamic_cast<const GCTAEventList*>(obs->events()) == NULL) {
521   - throw GException::no_list(G_BIN_EVENTS);
522   - }
523   -
524   - // Initialise empty Ebounds object
525   - GEbounds ebds;
526   -
527   - // First check wether we have FITS-file defining
528   - // the binning
529   - if (m_ebinalg == "FILE") {
530   -
531   - // Open fits file to check which extension is given
532   - GFits file(m_ebinfile);
533   -
534   - if (file.contains("EBOUNDS")) {
535   -
536   - // Close file and load ebounds
537   - file.close();
538   - ebds.load(m_ebinfile,"EBOUNDS");
539   -
540   - } // endif: EBOUNDS extension was given
541   -
542   - else if (file.contains("ENERGYBINS")) {
543   -
544   - // Close file and load ebounds
545   - file.close();
546   - ebds.load(m_ebinfile,"ENERGYBINS");
547   -
548   - } // endelse: ENERGYBINS extension was given
549   -
550   - else {
551   -
552   - // Close file
553   - file.close();
554   -
555   - // Signal that no energy boundaries have been found
556   - std::string msg = "No extension with name \"EBOUNDS\" or"
557   - " \"ENERGYBINS\" found in FITS file"
558   - " \""+m_ebinfile+"\".\n"
559   - "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
560   - " is required if the parameter \"ebinalg\""
561   - " is set to \"FILE\".";
562   - throw GException::invalid_value(G_BIN_EVENTS, msg);
563   -
564   - } // endelse: no suited extension found
565   -
566   - // Set enumbins parameter to number of ebounds
567   - m_enumbins = ebds.size();
568   -
569   - } // endif: ebinalg was "FILE"
570   -
571   - else {
572   -
573   - // Initialise log mode for ebinning
574   - bool log = true;
575   -
576   - // check if algorithm is linear
577   - if (m_ebinalg == "LIN") {
578   - log = false;
579   - }
580   -
581   - // todo: should we also check if m_ebinalg is "LOG"
582   - // and throw an exception if neither LIN/LOG/FILE
583   - // is given?
584   -
585   - // Setup energy range covered by data
586   - GEnergy emin(m_emin, "TeV");
587   - GEnergy emax(m_emax, "TeV");
588   - GEbounds ebounds(m_enumbins, emin, emax,log);
589   - ebds = ebounds;
590   -
591   - } //endif: ebinalg was not "FILE"
592   -
593   - // Get Good Time intervals
594   - GGti gti = obs->events()->gti();
595   -
596   - // Get map centre
597   - double xref = m_xref;
598   - double yref = m_yref;
599   - if (m_usepnt) {
600   -
601   - // Get pointer on CTA pointing
602   - const GCTAPointing& pnt = obs->pointing();
603   -
604   - // Set reference point to pointing
605   - if (gammalib::toupper(m_coordsys) == "GAL") {
606   - xref = pnt.dir().l_deg();
607   - yref = pnt.dir().b_deg();
608   - }
609   - else {
610   - xref = pnt.dir().ra_deg();
611   - yref = pnt.dir().dec_deg();
612   - }
613   -
614   - } // endif: used pointing
615   -
616   - // Create skymap
617   - GSkymap map = GSkymap(m_proj, m_coordsys,
618   - xref, yref, -m_binsz, m_binsz,
619   - m_nxpix, m_nypix, m_enumbins);
620   -
621   - // Initialise binning statistics
622   - int num_outside_map = 0;
623   - int num_outside_ebds = 0;
624   - int num_in_map = 0;
625   -
626   - // Fill sky map
627   - GCTAEventList* events =
628   - static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
629   - for (int i = 0; i < events->size(); ++i) {
630   -
631   - // Get event
632   - GCTAEventAtom* event = (*events)[i];
633   -
634   - // Determine sky pixel
635   - GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
636   - GSkyDir dir = inst->dir();
637   - GSkyPixel pixel = map.dir2pix(dir);
638   -
639   - // Skip if pixel is out of range
640   - if (pixel.x() < -0.5 || pixel.x() > (m_nxpix-0.5) ||
641   - pixel.y() < -0.5 || pixel.y() > (m_nypix-0.5)) {
642   - num_outside_map++;
643   - continue;
644   - }
645   -
646   - // Determine energy bin. Skip if we are outside the energy range
647   - int index = ebds.index(event->energy());
648   - if (index == -1) {
649   - num_outside_ebds++;
650   - continue;
651   - }
652   -
653   - // Fill event in skymap
654   - map(pixel, index) += 1.0;
655   - num_in_map++;
656   -
657   - } // endfor: looped over all events
658   -
659   - // Log binning results
660   - if (logTerse()) {
661   - log << std::endl;
662   - log.header1("Binning");
663   - log << gammalib::parformat("Events in list");
664   - log << obs->events()->size() << std::endl;
665   - log << gammalib::parformat("Events in map");
666   - log << num_in_map << std::endl;
667   - log << gammalib::parformat("Events outside map area");
668   - log << num_outside_map << std::endl;
669   - log << gammalib::parformat("Events outside energy bins");
670   - log << num_outside_ebds << std::endl;
671   - }
672   -
673   - // Log map
674   - if (logTerse()) {
675   - log << std::endl;
676   - log.header1("Counts map");
677   - log << map << std::endl;
678   - }
679   -
680   - // Create events cube from sky map
681   - GCTAEventCube cube(map, ebds, gti);
682   -
683   - // Replace event list by event cube in observation
684   - obs->events(cube);
685   -
686   - } // endif: observation was valid
687   -
688   - // Return
689   - return;
690   -}
691   -
692   -
693   -/***********************************************************************//**
694 535 * @brief Initialise counts cube information
695 536 *
696 537 * @exception GException::invalid_value
... ... @@ -764,7 +605,7 @@ void ctbin::init_cube(void)
764 605 *
765 606 * @param[in] obs CTA observation.
766 607 *
767   - * @exception GException::no_list
  608 + * @exception GException::invalid_value
768 609 * No event list found in observation.
769 610 *
770 611 * Fills the events from an event list in the counts cube setup by init_cube.
... ... @@ -776,28 +617,40 @@ void ctbin::fill_cube(GCTAObservation* obs)
776 617  
777 618 // Make sure that the observation holds a CTA event list. If this
778 619 // is not the case then throw an exception.
779   - if (dynamic_cast<const GCTAEventList*>(obs->events()) == NULL) {
780   - throw GException::no_list(G_FILL_CUBE);
  620 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs->events());
  621 + if (events == NULL) {
  622 + std::string msg = "CTA Observation does not contain an event "
  623 + "list. Event list information is needed to "
  624 + "fill the counts map.";
  625 + throw GException::invalid_value(G_FILL_CUBE, msg);
781 626 }
782 627  
  628 + // Get the RoI
  629 + const GCTARoi& roi = events->roi();
  630 +
783 631 // Initialise binning statistics
  632 + int num_outside_roi = 0;
784 633 int num_outside_map = 0;
785 634 int num_outside_ebds = 0;
786 635 int num_in_map = 0;
787 636  
788 637 // Fill sky map
789   - GCTAEventList* events =
790   - static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
791 638 for (int i = 0; i < events->size(); ++i) {
792 639  
793 640 // Get event
794   - GCTAEventAtom* event = (*events)[i];
  641 + const GCTAEventAtom* event = (*events)[i];
795 642  
796 643 // Determine sky pixel
797 644 GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
798 645 GSkyDir dir = inst->dir();
799 646 GSkyPixel pixel = m_cube.dir2pix(dir);
800 647  
  648 + // Skip is pixel is outside RoI
  649 + if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
  650 + num_outside_roi++;
  651 + continue;
  652 + }
  653 +
801 654 // Skip if pixel is out of range
802 655 if (pixel.x() < -0.5 || pixel.x() > (m_nxpix-0.5) ||
803 656 pixel.y() < -0.5 || pixel.y() > (m_nypix-0.5)) {
... ... @@ -831,6 +684,8 @@ void ctbin::fill_cube(GCTAObservation* obs)
831 684 log << obs->events()->size() << std::endl;
832 685 log << gammalib::parformat("Events in cube");
833 686 log << num_in_map << std::endl;
  687 + log << gammalib::parformat("Event bins outside RoI");
  688 + log << num_outside_roi << std::endl;
834 689 log << gammalib::parformat("Events outside cube area");
835 690 log << num_outside_map << std::endl;
836 691 log << gammalib::parformat("Events outside energy bins");
... ... @@ -906,7 +761,7 @@ void ctbin::obs_cube(void)
906 761 container.models(m_obs.models());
907 762  
908 763 // Append CTA observation
909   - m_obs.append(obs);
  764 + container.append(obs);
910 765  
911 766 // Copy over all remaining non-CTA observations
912 767 for (int i = 0; i < m_obs.size(); ++i) {
... ... @@ -926,110 +781,6 @@ void ctbin::obs_cube(void)
926 781 }
927 782  
928 783  
929   -/*==========================================================================
930   - = =
931   - = Private methods =
932   - = =
933   - ==========================================================================*/
934   -
935   -/***********************************************************************//**
936   - * @brief Initialise class members
937   - ***************************************************************************/
938   -void ctbin::init_members(void)
939   -{
940   - // Initialise members
941   - m_evfile.clear();
942   - m_outfile.clear();
943   - m_prefix.clear();
944   - m_proj.clear();
945   - m_coordsys.clear();
946   - m_ebinalg.clear();
947   - m_ebinfile.clear();
948   - m_usepnt = false;
949   - m_emin = 0.0;
950   - m_emax = 0.0;
951   - m_enumbins = 0;
952   - m_xref = 0.0;
953   - m_yref = 0.0;
954   - m_binsz = 0.0;
955   - m_nxpix = 0;
956   - m_nypix = 0;
957   -
958   - // Initialise protected members
959   - m_obs.clear();
960   - m_infiles.clear();
961   - m_use_xml = false;
962   - m_read_ahead = false;
963   - m_cube.clear();
964   - m_ebounds.clear();
965   - m_gti.clear();
966   - m_ontime = 0.0;
967   - m_livetime = 0.0;
968   -
969   - // Set logger properties
970   - log.date(true);
971   -
972   - // Return
973   - return;
974   -}
975   -
976   -
977   -/***********************************************************************//**
978   - * @brief Copy class members
979   - *
980   - * @param[in] app Application.
981   - ***************************************************************************/
982   -void ctbin::copy_members(const ctbin& app)
983   -{
984   - // Copy attributes
985   - m_evfile = app.m_evfile;
986   - m_outfile = app.m_outfile;
987   - m_prefix = app.m_prefix;
988   - m_proj = app.m_proj;
989   - m_coordsys = app.m_coordsys;
990   - m_ebinalg = app.m_ebinalg;
991   - m_ebinfile = app.m_ebinfile;
992   - m_usepnt = app.m_usepnt;
993   - m_emin = app.m_emin;
994   - m_emax = app.m_emax;
995   - m_enumbins = app.m_enumbins;
996   - m_xref = app.m_xref;
997   - m_yref = app.m_yref;
998   - m_binsz = app.m_binsz;
999   - m_nxpix = app.m_nxpix;
1000   - m_nypix = app.m_nypix;
1001   -
1002   - // Copy protected members
1003   - m_obs = app.m_obs;
1004   - m_infiles = app.m_infiles;
1005   - m_use_xml = app.m_use_xml;
1006   - m_read_ahead = app.m_read_ahead;
1007   - m_cube = app.m_cube;
1008   - m_ebounds = app.m_ebounds;
1009   - m_gti = app.m_gti;
1010   - m_ontime = app.m_ontime;
1011   - m_livetime = app.m_livetime;
1012   -
1013   - // Return
1014   - return;
1015   -}
1016   -
1017   -
1018   -/***********************************************************************//**
1019   - * @brief Delete class members
1020   - ***************************************************************************/
1021   -void ctbin::free_members(void)
1022   -{
1023   - // Write separator into logger
1024   - if (logTerse()) {
1025   - log << std::endl;
1026   - }
1027   -
1028   - // Return
1029   - return;
1030   -}
1031   -
1032   -
1033 784 /***********************************************************************//**
1034 785 * @brief Get the energy boundaries
1035 786 *
... ... @@ -1072,7 +823,7 @@ void ctbin::get_ebounds(void)
1072 823 "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
1073 824 " is required if the parameter \"ebinalg\""
1074 825 " is set to \"FILE\".";
1075   - throw GException::invalid_value(G_BIN_EVENTS, msg);
  826 + throw GException::invalid_value(G_GET_EBOUNDS, msg);
1076 827 }
1077 828  
1078 829 // Set enumbins parameter to number of ebounds
... ... @@ -1105,146 +856,3 @@ void ctbin::get_ebounds(void)
1105 856 // Return
1106 857 return;
1107 858 }
1108   -
1109   -
1110   -/***********************************************************************//**
1111   - * @brief Set output file name.
1112   - *
1113   - * @param[in] index Index of input file.
1114   - *
1115   - * Converts an input filename into an output filename by prepending the
1116   - * prefix stored in the member m_prefix to the input filename. Any path will
1117   - * be stripped from the input filename. Also a trailing ".gz" will be
1118   - * stripped. If no input file name exist, prefix and observation
1119   - * number are used.
1120   - ***************************************************************************/
1121   -std::string ctbin::set_outfile_name(const int index) const
1122   -{
1123   -
1124   - // Initialize output file name
1125   - std::string outname="";
1126   -
1127   - // Split input filename into path elements
1128   - std::vector<std::string> elements = gammalib::split(m_infiles[index], "/");
1129   -
1130   - // If there is an observation input file name
1131   - if (elements.size() > 0) {
1132   - outname = m_prefix + elements[elements.size()-1];
1133   - }
1134   - // ...else, use observation number
1135   - else {
1136   - outname = m_prefix + m_obs[index]->id() + ".fits";
1137   - }
1138   -
1139   - // Strip any ".gz"
1140   - outname = gammalib::strip_chars(outname, ".gz");
1141   -
1142   - // Return output filename
1143   - return outname;
1144   -}
1145   -
1146   -
1147   -/***********************************************************************//**
1148   - * @brief Save counts map in FITS format.
1149   - *
1150   - * Save the counts map as a FITS file. The filename of the FITS file is
1151   - * specified by the m_outfile member.
1152   - ***************************************************************************/
1153   -void ctbin::save_fits(void)
1154   -{
1155   - // Get output filename
1156   - m_outfile = (*this)["outfile"].filename();
1157   -
1158   - // Get CTA observation from observation container
1159   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
1160   -
1161   - // Save event list
1162   - save_counts_map(obs, m_outfile);
1163   -
1164   - // Return
1165   - return;
1166   -}
1167   -
1168   -
1169   -/***********************************************************************//**
1170   - * @brief Save counts map(s) in XML format.
1171   - *
1172   - * Save the counts map(s) into FITS files and write the file path information
1173   - * into a XML file. The filename of the XML file is specified by the
1174   - * m_outfile member, the filename(s) of the counts map(s) are built by
1175   - * prepending the prefix given by the m_prefix member to the input counts
1176   - * map(s) filenames. Any path present in the input filename will be stripped,
1177   - * i.e. the counts map(s) will be written in the local working directory
1178   - * (unless a path is specified in the m_prefix member).
1179   - ***************************************************************************/
1180   -void ctbin::save_xml(void)
1181   -{
1182   - // Get output filename and prefix
1183   - m_outfile = (*this)["outfile"].filename();
1184   - m_prefix = (*this)["prefix"].string();
1185   -
1186   - // Issue warning if output filename has no .xml suffix
1187   - std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4));
1188   - if (suffix != ".xml") {
1189   - log << "*** WARNING: Name of observation definition output file \""+
1190   - m_outfile+"\"" << std::endl;
1191   - log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
1192   - log << "*** WARNING: This is not an error, but might be misleading."
1193   - " It is recommended" << std::endl;
1194   - log << "*** WARNING: to use the suffix \".xml\" for observation"
1195   - " definition files." << std::endl;
1196   - }
1197   -
1198   - // Loop over all observation in the container
1199   - for (int i = 0; i < m_obs.size(); ++i) {
1200   -
1201   - // Get CTA observation
1202   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
1203   -
1204   - // Handle only CTA observations
1205   - if (obs != NULL) {
1206   -
1207   - // Set event output file name
1208   - std::string outfile = set_outfile_name(i);
1209   -
1210   - // Store output file name in observation
1211   - obs->eventfile(outfile);
1212   -
1213   - // Save event list
1214   - save_counts_map(obs, outfile);
1215   -
1216   - } // endif: observation was a CTA observations
1217   -
1218   - } // endfor: looped over observations
1219   -
1220   - // Save observations in XML file
1221   - m_obs.save(m_outfile);
1222   -
1223   - // Return
1224   - return;
1225   -}
1226   -
1227   -
1228   -/***********************************************************************//**
1229   - * @brief Save a single counts map into a FITS file
1230   - *
1231   - * @param[in] obs Pointer to CTA observation.
1232   - * @param[in] outfile Output file name.
1233   - *
1234   - * This method saves a single counts map into a FITS file. The method does
1235   - * nothing if the observation pointer is not valid.
1236   - ***************************************************************************/
1237   -void ctbin::save_counts_map(const GCTAObservation* obs,
1238   - const std::string& outfile) const
1239   -{
1240   - // Save only if observation is valid
1241   - if (obs != NULL) {
1242   -
1243   - // Save observation into FITS file
1244   - obs->save(outfile, clobber());
1245   -
1246   - } // endif: observation was valid
1247   -
1248   - // Return
1249   - return;
1250   -}
... ...
src/ctbin/ctbin.hpp
  View file @28b601b
... ... @@ -79,30 +79,21 @@ public:
79 79 void save(void);
80 80 const GObservations& obs(void) const;
81 81 GCTAEventCube cube(void) const;
82   - void get_parameters(void);
83   - void init_cube(void);
84   - void fill_cube(GCTAObservation* obs);
85   - void obs_cube(void);
86 82  
87 83 protected:
88 84 // Protected methods
89   - void init_members(void);
90   - void copy_members(const ctbin& app);
91   - void free_members(void);
92   - void get_ebounds(void);
93   - std::string set_outfile_name(const int index) const;
94   - void save_fits(void);
95   - void save_xml(void);
96   - void save_counts_map(const GCTAObservation* obs,
97   - const std::string& outfile) const;
98   -
99   - // Obsolete methods (to be reomved later)
100   - void bin_events(GCTAObservation* obs);
  85 + void init_members(void);
  86 + void copy_members(const ctbin& app);
  87 + void free_members(void);
  88 + void get_parameters(void);
  89 + void init_cube(void);
  90 + void fill_cube(GCTAObservation* obs);
  91 + void obs_cube(void);
  92 + void get_ebounds(void);
101 93  
102 94 // User parameters
103 95 std::string m_evfile; //!< Input event list or XML file
104 96 std::string m_outfile; //!< Output counts map or XML file
105   - std::string m_prefix; //!< Prefix for multiple counts maps
106 97 bool m_usepnt; //!< Use pointing instead of xref/yref parameters
107 98 double m_emin; //!< Lower energy
108 99 double m_emax; //!< Upper energy
... ... @@ -119,8 +110,6 @@ protected:
119 110  
120 111 // Protected members
121 112 GObservations m_obs; //!< Observation container
122   - std::vector<std::string> m_infiles; //!< Input event filenames
123   - bool m_use_xml; //!< Use XML file instead of FITS file
124 113 bool m_read_ahead; //!< Read ahead parameters
125 114 GSkymap m_cube; //!< Event cube
126 115 GEbounds m_ebounds; //!< Energy boundaries
... ...
src/ctbin/ctbin.par
  View file @28b601b
... ... @@ -19,7 +19,6 @@
19 19 #===================
20 20 evfile, f, a, "events.fits",,,"Input event list or observation definition file"
21 21 outfile, f, a, "cntmap.fits",,, "Output counts cube"
22   -prefix, s, h, "cntmap_",,,"Prefix (no longer used)"
23 22 #
24 23 # Energy binning parameters
25 24 #==========================
... ...