Eirik Hatlen / ctools

Commit 90c2e5d8972603060ae47fb0b749b6812890ada8

Authored by Jürgen Knödlseder
1 parent d135a149

Integrate ctfindvar tool

Inline Side-by-side
Showing 5 changed files with 364 additions and 336 deletions
.gitignore
  View file @90c2e5d
... ... @@ -57,6 +57,7 @@ pyext/ctools/tools.py
57 57 pyext/ctools/tools_wrap.cpp
58 58 cscripts/__init__.py
59 59 src/ctbin/ctbin
  60 +src/ctfindvar/ctfindvar
60 61 src/ctlike/ctlike
61 62 src/cttsmap/cttsmap
62 63 src/ctmodel/ctmodel
... ...
src/ctfindvar/ctfindvar.cpp
  View file @90c2e5d
1 1 /***************************************************************************
2   - * ctfindvar - search time variability tool *
  2 + * ctfindvar - Time variability search tool *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2018 by Simon Bonnefoy *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctfindvar.cpp
23   - * @brief search time variability tool implementation
  23 + * @brief Time variability search tool implementation
24 24 * @author Simon Bonnefoy
25 25 */
26 26  
... ... @@ -40,7 +40,8 @@
40 40  
41 41 /* __ Method name definitions ____________________________________________ */
42 42 #define G_FILL_CUBE "ctfindvar::fill_cube(GCTAObservation*)"
43   -#define G_FILL_ALPHA_VECTOR "ctfindvar::fill_alpha_vector(const int&, vector<double>&)"
  43 +#define G_FILL_ALPHA_VECTOR "ctfindvar::fill_alpha_vector(const int&, "\
  44 + "vector<double>&)"
44 45 #define G_INIT_CUBE "ctfindvar::init_cube(void)"
45 46 #define G_INIT_GTIS "ctfindvar::init_gtis(void)"
46 47  
... ... @@ -60,10 +61,9 @@
60 61 /***********************************************************************//**
61 62 * @brief Void constructor
62 63 *
63   - * Constructs empty search time variability tool.
  64 + * Constructs empty time variability search tool.
64 65 ***************************************************************************/
65   -ctfindvar::ctfindvar(void) :
66   - ctobservation(CTFINDVAR_NAME, VERSION)
  66 +ctfindvar::ctfindvar(void) : ctobservation(CTFINDVAR_NAME, VERSION)
