diff --git a/ChangeLog b/ChangeLog index 4805d13..706dd1d 100755 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2018-11-23 +2018-12-11 * Version 1.6.0 released ======================== @@ -7,6 +7,7 @@ CTA background rates are now per ontime (#2758) Apply event bin weighting in ctmodel for stacked computation (#2721) Enable ctbin to produce counts cubes for a joint binned analysis (#2714) + Add ctfindvar tool (#2713) Add Jupyter notebooks for H.E.S.S. DR1 (#2712) csphagen now creates appropriate model for likelihood fitting (#2711) csphagen now supports any kind of background model (#2710) diff --git a/NEWS b/NEWS index 60f4333..55437be 100755 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.6.0 -23 November 2018 +11 December 2018 Introduction @@ -126,6 +126,12 @@ ctexpcube - CTA exposure cube generation None +ctfindvar - Variability search +------------------------------ +New: the ctfindvar tool searches for source variability using an On/Off method +in the time domain (#2713). + + ctlike - Maximum likelihood fitting ----------------------------------- Add hidden "like_accuracy" parameter to allow control of the accuracy of diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 8dc4ea8..a75b110 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -14,6 +14,7 @@ In particular, this release provides: * Support for Xspec spectral analysis * Improved computation speed for :ref:`ctskymap` and :ref:`csphagen` * Parallelisation of various cscripts +* Support for source variability searches * Support for Python pickeling * A Code Generator to facilitate adding ctools and cscripts @@ -102,6 +103,8 @@ Improvements New features ------------ +* [`2713 `_] - + Add :ref:`ctfindvar` tool * [`2709 `_] - Add code generator * [`2694 `_] - diff --git a/doc/source/users/reference_manual/ctfindvar.rst b/doc/source/users/reference_manual/ctfindvar.rst index 1cc5204..d688200 100644 --- a/doc/source/users/reference_manual/ctfindvar.rst +++ b/doc/source/users/reference_manual/ctfindvar.rst @@ -1,50 +1,50 @@ .. _ctfindvar: ctfindvar -============================== +========= -ToDo: Describe in a one liner what the tool is doing. +Search for source variability. Synopsis -------- -This tool is used to search for source variability. -The variability is computed by binning the arrival time of the photons -and applying the on-off method to each bin, integrating the background -over all the other bins. -A Limit is set to include a bin in the background sample, i.e., its -significance with respect to all the others bins has to be lower than 4.5 -The significance of each bin is computed using Li&Ma eq. 17. -A limit on the background sample can also be included, taking only into -account bins that have a number of counts which is higher than a threshold -defined by the user via the "minoff" parameter. - -Variability is calculated for each pixel of the sky given in the input files. -The results of the variability is further stored in the srcsig.fits file, -where the GTIs and the significance of each bin are stored. -Only the sources with a significance > 4.5 sigmas are stored. - -:ref:`ctfindvar` takes in input some observation files. -It calculates a cube with several maps, each one with a count number integrated -over a time period defined by the tinterval parameter, for each pixel. -The tinterval parameter is the time scale over which the variability -will be searched. - -Inputs files: -- observation file xml, or event list .fits - -Output files: -- countscube.fits - provides all the count-skymaps contained in the cube. -- peaksigmap.fits - a skymap of the highest significance for each pixel -- srcsig.fits - a table containing the GTIs together with the significance - for all the sources with a significant variability detected. - - -The significance of a source can be displayed using the -``show_variability_evolution.py`` script in the example forlder + +This tool searches for source variability by applying the On/Off method in the +time domain. The variability of sky map pixels is computed by binning the events +as function of the photon arrival time and applying the On/Off method to each +time bin, integrating the background over all the other time bins. Only bins +with a variabilty significance below a given threshold, defined by the hidden +``threshold`` parameter, are included in the background estimate. Furthermore, +only time bins that exceed a minimum number of events, specified by the hidden +``minoff`` parameter, are considered for the background estimate. The significance +of each time bin is computed using Li & Ma formula. + +:ref:`ctfindvar` will search a grid of sky map pixels for source variability. Using +the hidden ``smooth_kernel`` and ``smooth_rad`` parameters, the sky map may be +smoothed before the variability search. The tool will compute sky maps for a +number of time bins. The duration of each time bin is specified by the ``tinterval`` +parameter, the start and stop time are defined by the ``tmin`` and ``tmax`` +parameters. If the start and stop time are ``NONE``, the start and stop times will +be extracted from the input event list(s). Note that ``tinterval`` specifies the +time scale over which the variability will be searched. + +The tool will determine the sky map pixel with the maximum variability and store +the variability histogram and pixel position as ``MAXSIGPIXEL`` in the output file, +specified using the ``outmap`` parameter. In addition, the tool will determine the +variability for a number of source positions. If an input model is specified +using the ``inmodel`` parameter, the source positions will be defined by the centre +positions of all gamma-ray source components in the input model. If no input +model is specified, the tool will use the ``xsrc`` and ``ysrc`` parameters to allow +the specification of a source position. The results will be stored in the output +file using either the source names found in the model definition file, or the +name ``SOURCE`` if the source position was specified by the ``xsrc`` and ``ysrc`` +parameters. + +All sky map pixels with a variability above the value specified by the ``threshold`` +will also added be as a point source model to the ``outmodel`` file. Each source will +have a power law as the spectral component with a total flux that corresponds to +1% of the Crab nebula flux. + General parameters ------------------ @@ -56,7 +56,17 @@ General parameters Input model definition file for extracting source positions. ``(outcube = NONE) [file]`` - Filename for saving counts cube, NONE will result in no cube being saved. + Filename for saving counts cube, ``NONE`` will result in no cube being saved. + +``outmap [file]`` + Filename for output sky map containing the maximum significance for each + spatial pixel. The file also contains the variability significance as function + of time for each source of interest, and the pixel with the maximum + significance. Furthermore, it contains the Right Ascension and Declination + of each source of interest and the pixel with the maximum significance. + +``outmodel [file]`` + Output model definition file containing sources for significant pixels. ``caldb [string]`` Calibration database (only required for IRF background subtraction if no @@ -66,40 +76,45 @@ General parameters Instrument response function (only required for IRF background subtraction if no response information is provided by ``inobs``). -``prefix [string]`` - Output file prefix. The method will save two files: - - Skymap containing the maximum significance for each spatial pixel - - File containing the significance as a function of time interval for each - source -``(minoff = 0) [real]`` - Minimum number of off counts to compute the significance +Variability search parameters +----------------------------- + +``coordsys [string]`` + Coordinate system (CEL - celestial, GAL - galactic). + +``xsrc [real]`` + Right Ascension / Galactic longitude of source of interest (J2000, in degrees) + +``ysrc [real]`` + Declination / Galactic latitude of source of interest (J2000, in degrees) ``emin [real]`` - Minimum energy (TeV) for extracting events. A value of 0 results in no cut - being applied. + Minimum energy (TeV) for extracting events. ``emax [real]`` - Maximum energy (TeV) for extracting events. A value of 0 results in no cut - being applied. + Maximum energy (TeV) for extracting events. -``(sig_thr = 4.5) [real]`` +``(threshold = 4.5) [real]`` Significance threshold for variability detection. -``(smoothkrnl = NONE) [string]`` +``(minoff = 0) [real]`` + Minimum number of off counts to compute the significance. + +``(smooth_kernel = NONE) [string]`` Kernel to be used in smoothing the counts maps before computing significances. -``(smoothpar = 0.05) [real]`` - Parameter for the smoothing kernel. This is either the sigma of the GAUSSIAN - or the radius of the DISK. +``(smooth_rad = 0.05) [real]`` + Radius of smoothing kernel in degrees. This is either the sigma of the + ``GAUSSIAN`` or the radius of the ``DISK``. Time binning parameters ----------------------- ``tinterval [real]`` - interval in seconds for each time step. + Interval in seconds for each time step. ``tmin [time]`` Start time for interval determination (UTC string, JD, MJD or MET in seconds). @@ -120,7 +135,7 @@ Spatial binning parameters -------------------------- ``(usepnt = no) [boolean]`` - Use CTA pointing direction for cube centre instead of xref/yref parameters? + Use CTA pointing direction for sky map centre instead of xref/yref parameters? ``nxpix [integer]`` Number of cube bins in Right Ascension or Galactic longitude. @@ -129,25 +144,16 @@ Spatial binning parameters Number of cube bins in Declination or Galactic latitude. ``binsz [real]`` - Cube bin size (in degrees/pixel). - -``coordsys [string]`` - Coordinate system (CEL - celestial, GAL - galactic). + Sky map pixel size (in degrees/pixel). ``proj [string]`` Projection method. ``xref [real]`` - Right Ascension / Galactic longitude of cube centre (J2000, in degrees). + Right Ascension / Galactic longitude of sky map centre (J2000, in degrees). ``yref [real]`` - Declination / Galactic latitude of cube centre (J2000, in degrees). - -``xsrc [real]`` - Right Ascension / Galactic longitude of source of interest (J2000, in degrees) - -``ysrc [real]`` - Declination / Galactic latitude of source of interest (J2000, in degrees) + Declination / Galactic latitude of sky map centre (J2000, in degrees). Standard parameters diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst index 1eb4d33..9159f89 100755 --- a/doc/source/users/reference_manual/index.rst +++ b/doc/source/users/reference_manual/index.rst @@ -17,7 +17,6 @@ ctools .. toctree:: :maxdepth: 1 - ctfindvar --- ToDo: Describe what tool is doing ctbin --- Generates counts cube ctbkgcube --- Generates background cube ctbutterfly --- Compute butterfly @@ -25,6 +24,7 @@ ctools ctedispcube --- Generates energy dispersion cube cterror --- Calculates likelihood profile errors ctexpcube --- Generates exposure cube + ctfindvar --- Search for source variability ctlike --- Performs maximum likelihood fitting ctmapcube --- Generates a map cube ctmodel --- Computes model counts cube diff --git a/examples/show_variability_evolution.py b/examples/show_variability_evolution.py index b4794eb..9f9adf1 100755 --- a/examples/show_variability_evolution.py +++ b/examples/show_variability_evolution.py @@ -1,8 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# Shows butterfly diagram created with ctbutterfly +# Shows source variability evolution created with ctfindvar # -# Copyright (C) 2014-2017 Simon Bonnefoy +# Copyright (C) 2018 Simon Bonnefoy # # 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 @@ -43,6 +43,8 @@ def plot_variability_evolution(filename, plotfile, srcname): Name of spectrum FITS file plotfile : str Plot file name + srcname : str + Source name """ # Read variability file fits = gammalib.GFits(filename) @@ -51,41 +53,45 @@ def plot_variability_evolution(filename, plotfile, srcname): c_gti_stop = table['TSTOP'] c_sig = table[srcname] - #convert the c_sig object to list - table_c_sig=[] - table_c_gti=[] - + # Convert the c_sig object to list + table_c_sig = [] + table_c_gti = [] for k in c_sig: table_c_sig.append(k) - for j,k in zip(c_gti_start, c_gti_stop): + for j, k in zip(c_gti_start, c_gti_stop): table_c_gti.append(0.5*(j+k)) - #Initializing the plots - axarr=[] - f = plt.figure(figsize=(12,7)) - ax1 = f.add_subplot(121) - ax2 = f.add_subplot(122) + # Initializing the plots + axarr = [] + f = plt.figure(figsize=(12,7)) + ax1 = f.add_subplot(121) + ax2 = f.add_subplot(122) axarr.append(ax1) axarr.append(ax2) + # Create first plot axarr[0].plot(table_c_gti, table_c_sig) axarr[0].set_title('Significance evolution') axarr[0].set_xlabel('MJD') axarr[0].set_ylabel('Significance') - axarr[1].hist(table_c_sig, 20, range=[min(table_c_sig), max(table_c_sig)], histtype='stepfilled', facecolor='g', alpha=0.75) + # Create second plot + axarr[1].hist(table_c_sig, 20, range=[min(table_c_sig), max(table_c_sig)], + histtype='stepfilled', facecolor='g', alpha=0.75) axarr[1].set_title('Significance distribution') axarr[1].set_xlabel('Significance') axarr[1].set_ylabel('Counts') - #IT SEEMS THAT THE LOG SCALE IS NOT WORKING WELL !!!!!! - #axarr[1].gca().set_yscale('log') - + # Optionally save result in file if plotfile != "": f.savefig(plotfile+".png") plt.show() + # Return + return + + # ========================== # # Show variability evolution # # ========================== # diff --git a/src/ctfindvar/ctfindvar.cpp b/src/ctfindvar/ctfindvar.cpp index 1da10a0..088bb62 100644 --- a/src/ctfindvar/ctfindvar.cpp +++ b/src/ctfindvar/ctfindvar.cpp @@ -39,6 +39,7 @@ #endif /* __ Method name definitions ____________________________________________ */ +#define G_INX2GTI "ctfindvar::inx2gti(int&)" #define G_FILL_CUBE "ctfindvar::fill_cube(GCTAObservation*)" #define G_FILL_ALPHA_VECTOR "ctfindvar::fill_alpha_vector(const int&, "\ "vector&)" @@ -235,146 +236,16 @@ void ctfindvar::run(void) // Write input observation container into logger log_observations(NORMAL, m_obs, "Input observation"); - // Loop over all unbinned CTA observations in the container - #pragma omp parallel for - for (int i = 0; i < m_obs.size(); ++i) { - - // Get pointer to observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); - - // Fill the cube - fill_cube(obs); - - // Dispose events to free memory - obs->dispose_events(); - - } // endfor: looped over observations - - // Smooth the maps if requested - if ((*this)["smoothkrnl"].is_valid() && - (*this)["smoothkrnl"].string() != "NONE") { - m_counts.smooth((*this)["smoothkrnl"].string(), - (*this)["smoothpar"].real()); - } - - // Get number of maps - const int nbins = m_counts.nmaps(); - - // creating GSkyDir to get the position of the source and each pixels - - // Extract which pixel in the map each source is located in - std::vector srcInxPix; - if (m_inmodel.size() > 0) { - for (int mod = 0; mod < m_inmodel.size(); ++mod) { - - // Extract the model - GModelSky* model = dynamic_cast(m_inmodel[mod]); - if (model != NULL) { - // Get the source position - GModelSpatial* model_spatial = model->spatial(); - GSkyRegion* model_region = model_spatial->region(); - GSkyDir srcSkyDir = - dynamic_cast(model_region)->centre(); - - // Store the map index of this source - int src_map_index = m_counts.dir2inx(srcSkyDir); - srcInxPix.push_back(src_map_index); - } else { - m_inmodel.remove(mod); - mod--; - } - } - } - // ... otherwise only one source position is queried - else { - GSkyDir srcSkyDir; - if ((*this)["coordsys"].string() == "CEL") { - srcSkyDir.radec_deg((*this)["xsrc"].real(),(*this)["ysrc"].real()); - } else { - srcSkyDir.lb_deg((*this)["xsrc"].real(),(*this)["ysrc"].real()); - } - - // Store the index of this source - int src_map_index = m_counts.dir2inx(srcSkyDir); - srcInxPix.push_back(src_map_index); - } - - //preparing the histogram with significance evolution - //for the source center and the highest sig-pix - m_pixsigsrc = GNdarray(srcInxPix.size(), nbins); - m_pixsigmax = GNdarray(nbins); - - //Prepare the final skymap with the max significance of each pixel - double max_sig = 0; - - log_header1(NORMAL, "Looping Over Pixels"); - - //looping over all the pixels in the cube - #pragma omp parallel for - for (int pix_number = 0; pix_number < m_counts.npix(); pix_number++) { - - //Preparing array to store significance evolution of - //each pixel - GNdarray pixSig(nbins); - double total_counts=0; - - //Checking the integrated count of each pixel - for (int k = 0; k < m_counts.nmaps(); k++) { - total_counts += m_counts(pix_number, k); - } - - //Getting the variability significance for the pixel - get_variability_sig(pix_number,nbins, pixSig); - - //Getting the significance evolution for the source - for (int src = 0; src < srcInxPix.size(); src++) { + // Create GTIs + init_gtis(); - // Store the distribution if the source is located at this position - if (srcInxPix[src] == pix_number) { + // Create counts cube + create_cube(); - #pragma omp critical(ctfind_run) - for (int i = 0; i < nbins; i++) { - m_pixsigsrc(src,i) = pixSig(i); - } - } - #ifdef G_DEBUG - std::cout << "checking pixel number of the source of interest" << pix_number << std::endl; - std::cout << "number of counts in pixel of the source of interest" << total_counts << std::endl; - pixSigSrc=pixSig; - std::cout << "checking pixel number " << pix_number << " with total number of counts of: " << total_counts << std::endl; - std::cin.ignore(); - std::cout << "checking pixel number of the source of interest: " << pix_number << " - with total number of counts of: " << total_counts << std::endl; - #endif - } - } - - // Getting the evolution for the pix with highest significance - #pragma omp critical(ctfind_run) - { - if (max(pixSig) > max_sig) { - // Update the pixel position with the maximum sigma - max_sig = max(pixSig); - m_max_sig_dir = m_counts.inx2dir(pix_number); - m_pixsigmax = pixSig; - } - } - //storing sig in skymap - m_peaksigmap(pix_number) = max(pixSig); - - //If the significance is greater than sig threshold, the position - //is stored in a model - if (max(pixSig) > m_sig_threshold) { - GSkyDir dir_pix = m_counts.inx2dir(pix_number); - fill_model_sig_pos(dir_pix); - } - - } - - log_header1(NORMAL, "Analysis finished"); - log_value(NORMAL, "Maximum sigma", max_sig); - log_value(NORMAL, "Max pixel RA", m_max_sig_dir.ra_deg()); - log_value(NORMAL, "Max pixel DEC", m_max_sig_dir.dec_deg()); + // Analyse cube + analyse_cube(); + // Return return; } @@ -388,37 +259,52 @@ void ctfindvar::run(void) ***************************************************************************/ void ctfindvar::save(void) { + // Write header log_header1(TERSE, "Saving results"); - // Filenames - std::string prefix((*this)["prefix"].string()); - // Write counts cube if ((*this)["outcube"].is_valid()) { + + // Get counts cube output filename GFilename outcube = (*this)["outcube"].filename(); - log_value(TERSE, "Saving counts cube", outcube); + + // Save counts cube m_counts.save(outcube, (*this)["clobber"].boolean()); + + // Log saving + log_value(TERSE, "Saving counts cube", outcube); } // Create a FITS file for storing the output - m_outfile.clear(); + GFits fits; // Write the most significant values for each pixel - log_string(TERSE, "Appending maximum significances map to output file"); - m_peaksigmap.write(m_outfile, "PEAKSIGMAP"); + m_peaksigmap.write(fits, "PEAKSIGMAP"); + + // Write source histograms + write_source_histograms(fits); - // Write individual source distributions - write_srchist(); + // Get output map filename + GFilename outmap = (*this)["outmap"].filename(); - // Write the output file - GFilename outfilename = prefix + "signifmap.fits"; - log_value(TERSE, "Saving results to ", outfilename); - m_outfile.saveto(outfilename, (*this)["clobber"].boolean()); + // Save output map + fits.saveto(outmap, (*this)["clobber"].boolean()); - //Save the model file with position with significance + // Log saving + log_value(TERSE, "Saving output map", outmap); + + // Save the model definition file if (m_model_above_thr.size() > 2) { - GFilename model_filename = prefix + "model_significant_directions.xml"; - m_model_above_thr.save(model_filename); + + // Get model output filename + GFilename outmodel = (*this)["outmodel"].filename(); + + // Save model + m_model_above_thr.save(outmodel); + + // Log saving + log_value(TERSE, "Saving output model", outmodel); + } // Return @@ -429,7 +315,7 @@ void ctfindvar::save(void) /***********************************************************************//** * @brief Get the map index associated with a given time * - * @param[in] time Time + * @param[in] time Time * @return Index of map in counts cube ***************************************************************************/ int ctfindvar::time2inx(const GTime& time) @@ -438,29 +324,35 @@ int ctfindvar::time2inx(const GTime& time) int map_index = -1; // Loop over all GTIs - for (int i = 0; i < m_gti.size(); i++) { - // Check if interval contains the time + for (int i = 0; i < m_gti.size(); ++i) { + + // If interval contains the time then break if (m_gti[i].contains(time)) { map_index = i; break; } + } + // Return map index return map_index; } /***********************************************************************//** - * @brief Return reference to significance cube object + * @brief Return Good Time Interval for index + * + * @param[in] index Index + * @return Good Time Interval ***************************************************************************/ GGti ctfindvar::inx2gti(const int& index) { + // Throw an exception if index is invalid if ((index < 0) || (index >= m_gti.size())) { - // Index is invalid, so throw an error - throw GException::invalid_value("ctfindvar::inx2gti(const int&)", - "'index' parameter out of range"); + throw GException::out_of_range(G_INX2GTI, "Time index", index, m_gti.size()); } + // Return Good Time Interval return m_gti[index]; } @@ -483,7 +375,6 @@ void ctfindvar::init_members(void) m_max_sig_dir.clear(); m_minoff = 0.0; m_sig_threshold = 0.0; - m_outfile.clear(); m_peaksigmap.clear(); m_pixsigsrc.clear(); m_tstart.clear(); @@ -511,7 +402,6 @@ void ctfindvar::copy_members(const ctfindvar& app) m_max_sig_dir = app.m_max_sig_dir; m_minoff = app.m_minoff; m_sig_threshold = app.m_sig_threshold; - m_outfile = app.m_outfile; m_peaksigmap = app.m_peaksigmap; m_pixsigsrc = app.m_pixsigsrc; m_tstart = app.m_tstart; @@ -537,39 +427,35 @@ void ctfindvar::free_members(void) /***********************************************************************//** * @brief Get application parameters - * - * @todo Implement method ***************************************************************************/ void ctfindvar::get_parameters(void) { - #ifdef _OPENMP - int nthreads = (*this)["nthreads"].integer(); - if (nthreads > 0) { - omp_set_num_threads(nthreads); - } - #endif - // Load the observations setup_observations(m_obs, true, true, false); - // Get the rest of the parameters - if (read_ahead()) { - (*this)["prefix"].string(); - - if ((*this)["inmodel"].is_valid()) { - m_inmodel = GModels((*this)["inmodel"].filename()); - } else { - // Get the desired source position from user - (*this)["xsrc"].real(); - (*this)["ysrc"].real(); - } + // Either load model or query source position parameters + if ((*this)["inmodel"].is_valid()) { + m_inmodel = GModels((*this)["inmodel"].filename()); + } + else { + (*this)["coordsys"].string(); + (*this)["xsrc"].real(); + (*this)["ysrc"].real(); } - // Create GTIs and counts cube - init_gtis(); - init_cube(); + // Query time limit parameters + (*this)["tinterval"].real(); + if ((*this)["tmin"].is_valid()) { + (*this)["tmin"].time(); + } + if ((*this)["tmax"].is_valid()) { + (*this)["tmax"].time(); + } + + // Create map to query spatial parameters + create_map(m_obs); - // Get the energy limits + // Read energy limits m_emin = GEnergy((*this)["emin"].real(), "TeV"); m_emax = GEnergy((*this)["emax"].real(), "TeV"); @@ -577,7 +463,28 @@ void ctfindvar::get_parameters(void) m_minoff = (*this)["minoff"].real(); // Get minimum significance to set a source as variable - m_sig_threshold = (*this)["sig_thr"].real(); + m_sig_threshold = (*this)["threshold"].real(); + + // Query smoothing parameters + if ((*this)["smooth_kernel"].is_valid()) { + (*this)["smooth_kernel"].string(); + (*this)["smooth_rad"].real(); + } + + // If needed later, query output filenames now + if (read_ahead()) { + (*this)["outcube"].query(); + (*this)["outmap"].query(); + (*this)["outmodel"].query(); + } + + // Set number of OpenMP threads + #ifdef _OPENMP + int nthreads = (*this)["nthreads"].integer(); + if (nthreads > 0) { + omp_set_num_threads(nthreads); + } + #endif // Write parameters into logger log_parameters(TERSE); @@ -588,117 +495,139 @@ void ctfindvar::get_parameters(void) /***********************************************************************//** - * @brief Initialize counts cube + * @brief Initialize Good Time Intervals * - * Sets up the required event cube and time interval storage objects - * associated with this object - ***************************************************************************/ -void ctfindvar::init_cube(void) -{ - // Create the basic skymap - m_counts = create_map(m_obs); - - // Resize to the appropriate number of time intervals - m_counts.nmaps(m_gti.size()); - - // Create the peaksigmap - m_peaksigmap = m_counts.extract(1); -} - - -/***********************************************************************//** - * @brief Initialize counts cube + * @exception GException::invalid_value + * Invalid start or stop time or insufficient time bins * - * Sets up the required event cube and time interval storage objects - * associated with this object + * Initialise Good Time Intervals for the variability search. ***************************************************************************/ void ctfindvar::init_gtis(void) { - double tinterval = (*this)["tinterval"].real(); - m_tstart = GTime("JD 999999999"); // large unreasonable Julian date - m_tstop = GTime("JD 0"); // small unreasonable Julian date + // Log header + log_header1(NORMAL, "Initialise Good Time Intervals"); - // Set the start time - if ((*this)["tmin"].is_valid()) { - log_string(NORMAL, "Setting start time from command line"); - m_tstart = (*this)["tmin"].time(); - } else { - log_string(NORMAL, "Setting start time from observations"); - for (int o=0; o(m_obs[o]); - - // Check that the start time of the run is less than tstart - if (obs->gti().tstart() < m_tstart) { - m_tstart = obs->gti().tstart(); - } - } - } + // Clear time intervals + m_gti.clear(); - // Set the stop time - if ((*this)["tmax"].is_valid()) { - log_string(NORMAL, "Setting stop time from command line"); - m_tstop = (*this)["tmax"].time(); - } else { + // Set start and stop times + m_tstart = get_tstart(); + m_tstop = get_tstop(); - log_string(NORMAL, "Setting stop time from observations"); - for (int o=0; o(m_obs[o]); - - // Check that the stop time of the run is larger than tstop - if (obs->gti().tstop() > m_tstop) { - m_tstop = obs->gti().tstop(); - } - } - } + // Query time interval + double tinterval = (*this)["tinterval"].real(); - // Fill the number of good time intervals + // Compute time bins double tstart_sec = m_tstart.secs(); double tstop_sec = m_tstop.secs(); int bins = (tstop_sec - tstart_sec) / tinterval + 0.5; // Log information - log_value(NORMAL, "treference (mjd)", m_tstart.reference().mjdref()); - log_value(NORMAL, "tstart (sec)", tstart_sec); - log_value(NORMAL, "tstop (sec)", tstop_sec); - log_value(NORMAL, "Total time", tstop_sec-tstart_sec); - log_value(NORMAL, "Time interval", tinterval); - log_value(NORMAL, "nbins", bins); - - // Throw if tstart == tstop + log_value(NORMAL, "Reference time (mjd)", m_tstart.reference().mjdref()); + log_value(NORMAL, "Start time (sec)", tstart_sec); + log_value(NORMAL, "Stop time (sec)", tstop_sec); + log_value(NORMAL, "Total time (sec)", tstop_sec-tstart_sec); + log_value(NORMAL, "Time interval (sec)", tinterval); + log_value(NORMAL, "Number of time bins", bins); + + // Throw exception if tstart == tstop if (m_tstart == m_tstop) { throw GException::invalid_value(G_INIT_GTIS, "Start time is equal to stop time"); - } else if (bins < 2) { + } + else if (bins < 2) { std::string msg = "Method requires at least two time bins (" + gammalib::str(bins) + " bins found)."; throw GException::invalid_value(G_INIT_GTIS, msg); } - for (int i=0; i(m_obs[o]); - - if (gti_overlap(obs->gti(), next_gti) > 0.0) { - m_gti.push_back(next_gti); - break; - } - } + for (int k = 0; k < m_obs.size(); ++k) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[k]); + + // Continue only if observation is valid + if (obs != NULL) { + + // Add GTI is the time intervals overlap + if (gti_overlap(obs->gti(), gti) > 0.0) { + m_gti.push_back(gti); + break; + } + + } // endif: observation was valid + + } // endfor: looped over observations + + } // endfor: looped over Good Time Intervals + + // Return + return; +} + + +/***********************************************************************//** + * @brief Create counts cube + * + * Creates a counts cube comprising a sky map for each Good Time Interval. + ***************************************************************************/ +void ctfindvar::create_cube(void) +{ + // Log header + log_header1(NORMAL, "Create counts cube"); + + // Create the basic skymap + m_counts = create_map(m_obs); + + // Resize to the appropriate number of time intervals + m_counts.nmaps(m_gti.size()); + + // Create the peaksigmap + m_peaksigmap = m_counts.extract(1); + + // Loop over all unbinned CTA observations in the container + #pragma omp parallel for + for (int i = 0; i < m_obs.size(); ++i) { + + // Get pointer to observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Fill the cube + fill_cube(obs); + + // Dispose events to free memory + obs->dispose_events(); + + } // endfor: looped over observations + + // Smooth the maps if requested + if ((*this)["smooth_kernel"].is_valid() && + (*this)["smooth_kernel"].string() != "NONE") { + m_counts.smooth((*this)["smooth_kernel"].string(), + (*this)["smooth_rad"].real()); } + + // Return + return; } /***********************************************************************//** - * @brief Fill the cube from the events in the map + * @brief Fill the cube from the events in an observation * - * Sets up the required event cube and time interval storage objects - * associated with this object + * @param[in] obs Pointer to CTA observation. + * + * Fills the cube from the events in an observation. ***************************************************************************/ void ctfindvar::fill_cube(GCTAObservation* obs) { @@ -733,10 +662,6 @@ void ctfindvar::fill_cube(GCTAObservation* obs) int num_outside_time = 0; int num_in_map = 0; - // Specify whether we are cutting on energy or not - bool cut_emin = (m_emin.TeV() > 0.0) ? true : false; - bool cut_emax = (m_emax.TeV() > 0.0) ? true : false; - // Fill counts sky map for (int i = 0; i < events->size(); ++i) { @@ -744,16 +669,14 @@ void ctfindvar::fill_cube(GCTAObservation* obs) const GCTAEventAtom* event = (*events)[i]; // Check that the energy is valid - if ((cut_emin && event->energy() < m_emin) || - (cut_emax && event->energy() > m_emax)) { + if ((event->energy() < m_emin) || (event->energy() > m_emax)) { num_outside_ecut++; continue; } - // Determine event time - GTime evnt_time = event->time(); - int time_indx = time2inx(evnt_time); - + // Determine event time index + int time_indx = time2inx(event->time()); + // Check if event is in the GTIs if ((time_indx > -1) && (time_indx < m_counts.nmaps())) { @@ -791,10 +714,11 @@ void ctfindvar::fill_cube(GCTAObservation* obs) // Increment number of maps num_in_map++; - } else { - // Event falls outside the specified time range + } + + // ... otherwise, event falls outside the specified time range + else { num_outside_time++; - continue; } } // endfor: looped over all events @@ -818,242 +742,557 @@ void ctfindvar::fill_cube(GCTAObservation* obs) /***********************************************************************//** - * @brief Save the current sky direction to model if significance > thr - * - * @param[in] dir_pix pixel sky direction - * @param[out] m_model_sig_above_thr model position with sig>thr + * @brief Analyse all pixels of counts cube ***************************************************************************/ -void ctfindvar::fill_model_sig_pos(const GSkyDir& dir_pix) +void ctfindvar::analyse_cube(void) { - //Set spatial components - GModelSpatialPointSource spatial(dir_pix); + // Log header + log_header1(NORMAL, "Analyse sky map pixels"); - //Set fake spectral component (powerlaw with 1% of Crab) - GModelSpectralPlaw spectral(5.7e-18, -2.48, GEnergy(0.3, "TeV")); + // Get number of GTIs + const int nbins = m_gti.size(); - //Creat a sky model - GModelSky model(spatial, spectral); - - //Extracting the coordinate to name the model - double ra = dir_pix.ra_deg(); - double dec = dir_pix.dec_deg(); - std::string name = "Sky_direction_" + gammalib::str(ra, 4) + - "_" + gammalib::str(dec, 4); - model.name(name); - - //Add the model for this pixel to the output mode - m_model_above_thr.append(model); + // Get relevant sky map pixels + std::vector srcInxPix = get_pixels(); + + // Preparing the histogram with significance evolution for the source + // center and the highest sig-pix + m_pixsigsrc = GNdarray(srcInxPix.size(), nbins); + m_pixsigmax = GNdarray(nbins); + + // Initialise maximum significance + double max_sig = 0; + + // looping over all the pixels in the cube + #pragma omp parallel for + for (int ipix = 0; ipix < m_counts.npix(); ++ipix) { + + // Compute total counts for pixel + double total_counts = 0.0; + for (int k = 0; k < nbins; ++k) { + total_counts += m_counts(ipix, k); + } + + // Getting the variability significance for the current pixel + GNdarray pixSig = get_variability_sig(ipix); + + // Getting the significance evolution for the source + for (int src = 0; src < srcInxPix.size(); ++src) { + + // Store the distribution if the source is located at this position + if (srcInxPix[src] == ipix) { + + #pragma omp critical(ctfindvar_analyse_cube) + for (int i = 0; i < nbins; i++) { + m_pixsigsrc(src,i) = pixSig(i); + } + } + + } // endfor: looped over sources + + // Make sure that this code is not executed in parallel + #pragma omp critical(ctfindvar_analyse_cube) + { + // Determine the pixel with the highest significance + if (max(pixSig) > max_sig) { + max_sig = max(pixSig); + m_max_sig_dir = m_counts.inx2dir(ipix); + m_pixsigmax = pixSig; + } + + // Storing sig in skymap + m_peaksigmap(ipix) = max(pixSig); + + // If the significance is greater than the significance threshold, + // the position is stored in a model + if (max(pixSig) > m_sig_threshold) { + GSkyDir dir_pix = m_counts.inx2dir(ipix); + m_model_above_thr.append(sky_model(dir_pix)); + } + } + } // endfor: looped over pixels + + // Log results + log_value(NORMAL, "Maximum significance", max_sig); + log_value(NORMAL, "Right Ascension of maximum", m_max_sig_dir.ra_deg()); + log_value(NORMAL, "Declination of maximum", m_max_sig_dir.dec_deg()); + + // Return + return; } + + /***********************************************************************//** - * @brief Get the map index associated with a given time + * @brief Return pixel vector + * + * @return Vector of sky map pixels * - * @param[in] pix_number Spatial pixel - * @param[out] alpha_vector Vector containing effective exposure + * Returns a vector of sky map pixels that corresponds to the source model + * positions. ***************************************************************************/ -void ctfindvar::fill_alpha_vector(const int& pix_number, - std::vector& alpha_vector) +std::vector ctfindvar::get_pixels(void) { - // Pixel position - GSkyDir pix_dir = m_counts.inx2dir(pix_number); + // Initialise pixel array + std::vector srcInxPix; - // Loop over all observations - for (int i = 0; i < m_obs.size(); i++) { - // Skip if observation does not overlap with this pixel position - GCTAObservation* obs = dynamic_cast(m_obs[i]); - GCTARoi roi = obs->roi(); - GSkyDir roi_centre = roi.centre().dir(); + // If there are input models then extract + if (m_inmodel.size() > 0) { - if (roi_centre.dist_deg(pix_dir) > roi.radius()) { - continue; - } - - // Convert sky direction to instrument direction - GCTAInstDir instdir = obs->pointing().instdir(pix_dir); + // Loop over models + for (int i = 0; i < m_inmodel.size(); ++i) { - // Extract the good time intervals of the observation - const GGti& obs_gti(obs->gti()); + // Extract the model + GModelSky* model = dynamic_cast(m_inmodel[i]); - // Loop over all time bins - for (int j = 0; j < m_gti.size(); j++) { + // If model is a sky model then extract the model centre. Note + // that this only works for a model that results a GSkyRegionCircle + // object as region. + if (model != NULL) { - // Make sure observation overlaps with this time interval - double exposure = gti_overlap(m_gti[j], obs_gti); - if (exposure > 0.0) { + // Get the source position + GModelSpatial* spatial = model->spatial(); + GSkyRegion* region = spatial->region(); + GSkyRegionCircle* circle = dynamic_cast(region); + + // If region circle is valid then + if (circle != NULL) { + int ipix = m_counts.dir2inx(circle->centre()); + srcInxPix.push_back(ipix); + continue; + } - // Get IRF response - const GCTAResponseIrf* rsp = dynamic_cast - (obs->response()); + } // endif: model was a sky model - // Throw an exception if observation has no instrument response function - if (rsp == NULL) { - std::string msg = - "No response information available for "+ - get_obs_header(obs)+" to compute IRF background. " - "Please specify response information or use " - "another background subtraction method."; - throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg); - } + // If we are still alive then remove the model since it is neither + // a sky model, nor does it have a region circle + m_inmodel.remove(i); + i--; - // Get IRF background template - const GCTABackground* bkg = rsp->background(); + } // endfor: looped over all models - // Throw an exception if observation has no IRF background template - if (bkg == NULL) { - std::string msg = - "No IRF background template found in instrument " - "response function for "+ - get_obs_header(obs)+". Please specify an instrument " - "response function containing a background template."; - throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg); - } + } // endif: there were input models + + // ... otherwise query the source position + else { - // Multiply background rate with livetime and solid angle - #pragma omp critical(ctfindvar_alpha_vector) - alpha_vector[j] += exposure * bkg->rate_ebin(instdir, m_emin, m_emax); - } + // Query source direction + GSkyDir dir; + if ((*this)["coordsys"].string() == "CEL") { + dir.radec_deg((*this)["xsrc"].real(),(*this)["ysrc"].real()); + } + else { + dir.lb_deg((*this)["xsrc"].real(),(*this)["ysrc"].real()); } - } - return; -} + // Store the index of this source + int ipix = m_counts.dir2inx(dir); + srcInxPix.push_back(ipix); + + } // endelse: queried source position + + // Return pixel array + return srcInxPix; +} + + +/***********************************************************************//** + * @brief Return sky model for a given sky direction + * + * @param[in] dir Sky direction + * @return Sky model for a given sky direction + * + * Returns a point source sky model at the position of the current pixel + * with a power law spectral model corresponding to 1% of the Crab flux. + * + * The model name is built from the celestial coordinates of the source. + ***************************************************************************/ +GModelSky ctfindvar::sky_model(const GSkyDir& dir) const +{ + // Set point source spatial component from sky direction + GModelSpatialPointSource spatial(dir); + // Set fake spectral component (powerlaw with 1% of Crab) + GModelSpectralPlaw spectral(5.7e-18, -2.48, GEnergy(0.3, "TeV")); + // Create a sky model + GModelSky model(spatial, spectral); + + // Build model name from celestial coordinates + double ra = dir.ra_deg(); + double dec = dir.dec_deg(); + std::string name = gammalib::str(ra, 4); + if (dec < 0.0) { + name += "-"; + } + else { + name += "+"; + } + name += gammalib::str(std::abs(dec), 4); + + // Set model name + model.name(name); + + // Return model + return model; +} /***********************************************************************//** - * @brief Obtain the significance of variability in a given bin + * @brief Get the significance of variability for a given skymap pixel * - * @param[in] pix_number Pixel in the skymap - * @param[in] nbins Number of bins in the skymap - * @param[out] sig_histogram Histogram for storing returned significances + * @param[in] ipix Sky map pixel index. + * @return Histogram of significances. * * Significance is computed according to Li & Ma equation 17: * see: https://doi.org/10.1086/161295 * Some modification is made in order to handle the case where Non = 0, * but Noff != 0. ***************************************************************************/ -void ctfindvar::get_variability_sig(const int& pix_number, - const int& nbins, - GNdarray& sig_histogram) +GNdarray ctfindvar::get_variability_sig(const int& ipix) { - std::vector accepted_bin_bckg_vector(nbins, true); - bool background_validated = false; - double non, noff; - double alpha, sig; + // Initialise result + GNdarray sig_histogram(m_gti.size()); + + // Initialise vectors + std::vector accepted_bin_bckg_vector(m_gti.size(), true); - // Fill the alpha values for each vector - std::vector alpha_vector(nbins, 0.0); - fill_alpha_vector(pix_number, alpha_vector); + // Get alpha values for specified pixel + std::vector alphas = get_alphas(ipix); - // Loop over pixels until all background pixls are removed + // Loop over pixels until all background pixels are removed + bool background_validated = false; while (background_validated == false) { + + // Signal that background was validated background_validated = true; // Loop over all the GTIs of the pixel - for (int i = 0; i < nbins; i++) { + for (int i = 0; i < m_gti.size(); ++i) { - // The GTI is discared from bckg calculation and not checked again. + // The GTI is discared from background calculation and not + // checked again if (!accepted_bin_bckg_vector[i]) { continue; - } - // Check if bin fails minoff check or no observations overlap it - else if (m_counts(pix_number, i) < m_minoff || alpha_vector[i] == 0) { + } + + // ... otherwise check if bin fails minoff check or no observations + // overlap it + else if (m_counts(ipix, i) < m_minoff || alphas[i] == 0) { accepted_bin_bckg_vector[i] = false; continue; } - noff = 0.0; - alpha = 0.0; - // Gor one GTI selected (i), looping over all the others (j). - for (int j = 0; j < nbins; j++) { + // ... otherwise + double noff = 0.0; + double alpha = 0.0; + + // Got one GTI selected (i), looping over all the others (j). + for (int j = 0; j < m_gti.size(); ++j) { if (j != i && accepted_bin_bckg_vector[j] == 1) { - noff += m_counts(pix_number, j); - alpha += alpha_vector[j]; + noff += m_counts(ipix, j); + alpha += alphas[j]; } } // The background is averaged on the number of bins -1 - non = m_counts(pix_number, i); - alpha = alpha_vector[i]/alpha; + double non = m_counts(ipix, i); + alpha = alphas[i]/alpha; // Compute sensitivity in Gaussian sigma using Li & Ma double alpha1 = alpha + 1.0; double ntotal = non+noff; double arg1 = non/ntotal; double arg2 = noff/ntotal; - if (noff == 0.0) { - sig = 0.0; - } else if (non > 0.0) { - double term1 = non * std::log((alpha1/alpha)*arg1); - double term2 = noff * std::log(alpha1*arg2); - sig = std::sqrt(2.0 * (term1 + term2)); - } else { - sig = std::sqrt(2.0 * noff * std::log(alpha1)); - } - - // Specify the sign of the significance - sig *= (non < alpha*noff) ? -1.0 : 1.0; - - if (alpha == 0.0) { - sig = 0.0; + double sig = 0.0; + if (alpha != 0.0) { + if (noff == 0.0) { + sig = 0.0; + } + else if (non > 0.0) { + double term1 = non * std::log((alpha1/alpha)*arg1); + double term2 = noff * std::log(alpha1*arg2); + sig = std::sqrt(2.0 * (term1 + term2)); + } + else { + sig = std::sqrt(2.0 * noff * std::log(alpha1)); + } + sig *= (non < alpha*noff) ? -1.0 : 1.0; } // Update the significance sig_histogram(i) = sig; - // If bin is significant, it's removed from bckg and we loop again. + // If bin is significant, it's removed from the background and + // we loop again if (sig > m_sig_threshold) { accepted_bin_bckg_vector[i] = false; background_validated = false; } - } - } + } // endfor: looped over all GTIs + + } // endwhile: iterate until converged + + // Return significance histogram + return sig_histogram; } /***********************************************************************//** + * @brief Get alpha vector + * + * @param[in] ipix Sky map pixel index + * @return Vector containing effective exposure + ***************************************************************************/ +std::vector ctfindvar::get_alphas(const int& ipix) const +{ + // Initialise alpha vector + std::vector alphas(m_gti.size(), 0.0); + + // Convert pixel index into sky direction + GSkyDir pix_dir = m_counts.inx2dir(ipix); + + // Loop over all observations + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + const GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Fall through if observation is not valid + if (obs == NULL) { + continue; + } + + // Get RoI. Fall through if observation does not overlap with this + // pixel position + GCTARoi roi = obs->roi(); + GSkyDir roi_centre = roi.centre().dir(); + if (roi_centre.dist_deg(pix_dir) > roi.radius()) { + continue; + } + + // Convert sky direction to instrument direction + GCTAInstDir instdir = obs->pointing().instdir(pix_dir); + + // Extract the good time intervals of the observation + const GGti& obs_gti(obs->gti()); + + // Loop over all time bins + for (int j = 0; j < m_gti.size(); ++j) { + + // Make sure observation overlaps with this time interval + double exposure = gti_overlap(m_gti[j], obs_gti); + if (exposure > 0.0) { + + // Get IRF response + const GCTAResponseIrf* rsp = + dynamic_cast(obs->response()); + + // Throw an exception if observation has no instrument response function + if (rsp == NULL) { + std::string msg = + "No response information available for "+ + get_obs_header(obs)+" to compute IRF background. " + "Please specify response information or use " + "another background subtraction method."; + throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg); + } + + // Get IRF background template + const GCTABackground* bkg = rsp->background(); + + // Throw an exception if observation has no IRF background template + if (bkg == NULL) { + std::string msg = + "No IRF background template found in instrument " + "response function for "+ + get_obs_header(obs)+". Please specify an instrument " + "response function containing a background template."; + throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg); + } + + // Add up background rate + alphas[j] += exposure * bkg->rate_ebin(instdir, m_emin, m_emax); + + } // endif: exposure was positive + + } // endfor: looped over all time bins + + } // endfor: looped over all observations + + // Return alphas + return alphas; +} + + +/***********************************************************************//** * @brief Returns number of seconds that two GTIs overlap * * @param[in] gti1 First GTI * @param[in] gti2 Second GTI * @return Number of seconds that @p gti1 and @p gti2 overlap ***************************************************************************/ -double ctfindvar::gti_overlap(const GGti& gti1, const GGti& gti2) +double ctfindvar::gti_overlap(const GGti& gti1, const GGti& gti2) const { + // Initialise no overlap double overlap = 0.0; + + // Set first and second GTI so that gti_1st starts before gti_2nd const GGti& gti_1st = (gti1.tstart() <= gti2.tstart()) ? gti1 : gti2; - const GGti& gti_2nd = (gti1.tstart() > gti2.tstart()) ? gti1 : gti2; + const GGti& gti_2nd = (gti1.tstart() > gti2.tstart()) ? gti1 : gti2; - // gti1 starts earlier + // Compute overlap if (gti_1st.tstop() > gti_2nd.tstart()) { if (gti_1st.tstop() <= gti_2nd.tstop()) { overlap = gti_1st.tstop() - gti_2nd.tstart(); - } else if (gti_1st.tstop() > gti_2nd.tstop()) { + } + else if (gti_1st.tstop() > gti_2nd.tstop()) { overlap = gti_2nd.tstop() - gti_2nd.tstart(); } } + // Return overlap return overlap; } /***********************************************************************//** - * @brief Write individual histograms in CSV format + * @brief Get start time + * + * @return Start time + ***************************************************************************/ +GTime ctfindvar::get_tstart(void) +{ + // Initialise start time with large unreasonable Julian date + GTime tstart("JD 999999999"); + + // If tmin parameter is valid then set the start time + if ((*this)["tmin"].is_valid()) { + + // Query tmin parameter + tstart = (*this)["tmin"].time(); + + // Log action + log_string(NORMAL, "Setting start time \"tmin\" parameter"); + + } // endif: tmin parameter was valid + + // ... otherwise set start time from observations + else { + + // Loop over all observations + for (int k = 0; k < m_obs.size(); ++k) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[k]); + + // Continue only if observation is valid + if (obs != NULL) { + + // Update start time if the start time of the observation is + // smaller than the current start time + if (obs->gti().tstart() < tstart) { + tstart = obs->gti().tstart(); + } + + } // endif: observation was valid + + } // endfor: looped over observations + + // Log action + log_string(NORMAL, "Setting start time from observations"); + + } // endelse: no valid time specified + + // Return start time + return tstart; +} + + +/***********************************************************************//** + * @brief Get stop time + * + * @return Stop time ***************************************************************************/ -void ctfindvar::write_srchist(void) +GTime ctfindvar::get_tstop(void) +{ + // Initialise stop time with small unreasonable Julian date + GTime tstop("JD 0"); + + // If tmax parameter is valid then set the stop time + if ((*this)["tmax"].is_valid()) { + + // Query tmax parameter + tstop = (*this)["tmax"].time(); + + // Log action + log_string(NORMAL, "Setting stop time \"tmax\" parameter"); + + } // endif: tmax parameter was valid + + // ... otherwise set stop time from observations + else { + + // Loop over all observations + for (int k = 0; k < m_obs.size(); ++k) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[k]); + + // Continue only if observation is valid + if (obs != NULL) { + + // Update stop time if the stop time of the observation is + // larger than the current stop time + if (obs->gti().tstop() > tstop) { + tstop = obs->gti().tstop(); + } + + } // endif: observation was valid + + } // endfor: looped over observations + + // Log action + log_string(NORMAL, "Setting stop time from observations"); + + } // endelse: no valid time specified + + // Return stop time + return tstop; +} + + +/***********************************************************************//** + * @brief Write source histograms in FITS format + * + * @param[in,out] fits FITS file + * + * This method creates two FITS binary table extensions and appends them to + * the FITS file. + * + * The first extension is named "SOURCE SIGNIFICANCE" and contains the + * significance of the variability as function of time bin for the pixel + * with the maximum significance and for each source. + * + * The second extension is named "SOURCE POSITION" and contains the Right + * Ascension and Declination for the maximum significance pixel and for each + * source. + ***************************************************************************/ +void ctfindvar::write_source_histograms(GFits& fits) { // Setup data storage for time and source significances double tinterval = (*this)["tinterval"].real(); - int nentries = (m_tstop-m_tstart) / tinterval + 0.5; - int nsources = m_pixsigsrc.shape()[0]; - GNdarray time_info(2, nentries); - GNdarray max_pixel_info(nentries); - GNdarray src_info(nsources, nentries); - - // Store the full set of information - for (int i=0; i= 0) { - src_info(src,i) = m_pixsigsrc(src,index); - } else { - src_info(src,i) = 0.0; - } + for (int k = 0; k < nsrc; ++k) { + src_info(k,i) = (index >= 0) ? m_pixsigsrc(k,index) : 0.0; } - // Fill max pixel information - if (index >= 0) { - max_pixel_info(i) = m_pixsigmax(index); - } + // Fill maximum pixel information + max_pixel_info(i) = (index >= 0) ? m_pixsigmax(index) : 0.0; + + } // endfor: looped over time bins + + // Create binary tables + GFitsBinTable table_signif(nbins); + GFitsBinTable table_pos(nsrc+1); + table_signif.extname("SOURCE SIGNIFICANCE"); + table_pos.extname("SOURCE POSITION"); + + // Create columns for time bins + GFitsTableDoubleCol time_start("TSTART", nbins); + GFitsTableDoubleCol time_stop("TSTOP", nbins); + GFitsTableDoubleCol time_sigma("MAXSIGPIXEL", nbins); + time_start.unit("MJD"); + time_stop.unit("MJD"); + time_sigma.unit("sigma"); + + // Create columns for source position + GFitsTableStringCol pos_name("SOURCE", nsrc+1, 20); + GFitsTableDoubleCol pos_ra("RA", nsrc+1); + GFitsTableDoubleCol pos_dec("DEC", nsrc+1); + pos_ra.unit("deg"); + pos_dec.unit("deg"); + + // Fill time bin columns + for (int i = 0; i < nbins; ++i) { + time_start(i) = time_info(0,i); + time_stop(i) = time_info(1,i); + time_sigma(i) = max_pixel_info(i); } - // Call the appropriate method for writing the data - write_srchist_fits(time_info, max_pixel_info, src_info); -} + // Fill source columns + pos_name(0) = "MAXSIGPIXEL"; + pos_ra(0) = m_max_sig_dir.ra_deg(); + pos_dec(0) = m_max_sig_dir.dec_deg(); - -/***********************************************************************//** - * @brief Write individual histograms in FITS format - * - * @param[in] time_info Start and stop time in MJD for each bin - * @param[in] src_info Significance distribution for each source - - * This method stores the extracted significance for each source at each time - * interval between tstart and tstop. Each time bin contains information on - * the start and stop time of that bin (in MJD) and the extracted variability - * significance of the source in that time bin. - ***************************************************************************/ -void ctfindvar::write_srchist_fits(const GNdarray& time_info, - const GNdarray& max_pixel_info, - const GNdarray& src_info) -{ - log_string(NORMAL, "Storing histograms in FITS format"); - - // Create the table - int nrows = time_info.shape()[1]; - GFitsBinTable table_signif(nrows); - GFitsBinTable table_dir(2); - table_signif.extname("SRCSIGNIF"); - table_dir.extname("SRCPOSITION"); - - // Append the start/stop time and max pixel columns - GFitsTableDoubleCol start_mjd("TSTART", nrows); - GFitsTableDoubleCol stop_mjd("TSTOP", nrows); - GFitsTableDoubleCol max_pix("MAXSIGPIXEL", nrows); - for (int i=0; i 0) { - src_name = m_inmodel[src]->name(); - - // Get source position - GModelSky* model = dynamic_cast(m_inmodel[src]); - GModelSpatial* model_spatial = model->spatial(); - GSkyRegion* model_region = model_spatial->region(); - src_dir = dynamic_cast(model_region)->centre(); - } - // ... otherwise set the position from xsrc,ysrc + if (m_inmodel.size() > k) { + + // Extract the model + GModelSky* model = dynamic_cast(m_inmodel[k]); + + // If model is a sky model then extract the model centre. Note + // that this only works for a model that results a GSkyRegionCircle + // object as region. + if (model != NULL) { + + // Get the source position + GModelSpatial* spatial = model->spatial(); + GSkyRegion* region = spatial->region(); + GSkyRegionCircle* circle = dynamic_cast(region); + + // If region circle is valid then extract source name and + // sky direction + if (circle != NULL) { + src_name = m_inmodel[k]->name(); + src_dir = circle->centre(); + } + + } // endif: model was a sky model + + } // endif: there were enough models in the input model + + // ... otherwise set the position from the xsrc and ysrc parameters else { if ((*this)["coordsys"].string() == "CEL") { - src_dir.radec_deg((*this)["xsrc"].real(), (*this)["yref"].real()); - } else { - src_dir.lb_deg((*this)["xsrc"].real(), (*this)["yref"].real()); + src_dir.radec_deg((*this)["xsrc"].real(), (*this)["ysrc"].real()); + } + else { + src_dir.lb_deg((*this)["xsrc"].real(), (*this)["ysrc"].real()); } } - - // Create and fill a new column - GFitsTableDoubleCol src_sig_col(src_name, nrows); - for (int i=0; i& alpha_vector); - void get_variability_sig(const int& pix_number, - const int& nbins, - GNdarray& sig_histogram); - double gti_overlap(const GGti& gti1, - const GGti& gti2); - - // Writing methods for individual source histograms - void write_srchist(void); - void write_srchist_fits(const GNdarray& time_info, - const GNdarray& max_pixel_info, - const GNdarray& src_info); + void init_members(void); + void copy_members(const ctfindvar& app); + void free_members(void); + void get_parameters(void); + void init_gtis(void); + void create_cube(void); + void fill_cube(GCTAObservation* obs); + void analyse_cube(void); + std::vector get_pixels(void); + GNdarray get_variability_sig(const int& ipix); + std::vector get_alphas(const int& ipix) const; + double gti_overlap(const GGti& gti1, const GGti& gti2) const; + GTime get_tstart(void); + GTime get_tstop(void); + GModelSky sky_model(const GSkyDir& dir) const; + void write_source_histograms(GFits& fits); // Protected members GSkyMap m_counts; //!< Counts for each time interval @@ -97,7 +92,6 @@ protected: GSkyDir m_max_sig_dir; //!< Sky direction associated with maximum significance double m_minoff; //!< Minimum counts for use in significance calculation double m_sig_threshold; //!< Minimum significance required to set source as variable - GFits m_outfile; //!< FITS object to store output GSkyMap m_peaksigmap; //!< Skymap holding the maximum significance GNdarray m_pixsigsrc; //!< Store distributions of the source significances GNdarray m_pixsigmax; //!< Store distribution for pixel with max significance diff --git a/src/ctfindvar/ctfindvar.par b/src/ctfindvar/ctfindvar.par index 40ce4d5..06a0da9 100644 --- a/src/ctfindvar/ctfindvar.par +++ b/src/ctfindvar/ctfindvar.par @@ -16,22 +16,26 @@ # # General parameters #=================== -inobs, f, a, events.fits,,, "Input event list, counts cube or observation definition XML file" +inobs, f, a, events.fits,,, "Input event list or observation definition XML file" inmodel, f, h, NONE,,, "Input model definition file for extracting source positions" -outcube, f, h, NONE,,, "Filename for counts cube" +outcube, f, h, NONE,,, "Output counts cube file" +outmap, f, a, map.fits,,, "Output map file" +outmodel, f, a, model.xml,,, "Output model definition file" caldb, s, a, prod2,,, "Calibration database" irf, s, a, South_0.5h,,, "Instrument response function" -prefix, s, a, findvar_,,, "Prefix for output files" # # Variability search parameters #============================== -minoff, r, h, 0,,, "Minimum required events in off bin for computing Noff and alpha" -emin, r, a, 0,,, "Lower energy limit (units of TeV; 0 => no cut)" -emax, r, a, 0,,, "Upper energy limit (units of TeV; 0 => no cut)" -sig_thr, r, h, 4.5,,, "Significance threshold for variability detection" -smoothkrnl, s, h, NONE, GAUSSIAN|DISK|NONE,, "Smoothing kernel for counts maps" -smoothpar, r, h, 0.05,,, "Parameter for the smoothing kernel" +coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +xsrc, r, a, 83.63,0,360, "First coordinate of the source of interest (RA or galactic l)" +ysrc, r, a, 22.01,-90,90, "Second coordinate of the source of interest (DEC or galactic b)" +emin, r, a, 0.1,,, "Lower energy limit (TeV)" +emax, r, a, 100.0,,, "Upper energy limit (TeV)" +threshold, r, h, 4.5,,, "Significance threshold for variability detection" +minoff, r, h, 0,,, "Minimum required events in off bin for computing Noff and alpha" +smooth_kernel, s, h, NONE, GAUSSIAN|DISK|NONE,, "Smoothing kernel for counts maps" +smooth_rad, r, h, 0.05,,, "Smoothing kernel radius (deg)" # # Time binning parameters @@ -47,12 +51,9 @@ usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" -coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" proj, s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method" xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" yref, r, a, 22.51,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" -xsrc, r, a, 83.63,0,360, "First coordinate of the source of interest (RA or galactic l)" -ysrc, r, a, 22.51,-90,90, "Second coordinate of the source of interest (DEC or galactic b)" # # Standard parameters diff --git a/test/test_ctfindvar.py b/test/test_ctfindvar.py index d20dfd0..0c0c4da 100755 --- a/test/test_ctfindvar.py +++ b/test/test_ctfindvar.py @@ -71,8 +71,9 @@ class Test(test): # Setup ctfindvar command cmd = ctfindvar+' inobs="'+self._events+'"'+\ + ' outmap="ctfindvar_cmd1.fits"'+\ + ' outmodel="ctfindvar_cmd1.xml"'+\ ' logfile="ctfindvar_cmd1.log" chatter=1'+\ - ' prefix=ctfindvar_test_ '+\ ' tmin=NONE tmax=NONE tinterval=20'+\ ' emin=0.2 emax=1.1'+\ ' nxpix=200 nypix=200 binsz=0.05'+\ @@ -80,12 +81,28 @@ class Test(test): ' xref=83.63 yref=22.51'+\ ' xsrc=83.63 ysrc=22.51'+\ ' caldb=prod2 irf=South_0.5h' - - + # Check if execution was successful self.test_value(self._execute(cmd), 0, 'Check successful execution from command line') + # Setup ctfindvar command + cmd = ctfindvar+' inobs="event_file_that_does_not_exist.fits"'+ \ + ' outmap="ctfindvar_cmd2.fits"'+ \ + ' outmodel="ctfindvar_cmd2.xml"'+ \ + ' logfile="ctfindvar_cmd2.log" chatter=1'+ \ + ' tmin=NONE tmax=NONE tinterval=20'+ \ + ' emin=0.2 emax=1.1'+ \ + ' nxpix=200 nypix=200 binsz=0.05'+ \ + ' proj=CAR coordsys=CEL'+ \ + ' xref=83.63 yref=22.51'+ \ + ' xsrc=83.63 ysrc=22.51'+ \ + ' caldb=prod2 irf=South_0.5h' + + # Check if execution was successful + self.test_assert(self._execute(cmd, success=False) != 0, + 'Check invalid input file when executed from command line') + # Check ctfindvar --help self._check_help(ctfindvar) @@ -101,11 +118,8 @@ class Test(test): # Allocate ctfindvar ctfindvar = ctools.ctfindvar() - #setup ctfindvar + # Setup ctfindvar ctfindvar['inobs'] = self._events - ctfindvar['logfile'] = 'ctfindvar_py0.log' - ctfindvar['chatter'] = 1 - ctfindvar['prefix'] = 'ctfindvar_test_python_' ctfindvar['tmin'] = 'NONE' ctfindvar['tmax'] = 'NONE' ctfindvar['tinterval'] = 20 @@ -122,32 +136,39 @@ class Test(test): ctfindvar['ysrc'] = 22.51 ctfindvar['caldb'] = 'prod2' ctfindvar['irf'] = 'South_0.5h' + ctfindvar['outmap'] = 'ctfindvar_py0.fits' + ctfindvar['outmodel'] = 'ctfindvar_py0.xml' + ctfindvar['logfile'] = 'ctfindvar_py0.log' + ctfindvar['chatter'] = 1 - ## Run ctfindvar tool + # Run ctfindvar tool ctfindvar.logFileOpen() ctfindvar.run() # Check that saving does not nothing - ctfindvar['logfile'] = 'ctfindvar_py1.log' + ctfindvar['outmap'] = 'ctfindvar_py1.fits' + ctfindvar['outmodel'] = 'ctfindvar_py1.xml' + ctfindvar['logfile'] = 'ctfindvar_py1.log' ctfindvar.logFileOpen() ctfindvar.save() # Check result - self._check_result_file('ctfindvar_test_python_signifmap.fits') + self._check_result_file('ctfindvar_py1.fits') # Copy ctfindvar tool cpy_tool = ctfindvar # Run copy of ctfindvar tool again - cpy_tool['logfile'] = 'ctfindvar_py2.log' - cpy_tool['chatter'] = 3 - cpy_tool['prefix'] = 'ctfindvar_test_copytool_' + cpy_tool['outmap'] = 'ctfindvar_py2.fits' + cpy_tool['outmodel'] = 'ctfindvar_py2.xml' + cpy_tool['logfile'] = 'ctfindvar_py2.log' + cpy_tool['chatter'] = 3 cpy_tool.logFileOpen() cpy_tool.run() cpy_tool.save() # Check result - self._check_result_file('ctfindvar_test_copytool_signifmap.fits') + self._check_result_file('ctfindvar_py2.fits') # Check that clearing does not lead to an exception or segfault ctfindvar.clear() @@ -158,17 +179,16 @@ class Test(test): # Check ctfindvar result def _check_result_file(self, filename): """ - Check content of tool + Check ctfindvar result Parameters ---------- - tool : `~ctools.ctfindvar` - ctfindvar instance + filename : str + Output file """ - print "checking if you amazing file is valid" - # Read variability file - fits = gammalib.GFits(filename) - self.test_value(fits.size(), 3, - 'Check for 3 rows in srcsig file') + # Read variability file + fits = gammalib.GFits(filename) + self.test_value(fits.size(), 3, 'Check for 3 extensions in output file') + # Return return