Karim Sabri / ctools

Commit fc5614e8a7710afc1af6e6000c1c19e6087548c5

Authored by Jürgen Knödlseder
1 parent 1a11125e

Integrate ctprob tool

Inline Side-by-side
Showing 12 changed files with 363 additions and 660 deletions
doc/source/users/reference_manual/ctprob.rst
  View file @fc5614e
1   - .. _ctprob:
  1 +.. _ctprob:
2 2  
3 3 ctprob
4 4 ========
5 5  
6   -Computes event probability for a given model.
  6 +Compute event probability for a given model.
7 7  
8 8  
9 9 Synopsis
10 10 --------
11 11  
12   -This tool takes an event list file (or an XML observation file) and a source
13   -model as input. For each event, it evalutes the probability for that event
14   -to come from each of the source in the model. This probability is obtained
15   -evaluating the differential expected counts for each source at the event
16   -direction and energy and normalizing it to the total differential expected
17   -counts for the given model. If only one source is provided in the model,
18   -the probability is set to 1 for each photon.
19   -A column for each source in the model is added to the event list in order to
20   -store these values. The name of each column is made by the source name with
21   -the prefix "PROB_".
  12 +This tool takes an event list file (or an XML observation definition file) and
  13 +a model as input and computes for each event the probability that it originates
  14 +from a specific model component. The probability is computed by by evaluating
  15 +the differential event probability for a given model and normalizing these
  16 +probabilities so that the sum of the probabilities for all model components is
  17 +unity for each event.
  18 +
  19 +The tool appends for each model component a single precision column to the event
  20 +list that contains the event probabilities. The name of the column is build from
  21 +the model name prefixed with ``PROB_``.
22 22  
23 23  
24 24 General parameters
... ... @@ -34,7 +34,7 @@ General parameters
34 34 Prefix for output event lists in observation definition file.
35 35  
36 36 ``inmodel [string]``
37   - Input model XML file.
  37 + Input model definition XML file.
38 38  
39 39 ``caldb [string]``
40 40 Calibration database.
... ... @@ -43,7 +43,7 @@ General parameters
43 43 Instrument response function.
44 44  
45 45 ``(edisp = no) [boolean]``
46   - Apply energy dispersion to response computation.
  46 + Applies energy dispersion to response computation.
47 47  
48 48  
49 49 Standard parameters
... ... @@ -73,7 +73,7 @@ Standard parameters
73 73 ``(mode = ql) [string]``
74 74 Mode of automatic parameters (default is "ql", i.e. "query and learn").
75 75  
76   -``(logfile = ctlike.log) [string]``
  76 +``(logfile = ctprob.log) [string]``
77 77 Name of log file.
78 78  
79 79  
... ...
doc/source/users/reference_manual/index.rst
  View file @fc5614e
... ... @@ -27,11 +27,11 @@ ctools
27 27 ctmodel --- Computes model counts cube <ctmodel>
28 28 ctobssim --- Simulate observations <ctobssim>
29 29 ctpsfcube --- Generates point spread function cube <ctpsfcube>
  30 + ctprob --- Computes event probability for a given model <ctprob>
30 31 ctselect --- Selects event data <ctselect>
31 32 ctskymap --- Generates sky map <ctskymap>
32 33 cttsmap --- Generates Test Statistic map <cttsmap>
33 34 ctulimit --- Calculates upper limit <ctulimit>
34   - ctprob --- Computes event probability for a given model <ctprob>
35 35  
36 36  
37 37 cscripts
... ...
src/ctprob/Makefile.am
  View file @fc5614e
1 1 #############################################################################
2   -# Makefile.am for ctprob tool #
  2 +# Makefile.am for ctprob tool #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016 Juergen Knoedlseder #
  4 +# Copyright (C) 2017 Leonardo Di Venere #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ...
src/ctprob/ctprob.cpp
  View file @fc5614e
1 1 /***************************************************************************
2   - * ctprob - Computes probability for a given model *
  2 + * ctprob - Computes event probability for a given model *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Leonardo Di Venere *
  4 + * copyright (C) 2017 by Leonardo Di Venere *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctprob.cpp
23   - * @brief Computes probability for a given model
  23 + * @brief Computes event probability for a given model
24 24 * @author Leonardo Di Venere
25 25 */
26 26  
... ... @@ -33,13 +33,12 @@
33 33 #include "GTools.hpp"
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
36   -#define G_RUN "ctprob::run()"
37   -#define G_GET_OBS "ctprob::get_obs()"
  36 +#define G_GET_OBS "ctprob::get_obs()"
  37 +#define G_EVALUATE_PROB "ctprob::evaluate_probability(GCTAObservation*)"
38 38  
39 39 /* __ Debug definitions __________________________________________________ */
40 40  
41 41 /* __ Coding definitions _________________________________________________ */
42   -//#define G_USE_MKSTEMP //!< Use mkstemp for temporary filename
43 42  
44 43  
45 44 /*==========================================================================
... ... @@ -51,7 +50,7 @@
51 50 /***********************************************************************//**
52 51 * @brief Void constructor
53 52 ***************************************************************************/
54   -ctprob::ctprob(void) : ctobservation(CTPROB_NAME, CTPROB_VERSION)
  53 +ctprob::ctprob(void) : ctobservation(CTPROB_NAME, VERSION)
