Luigi Tibaldo / ctools

Commit 0be0a1a4450dd64c86fc9b87dd08a291684a8716

Authored by Jürgen Knödlseder
1 parent 90c2e5d8

Integrate code

Inline Side-by-side
Showing 10 changed files with 940 additions and 650 deletions
ChangeLog
  View file @0be0a1a
1   -2018-11-23
  1 +2018-12-11
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
... ... @@ -7,6 +7,7 @@
7 7 CTA background rates are now per ontime (#2758)
8 8 Apply event bin weighting in ctmodel for stacked computation (#2721)
9 9 Enable ctbin to produce counts cubes for a joint binned analysis (#2714)
  10 + Add ctfindvar tool (#2713)
10 11 Add Jupyter notebooks for H.E.S.S. DR1 (#2712)
11 12 csphagen now creates appropriate model for likelihood fitting (#2711)
12 13 csphagen now supports any kind of background model (#2710)
... ...
1 1 New Features and Important Changes in ctools 1.6.0
2 2  
3   -23 November 2018
  3 +11 December 2018
4 4  
5 5  
6 6 Introduction
... ... @@ -126,6 +126,12 @@ ctexpcube - CTA exposure cube generation
126 126 None
127 127  
128 128  
  129 +ctfindvar - Variability search
  130 +------------------------------
  131 +New: the ctfindvar tool searches for source variability using an On/Off method
  132 +in the time domain (#2713).
  133 +
  134 +
129 135 ctlike - Maximum likelihood fitting
130 136 -----------------------------------
131 137 Add hidden "like_accuracy" parameter to allow control of the accuracy of
... ...
doc/source/admin/release_history/1.6.rst
  View file @0be0a1a
... ... @@ -14,6 +14,7 @@ In particular, this release provides:
14 14 * Support for Xspec spectral analysis
15 15 * Improved computation speed for :ref:`ctskymap` and :ref:`csphagen`
16 16 * Parallelisation of various cscripts
  17 +* Support for source variability searches
17 18 * Support for Python pickeling
18 19 * A Code Generator to facilitate adding ctools and cscripts
19 20  
... ... @@ -102,6 +103,8 @@ Improvements
102 103 New features
103 104 ------------
104 105  
  106 +* [`2713 <https://cta-redmine.irap.omp.eu/issues/2713>`_] -
  107 + Add :ref:`ctfindvar` tool
105 108 * [`2709 <https://cta-redmine.irap.omp.eu/issues/2709>`_] -
106 109 Add code generator
107 110 * [`2694 <https://cta-redmine.irap.omp.eu/issues/2694>`_] -
... ...
doc/source/users/reference_manual/ctfindvar.rst
  View file @0be0a1a
1 1 .. _ctfindvar:
2 2  
3 3 ctfindvar
4   -==============================
  4 +=========
5 5  
6   -ToDo: Describe in a one liner what the tool is doing.
  6 +Search for source variability.
7 7  
8 8  
9 9 Synopsis
10 10 --------
11   -This tool is used to search for source variability.
12   -The variability is computed by binning the arrival time of the photons
13   -and applying the on-off method to each bin, integrating the background
14   -over all the other bins.
15   -A Limit is set to include a bin in the background sample, i.e., its
16   -significance with respect to all the others bins has to be lower than 4.5
17   -The significance of each bin is computed using Li&Ma eq. 17.
18   -A limit on the background sample can also be included, taking only into
19   -account bins that have a number of counts which is higher than a threshold
20   -defined by the user via the "minoff" parameter.
21   -
22   -Variability is calculated for each pixel of the sky given in the input files.
23   -The results of the variability is further stored in the srcsig.fits file,
24   -where the GTIs and the significance of each bin are stored.
25   -Only the sources with a significance > 4.5 sigmas are stored.
26   -
27   -:ref:`ctfindvar` takes in input some observation files.
28   -It calculates a cube with several maps, each one with a count number integrated
29   -over a time period defined by the tinterval parameter, for each pixel.
30   -The tinterval parameter is the time scale over which the variability
31   -will be searched.
32   -
33   -Inputs files:
34   -- observation file xml, or event list .fits
35   -
36   -Output files:
37   -- countscube.fits
38   - provides all the count-skymaps contained in the cube.
39   -- peaksigmap.fits
40   - a skymap of the highest significance for each pixel
41   -- srcsig.fits
42   - a table containing the GTIs together with the significance
43   - for all the sources with a significant variability detected.
44   -
45   -
46   -The significance of a source can be displayed using the
47   -``show_variability_evolution.py`` script in the example forlder
  11 +
  12 +This tool searches for source variability by applying the On/Off method in the
  13 +time domain. The variability of sky map pixels is computed by binning the events
  14 +as function of the photon arrival time and applying the On/Off method to each
  15 +time bin, integrating the background over all the other time bins. Only bins
  16 +with a variabilty significance below a given threshold, defined by the hidden
  17 +``threshold`` parameter, are included in the background estimate. Furthermore,
  18 +only time bins that exceed a minimum number of events, specified by the hidden
  19 +``minoff`` parameter, are considered for the background estimate. The significance
  20 +of each time bin is computed using Li & Ma formula.
  21 +
  22 +:ref:`ctfindvar` will search a grid of sky map pixels for source variability. Using
  23 +the hidden ``smooth_kernel`` and ``smooth_rad`` parameters, the sky map may be
  24 +smoothed before the variability search. The tool will compute sky maps for a
  25 +number of time bins. The duration of each time bin is specified by the ``tinterval``
  26 +parameter, the start and stop time are defined by the ``tmin`` and ``tmax``
  27 +parameters. If the start and stop time are ``NONE``, the start and stop times will
  28 +be extracted from the input event list(s). Note that ``tinterval`` specifies the
  29 +time scale over which the variability will be searched.
  30 +
  31 +The tool will determine the sky map pixel with the maximum variability and store
  32 +the variability histogram and pixel position as ``MAXSIGPIXEL`` in the output file,
  33 +specified using the ``outmap`` parameter. In addition, the tool will determine the
  34 +variability for a number of source positions. If an input model is specified
  35 +using the ``inmodel`` parameter, the source positions will be defined by the centre
  36 +positions of all gamma-ray source components in the input model. If no input
  37 +model is specified, the tool will use the ``xsrc`` and ``ysrc`` parameters to allow
  38 +the specification of a source position. The results will be stored in the output
  39 +file using either the source names found in the model definition file, or the
  40 +name ``SOURCE`` if the source position was specified by the ``xsrc`` and ``ysrc``
  41 +parameters.
  42 +
  43 +All sky map pixels with a variability above the value specified by the ``threshold``
  44 +will also added be as a point source model to the ``outmodel`` file. Each source will
  45 +have a power law as the spectral component with a total flux that corresponds to
  46 +1% of the Crab nebula flux.
  47 +
48 48  
49 49 General parameters
50 50 ------------------
... ... @@ -56,7 +56,17 @@ General parameters
56 56 Input model definition file for extracting source positions.
57 57  
58 58 ``(outcube = NONE) [file]``
59   - Filename for saving counts cube, NONE will result in no cube being saved.
  59 + Filename for saving counts cube, ``NONE`` will result in no cube being saved.
  60 +
  61 +``outmap [file]``
  62 + Filename for output sky map containing the maximum significance for each
  63 + spatial pixel. The file also contains the variability significance as function
  64 + of time for each source of interest, and the pixel with the maximum
  65 + significance. Furthermore, it contains the Right Ascension and Declination
  66 + of each source of interest and the pixel with the maximum significance.
  67 +
  68 +``outmodel [file]``
  69 + Output model definition file containing sources for significant pixels.
60 70  
61 71 ``caldb [string]``
62 72 Calibration database (only required for IRF background subtraction if no
... ... @@ -66,40 +76,45 @@ General parameters
66 76 Instrument response function (only required for IRF background subtraction
67 77 if no response information is provided by ``inobs``).
68 78  
69   -``prefix [string]``
70   - Output file prefix. The method will save two files:
71   - - Skymap containing the maximum significance for each spatial pixel
72   - - File containing the significance as a function of time interval for each
73   - source
74 79  
75   -``(minoff = 0) [real]``
76   - Minimum number of off counts to compute the significance
  80 +Variability search parameters
  81 +-----------------------------
  82 +
  83 +``coordsys <CEL|GAL> [string]``
  84 + Coordinate system (CEL - celestial, GAL - galactic).
  85 +
  86 +``xsrc [real]``
  87 + Right Ascension / Galactic longitude of source of interest (J2000, in degrees)
  88 +
  89 +``ysrc [real]``
  90 + Declination / Galactic latitude of source of interest (J2000, in degrees)
77 91  
78 92 ``emin [real]``
79   - Minimum energy (TeV) for extracting events. A value of 0 results in no cut
80   - being applied.
  93 + Minimum energy (TeV) for extracting events.
81 94  
82 95 ``emax [real]``
83   - Maximum energy (TeV) for extracting events. A value of 0 results in no cut
84   - being applied.
  96 + Maximum energy (TeV) for extracting events.
85 97  
86   -``(sig_thr = 4.5) [real]``
  98 +``(threshold = 4.5) [real]``
87 99 Significance threshold for variability detection.
88 100  
89   -``(smoothkrnl = NONE) <GAUSSIAN|DISK|NONE> [string]``
  101 +``(minoff = 0) [real]``
  102 + Minimum number of off counts to compute the significance.
  103 +
  104 +``(smooth_kernel = NONE) <GAUSSIAN|DISK|NONE> [string]``
90 105 Kernel to be used in smoothing the counts maps before computing
91 106 significances.
92 107  
93   -``(smoothpar = 0.05) [real]``
94   - Parameter for the smoothing kernel. This is either the sigma of the GAUSSIAN
95   - or the radius of the DISK.
  108 +``(smooth_rad = 0.05) [real]``
  109 + Radius of smoothing kernel in degrees. This is either the sigma of the
  110 + ``GAUSSIAN`` or the radius of the ``DISK``.
96 111  
97 112  
98 113 Time binning parameters
99 114 -----------------------
100 115  
101 116 ``tinterval [real]``
102   - interval in seconds for each time step.
  117 + Interval in seconds for each time step.
103 118  
104 119 ``tmin [time]``
105 120 Start time for interval determination (UTC string, JD, MJD or MET in seconds).
... ... @@ -120,7 +135,7 @@ Spatial binning parameters
120 135 --------------------------
121 136  
122 137 ``(usepnt = no) [boolean]``
123   - Use CTA pointing direction for cube centre instead of xref/yref parameters?
  138 + Use CTA pointing direction for sky map centre instead of xref/yref parameters?
124 139  
125 140 ``nxpix [integer]``
126 141 Number of cube bins in Right Ascension or Galactic longitude.
... ... @@ -129,25 +144,16 @@ Spatial binning parameters
129 144 Number of cube bins in Declination or Galactic latitude.
130 145  
131 146 ``binsz [real]``
132   - Cube bin size (in degrees/pixel).
133   -
134   -``coordsys <CEL|GAL> [string]``
135   - Coordinate system (CEL - celestial, GAL - galactic).
  147 + Sky map pixel size (in degrees/pixel).
136 148  
137 149 ``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
138 150 Projection method.
139 151  
140 152 ``xref [real]``
141   - Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
  153 + Right Ascension / Galactic longitude of sky map centre (J2000, in degrees).
142 154  
143 155 ``yref [real]``
144   - Declination / Galactic latitude of cube centre (J2000, in degrees).
145   -
146   -``xsrc [real]``
147   - Right Ascension / Galactic longitude of source of interest (J2000, in degrees)
148   -
149   -``ysrc [real]``
150   - Declination / Galactic latitude of source of interest (J2000, in degrees)
  156 + Declination / Galactic latitude of sky map centre (J2000, in degrees).
151 157  
152 158  
153 159 Standard parameters
... ...
doc/source/users/reference_manual/index.rst
  View file @0be0a1a
... ... @@ -17,7 +17,6 @@ ctools
17 17 .. toctree::
18 18 :maxdepth: 1
19 19  
20   - ctfindvar --- ToDo: Describe what tool is doing <ctfindvar>
21 20 ctbin --- Generates counts cube <ctbin>
22 21 ctbkgcube --- Generates background cube <ctbkgcube>
23 22 ctbutterfly --- Compute butterfly <ctbutterfly>
... ... @@ -25,6 +24,7 @@ ctools
25 24 ctedispcube --- Generates energy dispersion cube <ctedispcube>
26 25 cterror --- Calculates likelihood profile errors <cterror>
27 26 ctexpcube --- Generates exposure cube <ctexpcube>
  27 + ctfindvar --- Search for source variability <ctfindvar>
28 28 ctlike --- Performs maximum likelihood fitting <ctlike>
29 29 ctmapcube --- Generates a map cube <ctmapcube>
30 30 ctmodel --- Computes model counts cube <ctmodel>
... ...
examples/show_variability_evolution.py
  View file @0be0a1a
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# Shows butterfly diagram created with ctbutterfly
  3 +# Shows source variability evolution created with ctfindvar
4 4 #
5   -# Copyright (C) 2014-2017 Simon Bonnefoy
  5 +# Copyright (C) 2018 Simon Bonnefoy
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -43,6 +43,8 @@ def plot_variability_evolution(filename, plotfile, srcname):
43 43 Name of spectrum FITS file
44 44 plotfile : str
45 45 Plot file name
  46 + srcname : str
  47 + Source name
46 48 """
47 49 # Read variability file
48 50 fits = gammalib.GFits(filename)
... ... @@ -51,41 +53,45 @@ def plot_variability_evolution(filename, plotfile, srcname):
51 53 c_gti_stop = table['TSTOP']
52 54 c_sig = table[srcname]
53 55  
54   - #convert the c_sig object to list
55   - table_c_sig=[]
56   - table_c_gti=[]
57   -
  56 + # Convert the c_sig object to list
  57 + table_c_sig = []
  58 + table_c_gti = []
58 59 for k in c_sig:
59 60 table_c_sig.append(k)
60   - for j,k in zip(c_gti_start, c_gti_stop):
  61 + for j, k in zip(c_gti_start, c_gti_stop):
61 62 table_c_gti.append(0.5*(j+k))
62 63  
63   - #Initializing the plots
64   - axarr=[]
65   - f = plt.figure(figsize=(12,7))
66   - ax1 = f.add_subplot(121)
67   - ax2 = f.add_subplot(122)
  64 + # Initializing the plots
  65 + axarr = []
  66 + f = plt.figure(figsize=(12,7))
  67 + ax1 = f.add_subplot(121)
  68 + ax2 = f.add_subplot(122)
68 69 axarr.append(ax1)
69 70 axarr.append(ax2)
70 71  
  72 + # Create first plot
71 73 axarr[0].plot(table_c_gti, table_c_sig)
72 74 axarr[0].set_title('Significance evolution')
73 75 axarr[0].set_xlabel('MJD')
74 76 axarr[0].set_ylabel('Significance')
75 77  
76   - axarr[1].hist(table_c_sig, 20, range=[min(table_c_sig), max(table_c_sig)], histtype='stepfilled', facecolor='g', alpha=0.75)
  78 + # Create second plot
  79 + axarr[1].hist(table_c_sig, 20, range=[min(table_c_sig), max(table_c_sig)],
  80 + histtype='stepfilled', facecolor='g', alpha=0.75)
77 81 axarr[1].set_title('Significance distribution')
78 82 axarr[1].set_xlabel('Significance')
79 83 axarr[1].set_ylabel('Counts')
80 84  
81   - #IT SEEMS THAT THE LOG SCALE IS NOT WORKING WELL !!!!!!
82   - #axarr[1].gca().set_yscale('log')
83   -
  85 + # Optionally save result in file
84 86 if plotfile != "":
85 87 f.savefig(plotfile+".png")
86 88  
87 89 plt.show()
88 90  
  91 + # Return
  92 + return
  93 +
  94 +
89 95 # ========================== #
90 96 # Show variability evolution #
91 97 # ========================== #
... ...
src/ctfindvar/ctfindvar.cpp
  View file @0be0a1a
... ... @@ -39,6 +39,7 @@
39 39 #endif
40 40  
41 41 /* __ Method name definitions ____________________________________________ */
  42 +#define G_INX2GTI "ctfindvar::inx2gti(int&)"
42 43 #define G_FILL_CUBE "ctfindvar::fill_cube(GCTAObservation*)"
43 44 #define G_FILL_ALPHA_VECTOR "ctfindvar::fill_alpha_vector(const int&, "\
44 45 "vector<double>&)"
... ... @@ -235,146 +236,16 @@ void ctfindvar::run(void)
235 236 // Write input observation container into logger
236 237 log_observations(NORMAL, m_obs, "Input observation");
237 238  
238   - // Loop over all unbinned CTA observations in the container
239   - #pragma omp parallel for
240   - for (int i = 0; i < m_obs.size(); ++i) {
241   -
242   - // Get pointer to observation
243   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
244   -
245   - // Fill the cube
246   - fill_cube(obs);
247   -
248   - // Dispose events to free memory
249   - obs->dispose_events();
250   -
251   - } // endfor: looped over observations
252   -
253   - // Smooth the maps if requested
254   - if ((*this)["smoothkrnl"].is_valid() &&
255   - (*this)["smoothkrnl"].string() != "NONE") {
256   - m_counts.smooth((*this)["smoothkrnl"].string(),
257   - (*this)["smoothpar"].real());
258   - }
259   -
260   - // Get number of maps
261   - const int nbins = m_counts.nmaps();
262   -
263   - // creating GSkyDir to get the position of the source and each pixels
264   -
265   - // Extract which pixel in the map each source is located in
266   - std::vector<int> srcInxPix;
267   - if (m_inmodel.size() > 0) {
268   - for (int mod = 0; mod < m_inmodel.size(); ++mod) {
269   -
270   - // Extract the model
271   - GModelSky* model = dynamic_cast<GModelSky*>(m_inmodel[mod]);
272   - if (model != NULL) {
273   - // Get the source position
274   - GModelSpatial* model_spatial = model->spatial();
275   - GSkyRegion* model_region = model_spatial->region();
276   - GSkyDir srcSkyDir =
277   - dynamic_cast<GSkyRegionCircle*>(model_region)->centre();
278   -
279   - // Store the map index of this source
280   - int src_map_index = m_counts.dir2inx(srcSkyDir);
281   - srcInxPix.push_back(src_map_index);
282   - } else {
283   - m_inmodel.remove(mod);
284   - mod--;
285   - }
286   - }
287   - }
288   - // ... otherwise only one source position is queried
289   - else {
290   - GSkyDir srcSkyDir;
291   - if ((*this)["coordsys"].string() == "CEL") {
292   - srcSkyDir.radec_deg((*this)["xsrc"].real(),(*this)["ysrc"].real());
293   - } else {
294   - srcSkyDir.lb_deg((*this)["xsrc"].real(),(*this)["ysrc"].real());
295   - }
296   -
297   - // Store the index of this source
298   - int src_map_index = m_counts.dir2inx(srcSkyDir);
299   - srcInxPix.push_back(src_map_index);
300   - }
301   -
302   - //preparing the histogram with significance evolution
303   - //for the source center and the highest sig-pix
304   - m_pixsigsrc = GNdarray(srcInxPix.size(), nbins);
305   - m_pixsigmax = GNdarray(nbins);
306   -
307   - //Prepare the final skymap with the max significance of each pixel
308   - double max_sig = 0;
309   -
310   - log_header1(NORMAL, "Looping Over Pixels");
311   -
312   - //looping over all the pixels in the cube
313   - #pragma omp parallel for
314   - for (int pix_number = 0; pix_number < m_counts.npix(); pix_number++) {
315   -
316   - //Preparing array to store significance evolution of
317   - //each pixel
318   - GNdarray pixSig(nbins);
319   - double total_counts=0;
320   -
321   - //Checking the integrated count of each pixel
322   - for (int k = 0; k < m_counts.nmaps(); k++) {
323   - total_counts += m_counts(pix_number, k);
324   - }
325   -
326   - //Getting the variability significance for the pixel
327   - get_variability_sig(pix_number,nbins, pixSig);
328   -
329   - //Getting the significance evolution for the source
330   - for (int src = 0; src < srcInxPix.size(); src++) {
  239 + // Create GTIs
  240 + init_gtis();
331 241  
332   - // Store the distribution if the source is located at this position
333   - if (srcInxPix[src] == pix_number) {
  242 + // Create counts cube
  243 + create_cube();
334 244  
335   - #pragma omp critical(ctfind_run)
336   - for (int i = 0; i < nbins; i++) {
337   - m_pixsigsrc(src,i) = pixSig(i);
338   - }
339   - }
340   - #ifdef G_DEBUG
341   - std::cout << "checking pixel number of the source of interest" << pix_number << std::endl;
342   - std::cout << "number of counts in pixel of the source of interest" << total_counts << std::endl;
343   - pixSigSrc=pixSig;
344   - std::cout << "checking pixel number " << pix_number << " with total number of counts of: " << total_counts << std::endl;
345   - std::cin.ignore();
346   - std::cout << "checking pixel number of the source of interest: " << pix_number << " - with total number of counts of: " << total_counts << std::endl;
347   - #endif
348   - }
349   - }
350   -
351   - // Getting the evolution for the pix with highest significance
352   - #pragma omp critical(ctfind_run)
353   - {
354   - if (max(pixSig) > max_sig) {
355   - // Update the pixel position with the maximum sigma
356   - max_sig = max(pixSig);
357   - m_max_sig_dir = m_counts.inx2dir(pix_number);
358   - m_pixsigmax = pixSig;
359   - }
360   - }
361   - //storing sig in skymap
362   - m_peaksigmap(pix_number) = max(pixSig);
363   -
364   - //If the significance is greater than sig threshold, the position
365   - //is stored in a model
366   - if (max(pixSig) > m_sig_threshold) {
367   - GSkyDir dir_pix = m_counts.inx2dir(pix_number);
368   - fill_model_sig_pos(dir_pix);
369   - }
370   -
371   - }
372   -
373   - log_header1(NORMAL, "Analysis finished");
374   - log_value(NORMAL, "Maximum sigma", max_sig);
375   - log_value(NORMAL, "Max pixel RA", m_max_sig_dir.ra_deg());
376   - log_value(NORMAL, "Max pixel DEC", m_max_sig_dir.dec_deg());
  245 + // Analyse cube
  246 + analyse_cube();
377 247  
  248 + // Return
378 249 return;
379 250 }
380 251  
... ... @@ -388,37 +259,52 @@ void ctfindvar::run(void)
388 259 ***************************************************************************/
389 260 void ctfindvar::save(void)
390 261 {
  262 + // Write header
391 263 log_header1(TERSE, "Saving results");
392 264  
393   - // Filenames
394   - std::string prefix((*this)["prefix"].string());
395   -
396 265 // Write counts cube
397 266 if ((*this)["outcube"].is_valid()) {
  267 +
  268 + // Get counts cube output filename
398 269 GFilename outcube = (*this)["outcube"].filename();
399   - log_value(TERSE, "Saving counts cube", outcube);
  270 +
  271 + // Save counts cube
400 272 m_counts.save(outcube, (*this)["clobber"].boolean());
  273 +
  274 + // Log saving
  275 + log_value(TERSE, "Saving counts cube", outcube);
401 276 }
402 277  
403 278 // Create a FITS file for storing the output
404   - m_outfile.clear();
  279 + GFits fits;
405 280  
406 281 // Write the most significant values for each pixel
407   - log_string(TERSE, "Appending maximum significances map to output file");
408   - m_peaksigmap.write(m_outfile, "PEAKSIGMAP");
  282 + m_peaksigmap.write(fits, "PEAKSIGMAP");
  283 +
  284 + // Write source histograms
  285 + write_source_histograms(fits);
409 286  
410   - // Write individual source distributions
411   - write_srchist();
  287 + // Get output map filename
  288 + GFilename outmap = (*this)["outmap"].filename();
412 289  
413   - // Write the output file
414   - GFilename outfilename = prefix + "signifmap.fits";
415   - log_value(TERSE, "Saving results to ", outfilename);
416   - m_outfile.saveto(outfilename, (*this)["clobber"].boolean());
  290 + // Save output map
  291 + fits.saveto(outmap, (*this)["clobber"].boolean());
417 292  
418   - //Save the model file with position with significance
  293 + // Log saving
  294 + log_value(TERSE, "Saving output map", outmap);
  295 +
  296 + // Save the model definition file
419 297 if (m_model_above_thr.size() > 2) {
420   - GFilename model_filename = prefix + "model_significant_directions.xml";
421   - m_model_above_thr.save(model_filename);
  298 +
  299 + // Get model output filename
  300 + GFilename outmodel = (*this)["outmodel"].filename();
  301 +
  302 + // Save model
  303 + m_model_above_thr.save(outmodel);
  304 +
  305 + // Log saving
  306 + log_value(TERSE, "Saving output model", outmodel);
  307 +
422 308 }
423 309  
424 310 // Return
... ... @@ -429,7 +315,7 @@ void ctfindvar::save(void)
429 315 /***********************************************************************//**
430 316 * @brief Get the map index associated with a given time
431 317 *
432   - * @param[in] time Time
  318 + * @param[in] time Time
433 319 * @return Index of map in counts cube
434 320 ***************************************************************************/
435 321 int ctfindvar::time2inx(const GTime& time)
... ... @@ -438,29 +324,35 @@ int ctfindvar::time2inx(const GTime&amp; time)
438 324 int map_index = -1;
439 325  
440 326 // Loop over all GTIs
441   - for (int i = 0; i < m_gti.size(); i++) {
442   - // Check if interval contains the time
  327 + for (int i = 0; i < m_gti.size(); ++i) {
  328 +
  329 + // If interval contains the time then break
443 330 if (m_gti[i].contains(time)) {
444 331 map_index = i;
445 332 break;
446 333 }
  334 +
447 335 }
448 336  
  337 + // Return map index
449 338 return map_index;
450 339 }
451 340  
452 341  
453 342 /***********************************************************************//**
454   - * @brief Return reference to significance cube object
  343 + * @brief Return Good Time Interval for index
  344 + *
  345 + * @param[in] index Index
  346 + * @return Good Time Interval
455 347 ***************************************************************************/
456 348 GGti ctfindvar::inx2gti(const int& index)
457 349 {
  350 + // Throw an exception if index is invalid
458 351 if ((index < 0) || (index >= m_gti.size())) {
459   - // Index is invalid, so throw an error
460   - throw GException::invalid_value("ctfindvar::inx2gti(const int&)",
461   - "'index' parameter out of range");
  352 + throw GException::out_of_range(G_INX2GTI, "Time index", index, m_gti.size());
462 353 }
463 354  
  355 + // Return Good Time Interval
464 356 return m_gti[index];
465 357 }
466 358  
... ... @@ -483,7 +375,6 @@ void ctfindvar::init_members(void)
483 375 m_max_sig_dir.clear();
484 376 m_minoff = 0.0;
485 377 m_sig_threshold = 0.0;
486   - m_outfile.clear();
487 378 m_peaksigmap.clear();
488 379 m_pixsigsrc.clear();
489 380 m_tstart.clear();
... ... @@ -511,7 +402,6 @@ void ctfindvar::copy_members(const ctfindvar&amp; app)
511 402 m_max_sig_dir = app.m_max_sig_dir;
512 403 m_minoff = app.m_minoff;
513 404 m_sig_threshold = app.m_sig_threshold;
514   - m_outfile = app.m_outfile;
515 405 m_peaksigmap = app.m_peaksigmap;
516 406 m_pixsigsrc = app.m_pixsigsrc;
517 407 m_tstart = app.m_tstart;
... ... @@ -537,39 +427,35 @@ void ctfindvar::free_members(void)
537 427  
538 428 /***********************************************************************//**
539 429 * @brief Get application parameters
540   - *
541   - * @todo Implement method
542 430 ***************************************************************************/
543 431 void ctfindvar::get_parameters(void)
544 432 {
545   - #ifdef _OPENMP
546   - int nthreads = (*this)["nthreads"].integer();
547   - if (nthreads > 0) {
548   - omp_set_num_threads(nthreads);
549   - }
550   - #endif
551   -
552 433 // Load the observations
553 434 setup_observations(m_obs, true, true, false);
554 435  
555   - // Get the rest of the parameters
556   - if (read_ahead()) {
557   - (*this)["prefix"].string();
558   -
559   - if ((*this)["inmodel"].is_valid()) {
560   - m_inmodel = GModels((*this)["inmodel"].filename());
561   - } else {
562   - // Get the desired source position from user
563   - (*this)["xsrc"].real();
564   - (*this)["ysrc"].real();
565   - }
  436 + // Either load model or query source position parameters
  437 + if ((*this)["inmodel"].is_valid()) {
  438 + m_inmodel = GModels((*this)["inmodel"].filename());
  439 + }
  440 + else {
  441 + (*this)["coordsys"].string();
  442 + (*this)["xsrc"].real();
  443 + (*this)["ysrc"].real();
566 444 }
567 445  
568   - // Create GTIs and counts cube
569   - init_gtis();
570   - init_cube();
  446 + // Query time limit parameters
  447 + (*this)["tinterval"].real();
  448 + if ((*this)["tmin"].is_valid()) {
  449 + (*this)["tmin"].time();
  450 + }
  451 + if ((*this)["tmax"].is_valid()) {
  452 + (*this)["tmax"].time();
  453 + }
  454 +
  455 + // Create map to query spatial parameters
  456 + create_map(m_obs);
571 457  
572   - // Get the energy limits
  458 + // Read energy limits
573 459 m_emin = GEnergy((*this)["emin"].real(), "TeV");
574 460 m_emax = GEnergy((*this)["emax"].real(), "TeV");
575 461  
... ... @@ -577,7 +463,28 @@ void ctfindvar::get_parameters(void)
577 463 m_minoff = (*this)["minoff"].real();
578 464  
579 465 // Get minimum significance to set a source as variable
580   - m_sig_threshold = (*this)["sig_thr"].real();
  466 + m_sig_threshold = (*this)["threshold"].real();
  467 +
  468 + // Query smoothing parameters
  469 + if ((*this)["smooth_kernel"].is_valid()) {
  470 + (*this)["smooth_kernel"].string();
  471 + (*this)["smooth_rad"].real();
  472 + }
  473 +
  474 + // If needed later, query output filenames now
  475 + if (read_ahead()) {
  476 + (*this)["outcube"].query();
  477 + (*this)["outmap"].query();
  478 + (*this)["outmodel"].query();
  479 + }
  480 +
  481 + // Set number of OpenMP threads
  482 + #ifdef _OPENMP
  483 + int nthreads = (*this)["nthreads"].integer();
  484 + if (nthreads > 0) {
  485 + omp_set_num_threads(nthreads);
  486 + }
  487 + #endif
581 488  
582 489 // Write parameters into logger
583 490 log_parameters(TERSE);
... ... @@ -588,117 +495,139 @@ void ctfindvar::get_parameters(void)
588 495  
589 496  
590 497 /***********************************************************************//**
591   - * @brief Initialize counts cube
  498 + * @brief Initialize Good Time Intervals
592 499 *
593   - * Sets up the required event cube and time interval storage objects
594   - * associated with this object
595   - ***************************************************************************/
596   -void ctfindvar::init_cube(void)
597   -{
598   - // Create the basic skymap
599   - m_counts = create_map(m_obs);
600   -
601   - // Resize to the appropriate number of time intervals
602   - m_counts.nmaps(m_gti.size());
603   -
604   - // Create the peaksigmap
605   - m_peaksigmap = m_counts.extract(1);
606   -}
607   -
608   -
609   -/***********************************************************************//**
610   - * @brief Initialize counts cube
  500 + * @exception GException::invalid_value
  501 + * Invalid start or stop time or insufficient time bins
611 502 *
612   - * Sets up the required event cube and time interval storage objects
613   - * associated with this object
  503 + * Initialise Good Time Intervals for the variability search.
614 504 ***************************************************************************/
615 505 void ctfindvar::init_gtis(void)
616 506 {
617   - double tinterval = (*this)["tinterval"].real();
618   - m_tstart = GTime("JD 999999999"); // large unreasonable Julian date
619   - m_tstop = GTime("JD 0"); // small unreasonable Julian date
  507 + // Log header
  508 + log_header1(NORMAL, "Initialise Good Time Intervals");
620 509  
621   - // Set the start time
622   - if ((*this)["tmin"].is_valid()) {
623   - log_string(NORMAL, "Setting start time from command line");
624   - m_tstart = (*this)["tmin"].time();
625   - } else {
626   - log_string(NORMAL, "Setting start time from observations");
627   - for (int o=0; o<m_obs.size(); o++) {
628   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[o]);
629   -
630   - // Check that the start time of the run is less than tstart
631   - if (obs->gti().tstart() < m_tstart) {
632   - m_tstart = obs->gti().tstart();
633   - }
634   - }
635   - }
  510 + // Clear time intervals
  511 + m_gti.clear();
636 512  
637   - // Set the stop time
638   - if ((*this)["tmax"].is_valid()) {
639   - log_string(NORMAL, "Setting stop time from command line");
640   - m_tstop = (*this)["tmax"].time();
641   - } else {
  513 + // Set start and stop times
  514 + m_tstart = get_tstart();
  515 + m_tstop = get_tstop();
642 516  
643   - log_string(NORMAL, "Setting stop time from observations");
644   - for (int o=0; o<m_obs.size(); o++) {
645   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[o]);
646   -
647   - // Check that the stop time of the run is larger than tstop
648   - if (obs->gti().tstop() > m_tstop) {
649   - m_tstop = obs->gti().tstop();
650   - }
651   - }
652   - }
  517 + // Query time interval
  518 + double tinterval = (*this)["tinterval"].real();
653 519  
654   - // Fill the number of good time intervals
  520 + // Compute time bins
655 521 double tstart_sec = m_tstart.secs();
656 522 double tstop_sec = m_tstop.secs();
657 523 int bins = (tstop_sec - tstart_sec) / tinterval + 0.5;
658 524  
659 525 // Log information
660   - log_value(NORMAL, "treference (mjd)", m_tstart.reference().mjdref());
661   - log_value(NORMAL, "tstart (sec)", tstart_sec);
662   - log_value(NORMAL, "tstop (sec)", tstop_sec);
663   - log_value(NORMAL, "Total time", tstop_sec-tstart_sec);
664   - log_value(NORMAL, "Time interval", tinterval);
665   - log_value(NORMAL, "nbins", bins);
666   -
667   - // Throw if tstart == tstop
  526 + log_value(NORMAL, "Reference time (mjd)", m_tstart.reference().mjdref());
  527 + log_value(NORMAL, "Start time (sec)", tstart_sec);
  528 + log_value(NORMAL, "Stop time (sec)", tstop_sec);
  529 + log_value(NORMAL, "Total time (sec)", tstop_sec-tstart_sec);
  530 + log_value(NORMAL, "Time interval (sec)", tinterval);
  531 + log_value(NORMAL, "Number of time bins", bins);
  532 +
  533 + // Throw exception if tstart == tstop
668 534 if (m_tstart == m_tstop) {
669 535 throw GException::invalid_value(G_INIT_GTIS,
670 536 "Start time is equal to stop time");
671   - } else if (bins < 2) {
  537 + }
  538 + else if (bins < 2) {
672 539 std::string msg = "Method requires at least two time bins (" +
673 540 gammalib::str(bins) + " bins found).";
674 541 throw GException::invalid_value(G_INIT_GTIS, msg);
675 542 }
676 543  
677   - for (int i=0; i<bins; i++) {
  544 + // Set up Good Time Intervals by appending all intervals that overlap
  545 + // with at least one of the observations
  546 + for (int i = 0; i < bins; ++i) {
678 547  
679   - // Get the next time interval
680   - GTime next_tstart(m_tstart + i*tinterval);
681   - GTime next_tstop(m_tstart + (i+1.0)*tinterval);
682   - GGti next_gti(next_tstart, next_tstop);
  548 + // Get time interval
  549 + GTime tstart(m_tstart + i*tinterval);
  550 + GTime tstop(m_tstart + (i+1.0)*tinterval);
  551 + GGti gti(tstart, tstop);
683 552  
684 553 // Make sure there's an observation that overlaps with this GTI
685   - for (int o=0; o<m_obs.size(); o++) {
686   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[o]);
687   -
688   - if (gti_overlap(obs->gti(), next_gti) > 0.0) {
689   - m_gti.push_back(next_gti);
690   - break;
691   - }
692   - }
  554 + for (int k = 0; k < m_obs.size(); ++k) {
  555 +
  556 + // Get CTA observation
  557 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[k]);
  558 +
  559 + // Continue only if observation is valid
  560 + if (obs != NULL) {
  561 +
  562 + // Add GTI is the time intervals overlap
  563 + if (gti_overlap(obs->gti(), gti) > 0.0) {
  564 + m_gti.push_back(gti);
  565 + break;
  566 + }
  567 +
  568 + } // endif: observation was valid
  569 +
  570 + } // endfor: looped over observations
  571 +
  572 + } // endfor: looped over Good Time Intervals
  573 +
  574 + // Return
  575 + return;
  576 +}
  577 +
  578 +
  579 +/***********************************************************************//**
  580 + * @brief Create counts cube
  581 + *
  582 + * Creates a counts cube comprising a sky map for each Good Time Interval.
  583 + ***************************************************************************/
  584 +void ctfindvar::create_cube(void)
  585 +{
  586 + // Log header
  587 + log_header1(NORMAL, "Create counts cube");
  588 +
  589 + // Create the basic skymap
  590 + m_counts = create_map(m_obs);
  591 +
  592 + // Resize to the appropriate number of time intervals
  593 + m_counts.nmaps(m_gti.size());
  594 +
  595 + // Create the peaksigmap
  596 + m_peaksigmap = m_counts.extract(1);
  597 +
  598 + // Loop over all unbinned CTA observations in the container
  599 + #pragma omp parallel for
  600 + for (int i = 0; i < m_obs.size(); ++i) {
  601 +
  602 + // Get pointer to observation
  603 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  604 +
  605 + // Fill the cube
  606 + fill_cube(obs);
  607 +
  608 + // Dispose events to free memory
  609 + obs->dispose_events();
  610 +
  611 + } // endfor: looped over observations
  612 +
  613 + // Smooth the maps if requested
  614 + if ((*this)["smooth_kernel"].is_valid() &&
  615 + (*this)["smooth_kernel"].string() != "NONE") {
  616 + m_counts.smooth((*this)["smooth_kernel"].string(),
  617 + (*this)["smooth_rad"].real());
693 618 }
  619 +
  620 + // Return
  621 + return;
694 622 }
695 623  
696 624  
697 625 /***********************************************************************//**
698   - * @brief Fill the cube from the events in the map
  626 + * @brief Fill the cube from the events in an observation
699 627 *
700   - * Sets up the required event cube and time interval storage objects
701   - * associated with this object
  628 + * @param[in] obs Pointer to CTA observation.
  629 + *
  630 + * Fills the cube from the events in an observation.
702 631 ***************************************************************************/
703 632 void ctfindvar::fill_cube(GCTAObservation* obs)
704 633 {
... ... @@ -733,10 +662,6 @@ void ctfindvar::fill_cube(GCTAObservation* obs)
733 662 int num_outside_time = 0;
734 663 int num_in_map = 0;
735 664  
736   - // Specify whether we are cutting on energy or not
737   - bool cut_emin = (m_emin.TeV() > 0.0) ? true : false;
738   - bool cut_emax = (m_emax.TeV() > 0.0) ? true : false;
739   -
740 665 // Fill counts sky map
741 666 for (int i = 0; i < events->size(); ++i) {
742 667  
... ... @@ -744,16 +669,14 @@ void ctfindvar::fill_cube(GCTAObservation* obs)
744 669 const GCTAEventAtom* event = (*events)[i];
745 670  
746 671 // Check that the energy is valid
747   - if ((cut_emin && event->energy() < m_emin) ||
748   - (cut_emax && event->energy() > m_emax)) {
  672 + if ((event->energy() < m_emin) || (event->energy() > m_emax)) {
749 673 num_outside_ecut++;
750 674 continue;
751 675 }
752 676  
753   - // Determine event time
754   - GTime evnt_time = event->time();
755   - int time_indx = time2inx(evnt_time);
756   -
  677 + // Determine event time index
  678 + int time_indx = time2inx(event->time());
  679 +
757 680 // Check if event is in the GTIs
758 681 if ((time_indx > -1) && (time_indx < m_counts.nmaps())) {
759 682  
... ... @@ -791,10 +714,11 @@ void ctfindvar::fill_cube(GCTAObservation* obs)
791 714 // Increment number of maps
792 715 num_in_map++;
793 716  
794   - } else {
795   - // Event falls outside the specified time range
  717 + }
  718 +
  719 + // ... otherwise, event falls outside the specified time range
  720 + else {
796 721 num_outside_time++;
797   - continue;
798 722 }
799 723  
800 724 } // endfor: looped over all events
... ... @@ -818,242 +742,557 @@ void ctfindvar::fill_cube(GCTAObservation* obs)
818 742  
819 743  
820 744 /***********************************************************************//**
821   - * @brief Save the current sky direction to model if significance > thr
822   - *
823   - * @param[in] dir_pix pixel sky direction
824   - * @param[out] m_model_sig_above_thr model position with sig>thr
  745 + * @brief Analyse all pixels of counts cube
825 746 ***************************************************************************/
826   -void ctfindvar::fill_model_sig_pos(const GSkyDir& dir_pix)
  747 +void ctfindvar::analyse_cube(void)
827 748 {
828   - //Set spatial components
829   - GModelSpatialPointSource spatial(dir_pix);
  749 + // Log header
  750 + log_header1(NORMAL, "Analyse sky map pixels");
830 751  
831   - //Set fake spectral component (powerlaw with 1% of Crab)
832   - GModelSpectralPlaw spectral(5.7e-18, -2.48, GEnergy(0.3, "TeV"));
  752 + // Get number of GTIs
  753 + const int nbins = m_gti.size();
833 754  
834   - //Creat a sky model
835   - GModelSky model(spatial, spectral);
836   -
837   - //Extracting the coordinate to name the model
838   - double ra = dir_pix.ra_deg();
839   - double dec = dir_pix.dec_deg();
840   - std::string name = "Sky_direction_" + gammalib::str(ra, 4) +
841   - "_" + gammalib::str(dec, 4);
842   - model.name(name);
843   -
844   - //Add the model for this pixel to the output mode
845   - m_model_above_thr.append(model);
  755 + // Get relevant sky map pixels
  756 + std::vector<int> srcInxPix = get_pixels();
  757 +
  758 + // Preparing the histogram with significance evolution for the source
  759 + // center and the highest sig-pix
  760 + m_pixsigsrc = GNdarray(srcInxPix.size(), nbins);
  761 + m_pixsigmax = GNdarray(nbins);
  762 +
  763 + // Initialise maximum significance
  764 + double max_sig = 0;
  765 +
  766 + // looping over all the pixels in the cube
  767 + #pragma omp parallel for
  768 + for (int ipix = 0; ipix < m_counts.npix(); ++ipix) {
  769 +
  770 + // Compute total counts for pixel
  771 + double total_counts = 0.0;
  772 + for (int k = 0; k < nbins; ++k) {
  773 + total_counts += m_counts(ipix, k);
  774 + }
  775 +
  776 + // Getting the variability significance for the current pixel
  777 + GNdarray pixSig = get_variability_sig(ipix);
  778 +
  779 + // Getting the significance evolution for the source
  780 + for (int src = 0; src < srcInxPix.size(); ++src) {
  781 +
  782 + // Store the distribution if the source is located at this position
  783 + if (srcInxPix[src] == ipix) {
  784 +
  785 + #pragma omp critical(ctfindvar_analyse_cube)
  786 + for (int i = 0; i < nbins; i++) {
  787 + m_pixsigsrc(src,i) = pixSig(i);
  788 + }
  789 + }
  790 +
  791 + } // endfor: looped over sources
  792 +
  793 + // Make sure that this code is not executed in parallel
  794 + #pragma omp critical(ctfindvar_analyse_cube)
  795 + {
  796 + // Determine the pixel with the highest significance
  797 + if (max(pixSig) > max_sig) {
  798 + max_sig = max(pixSig);
  799 + m_max_sig_dir = m_counts.inx2dir(ipix);
  800 + m_pixsigmax = pixSig;
  801 + }
  802 +
  803 + // Storing sig in skymap
  804 + m_peaksigmap(ipix) = max(pixSig);
  805 +
  806 + // If the significance is greater than the significance threshold,
  807 + // the position is stored in a model
  808 + if (max(pixSig) > m_sig_threshold) {
  809 + GSkyDir dir_pix = m_counts.inx2dir(ipix);
  810 + m_model_above_thr.append(sky_model(dir_pix));
  811 + }
  812 + }
846 813  
  814 + } // endfor: looped over pixels
  815 +
  816 + // Log results
  817 + log_value(NORMAL, "Maximum significance", max_sig);
  818 + log_value(NORMAL, "Right Ascension of maximum", m_max_sig_dir.ra_deg());
  819 + log_value(NORMAL, "Declination of maximum", m_max_sig_dir.dec_deg());
  820 +
  821 + // Return
  822 + return;
847 823 }
  824 +
  825 +
848 826 /***********************************************************************//**
849   - * @brief Get the map index associated with a given time
  827 + * @brief Return pixel vector
  828 + *
  829 + * @return Vector of sky map pixels
850 830 *
851   - * @param[in] pix_number Spatial pixel
852   - * @param[out] alpha_vector Vector containing effective exposure
  831 + * Returns a vector of sky map pixels that corresponds to the source model
  832 + * positions.
853 833 ***************************************************************************/
854   -void ctfindvar::fill_alpha_vector(const int& pix_number,
855   - std::vector<double>& alpha_vector)
  834 +std::vector<int> ctfindvar::get_pixels(void)
856 835 {
857   - // Pixel position
858   - GSkyDir pix_dir = m_counts.inx2dir(pix_number);
  836 + // Initialise pixel array
  837 + std::vector<int> srcInxPix;
859 838  
860   - // Loop over all observations
861   - for (int i = 0; i < m_obs.size(); i++) {
862   - // Skip if observation does not overlap with this pixel position
863   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
864   - GCTARoi roi = obs->roi();
865   - GSkyDir roi_centre = roi.centre().dir();
  839 + // If there are input models then extract
  840 + if (m_inmodel.size() > 0) {
866 841  
867   - if (roi_centre.dist_deg(pix_dir) > roi.radius()) {
868   - continue;
869   - }
870   -
871   - // Convert sky direction to instrument direction
872   - GCTAInstDir instdir = obs->pointing().instdir(pix_dir);
  842 + // Loop over models
  843 + for (int i = 0; i < m_inmodel.size(); ++i) {
873 844  
874   - // Extract the good time intervals of the observation
875   - const GGti& obs_gti(obs->gti());
  845 + // Extract the model
  846 + GModelSky* model = dynamic_cast<GModelSky*>(m_inmodel[i]);
876 847  
877   - // Loop over all time bins
878   - for (int j = 0; j < m_gti.size(); j++) {
  848 + // If model is a sky model then extract the model centre. Note
  849 + // that this only works for a model that results a GSkyRegionCircle
  850 + // object as region.
  851 + if (model != NULL) {
879 852  
880   - // Make sure observation overlaps with this time interval
881   - double exposure = gti_overlap(m_gti[j], obs_gti);
882   - if (exposure > 0.0) {
  853 + // Get the source position
  854 + GModelSpatial* spatial = model->spatial();
  855 + GSkyRegion* region = spatial->region();
  856 + GSkyRegionCircle* circle = dynamic_cast<GSkyRegionCircle*>(region);
  857 +
  858 + // If region circle is valid then
  859 + if (circle != NULL) {
  860 + int ipix = m_counts.dir2inx(circle->centre());
  861 + srcInxPix.push_back(ipix);
  862 + continue;
  863 + }
883 864  
884   - // Get IRF response
885   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
886   - (obs->response());
  865 + } // endif: model was a sky model
887 866  
888   - // Throw an exception if observation has no instrument response function
889   - if (rsp == NULL) {
890   - std::string msg =
891   - "No response information available for "+
892   - get_obs_header(obs)+" to compute IRF background. "
893   - "Please specify response information or use "
894   - "another background subtraction method.";
895   - throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg);
896   - }
  867 + // If we are still alive then remove the model since it is neither
  868 + // a sky model, nor does it have a region circle
  869 + m_inmodel.remove(i);
  870 + i--;
897 871  
898   - // Get IRF background template
899   - const GCTABackground* bkg = rsp->background();
  872 + } // endfor: looped over all models
900 873  
901   - // Throw an exception if observation has no IRF background template
902   - if (bkg == NULL) {
903   - std::string msg =
904   - "No IRF background template found in instrument "
905   - "response function for "+
906   - get_obs_header(obs)+". Please specify an instrument "
907   - "response function containing a background template.";
908   - throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg);
909   - }
  874 + } // endif: there were input models
  875 +
  876 + // ... otherwise query the source position
  877 + else {
910 878  
911   - // Multiply background rate with livetime and solid angle
912   - #pragma omp critical(ctfindvar_alpha_vector)
913   - alpha_vector[j] += exposure * bkg->rate_ebin(instdir, m_emin, m_emax);
914   - }
  879 + // Query source direction
  880 + GSkyDir dir;
  881 + if ((*this)["coordsys"].string() == "CEL") {
  882 + dir.radec_deg((*this)["xsrc"].real(),(*this)["ysrc"].real());
  883 + }
  884 + else {
  885 + dir.lb_deg((*this)["xsrc"].real(),(*this)["ysrc"].real());
915 886 }
916   - }
917 887  
918   - return;
919   -}
  888 + // Store the index of this source
  889 + int ipix = m_counts.dir2inx(dir);
  890 + srcInxPix.push_back(ipix);
  891 +
  892 + } // endelse: queried source position
  893 +
  894 + // Return pixel array
  895 + return srcInxPix;
  896 +}
  897 +
  898 +
  899 +/***********************************************************************//**
  900 + * @brief Return sky model for a given sky direction
  901 + *
  902 + * @param[in] dir Sky direction
  903 + * @return Sky model for a given sky direction
  904 + *
  905 + * Returns a point source sky model at the position of the current pixel
  906 + * with a power law spectral model corresponding to 1% of the Crab flux.
  907 + *
  908 + * The model name is built from the celestial coordinates of the source.
  909 + ***************************************************************************/
  910 +GModelSky ctfindvar::sky_model(const GSkyDir& dir) const
  911 +{
  912 + // Set point source spatial component from sky direction
  913 + GModelSpatialPointSource spatial(dir);
920 914  
  915 + // Set fake spectral component (powerlaw with 1% of Crab)
  916 + GModelSpectralPlaw spectral(5.7e-18, -2.48, GEnergy(0.3, "TeV"));
921 917  
  918 + // Create a sky model
  919 + GModelSky model(spatial, spectral);
  920 +
  921 + // Build model name from celestial coordinates
  922 + double ra = dir.ra_deg();
  923 + double dec = dir.dec_deg();
  924 + std::string name = gammalib::str(ra, 4);
  925 + if (dec < 0.0) {
  926 + name += "-";
  927 + }
  928 + else {
  929 + name += "+";
  930 + }
  931 + name += gammalib::str(std::abs(dec), 4);
  932 +
  933 + // Set model name
  934 + model.name(name);
  935 +
  936 + // Return model
  937 + return model;
  938 +}
922 939  
923 940  
924 941 /***********************************************************************//**
925   - * @brief Obtain the significance of variability in a given bin
  942 + * @brief Get the significance of variability for a given skymap pixel
926 943 *
927   - * @param[in] pix_number Pixel in the skymap
928   - * @param[in] nbins Number of bins in the skymap
929   - * @param[out] sig_histogram Histogram for storing returned significances
  944 + * @param[in] ipix Sky map pixel index.
  945 + * @return Histogram of significances.
930 946 *
931 947 * Significance is computed according to Li & Ma equation 17:
932 948 * see: https://doi.org/10.1086/161295
933 949 * Some modification is made in order to handle the case where Non = 0,
934 950 * but Noff != 0.
935 951 ***************************************************************************/
936   -void ctfindvar::get_variability_sig(const int& pix_number,
937   - const int& nbins,
938   - GNdarray& sig_histogram)
  952 +GNdarray ctfindvar::get_variability_sig(const int& ipix)
939 953 {
940   - std::vector<bool> accepted_bin_bckg_vector(nbins, true);
941   - bool background_validated = false;
942   - double non, noff;
943   - double alpha, sig;
  954 + // Initialise result
  955 + GNdarray sig_histogram(m_gti.size());
  956 +
  957 + // Initialise vectors
  958 + std::vector<bool> accepted_bin_bckg_vector(m_gti.size(), true);
944 959  
945   - // Fill the alpha values for each vector
946   - std::vector<double> alpha_vector(nbins, 0.0);
947   - fill_alpha_vector(pix_number, alpha_vector);
  960 + // Get alpha values for specified pixel
  961 + std::vector<double> alphas = get_alphas(ipix);
948 962  
949   - // Loop over pixels until all background pixls are removed
  963 + // Loop over pixels until all background pixels are removed
  964 + bool background_validated = false;
950 965 while (background_validated == false) {
  966 +
  967 + // Signal that background was validated
951 968 background_validated = true;
952 969  
953 970 // Loop over all the GTIs of the pixel
954   - for (int i = 0; i < nbins; i++) {
  971 + for (int i = 0; i < m_gti.size(); ++i) {
955 972  
956   - // The GTI is discared from bckg calculation and not checked again.
  973 + // The GTI is discared from background calculation and not
  974 + // checked again
957 975 if (!accepted_bin_bckg_vector[i]) {
958 976 continue;
959   - }
960   - // Check if bin fails minoff check or no observations overlap it
961   - else if (m_counts(pix_number, i) < m_minoff || alpha_vector[i] == 0) {
  977 + }
  978 +
  979 + // ... otherwise check if bin fails minoff check or no observations
  980 + // overlap it
  981 + else if (m_counts(ipix, i) < m_minoff || alphas[i] == 0) {
962 982 accepted_bin_bckg_vector[i] = false;
963 983 continue;
964 984 }
965 985  
966   - noff = 0.0;
967   - alpha = 0.0;
968   - // Gor one GTI selected (i), looping over all the others (j).
969   - for (int j = 0; j < nbins; j++) {
  986 + // ... otherwise
  987 + double noff = 0.0;
  988 + double alpha = 0.0;
  989 +
  990 + // Got one GTI selected (i), looping over all the others (j).
  991 + for (int j = 0; j < m_gti.size(); ++j) {
970 992 if (j != i && accepted_bin_bckg_vector[j] == 1) {
971   - noff += m_counts(pix_number, j);
972   - alpha += alpha_vector[j];
  993 + noff += m_counts(ipix, j);
  994 + alpha += alphas[j];
973 995 }
974 996 }
975 997  
976 998 // The background is averaged on the number of bins -1
977   - non = m_counts(pix_number, i);
978   - alpha = alpha_vector[i]/alpha;
  999 + double non = m_counts(ipix, i);
  1000 + alpha = alphas[i]/alpha;
979 1001  
980 1002 // Compute sensitivity in Gaussian sigma using Li & Ma
981 1003 double alpha1 = alpha + 1.0;
982 1004 double ntotal = non+noff;
983 1005 double arg1 = non/ntotal;
984 1006 double arg2 = noff/ntotal;
985   - if (noff == 0.0) {
986   - sig = 0.0;
987   - } else if (non > 0.0) {
988   - double term1 = non * std::log((alpha1/alpha)*arg1);
989   - double term2 = noff * std::log(alpha1*arg2);
990   - sig = std::sqrt(2.0 * (term1 + term2));
991   - } else {
992   - sig = std::sqrt(2.0 * noff * std::log(alpha1));
993   - }
994   -
995   - // Specify the sign of the significance
996   - sig *= (non < alpha*noff) ? -1.0 : 1.0;
997   -
998   - if (alpha == 0.0) {
999   - sig = 0.0;
  1007 + double sig = 0.0;
  1008 + if (alpha != 0.0) {
  1009 + if (noff == 0.0) {
  1010 + sig = 0.0;
  1011 + }
  1012 + else if (non > 0.0) {
  1013 + double term1 = non * std::log((alpha1/alpha)*arg1);
  1014 + double term2 = noff * std::log(alpha1*arg2);
  1015 + sig = std::sqrt(2.0 * (term1 + term2));
  1016 + }
  1017 + else {
  1018 + sig = std::sqrt(2.0 * noff * std::log(alpha1));
  1019 + }
  1020 + sig *= (non < alpha*noff) ? -1.0 : 1.0;
1000 1021 }
1001 1022  
1002 1023 // Update the significance
1003 1024 sig_histogram(i) = sig;
1004 1025  
1005   - // If bin is significant, it's removed from bckg and we loop again.
  1026 + // If bin is significant, it's removed from the background and
  1027 + // we loop again
1006 1028 if (sig > m_sig_threshold) {
1007 1029 accepted_bin_bckg_vector[i] = false;
1008 1030 background_validated = false;
1009 1031 }
1010   - }
1011 1032  
1012   - }
  1033 + } // endfor: looped over all GTIs
  1034 +
  1035 + } // endwhile: iterate until converged
  1036 +
  1037 + // Return significance histogram
  1038 + return sig_histogram;
1013 1039 }
1014 1040  
1015 1041  
1016 1042 /***********************************************************************//**
  1043 + * @brief Get alpha vector
  1044 + *
  1045 + * @param[in] ipix Sky map pixel index
  1046 + * @return Vector containing effective exposure
  1047 + ***************************************************************************/
  1048 +std::vector<double> ctfindvar::get_alphas(const int& ipix) const
  1049 +{
  1050 + // Initialise alpha vector
  1051 + std::vector<double> alphas(m_gti.size(), 0.0);
  1052 +
  1053 + // Convert pixel index into sky direction
  1054 + GSkyDir pix_dir = m_counts.inx2dir(ipix);
  1055 +
  1056 + // Loop over all observations
  1057 + for (int i = 0; i < m_obs.size(); ++i) {
  1058 +
  1059 + // Get CTA observation
  1060 + const GCTAObservation* obs = dynamic_cast<const GCTAObservation*>(m_obs[i]);
  1061 +
  1062 + // Fall through if observation is not valid
  1063 + if (obs == NULL) {
  1064 + continue;
  1065 + }
  1066 +
  1067 + // Get RoI. Fall through if observation does not overlap with this
  1068 + // pixel position
  1069 + GCTARoi roi = obs->roi();
  1070 + GSkyDir roi_centre = roi.centre().dir();
  1071 + if (roi_centre.dist_deg(pix_dir) > roi.radius()) {
  1072 + continue;
  1073 + }
  1074 +
  1075 + // Convert sky direction to instrument direction
  1076 + GCTAInstDir instdir = obs->pointing().instdir(pix_dir);
  1077 +
  1078 + // Extract the good time intervals of the observation
  1079 + const GGti& obs_gti(obs->gti());
  1080 +
  1081 + // Loop over all time bins
  1082 + for (int j = 0; j < m_gti.size(); ++j) {
  1083 +
  1084 + // Make sure observation overlaps with this time interval
  1085 + double exposure = gti_overlap(m_gti[j], obs_gti);
  1086 + if (exposure > 0.0) {
  1087 +
  1088 + // Get IRF response
  1089 + const GCTAResponseIrf* rsp =
  1090 + dynamic_cast<const GCTAResponseIrf*>(obs->response());
  1091 +
  1092 + // Throw an exception if observation has no instrument response function
  1093 + if (rsp == NULL) {
  1094 + std::string msg =
  1095 + "No response information available for "+
  1096 + get_obs_header(obs)+" to compute IRF background. "
  1097 + "Please specify response information or use "
  1098 + "another background subtraction method.";
  1099 + throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg);
  1100 + }
  1101 +
  1102 + // Get IRF background template
  1103 + const GCTABackground* bkg = rsp->background();
  1104 +
  1105 + // Throw an exception if observation has no IRF background template
  1106 + if (bkg == NULL) {
  1107 + std::string msg =
  1108 + "No IRF background template found in instrument "
  1109 + "response function for "+
  1110 + get_obs_header(obs)+". Please specify an instrument "
  1111 + "response function containing a background template.";
  1112 + throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg);
  1113 + }
  1114 +
  1115 + // Add up background rate
  1116 + alphas[j] += exposure * bkg->rate_ebin(instdir, m_emin, m_emax);
  1117 +
  1118 + } // endif: exposure was positive
  1119 +
  1120 + } // endfor: looped over all time bins
  1121 +
  1122 + } // endfor: looped over all observations
  1123 +
  1124 + // Return alphas
  1125 + return alphas;
  1126 +}
  1127 +
  1128 +
  1129 +/***********************************************************************//**
1017 1130 * @brief Returns number of seconds that two GTIs overlap
1018 1131 *
1019 1132 * @param[in] gti1 First GTI
1020 1133 * @param[in] gti2 Second GTI
1021 1134 * @return Number of seconds that @p gti1 and @p gti2 overlap
1022 1135 ***************************************************************************/
1023   -double ctfindvar::gti_overlap(const GGti& gti1, const GGti& gti2)
  1136 +double ctfindvar::gti_overlap(const GGti& gti1, const GGti& gti2) const
1024 1137 {
  1138 + // Initialise no overlap
1025 1139 double overlap = 0.0;
  1140 +
  1141 + // Set first and second GTI so that gti_1st starts before gti_2nd
1026 1142 const GGti& gti_1st = (gti1.tstart() <= gti2.tstart()) ? gti1 : gti2;
1027   - const GGti& gti_2nd = (gti1.tstart() > gti2.tstart()) ? gti1 : gti2;
  1143 + const GGti& gti_2nd = (gti1.tstart() > gti2.tstart()) ? gti1 : gti2;
1028 1144  
1029   - // gti1 starts earlier
  1145 + // Compute overlap
1030 1146 if (gti_1st.tstop() > gti_2nd.tstart()) {
1031 1147 if (gti_1st.tstop() <= gti_2nd.tstop()) {
1032 1148 overlap = gti_1st.tstop() - gti_2nd.tstart();
1033   - } else if (gti_1st.tstop() > gti_2nd.tstop()) {
  1149 + }
  1150 + else if (gti_1st.tstop() > gti_2nd.tstop()) {
1034 1151 overlap = gti_2nd.tstop() - gti_2nd.tstart();
1035 1152 }
1036 1153 }
1037 1154  
  1155 + // Return overlap
1038 1156 return overlap;
1039 1157 }
1040 1158  
1041 1159  
1042 1160 /***********************************************************************//**
1043   - * @brief Write individual histograms in CSV format
  1161 + * @brief Get start time
  1162 + *
  1163 + * @return Start time
  1164 + ***************************************************************************/
  1165 +GTime ctfindvar::get_tstart(void)
  1166 +{
  1167 + // Initialise start time with large unreasonable Julian date
  1168 + GTime tstart("JD 999999999");
  1169 +
  1170 + // If tmin parameter is valid then set the start time
  1171 + if ((*this)["tmin"].is_valid()) {
  1172 +
  1173 + // Query tmin parameter
  1174 + tstart = (*this)["tmin"].time();
  1175 +
  1176 + // Log action
  1177 + log_string(NORMAL, "Setting start time \"tmin\" parameter");
  1178 +
  1179 + } // endif: tmin parameter was valid
  1180 +
  1181 + // ... otherwise set start time from observations
  1182 + else {
  1183 +
  1184 + // Loop over all observations
  1185 + for (int k = 0; k < m_obs.size(); ++k) {
  1186 +
  1187 + // Get CTA observation
  1188 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[k]);
  1189 +
  1190 + // Continue only if observation is valid
  1191 + if (obs != NULL) {
  1192 +
  1193 + // Update start time if the start time of the observation is
  1194 + // smaller than the current start time
  1195 + if (obs->gti().tstart() < tstart) {
  1196 + tstart = obs->gti().tstart();
  1197 + }
  1198 +
  1199 + } // endif: observation was valid
  1200 +
  1201 + } // endfor: looped over observations
  1202 +
  1203 + // Log action
  1204 + log_string(NORMAL, "Setting start time from observations");
  1205 +
  1206 + } // endelse: no valid time specified
  1207 +
  1208 + // Return start time
  1209 + return tstart;
  1210 +}
  1211 +
  1212 +
  1213 +/***********************************************************************//**
  1214 + * @brief Get stop time
  1215 + *
  1216 + * @return Stop time
1044 1217 ***************************************************************************/
1045   -void ctfindvar::write_srchist(void)
  1218 +GTime ctfindvar::get_tstop(void)
  1219 +{
  1220 + // Initialise stop time with small unreasonable Julian date
  1221 + GTime tstop("JD 0");
  1222 +
  1223 + // If tmax parameter is valid then set the stop time
  1224 + if ((*this)["tmax"].is_valid()) {
  1225 +
  1226 + // Query tmax parameter
  1227 + tstop = (*this)["tmax"].time();
  1228 +
  1229 + // Log action
  1230 + log_string(NORMAL, "Setting stop time \"tmax\" parameter");
  1231 +
  1232 + } // endif: tmax parameter was valid
  1233 +
  1234 + // ... otherwise set stop time from observations
  1235 + else {
  1236 +
  1237 + // Loop over all observations
  1238 + for (int k = 0; k < m_obs.size(); ++k) {
  1239 +
  1240 + // Get CTA observation
  1241 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[k]);
  1242 +
  1243 + // Continue only if observation is valid
  1244 + if (obs != NULL) {
  1245 +
  1246 + // Update stop time if the stop time of the observation is
  1247 + // larger than the current stop time
  1248 + if (obs->gti().tstop() > tstop) {
  1249 + tstop = obs->gti().tstop();
  1250 + }
  1251 +
  1252 + } // endif: observation was valid
  1253 +
  1254 + } // endfor: looped over observations
  1255 +
  1256 + // Log action
  1257 + log_string(NORMAL, "Setting stop time from observations");
  1258 +
  1259 + } // endelse: no valid time specified
  1260 +
  1261 + // Return stop time
  1262 + return tstop;
  1263 +}
  1264 +
  1265 +
  1266 +/***********************************************************************//**
  1267 + * @brief Write source histograms in FITS format
  1268 + *
  1269 + * @param[in,out] fits FITS file
  1270 + *
  1271 + * This method creates two FITS binary table extensions and appends them to
  1272 + * the FITS file.
  1273 + *
  1274 + * The first extension is named "SOURCE SIGNIFICANCE" and contains the
  1275 + * significance of the variability as function of time bin for the pixel
  1276 + * with the maximum significance and for each source.
  1277 + *
  1278 + * The second extension is named "SOURCE POSITION" and contains the Right
  1279 + * Ascension and Declination for the maximum significance pixel and for each
  1280 + * source.
  1281 + ***************************************************************************/
  1282 +void ctfindvar::write_source_histograms(GFits& fits)
1046 1283 {
1047 1284 // Setup data storage for time and source significances
1048 1285 double tinterval = (*this)["tinterval"].real();
1049   - int nentries = (m_tstop-m_tstart) / tinterval + 0.5;
1050   - int nsources = m_pixsigsrc.shape()[0];
1051   - GNdarray time_info(2, nentries);
1052   - GNdarray max_pixel_info(nentries);
1053   - GNdarray src_info(nsources, nentries);
1054   -
1055   - // Store the full set of information
1056   - for (int i=0; i<nentries; i++) {
  1286 + int nbins = (m_tstop-m_tstart) / tinterval + 0.5;
  1287 + int nsrc = m_pixsigsrc.shape()[0];
  1288 + GNdarray time_info(2, nbins);
  1289 + GNdarray max_pixel_info(nbins);
  1290 + GNdarray src_info(nsrc, nbins);
  1291 +
  1292 + // Loop over time bins
  1293 + for (int i = 0; i < nbins; ++i) {
  1294 +
  1295 + // Set time bin mid point
1057 1296 GTime tstart_bin(m_tstart + i * tinterval);
1058 1297 GTime tstop_bin(m_tstart + (i+1) * tinterval);
1059 1298 GTime midpoint((tstop_bin.secs()+tstart_bin.secs())*0.5);
... ... @@ -1066,111 +1305,125 @@ void ctfindvar::write_srchist(void)
1066 1305 int index = time2inx(midpoint);
1067 1306  
1068 1307 // Loop over every source
1069   - for (int src=0; src<nsources; src++) {
1070   - if (index >= 0) {
1071   - src_info(src,i) = m_pixsigsrc(src,index);
1072   - } else {
1073   - src_info(src,i) = 0.0;
1074   - }
  1308 + for (int k = 0; k < nsrc; ++k) {
  1309 + src_info(k,i) = (index >= 0) ? m_pixsigsrc(k,index) : 0.0;
1075 1310 }
1076 1311  
1077   - // Fill max pixel information
1078   - if (index >= 0) {
1079   - max_pixel_info(i) = m_pixsigmax(index);
1080   - }
  1312 + // Fill maximum pixel information
  1313 + max_pixel_info(i) = (index >= 0) ? m_pixsigmax(index) : 0.0;
  1314 +
  1315 + } // endfor: looped over time bins
  1316 +
  1317 + // Create binary tables
  1318 + GFitsBinTable table_signif(nbins);
  1319 + GFitsBinTable table_pos(nsrc+1);
  1320 + table_signif.extname("SOURCE SIGNIFICANCE");
  1321 + table_pos.extname("SOURCE POSITION");
  1322 +
  1323 + // Create columns for time bins
  1324 + GFitsTableDoubleCol time_start("TSTART", nbins);
  1325 + GFitsTableDoubleCol time_stop("TSTOP", nbins);
  1326 + GFitsTableDoubleCol time_sigma("MAXSIGPIXEL", nbins);
  1327 + time_start.unit("MJD");
  1328 + time_stop.unit("MJD");
  1329 + time_sigma.unit("sigma");
  1330 +
  1331 + // Create columns for source position
  1332 + GFitsTableStringCol pos_name("SOURCE", nsrc+1, 20);
  1333 + GFitsTableDoubleCol pos_ra("RA", nsrc+1);
  1334 + GFitsTableDoubleCol pos_dec("DEC", nsrc+1);
  1335 + pos_ra.unit("deg");
  1336 + pos_dec.unit("deg");
  1337 +
  1338 + // Fill time bin columns
  1339 + for (int i = 0; i < nbins; ++i) {
  1340 + time_start(i) = time_info(0,i);
  1341 + time_stop(i) = time_info(1,i);
  1342 + time_sigma(i) = max_pixel_info(i);
1081 1343 }
1082 1344  
1083   - // Call the appropriate method for writing the data
1084   - write_srchist_fits(time_info, max_pixel_info, src_info);
1085   -}
  1345 + // Fill source columns
  1346 + pos_name(0) = "MAXSIGPIXEL";
  1347 + pos_ra(0) = m_max_sig_dir.ra_deg();
  1348 + pos_dec(0) = m_max_sig_dir.dec_deg();
1086 1349  
1087   -
1088   -/***********************************************************************//**
1089   - * @brief Write individual histograms in FITS format
1090   - *
1091   - * @param[in] time_info Start and stop time in MJD for each bin
1092   - * @param[in] src_info Significance distribution for each source
1093   -
1094   - * This method stores the extracted significance for each source at each time
1095   - * interval between tstart and tstop. Each time bin contains information on
1096   - * the start and stop time of that bin (in MJD) and the extracted variability
1097   - * significance of the source in that time bin.
1098   - ***************************************************************************/
1099   -void ctfindvar::write_srchist_fits(const GNdarray& time_info,
1100   - const GNdarray& max_pixel_info,
1101   - const GNdarray& src_info)
1102   -{
1103   - log_string(NORMAL, "Storing histograms in FITS format");
1104   -
1105   - // Create the table
1106   - int nrows = time_info.shape()[1];
1107   - GFitsBinTable table_signif(nrows);
1108   - GFitsBinTable table_dir(2);
1109   - table_signif.extname("SRCSIGNIF");
1110   - table_dir.extname("SRCPOSITION");
1111   -
1112   - // Append the start/stop time and max pixel columns
1113   - GFitsTableDoubleCol start_mjd("TSTART", nrows);
1114   - GFitsTableDoubleCol stop_mjd("TSTOP", nrows);
1115   - GFitsTableDoubleCol max_pix("MAXSIGPIXEL", nrows);
1116   - for (int i=0; i<nrows; i++) {
1117   - start_mjd(i) = time_info(0,i);
1118   - stop_mjd(i) = time_info(1,i);
1119   - max_pix(i) = max_pixel_info(i);
1120   - }
1121   - table_signif.append(start_mjd);
1122   - table_signif.append(stop_mjd);
1123   - table_signif.append(max_pix);
1124   -
1125   - // Store the position of the peak pixel
1126   - GFitsTableDoubleCol max_pix_dir("MAXSIGPIXEL",2);
1127   - max_pix_dir(0) = m_max_sig_dir.ra_deg();
1128   - max_pix_dir(1) = m_max_sig_dir.dec_deg();
1129   - table_dir.append(max_pix_dir);
  1350 + // Append columns to tables
  1351 + table_signif.append(time_start);
  1352 + table_signif.append(time_stop);
  1353 + table_signif.append(time_sigma);
1130 1354  
1131 1355 // Append the distribution for each source
1132   - int src_count = src_info.shape()[0];
1133   - for (int src=0; src<src_count; src++) {
1134   -
  1356 + for (int k = 0; k < nsrc; ++k) {
  1357 +
  1358 + // Initialise source name and direction
1135 1359 std::string src_name("SOURCE");
1136 1360 GSkyDir src_dir;
1137 1361  
1138 1362 // Get the name from the input model file
1139   - if (m_inmodel.size() > 0) {
1140   - src_name = m_inmodel[src]->name();
1141   -
1142   - // Get source position
1143   - GModelSky* model = dynamic_cast<GModelSky*>(m_inmodel[src]);
1144   - GModelSpatial* model_spatial = model->spatial();
1145   - GSkyRegion* model_region = model_spatial->region();
1146   - src_dir = dynamic_cast<GSkyRegionCircle*>(model_region)->centre();
1147   - }
1148   - // ... otherwise set the position from xsrc,ysrc
  1363 + if (m_inmodel.size() > k) {
  1364 +
  1365 + // Extract the model
  1366 + GModelSky* model = dynamic_cast<GModelSky*>(m_inmodel[k]);
  1367 +
  1368 + // If model is a sky model then extract the model centre. Note
  1369 + // that this only works for a model that results a GSkyRegionCircle
  1370 + // object as region.
  1371 + if (model != NULL) {
  1372 +
  1373 + // Get the source position
  1374 + GModelSpatial* spatial = model->spatial();
  1375 + GSkyRegion* region = spatial->region();
  1376 + GSkyRegionCircle* circle = dynamic_cast<GSkyRegionCircle*>(region);
  1377 +
  1378 + // If region circle is valid then extract source name and
  1379 + // sky direction
  1380 + if (circle != NULL) {
  1381 + src_name = m_inmodel[k]->name();
  1382 + src_dir = circle->centre();
  1383 + }
  1384 +
  1385 + } // endif: model was a sky model
  1386 +
  1387 + } // endif: there were enough models in the input model
  1388 +
  1389 + // ... otherwise set the position from the xsrc and ysrc parameters
1149 1390 else {
1150 1391 if ((*this)["coordsys"].string() == "CEL") {
1151   - src_dir.radec_deg((*this)["xsrc"].real(), (*this)["yref"].real());
1152   - } else {
1153   - src_dir.lb_deg((*this)["xsrc"].real(), (*this)["yref"].real());
  1392 + src_dir.radec_deg((*this)["xsrc"].real(), (*this)["ysrc"].real());
  1393 + }
  1394 + else {
  1395 + src_dir.lb_deg((*this)["xsrc"].real(), (*this)["ysrc"].real());
1154 1396 }
1155 1397 }
1156   -
1157   - // Create and fill a new column
1158   - GFitsTableDoubleCol src_sig_col(src_name, nrows);
1159   - for (int i=0; i<nrows; i++) {
1160   - src_sig_col(i) = src_info(src, i);
  1398 +
  1399 + // Create new time column for source
  1400 + GFitsTableDoubleCol column(src_name, nbins);
  1401 + column.unit("sigma");
  1402 +
  1403 + // Fill column
  1404 + for (int i = 0; i < nbins; ++i) {
  1405 + column(i) = src_info(k, i);
1161 1406 }
1162   - table_signif.append(src_sig_col);
1163 1407  
1164   - // Store the positions of the source
1165   - GFitsTableDoubleCol src_dir_col(src_name,2);
1166   - src_dir_col(0) = src_dir.ra_deg();
1167   - src_dir_col(1) = src_dir.dec_deg();
1168   - table_dir.append(src_dir_col);
1169   - }
  1408 + // Append column to time table
  1409 + table_signif.append(column);
  1410 +
  1411 + // Store source infirmation
  1412 + pos_name(k+1) = src_name;
  1413 + pos_ra(k+1) = src_dir.ra_deg();
  1414 + pos_dec(k+1) = src_dir.dec_deg();
  1415 +
  1416 + } // endfor: looped over sources
1170 1417  
1171   - // Store information regarding the source positions
  1418 + // Append columns to tables
  1419 + table_pos.append(pos_name);
  1420 + table_pos.append(pos_ra);
  1421 + table_pos.append(pos_dec);
1172 1422  
1173   - // Write the tables to a file
1174   - m_outfile.append(table_signif);
1175   - m_outfile.append(table_dir);
  1423 + // Append tables to FITS file
  1424 + fits.append(table_signif);
  1425 + fits.append(table_pos);
  1426 +
  1427 + // Return
  1428 + return;
1176 1429 }
... ...
src/ctfindvar/ctfindvar.hpp
  View file @0be0a1a
... ... @@ -68,27 +68,22 @@ public:
68 68  
69 69 protected:
70 70 // Protected methods
71   - void init_members(void);
72   - void copy_members(const ctfindvar& app);
73   - void free_members(void);
74   - void get_parameters(void);
75   - void init_cube(void);
76   - void init_gtis(void);
77   - void fill_cube(GCTAObservation* obs);
78   - void fill_model_sig_pos(const GSkyDir& dir_pix);
79   - void fill_alpha_vector(const int& pix_number,
80   - std::vector<double>& alpha_vector);
81   - void get_variability_sig(const int& pix_number,
82   - const int& nbins,
83   - GNdarray& sig_histogram);
84   - double gti_overlap(const GGti& gti1,
85   - const GGti& gti2);
86   -
87   - // Writing methods for individual source histograms
88   - void write_srchist(void);
89   - void write_srchist_fits(const GNdarray& time_info,
90   - const GNdarray& max_pixel_info,
91   - const GNdarray& src_info);
  71 + void init_members(void);
  72 + void copy_members(const ctfindvar& app);
  73 + void free_members(void);
  74 + void get_parameters(void);
  75 + void init_gtis(void);
  76 + void create_cube(void);
  77 + void fill_cube(GCTAObservation* obs);
  78 + void analyse_cube(void);
  79 + std::vector<int> get_pixels(void);
  80 + GNdarray get_variability_sig(const int& ipix);
  81 + std::vector<double> get_alphas(const int& ipix) const;
  82 + double gti_overlap(const GGti& gti1, const GGti& gti2) const;
  83 + GTime get_tstart(void);
  84 + GTime get_tstop(void);
  85 + GModelSky sky_model(const GSkyDir& dir) const;
  86 + void write_source_histograms(GFits& fits);
92 87  
93 88 // Protected members
94 89 GSkyMap m_counts; //!< Counts for each time interval
... ... @@ -97,7 +92,6 @@ protected:
97 92 GSkyDir m_max_sig_dir; //!< Sky direction associated with maximum significance
98 93 double m_minoff; //!< Minimum counts for use in significance calculation
99 94 double m_sig_threshold; //!< Minimum significance required to set source as variable
100   - GFits m_outfile; //!< FITS object to store output
101 95 GSkyMap m_peaksigmap; //!< Skymap holding the maximum significance
102 96 GNdarray m_pixsigsrc; //!< Store distributions of the source significances
103 97 GNdarray m_pixsigmax; //!< Store distribution for pixel with max significance
... ...
src/ctfindvar/ctfindvar.par
  View file @0be0a1a
... ... @@ -16,22 +16,26 @@
16 16 #
17 17 # General parameters
18 18 #===================
19   -inobs, f, a, events.fits,,, "Input event list, counts cube or observation definition XML file"
  19 +inobs, f, a, events.fits,,, "Input event list or observation definition XML file"
20 20 inmodel, f, h, NONE,,, "Input model definition file for extracting source positions"
21   -outcube, f, h, NONE,,, "Filename for counts cube"
  21 +outcube, f, h, NONE,,, "Output counts cube file"
  22 +outmap, f, a, map.fits,,, "Output map file"
  23 +outmodel, f, a, model.xml,,, "Output model definition file"
22 24 caldb, s, a, prod2,,, "Calibration database"
23 25 irf, s, a, South_0.5h,,, "Instrument response function"
24   -prefix, s, a, findvar_,,, "Prefix for output files"
25 26  
26 27 #
27 28 # Variability search parameters
28 29 #==============================
29   -minoff, r, h, 0,,, "Minimum required events in off bin for computing Noff and alpha"
30   -emin, r, a, 0,,, "Lower energy limit (units of TeV; 0 => no cut)"
31   -emax, r, a, 0,,, "Upper energy limit (units of TeV; 0 => no cut)"
32   -sig_thr, r, h, 4.5,,, "Significance threshold for variability detection"
33   -smoothkrnl, s, h, NONE, GAUSSIAN|DISK|NONE,, "Smoothing kernel for counts maps"
34   -smoothpar, r, h, 0.05,,, "Parameter for the smoothing kernel"
  30 +coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
  31 +xsrc, r, a, 83.63,0,360, "First coordinate of the source of interest (RA or galactic l)"
  32 +ysrc, r, a, 22.01,-90,90, "Second coordinate of the source of interest (DEC or galactic b)"
  33 +emin, r, a, 0.1,,, "Lower energy limit (TeV)"
  34 +emax, r, a, 100.0,,, "Upper energy limit (TeV)"
  35 +threshold, r, h, 4.5,,, "Significance threshold for variability detection"
  36 +minoff, r, h, 0,,, "Minimum required events in off bin for computing Noff and alpha"
  37 +smooth_kernel, s, h, NONE, GAUSSIAN|DISK|NONE,, "Smoothing kernel for counts maps"
  38 +smooth_rad, r, h, 0.05,,, "Smoothing kernel radius (deg)"
35 39  
36 40 #
37 41 # Time binning parameters
... ... @@ -47,12 +51,9 @@ usepnt, b, h, no,,, &quot;Use pointing instead of xref/yref parameters?&quot;
47 51 nxpix, i, a, 200,,, "Size of the X axis in pixels"
48 52 nypix, i, a, 200,,, "Size of the Y axis in pixels"
49 53 binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
50   -coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
51 54 proj, s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
52 55 xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
53 56 yref, r, a, 22.51,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
54   -xsrc, r, a, 83.63,0,360, "First coordinate of the source of interest (RA or galactic l)"
55   -ysrc, r, a, 22.51,-90,90, "Second coordinate of the source of interest (DEC or galactic b)"
56 57  
57 58 #
58 59 # Standard parameters
... ...
test/test_ctfindvar.py
  View file @0be0a1a
... ... @@ -71,8 +71,9 @@ class Test(test):
71 71  
72 72 # Setup ctfindvar command
73 73 cmd = ctfindvar+' inobs="'+self._events+'"'+\
  74 + ' outmap="ctfindvar_cmd1.fits"'+\
  75 + ' outmodel="ctfindvar_cmd1.xml"'+\
74 76 ' logfile="ctfindvar_cmd1.log" chatter=1'+\
75   - ' prefix=ctfindvar_test_ '+\
76 77 ' tmin=NONE tmax=NONE tinterval=20'+\
77 78 ' emin=0.2 emax=1.1'+\
78 79 ' nxpix=200 nypix=200 binsz=0.05'+\
... ... @@ -80,12 +81,28 @@ class Test(test):
80 81 ' xref=83.63 yref=22.51'+\
81 82 ' xsrc=83.63 ysrc=22.51'+\
82 83 ' caldb=prod2 irf=South_0.5h'
83   -
84   -
  84 +
85 85 # Check if execution was successful
86 86 self.test_value(self._execute(cmd), 0,
87 87 'Check successful execution from command line')
88 88  
  89 + # Setup ctfindvar command
  90 + cmd = ctfindvar+' inobs="event_file_that_does_not_exist.fits"'+ \
  91 + ' outmap="ctfindvar_cmd2.fits"'+ \
  92 + ' outmodel="ctfindvar_cmd2.xml"'+ \
  93 + ' logfile="ctfindvar_cmd2.log" chatter=1'+ \
  94 + ' tmin=NONE tmax=NONE tinterval=20'+ \
  95 + ' emin=0.2 emax=1.1'+ \
  96 + ' nxpix=200 nypix=200 binsz=0.05'+ \
  97 + ' proj=CAR coordsys=CEL'+ \
  98 + ' xref=83.63 yref=22.51'+ \
  99 + ' xsrc=83.63 ysrc=22.51'+ \
  100 + ' caldb=prod2 irf=South_0.5h'
  101 +
  102 + # Check if execution was successful
  103 + self.test_assert(self._execute(cmd, success=False) != 0,
  104 + 'Check invalid input file when executed from command line')
  105 +
89 106 # Check ctfindvar --help
90 107 self._check_help(ctfindvar)
91 108  
... ... @@ -101,11 +118,8 @@ class Test(test):
101 118 # Allocate ctfindvar
102 119 ctfindvar = ctools.ctfindvar()
103 120  
104   - #setup ctfindvar
  121 + # Setup ctfindvar
105 122 ctfindvar['inobs'] = self._events
106   - ctfindvar['logfile'] = 'ctfindvar_py0.log'
107   - ctfindvar['chatter'] = 1
108   - ctfindvar['prefix'] = 'ctfindvar_test_python_'
109 123 ctfindvar['tmin'] = 'NONE'
110 124 ctfindvar['tmax'] = 'NONE'
111 125 ctfindvar['tinterval'] = 20
... ... @@ -122,32 +136,39 @@ class Test(test):
122 136 ctfindvar['ysrc'] = 22.51
123 137 ctfindvar['caldb'] = 'prod2'
124 138 ctfindvar['irf'] = 'South_0.5h'
  139 + ctfindvar['outmap'] = 'ctfindvar_py0.fits'
  140 + ctfindvar['outmodel'] = 'ctfindvar_py0.xml'
  141 + ctfindvar['logfile'] = 'ctfindvar_py0.log'
  142 + ctfindvar['chatter'] = 1
125 143  
126   - ## Run ctfindvar tool
  144 + # Run ctfindvar tool
127 145 ctfindvar.logFileOpen()
128 146 ctfindvar.run()
129 147  
130 148 # Check that saving does not nothing
131   - ctfindvar['logfile'] = 'ctfindvar_py1.log'
  149 + ctfindvar['outmap'] = 'ctfindvar_py1.fits'
  150 + ctfindvar['outmodel'] = 'ctfindvar_py1.xml'
  151 + ctfindvar['logfile'] = 'ctfindvar_py1.log'
132 152 ctfindvar.logFileOpen()
133 153 ctfindvar.save()
134 154  
135 155 # Check result
136   - self._check_result_file('ctfindvar_test_python_signifmap.fits')
  156 + self._check_result_file('ctfindvar_py1.fits')
137 157  
138 158 # Copy ctfindvar tool
139 159 cpy_tool = ctfindvar
140 160  
141 161 # Run copy of ctfindvar tool again
142   - cpy_tool['logfile'] = 'ctfindvar_py2.log'
143   - cpy_tool['chatter'] = 3
144   - cpy_tool['prefix'] = 'ctfindvar_test_copytool_'
  162 + cpy_tool['outmap'] = 'ctfindvar_py2.fits'
  163 + cpy_tool['outmodel'] = 'ctfindvar_py2.xml'
  164 + cpy_tool['logfile'] = 'ctfindvar_py2.log'
  165 + cpy_tool['chatter'] = 3
145 166 cpy_tool.logFileOpen()
146 167 cpy_tool.run()
147 168 cpy_tool.save()
148 169  
149 170 # Check result
150   - self._check_result_file('ctfindvar_test_copytool_signifmap.fits')
  171 + self._check_result_file('ctfindvar_py2.fits')
151 172  
152 173 # Check that clearing does not lead to an exception or segfault
153 174 ctfindvar.clear()
... ... @@ -158,17 +179,16 @@ class Test(test):
158 179 # Check ctfindvar result
159 180 def _check_result_file(self, filename):
160 181 """
161   - Check content of tool
  182 + Check ctfindvar result
162 183  
163 184 Parameters
164 185 ----------
165   - tool : `~ctools.ctfindvar`
166   - ctfindvar instance
  186 + filename : str
  187 + Output file
167 188 """
168   - print "checking if you amazing file is valid"
169   - # Read variability file
170   - fits = gammalib.GFits(filename)
171   - self.test_value(fits.size(), 3,
172   - 'Check for 3 rows in srcsig file')
  189 + # Read variability file
  190 + fits = gammalib.GFits(filename)
  191 + self.test_value(fits.size(), 3, 'Check for 3 extensions in output file')
  192 +
173 193 # Return
174 194 return
... ...