Line data    Source code

       1             : //# SimpleSubMS.cc 
       2             : //# Copyright (C) 2010
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify
       6             : //# it under the terms of the GNU General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or
       8             : //# (at your option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful,
      11             : //# but WITHOUT ANY WARRANTY; without even the implied warranty of
      12             : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      13             : //# GNU General Public License for more details.
      14             : //# 
      15             : //# You should have received a copy of the GNU General Public License
      16             : //# along with this library; if not, write to the Free Software
      17             : //# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id: $
      27             : #include <casacore/casa/System/AppInfo.h>
      28             : #include <casacore/casa/Utilities/GenSort.h>
      29             : #include <msvis/MSVis/VisBuffer.h>
      30             : #include <msvis/MSVis/VisibilityIterator.h>
      31             : #include <msvis/MSVis/SimpleSubMS.h>
      32             : 
      33             : using namespace casacore;
      34             : namespace casa {
      35             : 
      36           0 :   SimpleSubMS::SimpleSubMS(MeasurementSet& ms) : SubMS(ms){
      37             : 
      38           0 :   }
      39             : 
      40           0 :   SimpleSubMS::~SimpleSubMS(){
      41             : 
      42           0 :   }
      43             : 
      44           0 :   MeasurementSet* SimpleSubMS::makeMemSubMS(const MS::PredefinedColumns& whichcol, const String& name){
      45           0 :     LogIO os(LogOrigin("SimpleSubMS", "makeMemSubMS()"));
      46             :     
      47           0 :     if(max(fieldid_p) >= Int(ms_p.field().nrow())){
      48             :       os << LogIO::SEVERE 
      49             :          << "Field selection contains elements that do not exist in "
      50             :          << "this MS"
      51           0 :          << LogIO::POST;
      52           0 :       ms_p=MeasurementSet();
      53           0 :       return 0;  
      54             :     }
      55           0 :     if(max(spw_p) >= Int(ms_p.spectralWindow().nrow())){
      56             :       os << LogIO::SEVERE 
      57             :          << "SpectralWindow selection contains elements that do not exist in "
      58             :          << "this MS"
      59           0 :          << LogIO::POST;
      60           0 :       ms_p=MeasurementSet();
      61           0 :       return 0;
      62             :     }
      63             :     
      64           0 :     if(!makeSelection()){
      65             :       os << LogIO::SEVERE 
      66             :          << "Failed on selection: combination of spw and/or field and/or time "
      67             :          << "chosen may be invalid."
      68           0 :          << LogIO::POST;
      69           0 :       ms_p=MeasurementSet();
      70           0 :       return 0;
      71             :     }
      72             :     //Make sure the damn ms is sorted the way we want ...to copy the meta-columns as is.
      73           0 :     Block<String> sortcol(4);
      74           0 :     sortcol[0]=MS::columnName(MS::ARRAY_ID);
      75           0 :     sortcol[1]=MS::columnName(MS::FIELD_ID);
      76           0 :     sortcol[2]=MS::columnName(MS::DATA_DESC_ID);
      77           0 :     sortcol[3]=MS::columnName(MS::TIME);
      78           0 :     MeasurementSet msselsort(mssel_p.sort(sortcol, Sort::Ascending, Sort::QuickSort));
      79             : 
      80           0 :     mscIn_p=new MSColumns(msselsort);
      81           0 :     String msname=name;
      82           0 :     if(msname==""){
      83           0 :       msname=AppInfo::workFileName(100, "TempSubMS");
      84             :     }
      85           0 :     MeasurementSet* ah=setupMS(msname, nchan_p[0], ncorr_p[0],  
      86           0 :                                mscIn_p->observation().telescopeName()(0),
      87           0 :                                Vector<MS::PredefinedColumns>(1,whichcol));
      88           0 :     if(!ah)
      89           0 :       return 0;
      90             :     
      91             :     MeasurementSet * outpointer;
      92           0 :     if(name==""){
      93           0 :       ah->markForDelete();
      94           0 :       outpointer= new MeasurementSet(ah->copyToMemoryTable("TmpMemoryMS"));
      95           0 :       delete ah;
      96             :     }
      97             :     else{
      98           0 :       outpointer=ah;
      99             :     }
     100           0 :     if(!outpointer)
     101           0 :       return 0;
     102           0 :     outpointer->initRefs();
     103             :     
     104           0 :     msOut_p = *outpointer;
     105           0 :     msc_p = new MSColumns(msOut_p);
     106             :     
     107             :     //Do the subtables
     108           0 :     msc_p->setEpochRef(MEpoch::castType(mscIn_p->timeMeas().getMeasRef().getType()));
     109             :     
     110             :     // UVW is the only other Measures column in the main table.
     111           0 :     msc_p->uvwMeas().setDescRefCode(Muvw::castType(mscIn_p->uvwMeas().getMeasRef().getType()));
     112             :     
     113             :     // fill or update
     114           0 :     if(!fillDDTables()){
     115           0 :       return 0;
     116             :     }
     117             : 
     118             :     // SourceIDs need to be remapped around here.  It could not be done in
     119             :     // selectSource() because mssel_p was not setup yet.
     120           0 :     relabelSources();
     121             :     
     122           0 :     fillFieldTable();
     123           0 :     copySource();
     124             :     
     125           0 :     copyAntenna();
     126           0 :     if(!copyFeed())         // Feed table writing has to be after antenna 
     127           0 :       return 0;
     128             :     
     129           0 :     copyObservation();
     130           0 :     copyPointing();
     131             :     // copyWeather();
     132             :     // copySyscal();   // ?
     133             :     /////////////Done with the sub tables...now to the main
     134           0 :     fillAccessoryMainCols();
     135           0 :     msc_p->weight().putColumn(mscIn_p->weight());
     136           0 :     msc_p->sigma().putColumn(mscIn_p->sigma());
     137           0 :     Block<Int> sort(4);
     138           0 :     sort[0]=MS::ARRAY_ID;
     139           0 :     sort[1]=MS::FIELD_ID;
     140           0 :     sort[2]=MS::DATA_DESC_ID;
     141           0 :     sort[3]=MS::TIME;
     142           0 :     ROVisibilityIterator vi(msselsort, sort);
     143           0 :     for (Int k=0; k < spw_p.shape()(0) ; ++k){ 
     144             :        os << LogIO::NORMAL
     145           0 :           << "Selecting "<< nchan_p[k] << " channels, starting at " <<chanStart_p[k] << " for spw " << spw_p[k] << LogIO::POST; ;
     146             :        
     147           0 :        vi.selectChannel(1, chanStart_p[k], nchan_p[k],
     148           0 :                         chanStep_p[k], spw_p[k]);
     149             :     }
     150             :     /////slurp test
     151             :     //vi.slurp();
     152             :     ///
     153           0 :     Bool doSpWeight = !(mscIn_p->weightSpectrum().isNull()) &&
     154           0 :                       mscIn_p->weightSpectrum().isDefined(0);
     155           0 :     Int rowsdone=0;
     156           0 :     Int rowsnow=0;
     157           0 :     VisBuffer vb(vi);
     158           0 :     Vector<Int> spwindex(max(spw_p)+1);
     159           0 :     spwindex.set(-1);
     160           0 :     for(uInt k = 0; k < spw_p.nelements(); ++k)
     161           0 :       spwindex[spw_p[k]] = k;
     162           0 :     for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
     163           0 :       for (vi.origin(); vi.more(); vi++) {
     164           0 :         Int spw=spwindex[vb.spectralWindow()];
     165           0 :         if(spw <0){
     166           0 :           cerr << "Huh ?: The programmer calling this code is useless.." << endl;
     167           0 :           return 0;
     168             :         }
     169           0 :         rowsnow=vb.nRow();      
     170           0 :         RefRows rowstoadd(rowsdone, rowsdone+rowsnow-1);
     171           0 :         if(whichcol==MS::DATA)
     172           0 :           msc_p->data().putColumnCells(rowstoadd, vb.visCube());
     173           0 :         else if(whichcol==MS::MODEL_DATA)
     174           0 :           msc_p->data().putColumnCells(rowstoadd, vb.modelVisCube());
     175           0 :         else if(whichcol==MS::CORRECTED_DATA)
     176           0 :           msc_p->data().putColumnCells(rowstoadd, vb.correctedVisCube());
     177             :         else
     178           0 :           cerr << "Column not requested not yet supported to be loaded to memory" << endl;
     179           0 :         msc_p->flag().putColumnCells(rowstoadd, vb.flagCube());
     180           0 :         if(doSpWeight)
     181           0 :           msc_p->weightSpectrum().putColumnCells(rowstoadd, vb.weightSpectrum());
     182           0 :         rowsdone += rowsnow;
     183             :       }
     184             :     }
     185           0 :     return outpointer;
     186             :   }
     187             : 
     188             : 
     189             : } //#End casa namespace