55 54 {
56 55 // Initialise members
57 56 init_members();
... ... @@ -70,7 +69,7 @@ ctprob::ctprob(void) : ctobservation(CTPROB_NAME, CTPROB_VERSION)
70 69 * provided in an observation container.
71 70 ***************************************************************************/
72 71 ctprob::ctprob(const GObservations& obs) :
73   - ctobservation(CTPROB_NAME, CTPROB_VERSION, obs)
  72 + ctobservation(CTPROB_NAME, VERSION, obs)
74 73 {
75 74 // Initialise members
76 75 init_members();
... ... @@ -89,7 +88,7 @@ ctprob::ctprob(const GObservations&amp; obs) :
89 88 * @param[in] argv Array of command line arguments.
90 89 ***************************************************************************/
91 90 ctprob::ctprob(int argc, char *argv[]) :
92   - ctobservation(CTPROB_NAME, CTPROB_VERSION, argc, argv)
  91 + ctobservation(CTPROB_NAME, VERSION, argc, argv)
93 92 {
94 93 // Initialise members
95 94 init_members();
... ... @@ -201,13 +200,11 @@ void ctprob::clear(void)
201 200  
202 201  
203 202 /***********************************************************************//**
204   - * @brief Append probability columns event data
  203 + * @brief Runs the ctprob tool
205 204 *
206 205 * This method reads in the application parameters and loops over all
207   - * observations that were found to perform the probability calculation.
208   - * Each observation is written to a temporary file, which is then
209   - * re-opened and used for the actual calculation. The temporary file
210   - * is deleted after this action so that no disk overflow will occur.
  206 + * unbinned observations to compute the probability for each event that it
  207 + * comes from one of the model components.
211 208 ***************************************************************************/
212 209 void ctprob::run(void)
213 210 {
... ... @@ -219,123 +216,30 @@ void ctprob::run(void)
219 216 // Get parameters
220 217 get_parameters();
221 218  
  219 + // Set energy dispersion flags of all CTA observations and save old
  220 + // values in save_edisp vector
  221 + std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
  222 +
222 223 // Write input observation container into logger
223 224 log_observations(NORMAL, m_obs, "Input observation");
224 225  
225 226 // Write header into logger
226   - log_header1(TERSE, "Event probability evaluation");
  227 + log_header1(TERSE, "Compute event probabilities");
227 228  
228 229 // Initialise counters
229 230 int n_observations = 0;
230 231  
231   - // Loop over all observation in the container
232   - for (int i = 0; i < m_obs.size(); ++i) {
233   -
234   - // Write header for the current observation
235   - log_header3(TERSE, get_obs_header(m_obs[i]));
236   -
237   - // Initialise event input and output filenames and the event
238   - // and GTI extension names
239   - m_infiles.push_back("");
240   - m_evtname.push_back(gammalib::extname_cta_events);
241   - m_gtiname.push_back(gammalib::extname_gti);
242   -
243   - // Get CTA observation
244   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
245   -
246   - // Skip observation if it's not CTA
247   - if (obs == NULL) {
248   - std::string msg = " Skipping "+m_obs[i]->instrument()+
249   - " observation";
250   - log_string(NORMAL, msg);
251   - continue;
252   - }
253   -
254   - // Skip observation if we have a binned observation
255   - if (obs->eventtype() == "CountsCube") {
256   - std::string msg = " Skipping binned "+m_obs[i]->instrument()+
257   - " observation";
258   - log_string(NORMAL, msg);
259   - continue;
260   - }
  232 + // Loop over all unbinned CTA observations in the container
  233 + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
  234 + obs = next_unbinned_observation()) {
261 235  
262 236 // Increment counter
263 237 n_observations++;
264 238  
265   - // Save event file name (for possible saving)
266   - m_infiles[i] = obs->eventfile();
267   -
268   - // Extract event and GTI extension names from input FITS file
269   - GFilename fname(m_infiles[i]);
270   - if (fname.has_extname()) {
271   - m_evtname[i] = fname.extname();
272   - }
273   - m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]);
274   -
275   - // Write input file information into logger
276   - log_value(NORMAL, "Input filename", m_infiles[i]);
277   - log_value(NORMAL, "Event extension name", m_evtname[i]);
278   - log_value(NORMAL, "GTI extension name", m_gtiname[i]);
279   -
280   - // Fall through in case that the event file is empty
281   - if (obs->events()->size() == 0) {
282   - log_string(NORMAL, " Warning: No events in event file \""+
283   - m_infiles[i]+"\". Event probability evaluation skipped.");
284   - continue;
285   - }
286   -
287   - // If we have an input file then check it
288   - if (!m_infiles[i].empty()) {
289   - std::string message = check_infile(m_infiles[i], m_evtname[i]);
290   - if (!message.empty()) {
291   - throw GException::invalid_value(G_RUN, message);
292   - }
293   - }
294   -
295   - // Get temporary file name
296   - #if G_USE_MKSTEMP
297   - char tpl[] = "ctprobXXXXXX";
298   - int fileid = mkstemp(tpl);
299   - std::string filename(tpl);
300   - #else
301   - std::string filename = std::tmpnam(NULL);
302   - #endif
303   -
304   - // Save observation in temporary file. We add here the events and
305   - // GTI extension name so that the GCTAObservation::save method can
306   - // use this information for writing the proper extension names into
307   - // the temporary file
308   - obs->save(filename+"["+m_evtname[i]+";"+ m_gtiname[i]+"]", true);
309   -
310   - // Log saved FITS file.
311   - if (logExplicit()) {
312   - GFits tmpfile(filename);
313   - log.header3("FITS file content of temporary file");
314   - log << tmpfile << std::endl;
315   - tmpfile.close();
316   - }
317   -
318   - // If we have a temporary file then check it
319   - if (!filename.empty()) {
320   - std::string message = check_infile(filename, m_evtname[i]);
321   - if (!message.empty()) {
322   - throw GException::invalid_value(G_RUN, message);
323   - }
324   - }
325   -
326   - // Load observation from temporary file, including event selection
  239 + // Evaluate event probabilities
327 240 evaluate_probability(obs);
328   - //select_events(obs, filename, m_evtname[i], m_gtiname[i]);
329   -
330   - // Close temporary file
331   - #if G_USE_MKSTEMP
332   - close(fileid);
333   - #endif
334 241  
335   - // Remove temporary file
336   - std::remove(filename.c_str());
337   -
338   - } // endfor: looped over all observations
  242 + } // endfor: looped over observations
339 243  
340 244 // If more than a single observation has been handled then make sure that
341 245 // an XML file will be used for storage
... ... @@ -346,6 +250,9 @@ void ctprob::run(void)
346 250 // Write observation(s) into logger
347 251 log_observations(NORMAL, m_obs, "Output observation");
348 252  
  253 + // Restore energy dispersion flags of all CTA observations
  254 + restore_edisp(m_obs, save_edisp);
  255 +
349 256 // Optionally publish event list(s)
350 257 if ((*this)["publish"].boolean()) {
351 258 publish();
... ... @@ -379,6 +286,7 @@ void ctprob::save(void)
379 286 {
380 287 // Write header into logger
381 288 log_header1(TERSE, gammalib::number("Save event list", m_obs.size()));
  289 +
382 290 // Case A: Save event file(s) and XML metadata information
383 291 if (m_use_xml) {
384 292 save_xml();
... ... @@ -544,11 +452,6 @@ void ctprob::free_members(void)
544 452 ***************************************************************************/
545 453 void ctprob::get_parameters(void)
546 454 {
547   - // Initialise selection flags
548   - /* m_select_energy = true;
549   - m_select_roi = true;
550   - m_select_time = true;
551   - */
552 455 // Setup observations from "inobs" parameter. Do not request response
553 456 // information and do not accept counts cubes.
554 457 if (m_obs.size() == 0) {
... ... @@ -559,7 +462,6 @@ void ctprob::get_parameters(void)
559 462 else {
560 463 setup_observations(m_obs);
561 464 }
562   - //setup_observations(m_obs, false, true, false);
563 465  
564 466 // Read model definition file if required
565 467 if (m_obs.models().size() == 0) {
... ... @@ -671,144 +573,90 @@ void ctprob::get_obs(void)
671 573 *
672 574 * @param[in,out] obs CTA observation.
673 575 *
  576 + * @exception GException::invalid_value
  577 + * Observation does not contain an event list.
674 578 *
675   - * For each event from the input event file evaluates the probability that
676   - * the event comes from any of the source in the input model. This is done
677   - * by evaluating the differential expected counts for each source at the event
678   - * direction and energy and normalizing it to the total differential expected
679   - * counts for the given model. If only one source is provided in the model,
680   - * the probability is set to 1 for each photon.
  579 + * Computes for all events the probability that an event arises from a
  580 + * specific model component. The computation is done by by evaluating the
  581 + * differential event probability for a given model and normalizing these
  582 + * probabilities so that the sum of the probabilities for all model
  583 + * components is unity for each event.
681 584 *
682   - * A FITS column for each source in the model is created and added to
683   - * the event list. The name of each column is made by the source name with
684   - * the prefix "PROB_".
  585 + * The method appends for each model component a single precision column to
  586 + * the event list that contains the event probabilities. The name of the
  587 + * column is build from the model name prefixed with `PROB_`.
685 588 ***************************************************************************/
686   -void ctprob::evaluate_probability(GCTAObservation* obs)
  589 +void ctprob::evaluate_probability(GCTAObservation* obs)
687 590 {
688 591 // Write header into logger
689   - log_header3(NORMAL, "Evaluating probability for events");
690   -
691   - GCTAEventList* evts = static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
692   - const GEvent* evt;
693   -
694   - double total=0.;
695   - double value;
696   -
697   - std::vector<GFitsTableFloatCol*> columns;
698   -
699   - //define colums
700   - for (int j = 0; j<m_obs.models().size(); ++j) {
701   - const std::string mdl_name = m_obs.models()[j] ->name();
702   - GFitsTableFloatCol* col = new GFitsTableFloatCol("PROB_"+mdl_name, evts->size());
703   - columns.push_back(col);
704   - }
705   -
706   - // Loop over events
707   - for (int i = 0; i < evts->size(); ++i) {
708   - total = 0.;
709   - evt = (*evts)[i];
710   - std::vector<double> values;
711   - //Check if you have more than one model
712   - if (m_obs.models().size()>1) {
713   - //Loop over models
714   - for (int j = 0; j<m_obs.models().size(); ++j) {
715   - value = m_obs.models()[j]->eval( *evt, *obs);
716   - values.push_back(value);
717   - total += value;
718   - }
719   - }
720   - // if there is only one model a constant value equal to 1. is set for each photon
721   - else {
722   - value = 1.;
723   - values.push_back(value);
724   - total += value;
725   - } // end if
726   - for (int j = 0; j<m_obs.models().size(); ++j) {
727   - values[j] /= total;
728   - (*(columns[j]))(i) = values[j];
729   - }
730   - }
731   -
732   - //append columns to observation event list
733   - for (int j = 0; j<m_obs.models().size(); ++j) {
734   - evts->append_column(*columns[j]);
735   - }
736   -
737   - // Return
738   - return;
739   -}
740   -
741   -
742   -
743   -/***********************************************************************//**
744   - * @brief Check input filename
745   - *
746   - * @param[in] filename File name.
747   - * @param[in] evtname Event extension name.
748   - *
749   - * This method checks if the input FITS file is correct.
750   - ***************************************************************************/
751   -std::string ctprob::check_infile(const std::string& filename,
752   - const std::string& evtname) const
753   -{
754   - // Initialise message string
755   - std::string message = "";
756   -
757   - // Open FITS file
758   - GFits fits(filename);
759   -
760   - // Check for existence of events extensions
761   - if (!fits.contains(evtname)) {
762   - message = "No \""+evtname+"\" extension found in input file \""+
763   - filename + "\".";
  592 + log_header3(NORMAL, "Evaluate events probabilities");
  593 +
  594 + // Make sure that the observation holds a CTA event list. If this
  595 + // is not the case then throw an exception.
  596 + GCTAEventList* events = dynamic_cast<GCTAEventList*>(const_cast<GEvents*>
  597 + (obs->events()));
  598 + if (events == NULL) {
  599 + std::string msg = "CTA Observation does not contain an event list. An "
  600 + "event list is needed to compute the event "
  601 + "probabilities.";
  602 + throw GException::invalid_value(G_EVALUATE_PROB, msg);
764 603 }
765 604  
766   - // ... otherwise check column names
767   - else {
768   -
769   - // Get pointer to FITS table
770   - GFitsTable* table = fits.table(evtname);
  605 + // Continue only if there are model components
  606 + int nmodels = m_obs.models().size();
  607 + if (nmodels > 0) {
  608 +
  609 + // Allocate single precision FITS table columns for all model components
  610 + std::vector<GFitsTableFloatCol*> columns;
  611 + for (int k = 0; k < nmodels; ++k) {
  612 + const std::string colname = "PROB_" + m_obs.models()[k] ->name();
  613 + GFitsTableFloatCol* col =
  614 + new GFitsTableFloatCol(colname, events->size());
  615 + columns.push_back(col);
  616 + }
771 617  
772   - // Initialise list of missing columns
773   - std::vector<std::string> missing;
  618 + // Loop over all events in the event list
  619 + for (int i = 0; i < events->size(); ++i) {
774 620  
775   - // Check for existence of TIME column
776   - if (!table->contains("TIME")) {
777   - missing.push_back("TIME");
778   - }
  621 + // Initialise total probability
  622 + double total = 0.0;
779 623  
780   - // Check for existence of ENERGY column
781   - if (!table->contains("ENERGY")) {
782   - missing.push_back("ENERGY");
783   - }
  624 + // Get point to event
  625 + const GEvent* event = (*events)[i];
784 626  
785   - // Check for existence of RA column
786   - if (!table->contains("RA")) {
787   - missing.push_back("RA");
788   - }
  627 + // Initialise probability vector
  628 + std::vector<double> probabilities(nmodels, 0.0);
789 629  
790   - // Check for existence of DEC column
791   - if (!table->contains("DEC")) {
792   - missing.push_back("DEC");
793   - }
  630 + // Loop over models
  631 + for (int k = 0; k < nmodels; ++k) {
  632 + double value = m_obs.models()[k]->eval(*event, *obs);
  633 + probabilities[k] = value;
  634 + total += value;
  635 + }
794 636  
795   - // Set error message for missing columns
796   - if (!missing.empty()) {
797   - message = "The following columns are missing in the "
798   - "\""+evtname+"\" extension of input file \""+
799   - filename + "\": ";
800   - for (int i = 0; i < missing.size(); ++i) {
801   - message += "\"" + missing[i] + "\"";
802   - if (i < missing.size()-1) {
803   - message += ", ";
  637 + // Normalize probabilities to unity
  638 + if (total > 0.0) {
  639 + for (int k = 0; k < nmodels; ++k) {
  640 + probabilities[k] /= total;
804 641 }
805 642 }
  643 +
  644 + // Store probabilities in FITS columns
  645 + for (int k = 0; k < nmodels; ++k) {
  646 + (*(columns[k]))(i) = probabilities[k];
  647 + }
  648 +
  649 + } // endfor: looped over all events
  650 +
  651 + // Append columns to event list
  652 + for (int k = 0; k < nmodels; ++k) {
  653 + events->append_column(*columns[k]);
806 654 }
807 655  
808   - } // endelse: checked column names
  656 + } // endif: there were model components
809 657  
810 658 // Return
811   - return message;
  659 + return;
812 660 }
813 661  
814 662  
... ... @@ -842,43 +690,10 @@ std::string ctprob::set_outfile_name(const std::string&amp; filename) const
842 690  
843 691  
844 692 /***********************************************************************//**
845   - * @brief Get Good Time Intervals extension name
846   - *
847   - * @param[in] filename Input file name.
848   - * @param[in] evtname Events extension name.
849   - *
850   - * Extracts the Good Time Intervals extension name from the event file. We
851   - * do this by loading the events and accessing the Good Time Intervals
852   - * extension name using the GCTAEventList::gtiname() method. If the file name
853   - * is empty, the method returns `GTI`.
854   - ***************************************************************************/
855   -std::string ctprob::get_gtiname(const std::string& filename,
856   - const std::string& evtname) const
857   -{
858   - // Initialise GTI name
859   - std::string gtiname = gammalib::extname_gti;
860   -
861   - // Continue only if the filename is not empty
862   - if (!filename.empty()) {
863   -
864   - // Load events
865   - GCTAEventList events(filename+"["+evtname+"]");
866   -
867   - // Get GTI name
868   - gtiname = events.gtiname();
869   -
870   - }
871   -
872   - // Return GTI name
873   - return (gtiname);
874   -}
875   -
876   -
877   -/***********************************************************************//**
878 693 * @brief Save event list in FITS format.
879 694 *
880 695 * Save the event list as a FITS file. The file name of the FITS file is
881   - * specified by the "outobs" parameter.
  696 + * specified by the `outobs` parameter.
882 697 ***************************************************************************/
883 698 void ctprob::save_fits(void)
884 699 {
... ... @@ -888,39 +703,42 @@ void ctprob::save_fits(void)
888 703 // Get output filename
889 704 m_outobs = (*this)["outobs"].filename();
890 705  
891   - // Get CTA observation from observation container
892   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
  706 + // Loop over all unbinned CTA observations in the container
  707 + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
  708 + obs = next_unbinned_observation()) {
893 709  
894   - // Save only if it's a CTA observation
895   - if (obs != NULL) {
896   -
897   - // Save only if file name is non-empty
898   - if (m_infiles[0].length() > 0) {
  710 + // Get input file name and default extension names
  711 + std::string infile = obs->eventfile();
  712 + std::string evtname = gammalib::extname_cta_events;
  713 + std::string gtiname = gammalib::extname_gti;
899 714  
900   - // Create file name object
901   - GFilename fname(m_outobs);
  715 + // Create file name object
  716 + GFilename fname(m_outobs);
902 717  
903   - // Extract filename and event extension name
904   - std::string outfile = fname.url();
  718 + // Extract filename and event extension name
  719 + std::string outfile = fname.url();
905 720  
906   - // Append event extension name. We handle here the possibility
907   - // to write the events into a different extension.
908   - if (fname.has_extname()) {
909   - outfile += "["+fname.extname()+"]";
910   - }
911   - else {
912   - outfile += "["+m_evtname[0]+"]";
913   - }
914   -
915   - // Log filename
916   - log_value(NORMAL, "Event list file", outfile);
917   - // Save event list
918   - save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0],
919   - outfile);
  721 + // Append event extension name. We handle here the possibility
  722 + // to write the events into a different extension.
  723 + if (fname.has_extname()) {
  724 + outfile += "["+fname.extname()+"]";
  725 + }
  726 + else {
  727 + outfile += "["+gammalib::extname_cta_events+"]";
  728 + }
  729 +
  730 + // Log filename
  731 + log_value(NORMAL, "Event list file", outfile);
  732 +
  733 + // Save event list
  734 + save_event_list(obs, infile, evtname, gtiname, outfile);
920 735  
921   - } // endif: filename was non empty
  736 + // Exit the loop (if this method is called we should only
  737 + // have a single unbinned observation in the observation
  738 + // container
  739 + break;
922 740  
923   - } // endif: observation was CTA observation
  741 + } // endfor: looped over all unbinned observations
924 742  
925 743 } // endif: there were observations
926 744  
... ... @@ -933,7 +751,7 @@ void ctprob::save_fits(void)
933 751 * @brief Save event list(s) in XML format.
934 752 *
935 753 * Save the event list(s) into FITS files and write the file path information
936   - * into a XML file. The filename of the XML file is specified by the outfile
  754 + * into a XML file. The filename of the XML file is specified by the `outobs`
937 755 * parameter, the filename(s) of the event lists are built by prepending a
938 756 * prefix to the input event list filenames. Any path present in the input
939 757 * filename will be stripped, i.e. the event list(s) will be written in the
... ... @@ -948,27 +766,27 @@ void ctprob::save_xml(void)
948 766 // Issue warning if output filename has no .xml suffix
949 767 log_string(TERSE, warn_xml_suffix(m_outobs));
950 768  
951   - // Loop over all observation in the container
952   - for (int i = 0; i < m_obs.size(); ++i) {
953   -
954   - // Get CTA observation
955   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  769 + // Loop over all unbinned CTA observations in the container
  770 + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
  771 + obs = next_unbinned_observation()) {
956 772  
957   - // Skip observations that are no CTA observations
958   - if (obs == NULL) {
959   - continue;
960   - }
  773 + // Get input file name and default extension names
  774 + std::string infile = obs->eventfile();
  775 + std::string evtname = gammalib::extname_cta_events;
  776 + std::string gtiname = gammalib::extname_gti;
961 777  
962   - // Skip observations that have empty names
963   - if (m_infiles[i].length() == 0) {
964   - continue;
  778 + // Extract event and GTI extension names from input FITS file
  779 + GFilename fname(infile);
  780 + if (fname.has_extname()) {
  781 + evtname = fname.extname();
965 782 }
  783 + gtiname = get_gtiname(fname.url(), evtname);
966 784  
967 785 // Set event output file name
968   - std::string outfile = set_outfile_name(m_infiles[i]);
  786 + std::string outfile = set_outfile_name(infile);
969 787  
970 788 // Append event extension name
971   - outfile += "["+m_evtname[i]+"]";
  789 + outfile += "["+evtname+"]";
972 790  
973 791 // Log filename
974 792 log_value(NORMAL, "Event list file", outfile);
... ... @@ -977,8 +795,7 @@ void ctprob::save_xml(void)
977 795 obs->eventfile(outfile);
978 796  
979 797 // Save event list
980   - save_event_list(obs, m_infiles[i], m_evtname[i], m_gtiname[i],
981   - outfile);
  798 + save_event_list(obs, infile, evtname, gtiname, outfile);
982 799  
983 800 } // endfor: looped over observations
984 801  
... ... @@ -991,109 +808,3 @@ void ctprob::save_xml(void)
991 808 // Return
992 809 return;
993 810 }
994   -
995   -
996   -/***********************************************************************//**
997   - * @brief Save event list into FITS file
998   - *
999   - * @param[in] obs Pointer to CTA observation.
1000   - * @param[in] infile Input file name.
1001   - * @param[in] evtname Event extension name.
1002   - * @param[in] gtiname GTI extension name.
1003   - * @param[in] outfile Output file name.
1004   - *
1005   - * Saves an event list and the corresponding Good Time Intervals into a FITS
1006   - * file and copy all others extensions from the input file to the output
1007   - * file.
1008   - *
1009   - * If an extension name is specified in the @p outfile argument, the events
1010   - * and eventually also the Good Time Intervals will be extracted from the
1011   - * argument and used for writing the events. The format is
1012   - *
1013   - * <filename>[<event extension name;GTI extension name>]
1014   - *
1015   - * where <filename> needs to be replaced by the name of the FITS file,
1016   - * and <event extension name;GTI extension name> by the name of the events
1017   - * and Good Time Intervals extensions. For example
1018   - *
1019   - * myfits.fits[EVENTS1;GTI1]
1020   - *
1021   - * will write the selected events into the "EVENTS1" extension and the
1022   - * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS
1023   - * file. If the Good Time Intervals extension name is skipped, e.g.
1024   - *
1025   - * myfits.fits[EVENTS1]
1026   - *
1027   - * the original extension name for the Good Time Intervals will be kept.
1028   - * Analogously, only the Good Time Intervals extension name can be changed
1029   - * by specifying
1030   - *
1031   - * myfits.fits[;GTI1]
1032   - *
1033   - * In none of the cases will the original events and Good Time Intervals be
1034   - * copied over to the output file.
1035   - ***************************************************************************/
1036   -void ctprob::save_event_list(const GCTAObservation* obs,
1037   - const std::string& infile,
1038   - const std::string& evtname,
1039   - const std::string& gtiname,
1040   - const std::string& outfile) const
1041   -{
1042   - // Save only if we have an event list
1043   - if (obs->eventtype() == "EventList") {
1044   -
1045   - // Set output FITS file event extension names
1046   - GFilename outname(outfile);
1047   - std::string outevt = evtname;
1048   - std::string outgti = gtiname;
1049   - if (outname.has_extname()) {
1050   - std::vector<std::string> extnames =
1051   - gammalib::split(outname.extname(), ";");
1052   - if (extnames.size() > 0) {
1053   - std::string extname = gammalib::strip_whitespace(extnames[0]);
1054   - if (!extname.empty()) {
1055   - outevt = extname;
1056   - }
1057   - }
1058   - if (extnames.size() > 1) {
1059   - std::string extname = gammalib::strip_whitespace(extnames[1]);
1060   - if (!extname.empty()) {
1061   - outgti = extname;
1062   - }
1063   - }
1064   - }
1065   -
1066   - // Create output FITS file
1067   - GFits outfits;
1068   -
1069   - // Write observation into FITS file
1070   - obs->write(outfits, outevt, outgti);
1071   -
1072   - // Copy all extensions other than evtname and gtiname extensions
1073   - // from the input to the output event list. The evtname and
1074   - // gtiname extensions are written by the save method, all others
1075   - // that may eventually be present have to be copied over
1076   - // explicitly.
1077   - GFits infits(infile);
1078   - for (int extno = 1; extno < infits.size(); ++extno) {
1079   - GFitsHDU* hdu = infits.at(extno);
1080   - if (hdu->extname() != evtname &&
1081   - hdu->extname() != gtiname &&
1082   - hdu->extname() != outevt &&
1083   - hdu->extname() != outgti) {
1084   - outfits.append(*hdu);
1085   - }
1086   - }
1087   -
1088   - // Close input file
1089   - infits.close();
1090   -
1091   - // Save file to disk and close it (we need both operations)
1092   - outfits.saveto(outname.url(), clobber());
1093   - outfits.close();
1094   -
1095   - } // endif: observation was unbinned
1096   -
1097   - // Return
1098   - return;
1099   -}
... ...
src/ctprob/ctprob.hpp
  View file @fc5614e
1 1 /***************************************************************************
2   - * ctprob - Computes probability for a given model *
  2 + * ctprob - Computes event probability for a given model *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Leonardo Di Venere *
  4 + * copyright (C) 2017 by Leonardo Di Venere *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctprob.hpp
23   - * @brief Computes probability for a given model
  23 + * @brief Computes event probability for a given model
24 24 * @author Leonardo Di Venere
25 25 */
26 26  
... ... @@ -35,8 +35,7 @@
35 35 #include "ctobservation.hpp"
36 36  
37 37 /* __Definitions _________________________________________________________ */
38   -#define CTPROB_NAME "ctprob"
39   -#define CTPROB_VERSION "0.0.1"
  38 +#define CTPROB_NAME "ctprob"
40 39  
41 40  
42 41 /***********************************************************************//**
... ... @@ -70,32 +69,22 @@ protected:
70 69 void free_members(void);
71 70 void get_parameters(void);
72 71 void get_obs(void);
73   - void evaluate_probability(GCTAObservation* obs);
74   - std::string check_infile(const std::string& filename,
75   - const std::string& evtname) const;
  72 + void evaluate_probability(GCTAObservation* obs);
76 73 std::string set_outfile_name(const std::string& filename) const;
77   - std::string get_gtiname(const std::string& filename,
78   - const std::string& evtname) const;
79 74 void save_fits(void);
80 75 void save_xml(void);
81   - void save_event_list(const GCTAObservation* obs,
82   - const std::string& infile,
83   - const std::string& evtname,
84   - const std::string& gtiname,
85   - const std::string& outfile) const;
86 76  
87 77 // User parameters
88   - std::string m_outobs; //!< Output event list or XML file
89   - std::string m_prefix; //!< Prefix for multiple event lists
90   - bool m_apply_edisp; //!< Apply energy dispersion?
91   - bool m_publish; //!< Publish model cube?
92   - GChatter m_chatter; //!< Chattiness
  78 + std::string m_outobs; //!< Output event list or XML file
  79 + std::string m_prefix; //!< Prefix for multiple event lists
  80 + bool m_apply_edisp; //!< Apply energy dispersion?
  81 + bool m_publish; //!< Publish event list?
  82 + GChatter m_chatter; //!< Chattiness
93 83  
94 84 // Protected members
95   - std::vector<std::string> m_infiles; //!< Input event filenames
96   - std::vector<std::string> m_evtname; //!< Event extension names
97   - std::vector<std::string> m_gtiname; //!< GTI extension names
  85 + std::vector<std::string> m_infiles; //!< Input event filenames
  86 + std::vector<std::string> m_evtname; //!< Event extension names
  87 + std::vector<std::string> m_gtiname; //!< GTI extension names
98 88 };
99 89  
100   -
101 90 #endif /* CTPROB_HPP */
... ...
src/ctprob/ctprob.par
  View file @fc5614e
... ... @@ -8,7 +8,7 @@
8 8 # #
9 9 # Author: Leonardo Di Venere #
10 10 # leonardo.divenere@ba.infn.it #
11   -# INFN #
  11 +# INFN #
12 12 # #
13 13 # Purpose: Parameter file for the ctprob application #
14 14 # #
... ... @@ -16,13 +16,13 @@
16 16 #
17 17 # General parameters
18 18 #===================
19   -inobs, f, a, "events.fits",,, "Input event list or observation definition XML file"
20   -outobs, f, a, "prob_events.fits",,, "Output event list or observation definition XML file"
21   -prefix, s, h, "prob_",,, "Prefix for event lists in observation definition XML file"
22   -inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file"
23   -caldb, s, a, "prod2",,, "Calibration database"
24   -irf, s, a, "South_0.5h",,, "Instrument response function"
25   -edisp, b, h, no,,, "Apply energy dispersion?"
  19 +inobs, f, a, "events.fits",,, "Input event list or observation definition XML file"
  20 +outobs, f, a, "prob_events.fits",,, "Output event list or observation definition XML file"
  21 +prefix, s, h, "prob_",,, "Prefix for event lists in observation definition XML file"
  22 +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file"
  23 +caldb, s, a, "prod2",,, "Calibration database"
  24 +irf, s, a, "South_0.5h",,, "Instrument response function"
  25 +edisp, b, h, no,,, "Apply energy dispersion?"
26 26  
27 27 #
28 28 # Standard parameters
... ...
src/ctprob/main.cpp
  View file @fc5614e
1 1 /***************************************************************************
2   - * ctprob - Computes probability for a given model *
  2 + * ctprob - Computes event probability for a given model *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Leonardo Di Venere *
  4 + * copyright (C) 2017 by Leonardo Di Venere *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file main.cpp
23   - * @brief Computes probability for a given model
  23 + * @brief Computes event probability for a given model
24 24 * @author Leonardo Di Venere
25 25 */
26 26  
... ...
src/ctselect/ctselect.cpp
  View file @fc5614e
... ... @@ -1175,39 +1175,6 @@ std::string ctselect::set_outfile_name(const std::string&amp; filename) const
1175 1175  
1176 1176  
1177 1177 /***********************************************************************//**
1178   - * @brief Get Good Time Intervals extension name
1179   - *
1180   - * @param[in] filename Input file name.
1181   - * @param[in] evtname Events extension name.
1182   - *
1183   - * Extracts the Good Time Intervals extension name from the event file. We
1184   - * do this by loading the events and accessing the Good Time Intervals
1185   - * extension name using the GCTAEventList::gtiname() method. If the file name
1186   - * is empty, the method returns `GTI`.
1187   - ***************************************************************************/
1188   -std::string ctselect::get_gtiname(const std::string& filename,
1189   - const std::string& evtname) const
1190   -{
1191   - // Initialise GTI name
1192   - std::string gtiname = gammalib::extname_gti;
1193   -
1194   - // Continue only if the filename is not empty
1195   - if (!filename.empty()) {
1196   -
1197   - // Load events
1198   - GCTAEventList events(filename+"["+evtname+"]");
1199   -
1200   - // Get GTI name
1201   - gtiname = events.gtiname();
1202   -
1203   - }
1204   -
1205   - // Return GTI name
1206   - return (gtiname);
1207   -}
1208   -
1209   -
1210   -/***********************************************************************//**
1211 1178 * @brief Save event list in FITS format.
1212 1179 *
1213 1180 * Save the event list as a FITS file. The file name of the FITS file is
... ... @@ -1325,109 +1292,3 @@ void ctselect::save_xml(void)
1325 1292 // Return
1326 1293 return;
1327 1294 }
1328   -
1329   -
1330   -/***********************************************************************//**
1331   - * @brief Save event list into FITS file
1332   - *
1333   - * @param[in] obs Pointer to CTA observation.
1334   - * @param[in] infile Input file name.
1335   - * @param[in] evtname Event extension name.
1336   - * @param[in] gtiname GTI extension name.
1337   - * @param[in] outfile Output file name.
1338   - *
1339   - * Saves an event list and the corresponding Good Time Intervals into a FITS
1340   - * file and copy all others extensions from the input file to the output
1341   - * file.
1342   - *
1343   - * If an extension name is specified in the @p outfile argument, the events
1344   - * and eventually also the Good Time Intervals will be extracted from the
1345   - * argument and used for writing the events. The format is
1346   - *
1347   - * <filename>[<event extension name;GTI extension name>]
1348   - *
1349   - * where <filename> needs to be replaced by the name of the FITS file,
1350   - * and <event extension name;GTI extension name> by the name of the events
1351   - * and Good Time Intervals extensions. For example
1352   - *
1353   - * myfits.fits[EVENTS1;GTI1]
1354   - *
1355   - * will write the selected events into the "EVENTS1" extension and the
1356   - * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS
1357   - * file. If the Good Time Intervals extension name is skipped, e.g.
1358   - *
1359   - * myfits.fits[EVENTS1]
1360   - *
1361   - * the original extension name for the Good Time Intervals will be kept.
1362   - * Analogously, only the Good Time Intervals extension name can be changed
1363   - * by specifying
1364   - *
1365   - * myfits.fits[;GTI1]
1366   - *
1367   - * In none of the cases will the original events and Good Time Intervals be
1368   - * copied over to the output file.
1369   - ***************************************************************************/
1370   -void ctselect::save_event_list(const GCTAObservation* obs,
1371   - const std::string& infile,
1372   - const std::string& evtname,
1373   - const std::string& gtiname,
1374   - const std::string& outfile) const
1375   -{
1376   - // Save only if we have an event list
1377   - if (obs->eventtype() == "EventList") {
1378   -
1379   - // Set output FITS file event extension names
1380   - GFilename outname(outfile);
1381   - std::string outevt = evtname;
1382   - std::string outgti = gtiname;
1383   - if (outname.has_extname()) {
1384   - std::vector<std::string> extnames =
1385   - gammalib::split(outname.extname(), ";");
1386   - if (extnames.size() > 0) {
1387   - std::string extname = gammalib::strip_whitespace(extnames[0]);
1388   - if (!extname.empty()) {
1389   - outevt = extname;
1390   - }
1391   - }
1392   - if (extnames.size() > 1) {
1393   - std::string extname = gammalib::strip_whitespace(extnames[1]);
1394   - if (!extname.empty()) {
1395   - outgti = extname;
1396   - }
1397   - }
1398   - }
1399   -
1400   - // Create output FITS file
1401   - GFits outfits;
1402   -
1403   - // Write observation into FITS file
1404   - obs->write(outfits, outevt, outgti);
1405   -
1406   - // Copy all extensions other than evtname and gtiname extensions
1407   - // from the input to the output event list. The evtname and
1408   - // gtiname extensions are written by the save method, all others
1409   - // that may eventually be present have to be copied over
1410   - // explicitly.
1411   - GFits infits(infile);
1412   - for (int extno = 1; extno < infits.size(); ++extno) {
1413   - GFitsHDU* hdu = infits.at(extno);
1414   - if (hdu->extname() != evtname &&
1415   - hdu->extname() != gtiname &&
1416   - hdu->extname() != outevt &&
1417   - hdu->extname() != outgti) {
1418   - outfits.append(*hdu);
1419   - }
1420   - }
1421   -
1422   - // Close input file
1423   - infits.close();
1424   -
1425   - // Save file to disk and close it (we need both operations)
1426   - outfits.saveto(outname.url(), clobber());
1427   - outfits.close();
1428   -
1429   - } // endif: observation was unbinned
1430   -
1431   - // Return
1432   - return;
1433   -}
... ...
src/ctselect/ctselect.hpp
  View file @fc5614e
1 1 /***************************************************************************
2 2 * ctselect - Data selection tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -78,15 +78,8 @@ protected:
78 78 std::string check_infile(const std::string& filename,
79 79 const std::string& evtname) const;
80 80 std::string set_outfile_name(const std::string& filename) const;
81   - std::string get_gtiname(const std::string& filename,
82   - const std::string& evtname) const;
83 81 void save_fits(void);
84 82 void save_xml(void);
85   - void save_event_list(const GCTAObservation* obs,
86   - const std::string& infile,
87   - const std::string& evtname,
88   - const std::string& gtiname,
89   - const std::string& outfile) const;
90 83  
91 84 // User parameters
92 85 std::string m_outobs; //!< Output event list or XML file
... ...
src/support/ctool.cpp
  View file @fc5614e
... ... @@ -1730,6 +1730,145 @@ std::vector&lt;bool&gt; ctool::cube_layer_usage(const GEbounds&amp; cube_ebounds,
1730 1730  
1731 1731  
1732 1732 /***********************************************************************//**
  1733 + * @brief Save event list into FITS file
  1734 + *
  1735 + * @param[in] obs Pointer to CTA observation.
  1736 + * @param[in] infile Input file name.
  1737 + * @param[in] evtname Event extension name.
  1738 + * @param[in] gtiname GTI extension name.
  1739 + * @param[in] outfile Output file name.
  1740 + *
  1741 + * Saves an event list and the corresponding Good Time Intervals into a FITS
  1742 + * file and copy all others extensions from the input file to the output
  1743 + * file.
  1744 + *
  1745 + * If an extension name is specified in the @p outfile argument, the events
  1746 + * and eventually also the Good Time Intervals will be extracted from the
  1747 + * argument and used for writing the events. The format is
  1748 + *
  1749 + * <filename>[<event extension name;GTI extension name>]
  1750 + *
  1751 + * where <filename> needs to be replaced by the name of the FITS file,
  1752 + * and <event extension name;GTI extension name> by the name of the events
  1753 + * and Good Time Intervals extensions. For example
  1754 + *
  1755 + * myfits.fits[EVENTS1;GTI1]
  1756 + *
  1757 + * will write the selected events into the `EVENTS1` extension and the
  1758 + * Good Time Intervals into the `GTI1` extension of the `myfits.fits` FITS
  1759 + * file. If the Good Time Intervals extension name is skipped, e.g.
  1760 + *
  1761 + * myfits.fits[EVENTS1]
  1762 + *
  1763 + * the original extension name for the Good Time Intervals will be kept.
  1764 + * Analogously, only the Good Time Intervals extension name can be changed
  1765 + * by specifying
  1766 + *
  1767 + * myfits.fits[;GTI1]
  1768 + *
  1769 + * In none of the cases will the original events and Good Time Intervals be
  1770 + * copied over to the output file.
  1771 + ***************************************************************************/
  1772 +void ctool::save_event_list(const GCTAObservation* obs,
  1773 + const std::string& infile,
  1774 + const std::string& evtname,
  1775 + const std::string& gtiname,
  1776 + const std::string& outfile) const
  1777 +{
  1778 + // Save only if we have an event list
  1779 + if (obs->eventtype() == "EventList") {
  1780 +
  1781 + // Set output FITS file event extension names
  1782 + GFilename outname(outfile);
  1783 + std::string outevt = evtname;
  1784 + std::string outgti = gtiname;
  1785 + if (outname.has_extname()) {
  1786 + std::vector<std::string> extnames =
  1787 + gammalib::split(outname.extname(), ";");
  1788 + if (extnames.size() > 0) {
  1789 + std::string extname = gammalib::strip_whitespace(extnames[0]);
  1790 + if (!extname.empty()) {
  1791 + outevt = extname;
  1792 + }
  1793 + }
  1794 + if (extnames.size() > 1) {
  1795 + std::string extname = gammalib::strip_whitespace(extnames[1]);
  1796 + if (!extname.empty()) {
  1797 + outgti = extname;
  1798 + }
  1799 + }
  1800 + }
  1801 +
  1802 + // Create output FITS file
  1803 + GFits outfits;
  1804 +
  1805 + // Write observation into FITS file
  1806 + obs->write(outfits, outevt, outgti);
  1807 +
  1808 + // Copy all extensions other than evtname and gtiname extensions
  1809 + // from the input to the output event list. The evtname and
  1810 + // gtiname extensions are written by the save method, all others
  1811 + // that may eventually be present have to be copied over
  1812 + // explicitly.
  1813 + GFits infits(infile);
  1814 + for (int extno = 1; extno < infits.size(); ++extno) {
  1815 + GFitsHDU* hdu = infits.at(extno);
  1816 + if (hdu->extname() != evtname &&
  1817 + hdu->extname() != gtiname &&
  1818 + hdu->extname() != outevt &&
  1819 + hdu->extname() != outgti) {
  1820 + outfits.append(*hdu);
  1821 + }
  1822 + }
  1823 +
  1824 + // Close input file
  1825 + infits.close();
  1826 +
  1827 + // Save file to disk and close it (we need both operations)
  1828 + outfits.saveto(outname.url(), clobber());
  1829 + outfits.close();
  1830 +
  1831 + } // endif: observation was unbinned
  1832 +
  1833 + // Return
  1834 + return;
  1835 +}
  1836 +
  1837 +
  1838 +/***********************************************************************//**
  1839 + * @brief Get Good Time Intervals extension name
  1840 + *
  1841 + * @param[in] filename Input file name.
  1842 + * @param[in] evtname Events extension name.
  1843 + *
  1844 + * Extracts the Good Time Intervals extension name from the event file. We
  1845 + * do this by loading the events and accessing the Good Time Intervals
  1846 + * extension name using the GCTAEventList::gtiname() method. If the file name
  1847 + * is empty, the method returns `GTI`.
  1848 + ***************************************************************************/
  1849 +std::string ctool::get_gtiname(const std::string& filename,
  1850 + const std::string& evtname) const
  1851 +{
  1852 + // Initialise GTI name
  1853 + std::string gtiname = gammalib::extname_gti;
  1854 +
  1855 + // Continue only if the filename is not empty
  1856 + if (!filename.empty()) {
  1857 +
  1858 + // Load events
  1859 + GCTAEventList events(filename+"["+evtname+"]");
  1860 +
  1861 + // Get GTI name
  1862 + gtiname = events.gtiname();
  1863 +
  1864 + }
  1865 +
  1866 + // Return GTI name
  1867 + return (gtiname);
  1868 +}
  1869 +
  1870 +
  1871 +/***********************************************************************//**
1733 1872 * @brief Set warning string if there are too few energies
1734 1873 *
1735 1874 * @param[in] energies Energies.
... ...
src/support/ctool.hpp
  View file @fc5614e
1 1 /***************************************************************************
2 2 * ctool - ctool base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -121,6 +121,13 @@ protected:
121 121 std::vector<bool> cube_layer_usage(const GEbounds& cube_ebounds,
122 122 const GEbounds& list_ebounds) const;
123 123 bool is_valid_filename(const GFilename& filename) const;
  124 + std::string get_gtiname(const std::string& filename,
  125 + const std::string& evtname) const;
  126 + void save_event_list(const GCTAObservation* obs,
  127 + const std::string& infile,
  128 + const std::string& evtname,
  129 + const std::string& gtiname,
  130 + const std::string& outfile) const;
124 131  
125 132 // Protected warning strings
126 133 std::string warn_too_few_energies(const GEnergies& energies) const;
... ...
test/test_ctprob.py
  View file @fc5614e
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # This scripts performs unit tests for the ctprob tool.
4 4 #
5   -# Copyright (C) 2014-2016 Juergen Knoedlseder
  5 +# Copyright (C) 2017 by Leonardo Di Venere
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -73,7 +73,7 @@ class Test(test):
73 73 ' outobs="ctprob_cmd1.fits"'+ \
74 74 ' inmodel='+self._model+ \
75 75 ' caldb='+self._caldb+' irf='+self._irf+ \
76   - ' edisp=yes'+ \
  76 + ' edisp=no'+ \
77 77 ' logfile="ctprob_cmd1.log" chatter=1'
78 78  
79 79 # Check if execution was successful
... ... @@ -88,7 +88,7 @@ class Test(test):
88 88 ' outobs="ctprob_cmd2.fits"'+ \
89 89 ' inmodel='+self._model+ \
90 90 ' caldb='+self._caldb+' irf='+self._irf+ \
91   - ' edisp=yes'+ \
  91 + ' edisp=no'+ \
92 92 ' logfile="ctprob_cmd2.log" chatter=1'
93 93  
94 94 # Check if execution failed
... ... @@ -135,12 +135,13 @@ class Test(test):
135 135 # Now set ctprob parameters
136 136 prob['inobs'] = self._events
137 137 prob['outobs'] = 'ctprob_py1.fits'
138   - prob['inmodel'] = self._model
139   - prob['caldb'] = self._caldb
140   - prob['irf'] = self._irf
141   - prob['edisp'] = True
  138 + prob['inmodel'] = self._model
  139 + prob['caldb'] = self._caldb
  140 + prob['irf'] = self._irf
  141 + prob['edisp'] = False
142 142 prob['logfile'] = 'ctprob_py1.log'
143 143 prob['chatter'] = 2
  144 +
144 145 # Run ctprob tool
145 146 prob.logFileOpen() # Make sure we get a log file
146 147 prob.run()
... ... @@ -149,15 +150,14 @@ class Test(test):
149 150 # Check result file
150 151 self._check_result_file('ctprob_py1.fits')
151 152  
152   -
153 153 # Setup ctprob tool for an invalid event file
154 154 prob = ctools.ctprob()
155 155 prob['inobs'] = self._invalid_events
156 156 prob['outobs'] = 'ctprob_py2.fits'
157   - prob['inmodel'] = self._model
158   - prob['caldb'] = self._caldb
159   - prob['irf'] = self._irf
160   - prob['edisp'] = True
  157 + prob['inmodel'] = self._model
  158 + prob['caldb'] = self._caldb
  159 + prob['irf'] = self._irf
  160 + prob['edisp'] = False
161 161 prob['logfile'] = 'ctprob_py2.log'
162 162 prob['chatter'] = 3
163 163  
... ... @@ -167,26 +167,24 @@ class Test(test):
167 167 try:
168 168 prob.execute()
169 169 self.test_try_failure('Exception not thrown')
170   - except ValueError:
  170 + except RuntimeError:
171 171 self.test_try_success()
172 172  
173   -
174 173 # Setup ctprob with a list of observations and test prefix
175 174 prob = ctools.ctprob()
176   - prob['inobs'] = self._datadir + '/obs_unbinned.xml'
177   - prob['prefix'] = 'ctprob_py3_'
  175 + prob['inobs'] = self._datadir + '/obs_unbinned.xml'
  176 + prob['prefix'] = 'ctprob_py3_'
178 177 prob['outobs'] = 'ctprob_py3_obs.xml'
179   - prob['inmodel'] = self._model
180   - prob['caldb'] = self._caldb
181   - prob['irf'] = self._irf
182   - prob['edisp'] = True
  178 + prob['inmodel'] = self._model
  179 + prob['caldb'] = self._caldb
  180 + prob['irf'] = self._irf
  181 + prob['edisp'] = True
183 182 prob['logfile'] = 'ctprob_py3.log'
184   - prob['chatter'] = 3
  183 + prob['chatter'] = 4
185 184  
186 185 # Run ctprob tool
187 186 prob.logFileOpen() # Make sure we get a log file
188   - prob.run()
189   - prob.save()
  187 + prob.execute()
190 188  
191 189 # Check result file
192 190 self._check_result_file('ctprob_py3_crab_events.fits')
... ... @@ -213,8 +211,8 @@ class Test(test):
213 211 self._check_events(events, nevents=nevents)
214 212  
215 213 # Check that probability columns exist
216   - colname1 = "PROB_Background"
217   - colname2 = "PROB_Crab"
  214 + colname1 = 'PROB_Background'
  215 + colname2 = 'PROB_Crab'
218 216 self._check_column(filename, colname1)
219 217 self._check_column(filename, colname2)
220 218  
... ... @@ -291,8 +289,8 @@ class Test(test):
291 289  
292 290 def _check_normalization(self, filename, colnames):
293 291 """
294   - Check that the probability values in columns identified
295   - by colnames are correctly normalized.
  292 + Check that the probability values in columns identified by colnames are
  293 + correctly normalized.
296 294  
297 295 Parameters
298 296 ----------
... ... @@ -304,26 +302,31 @@ class Test(test):
304 302 # Check that the probability values are correctly normalized
305 303  
306 304 # Open fits table
307   - gfits = gammalib.GFits(filename)
  305 + gfits = gammalib.GFits(filename)
308 306 gtable = gfits[1]
309   - ncols = gtable.ncols()
310   - nevt = gtable.nrows()
311   - srcs=[]
  307 + ncols = gtable.ncols()
  308 + nevt = gtable.nrows()
  309 + srcs = []
  310 +
312 311 # Loop over colnames
313 312 for colname in colnames:
314   - # self._check_column(filename, colnames) # Check if column exists
  313 +
315 314 # Get the correct column corresponding to colname from the table
316 315 for icol in range(ncols):
317 316 gcol = gtable[icol]
318 317 if gcol.name()==colname:
319 318 break
  319 +
320 320 # Read and save column content
321 321 src1 = []
322 322 for ievt in range(nevt):
323 323 src1.append(gcol[ievt])
324 324 srcs.extend(src1)
  325 +
  326 + # Compute sum
325 327 tot = sum(srcs)
326   - self.test_value(int(tot+0.5),nevt, 'Check that probability columns are normalized')
  328 + self.test_value(int(tot+0.5), nevt, 'Check that probability columns '
  329 + 'are normalized')
327 330  
328 331 # Return
329 332 return
... ...