Eirik Hatlen / ctools

Commit 42ccfdb30c4258e667436f0ae3fd3a369c6e15f7

Authored by Josh Cardenzana
Committed by Jürgen Knödlseder
1 parent bf9f1f3e

Add output for significance distribution of pixel with maximum significance (#2713)

Inline Side-by-side
Showing 2 changed files with 67 additions and 27 deletions
src/ctfindvar/ctfindvar.cpp
  View file @42ccfdb
... ... @@ -303,11 +303,10 @@ void ctfindvar::run(void)
303 303  
304 304 //preparing the histogram with significance evolution for the source center and the highest sig-pix
305 305 m_pixsigsrc = GNdarray(srcInxPix.size(), nbins);
306   - m_pixsigmax = GNdarray(1,nbins);
  306 + m_pixsigmax = GNdarray(nbins);
307 307  
308 308 //Prepare the final skymap with the max significance of each pixel
309 309 double max_sig=0;
310   - GSkyDir max_sig_dir;
311 310  
312 311 //looping over all the pixels in the cube
313 312 #pragma omp parallel for
... ... @@ -365,15 +364,10 @@ void ctfindvar::run(void)
365 364 {
366 365 if (max(pixSig) > max_sig)
367 366 {
368   - max_sig = max(pixSig);
369   - // Store information for pixel with the maximum significance
370   - for (int i=0; i<nbins; i++) {
371   - m_pixsigmax(0,i) = pixSig(i);
372   - }
373   -
374 367 // Update the pixel position with the maximum sigma
375   - max_sig = max(pixSig);
376   - max_sig_dir = m_counts.inx2dir(pix_number);
  368 + max_sig = max(pixSig);
  369 + m_max_sig_dir = m_counts.inx2dir(pix_number);
  370 + m_pixsigmax = pixSig;
377 371 }
378 372 }
379 373 //storing sig in skymap
... ... @@ -382,8 +376,8 @@ void ctfindvar::run(void)
382 376  
383 377 log_header1(NORMAL, "Analysis finished");
384 378 log_value(NORMAL, "Maximum sigma", max_sig);
385   - log_value(NORMAL, "Max pixel RA", max_sig_dir.ra_deg());
386   - log_value(NORMAL, "Max pixel DEC", max_sig_dir.dec_deg());
  379 + log_value(NORMAL, "Max pixel RA", m_max_sig_dir.ra_deg());
  380 + log_value(NORMAL, "Max pixel DEC", m_max_sig_dir.dec_deg());
387 381  
388 382 return;
389 383 }
... ... @@ -925,7 +919,8 @@ void ctfindvar::write_srchist(void)
925 919 double tinterval = (*this)["tinterval"].real();
926 920 int nentries = (m_tstop-m_tstart) / tinterval + 0.5;
927 921 int nsources = m_pixsigsrc.shape()[0];
928   - GNdarray time_info(2,nentries);
  922 + GNdarray time_info(2, nentries);
  923 + GNdarray max_pixel_info(nentries);
929 924 GNdarray src_info(nsources, nentries);
930 925  
931 926 // Store the full set of information
... ... @@ -949,13 +944,18 @@ void ctfindvar::write_srchist(void)
949 944 src_info(src,i) = 0.0;
950 945 }
951 946 }
  947 +
  948 + // Fill max pixel information
  949 + if (index >= 0) {
  950 + max_pixel_info(i) = m_pixsigmax(index);
  951 + }
