From a35f9317adc1b83917b16433b4aaa3bd2f60114a Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 12 Aug 2016 00:19:40 +0200 Subject: [PATCH] Enhance ctools code coverage (#1828) --- ChangeLog | 3 ++- NEWS | 5 ++++- README.md | 2 +- pyext/Makefile.am | 1 + pyext/ctselect.i | 9 ++++----- src/ctbkgcube/ctbkgcube.cpp | 20 ++++++-------------- src/ctselect/ctselect.cpp | 167 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------- src/ctselect/ctselect.hpp | 13 ------------- src/cttsmap/cttsmap.cpp | 113 ++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------- src/ctulimit/ctulimit.cpp | 8 ++++---- test/data/invalid_event_list.fits.gz | Bin 0 -> 15054 bytes test/test_ctselect.py | 69 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------ test/test_cttsmap.py | 135 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------- test/test_ctulimit.py | 57 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 14 files changed, 387 insertions(+), 215 deletions(-) create mode 100644 test/data/invalid_event_list.fits.gz diff --git a/ChangeLog b/ChangeLog index 4295119..b8867aa 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2016-08-06 +2016-08-12 * Version 1.2.0 released ======================== + Remove ctselect::time_reference() method Add ctool and cscript pars() methods (#1838) Fix ctbin::save() bug for empty observation container (#1835) Fix exception when calling ctools clear() methods (#1835) diff --git a/NEWS b/NEWS index 3d2ad13..3540a36 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.2.0 -6 August 2016 +12 August 2016 Introduction @@ -38,6 +38,9 @@ ctselect - CTA event selection The clear() method now calls GApplication::clear() to conserve the properties of the tool (#1835). +The time_reference() method has been removed since the method is part of the +ctool base class, and exposed in Python as _time_reference(). + ctbin - CTA event binning ------------------------- diff --git a/README.md b/README.md index 357dae6..f29c028 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ctools information ================== -* Version: 1.2.0.dev1 (6 August 2016) +* Version: 1.2.0.dev1 (12 August 2016) * Author: Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu) * GammaLib dependency: 1.2.0.dev1 diff --git a/pyext/Makefile.am b/pyext/Makefile.am index 5ac7a53..770b777 100755 --- a/pyext/Makefile.am +++ b/pyext/Makefile.am @@ -64,6 +64,7 @@ testscripts = ctools/tests/__init__.py \ testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/crab_events.fits.gz \ $(top_srcdir)/test/data/crab_cntmap.fits.gz \ + $(top_srcdir)/test/data/invalid_event_list.fits.gz \ $(top_srcdir)/test/data/exclusion.reg # Install test scripts and test data. Note that we don't need to prepend diff --git a/pyext/ctselect.i b/pyext/ctselect.i index fa1eff9..2b5a236 100644 --- a/pyext/ctselect.i +++ b/pyext/ctselect.i @@ -44,11 +44,10 @@ public: virtual ~ctselect(void); // Methods - void clear(void); - void run(void); - void save(void); - const GObservations& obs(void) const; - const GTimeReference& time_reference(void) const; + void clear(void); + void run(void); + void save(void); + const GObservations& obs(void) const; }; diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp index e33af8b..ff5443f 100644 --- a/src/ctbkgcube/ctbkgcube.cpp +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -246,20 +246,12 @@ void ctbkgcube::run(void) // Flag removal bool remove = true; - // Do we have a CTA, HESS, MAGIC or VERITAS specific model? - if (m_bkgmdl[i]->is_valid("CTA", "") || - m_bkgmdl[i]->is_valid("HESS", "") || - m_bkgmdl[i]->is_valid("MAGIC", "") || - m_bkgmdl[i]->is_valid("VERITAS", "")) { - - // Do we have a background model? - if (dynamic_cast(m_bkgmdl[i]) != NULL) { - - // ... then keep model - remove = false; - - } // endif: had a background model - } // endif: had a CTA, HESS, MAGIC or VERITAS model + // If we have a data space model with "GCTA" classname then we + // have a CTA background model and we want to keep the model + if ((dynamic_cast(m_bkgmdl[i]) != NULL) && + (m_bkgmdl[i]->classname().substr(0,4) == "GCTA")) { + remove = false; + } // Log model removal or keeping std::string what = (remove) ? "Remove model" : "Keep model"; diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 86a2503..fff9a3c 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -291,8 +291,8 @@ void ctselect::run(void) // 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.length() > 0) { - throw GException::app_error(G_RUN, message); + if (!message.empty()) { + throw GException::invalid_value(G_RUN, message); } } @@ -322,8 +322,8 @@ void ctselect::run(void) // If we have a temporary file then check it if (!filename.empty()) { std::string message = check_infile(filename, m_evtname[i]); - if (message.length() > 0) { - throw GException::app_error(G_RUN, message); + if (!message.empty()) { + throw GException::invalid_value(G_RUN, message); } } @@ -611,8 +611,9 @@ void ctselect::select_events(GCTAObservation* obs, // Write header into logger log_header3(NORMAL, "Events selection"); - // Initialise selection string + // Initialise selection and addition strings std::string selection; + std::string add; // Initialise selection flags bool remove_all = false; @@ -679,6 +680,7 @@ void ctselect::select_events(GCTAObservation* obs, sprintf(cmin, "%.8f", tmin); sprintf(cmax, "%.8f", tmax); selection = "TIME >= "+std::string(cmin)+" && TIME <= "+std::string(cmax); + add = " && "; log_value(NORMAL, "Time range", gammalib::str(tmin)+" - "+gammalib::str(tmax)+" s"); @@ -687,11 +689,6 @@ void ctselect::select_events(GCTAObservation* obs, // Make energy selection if (m_select_energy) { - // If we have aready a selection than add an "&&" operator - if (selection.length() > 0) { - selection += " && "; - } - // Log the requested energy selection if (logNormal()) { log << gammalib::parformat("Selected energy range"); @@ -725,43 +722,40 @@ void ctselect::select_events(GCTAObservation* obs, char cmax[80]; sprintf(cmin, "%.8f", emin); sprintf(cmax, "%.8f", emax); - selection += "ENERGY >= "+std::string(cmin)+" && ENERGY <= "+ - std::string(cmax); + selection += add + "ENERGY >= " + std::string(cmin)+ + " && ENERGY <= " + std::string(cmax); + add = " && "; } } else if (emin > 0.0) { char cmin[80]; sprintf(cmin, "%.8f", emin); - selection += "ENERGY >= "+std::string(cmin); + selection += add + "ENERGY >= " + std::string(cmin); + add = " && "; } else if (emax > 0.0) { char cmax[80]; sprintf(cmax, "%.8f", emax); - selection += "ENERGY <= "+std::string(cmax); + selection += add + "ENERGY <= " + std::string(cmax); + add = " && "; } else { - log << "None" << std::endl; remove_all = true; } } // endif: made energy selection - // Make ROI selection + // Make RoI selection if (m_select_roi) { - // If we have aready a selection than add an "&&" operator - if (selection.length() > 0) { - selection += " && "; - } - - // Log the requested ROI + // Log the requested RoI log_value(NORMAL, "Requested RoI", "Centre(RA,DEC)=("+gammalib::str(ra)+", "+ gammalib::str(dec)+") deg, Radius="+gammalib::str(m_rad)+ " deg"); - // If we have already an ROI then make sure that the selected - // ROI overlaps with the existing ROI + // If we have already an RoI then make sure that the selected + // ROI overlaps with the existing RoI double roi_radius = list->roi().radius(); if (roi_radius > 0.0) { GSkyDir roi_centre = list->roi().centre().dir(); @@ -777,7 +771,7 @@ void ctselect::select_events(GCTAObservation* obs, } } - // Log the selected ROI + // Log the selected RoI if (logNormal()) { log << gammalib::parformat("Selected RoI"); if (rad <= 0.0) { @@ -790,29 +784,25 @@ void ctselect::select_events(GCTAObservation* obs, } } - // Format ROI selection + // Format RoI selection char cra[80]; char cdec[80]; char crad[80]; sprintf(cra, "%.6f", ra); sprintf(cdec, "%.6f", dec); sprintf(crad, "%.6f", rad); - selection += "ANGSEP("+std::string(cra)+"," + - std::string(cdec)+",RA,DEC) <= " + - std::string(crad); + selection += add + "ANGSEP("+std::string(cra)+"," + + std::string(cdec)+",RA,DEC) <= " + + std::string(crad); + add = " && "; } // endif: made ROI selection // Make an expression selection if (select_expr) { - // If we have aready a selection than add an "&&" operator - if (selection.length() > 0) { - selection += " && "; - } - // Append the expression - selection += "("+gammalib::strip_whitespace(m_expr)+")"; + selection += add + "("+gammalib::strip_whitespace(m_expr)+")"; } // endif: made expression selection @@ -875,11 +865,10 @@ void ctselect::select_events(GCTAObservation* obs, GCTAInstDir instdir; instdir.dir().radec_deg(ra, dec); list->roi(GCTARoi(instdir, rad)); - } // endif: Roi selection was performed + } + // ... otherwise restore old RoI information if it existed else if (old_roi.is_valid()) { - // Restore old RoI information in case no selection was performed - // and RoI was existing before list->roi(old_roi); } @@ -891,11 +880,10 @@ void ctselect::select_events(GCTAObservation* obs, GEbounds ebounds; ebounds.append(GEnergy(emin, "TeV"), GEnergy(emax, "TeV")); list->ebounds(ebounds); - } //endif: energy selection was performed + } + // ... otherwise restore old Ebounds if they existed else if (old_ebounds.size() > 0) { - // Restore old Ebounds in case no energy selection was performed - // and observation already had valid Ebounds list->ebounds(old_ebounds); } @@ -1342,65 +1330,60 @@ void ctselect::save_event_list(const GCTAObservation* obs, const std::string& gtiname, const std::string& outfile) const { - // Save only if observation is valid - if (obs != NULL) { - - // 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; - } + // 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; } } - - // 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); + if (extnames.size() > 1) { + std::string extname = gammalib::strip_whitespace(extnames[1]); + if (!extname.empty()) { + outgti = extname; } } + } - // Close input file - infits.close(); + // 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); + } + } - // Save file to disk and close it (we need both operations) - outfits.saveto(outname.url(), clobber()); - outfits.close(); + // Close input file + infits.close(); - } // endif: observation was unbinned + // Save file to disk and close it (we need both operations) + outfits.saveto(outname.url(), clobber()); + outfits.close(); - } // endif: observation was valid + } // endif: observation was unbinned // Return return; diff --git a/src/ctselect/ctselect.hpp b/src/ctselect/ctselect.hpp index e73d8c3..2f9cda5 100644 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -62,7 +62,6 @@ public: void run(void); void save(void); const GObservations& obs(void) const; - const GTimeReference& time_reference(void) const; protected: // Protected methods @@ -128,16 +127,4 @@ const GObservations& ctselect::obs(void) const return m_obs; } - -/***********************************************************************//** - * @brief Return CTA time reference - * - * @return CTA time reference - ***************************************************************************/ -inline -const GTimeReference& ctselect::time_reference(void) const -{ - return m_cta_ref; -} - #endif /* CTSELECT_HPP */ diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp index db730cb..664d983 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -234,7 +234,7 @@ void cttsmap::run(void) // Initialise optimizer GOptimizerLM* opt = (logExplicit()) ? new GOptimizerLM(&log) - : new GOptimizerLM(); + : new GOptimizerLM(); // Store initial models GModels models_orig = m_obs.models(); @@ -246,15 +246,15 @@ void cttsmap::run(void) models.remove(m_srcname); // Get likelihood of null hypothesis if not given before - if (m_logL0 == 0.0) { + if (m_logL0 == 0.0) { // Write header log_header1(TERSE, "Compute NULL Hypothesis for TS computation"); - // Compute likelihood without the test source - m_obs.models(models); - m_obs.optimize(*opt); - m_logL0 = -(opt->value()); + // Compute likelihood without the test source + m_obs.models(models); + m_obs.optimize(*opt); + m_logL0 = -(opt->value()); } @@ -272,11 +272,9 @@ void cttsmap::run(void) " at "+bincentre.print()); // Add test source at current bin position - if (m_testsource != NULL) { - (*m_testsource)["RA"].value(bincentre.ra_deg()); - (*m_testsource)["DEC"].value(bincentre.dec_deg()); - models.append(*m_testsource); - } + (*m_testsource)["RA"].value(bincentre.ra_deg()); + (*m_testsource)["DEC"].value(bincentre.dec_deg()); + models.append(*m_testsource); // Assign models to observations m_obs.models(models); @@ -292,16 +290,11 @@ void cttsmap::run(void) // Log information if (logExplicit()) { - log << " TS value ....: "; - log << ts << std::endl; + log_value(EXPLICIT, "TS value", ts); } - else if (logTerse()) { - log << "TS for bin number "; - log << i; - log << " at "; - log << bincentre.print(); - log << ": "; - log << ts << std::endl; + else if (logNormal()) { + log_value(NORMAL, "TS value (bin "+gammalib::str(i)+")", + gammalib::str(ts)+" ("+bincentre.print()+")"); } // Get test source model instance @@ -312,11 +305,9 @@ void cttsmap::run(void) m_tsmap(i) = ts; // Extract fitted test source parameters - if (m_testsource != NULL) { - for (int j = 0; j < m_mapnames.size(); ++j) { - m_maps[j](i) = (*testsource)[m_mapnames[j]].value(); - } - } + for (int j = 0; j < m_mapnames.size(); ++j) { + m_maps[j](i) = (*testsource)[m_mapnames[j]].value(); + } // Set status of bin to true m_statusmap(i) = 1.0; @@ -479,7 +470,8 @@ void cttsmap::copy_members(const cttsmap& app) m_statusmap = app.m_statusmap; // Clone protected members - m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() : NULL; + m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone() + : NULL; // Return return; @@ -491,6 +483,12 @@ void cttsmap::copy_members(const cttsmap& app) ***************************************************************************/ void cttsmap::free_members(void) { + // Free test source + if (m_testsource != NULL) { + delete m_testsource; + m_testsource = NULL; + } + // Return return; } @@ -519,26 +517,24 @@ void cttsmap::get_parameters(void) // Get name of test source m_srcname = (*this)["srcname"].string(); - // Get Model and store it as protected member + // Get model and store it as protected member. Note that the clone + // method will never return a NULL pointer, hence we do not need + // to check the pointer. m_testsource = m_obs.models()[m_srcname]->clone(); // Check if model has RA and DEC parameters which are necessary to // continue - if (m_testsource != NULL) { - if (!m_testsource->has_par("RA") || !m_testsource->has_par("DEC")) { - std::string msg = "Source \""+m_srcname+"\" has no \"RA\" and " - "\"DEC\" parameters. Only sources with \"RA\" " - " and \"DEC\" parameters can be used as test " - "sources."; - throw GException::invalid_value(G_GET_PARAMETERS, msg); - } + if (!m_testsource->has_par("RA") || !m_testsource->has_par("DEC")) { + std::string msg = "Source \""+m_srcname+"\" has no \"RA\" and " + "\"DEC\" parameters. Only sources with \"RA\" " + "and \"DEC\" parameters can be used as test " + "sources."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); } // Fix the spatial parameters of the test source - if (m_testsource != NULL) { - (*m_testsource)["RA"].fix(); - (*m_testsource)["DEC"].fix(); - } + (*m_testsource)["RA"].fix(); + (*m_testsource)["DEC"].fix(); // Create sky map based on task parameters GSkyMap map = create_map(m_obs); @@ -583,32 +579,27 @@ void cttsmap::init_maps(const GSkyMap& map) m_tsmap = GSkyMap(map); // Initialise map information - m_statusmap.clear(); - - // Create status map - m_statusmap = GSkyMap(map); - - // Initialise maps of free parameters - if (m_testsource != NULL) { + m_statusmap.clear(); - // Loop over all model parameters - for (int i = 0; i < m_testsource->size(); i++) { - - // Skip fixed parameters - if ((*m_testsource)[i].is_fixed()) { - continue; - } + // Create status map + m_statusmap = GSkyMap(map); - // Add sky map for free parameters. Note that push_back will - // create a copy. - m_maps.push_back(map); + // Loop over all model parameters + for (int i = 0; i < m_testsource->size(); i++) { + + // Skip fixed parameters + if ((*m_testsource)[i].is_fixed()) { + continue; + } - // Store parameter name - m_mapnames.push_back((*m_testsource)[i].name()); - - } // endfor: looped over all model parameters + // Add sky map for free parameters. Note that push_back will + // create a copy. + m_maps.push_back(map); - } // endif: test source was valid + // Store parameter name + m_mapnames.push_back((*m_testsource)[i].name()); + + } // endfor: looped over all model parameters // Return return; diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp index 887d2db..654e69b 100644 --- a/src/ctulimit/ctulimit.cpp +++ b/src/ctulimit/ctulimit.cpp @@ -618,10 +618,10 @@ void ctulimit::ulimit_bisection(const double& min, const double& max) // Change boundaries for further iteration if (eval_mid > 0.0) { - wrk_max = mid; + wrk_max = mid; // logL too large, new range = [wrk_min, mid] } - else if (eval_mid < 0.0) { - wrk_min = mid; + else { + wrk_min = mid; // logL too small, new range = [mid, wrk_max] } // Increment counter @@ -672,7 +672,7 @@ double ctulimit::evaluate(const double& value) else { std::string msg = "Value of parameter \""+m_model_par->name()+"\" " "outside of validity range requested. To omit " - "this error please enlarge the parameter range " + "this error please enlarge the parameter range " "in the model XML file."; throw GException::invalid_value(G_EVALUATE, msg); } diff --git a/test/data/invalid_event_list.fits.gz b/test/data/invalid_event_list.fits.gz new file mode 100644 index 0000000..c763170 Binary files /dev/null and b/test/data/invalid_event_list.fits.gz differ diff --git a/test/test_ctselect.py b/test/test_ctselect.py index b793e1d..34efcaa 100755 --- a/test/test_ctselect.py +++ b/test/test_ctselect.py @@ -39,6 +39,9 @@ class Test(test): # Call base class constructor test.__init__(self) + # Set test data + self._invalid_events = self._datadir + '/invalid_event_list.fits' + # Return return @@ -104,7 +107,7 @@ class Test(test): """ Test ctselect from Python """ - # Allocate ctselect + # Allocate empty ctselect tool select = ctools.ctselect() # Check that empty ctselect tool holds an empty observation @@ -122,7 +125,6 @@ class Test(test): select.clear() # Now set ctselect parameters - select = ctools.ctselect() select['inobs'] = self._events select['ra'] = 83.63 select['dec'] = 22.01 @@ -130,7 +132,7 @@ class Test(test): select['tmin'] = 0 select['tmax'] = 1800 select['emin'] = 0.1 - select['emax'] = 100 + select['emax'] = 100.0 select['outobs'] = 'ctselect_py1.fits' select['logfile'] = 'ctselect_py1.log' select['chatter'] = 2 @@ -166,7 +168,6 @@ class Test(test): # Check result file self._check_result_file('ctselect_py2.fits') - print('Step3') # Execute again the copy of ctselect tool again, now using the pointing # information in the input observation @@ -186,19 +187,75 @@ class Test(test): # Check that cleared ctselect tool holds no observations and events self._check_obs(cpy_select.obs(), nobs=0) + # Get mixed observation container + obs = self._obs_mixed() + + # Now setup ctselect tool from observation container. We make an + # event selection that will result in zero events in the file. + select = ctools.ctselect(obs) + select['ra'] = 83.63 + select['dec'] = 22.01 + select['rad'] = 3.0 + select['tmin'] = 0 + select['tmax'] = 1800 + select['emin'] = 120.0 + select['emax'] = 130.0 + select['expr'] = 'DETX == 0' + select['outobs'] = 'ctselect_py4.fits' + select['logfile'] = 'ctselect_py4.log' + select['chatter'] = 3 + + # Execute tool + select.logFileOpen() # Needed to get a new log file + select.execute() + + # Check result file + self._check_result_file('ctselect_py4.fits', nevents=0) + + # Setup ctselect tool for an invalid event file + select = ctools.ctselect() + select['inobs'] = self._invalid_events + select['ra'] = 83.63 + select['dec'] = 22.01 + select['rad'] = 3 + select['tmin'] = 0 + select['tmax'] = 1800 + select['emin'] = 0.1 + select['emax'] = 100.0 + select['outobs'] = 'ctselect_py5.fits' + select['logfile'] = 'ctselect_py5.log' + select['chatter'] = 3 + + # Execute tool + select.logFileOpen() # Needed to get a new log file + self.test_try('Test invalid event file') + try: + select.execute() + self.test_try_failure('Exception not thrown') + except ValueError: + self.test_try_success() + + # Return return # Check result file - def _check_result_file(self, filename): + def _check_result_file(self, filename, nevents=6127): """ Check result file + + Parameters + ---------- + filename : str + Event list file name + nevents : int, optional + Expected number of events """ # Open result file events = gammalib.GCTAEventList(filename) # Check event list - self._check_events(events, nevents=6127) + self._check_events(events, nevents=nevents) # Return return diff --git a/test/test_cttsmap.py b/test/test_cttsmap.py index e81ad49..5b3890d 100755 --- a/test/test_cttsmap.py +++ b/test/test_cttsmap.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# This scripts performs unit tests for the cttsmap tool. +# This scripts performs unit tests for the cttsmap tool # # Copyright (C) 2014-2016 Juergen Knoedlseder # @@ -18,6 +18,7 @@ # along with this program. If not, see . # # ========================================================================== +import os import gammalib import ctools from testing import test @@ -105,14 +106,31 @@ class Test(test): """ Test cttsmap from Python """ - # Set-up cttsmap + # Allocate empty cttsmap tool tsmap = ctools.cttsmap() + + # Check that empty cttsmap tool holds an empty observation and an + # empty TS map + self._check_obs(tsmap.obs(), nobs=0) + self._check_ts_map(tsmap.tsmap(), nx=0, ny=0) + + # Check that saving does not nothing + tsmap['outmap'] = 'cttsmap_py0.fits' + tsmap['logfile'] = 'cttsmap_py0.log' + tsmap.logFileOpen() + tsmap.save() + self.test_assert(not os.path.isfile('cttsmap_py0.fits'), + 'Check that no TS map has been created') + + # Check that clearing does not lead to an exception or segfault + tsmap.clear() + + # Now set cttsmap parameters tsmap['inobs'] = self._events tsmap['inmodel'] = self._model tsmap['srcname'] = 'Crab' tsmap['caldb'] = self._caldb tsmap['irf'] = self._irf - tsmap['outmap'] = 'cttsmap_py1.fits' tsmap['nxpix'] = 3 tsmap['nypix'] = 3 tsmap['binsz'] = 0.02 @@ -120,6 +138,7 @@ class Test(test): tsmap['proj'] = 'CAR' tsmap['xref'] = 83.63 tsmap['yref'] = 22.01 + tsmap['outmap'] = 'cttsmap_py1.fits' tsmap['logfile'] = 'cttsmap_py1.log' tsmap['chatter'] = 2 @@ -131,32 +150,116 @@ class Test(test): # Check result file self._check_result_file('cttsmap_py1.fits') + # Copy cttsmap tool + cpy_tsmap = tsmap.copy() + + # Check observation container and TS map of copy + self._check_obs(cpy_tsmap.obs()) + self._check_ts_map(cpy_tsmap.tsmap()) + + # Execute copy of cttsmap tool again, now with a higher chatter + # level than before + cpy_tsmap['binmin'] = 0 + cpy_tsmap['binmax'] = 1 + cpy_tsmap['outmap'] = 'cttsmap_py2.fits' + cpy_tsmap['logfile'] = 'cttsmap_py2.log' + cpy_tsmap['chatter'] = 3 + cpy_tsmap['publish'] = True + cpy_tsmap.logFileOpen() # Needed to get a new log file + cpy_tsmap.execute() + + # Check result file + self._check_result_file('cttsmap_py2.fits') + + # Now clear copy of cttsmap tool + cpy_tsmap.clear() + + # Check that the cleared copy has also cleared the observation + # container and TS map + self._check_obs(cpy_tsmap.obs(), nobs=0) + self._check_ts_map(cpy_tsmap.tsmap(), nx=0, ny=0) + # Return return # Check result file - def _check_result_file(self, filename): + def _check_result_file(self, filename, nx=3, ny=3): """ Check result file + + Parameters + ---------- + filename : str + TS map file name + nx : int, optional + Number of X pixels + ny : int, optional + Number of Y pixels """ # Open result file fits = gammalib.GFits(filename) # Get HDUs - ts = fits['Primary'] - prefactor = fits['Prefactor'] - index = fits['Index'] + ts = fits['Primary'] + pre = fits['Prefactor'] + index = fits['Index'] # Check dimensions - self.test_value(ts.naxis(), 2, 'Check for 2 dimensions') - self.test_value(ts.naxes(0), 3, 'Check for 3 pixels in X') - self.test_value(ts.naxes(1), 3, 'Check for 3 pixels in Y') - self.test_value(prefactor.naxis(), 2, 'Check for 2 dimensions') - self.test_value(prefactor.naxes(0), 3, 'Check for 3 pixels in X') - self.test_value(prefactor.naxes(1), 3, 'Check for 3 pixels in Y') - self.test_value(index.naxis(), 2, 'Check for 2 dimensions') - self.test_value(index.naxes(0), 3, 'Check for 3 pixels in X') - self.test_value(index.naxes(1), 3, 'Check for 3 pixels in Y') + self.test_value(ts.naxis(), 2, 'Check for TS map dimensions') + self.test_value(ts.naxes(0), nx, 'Check for TS map pixels in X') + self.test_value(ts.naxes(1), ny, 'Check for TS map pixels in Y') + self.test_value(pre.naxis(), 2, 'Check for Prefactor map dimensions') + self.test_value(pre.naxes(0), nx, 'Check for Prefactor map pixels in X') + self.test_value(pre.naxes(1), ny, 'Check for Prefactor map pixels in Y') + self.test_value(index.naxis(), 2, 'Check for Index map dimensions') + self.test_value(index.naxes(0), nx, 'Check for Index map pixels in X') + self.test_value(index.naxes(1), ny, 'Check for Index map pixels in Y') + + # Return + return + + # Check observation + def _check_obs(self, obs, nobs=1): + """ + Check observation + + Parameters + ---------- + obs : `~gammalib.GObservations` + Models + nobs : int, optional + Expected number of observations + """ + # Check number of observations + self.test_value(obs.size(), nobs, 'Check number of observations') + + # Return + return + + # Check TS map + def _check_ts_map(self, tsmap, nx=3, ny=3): + """ + Check TS map + + Parameters + ---------- + tsmap : `~gammalib.GSkyMap` + TS map + nx : int, optional + Number of X pixels + ny : int, optional + Number of Y pixels + """ + # Determine number of maps + if nx > 0 and ny > 0: + nmaps = 1 + else: + nmaps = 0 + + # Check dimensions + self.test_value(tsmap.nmaps(), nmaps, 'Check number of maps') + self.test_value(tsmap.nx(), nx, 'Check for number of X pixels') + self.test_value(tsmap.ny(), ny, 'Check for number of Y pixels') # Return return diff --git a/test/test_ctulimit.py b/test/test_ctulimit.py index 81ea41e..66705b9 100755 --- a/test/test_ctulimit.py +++ b/test/test_ctulimit.py @@ -118,7 +118,6 @@ class Test(test): ulimit.clear() # Now set ctulimit parameters - ulimit = ctools.ctulimit() ulimit['inobs'] = self._events ulimit['inmodel'] = self._model ulimit['srcname'] = 'Crab' @@ -140,6 +139,13 @@ class Test(test): self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16, 'Check upper limit on energy flux') + # Check obs() method + self.test_value(ulimit.obs().size(), 1, + 'Check number of observations in container') + + # Check opt() method + self.test_value(ulimit.opt().status(), 0, 'Check optimizer status') + # Copy ctulimit tool cpy_ulimit = ulimit.copy() @@ -166,6 +172,9 @@ class Test(test): self.test_value(cpy_ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16, 'Check upper limit on energy flux') + # Save copy of observation container for later + obs = cpy_ulimit.obs().copy() + # Now clear copy of ctulimit tool cpy_ulimit.clear() @@ -177,6 +186,52 @@ class Test(test): self.test_value(cpy_ulimit.eflux_ulimit(), 0.0, 1.0e-16, 'Check upper limit on energy flux') + # Now set ctulimit tool using the observation container from last + # time. This should avoid the necessity to recompute the maximum + # likelihood + ulimit = ctools.ctulimit(obs) + #ulimit['inobs'] = self._events + #ulimit['inmodel'] = self._model + ulimit['srcname'] = 'Crab' + #ulimit['caldb'] = self._caldb + #ulimit['irf'] = self._irf + ulimit['logfile'] = 'ctulimit_py3.log' + ulimit['chatter'] = 4 + + # Execute ctulimit tool + ulimit.logFileOpen() # Needed to get a new log file + ulimit.execute() + + # Check results + self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23, + 'Check differential upper limit') + self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17, + 'ulimit upper limit on photon flux') + self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16, + 'Check upper limit on energy flux') + + # Test invalid model name + self.test_try('Test invalid model name') + try: + ulimit['srcname'] = 'Weihnachtsstern' + ulimit['logfile'] = 'ctulimit_py4.log' + ulimit.logFileOpen() + ulimit.execute() + self.test_try_failure('Exception not thrown') + except ValueError: + self.test_try_success() + + # Test specification of background model + ulimit['srcname'] = 'Background' + ulimit['logfile'] = 'ctulimit_py5.log' + ulimit.logFileOpen() + self.test_try('Test invalid model name') + try: + ulimit.execute() + self.test_try_failure('Exception not thrown') + except ValueError: + self.test_try_success() + # Return return -- libgit2 0.21.4