Line data    Source code

       1             : /**
       2             :    Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
       3             :    Initial version January 2010.
       4             :    Maintained by ESO since 2013. 
       5             :    
       6             :    This file is part of LibAIR and is licensed under GNU Public
       7             :    License Version 2
       8             :    
       9             :    \file msutils.cpp
      10             :    Renamed msutils.cc 2023
      11             : 
      12             : */
      13             : 
      14             : #include <stdexcept>
      15             : #include "msutils.h"
      16             : 
      17             : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
      18             : #include <casacore/tables/Tables/Table.h>
      19             : #include <casacore/ms/MeasurementSets/MSFieldColumns.h>
      20             : #include <casacore/ms/MeasurementSets/MSField.h>
      21             : #include <casacore/ms/MeasurementSets/MSColumns.h>
      22             : 
      23             : namespace LibAIR2 {
      24             : 
      25             :   /** Channel frequencies for spectral window spw     
      26             :    */
      27           0 :   void spwChannelFreq(const casacore::MeasurementSet &ms,
      28             :                       size_t spw,
      29             :                       std::vector<double> &fres)
      30             :   {
      31           0 :     casacore::MSSpectralWindow specTable(ms.spectralWindow());
      32             : 
      33             :     // Number of channels 
      34             :     casacore::ScalarColumn<casacore::Int> nc
      35             :       (specTable,
      36           0 :        casacore::MSSpectralWindow::columnName(casacore::MSSpectralWindow::NUM_CHAN));
      37             : 
      38             :     // Frequencies of the channels
      39             :     casacore::ArrayColumn<casacore::Double>
      40             :       chfreq(specTable,
      41           0 :              casacore::MSSpectralWindow::columnName(casacore::MSSpectralWindow::CHAN_FREQ));
      42             : 
      43           0 :     fres.resize(nc(spw));
      44           0 :     casacore::Array<casacore::Double> freq;
      45           0 :     chfreq.get(spw, freq, casacore::True);
      46             :     
      47           0 :     for(size_t i=0; i< static_cast<size_t>(nc(spw)); ++i)
      48             :     {
      49           0 :       fres[i]=freq(casacore::IPosition(1,i));
      50             :     }
      51           0 :   }
      52             : 
      53           0 :   void fieldIDs(const casacore::MeasurementSet &ms,
      54             :                 std::vector<double> &time,
      55             :                 std::vector<int> &fieldID,
      56             :                 std::vector<int> &sourceID,
      57             :                 const std::vector<size_t> &sortedI)
      58             :   {
      59           0 :     const casacore::MSMainColumns cols(ms);
      60           0 :     const casacore::ScalarColumn<casacore::Int> &f= cols.fieldId();
      61           0 :     const casacore::ScalarColumn<casacore::Double> &t= cols.time();
      62             : 
      63           0 :     std::map<size_t, size_t> srcmap=getFieldSrcMap(ms);
      64             : 
      65           0 :     const size_t nrows=f.nrow();    
      66             : 
      67           0 :     if(sortedI.size()!= nrows){
      68           0 :       throw std::runtime_error("Time-sorted row vector must have same size as MS.");
      69             :     }
      70             : 
      71           0 :     time.resize(nrows);
      72           0 :     fieldID.resize(nrows);
      73           0 :     sourceID.resize(nrows);
      74           0 :     for(size_t ii=0; ii<nrows; ++ii)
      75             :     {
      76           0 :       size_t i = sortedI[ii];
      77             : 
      78           0 :       time[ii]=t(i);
      79           0 :       fieldID[ii]=f(i);
      80           0 :       if (srcmap.count(f(i)) == 0)
      81             :       {
      82           0 :         throw std::runtime_error("Encountered data without associated source");
      83             :       }
      84             :       /**checking if compiler version has map::at*/
      85             : #if __GNUC__ <= 4 and __GNUC_MINOR__ < 1
      86             :       sourceID[ii]=srcmap[(f(i))];
      87             : #else
      88           0 :       sourceID[ii]=srcmap.at(f(i));
      89             : #endif
      90             : 
      91             :     }
      92           0 :   }
      93             : 
      94             : 
      95           0 :   std::ostream & operator<<(std::ostream &o, 
      96             :                             const StateIntentMap &t)
      97             :   {
      98           0 :     o<<"     Map between State_ID and Scan Intents "<<std::endl
      99           0 :      <<"----------------------------------------------"<<std::endl;
     100           0 :     for(StateIntentMap::const_iterator i=t.begin(); 
     101           0 :         i!= t.end(); 
     102           0 :         ++i)
     103             :     {
     104           0 :       o<<i->first<<": "<<i->second
     105           0 :        <<std::endl;
     106             :     }
     107           0 :     return o;
     108             :   }
     109             : 
     110             : 
     111           0 :   void scanIntents(const casacore::MeasurementSet &ms,
     112             :                    StateIntentMap &mi)
     113             :   {
     114           0 :     mi.clear();
     115           0 :     casacore::MSState state(ms.state());
     116             :     const casacore::ScalarColumn<casacore::String> modes(state,
     117           0 :                                                    casacore::MSState::columnName(casacore::MSState::OBS_MODE));
     118           0 :     for(size_t i=0; i<state.nrow(); ++i)
     119             :     {
     120           0 :       casacore::String t;
     121           0 :       modes.get(i,t);
     122           0 :       mi.insert(std::pair<size_t, std::string>(i, t));
     123           0 :     }
     124           0 :   }
     125             : 
     126           0 :   std::set<size_t> skyStateIDs(const casacore::MeasurementSet &ms)
     127             :   {
     128           0 :     std::set<size_t> res;
     129           0 :     StateIntentMap mi;
     130           0 :     scanIntents(ms, mi);
     131           0 :     for(StateIntentMap::const_iterator i=mi.begin();
     132           0 :         i!=mi.end();
     133           0 :         ++i)
     134             :     {
     135           0 :       if((i->second.find("ON_SOURCE") != i->second.npos
     136           0 :           || i->second.find("REFERENCE") != i->second.npos
     137           0 :           || i->second.find("SIGNAL") != i->second.npos)
     138           0 :          && (i->second.find("CALIBRATE_ATMOSPHERE") == i->second.npos))
     139             :         {
     140           0 :           res.insert(i->first);
     141             :         }
     142             :     }
     143           0 :     return res;
     144           0 :   }
     145             : 
     146           0 :   field_t getFieldNames(const casacore::MeasurementSet &ms)
     147             :   {
     148           0 :     field_t res;
     149           0 :     const casacore::MSField & fieldT(ms.field());
     150           0 :     const size_t nfields=fieldT.nrow();
     151             : 
     152           0 :     casacore::MSFieldColumns fcols(fieldT);
     153           0 :     const casacore::ScalarColumn<casacore::String> &names (fcols.name());
     154           0 :     for(size_t i=0; i<nfields; ++i)
     155             :     {
     156           0 :       res.insert(field_t::value_type(i, std::string(names(i))));
     157             :     }
     158           0 :     return res;
     159           0 :   }
     160             : 
     161           0 :   field_t getSourceNames(const casacore::MeasurementSet &ms)
     162             :   {
     163           0 :     field_t res;
     164           0 :     const casacore::MSSource & srcTab(ms.source());
     165           0 :     const size_t nsource=srcTab.nrow();
     166             : 
     167           0 :     casacore::MSSourceColumns scols(srcTab);
     168           0 :     const casacore::ScalarColumn<casacore::String> &names (scols.name());
     169           0 :     const casacore::ScalarColumn<casacore::Int> & ids (scols.sourceId());
     170           0 :     for(size_t i=0; i<nsource; ++i)
     171             :     {
     172           0 :       res.insert(field_t::value_type(ids(i), std::string(names(i))));
     173             :     }
     174           0 :     return res;
     175             :     
     176           0 :   }
     177             : 
     178           0 :   std::set<size_t> getSrcFields(const casacore::MeasurementSet &ms,
     179             :                                 const std::string &source)
     180             :   {
     181           0 :     std::set<size_t> res;
     182           0 :     const casacore::MSField & fieldT(ms.field());
     183           0 :     const size_t nfields=fieldT.nrow();
     184             : 
     185           0 :     field_t sources=getSourceNames(ms);
     186             : 
     187           0 :     casacore::MSFieldColumns fcols(fieldT);
     188           0 :     const casacore::ScalarColumn<casacore::Int> &fieldsrc (fcols.sourceId());
     189           0 :     for(size_t i=0; i<nfields; ++i)
     190             :     {
     191             :       try {
     192           0 :         if(sources.at(fieldsrc(i))==source)
     193           0 :           res.insert(i);
     194             :       }
     195           0 :       catch (const std::out_of_range &e)
     196             :       {
     197           0 :         std::cout<<"ERROR: Could not find source name for source ID " <<fieldsrc(i)<<std::endl
     198           0 :                  <<"       This could be related to CSV-910"<<std::endl
     199           0 :                  <<"       Continuing with processing but the statistics might be affected by this"<<std::endl
     200           0 :                  <<std::endl;
     201           0 :       }
     202             :     }
     203           0 :     return res;
     204           0 :   }
     205             : 
     206           0 :   std::map<size_t, size_t> getFieldSrcMap(const casacore::MeasurementSet &ms)
     207             :   {
     208           0 :     std::map<size_t, size_t> res;
     209           0 :     const casacore::MSField & fieldT(ms.field());
     210           0 :     const size_t nfields=fieldT.nrow();
     211           0 :     casacore::MSFieldColumns fcols(fieldT);
     212           0 :     const casacore::ScalarColumn<casacore::Int> &fieldsrc (fcols.sourceId());
     213           0 :     for(size_t i=0; i<nfields; ++i)
     214             :     {
     215           0 :       res.insert(std::pair<size_t, size_t>(i, fieldsrc(i)));
     216             :     }
     217           0 :     return res;
     218           0 :   }
     219             : }
     220             : 
     221             :