diff --git a/doc/source/users/reference_manual/ctprob.rst b/doc/source/users/reference_manual/ctprob.rst index d7fcf0b..5da08e3 100644 --- a/doc/source/users/reference_manual/ctprob.rst +++ b/doc/source/users/reference_manual/ctprob.rst @@ -1,24 +1,24 @@ - .. _ctprob: +.. _ctprob: ctprob ======== -Computes event probability for a given model. +Compute event probability for a given model. Synopsis -------- -This tool takes an event list file (or an XML observation file) and a source -model as input. For each event, it evalutes the probability for that event -to come from each of the source in the model. This probability is obtained -evaluating the differential expected counts for each source at the event -direction and energy and normalizing it to the total differential expected -counts for the given model. If only one source is provided in the model, -the probability is set to 1 for each photon. -A column for each source in the model is added to the event list in order to -store these values. The name of each column is made by the source name with -the prefix "PROB_". +This tool takes an event list file (or an XML observation definition file) and +a model as input and computes for each event the probability that it originates +from a specific model component. The probability is computed by by evaluating +the differential event probability for a given model and normalizing these +probabilities so that the sum of the probabilities for all model components is +unity for each event. + +The tool appends for each model component a single precision column to the event +list that contains the event probabilities. The name of the column is build from +the model name prefixed with ``PROB_``. General parameters @@ -34,7 +34,7 @@ General parameters Prefix for output event lists in observation definition file. ``inmodel [string]`` - Input model XML file. + Input model definition XML file. ``caldb [string]`` Calibration database. @@ -43,7 +43,7 @@ General parameters Instrument response function. ``(edisp = no) [boolean]`` - Apply energy dispersion to response computation. + Applies energy dispersion to response computation. Standard parameters @@ -73,7 +73,7 @@ Standard parameters ``(mode = ql) [string]`` Mode of automatic parameters (default is "ql", i.e. "query and learn"). -``(logfile = ctlike.log) [string]`` +``(logfile = ctprob.log) [string]`` Name of log file. diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst index a67845e..d4a2c2c 100644 --- a/doc/source/users/reference_manual/index.rst +++ b/doc/source/users/reference_manual/index.rst @@ -27,11 +27,11 @@ ctools ctmodel --- Computes model counts cube ctobssim --- Simulate observations ctpsfcube --- Generates point spread function cube + ctprob --- Computes event probability for a given model ctselect --- Selects event data ctskymap --- Generates sky map cttsmap --- Generates Test Statistic map ctulimit --- Calculates upper limit - ctprob --- Computes event probability for a given model cscripts diff --git a/src/ctprob/Makefile.am b/src/ctprob/Makefile.am index dc0bdac..1ef279c 100644 --- a/src/ctprob/Makefile.am +++ b/src/ctprob/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# -# Makefile.am for ctprob tool # +# Makefile.am for ctprob tool # # ------------------------------------------------------------------------- # -# Copyright (C) 2016 Juergen Knoedlseder # +# Copyright (C) 2017 Leonardo Di Venere # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # diff --git a/src/ctprob/ctprob.cpp b/src/ctprob/ctprob.cpp index b574601..7e77765 100644 --- a/src/ctprob/ctprob.cpp +++ b/src/ctprob/ctprob.cpp @@ -1,7 +1,7 @@ /*************************************************************************** - * ctprob - Computes probability for a given model * + * ctprob - Computes event probability for a given model * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2016 by Leonardo Di Venere * + * copyright (C) 2017 by Leonardo Di Venere * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file ctprob.cpp - * @brief Computes probability for a given model + * @brief Computes event probability for a given model * @author Leonardo Di Venere */ @@ -33,13 +33,12 @@ #include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ -#define G_RUN "ctprob::run()" -#define G_GET_OBS "ctprob::get_obs()" +#define G_GET_OBS "ctprob::get_obs()" +#define G_EVALUATE_PROB "ctprob::evaluate_probability(GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ /* __ Coding definitions _________________________________________________ */ -//#define G_USE_MKSTEMP //!< Use mkstemp for temporary filename /*========================================================================== @@ -51,7 +50,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctprob::ctprob(void) : ctobservation(CTPROB_NAME, CTPROB_VERSION) +ctprob::ctprob(void) : ctobservation(CTPROB_NAME, VERSION) { // Initialise members init_members(); @@ -70,7 +69,7 @@ ctprob::ctprob(void) : ctobservation(CTPROB_NAME, CTPROB_VERSION) * provided in an observation container. ***************************************************************************/ ctprob::ctprob(const GObservations& obs) : - ctobservation(CTPROB_NAME, CTPROB_VERSION, obs) + ctobservation(CTPROB_NAME, VERSION, obs) { // Initialise members init_members(); @@ -89,7 +88,7 @@ ctprob::ctprob(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctprob::ctprob(int argc, char *argv[]) : - ctobservation(CTPROB_NAME, CTPROB_VERSION, argc, argv) + ctobservation(CTPROB_NAME, VERSION, argc, argv) { // Initialise members init_members(); @@ -201,13 +200,11 @@ void ctprob::clear(void) /***********************************************************************//** - * @brief Append probability columns event data + * @brief Runs the ctprob tool * * This method reads in the application parameters and loops over all - * observations that were found to perform the probability calculation. - * Each observation is written to a temporary file, which is then - * re-opened and used for the actual calculation. The temporary file - * is deleted after this action so that no disk overflow will occur. + * unbinned observations to compute the probability for each event that it + * comes from one of the model components. ***************************************************************************/ void ctprob::run(void) { @@ -219,123 +216,30 @@ void ctprob::run(void) // Get parameters get_parameters(); + // Set energy dispersion flags of all CTA observations and save old + // values in save_edisp vector + std::vector save_edisp = set_edisp(m_obs, m_apply_edisp); + // Write input observation container into logger log_observations(NORMAL, m_obs, "Input observation"); // Write header into logger - log_header1(TERSE, "Event probability evaluation"); + log_header1(TERSE, "Compute event probabilities"); // Initialise counters int n_observations = 0; - // Loop over all observation in the container - for (int i = 0; i < m_obs.size(); ++i) { - - // Write header for the current observation - log_header3(TERSE, get_obs_header(m_obs[i])); - - // Initialise event input and output filenames and the event - // and GTI extension names - m_infiles.push_back(""); - m_evtname.push_back(gammalib::extname_cta_events); - m_gtiname.push_back(gammalib::extname_gti); - - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); - - // Skip observation if it's not CTA - if (obs == NULL) { - std::string msg = " Skipping "+m_obs[i]->instrument()+ - " observation"; - log_string(NORMAL, msg); - continue; - } - - // Skip observation if we have a binned observation - if (obs->eventtype() == "CountsCube") { - std::string msg = " Skipping binned "+m_obs[i]->instrument()+ - " observation"; - log_string(NORMAL, msg); - continue; - } + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { // Increment counter n_observations++; - // Save event file name (for possible saving) - m_infiles[i] = obs->eventfile(); - - // Extract event and GTI extension names from input FITS file - GFilename fname(m_infiles[i]); - if (fname.has_extname()) { - m_evtname[i] = fname.extname(); - } - m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]); - - // Write input file information into logger - log_value(NORMAL, "Input filename", m_infiles[i]); - log_value(NORMAL, "Event extension name", m_evtname[i]); - log_value(NORMAL, "GTI extension name", m_gtiname[i]); - - // Fall through in case that the event file is empty - if (obs->events()->size() == 0) { - log_string(NORMAL, " Warning: No events in event file \""+ - m_infiles[i]+"\". Event probability evaluation skipped."); - continue; - } - - // If we have an input file then check it - if (!m_infiles[i].empty()) { - std::string message = check_infile(m_infiles[i], m_evtname[i]); - if (!message.empty()) { - throw GException::invalid_value(G_RUN, message); - } - } - - // Get temporary file name - #if G_USE_MKSTEMP - char tpl[] = "ctprobXXXXXX"; - int fileid = mkstemp(tpl); - std::string filename(tpl); - #else - std::string filename = std::tmpnam(NULL); - #endif - - // Save observation in temporary file. We add here the events and - // GTI extension name so that the GCTAObservation::save method can - // use this information for writing the proper extension names into - // the temporary file - obs->save(filename+"["+m_evtname[i]+";"+ m_gtiname[i]+"]", true); - - // Log saved FITS file. - if (logExplicit()) { - GFits tmpfile(filename); - log.header3("FITS file content of temporary file"); - log << tmpfile << std::endl; - tmpfile.close(); - } - - // If we have a temporary file then check it - if (!filename.empty()) { - std::string message = check_infile(filename, m_evtname[i]); - if (!message.empty()) { - throw GException::invalid_value(G_RUN, message); - } - } - - // Load observation from temporary file, including event selection + // Evaluate event probabilities evaluate_probability(obs); - //select_events(obs, filename, m_evtname[i], m_gtiname[i]); - - // Close temporary file - #if G_USE_MKSTEMP - close(fileid); - #endif - // Remove temporary file - std::remove(filename.c_str()); - - } // endfor: looped over all observations + } // endfor: looped over observations // If more than a single observation has been handled then make sure that // an XML file will be used for storage @@ -346,6 +250,9 @@ void ctprob::run(void) // Write observation(s) into logger log_observations(NORMAL, m_obs, "Output observation"); + // Restore energy dispersion flags of all CTA observations + restore_edisp(m_obs, save_edisp); + // Optionally publish event list(s) if ((*this)["publish"].boolean()) { publish(); @@ -379,6 +286,7 @@ void ctprob::save(void) { // Write header into logger log_header1(TERSE, gammalib::number("Save event list", m_obs.size())); + // Case A: Save event file(s) and XML metadata information if (m_use_xml) { save_xml(); @@ -544,11 +452,6 @@ void ctprob::free_members(void) ***************************************************************************/ void ctprob::get_parameters(void) { - // Initialise selection flags - /* m_select_energy = true; - m_select_roi = true; - m_select_time = true; - */ // Setup observations from "inobs" parameter. Do not request response // information and do not accept counts cubes. if (m_obs.size() == 0) { @@ -559,7 +462,6 @@ void ctprob::get_parameters(void) else { setup_observations(m_obs); } - //setup_observations(m_obs, false, true, false); // Read model definition file if required if (m_obs.models().size() == 0) { @@ -671,144 +573,90 @@ void ctprob::get_obs(void) * * @param[in,out] obs CTA observation. * + * @exception GException::invalid_value + * Observation does not contain an event list. * - * For each event from the input event file evaluates the probability that - * the event comes from any of the source in the input model. This is done - * by evaluating the differential expected counts for each source at the event - * direction and energy and normalizing it to the total differential expected - * counts for the given model. If only one source is provided in the model, - * the probability is set to 1 for each photon. + * Computes for all events the probability that an event arises from a + * specific model component. The computation is done by by evaluating the + * differential event probability for a given model and normalizing these + * probabilities so that the sum of the probabilities for all model + * components is unity for each event. * - * A FITS column for each source in the model is created and added to - * the event list. The name of each column is made by the source name with - * the prefix "PROB_". + * The method appends for each model component a single precision column to + * the event list that contains the event probabilities. The name of the + * column is build from the model name prefixed with `PROB_`. ***************************************************************************/ -void ctprob::evaluate_probability(GCTAObservation* obs) +void ctprob::evaluate_probability(GCTAObservation* obs) { // Write header into logger - log_header3(NORMAL, "Evaluating probability for events"); - - GCTAEventList* evts = static_cast(const_cast(obs->events())); - const GEvent* evt; - - double total=0.; - double value; - - std::vector columns; - - //define colums - for (int j = 0; jname(); - GFitsTableFloatCol* col = new GFitsTableFloatCol("PROB_"+mdl_name, evts->size()); - columns.push_back(col); - } - - // Loop over events - for (int i = 0; i < evts->size(); ++i) { - total = 0.; - evt = (*evts)[i]; - std::vector values; - //Check if you have more than one model - if (m_obs.models().size()>1) { - //Loop over models - for (int j = 0; jeval( *evt, *obs); - values.push_back(value); - total += value; - } - } - // if there is only one model a constant value equal to 1. is set for each photon - else { - value = 1.; - values.push_back(value); - total += value; - } // end if - for (int j = 0; jappend_column(*columns[j]); - } - - // Return - return; -} - - - -/***********************************************************************//** - * @brief Check input filename - * - * @param[in] filename File name. - * @param[in] evtname Event extension name. - * - * This method checks if the input FITS file is correct. - ***************************************************************************/ -std::string ctprob::check_infile(const std::string& filename, - const std::string& evtname) const -{ - // Initialise message string - std::string message = ""; - - // Open FITS file - GFits fits(filename); - - // Check for existence of events extensions - if (!fits.contains(evtname)) { - message = "No \""+evtname+"\" extension found in input file \""+ - filename + "\"."; + log_header3(NORMAL, "Evaluate events probabilities"); + + // Make sure that the observation holds a CTA event list. If this + // is not the case then throw an exception. + GCTAEventList* events = dynamic_cast(const_cast + (obs->events())); + if (events == NULL) { + std::string msg = "CTA Observation does not contain an event list. An " + "event list is needed to compute the event " + "probabilities."; + throw GException::invalid_value(G_EVALUATE_PROB, msg); } - // ... otherwise check column names - else { - - // Get pointer to FITS table - GFitsTable* table = fits.table(evtname); + // Continue only if there are model components + int nmodels = m_obs.models().size(); + if (nmodels > 0) { + + // Allocate single precision FITS table columns for all model components + std::vector columns; + for (int k = 0; k < nmodels; ++k) { + const std::string colname = "PROB_" + m_obs.models()[k] ->name(); + GFitsTableFloatCol* col = + new GFitsTableFloatCol(colname, events->size()); + columns.push_back(col); + } - // Initialise list of missing columns - std::vector missing; + // Loop over all events in the event list + for (int i = 0; i < events->size(); ++i) { - // Check for existence of TIME column - if (!table->contains("TIME")) { - missing.push_back("TIME"); - } + // Initialise total probability + double total = 0.0; - // Check for existence of ENERGY column - if (!table->contains("ENERGY")) { - missing.push_back("ENERGY"); - } + // Get point to event + const GEvent* event = (*events)[i]; - // Check for existence of RA column - if (!table->contains("RA")) { - missing.push_back("RA"); - } + // Initialise probability vector + std::vector probabilities(nmodels, 0.0); - // Check for existence of DEC column - if (!table->contains("DEC")) { - missing.push_back("DEC"); - } + // Loop over models + for (int k = 0; k < nmodels; ++k) { + double value = m_obs.models()[k]->eval(*event, *obs); + probabilities[k] = value; + total += value; + } - // Set error message for missing columns - if (!missing.empty()) { - message = "The following columns are missing in the " - "\""+evtname+"\" extension of input file \""+ - filename + "\": "; - for (int i = 0; i < missing.size(); ++i) { - message += "\"" + missing[i] + "\""; - if (i < missing.size()-1) { - message += ", "; + // Normalize probabilities to unity + if (total > 0.0) { + for (int k = 0; k < nmodels; ++k) { + probabilities[k] /= total; } } + + // Store probabilities in FITS columns + for (int k = 0; k < nmodels; ++k) { + (*(columns[k]))(i) = probabilities[k]; + } + + } // endfor: looped over all events + + // Append columns to event list + for (int k = 0; k < nmodels; ++k) { + events->append_column(*columns[k]); } - } // endelse: checked column names + } // endif: there were model components // Return - return message; + return; } @@ -842,43 +690,10 @@ std::string ctprob::set_outfile_name(const std::string& filename) const /***********************************************************************//** - * @brief Get Good Time Intervals extension name - * - * @param[in] filename Input file name. - * @param[in] evtname Events extension name. - * - * Extracts the Good Time Intervals extension name from the event file. We - * do this by loading the events and accessing the Good Time Intervals - * extension name using the GCTAEventList::gtiname() method. If the file name - * is empty, the method returns `GTI`. - ***************************************************************************/ -std::string ctprob::get_gtiname(const std::string& filename, - const std::string& evtname) const -{ - // Initialise GTI name - std::string gtiname = gammalib::extname_gti; - - // Continue only if the filename is not empty - if (!filename.empty()) { - - // Load events - GCTAEventList events(filename+"["+evtname+"]"); - - // Get GTI name - gtiname = events.gtiname(); - - } - - // Return GTI name - return (gtiname); -} - - -/***********************************************************************//** * @brief Save event list in FITS format. * * Save the event list as a FITS file. The file name of the FITS file is - * specified by the "outobs" parameter. + * specified by the `outobs` parameter. ***************************************************************************/ void ctprob::save_fits(void) { @@ -888,39 +703,42 @@ void ctprob::save_fits(void) // Get output filename m_outobs = (*this)["outobs"].filename(); - // Get CTA observation from observation container - GCTAObservation* obs = dynamic_cast(m_obs[0]); + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { - // Save only if it's a CTA observation - if (obs != NULL) { - - // Save only if file name is non-empty - if (m_infiles[0].length() > 0) { + // Get input file name and default extension names + std::string infile = obs->eventfile(); + std::string evtname = gammalib::extname_cta_events; + std::string gtiname = gammalib::extname_gti; - // Create file name object - GFilename fname(m_outobs); + // Create file name object + GFilename fname(m_outobs); - // Extract filename and event extension name - std::string outfile = fname.url(); + // Extract filename and event extension name + std::string outfile = fname.url(); - // Append event extension name. We handle here the possibility - // to write the events into a different extension. - if (fname.has_extname()) { - outfile += "["+fname.extname()+"]"; - } - else { - outfile += "["+m_evtname[0]+"]"; - } - - // Log filename - log_value(NORMAL, "Event list file", outfile); - // Save event list - save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0], - outfile); + // Append event extension name. We handle here the possibility + // to write the events into a different extension. + if (fname.has_extname()) { + outfile += "["+fname.extname()+"]"; + } + else { + outfile += "["+gammalib::extname_cta_events+"]"; + } + + // Log filename + log_value(NORMAL, "Event list file", outfile); + + // Save event list + save_event_list(obs, infile, evtname, gtiname, outfile); - } // endif: filename was non empty + // Exit the loop (if this method is called we should only + // have a single unbinned observation in the observation + // container + break; - } // endif: observation was CTA observation + } // endfor: looped over all unbinned observations } // endif: there were observations @@ -933,7 +751,7 @@ void ctprob::save_fits(void) * @brief Save event list(s) in XML format. * * Save the event list(s) into FITS files and write the file path information - * into a XML file. The filename of the XML file is specified by the outfile + * into a XML file. The filename of the XML file is specified by the `outobs` * parameter, the filename(s) of the event lists are built by prepending a * prefix to the input event list filenames. Any path present in the input * filename will be stripped, i.e. the event list(s) will be written in the @@ -948,27 +766,27 @@ void ctprob::save_xml(void) // Issue warning if output filename has no .xml suffix log_string(TERSE, warn_xml_suffix(m_outobs)); - // Loop over all observation in the container - for (int i = 0; i < m_obs.size(); ++i) { - - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { - // Skip observations that are no CTA observations - if (obs == NULL) { - continue; - } + // Get input file name and default extension names + std::string infile = obs->eventfile(); + std::string evtname = gammalib::extname_cta_events; + std::string gtiname = gammalib::extname_gti; - // Skip observations that have empty names - if (m_infiles[i].length() == 0) { - continue; + // Extract event and GTI extension names from input FITS file + GFilename fname(infile); + if (fname.has_extname()) { + evtname = fname.extname(); } + gtiname = get_gtiname(fname.url(), evtname); // Set event output file name - std::string outfile = set_outfile_name(m_infiles[i]); + std::string outfile = set_outfile_name(infile); // Append event extension name - outfile += "["+m_evtname[i]+"]"; + outfile += "["+evtname+"]"; // Log filename log_value(NORMAL, "Event list file", outfile); @@ -977,8 +795,7 @@ void ctprob::save_xml(void) obs->eventfile(outfile); // Save event list - save_event_list(obs, m_infiles[i], m_evtname[i], m_gtiname[i], - outfile); + save_event_list(obs, infile, evtname, gtiname, outfile); } // endfor: looped over observations @@ -991,109 +808,3 @@ void ctprob::save_xml(void) // Return return; } - - -/***********************************************************************//** - * @brief Save event list into FITS file - * - * @param[in] obs Pointer to CTA observation. - * @param[in] infile Input file name. - * @param[in] evtname Event extension name. - * @param[in] gtiname GTI extension name. - * @param[in] outfile Output file name. - * - * Saves an event list and the corresponding Good Time Intervals into a FITS - * file and copy all others extensions from the input file to the output - * file. - * - * If an extension name is specified in the @p outfile argument, the events - * and eventually also the Good Time Intervals will be extracted from the - * argument and used for writing the events. The format is - * - * [] - * - * where needs to be replaced by the name of the FITS file, - * and by the name of the events - * and Good Time Intervals extensions. For example - * - * myfits.fits[EVENTS1;GTI1] - * - * will write the selected events into the "EVENTS1" extension and the - * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS - * file. If the Good Time Intervals extension name is skipped, e.g. - * - * myfits.fits[EVENTS1] - * - * the original extension name for the Good Time Intervals will be kept. - * Analogously, only the Good Time Intervals extension name can be changed - * by specifying - * - * myfits.fits[;GTI1] - * - * In none of the cases will the original events and Good Time Intervals be - * copied over to the output file. - ***************************************************************************/ -void ctprob::save_event_list(const GCTAObservation* obs, - const std::string& infile, - const std::string& evtname, - const std::string& gtiname, - const std::string& outfile) const -{ - // Save only if we have an event list - if (obs->eventtype() == "EventList") { - - // Set output FITS file event extension names - GFilename outname(outfile); - std::string outevt = evtname; - std::string outgti = gtiname; - if (outname.has_extname()) { - std::vector extnames = - gammalib::split(outname.extname(), ";"); - if (extnames.size() > 0) { - std::string extname = gammalib::strip_whitespace(extnames[0]); - if (!extname.empty()) { - outevt = extname; - } - } - if (extnames.size() > 1) { - std::string extname = gammalib::strip_whitespace(extnames[1]); - if (!extname.empty()) { - outgti = extname; - } - } - } - - // Create output FITS file - GFits outfits; - - // Write observation into FITS file - obs->write(outfits, outevt, outgti); - - // Copy all extensions other than evtname and gtiname extensions - // from the input to the output event list. The evtname and - // gtiname extensions are written by the save method, all others - // that may eventually be present have to be copied over - // explicitly. - GFits infits(infile); - for (int extno = 1; extno < infits.size(); ++extno) { - GFitsHDU* hdu = infits.at(extno); - if (hdu->extname() != evtname && - hdu->extname() != gtiname && - hdu->extname() != outevt && - hdu->extname() != outgti) { - outfits.append(*hdu); - } - } - - // Close input file - infits.close(); - - // Save file to disk and close it (we need both operations) - outfits.saveto(outname.url(), clobber()); - outfits.close(); - - } // endif: observation was unbinned - - // Return - return; -} diff --git a/src/ctprob/ctprob.hpp b/src/ctprob/ctprob.hpp index bdf2dcf..607d2e8 100644 --- a/src/ctprob/ctprob.hpp +++ b/src/ctprob/ctprob.hpp @@ -1,7 +1,7 @@ /*************************************************************************** - * ctprob - Computes probability for a given model * + * ctprob - Computes event probability for a given model * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2016 by Leonardo Di Venere * + * copyright (C) 2017 by Leonardo Di Venere * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file ctprob.hpp - * @brief Computes probability for a given model + * @brief Computes event probability for a given model * @author Leonardo Di Venere */ @@ -35,8 +35,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTPROB_NAME "ctprob" -#define CTPROB_VERSION "0.0.1" +#define CTPROB_NAME "ctprob" /***********************************************************************//** @@ -70,32 +69,22 @@ protected: void free_members(void); void get_parameters(void); void get_obs(void); - void evaluate_probability(GCTAObservation* obs); - std::string check_infile(const std::string& filename, - const std::string& evtname) const; + void evaluate_probability(GCTAObservation* obs); std::string set_outfile_name(const std::string& filename) const; - std::string get_gtiname(const std::string& filename, - const std::string& evtname) const; void save_fits(void); void save_xml(void); - void save_event_list(const GCTAObservation* obs, - const std::string& infile, - const std::string& evtname, - const std::string& gtiname, - const std::string& outfile) const; // User parameters - std::string m_outobs; //!< Output event list or XML file - std::string m_prefix; //!< Prefix for multiple event lists - bool m_apply_edisp; //!< Apply energy dispersion? - bool m_publish; //!< Publish model cube? - GChatter m_chatter; //!< Chattiness + std::string m_outobs; //!< Output event list or XML file + std::string m_prefix; //!< Prefix for multiple event lists + bool m_apply_edisp; //!< Apply energy dispersion? + bool m_publish; //!< Publish event list? + GChatter m_chatter; //!< Chattiness // Protected members - std::vector m_infiles; //!< Input event filenames - std::vector m_evtname; //!< Event extension names - std::vector m_gtiname; //!< GTI extension names + std::vector m_infiles; //!< Input event filenames + std::vector m_evtname; //!< Event extension names + std::vector m_gtiname; //!< GTI extension names }; - #endif /* CTPROB_HPP */ diff --git a/src/ctprob/ctprob.par b/src/ctprob/ctprob.par index edbb1dd..39ee4ec 100644 --- a/src/ctprob/ctprob.par +++ b/src/ctprob/ctprob.par @@ -8,7 +8,7 @@ # # # Author: Leonardo Di Venere # # leonardo.divenere@ba.infn.it # -# INFN # +# INFN # # # # Purpose: Parameter file for the ctprob application # # # @@ -16,13 +16,13 @@ # # General parameters #=================== -inobs, f, a, "events.fits",,, "Input event list or observation definition XML file" -outobs, f, a, "prob_events.fits",,, "Output event list or observation definition XML file" -prefix, s, h, "prob_",,, "Prefix for event lists in observation definition XML file" -inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file" -caldb, s, a, "prod2",,, "Calibration database" -irf, s, a, "South_0.5h",,, "Instrument response function" -edisp, b, h, no,,, "Apply energy dispersion?" +inobs, f, a, "events.fits",,, "Input event list or observation definition XML file" +outobs, f, a, "prob_events.fits",,, "Output event list or observation definition XML file" +prefix, s, h, "prob_",,, "Prefix for event lists in observation definition XML file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file" +caldb, s, a, "prod2",,, "Calibration database" +irf, s, a, "South_0.5h",,, "Instrument response function" +edisp, b, h, no,,, "Apply energy dispersion?" # # Standard parameters diff --git a/src/ctprob/main.cpp b/src/ctprob/main.cpp index 5b7a804..05257c9 100644 --- a/src/ctprob/main.cpp +++ b/src/ctprob/main.cpp @@ -1,7 +1,7 @@ /*************************************************************************** - * ctprob - Computes probability for a given model * + * ctprob - Computes event probability for a given model * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2016 by Leonardo Di Venere * + * copyright (C) 2017 by Leonardo Di Venere * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file main.cpp - * @brief Computes probability for a given model + * @brief Computes event probability for a given model * @author Leonardo Di Venere */ diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 870d8f0..7b7afdf 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -1175,39 +1175,6 @@ std::string ctselect::set_outfile_name(const std::string& filename) const /***********************************************************************//** - * @brief Get Good Time Intervals extension name - * - * @param[in] filename Input file name. - * @param[in] evtname Events extension name. - * - * Extracts the Good Time Intervals extension name from the event file. We - * do this by loading the events and accessing the Good Time Intervals - * extension name using the GCTAEventList::gtiname() method. If the file name - * is empty, the method returns `GTI`. - ***************************************************************************/ -std::string ctselect::get_gtiname(const std::string& filename, - const std::string& evtname) const -{ - // Initialise GTI name - std::string gtiname = gammalib::extname_gti; - - // Continue only if the filename is not empty - if (!filename.empty()) { - - // Load events - GCTAEventList events(filename+"["+evtname+"]"); - - // Get GTI name - gtiname = events.gtiname(); - - } - - // Return GTI name - return (gtiname); -} - - -/***********************************************************************//** * @brief Save event list in FITS format. * * 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) // Return return; } - - -/***********************************************************************//** - * @brief Save event list into FITS file - * - * @param[in] obs Pointer to CTA observation. - * @param[in] infile Input file name. - * @param[in] evtname Event extension name. - * @param[in] gtiname GTI extension name. - * @param[in] outfile Output file name. - * - * Saves an event list and the corresponding Good Time Intervals into a FITS - * file and copy all others extensions from the input file to the output - * file. - * - * If an extension name is specified in the @p outfile argument, the events - * and eventually also the Good Time Intervals will be extracted from the - * argument and used for writing the events. The format is - * - * [] - * - * where needs to be replaced by the name of the FITS file, - * and by the name of the events - * and Good Time Intervals extensions. For example - * - * myfits.fits[EVENTS1;GTI1] - * - * will write the selected events into the "EVENTS1" extension and the - * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS - * file. If the Good Time Intervals extension name is skipped, e.g. - * - * myfits.fits[EVENTS1] - * - * the original extension name for the Good Time Intervals will be kept. - * Analogously, only the Good Time Intervals extension name can be changed - * by specifying - * - * myfits.fits[;GTI1] - * - * In none of the cases will the original events and Good Time Intervals be - * copied over to the output file. - ***************************************************************************/ -void ctselect::save_event_list(const GCTAObservation* obs, - const std::string& infile, - const std::string& evtname, - const std::string& gtiname, - const std::string& outfile) const -{ - // Save only if we have an event list - if (obs->eventtype() == "EventList") { - - // Set output FITS file event extension names - GFilename outname(outfile); - std::string outevt = evtname; - std::string outgti = gtiname; - if (outname.has_extname()) { - std::vector extnames = - gammalib::split(outname.extname(), ";"); - if (extnames.size() > 0) { - std::string extname = gammalib::strip_whitespace(extnames[0]); - if (!extname.empty()) { - outevt = extname; - } - } - if (extnames.size() > 1) { - std::string extname = gammalib::strip_whitespace(extnames[1]); - if (!extname.empty()) { - outgti = extname; - } - } - } - - // Create output FITS file - GFits outfits; - - // Write observation into FITS file - obs->write(outfits, outevt, outgti); - - // Copy all extensions other than evtname and gtiname extensions - // from the input to the output event list. The evtname and - // gtiname extensions are written by the save method, all others - // that may eventually be present have to be copied over - // explicitly. - GFits infits(infile); - for (int extno = 1; extno < infits.size(); ++extno) { - GFitsHDU* hdu = infits.at(extno); - if (hdu->extname() != evtname && - hdu->extname() != gtiname && - hdu->extname() != outevt && - hdu->extname() != outgti) { - outfits.append(*hdu); - } - } - - // Close input file - infits.close(); - - // Save file to disk and close it (we need both operations) - outfits.saveto(outname.url(), clobber()); - outfits.close(); - - } // endif: observation was unbinned - - // Return - return; -} diff --git a/src/ctselect/ctselect.hpp b/src/ctselect/ctselect.hpp index f41e209..1a1df7e 100644 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctselect - Data selection tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2016 by Juergen Knoedlseder * + * copyright (C) 2010-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -78,15 +78,8 @@ protected: std::string check_infile(const std::string& filename, const std::string& evtname) const; std::string set_outfile_name(const std::string& filename) const; - std::string get_gtiname(const std::string& filename, - const std::string& evtname) const; void save_fits(void); void save_xml(void); - void save_event_list(const GCTAObservation* obs, - const std::string& infile, - const std::string& evtname, - const std::string& gtiname, - const std::string& outfile) const; // User parameters std::string m_outobs; //!< Output event list or XML file diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index d71b914..2dbc5c9 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -1730,6 +1730,145 @@ std::vector ctool::cube_layer_usage(const GEbounds& cube_ebounds, /***********************************************************************//** + * @brief Save event list into FITS file + * + * @param[in] obs Pointer to CTA observation. + * @param[in] infile Input file name. + * @param[in] evtname Event extension name. + * @param[in] gtiname GTI extension name. + * @param[in] outfile Output file name. + * + * Saves an event list and the corresponding Good Time Intervals into a FITS + * file and copy all others extensions from the input file to the output + * file. + * + * If an extension name is specified in the @p outfile argument, the events + * and eventually also the Good Time Intervals will be extracted from the + * argument and used for writing the events. The format is + * + * [] + * + * where needs to be replaced by the name of the FITS file, + * and by the name of the events + * and Good Time Intervals extensions. For example + * + * myfits.fits[EVENTS1;GTI1] + * + * will write the selected events into the `EVENTS1` extension and the + * Good Time Intervals into the `GTI1` extension of the `myfits.fits` FITS + * file. If the Good Time Intervals extension name is skipped, e.g. + * + * myfits.fits[EVENTS1] + * + * the original extension name for the Good Time Intervals will be kept. + * Analogously, only the Good Time Intervals extension name can be changed + * by specifying + * + * myfits.fits[;GTI1] + * + * In none of the cases will the original events and Good Time Intervals be + * copied over to the output file. + ***************************************************************************/ +void ctool::save_event_list(const GCTAObservation* obs, + const std::string& infile, + const std::string& evtname, + const std::string& gtiname, + const std::string& outfile) const +{ + // Save only if we have an event list + if (obs->eventtype() == "EventList") { + + // Set output FITS file event extension names + GFilename outname(outfile); + std::string outevt = evtname; + std::string outgti = gtiname; + if (outname.has_extname()) { + std::vector extnames = + gammalib::split(outname.extname(), ";"); + if (extnames.size() > 0) { + std::string extname = gammalib::strip_whitespace(extnames[0]); + if (!extname.empty()) { + outevt = extname; + } + } + if (extnames.size() > 1) { + std::string extname = gammalib::strip_whitespace(extnames[1]); + if (!extname.empty()) { + outgti = extname; + } + } + } + + // Create output FITS file + GFits outfits; + + // Write observation into FITS file + obs->write(outfits, outevt, outgti); + + // Copy all extensions other than evtname and gtiname extensions + // from the input to the output event list. The evtname and + // gtiname extensions are written by the save method, all others + // that may eventually be present have to be copied over + // explicitly. + GFits infits(infile); + for (int extno = 1; extno < infits.size(); ++extno) { + GFitsHDU* hdu = infits.at(extno); + if (hdu->extname() != evtname && + hdu->extname() != gtiname && + hdu->extname() != outevt && + hdu->extname() != outgti) { + outfits.append(*hdu); + } + } + + // Close input file + infits.close(); + + // Save file to disk and close it (we need both operations) + outfits.saveto(outname.url(), clobber()); + outfits.close(); + + } // endif: observation was unbinned + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get Good Time Intervals extension name + * + * @param[in] filename Input file name. + * @param[in] evtname Events extension name. + * + * Extracts the Good Time Intervals extension name from the event file. We + * do this by loading the events and accessing the Good Time Intervals + * extension name using the GCTAEventList::gtiname() method. If the file name + * is empty, the method returns `GTI`. + ***************************************************************************/ +std::string ctool::get_gtiname(const std::string& filename, + const std::string& evtname) const +{ + // Initialise GTI name + std::string gtiname = gammalib::extname_gti; + + // Continue only if the filename is not empty + if (!filename.empty()) { + + // Load events + GCTAEventList events(filename+"["+evtname+"]"); + + // Get GTI name + gtiname = events.gtiname(); + + } + + // Return GTI name + return (gtiname); +} + + +/***********************************************************************//** * @brief Set warning string if there are too few energies * * @param[in] energies Energies. diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index 5a6dde5..6e0d30d 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctool - ctool base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -121,6 +121,13 @@ protected: std::vector cube_layer_usage(const GEbounds& cube_ebounds, const GEbounds& list_ebounds) const; bool is_valid_filename(const GFilename& filename) const; + std::string get_gtiname(const std::string& filename, + const std::string& evtname) const; + void save_event_list(const GCTAObservation* obs, + const std::string& infile, + const std::string& evtname, + const std::string& gtiname, + const std::string& outfile) const; // Protected warning strings std::string warn_too_few_energies(const GEnergies& energies) const; diff --git a/test/test_ctprob.py b/test/test_ctprob.py index 390add5..0ad8306 100755 --- a/test/test_ctprob.py +++ b/test/test_ctprob.py @@ -2,7 +2,7 @@ # ========================================================================== # This scripts performs unit tests for the ctprob tool. # -# Copyright (C) 2014-2016 Juergen Knoedlseder +# Copyright (C) 2017 by Leonardo Di Venere # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -73,7 +73,7 @@ class Test(test): ' outobs="ctprob_cmd1.fits"'+ \ ' inmodel='+self._model+ \ ' caldb='+self._caldb+' irf='+self._irf+ \ - ' edisp=yes'+ \ + ' edisp=no'+ \ ' logfile="ctprob_cmd1.log" chatter=1' # Check if execution was successful @@ -88,7 +88,7 @@ class Test(test): ' outobs="ctprob_cmd2.fits"'+ \ ' inmodel='+self._model+ \ ' caldb='+self._caldb+' irf='+self._irf+ \ - ' edisp=yes'+ \ + ' edisp=no'+ \ ' logfile="ctprob_cmd2.log" chatter=1' # Check if execution failed @@ -135,12 +135,13 @@ class Test(test): # Now set ctprob parameters prob['inobs'] = self._events prob['outobs'] = 'ctprob_py1.fits' - prob['inmodel'] = self._model - prob['caldb'] = self._caldb - prob['irf'] = self._irf - prob['edisp'] = True + prob['inmodel'] = self._model + prob['caldb'] = self._caldb + prob['irf'] = self._irf + prob['edisp'] = False prob['logfile'] = 'ctprob_py1.log' prob['chatter'] = 2 + # Run ctprob tool prob.logFileOpen() # Make sure we get a log file prob.run() @@ -149,15 +150,14 @@ class Test(test): # Check result file self._check_result_file('ctprob_py1.fits') - # Setup ctprob tool for an invalid event file prob = ctools.ctprob() prob['inobs'] = self._invalid_events prob['outobs'] = 'ctprob_py2.fits' - prob['inmodel'] = self._model - prob['caldb'] = self._caldb - prob['irf'] = self._irf - prob['edisp'] = True + prob['inmodel'] = self._model + prob['caldb'] = self._caldb + prob['irf'] = self._irf + prob['edisp'] = False prob['logfile'] = 'ctprob_py2.log' prob['chatter'] = 3 @@ -167,26 +167,24 @@ class Test(test): try: prob.execute() self.test_try_failure('Exception not thrown') - except ValueError: + except RuntimeError: self.test_try_success() - # Setup ctprob with a list of observations and test prefix prob = ctools.ctprob() - prob['inobs'] = self._datadir + '/obs_unbinned.xml' - prob['prefix'] = 'ctprob_py3_' + prob['inobs'] = self._datadir + '/obs_unbinned.xml' + prob['prefix'] = 'ctprob_py3_' prob['outobs'] = 'ctprob_py3_obs.xml' - prob['inmodel'] = self._model - prob['caldb'] = self._caldb - prob['irf'] = self._irf - prob['edisp'] = True + prob['inmodel'] = self._model + prob['caldb'] = self._caldb + prob['irf'] = self._irf + prob['edisp'] = True prob['logfile'] = 'ctprob_py3.log' - prob['chatter'] = 3 + prob['chatter'] = 4 # Run ctprob tool prob.logFileOpen() # Make sure we get a log file - prob.run() - prob.save() + prob.execute() # Check result file self._check_result_file('ctprob_py3_crab_events.fits') @@ -213,8 +211,8 @@ class Test(test): self._check_events(events, nevents=nevents) # Check that probability columns exist - colname1 = "PROB_Background" - colname2 = "PROB_Crab" + colname1 = 'PROB_Background' + colname2 = 'PROB_Crab' self._check_column(filename, colname1) self._check_column(filename, colname2) @@ -291,8 +289,8 @@ class Test(test): def _check_normalization(self, filename, colnames): """ - Check that the probability values in columns identified - by colnames are correctly normalized. + Check that the probability values in columns identified by colnames are + correctly normalized. Parameters ---------- @@ -304,26 +302,31 @@ class Test(test): # Check that the probability values are correctly normalized # Open fits table - gfits = gammalib.GFits(filename) + gfits = gammalib.GFits(filename) gtable = gfits[1] - ncols = gtable.ncols() - nevt = gtable.nrows() - srcs=[] + ncols = gtable.ncols() + nevt = gtable.nrows() + srcs = [] + # Loop over colnames for colname in colnames: - # self._check_column(filename, colnames) # Check if column exists + # Get the correct column corresponding to colname from the table for icol in range(ncols): gcol = gtable[icol] if gcol.name()==colname: break + # Read and save column content src1 = [] for ievt in range(nevt): src1.append(gcol[ievt]) srcs.extend(src1) + + # Compute sum tot = sum(srcs) - self.test_value(int(tot+0.5),nevt, 'Check that probability columns are normalized') + self.test_value(int(tot+0.5), nevt, 'Check that probability columns ' + 'are normalized') # Return return