67 67 {
68 68 // Initialise members
69 69 init_members();
... ... @@ -78,10 +78,10 @@ ctfindvar::ctfindvar(void) :
78 78 *
79 79 * param[in] obs Observation container.
80 80 *
81   - * Constructs search time variability tool from an observation container.
  81 + * Constructs time variability search tool from an observation container.
82 82 ***************************************************************************/
83 83 ctfindvar::ctfindvar(const GObservations& obs) :
84   - ctobservation(CTFINDVAR_NAME, VERSION, obs)
  84 + ctobservation(CTFINDVAR_NAME, VERSION, obs)
85 85 {
86 86 // Initialise members
87 87 init_members();
... ... @@ -97,11 +97,11 @@ ctfindvar::ctfindvar(const GObservations&amp; obs) :
97 97 * @param[in] argc Number of arguments in command line.
98 98 * @param[in] argv Array of command line arguments.
99 99 *
100   - * Constructs search time variability tool using command line arguments for user
101   - * parameter setting.
  100 + * Constructs time variability search tool using command line arguments for
  101 + * user parameter setting.
102 102 ***************************************************************************/
103 103 ctfindvar::ctfindvar(int argc, char *argv[]) :
104   - ctobservation(CTFINDVAR_NAME, VERSION, argc, argv)
  104 + ctobservation(CTFINDVAR_NAME, VERSION, argc, argv)
105 105 {
106 106 // Initialise members
107 107 init_members();
... ... @@ -116,11 +116,10 @@ ctfindvar::ctfindvar(int argc, char *argv[]) :
116 116 *
117 117 * @param[in] app search time variability tool.
118 118 *
119   - * Constructs search time variability tool from another search time variability
120   - * tool.
  119 + * Constructs time variability search tool from another time variability
  120 + * search tool.
121 121 ***************************************************************************/
122   -ctfindvar::ctfindvar(const ctfindvar& app) :
123   - ctobservation(app)
  122 +ctfindvar::ctfindvar(const ctfindvar& app) : ctobservation(app)
124 123 {
125 124 // Initialise members
126 125 init_members();
... ... @@ -157,10 +156,10 @@ ctfindvar::~ctfindvar(void)
157 156 /***********************************************************************//**
158 157 * @brief Assignment operator
159 158 *
160   - * @param[in] app search time variability tool.
161   - * @return search time variability tool.
  159 + * @param[in] app Time variability search tool.
  160 + * @return Time variability search tool.
162 161 *
163   - * Assigns search time variability tool.
  162 + * Assigns time variability search tool.
164 163 ***************************************************************************/
165 164 ctfindvar& ctfindvar::operator=(const ctfindvar& app)
166 165 {
... ... @@ -193,9 +192,9 @@ ctfindvar&amp; ctfindvar::operator=(const ctfindvar&amp; app)
193 192 ==========================================================================*/
194 193  
195 194 /***********************************************************************//**
196   - * @brief Clear search time variability tool
  195 + * @brief Clear time variability search tool
197 196 *
198   - * Clears search time variability tool.
  197 + * Clears time variability search tool.
199 198 ***************************************************************************/
200 199 void ctfindvar::clear(void)
201 200 {
... ... @@ -221,7 +220,7 @@ void ctfindvar::clear(void)
221 220  
222 221  
223 222 /***********************************************************************//**
224   - * @brief Run search time variability tool
  223 + * @brief Run time variability search tool
225 224 ***************************************************************************/
226 225 void ctfindvar::run(void)
227 226 {
... ... @@ -258,6 +257,7 @@ void ctfindvar::run(void)
258 257 (*this)["smoothpar"].real());
259 258 }
260 259  
  260 + // Get number of maps
261 261 const int nbins = m_counts.nmaps();
262 262  
263 263 // creating GSkyDir to get the position of the source and each pixels
... ... @@ -265,7 +265,7 @@ void ctfindvar::run(void)
265 265 // Extract which pixel in the map each source is located in
266 266 std::vector<int> srcInxPix;
267 267 if (m_inmodel.size() > 0) {
268   - for (int mod=0; mod<m_inmodel.size(); mod++) {
  268 + for (int mod = 0; mod < m_inmodel.size(); ++mod) {
269 269  
270 270 // Extract the model
271 271 GModelSky* model = dynamic_cast<GModelSky*>(m_inmodel[mod]);
... ... @@ -380,196 +380,50 @@ void ctfindvar::run(void)
380 380  
381 381  
382 382 /***********************************************************************//**
383   - * @brief Obtain the significance of variability in a given bin
384   - *
385   - * @param[in] pix_number Pixel in the skymap
386   - * @param[in] nbins Number of bins in the skymap
387   - * @param[out] sig_histogram Histogram for storing returned significances
388   - *
389   - * Significance is computed according to Li & Ma equation 17:
390   - * see: https://doi.org/10.1086/161295
391   - * Some modification is made in order to handle the case where Non = 0,
392   - * but Noff != 0.
393   - ***************************************************************************/
394   -void ctfindvar::get_variability_sig(const int& pix_number,
395   - const int& nbins,
396   - GNdarray& sig_histogram)
397   -{
398   - std::vector<bool> accepted_bin_bckg_vector(nbins, true);
399   - bool background_validated = false;
400   - double non, noff;
401   - double alpha, sig;
402   -
403   - // Fill the alpha values for each vector
404   - std::vector<double> alpha_vector(nbins, 0.0);
405   - fill_alpha_vector(pix_number, alpha_vector);
406   -
407   - // Loop over pixels until all background pixls are removed
408   - while (background_validated == false) {
409   - background_validated = true;
410   -
411   - // Loop over all the GTIs of the pixel
412   - for (int i = 0; i < nbins; i++) {
413   -
414   - // The GTI is discared from bckg calculation and not checked again.
415   - if (!accepted_bin_bckg_vector[i]) {
416   - continue;
417   - }
418   - // Check if bin fails minoff check or no observations overlap it
419   - else if (m_counts(pix_number, i) < m_minoff || alpha_vector[i] == 0) {
420   - accepted_bin_bckg_vector[i] = false;
421   - continue;
422   - }
423   -
424   - noff = 0.0;
425   - alpha = 0.0;
426   - // Gor one GTI selected (i), looping over all the others (j).
427   - for (int j = 0; j < nbins; j++) {
428   - if (j != i && accepted_bin_bckg_vector[j] == 1) {
429   - noff += m_counts(pix_number, j);
430   - alpha += alpha_vector[j];
431   - }
432   - }
433   -
434   - // The background is averaged on the number of bins -1
435   - non = m_counts(pix_number, i);
436   - alpha = alpha_vector[i]/alpha;
437   -
438   - // Compute sensitivity in Gaussian sigma using Li & Ma
439   - double alpha1 = alpha + 1.0;
440   - double ntotal = non+noff;
441   - double arg1 = non/ntotal;
442   - double arg2 = noff/ntotal;
443   - if (noff == 0.0) {
444   - sig = 0.0;
445   - } else if (non > 0.0) {
446   - double term1 = non * std::log((alpha1/alpha)*arg1);
447   - double term2 = noff * std::log(alpha1*arg2);
448   - sig = std::sqrt(2.0 * (term1 + term2));
449   - } else {
450   - sig = std::sqrt(2.0 * noff * std::log(alpha1));
451   - }
452   -
453   - // Specify the sign of the significance
454   - sig *= (non < alpha*noff) ? -1.0 : 1.0;
455   -
456   - if (alpha == 0.0) {
457   - sig = 0.0;
458   - }
459   -
460   - // Update the significance
461   - sig_histogram(i) = sig;
462   -
463   - // If bin is significant, it's removed from bckg and we loop again.
464   - if (sig > m_sig_threshold) {
465   - accepted_bin_bckg_vector[i] = false;
466   - background_validated = false;
467   - }
468   - }
469   -
470   - }
471   -}
472   -
473   -/***********************************************************************//**
474   - * @brief Save the current sky direction to model if significance > thr
475   - *
476   - * @param[in] dir_pix pixel sky direction
477   - * @param[out] m_model_sig_above_thr model position with sig>thr
478   - ***************************************************************************/
479   -void ctfindvar::fill_model_sig_pos(const GSkyDir& dir_pix)
480   -{
481   - //Set spatial components
482   - GModelSpatialPointSource spatial(dir_pix);
483   -
484   - //Set fake spectral component (powerlaw with 1% of Crab)
485   - GModelSpectralPlaw spectral(5.7e-18, -2.48, GEnergy(0.3, "TeV"));
486   -
487   - //Creat a sky model
488   - GModelSky model(spatial, spectral);
489   -
490   - //Extracting the coordinate to name the model
491   - double ra = dir_pix.ra_deg();
492   - double dec = dir_pix.dec_deg();
493   - std::string name = "Sky_direction_" + gammalib::str(ra, 4) +
494   - "_" + gammalib::str(dec, 4);
495   - model.name(name);
496   -
497   - //Add the model for this pixel to the output mode
498   - m_model_above_thr.append(model);
499   -
500   -}
501   -/***********************************************************************//**
502   - * @brief Get the map index associated with a given time
  383 + * @brief Save peak significance map and significance distributions
503 384 *
504   - * @param[in] pix_number Spatial pixel
505   - * @param[out] alpha_vector Vector containing effective exposure
  385 + * Saves the peak significance and source signficance distributions.
  386 + * Optionally, the generated counts cube can also be saved if 'outcube' has
  387 + * been specified.
506 388 ***************************************************************************/
507   -void ctfindvar::fill_alpha_vector(const int& pix_number,
508   - std::vector<double>& alpha_vector)
  389 +void ctfindvar::save(void)
509 390 {
510   - // Pixel position
511   - GSkyDir pix_dir = m_counts.inx2dir(pix_number);
512   -
513   - // Loop over all observations
514   - for (int i = 0; i < m_obs.size(); i++) {
515   - // Skip if observation does not overlap with this pixel position
516   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
517   - GCTARoi roi = obs->roi();
518   - GSkyDir roi_centre = roi.centre().dir();
519   -
520   - if (roi_centre.dist_deg(pix_dir) > roi.radius()) {
521   - continue;
522   - }
523   -
524   - // Convert sky direction to instrument direction
525   - GCTAInstDir instdir = obs->pointing().instdir(pix_dir);
526   -
527   - // Extract the good time intervals of the observation
528   - const GGti& obs_gti(obs->gti());
  391 + log_header1(TERSE, "Saving results");
529 392  
530   - // Loop over all time bins
531   - for (int j = 0; j < m_gti.size(); j++) {
  393 + // Filenames
  394 + std::string prefix((*this)["prefix"].string());
532 395  
533   - // Make sure observation overlaps with this time interval
534   - double exposure = gti_overlap(m_gti[j], obs_gti);
535   - if (exposure > 0.0) {
  396 + // Write counts cube
  397 + if ((*this)["outcube"].is_valid()) {
  398 + GFilename outcube = (*this)["outcube"].filename();
  399 + log_value(TERSE, "Saving counts cube", outcube);
  400 + m_counts.save(outcube, (*this)["clobber"].boolean());
  401 + }
536 402  
537   - // Get IRF response
538   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
539   - (obs->response());
  403 + // Create a FITS file for storing the output
  404 + m_outfile.clear();
540 405  
541   - // Throw an exception if observation has no instrument response function
542   - if (rsp == NULL) {
543   - std::string msg =
544   - "No response information available for "+
545   - get_obs_header(obs)+" to compute IRF background. "
546   - "Please specify response information or use "
547   - "another background subtraction method.";
548   - throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg);
549   - }
  406 + // 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");
550 409  
551   - // Get IRF background template
552   - const GCTABackground* bkg = rsp->background();
  410 + // Write individual source distributions
  411 + write_srchist();
553 412  
554   - // Throw an exception if observation has no IRF background template
555   - if (bkg == NULL) {
556   - std::string msg =
557   - "No IRF background template found in instrument "
558   - "response function for "+
559   - get_obs_header(obs)+". Please specify an instrument "
560   - "response function containing a background template.";
561   - throw GException::invalid_value(G_FILL_ALPHA_VECTOR, msg);
562   - }
  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());
563 417  
564   - // Multiply background rate with livetime and solid angle
565   - #pragma omp critical(ctfindvar_alpha_vector)
566   - alpha_vector[j] += exposure * bkg->rate_ebin(instdir, m_emin, m_emax);
567   - }
568   - }
  418 + //Save the model file with position with significance
  419 + if (m_model_above_thr.size() > 2) {
  420 + GFilename model_filename = prefix + "model_significant_directions.xml";
  421 + m_model_above_thr.save(model_filename);
569 422 }
570 423  
  424 + // Return
571 425 return;
572   -}
  426 +}
573 427  
574 428  
575 429 /***********************************************************************//**
... ... @@ -611,53 +465,6 @@ GGti ctfindvar::inx2gti(const int&amp; index)
611 465 }
612 466  
613 467  
614   -/***********************************************************************//**
615   - * @brief Save peak significance map and significance distributions
616   - *
617   - * Saves the peak significance and source signficance distributions.
618   - * Optionally, the generated counts cube can also be saved if 'outcube' has
619   - * been specified.
620   - ***************************************************************************/
621   -void ctfindvar::save(void)
622   -{
623   - log_header1(TERSE, "Saving results");
624   -
625   - // Filenames
626   - std::string prefix((*this)["prefix"].string());
627   -
628   - // Write counts cube
629   - if ((*this)["outcube"].is_valid()) {
630   - GFilename outcube = (*this)["outcube"].filename();
631   - log_value(TERSE, "Saving counts cube", outcube);
632   - m_counts.save(outcube, (*this)["clobber"].boolean());
633   - }
634   -
635   - // Create a FITS file for storing the output
636   - m_outfile.clear();
637   -
638   - // Write the most significant values for each pixel
639   - log_string(TERSE, "Appending maximum significances map to output file");
640   - m_peaksigmap.write(m_outfile, "PEAKSIGMAP");
641   -
642   - // Write individual source distributions
643   - write_srchist();
644   -
645   - // Write the output file
646   - GFilename outfilename = prefix + "signifmap.fits";
647   - log_value(TERSE, "Saving results to ", outfilename);
648   - m_outfile.saveto(outfilename, (*this)["clobber"].boolean());
649   -
650   - //Save the model file with position with significance
651   - if (m_model_above_thr.size() > 2) {
652   - GFilename model_filename = prefix + "model_significant_directions.xml";
653   - m_model_above_thr.save(model_filename);
654   - }
655   -
656   - // Return
657   - return;
658   -}
659   -
660   -
661 468 /*==========================================================================
662 469 = =
663 470 = Private methods =
... ... @@ -673,10 +480,17 @@ void ctfindvar::init_members(void)
673 480 m_counts.clear();
674 481 m_gti.clear();
675 482 m_inmodel.clear();
  483 + m_max_sig_dir.clear();
  484 + m_minoff = 0.0;
  485 + m_sig_threshold = 0.0;
  486 + m_outfile.clear();
676 487 m_peaksigmap.clear();
677 488 m_pixsigsrc.clear();
678 489 m_tstart.clear();
679 490 m_tstop.clear();
  491 + m_emin.clear();
  492 + m_emax.clear();
  493 + m_model_above_thr.clear();
680 494  
681 495 // Return
682 496 return;
... ... @@ -690,8 +504,21 @@ void ctfindvar::init_members(void)
690 504 ***************************************************************************/
691 505 void ctfindvar::copy_members(const ctfindvar& app)
692 506 {
693   - // Copy attributes
694   - // TODO: Your code goes here
  507 + // Copy members
  508 + m_counts = app.m_counts;
  509 + m_gti = app.m_gti;
  510 + m_inmodel = app.m_inmodel;
  511 + m_max_sig_dir = app.m_max_sig_dir;
  512 + m_minoff = app.m_minoff;
  513 + m_sig_threshold = app.m_sig_threshold;
  514 + m_outfile = app.m_outfile;
  515 + m_peaksigmap = app.m_peaksigmap;
  516 + m_pixsigsrc = app.m_pixsigsrc;
  517 + m_tstart = app.m_tstart;
  518 + m_tstop = app.m_tstop;
  519 + m_emin = app.m_emin;
  520 + m_emax = app.m_emax;
  521 + m_model_above_thr = app.m_model_above_thr;
695 522  
696 523 // Return
697 524 return;
... ... @@ -761,6 +588,113 @@ void ctfindvar::get_parameters(void)
761 588  
762 589  
763 590 /***********************************************************************//**
  591 + * @brief Initialize counts cube
  592 + *
  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
  611 + *
  612 + * Sets up the required event cube and time interval storage objects
  613 + * associated with this object
  614 + ***************************************************************************/
  615 +void ctfindvar::init_gtis(void)
  616 +{
  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
  620 +
  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 + }
  636 +
  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 {
  642 +
  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 + }
  653 +
  654 + // Fill the number of good time intervals
  655 + double tstart_sec = m_tstart.secs();
  656 + double tstop_sec = m_tstop.secs();
  657 + int bins = (tstop_sec - tstart_sec) / tinterval + 0.5;
  658 +
  659 + // 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
  668 + if (m_tstart == m_tstop) {
  669 + throw GException::invalid_value(G_INIT_GTIS,
  670 + "Start time is equal to stop time");
  671 + } else if (bins < 2) {
  672 + std::string msg = "Method requires at least two time bins (" +
  673 + gammalib::str(bins) + " bins found).";
  674 + throw GException::invalid_value(G_INIT_GTIS, msg);
  675 + }
  676 +
  677 + for (int i=0; i<bins; i++) {
  678 +
  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);
  683 +
  684 + // 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 + }
  693 + }
  694 +}
  695 +
  696 +
  697 +/***********************************************************************//**
764 698 * @brief Fill the cube from the events in the map
765 699 *
766 700 * Sets up the required event cube and time interval storage objects
... ... @@ -884,108 +818,197 @@ void ctfindvar::fill_cube(GCTAObservation* obs)
884 818  
885 819  
886 820 /***********************************************************************//**
887   - * @brief Initialize counts cube
  821 + * @brief Save the current sky direction to model if significance > thr
888 822 *
889   - * Sets up the required event cube and time interval storage objects
890   - * associated with this object
  823 + * @param[in] dir_pix pixel sky direction
  824 + * @param[out] m_model_sig_above_thr model position with sig>thr
891 825 ***************************************************************************/
892   -void ctfindvar::init_cube(void)
  826 +void ctfindvar::fill_model_sig_pos(const GSkyDir& dir_pix)
893 827 {
894   - // Create the basic skymap
895   - m_counts = create_map(m_obs);
  828 + //Set spatial components
  829 + GModelSpatialPointSource spatial(dir_pix);
896 830  
897   - // Resize to the appropriate number of time intervals
898   - m_counts.nmaps(m_gti.size());
  831 + //Set fake spectral component (powerlaw with 1% of Crab)
  832 + GModelSpectralPlaw spectral(5.7e-18, -2.48, GEnergy(0.3, "TeV"));
899 833  
900   - // Create the peaksigmap
901   - m_peaksigmap = m_counts.extract(1);
  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);
  846 +
902 847 }
903   -
904   -
905 848 /***********************************************************************//**
906   - * @brief Initialize counts cube
  849 + * @brief Get the map index associated with a given time
907 850 *
908   - * Sets up the required event cube and time interval storage objects
909   - * associated with this object
  851 + * @param[in] pix_number Spatial pixel
  852 + * @param[out] alpha_vector Vector containing effective exposure
910 853 ***************************************************************************/
911   -void ctfindvar::init_gtis(void)
  854 +void ctfindvar::fill_alpha_vector(const int& pix_number,
  855 + std::vector<double>& alpha_vector)
912 856 {
913   - double tinterval = (*this)["tinterval"].real();
914   - m_tstart = GTime("JD 999999999"); // large unreasonable Julian date
915   - m_tstop = GTime("JD 0"); // small unreasonable Julian date
  857 + // Pixel position
  858 + GSkyDir pix_dir = m_counts.inx2dir(pix_number);
916 859  
917   - // Set the start time
918   - if ((*this)["tmin"].is_valid()) {
919   - log_string(NORMAL, "Setting start time from command line");
920   - m_tstart = (*this)["tmin"].time();
921   - } else {
922   - log_string(NORMAL, "Setting start time from observations");
923   - for (int o=0; o<m_obs.size(); o++) {
924   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[o]);
925   -
926   - // Check that the start time of the run is less than tstart
927   - if (obs->gti().tstart() < m_tstart) {
928   - m_tstart = obs->gti().tstart();
929   - }
  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();
  866 +
  867 + if (roi_centre.dist_deg(pix_dir) > roi.radius()) {
  868 + continue;
930 869 }
931   - }
  870 +
  871 + // Convert sky direction to instrument direction
  872 + GCTAInstDir instdir = obs->pointing().instdir(pix_dir);
932 873  
933   - // Set the stop time
934   - if ((*this)["tmax"].is_valid()) {
935   - log_string(NORMAL, "Setting stop time from command line");
936   - m_tstop = (*this)["tmax"].time();
937   - } else {
  874 + // Extract the good time intervals of the observation
  875 + const GGti& obs_gti(obs->gti());
938 876  
939   - log_string(NORMAL, "Setting stop time from observations");
940   - for (int o=0; o<m_obs.size(); o++) {
941   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[o]);
942   -
943   - // Check that the stop time of the run is larger than tstop
944   - if (obs->gti().tstop() > m_tstop) {
945   - m_tstop = obs->gti().tstop();
  877 + // Loop over all time bins
  878 + for (int j = 0; j < m_gti.size(); j++) {
  879 +
  880 + // Make sure observation overlaps with this time interval
  881 + double exposure = gti_overlap(m_gti[j], obs_gti);
  882 + if (exposure > 0.0) {
  883 +
  884 + // Get IRF response
  885 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>
  886 + (obs->response());
  887 +
  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 + }
  897 +
  898 + // Get IRF background template
  899 + const GCTABackground* bkg = rsp->background();
  900 +
  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 + }
  910 +
  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);
946 914 }
947 915 }
948 916 }
949 917  
950   - // Fill the number of good time intervals
951   - double tstart_sec = m_tstart.secs();
952   - double tstop_sec = m_tstop.secs();
953   - int bins = (tstop_sec - tstart_sec) / tinterval + 0.5;
  918 + return;
  919 +}
954 920  
955   - // Log information
956   - log_value(NORMAL, "treference (mjd)", m_tstart.reference().mjdref());
957   - log_value(NORMAL, "tstart (sec)", tstart_sec);
958   - log_value(NORMAL, "tstop (sec)", tstop_sec);
959   - log_value(NORMAL, "Total time", tstop_sec-tstart_sec);
960   - log_value(NORMAL, "Time interval", tinterval);
961   - log_value(NORMAL, "nbins", bins);
962 921  
963   - // Throw if tstart == tstop
964   - if (m_tstart == m_tstop) {
965   - throw GException::invalid_value(G_INIT_GTIS,
966   - "Start time is equal to stop time");
967   - } else if (bins < 2) {
968   - std::string msg = "Method requires at least two time bins (" +
969   - gammalib::str(bins) + " bins found).";
970   - throw GException::invalid_value(G_INIT_GTIS, msg);
971   - }
972 922  
973   - for (int i=0; i<bins; i++) {
974   -
975   - // Get the next time interval
976   - GTime next_tstart(m_tstart + i*tinterval);
977   - GTime next_tstop(m_tstart + (i+1.0)*tinterval);
978   - GGti next_gti(next_tstart, next_tstop);
979 923  
980   - // Make sure there's an observation that overlaps with this GTI
981   - for (int o=0; o<m_obs.size(); o++) {
982   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[o]);
983   -
984   - if (gti_overlap(obs->gti(), next_gti) > 0.0) {
985   - m_gti.push_back(next_gti);
986   - break;
  924 +/***********************************************************************//**
  925 + * @brief Obtain the significance of variability in a given bin
  926 + *
  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
  930 + *
  931 + * Significance is computed according to Li & Ma equation 17:
  932 + * see: https://doi.org/10.1086/161295
  933 + * Some modification is made in order to handle the case where Non = 0,
  934 + * but Noff != 0.
  935 + ***************************************************************************/
  936 +void ctfindvar::get_variability_sig(const int& pix_number,
  937 + const int& nbins,
  938 + GNdarray& sig_histogram)
  939 +{
  940 + std::vector<bool> accepted_bin_bckg_vector(nbins, true);
  941 + bool background_validated = false;
  942 + double non, noff;
  943 + double alpha, sig;
  944 +
  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);
  948 +
  949 + // Loop over pixels until all background pixls are removed
  950 + while (background_validated == false) {
  951 + background_validated = true;
  952 +
  953 + // Loop over all the GTIs of the pixel
  954 + for (int i = 0; i < nbins; i++) {
  955 +
  956 + // The GTI is discared from bckg calculation and not checked again.
  957 + if (!accepted_bin_bckg_vector[i]) {
  958 + 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) {
  962 + accepted_bin_bckg_vector[i] = false;
  963 + continue;
  964 + }
  965 +
  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++) {
  970 + if (j != i && accepted_bin_bckg_vector[j] == 1) {
  971 + noff += m_counts(pix_number, j);
  972 + alpha += alpha_vector[j];
  973 + }
  974 + }
  975 +
  976 + // The background is averaged on the number of bins -1
  977 + non = m_counts(pix_number, i);
  978 + alpha = alpha_vector[i]/alpha;
  979 +
  980 + // Compute sensitivity in Gaussian sigma using Li & Ma
  981 + double alpha1 = alpha + 1.0;
  982 + double ntotal = non+noff;
  983 + double arg1 = non/ntotal;
  984 + 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;
  1000 + }
  1001 +
  1002 + // Update the significance
  1003 + sig_histogram(i) = sig;
  1004 +
  1005 + // If bin is significant, it's removed from bckg and we loop again.
  1006 + if (sig > m_sig_threshold) {
  1007 + accepted_bin_bckg_vector[i] = false;
  1008 + background_validated = false;
987 1009 }
988 1010 }
  1011 +
989 1012 }
990 1013 }
991 1014  
... ...
src/ctfindvar/ctfindvar.hpp
  View file @90c2e5d
1 1 /***************************************************************************
2   - * ctfindvar - search time variability tool *
  2 + * ctfindvar - Time variability search tool *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2018 by Simon Bonnefoy *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctfindvar.hpp
23   - * @brief search time variability tool definition
  23 + * @brief Time variability search tool definition
24 24 * @author Simon Bonnefoy
25 25 */
26 26  
... ... @@ -40,7 +40,7 @@
40 40 /***********************************************************************//**
41 41 * @class ctfindvar
42 42 *
43   - * @brief search time variability tool
  43 + * @brief Time variability search tool
44 44 *
45 45 * @todo Add tool description.
46 46 ***************************************************************************/
... ... @@ -68,17 +68,17 @@ public:
68 68  
69 69 protected:
70 70 // Protected methods
71   - void init_cube(void);
72   - void fill_model_sig_pos(const GSkyDir& dir_pix);
73   - void init_gtis(void);
74 71 void init_members(void);
75 72 void copy_members(const ctfindvar& app);
76 73 void free_members(void);
  74 + void get_parameters(void);
  75 + void init_cube(void);
  76 + void init_gtis(void);
77 77 void fill_cube(GCTAObservation* obs);
  78 + void fill_model_sig_pos(const GSkyDir& dir_pix);
78 79 void fill_alpha_vector(const int& pix_number,
79 80 std::vector<double>& alpha_vector);
80   - void get_parameters(void);
81   - void get_variability_sig(const int& pix_number,
  81 + void get_variability_sig(const int& pix_number,
82 82 const int& nbins,
83 83 GNdarray& sig_histogram);
84 84 double gti_overlap(const GGti& gti1,
... ...
src/ctfindvar/ctfindvar.par
  View file @90c2e5d
... ... @@ -22,6 +22,10 @@ outcube, f, h, NONE,,, &quot;Filename for counts cube&quot;
22 22 caldb, s, a, prod2,,, "Calibration database"
23 23 irf, s, a, South_0.5h,,, "Instrument response function"
24 24 prefix, s, a, findvar_,,, "Prefix for output files"
  25 +
  26 +#
  27 +# Variability search parameters
  28 +#==============================
25 29 minoff, r, h, 0,,, "Minimum required events in off bin for computing Noff and alpha"
26 30 emin, r, a, 0,,, "Lower energy limit (units of TeV; 0 => no cut)"
27 31 emax, r, a, 0,,, "Upper energy limit (units of TeV; 0 => no cut)"
... ...
src/ctfindvar/main.cpp
  View file @90c2e5d
1 1 /***************************************************************************
2   - *ctfindvar - search time variability tool main code *
  2 + * ctfindvar - Time variability search tool main code *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2018 by Simon Bonnefoy *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctfindvar/main.cpp
23   - * @brief search time variability tool main code
  23 + * @brief Time variability search tool main code
24 24 * @author Simon Bonnefoy
25 25 */
26 26  
... ...