Iftach Sadeh / ctools

Commit 6985d99f88789be96a27619f22cf2da28dcaab8e

Authored by ldivenere
1 parent a21ba927

developed ctprob app - works properly

Inline Side-by-side
Showing 5 changed files with 708 additions and 471 deletions
configure.ac
  View file @6985d99
... ... @@ -1035,6 +1035,7 @@ AC_CONFIG_FILES([Makefile
1035 1035 src/ctbutterfly/Makefile
1036 1036 src/ctulimit/Makefile
1037 1037 src/cterror/Makefile
  1038 + src/ctprob/Makefile
1038 1039 src/support/Makefile
1039 1040 src/ctools-setup
1040 1041 test/Makefile
... ...
src/ctprob/ctprob.cpp
  View file @6985d99
... ... @@ -28,20 +28,18 @@
28 28 #ifdef HAVE_CONFIG_H
29 29 #include <config.h>
30 30 #endif
31   -#include <cstdio>
  31 +#include <cstdlib>
32 32 #include "ctprob.hpp"
33 33 #include "GTools.hpp"
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
36   -#define G_GET_PARAMETERS "ctprob::get_parameters()"
37   -#define G_FILL_CUBE "ctprob::fill_cube(GCTAObservation*)"
  36 +#define G_RUN "ctprob::run()"
  37 +#define G_GET_OBS "ctprob::get_obs()"
38 38  
39 39 /* __ Debug definitions __________________________________________________ */
40 40  
41 41 /* __ Coding definitions _________________________________________________ */
42   -
43   -/* __ Constants __________________________________________________________ */
44   -const GEnergy g_energy_margin(1.0e-12, "TeV");
  42 +//#define G_USE_MKSTEMP //!< Use mkstemp for temporary filename
45 43  
46 44  
47 45 /*==========================================================================
... ... @@ -68,11 +66,11 @@ ctprob::ctprob(void) : ctobservation(CTPROB_NAME, CTPROB_VERSION)
68 66 *
69 67 * param[in] obs Observation container.
70 68 *
71   - * This method creates an instance of the class by copying an existing
72   - * observations container.
  69 + * Creates an instance of the class that is initialised using the information
  70 + * provided in an observation container.
73 71 ***************************************************************************/
74 72 ctprob::ctprob(const GObservations& obs) :
75   - ctobservation(CTPROB_NAME, CTPROB_VERSION, obs)
  73 + ctobservation(CTPROB_NAME, CTPROB_VERSION, obs)
76 74 {
77 75 // Initialise members
78 76 init_members();
... ... @@ -83,14 +81,15 @@ ctprob::ctprob(const GObservations&amp; obs) :
83 81  
84 82  
85 83  
  84 +
86 85 /***********************************************************************//**
87 86 * @brief Command line constructor
88 87 *
89 88 * @param[in] argc Number of arguments in command line.
90 89 * @param[in] argv Array of command line arguments.
91 90 ***************************************************************************/
92   -ctprob::ctprob(int argc, char *argv[]) :
93   - ctobservation(CTPROB_NAME, CTPROB_VERSION, argc, argv)
  91 +ctprob::ctprob(int argc, char *argv[]) :
  92 + ctobservation(CTPROB_NAME, CTPROB_VERSION, argc, argv)
94 93 {
95 94 // Initialise members
96 95 init_members();
... ... @@ -202,41 +201,45 @@ void ctprob::clear(void)
202 201  
203 202  
204 203 /***********************************************************************//**
205   - * @brief Generate the model map(s)
  204 + * @brief Append probability columns event data
206 205 *
207   - * This method reads the task parameters from the parfile, sets up the
208   - * observation container, loops over all CTA observations in the container
209   - * and generates a model map for each CTA observation.
  206 + * 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.
210 211 ***************************************************************************/
211 212 void ctprob::run(void)
212 213 {
213   - // If we're in debug mode then all output is also dumped on the screen
  214 + // Switch screen logging on in debug mode
214 215 if (logDebug()) {
215 216 log.cout(true);
216 217 }
217 218  
218   - // Get task parameters
  219 + // Get parameters
219 220 get_parameters();
220 221  
221   - // Set energy dispersion flags of all CTA observations and save old
222   - // values in save_edisp vector
223   - std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
224   -
225 222 // Write input observation container into logger
226 223 log_observations(NORMAL, m_obs, "Input observation");
227 224  
228   - // Write input model container into logger
229   - log_models(NORMAL, m_obs.models(), "Input model");
  225 + // Write header into logger
  226 + log_header1(TERSE, "Event selection");
230 227  
231   - // Write header
232   - log_header1(TERSE, "Generate model cube");
  228 + // Initialise counters
  229 + int n_observations = 0;
233 230  
234   - // Loop over all observations in the container
  231 + // Loop over all observation in the container
235 232 for (int i = 0; i < m_obs.size(); ++i) {
236 233  
237 234 // Write header for the current observation
238 235 log_header3(TERSE, get_obs_header(m_obs[i]));
239 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 +
240 243 // Get CTA observation
241 244 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
242 245  
... ... @@ -248,13 +251,6 @@ void ctprob::run(void)
248 251 continue;
249 252 }
250 253  
251   - // Fill cube and leave loop if we are binned mode (meaning we
252   - // only have one binned observation)
253   - if (m_binned) {
254   - fill_cube(obs);
255   - break;
256   - }
257   -
258 254 // Skip observation if we have a binned observation
259 255 if (obs->eventtype() == "CountsCube") {
260 256 std::string msg = " Skipping binned "+m_obs[i]->instrument()+
... ... @@ -263,25 +259,95 @@ void ctprob::run(void)
263 259 continue;
264 260 }
265 261  
266   - // Fill the cube
267   - fill_cube(obs);
  262 + // Increment counter
  263 + n_observations++;
  264 +
  265 + // Save event file name (for possible saving)
  266 + m_infiles[i] = obs->eventfile();
268 267  
269   - // Dispose events to free memory if event file exists on disk
270   - if (obs->eventfile().length() > 0 && obs->eventfile().exists()) {
271   - obs->dispose_events();
  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 272 }
  273 + m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]);
273 274  
274   - } // endfor: looped over observations
  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 selection skipped.");
  284 + continue;
  285 + }
275 286  
276   - // Write model cube into header
277   - log_header1(NORMAL, "Model cube");
278   - log_string(NORMAL, m_cube.print());
  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 + }
279 294  
280   - // Restore energy dispersion flags of all CTA observations
281   - restore_edisp(m_obs, save_edisp);
  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 + }
