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