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           8 : void smoothCT(NewCalTable ct,
      50             :               const String& smtype,
      51             :               const Double& smtime,
      52             :               Vector<Int> selfields) {
      53             : 
      54             :   // Complex parameters?
      55           8 :   Bool cmplx=ct.isComplex();
      56             : 
      57             :   // half-width
      58           8 :   Double thw(smtime/2.0);
      59             : 
      60             :   // Workspace
      61           8 :   Vector<Double> times;
      62           8 :   Vector<Float> p,newp;
      63           8 :   Vector<Bool> pOK, newpOK;
      64             : 
      65           8 :   Cube<Float> fpar;
      66           8 :   Cube<Bool> fparok,newfparok;
      67             : 
      68           8 :   Vector<Bool> mask;
      69             : 
      70           8 :   IPosition blc(3,0,0,0), fblc(3,0,0,0);
      71           8 :   IPosition trc(3,0,0,0), ftrc(3,0,0,0);
      72           8 :   IPosition vec(1,0);
      73             : 
      74           8 :   Block<String> cols(4);
      75           8 :   cols[0]="SPECTRAL_WINDOW_ID";
      76           8 :   cols[1]="FIELD_ID";
      77           8 :   cols[2]="ANTENNA1";
      78           8 :   cols[3]="ANTENNA2";
      79           8 :   CTIter ctiter(ct,cols);
      80             : 
      81        4391 :   while (!ctiter.pastEnd()) {
      82             : 
      83        4383 :     Int nSlot=ctiter.nrow();
      84        4383 :     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        8793 :     if (nSlot>1 && 
      89        4410 :         (selfields.nelements()<1 || anyEQ(selfields,ifld))) {
      90             : 
      91             :       //UNUSED: Int ispw=ctiter.thisSpw();
      92             :       
      93        4236 :       vec(0)=nSlot;
      94        4236 :       trc(2)=ftrc(2)=nSlot-1;
      95             : 
      96        4236 :       times.assign(ctiter.time());
      97             : 
      98             :       // Extract Float info
      99        4236 :       if (cmplx)
     100        4236 :         fpar.assign(ctiter.casfparam("AP"));
     101             :       else
     102           0 :         fpar.assign(ctiter.fparam());
     103             : 
     104        4236 :       fparok.assign(!ctiter.flag());
     105        4236 :       newfparok.assign(fparok);
     106        4236 :       IPosition fsh(fpar.shape());
     107             : 
     108             :       // For each channel
     109        8472 :       for (int ichan=0;ichan<fsh(1);++ichan) {
     110        4236 :         blc(1)=trc(1)=fblc(1)=ftrc(1)=ichan;
     111             :         // For each param (pol)
     112       13404 :         for (Int ipar=0;ipar<fsh(0);++ipar) {
     113        9168 :           blc(0)=trc(0)=ipar;
     114        9168 :           fblc(0)=ftrc(0)=ipar/(cmplx?2:1);
     115             :           
     116             :           // Reference slices of par/parOK
     117        9168 :           p.reference(fpar(blc,trc).reform(vec));
     118        9168 :           newp.assign(p);
     119        9168 :           pOK.reference(fparok(fblc,ftrc).reform(vec));
     120        9168 :           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        9168 :           Vector<Bool> mask;
     133      191264 :           for (Int i=0;i<nSlot;++i) {
     134             :             // Make mask
     135      182096 :             mask = pOK;
     136      546288 :             mask = (mask && ( (times >  (times(i)-thw)) && 
     137      546288 :                               (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      182096 :             if (ntrue(mask)>0) {
     149      169512 :               if (smtype=="mean") {
     150      143804 :                 newp(i)=mean(p(mask));
     151             :               }
     152       25708 :               else if (smtype=="median") {
     153       25708 :                 newp(i)= median(p(mask),false);
     154             :               }
     155      169512 :               newpOK(i)=true;
     156             :             }
     157             :             else 
     158       12584 :               newpOK(i)=false;
     159             :             
     160             :           } // i
     161             :           // keep new ok info
     162        9168 :           p=newp;
     163        9168 :         } // ipar
     164             :       } // ichan
     165             : 
     166             :       // Put info back
     167        4236 :       if (cmplx)
     168        4236 :         ctiter.setcparam(RIorAPArray(fpar).c());
     169             :       else
     170           0 :         ctiter.setfparam(fpar);
     171             : 
     172        4236 :       ctiter.setflag(!newfparok);
     173             : 
     174        4236 :     } // nSlot>1
     175             : 
     176        4383 :     ctiter.next();
     177             :   } // ispw
     178             :         
     179           8 : }
     180             : 
     181          46 : 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          46 :   Table mstab(msName,Table::Old);
     188             : 
     189             :   // How many scans in total?
     190          46 :   TableVector<Int> allscansTV(mstab,"SCAN_NUMBER");
     191          46 :   Vector<Int> allscans=allscansTV.makeVector();
     192          46 :   Int nScan=genSort(allscans,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
     193             : 
     194             :   //  cout << "Found " << nScan << " scans in " << msName << "." << endl;
     195             : 
     196             :   // Workspace
     197          46 :   Vector<Int> scanlist(nScan,-1);
     198          46 :   Vector<Double> timelo(nScan,DBL_MIN);
     199          46 :   Vector<Double> timehi(nScan,DBL_MAX);
     200          46 :   Vector<Int> fieldlist(nScan,-1);
     201          46 :   Vector<Int> obslist(nScan,-1);
     202          46 :   Vector<uInt> ord;
     203             :   {
     204          46 :     Block<String> cols(1);
     205          46 :     cols[0]="SCAN_NUMBER";
     206          46 :     TableIterator mstiter(mstab,cols);
     207             :     
     208             :     // Get time boundares and fields for each scan
     209          46 :     Int iscan(0);
     210         176 :     while (!mstiter.pastEnd()) {
     211         130 :       Table thistab(mstiter.table());
     212             :       
     213         130 :       Int scan=TableVector<Int>(thistab,"SCAN_NUMBER")(0);
     214         130 :       scanlist(iscan)=scan;
     215             :       
     216         130 :       fieldlist(iscan)=TableVector<Int>(thistab,"FIELD_ID")(0);
     217         130 :       obslist(iscan)=TableVector<Int>(thistab,"OBSERVATION_ID")(0);
     218             :       
     219         260 :       Vector<Double> times=TableVector<Double>(thistab,"TIME").makeVector();
     220         130 :       timelo(iscan)=min(times)-1e-5;
     221         130 :       timehi(iscan)=max(times)+1e-5;
     222             :       
     223         130 :       mstiter.next();
     224         130 :       ++iscan;
     225         130 :     }
     226             : 
     227             : 
     228             :     // Ensure time orderliness
     229          46 :     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          46 :   }
     241             : 
     242             :   //UNUSED: Double rTime=timelo(ord(0));
     243             : 
     244             :   // Now iterate throught the NCT and set field and scan according to time
     245          46 :   Block<String> cols(1);
     246          46 :   cols[0]="TIME";
     247          46 :   CTIter ctiter(ct,cols);
     248             : 
     249          46 :   Int itime(0);
     250          46 :   Int thisObs(0);
     251          46 :   Int thisScan(0);
     252          46 :   Int thisField(0);
     253        7814 :   while (!ctiter.pastEnd()) {
     254        7768 :     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        7768 :     if (thisTime<timehi(ord(0))) {
     261         384 :       itime=0;
     262             :       //      cout << " Pre: ";
     263             :     }
     264             :     // If time after last MS time, use last
     265        7384 :     else if (thisTime>timelo(ord(nScan-1))) {
     266        7325 :       itime=nScan-1;
     267             :       //      cout << " Post: ";
     268             :     }
     269          59 :     else if (thisTime>timehi(ord(itime))) {
     270             : 
     271             :       // Isolate correct time index
     272          41 :       while (thisTime>timehi(ord(itime))&& itime<nScan) {
     273             :         //      cout << itime << " " << thisTime-rTime << ">" << timehi(ord(itime))-rTime << endl;
     274          28 :         ++itime;
     275             :       }
     276             :       //      cout << " Found: ";
     277             :     }
     278             :     //else 
     279             :       //      cout << " Still: ";
     280             : 
     281        7768 :     thisObs=obslist(ord(itime));
     282        7768 :     thisScan=scanlist(ord(itime));
     283        7768 :     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        7768 :     if (doField) 
     295        7768 :       ctiter.setfield(thisField);
     296        7768 :     if (doScan) 
     297        7768 :       ctiter.setscan(thisScan);
     298        7768 :     if (doObs) 
     299        7768 :       ctiter.setobs(thisObs);
     300             :     
     301        7768 :     ctiter.next();
     302             :   }
     303             : 
     304          46 : }
     305             : 
     306             : } //# NAMESPACE CASA - END