282 317  
283   - // Optionally publish model cube
284   - if (m_publish) {
  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
  327 + 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 +
  335 + // Remove temporary file
  336 + std::remove(filename.c_str());
  337 +
  338 + } // endfor: looped over all observations
  339 +
  340 + // If more than a single observation has been handled then make sure that
  341 + // an XML file will be used for storage
  342 + if (n_observations > 1) {
  343 + m_use_xml = true;
  344 + }
  345 +
  346 + // Write observation(s) into logger
  347 + log_observations(NORMAL, m_obs, "Output observation");
  348 +
  349 + // Optionally publish event list(s)
  350 + if ((*this)["publish"].boolean()) {
285 351 publish();
286 352 }
287 353  
... ... @@ -291,30 +357,37 @@ void ctprob::run(void)
291 357  
292 358  
293 359 /***********************************************************************//**
294   - * @brief Save model cube
  360 + * @brief Save the selected event list(s)
  361 + *
  362 + * This method saves the selected event list(s) into FITS file(s). There are
  363 + * two modes, depending on the m_use_xml flag.
295 364 *
296   - * Saves the model cube into a FITS file specified using the "outfile"
297   - * task parameter.
  365 + * If m_use_xml is true, all selected event list(s) will be saved into FITS
  366 + * files, where the output filenames are constructued from the input
  367 + * filenames by prepending the m_prefix string to name. Any path information
  368 + * will be stripped form the input name, hence event files will be written
  369 + * into the local working directory (unless some path information is present
  370 + * in the prefix). In addition, an XML file will be created that gathers
  371 + * the filename information for the selected event list(s). If an XML file
  372 + * was present on input, all metadata information will be copied from this
  373 + * input file.
  374 + *
  375 + * If m_use_xml is false, the selected event list will be saved into a FITS
  376 + * file.
298 377 ***************************************************************************/
299 378 void ctprob::save(void)
300 379 {
301   - // Write header
302   - log_header1(TERSE, "Save model cube");
303   -
304   - // Get model cube filename
305   - m_outcube = (*this)["outcube"].filename();
306   -
307   - // Save only if filename is non-empty and cube has some size
308   - if (!m_outcube.is_empty() && m_cube.size() > 0) {
309   - m_cube.save(m_outcube, clobber());
  380 + // Write header into logger
  381 + log_header1(TERSE, gammalib::number("Save event list", m_obs.size()));
  382 + // Case A: Save event file(s) and XML metadata information
  383 + if (m_use_xml) {
  384 + save_xml();
310 385 }
311 386  
312   - // Write into logger what has been done
313   - std::string fname = (m_outcube.is_empty()) ? "NONE" : m_outcube.url();
314   - if (m_cube.size() == 0) {
315   - fname.append(" (cube is empty, no file created)");
  387 + // Case B: Save event file as FITS file
  388 + else {
  389 + save_fits();
316 390 }
317   - log_value(NORMAL, "Model cube file", fname);
318 391  
319 392 // Return
320 393 return;
... ... @@ -322,51 +395,53 @@ void ctprob::save(void)
322 395  
323 396  
324 397 /***********************************************************************//**
325   - * @brief Publish model cube
  398 + * @brief Publish event lists
326 399 *
327   - * @param[in] name Model cube name.
  400 + * @param[in] name Event list name.
328 401 ***************************************************************************/
329 402 void ctprob::publish(const std::string& name)
330 403 {
331 404 // Write header into logger
332   - log_header1(TERSE, "Publish model cube");
  405 + log_header1(TERSE, gammalib::number("Publish event list", m_obs.size()));
333 406  
334   - // Set default name is user name is empty
335   - std::string user_name(name);
336   - if (user_name.empty()) {
337   - user_name = CTPROB_NAME;
338   - }
  407 + // Loop over all observation in the container
  408 + for (int i = 0; i < m_obs.size(); ++i) {
339 409  
340   - // Log filename
341   - log_value(NORMAL, "Model cube name", user_name);
  410 + // Get CTA observation
  411 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
342 412  
343   - // Publish model cube
344   - m_cube.counts().publish(user_name);
  413 + // Handle only CTA observations
  414 + if (obs != NULL) {
345 415  
346   - // Return
347   - return;
348   -}
  416 + // Continue only if there is an event list
  417 + if (obs->events()->size() != 0) {
349 418  
  419 + // Set default name if user name is empty
  420 + std::string user_name(name);
  421 + if (user_name.empty()) {
  422 + user_name = CTPROB_NAME;
  423 + }
350 424  
351   -/***********************************************************************//**
352   - * @brief Set model cube
353   - *
354   - * @param[in] cube Model cube.
355   - *
356   - * Set model cube and set all cube bins to zero.
357   - ***************************************************************************/
358   -void ctprob::cube(const GCTAEventCube& cube)
359   -{
360   - // Set cube
361   - m_cube = cube;
  425 + // If there are several event lists then add an index
  426 + if (m_use_xml) {
  427 + user_name += gammalib::str(i);
  428 + }
362 429  
363   - // Set all cube bins to zero
364   - for (int i = 0; i < m_cube.size(); ++i) {
365   - m_cube[i]->counts(0.0);
366   - }
  430 + // Write event list name into logger
  431 + log_value(NORMAL, "Event list name", user_name);
  432 +
  433 + // Write events into in-memory FITS file
  434 + GFits fits;
  435 + obs->write(fits);
  436 +
  437 + // Publish
  438 + fits.publish(gammalib::extname_cta_events, user_name);
367 439  
368   - // Signal that cube has been set
369   - m_has_cube = true;
  440 + } // endif: there were events
  441 +
  442 + } // endif: observation was a CTA observation
  443 +
  444 + } // endfor: looped over observations
370 445  
371 446 // Return
372 447 return;
... ... @@ -384,18 +459,28 @@ void ctprob::cube(const GCTAEventCube&amp; cube)
384 459 ***************************************************************************/
385 460 void ctprob::init_members(void)
386 461 {
387   - // Initialise members
388   - m_outcube.clear();
389   - m_apply_edisp = false;
390   - m_publish = false;
391   - m_chatter = static_cast<GChatter>(2);
  462 + // Initialise parameters
  463 + m_outobs.clear();
  464 + m_prefix.clear();
  465 + /* m_usepnt = false;
  466 + m_roi.clear();
  467 + m_tmin = 0.0;
  468 + m_tmax = 0.0;
  469 + m_emin = 0.0;
  470 + m_emax = 0.0;
  471 + m_expr.clear();
  472 + m_usethres.clear();*/
  473 + m_chatter = static_cast<GChatter>(2);
392 474  
393 475 // Initialise protected members
394   - m_cube.clear();
395   - m_gti.clear();
396   - m_has_cube = false;
397   - m_append_cube = false;
398   - m_binned = false;
  476 + m_infiles.clear();
  477 + m_evtname.clear();
  478 + m_gtiname.clear();
  479 + // m_timemin.clear();
  480 + //m_timemax.clear();
  481 + //m_select_energy = true;
  482 + //m_select_roi = true;
  483 + //m_select_time = true;
399 484  
400 485 // Return
401 486 return;
... ... @@ -409,18 +494,28 @@ void ctprob::init_members(void)
409 494 ***************************************************************************/
410 495 void ctprob::copy_members(const ctprob& app)
411 496 {
412   - // Copy attributes
413   - m_outcube = app.m_outcube;
414   - m_apply_edisp = app.m_apply_edisp;
415   - m_publish = app.m_publish;
416   - m_chatter = app.m_chatter;
  497 + // Copy parameters
  498 + m_outobs = app.m_outobs;
  499 + m_prefix = app.m_prefix;
  500 + /* m_usepnt = app.m_usepnt;
  501 + m_roi = app.m_roi;
  502 + m_tmin = app.m_tmin;
  503 + m_tmax = app.m_tmax;
  504 + m_emin = app.m_emin;
  505 + m_emax = app.m_emax;
  506 + m_expr = app.m_expr;
  507 + m_usethres = app.m_usethres;*/
  508 + m_chatter = app.m_chatter;
417 509  
418 510 // Copy protected members
419   - m_cube = app.m_cube;
420   - m_gti = app.m_gti;
421   - m_has_cube = app.m_has_cube;
422   - m_append_cube = app.m_append_cube;
423   - m_binned = app.m_binned;
  511 + m_infiles = app.m_infiles;
  512 + m_evtname = app.m_evtname;
  513 + m_gtiname = app.m_gtiname;
  514 + /* m_timemin = app.m_timemin;
  515 + m_timemax = app.m_timemax;
  516 + m_select_energy = app.m_select_energy;
  517 + m_select_roi = app.m_select_roi;
  518 + m_select_time = app.m_select_time;*/
424 519  
425 520 // Return
426 521 return;
... ... @@ -440,58 +535,31 @@ void ctprob::free_members(void)
440 535 /***********************************************************************//**
441 536 * @brief Get application parameters
442 537 *
443   - * Get all task parameters from parameter file or (if required) by querying
444   - * the user. The parameters are read in the correct order.
  538 + * Get all user parameters from parameter file or (if required) by querying
  539 + * the user. Times are assumed to be in the native CTA MJD format.
  540 + *
  541 + * This method also loads observations if no observations are yet allocated.
  542 + * Observations are either loaded from a single CTA even list, or from a
  543 + * XML file using the metadata information that is stored in that file.
445 544 ***************************************************************************/
446 545 void ctprob::get_parameters(void)
447 546 {
448   - // Reset cube append flag
449   - m_append_cube = false;
450   -
451   - // If there are no observations in container then load them via user
452   - // parameters.
  547 + // Initialise selection flags
  548 + /* m_select_energy = true;
  549 + m_select_roi = true;
  550 + m_select_time = true;
  551 + */
  552 + // Setup observations from "inobs" parameter. Do not request response
  553 + // information and do not accept counts cubes.
453 554 if (m_obs.size() == 0) {
454 555 get_obs();
455 556 }
456   -
457 557 // ... otherwise add response information and energy boundaries in case
458 558 // that they are missing
459 559 else {
460 560 setup_observations(m_obs);
461 561 }
462   -
463   - // If we have now excactly one CTA observation (but no cube has yet been
464   - // appended to the observation) then check whether this observation
465   - // is a binned observation, and if yes, extract the counts cube for
466   - // model generation
467   - if ((m_obs.size() == 1) && (m_append_cube == false)) {
468   -
469   - // Get CTA observation
470   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
471   -
472   - // Continue only if observation is a CTA observation
473   - if (obs != NULL) {
474   -
475   - // Check for binned observation
476   - if (obs->eventtype() == "CountsCube") {
477   -
478   - // Set cube from binned observation
479   - GCTAEventCube* evtcube = dynamic_cast<GCTAEventCube*>
480   - (const_cast<GEvents*>(obs->events()));
481   -
482   - cube(*evtcube);
483   -
484   - // Signal that cube has been set
485   - m_has_cube = true;
486   -
487   - // Signal that we are in binned mode
488   - m_binned = true;
489   -
490   - } // endif: observation was binned
491   -
492   - } // endif: observation was CTA
493   -
494   - } // endif: had exactly one observation
  562 + //setup_observations(m_obs, false, true, false);
495 563  
496 564 // Read model definition file if required
497 565 if (m_obs.models().size() == 0) {
... ... @@ -507,66 +575,18 @@ void ctprob::get_parameters(void)
507 575 // Get energy dispersion flag parameters
508 576 m_apply_edisp = (*this)["edisp"].boolean();
509 577  
510   - // If we do not have yet a counts cube for model computation then check
511   - // whether we should read it from the "incube" parameter or whether we
512   - // should create it from scratch using the task parameters
513   - if (!m_has_cube) {
514   -
515   - // Read cube definition file
516   - std::string incube = (*this)["incube"].filename();
517   -
518   - // If the cube filename is valid the load the cube and set all cube
519   - // bins to zero
520   - if (is_valid_filename(incube)) {
521   -
522   - // Load cube from given file
523   - m_cube.load(incube);
524   -
525   - // Set all cube bins to zero
526   - for (int i = 0; i < m_cube.size(); ++i) {
527   - m_cube[i]->counts(0.0);
528   - }
529   -
530   - } // endif: cube filename was valid
531   -
532   - // ... otherwise create a cube from the user parameters
533   - else {
534   - m_cube = create_cube(m_obs);
535   - }
536   -
537   - // Signal that cube has been set
538   - m_has_cube = true;
539   -
540   - } // endif: we had no cube yet
541 578  
542 579 // Get remaining parameters
543 580 m_publish = (*this)["publish"].boolean();
544 581 m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
545 582  
546   - // Read optionally output cube filenames
  583 + // Optionally read ahead parameters so that they get correctly
  584 + // dumped into the log file
547 585 if (read_ahead()) {
548   - m_outcube = (*this)["outcube"].filename();
  586 + m_outobs = (*this)["outobs"].filename();
  587 + m_prefix = (*this)["prefix"].string();
549 588 }
550 589  
551   - // If cube should be appended to first observation then do that now.
552   - // This is a kluge that makes sure that the cube is passed as part
553   - // of the observation in case that a cube response is used. The kluge
554   - // is needed because the GCTACubeSourceDiffuse::set method needs to
555   - // get the full event cube from the observation. It is also at this
556   - // step that the GTI, which may just be a dummy GTI when create_cube()
557   - // has been used, will be set.
558   - if (m_append_cube) {
559   -
560   - //TODO: Check that energy boundaries are compatible
561   -
562   - // Attach GTI of observations to model cube
563   - m_cube.gti(m_obs[0]->events()->gti());
564   -
565   - // Attach model cube to observations
566   - m_obs[0]->events(m_cube);
567   -
568   - } // endif: cube was scheduled for appending
569   -
570 590 // Write parameters into logger
571 591 log_parameters(TERSE);
572 592  
... ... @@ -578,89 +598,21 @@ void ctprob::get_parameters(void)
578 598 /***********************************************************************//**
579 599 * @brief Get observation container
580 600 *
  601 + * @exception GException::invalid_value
  602 + * Invalid FITS file.
  603 + *
581 604 * Get an observation container according to the user parameters. The method
582 605 * supports loading of a individual FITS file or an observation definition
583 606 * file in XML format.
584   - *
585   - * If the input filename is empty, the method checks for the existence of the
586   - * "expcube", "psfcube" and "bkgcube" parameters. If file names have been
587   - * specified, the method loads the files and creates a dummy events cube that
588   - * is appended to the observation container.
589   - *
590   - * If no file names are specified for the "expcube", "psfcube" or "bkgcube"
591   - * parameters, the method reads the necessary parameters to build a CTA
592   - * observation from scratch.
593   - *
594   - * The method sets m_append_cube = true and m_binned = true in case that
595   - * a stacked observation is requested (as detected by the presence of the
596   - * "expcube", "psfcube", and "bkgcube" parameters). In that case, it appended
597   - * a dummy event cube to the observation.
598   - *
599   - * @todo Support stacked energy dispersion
600 607 ***************************************************************************/
601 608 void ctprob::get_obs(void)
602 609 {
603 610 // Get the filename from the input parameters
604 611 std::string filename = (*this)["inobs"].filename();
605 612  
606   - // If no observation definition file has been specified then read all
607   - // parameters that are necessary to create an observation from scratch
  613 + // If no observation definition file throw an exception
608 614 if (!is_valid_filename(filename)) {
609   -
610   - // Get response cube filenames
611   - std::string expcube = (*this)["expcube"].filename();
612   - std::string psfcube = (*this)["psfcube"].filename();
613   - std::string bkgcube = (*this)["bkgcube"].filename();
614   -
615   - // If the filenames are valid then build an observation from cube
616   - // response information
617   - if (is_valid_filename(expcube) && is_valid_filename(psfcube) &&
618   - is_valid_filename(bkgcube)) {
619   -
620   - // Get exposure, PSF and background cubes
621   - GCTACubeExposure exposure(expcube);
622   - GCTACubePsf psf(psfcube);
623   - GCTACubeBackground background(bkgcube);
624   -
625   - // Create energy boundaries
626   - GEbounds ebounds = create_ebounds();
627   -
628   - // Create dummy sky map cube
629   - GSkyMap map("CAR","GAL",0.0,0.0,1.0,1.0,1,1,ebounds.size());
630   -
631   - // Create event cube
632   - GCTAEventCube cube(map, ebounds, exposure.gti());
633   -
634   - // Create CTA observation
635   - GCTAObservation cta;
636   - cta.events(cube);
637   - cta.response(exposure, psf, background);
638   -
639   - // Append observation to container
640   - m_obs.append(cta);
641   -
642   - // Signal that we are in binned mode
643   - m_binned = true;
644   -
645   - // Signal that we appended a cube
646   - m_append_cube = true;
647   -
648   - } // endif: cube response information was available
649   -
650   - // ... otherwise build an observation from IRF response information
651   - else {
652   -
653   - // Create CTA observation
654   - GCTAObservation cta = create_cta_obs();
655   -
656   - // Set response
657   - set_obs_response(&cta);
658   -
659   - // Append observation to container
660   - m_obs.append(cta);
661   -
662   - }
663   -
  615 + throw GException::invalid_value(G_GET_OBS, "Input event list is not valid.");
664 616 } // endif: filename was not valid
665 617  
666 618 // ... otherwise we have a file name
... ... @@ -713,99 +665,426 @@ void ctprob::get_obs(void)
713 665 }
714 666  
715 667  
  668 +
716 669 /***********************************************************************//**
717   - * @brief Fill model into model cube
  670 + * @brief Evaluate probability for events
718 671 *
719   - * @param[in] obs CTA observation.
  672 + * @param[in,out] obs CTA observation.
720 673 *
721   - * Adds the expected number of events for a given observation to the events
722   - * that are already found in the model cube. The method also updates the
723   - * GTI of the model cube so that cube GTI is a list of the GTIs of all
724   - * observations that were used to generate the model cube.
  674 + *
  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 differential expected counts for the event direction and
  678 + * energy and normalizing it to the total differential expected counts for
  679 + * given model.
  680 + * observation is replaced by selected event list read from the FITS file.
  681 + *
  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_".
725 685 ***************************************************************************/
726   -void ctprob::fill_cube(const GCTAObservation* obs)
  686 +void ctprob::evaluate_probability(GCTAObservation* obs)
727 687 {
728   - // Get references to GTI and energy boundaries for the event list
729   - const GGti& gti = obs->events()->gti();
730   - const GEbounds& obs_ebounds = obs->ebounds();
  688 + // Write header into logger
  689 + log_header3(NORMAL, "Evaluating probability for events");
731 690  
732   - // Get cube energy boundaries
733   - const GEbounds& cube_ebounds = m_cube.ebounds();
  691 + GCTAEventList* evts = static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
  692 + const GEvent* evt;
734 693  
735   - // Get counts cube usage flags
736   - std::vector<bool> usage = cube_layer_usage(cube_ebounds, obs_ebounds);
  694 + double total=0.;
  695 + double value;
737 696  
738   - // Initialise empty, invalid RoI
739   - GCTARoi roi;
  697 + std::vector<GFitsTableFloatCol*> columns;
740 698  
741   - // Retrieve RoI in case we have an unbinned observation
742   - if (obs->eventtype() == "EventList") {
743   - roi = obs->roi();
  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);
744 704 }
745 705  
746   - // Initialise statistics
747   - double sum = 0.0;
748   - int num_outside_ebds = 0;
749   - int num_outside_roi = 0;
  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 + //Loop over models
  712 + for (int j = 0; j<m_obs.models().size(); ++j) {
  713 + value = m_obs.models()[j]->eval( *evt, *obs);
  714 + values.push_back(value);
  715 + total += value;
  716 + }
  717 + for (int j = 0; j<m_obs.models().size(); ++j) {
  718 + values[j] /= total;
  719 + (*(columns[j]))(i) = values[j];
  720 + }
  721 + }
750 722  
751   - // Setup cube GTIs for this observation
752   - m_cube.gti(obs->events()->gti());
  723 + //append columns to observation event list
  724 + for (int j = 0; j<m_obs.models().size(); ++j) {
  725 + evts->append_column(*columns[j]);
  726 + }
753 727  
754   - // Loop over all cube bins
755   - for (int i = 0; i < m_cube.size(); ++i) {
  728 + // Return
  729 + return;
  730 +}
  731 +
  732 +
  733 +
  734 +/***********************************************************************//**
  735 + * @brief Check input filename
  736 + *
  737 + * @param[in] filename File name.
  738 + * @param[in] evtname Event extension name.
  739 + *
  740 + * This method checks if the input FITS file is correct.
  741 + ***************************************************************************/
  742 +std::string ctprob::check_infile(const std::string& filename,
  743 + const std::string& evtname) const
  744 +{
  745 + // Initialise message string
  746 + std::string message = "";
  747 +
  748 + // Open FITS file
  749 + GFits fits(filename);
  750 +
  751 + // Check for existence of events extensions
  752 + if (!fits.contains(evtname)) {
  753 + message = "No \""+evtname+"\" extension found in input file \""+
  754 + filename + "\".";
  755 + }
  756 +
  757 + // ... otherwise check column names
  758 + else {
  759 +
  760 + // Get pointer to FITS table
  761 + GFitsTable* table = fits.table(evtname);
  762 +
  763 + // Initialise list of missing columns
  764 + std::vector<std::string> missing;
  765 +
  766 + // Check for existence of TIME column
  767 + if (!table->contains("TIME")) {
  768 + missing.push_back("TIME");
  769 + }
  770 +
  771 + // Check for existence of ENERGY column
  772 + if (!table->contains("ENERGY")) {
  773 + missing.push_back("ENERGY");
  774 + }
  775 +
  776 + // Check for existence of RA column
  777 + if (!table->contains("RA")) {
  778 + missing.push_back("RA");
  779 + }
  780 +
  781 + // Check for existence of DEC column
  782 + if (!table->contains("DEC")) {
  783 + missing.push_back("DEC");
  784 + }
  785 +
  786 + // Set error message for missing columns
  787 + if (!missing.empty()) {
  788 + message = "The following columns are missing in the "
  789 + "\""+evtname+"\" extension of input file \""+
  790 + filename + "\": ";
  791 + for (int i = 0; i < missing.size(); ++i) {
  792 + message += "\"" + missing[i] + "\"";
  793 + if (i < missing.size()-1) {
  794 + message += ", ";
  795 + }
  796 + }
  797 + }
  798 +
  799 + } // endelse: checked column names
  800 +
  801 + // Return
  802 + return message;
  803 +}
  804 +
  805 +
  806 +/***********************************************************************//**
  807 + * @brief Set output file name.
  808 + *
  809 + * @param[in] filename Input file name.
  810 + *
  811 + * Converts an input file name into an output filename by prepending the
  812 + * prefix stored in the member m_prefix to the input file name. Any path as
  813 + * well as extension will be stripped from the input file name. Also a
  814 + * trailing ".gz" will be stripped as one cannot write into gzipped files.
  815 + ***************************************************************************/
  816 +std::string ctprob::set_outfile_name(const std::string& filename) const
  817 +{
  818 + // Create filename
  819 + GFilename fname(filename);
  820 +
  821 + // Split input filename without any extensions into path elements
  822 + std::vector<std::string> elements = gammalib::split(fname.url(), "/");
  823 +
  824 + // The last path element is the filename
  825 + std::string outname = m_prefix + elements[elements.size()-1];
  826 +
  827 + // Strip any ".gz"
  828 + outname = gammalib::strip_chars(outname, ".gz");
  829 +
  830 + // Return output filename
  831 + return outname;
  832 +}
  833 +
  834 +
  835 +/***********************************************************************//**
  836 + * @brief Get Good Time Intervals extension name
  837 + *
  838 + * @param[in] filename Input file name.
  839 + * @param[in] evtname Events extension name.
  840 + *
  841 + * Extracts the Good Time Intervals extension name from the event file. We
  842 + * do this by loading the events and accessing the Good Time Intervals
  843 + * extension name using the GCTAEventList::gtiname() method. If the file name
  844 + * is empty, the method returns `GTI`.
  845 + ***************************************************************************/
  846 +std::string ctprob::get_gtiname(const std::string& filename,
  847 + const std::string& evtname) const
  848 +{
  849 + // Initialise GTI name
  850 + std::string gtiname = gammalib::extname_gti;
  851 +
  852 + // Continue only if the filename is not empty
  853 + if (!filename.empty()) {
  854 +
  855 + // Load events
  856 + GCTAEventList events(filename+"["+evtname+"]");
  857 +
  858 + // Get GTI name
  859 + gtiname = events.gtiname();
  860 +
  861 + }
  862 +
  863 + // Return GTI name
  864 + return (gtiname);
  865 +}
  866 +
  867 +
  868 +/***********************************************************************//**
  869 + * @brief Save event list in FITS format.
  870 + *
  871 + * Save the event list as a FITS file. The file name of the FITS file is
  872 + * specified by the "outobs" parameter.
  873 + ***************************************************************************/
  874 +void ctprob::save_fits(void)
  875 +{
  876 + // Save only if there are observations
  877 + if (m_obs.size() > 0) {
  878 +
  879 + // Get output filename
  880 + m_outobs = (*this)["outobs"].filename();
  881 +
  882 + // Get CTA observation from observation container
  883 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
  884 +
  885 + // Save only if it's a CTA observation
  886 + if (obs != NULL) {
  887 +
  888 + // Save only if file name is non-empty
  889 + if (m_infiles[0].length() > 0) {
  890 +
  891 + // Create file name object
  892 + GFilename fname(m_outobs);
  893 +
  894 + // Extract filename and event extension name
  895 + std::string outfile = fname.url();
  896 +
  897 + // Append event extension name. We handle here the possibility
  898 + // to write the events into a different extension.
  899 + if (fname.has_extname()) {
  900 + outfile += "["+fname.extname()+"]";
  901 + }
  902 + else {
  903 + outfile += "["+m_evtname[0]+"]";
  904 + }
  905 +
  906 + // Log filename
  907 + log_value(NORMAL, "Event list file", outfile);
  908 + // Save event list
  909 + save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0],
  910 + outfile);
756 911  
757   - // Get cube bin
758   - GCTAEventBin* bin = m_cube[i];
  912 + } // endif: filename was non empty
759 913  
760   - // Determine counts cube energy bin
761   - int iebin = cube_ebounds.index(bin->energy());
  914 + } // endif: observation was CTA observation
762 915  
763   - // Skip bin if the corresponding counts cube energy bin is not fully
764   - // contained in the event list energy range. This avoids having
765   - // partially filled bins.
766   - if (!usage[iebin]) {
767   - num_outside_ebds++;
  916 + } // endif: there were observations
  917 +
  918 + // Return
  919 + return;
  920 +}
  921 +
  922 +
  923 +/***********************************************************************//**
  924 + * @brief Save event list(s) in XML format.
  925 + *
  926 + * Save the event list(s) into FITS files and write the file path information
  927 + * into a XML file. The filename of the XML file is specified by the outfile
  928 + * parameter, the filename(s) of the event lists are built by prepending a
  929 + * prefix to the input event list filenames. Any path present in the input
  930 + * filename will be stripped, i.e. the event list(s) will be written in the
  931 + * local working directory (unless a path is specified in the prefix).
  932 + ***************************************************************************/
  933 +void ctprob::save_xml(void)
  934 +{
  935 + // Get output filename and prefix
  936 + m_outobs = (*this)["outobs"].filename();
  937 + m_prefix = (*this)["prefix"].string();
  938 +
  939 + // Issue warning if output filename has no .xml suffix
  940 + log_string(TERSE, warn_xml_suffix(m_outobs));
  941 +
  942 + // Loop over all observation in the container
  943 + for (int i = 0; i < m_obs.size(); ++i) {
  944 +
  945 + // Get CTA observation
  946 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  947 +
  948 + // Skip observations that are no CTA observations
  949 + if (obs == NULL) {
768 950 continue;
769 951 }
770 952  
771   - // If RoI is valid then skip bin if it is outside the RoI of the
772   - // observation
773   - if (roi.is_valid() && !roi.contains(*bin)) {
774   - num_outside_roi++;
  953 + // Skip observations that have empty names
  954 + if (m_infiles[i].length() == 0) {
775 955 continue;
776 956 }
777 957  
778   - // Get actual bin value
779   - double value = bin->counts();
780   -
781   - // Compute model value for cube bin
782   - double model = m_obs.models().eval(*bin, *obs) * bin->size();
  958 + // Set event output file name
  959 + std::string outfile = set_outfile_name(m_infiles[i]);
783 960  
784   - // Add model to actual value
785   - value += model;
786   - sum += model;
  961 + // Append event extension name
  962 + outfile += "["+m_evtname[i]+"]";
787 963  
788   - // Store value
789   - bin->counts(value);
  964 + // Log filename
  965 + log_value(NORMAL, "Event list file", outfile);
790 966  
791   - } // endfor: looped over all cube bins
  967 + // Store output file name in observation
  968 + obs->eventfile(outfile);
792 969  
793   - // Append GTIs of observation to list of GTIs
794   - m_gti.extend(gti);
  970 + // Save event list
  971 + save_event_list(obs, m_infiles[i], m_evtname[i], m_gtiname[i],
  972 + outfile);
795 973  
796   - // Update GTIs
797   - m_cube.gti(m_gti);
  974 + } // endfor: looped over observations
798 975  
799   - // Log filling results
800   - log_value(NORMAL, "Model events in cube", sum);
801   - log_value(NORMAL, "Bins outside energy range", num_outside_ebds);
802   - log_value(NORMAL, "Bins outside RoI", num_outside_roi);
  976 + // Write observation definition XML file name into logger
  977 + log_value(NORMAL, "Obs. definition file", m_outobs);
803 978  
804   - // Write model cube into header
805   - log_header2(EXPLICIT, "Model cube");
806   - log_string(EXPLICIT, m_cube.print(m_chatter));
  979 + // Save observations in XML file
  980 + m_obs.save(m_outobs);
807 981  
808 982 // Return
809 983 return;
810 984 }
811 985  
  986 +
  987 +/***********************************************************************//**
  988 + * @brief Save event list into FITS file
  989 + *
  990 + * @param[in] obs Pointer to CTA observation.
  991 + * @param[in] infile Input file name.
  992 + * @param[in] evtname Event extension name.
  993 + * @param[in] gtiname GTI extension name.
  994 + * @param[in] outfile Output file name.
  995 + *
  996 + * Saves an event list and the corresponding Good Time Intervals into a FITS
  997 + * file and copy all others extensions from the input file to the output
  998 + * file.
  999 + *
  1000 + * If an extension name is specified in the @p outfile argument, the events
  1001 + * and eventually also the Good Time Intervals will be extracted from the
  1002 + * argument and used for writing the events. The format is
  1003 + *
  1004 + * <filename>[<event extension name;GTI extension name>]
  1005 + *
  1006 + * where <filename> needs to be replaced by the name of the FITS file,
  1007 + * and <event extension name;GTI extension name> by the name of the events
  1008 + * and Good Time Intervals extensions. For example
  1009 + *
  1010 + * myfits.fits[EVENTS1;GTI1]
  1011 + *
  1012 + * will write the selected events into the "EVENTS1" extension and the
  1013 + * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS
  1014 + * file. If the Good Time Intervals extension name is skipped, e.g.
  1015 + *
  1016 + * myfits.fits[EVENTS1]
  1017 + *
  1018 + * the original extension name for the Good Time Intervals will be kept.
  1019 + * Analogously, only the Good Time Intervals extension name can be changed
  1020 + * by specifying
  1021 + *
  1022 + * myfits.fits[;GTI1]
  1023 + *
  1024 + * In none of the cases will the original events and Good Time Intervals be
  1025 + * copied over to the output file.
  1026 + ***************************************************************************/
  1027 +void ctprob::save_event_list(const GCTAObservation* obs,
  1028 + const std::string& infile,
  1029 + const std::string& evtname,
  1030 + const std::string& gtiname,
  1031 + const std::string& outfile) const
  1032 +{
  1033 + // Save only if we have an event list
  1034 + if (obs->eventtype() == "EventList") {
  1035 +
  1036 + // Set output FITS file event extension names
  1037 + GFilename outname(outfile);
  1038 + std::string outevt = evtname;
  1039 + std::string outgti = gtiname;
  1040 + if (outname.has_extname()) {
  1041 + std::vector<std::string> extnames =
  1042 + gammalib::split(outname.extname(), ";");
  1043 + if (extnames.size() > 0) {
  1044 + std::string extname = gammalib::strip_whitespace(extnames[0]);
  1045 + if (!extname.empty()) {
  1046 + outevt = extname;
  1047 + }
  1048 + }
  1049 + if (extnames.size() > 1) {
  1050 + std::string extname = gammalib::strip_whitespace(extnames[1]);
  1051 + if (!extname.empty()) {
  1052 + outgti = extname;
  1053 + }
  1054 + }
  1055 + }
  1056 +
  1057 + // Create output FITS file
  1058 + GFits outfits;
  1059 +
  1060 + // Write observation into FITS file
  1061 + obs->write(outfits, outevt, outgti);
  1062 +
  1063 + // Copy all extensions other than evtname and gtiname extensions
  1064 + // from the input to the output event list. The evtname and
  1065 + // gtiname extensions are written by the save method, all others
  1066 + // that may eventually be present have to be copied over
  1067 + // explicitly.
  1068 + GFits infits(infile);
  1069 + for (int extno = 1; extno < infits.size(); ++extno) {
  1070 + GFitsHDU* hdu = infits.at(extno);
  1071 + if (hdu->extname() != evtname &&
  1072 + hdu->extname() != gtiname &&
  1073 + hdu->extname() != outevt &&
  1074 + hdu->extname() != outgti) {
  1075 + outfits.append(*hdu);
  1076 + }
  1077 + }
  1078 +
  1079 + // Close input file
  1080 + infits.close();
  1081 +
  1082 + // Save file to disk and close it (we need both operations)
  1083 + outfits.saveto(outname.url(), clobber());
  1084 + outfits.close();
  1085 +
  1086 + } // endif: observation was unbinned
  1087 +
  1088 + // Return
  1089 + return;
  1090 +}
... ...
src/ctprob/ctprob.hpp
  View file @6985d99
... ... @@ -24,117 +24,91 @@
24 24 * @author Leonardo Di Venere
25 25 */
26 26  
27   -
28   -#ifndef CTMODEL_HPP
29   -#define CTMODEL_HPP
  27 +#ifndef CTPROB_HPP
  28 +#define CTPROB_HPP
30 29  
31 30 /* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
32 33 #include "GammaLib.hpp"
33 34 #include "GCTALib.hpp"
34 35 #include "ctobservation.hpp"
35 36  
36 37 /* __Definitions _________________________________________________________ */
37   -#define CTMODEL_NAME "ctmodel"
38   -#define CTMODEL_VERSION "1.1.0"
  38 +#define CTPROB_NAME "ctprob"
  39 +#define CTPROB_VERSION "0.0.1"
39 40  
40 41  
41 42 /***********************************************************************//**
42   - * @class ctmodel
43   - *
44   - * @brief Model cube generation tool
45   - *
46   - * This class creates counts model map(s). The definition of the counts model
47   - * can be taken from a predefined observation container, from the counts maps
48   - * found in an observation definition XML file, from an individual counts
49   - * map, or from the parameter definition.
  43 + * @class ctprob
50 44 *
51   - * Results are stored in an observation container that can be written to disk
52   - * in form of FITS files (model maps) and an updated observation definition
53   - * XML file.
  45 + * @brief Data selection tool
54 46 ***************************************************************************/
55   -class ctmodel : public ctobservation {
  47 +class ctprob : public ctobservation {
56 48  
57 49 public:
58 50 // Constructors and destructors
59   - ctmodel(void);
60   - explicit ctmodel(const GObservations& obs);
61   - ctmodel(int argc, char *argv[]);
62   - ctmodel(const ctmodel& app);
63   - virtual ~ctmodel(void);
  51 + ctprob(void);
  52 + explicit ctprob(const GObservations& obs);
  53 + ctprob(int argc, char *argv[]);
  54 + ctprob(const ctprob& app);
  55 + virtual ~ctprob(void);
64 56  
65 57 // Operators
66   - ctmodel& operator=(const ctmodel& app);
  58 + ctprob& operator=(const ctprob& app);
67 59  
68 60 // Methods
69   - void clear(void);
70   - void run(void);
71   - void save(void);
72   - void publish(const std::string& name = "");
73   - const GCTAEventCube& cube(void) const;
74   - void cube(const GCTAEventCube& cube);
75   - void models(const GModels& models);
  61 + void clear(void);
  62 + void run(void);
  63 + void save(void);
  64 + void publish(const std::string& name = "");
76 65  
77 66 protected:
78 67 // Protected methods
79   - void init_members(void);
80   - void copy_members(const ctmodel& app);
81   - void free_members(void);
82   - void get_parameters(void);
83   - void get_obs(void);
84   - void fill_cube(const GCTAObservation* obs);
85   - bool has_cube(void) const;
  68 + void init_members(void);
  69 + void copy_members(const ctprob& app);
  70 + void free_members(void);
  71 + void get_parameters(void);
  72 + 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;
  76 + 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 + void save_fits(void);
  80 + 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 86  
87 87 // User parameters
88   - GFilename m_outcube; //!< Output model cube
89   - bool m_apply_edisp; //!< Apply energy dispersion?
90   - bool m_publish; //!< Publish model cube?
91   - GChatter m_chatter; //!< Chattiness
  88 + std::string m_outobs; //!< Output event list or XML file
  89 + std::string m_prefix; //!< Prefix for multiple event lists
  90 + // bool m_usepnt; //!< Use pointing instead of RA/DEC parameters
  91 + //GCTARoi m_roi; //!< RoI selection
  92 + //double m_tmin; //!< Start time
  93 + //double m_tmax; //!< Stop time
  94 + //double m_emin; //!< Lower energy
  95 + //double m_emax; //!< Upper energy
  96 + //std::string m_expr; //!< Selection expression
  97 + //std::string m_usethres; //!< Energy threshold type
  98 + bool m_apply_edisp; //!< Apply energy dispersion?
  99 + bool m_publish; //!< Publish model cube?
  100 + GChatter m_chatter; //!< Chattiness
92 101  
93 102 // Protected members
94   - GCTAEventCube m_cube; //!< Model cube
95   - GGti m_gti; //!< Model cube GTIs
96   - bool m_has_cube; //!< Signal if cube has been set or loaded
97   - bool m_append_cube; //!< Signal that cube should be appended
98   - bool m_binned; //!< Signals that we are in binned mode
  103 + std::vector<std::string> m_infiles; //!< Input event filenames
  104 + std::vector<std::string> m_evtname; //!< Event extension names
  105 + std::vector<std::string> m_gtiname; //!< GTI extension names
  106 + //GTime m_timemin; //!< Earliest time
  107 + //GTime m_timemax; //!< Latest time
  108 + //bool m_select_energy; //!< Perform energy selection
  109 + //bool m_select_roi; //!< Perform RoI selection
  110 + //bool m_select_time; //!< Perform time selection
99 111 };
100 112  
101 113  
102   -/***********************************************************************//**
103   - * @brief Return model cube
104   - *
105   - * @return Reference to model cube.
106   - ***************************************************************************/
107   -inline
108   -const GCTAEventCube& ctmodel::cube(void) const
109   -{
110   - return m_cube;
111   -}
112   -
113   -
114   -/***********************************************************************//**
115   - * @brief Signal if cube has been set or loaded
116   - *
117   - * @return True if cube has been set or loaded.
118   - ***************************************************************************/
119   -inline
120   -bool ctmodel::has_cube(void) const
121   -{
122   - return m_has_cube;
123   -}
124   -
125   -
126   -/***********************************************************************//**
127   - * @brief Set models
128   - *
129   - * @param[in] models Model container.
130   - *
131   - * Set model container that should be used for model generation.
132   - ***************************************************************************/
133   -inline
134   -void ctmodel::models(const GModels& models)
135   -{
136   - m_obs.models(models);
137   - return;
138   -}
139   -
140   -#endif /* CTMODEL_HPP */
  114 +#endif /* CTPROB_HPP */
... ...
src/ctprob/ctprob.par
  View file @6985d99
... ... @@ -16,53 +16,36 @@
16 16 #
17 17 # General parameters
18 18 #===================
19   -inobs, f, a, "NONE",,, "Input event list, counts cube or observation definition XML file"
20   -incube, f, a, "NONE",,, "Input counts cube file to extract model cube definition"
  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"
21 22 inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file"
22   -expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
23   -psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
24   -edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
25   -bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
  23 +#expcube, f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
  24 +#psfcube, f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
  25 +#edispcube,f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
  26 +#bkgcube, f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
26 27 caldb, s, a, "prod2",,, "Calibration database"
27 28 irf, s, a, "South_0.5h",,, "Instrument response function"
28 29 edisp, b, h, no,,, "Apply energy dispersion?"
29   -outcube, f, a, "modcube.fits",,, "Output model cube file"
30 30  
31 31 #
32   -# Observation parameters
33   -#=======================
34   -ra, r, a, 83.63,0,360, "RA of pointing (degrees)"
35   -dec, r, a, 22.01,-90,90, "Dec of pointing (degrees)"
36   -rad, r, a, 5.0,0,180, "Radius of FOV (degrees)"
37   -tmin, r, a, 0.0,,, "Start time (MET in s)"
38   -tmax, r, a, 1800.0,,, "End time (MET in s)"
39   -deadc, r, h, 0.98,0,1, "Average deadtime correction factor"
40   -
41   -#
42   -# Energy binning parameters
43   -#==========================
44   -ebinalg, s, h, "LOG",FILE|LIN|LOG,, "Algorithm for defining energy bins"
45   -emin, r, a, 0.1,,, "Start value for first energy bin in TeV"
46   -emax, r, a, 100.0,,, "Stop value for last energy bin in TeV"
47   -enumbins, i, a, 20,,, "Number of energy bins"
48   -ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition"
49   -
50   -#
51   -# Spatial binning parameters
52   -#===========================
53   -usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
54   -nxpix, i, a, 200,,, "Size of the X axis in pixels"
55   -nypix, i, a, 200,,, "Size of the Y axis in pixels"
56   -binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
57   -coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
58   -xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
59   -yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
60   -proj, s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
  32 +# Selection parameters
  33 +#=====================
  34 +#usepnt, b, h, no,,, "Use pointing instead of RA/DEC parameters?"
  35 +#ra, r, a, 83.63,0,360, "RA for ROI centre (degrees)"
  36 +#dec, r, a, 22.01,-90,90, "Dec for ROI centre (degrees)"
  37 +#rad, r, a, 3.0,0,180, "Radius of ROI (degrees)"
  38 +#tmin, r, a, 0.0,,, "Start time (CTA MET in seconds)"
  39 +#tmax, r, a, 0.0,,, "End time (CTA MET in seconds)"
  40 +#emin, r, a, 0.1,,, "Lower energy limit (TeV)"
  41 +#emax, r, a, 100.0,,, "Upper energy limit (TeV)"
  42 +#expr, s, h, "",,, "Additional selection expression"
  43 +#usethres,s, h, "NONE",NONE|DEFAULT|USER,, "Energy threshold type"
61 44  
62 45 #
63 46 # Standard parameters
64 47 #====================
65   -publish, b, h, no,,, "Publish model cube on VO Hub?"
  48 +publish, b, h, no,,, "Publish event list on VO Hub?"
66 49 chatter, i, h, 2,0,4, "Chattiness of output"
67 50 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
68 51 debug, b, h, no,,, "Debugging mode activated"
... ...
src/ctprob/main.cpp
  View file @6985d99
... ... @@ -29,7 +29,7 @@
29 29 #include <config.h>
30 30 #endif
31 31 #include "support.hpp"
32   -#include "ctmodel.hpp"
  32 +#include "ctprob.hpp"
33 33  
34 34  
35 35 /***********************************************************************//**
... ... @@ -38,8 +38,8 @@
38 38 * @param[in] argc Number of command line arguments.
39 39 * @param[in] argv Command line arguments.
40 40 *
41   - * This is the main entry point of the ctmodel application. It allocates a
42   - * ctmodel object and executes the application. Any exceptions that occur
  41 + * This is the main entry point of the ctprob application. It allocates a
  42 + * ctprob object and executes the application. Any exceptions that occur
43 43 * will be catched and corresponding error messages written in the
44 44 * application logger and into the standard output.
45 45 ***************************************************************************/
... ... @@ -49,13 +49,13 @@ int main (int argc, char *argv[])
49 49 int rc = 1;
50 50  
51 51 // Initialise pointer on application
52   - ctmodel* application = NULL;
  52 + ctprob* application = NULL;
53 53  
54 54 // Execute application
55 55 try {
56 56  
57 57 // Create instance of application
58   - application = new ctmodel(argc, argv);
  58 + application = new ctprob(argc, argv);
59 59  
60 60 // Execute ctool
61 61 rc = execute_ctool(application);
... ... @@ -65,7 +65,7 @@ int main (int argc, char *argv[])
65 65 catch (std::exception &e) {
66 66  
67 67 // Report exception
68   - report_ctool_failure("ctmodel", e.what());
  68 + report_ctool_failure("ctprob", e.what());
69 69  
70 70 }
71 71  
... ...