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