952 952 }
953 953  
954 954 // Call the appropriate method for writing the data
955 955 if ((*this)["histtype"].string() == "CSV") {
956   - write_srchist_csv(time_info, src_info);
  956 + write_srchist_csv(time_info, max_pixel_info, src_info);
957 957 } else {
958   - write_srchist_fits(time_info, src_info);
  958 + write_srchist_fits(time_info, max_pixel_info, src_info);
959 959 }
960 960 }
961 961  
... ... @@ -963,6 +963,7 @@ void ctfindvar::write_srchist(void)
963 963 * @brief Write individual histograms in CSV format
964 964 ***************************************************************************/
965 965 void ctfindvar::write_srchist_csv(const GNdarray& time_info,
  966 + const GNdarray& max_pixel_info,
966 967 const GNdarray& src_info)
967 968 {
968 969 log_string(NORMAL,"Storing histograms in CSV format is not currently supported.");
... ... @@ -981,47 +982,83 @@ void ctfindvar::write_srchist_csv(const GNdarray&amp; time_info,
981 982 * significance of the source in that time bin.
982 983 ***************************************************************************/
983 984 void ctfindvar::write_srchist_fits(const GNdarray& time_info,
  985 + const GNdarray& max_pixel_info,
984 986 const GNdarray& src_info)
985 987 {
986 988 log_string(NORMAL, "Storing histograms in FITS format");
987 989  
988 990 // Create the table
989 991 int nrows = time_info.shape()[1];
990   - GFitsBinTable table(nrows);
991   - table.extname("SRCSIGNIF");
  992 + GFitsBinTable table_signif(nrows);
  993 + GFitsBinTable table_dir(2);
  994 + table_signif.extname("SRCSIGNIF");
  995 + table_dir.extname("SRCPOSITION");
992 996  
993   - // Append the start and stop time columns
  997 + // Append the start/stop time and max pixel columns
994 998 GFitsTableFloatCol start_mjd("TSTART", nrows);
995 999 GFitsTableFloatCol stop_mjd("TSTOP", nrows);
996   -
  1000 + GFitsTableFloatCol max_pix("MAXSIGPIXEL", nrows);
997 1001 for (int i=0; i<nrows; i++) {
998 1002 start_mjd(i) = time_info(0,i);
999 1003 stop_mjd(i) = time_info(1,i);
  1004 + max_pix(i) = max_pixel_info(i);
1000 1005 }
1001   - table.append(start_mjd);
1002   - table.append(stop_mjd);
  1006 + table_signif.append(start_mjd);
  1007 + table_signif.append(stop_mjd);
  1008 + table_signif.append(max_pix);
  1009 +
  1010 + // Store the position of the peak pixel
  1011 + GFitsTableFloatCol max_pix_dir("MAXSIGPIXEL",2);
  1012 + max_pix_dir(0) = m_max_sig_dir.ra_deg();
  1013 + max_pix_dir(1) = m_max_sig_dir.dec_deg();
  1014 + table_dir.append(max_pix_dir);
1003 1015  
1004 1016 // Append the distribution for each source
1005 1017 int src_count = src_info.shape()[0];
1006 1018 for (int src=0; src<src_count; src++) {
  1019 +
1007 1020 std::string src_name("SOURCE");
  1021 + GSkyDir src_dir;
  1022 +
  1023 + // Get the name from the input model file
1008 1024 if (m_inmodel.size() > 0) {
1009 1025 src_name = m_inmodel[src]->name();
  1026 +
  1027 + // Get source position
  1028 + GModelSky* model = dynamic_cast<GModelSky*>(m_inmodel[src]);
  1029 + GModelSpatial* model_spatial = model->spatial();
  1030 + GSkyRegion* model_region = model_spatial->region();
  1031 + src_dir = dynamic_cast<GSkyRegionCircle*>(model_region)->centre();
  1032 + }
  1033 + // ... otherwise set the position from xsrc,ysrc
  1034 + else {
  1035 + if ((*this)["coordsys"].string() == "CEL") {
  1036 + src_dir.radec_deg((*this)["xsrc"].real(), (*this)["yref"].real());
  1037 + } else {
  1038 + src_dir.lb_deg((*this)["xsrc"].real(), (*this)["yref"].real());
  1039 + }
1010 1040 }
1011 1041  
1012   - GFitsTableFloatCol src_sig(src_name, nrows);
1013   -
  1042 + // Create and fill a new column
  1043 + GFitsTableFloatCol src_sig_col(src_name, nrows);
1014 1044 for (int i=0; i<nrows; i++) {
1015   - src_sig(i) = src_info(src, i);
  1045 + src_sig_col(i) = src_info(src, i);
1016 1046 }
1017   - table.append(src_sig);
  1047 + table_signif.append(src_sig_col);
  1048 +
  1049 + // Store the positions of the source
  1050 + GFitsTableFloatCol src_dir_col(src_name,2);
  1051 + src_dir_col(0) = src_dir.ra_deg();
  1052 + src_dir_col(1) = src_dir.dec_deg();
  1053 + table_dir.append(src_dir_col);
1018 1054 }
1019 1055  
1020 1056 // Store information regarding the source positions
1021 1057  
1022   - // Write the table to a file
  1058 + // Write the tables to a file
1023 1059 GFilename outfile((*this)["prefix"].string() + "srcsig.fits");
1024 1060 GFits fitsfile;
1025   - fitsfile.append(table);
  1061 + fitsfile.append(table_signif);
  1062 + fitsfile.append(table_dir);
1026 1063 fitsfile.saveto(outfile, (*this)["clobber"].boolean() );
1027 1064 }
... ...
src/ctfindvar/ctfindvar.hpp
  View file @42ccfdb
... ... @@ -81,14 +81,17 @@ protected:
81 81 // Writing methods for individual source histograms
82 82 void write_srchist(void);
83 83 void write_srchist_csv(const GNdarray& time_info,
  84 + const GNdarray& max_pixel_info,
84 85 const GNdarray& src_info);
85 86 void write_srchist_fits(const GNdarray& time_info,
  87 + const GNdarray& max_pixel_info,
86 88 const GNdarray& src_info);
87 89  
88 90 // Protected members
89 91 GSkyMap m_counts; //!< Counts for each time interval
90 92 GGti m_gti; //!< List of time intervals
91 93 GModels m_inmodel; //!< List of models for source positions
  94 + GSkyDir m_max_sig_dir; //!< Sky direction associated with maximum significance
92 95 double m_minoff; //!< Minimum counts for use in significance calculation
93 96 GSkyMap m_peaksigmap; //!< Skymap holding the maximum significance
94 97 GNdarray m_pixsigsrc; //!< Store distributions of the source significances
... ...