Line data    Source code

       1             : //# CalSet.cc: Implementation of Calibration parameter cache
       2             : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: casa-feedback@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             : 
      27             : #include <synthesis/CalTables/CTGlobals.h>
      28             : #include <synthesis/CalTables/NewCalTable.h>
      29             : #include <synthesis/CalTables/CTIter.h>
      30             : #include <synthesis/CalTables/RIorAParray.h>
      31             : 
      32             : #include <casacore/casa/Arrays.h>
      33             : #include <casacore/casa/BasicSL/String.h>
      34             : #include <casacore/casa/Utilities/GenSort.h>
      35             : #include <casacore/casa/Exceptions/Error.h>
      36             : #include <casacore/tables/Tables/Table.h>
      37             : #include <casacore/tables/Tables/TableIter.h>
      38             : #include <casacore/tables/Tables/TableVector.h>
      39             : 
      40             : #include <sstream>
      41             : 
      42             : #include <casacore/casa/Logging/LogMessage.h>
      43             : #include <casacore/casa/Logging/LogSink.h>
      44             : 
      45             : using namespace casacore;
      46             : namespace casa { //# NAMESPACE CASA - BEGIN
      47             : 
      48             : 
      49           0 : void smoothCT(NewCalTable ct,
      50             :               const String& smtype,
      51             :               const Double& smtime,
      52             :               Vector<Int> selfields) {
      53             : 
      54             :   // Complex parameters?
      55           0 :   Bool cmplx=ct.isComplex();
      56             : 
      57             :   // half-width
      58           0 :   Double thw(smtime/2.0);
      59             : 
      60             :   // Workspace
      61           0 :   Vector<Double> times;
      62           0 :   Vector<Float> p,newp;
      63           0 :   Vector<Bool> pOK, newpOK;
      64             : 
      65           0 :   Cube<Float> fpar;
      66           0 :   Cube<Bool> fparok,newfparok;
      67             : 
      68           0 :   Vector<Bool> mask;
      69             : 
      70           0 :   IPosition blc(3,0,0,0), fblc(3,0,0,0);
      71           0 :   IPosition trc(3,0,0,0), ftrc(3,0,0,0);
      72           0 :   IPosition vec(1,0);
      73             : 
      74           0 :   Block<String> cols(4);
      75           0 :   cols[0]="SPECTRAL_WINDOW_ID";
      76           0 :   cols[1]="FIELD_ID";
      77           0 :   cols[2]="ANTENNA1";
      78           0 :   cols[3]="ANTENNA2";
      79           0 :   CTIter ctiter(ct,cols);
      80             : 
      81           0 :   while (!ctiter.pastEnd()) {
      82             : 
      83           0 :     Int nSlot=ctiter.nrow();
      84           0 :     Int ifld=ctiter.thisField();
      85             : 
      86             :     // Only if more than one slot in this spw _AND_
      87             :     //  field is among those requested (if any)
      88           0 :     if (nSlot>1 && 
      89           0 :         (selfields.nelements()<1 || anyEQ(selfields,ifld))) {
      90             : 
      91             :       //UNUSED: Int ispw=ctiter.thisSpw();
      92             :       
      93           0 :       vec(0)=nSlot;
      94           0 :       trc(2)=ftrc(2)=nSlot-1;
      95             : 
      96           0 :       times.assign(ctiter.time());
      97             : 
      98             :       // Extract Float info
      99           0 :       if (cmplx)
     100           0 :         fpar.assign(ctiter.casfparam("AP"));
     101             :       else
     102           0 :         fpar.assign(ctiter.fparam());
     103             : 
     104           0 :       fparok.assign(!ctiter.flag());
     105           0 :       newfparok.assign(fparok);
     106           0 :       IPosition fsh(fpar.shape());
     107             : 
     108             :       // For each channel
     109           0 :       for (int ichan=0;ichan<fsh(1);++ichan) {
     110           0 :         blc(1)=trc(1)=fblc(1)=ftrc(1)=ichan;
     111             :         // For each param (pol)
     112           0 :         for (Int ipar=0;ipar<fsh(0);++ipar) {
     113           0 :           blc(0)=trc(0)=ipar;
     114           0 :           fblc(0)=ftrc(0)=ipar/(cmplx?2:1);
     115             :           
     116             :           // Reference slices of par/parOK
     117           0 :           p.reference(fpar(blc,trc).reform(vec));
     118           0 :           newp.assign(p);
     119           0 :           pOK.reference(fparok(fblc,ftrc).reform(vec));
     120           0 :           newpOK.reference(newfparok(fblc,ftrc).reform(vec));
     121             : 
     122             :        /*
     123             :             cout << ispw << " "
     124             :                  << ichan << " "
     125             :                  << ipar << " "
     126             :                  << "p.shape() = " << p.shape() << " "
     127             :                  << "pOK.shape() = " << pOK.shape() << " "
     128             :                  << endl;
     129             :        */
     130             : 
     131             :             
     132           0 :           Vector<Bool> mask;
     133           0 :           for (Int i=0;i<nSlot;++i) {
     134             :             // Make mask
     135           0 :             mask = pOK;
     136           0 :             mask = (mask && ( (times >  (times(i)-thw)) && 
     137           0 :                               (times <= (times(i)+thw)) ) );
     138             : 
     139             :             // Avoid explicit zeros, for now
     140             :             //      mask = (mask && amp>=FLT_MIN);
     141             : 
     142             : 
     143             :             //cout << "    " << ifld << " " << i << " " << idx(i) << " ";
     144             :             //for (Int j=0;j<mask.nelements();++j)
     145             :             //  cout << mask(j);
     146             :             //cout << endl;
     147             :             
     148           0 :             if (ntrue(mask)>0) {
     149           0 :               if (smtype=="mean") {
     150           0 :                 newp(i)=mean(p(mask));
     151             :               }
     152           0 :               else if (smtype=="median") {
     153           0 :                 newp(i)= median(p(mask),false);
     154             :               }
     155           0 :               newpOK(i)=true;
     156             :             }
     157             :             else 
     158           0 :               newpOK(i)=false;
     159             :             
     160             :           } // i
     161             :           // keep new ok info
     162           0 :           p=newp;
     163           0 :         } // ipar
     164             :       } // ichan
     165             : 
     166             :       // Put info back
     167           0 :       if (cmplx)
     168           0 :         ctiter.setcparam(RIorAPArray(fpar).c());
     169             :       else
     170           0 :         ctiter.setfparam(fpar);
     171             : 
     172           0 :       ctiter.setflag(!newfparok);
     173             : 
     174           0 :     } // nSlot>1
     175             : 
     176           0 :     ctiter.next();
     177             :   } // ispw
     178             :         
     179           0 : }
     180             : 
     181           0 : void assignCTScanField(NewCalTable& ct, String msName, 
     182             :                        Bool doField, Bool doScan, Bool doObs) {
     183             : 
     184             :   // TBD: verify msName is present and is an MS
     185             : 
     186             :   // Arrange to iterate only on SCAN  (and SPW?)
     187           0 :   Table mstab(msName,Table::Old);
     188             : 
     189             :   // How many scans in total?
     190           0 :   TableVector<Int> allscansTV(mstab,"SCAN_NUMBER");
     191           0 :   Vector<Int> allscans=allscansTV.makeVector();
     192           0 :   Int nScan=genSort(allscans,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
     193             : 
     194             :   //  cout << "Found " << nScan << " scans in " << msName << "." << endl;
     195             : 
     196             :   // Workspace
     197           0 :   Vector<Int> scanlist(nScan,-1);
     198           0 :   Vector<Double> timelo(nScan,DBL_MIN);
     199           0 :   Vector<Double> timehi(nScan,DBL_MAX);
     200           0 :   Vector<Int> fieldlist(nScan,-1);
     201           0 :   Vector<Int> obslist(nScan,-1);
     202           0 :   Vector<uInt> ord;
     203             :   {
     204           0 :     Block<String> cols(1);
     205           0 :     cols[0]="SCAN_NUMBER";
     206           0 :     TableIterator mstiter(mstab,cols);
     207             :     
     208             :     // Get time boundares and fields for each scan
     209           0 :     Int iscan(0);
     210           0 :     while (!mstiter.pastEnd()) {
     211           0 :       Table thistab(mstiter.table());
     212             :       
     213           0 :       Int scan=TableVector<Int>(thistab,"SCAN_NUMBER")(0);
     214           0 :       scanlist(iscan)=scan;
     215             :       
     216           0 :       fieldlist(iscan)=TableVector<Int>(thistab,"FIELD_ID")(0);
     217           0 :       obslist(iscan)=TableVector<Int>(thistab,"OBSERVATION_ID")(0);
     218             :       
     219           0 :       Vector<Double> times=TableVector<Double>(thistab,"TIME").makeVector();
     220           0 :       timelo(iscan)=min(times)-1e-5;
     221           0 :       timehi(iscan)=max(times)+1e-5;
     222             :       
     223           0 :       mstiter.next();
     224           0 :       ++iscan;
     225           0 :     }
     226             : 
     227             : 
     228             :     // Ensure time orderliness
     229           0 :     genSort(ord,timehi);
     230             : 
     231             :     /*
     232             :     cout << "scanlist=" << scanlist << endl;
     233             :     cout << "fieldlist=" << fieldlist << endl;
     234             :     cout << "timelo=" << timelo-timelo(0) << endl;
     235             :     cout << "timehi=" << timehi-timelo(0) << endl;
     236             :     cout << "ord.nelements() = " << ord.nelements() << endl;
     237             :     cout << "ord = " << ord << endl;
     238             :     */
     239             : 
     240           0 :   }
     241             : 
     242             :   //UNUSED: Double rTime=timelo(ord(0));
     243             : 
     244             :   // Now iterate throught the NCT and set field and scan according to time
     245           0 :   Block<String> cols(1);
     246           0 :   cols[0]="TIME";
     247           0 :   CTIter ctiter(ct,cols);
     248             : 
     249           0 :   Int itime(0);
     250           0 :   Int thisObs(0);
     251           0 :   Int thisScan(0);
     252           0 :   Int thisField(0);
     253           0 :   while (!ctiter.pastEnd()) {
     254           0 :     Double thisTime=ctiter.thisTime();
     255             : 
     256             :     //    cout.precision(12);
     257             :     //    cout << "thisTime = " << thisTime-rTime << endl;
     258             : 
     259             :     // If time before first MS time, just use first
     260           0 :     if (thisTime<timehi(ord(0))) {
     261           0 :       itime=0;
     262             :       //      cout << " Pre: ";
     263             :     }
     264             :     // If time after last MS time, use last
     265           0 :     else if (thisTime>timelo(ord(nScan-1))) {
     266           0 :       itime=nScan-1;
     267             :       //      cout << " Post: ";
     268             :     }
     269           0 :     else if (thisTime>timehi(ord(itime))) {
     270             : 
     271             :       // Isolate correct time index
     272           0 :       while (thisTime>timehi(ord(itime))&& itime<nScan) {
     273             :         //      cout << itime << " " << thisTime-rTime << ">" << timehi(ord(itime))-rTime << endl;
     274           0 :         ++itime;
     275             :       }
     276             :       //      cout << " Found: ";
     277             :     }
     278             :     //else 
     279             :       //      cout << " Still: ";
     280             : 
     281           0 :     thisObs=obslist(ord(itime));
     282           0 :     thisScan=scanlist(ord(itime));
     283           0 :     thisField=fieldlist(ord(itime));
     284             : 
     285             :     /*
     286             :     cout << " itime=" << itime << " "
     287             :          << timelo(ord(itime))-rTime << " < "
     288             :          << thisTime-rTime << " < "
     289             :          << timehi(ord(itime))-rTime
     290             :          << " s=" << thisScan << " f=" << thisField << endl;
     291             :     */
     292             : 
     293             :     // Set the field and scan
     294           0 :     if (doField) 
     295           0 :       ctiter.setfield(thisField);
     296           0 :     if (doScan) 
     297           0 :       ctiter.setscan(thisScan);
     298           0 :     if (doObs) 
     299           0 :       ctiter.setobs(thisObs);
     300             :     
     301           0 :     ctiter.next();
     302             :   }
     303             : 
     304           0 : }
     305             : 
     306             : } //# NAMESPACE CASA - END