Ignore:
Timestamp:
08/01/07 02:54:51 (17 years ago)
Author:
TakTsutsumi
Message:

Changes to handle GBT MS data

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/external/atnf/PKSIO/PKSMS2reader.cc

    r1325 r1393  
    2626//#                        Charlottesville, VA 22903-2475 USA
    2727//#
    28 //# $Id: PKSMS2reader.cc,v 19.11 2006/07/05 04:59:20 mcalabre Exp $
     28//# $Id$
    2929//#---------------------------------------------------------------------------
    3030//# Original: 2000/08/03, Mark Calabretta, ATNF
    … …  
    3838#include <ms/MeasurementSets/MSColumns.h>
    3939#include <tables/Tables.h>
     40#include <casa/Quanta/MVTime.h>
     41#include <casa/Quanta/MVAngle.h>
     42#include <casa/BasicMath/Math.h>
     43#include <measures/Measures/MeasConvert.h>
     44#include <measures/Measures/MEpoch.h>
     45#include <measures/Measures/MeasRef.h>
     46
    4047
    4148// Parkes includes.
    … …  
    8491
    8592  cPKSMS  = MeasurementSet(msName);
     93  // taql access to the syscal table
     94  String tmpcalName = msName;
     95  tmpcalName.append("/SYSCAL");
     96  cSysCalTab = Table(tmpcalName);
     97
    8698  cIdx    = 0;
     99  lastmjd = 0.0;
    87100  cNRow   = cPKSMS.nrow();
    88101  cMSopen = True;
    … …  
    111124  cFieldIdCol.reference(msCols.fieldId());
    112125  cFieldNameCol.reference(fieldCols.name());
     126  cFieldDelayDirCol.reference(fieldCols.delayDir());
    113127
    114128  cSrcIdCol.reference(fieldCols.sourceId());
     129  cSrcId2Col.reference(sourceCols.sourceId());
    115130  cSrcNameCol.reference(sourceCols.name());
    116131  cSrcDirCol.reference(sourceCols.direction());
    … …  
    120135  cStateIdCol.reference(msCols.stateId());
    121136  cObsModeCol.reference(stateCols.obsMode());
    122 
     137  cCalCol.reference(stateCols.cal());
     138  cSigStateCol.reference(stateCols.sig());
     139  cRefStateCol.reference(stateCols.ref());
    123140  cDataDescIdCol.reference(msCols.dataDescId());
     141  cSpWinIdCol.reference(dataDescCols.spectralWindowId());
    124142  cChanFreqCol.reference(spWinCols.chanFreq());
    125143
    … …  
    142160    cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
    143161  }
     162 
     163  if ((cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {
     164    cTcalCol.attach(cPKSMS.sysCal(), "TCAL");
     165  }
    144166
    145167  if ((cHaveCalFctr = cPKSMS.tableDesc().isColumn("CALFCTR"))) {
    … …  
    157179  cFlagCol.reference(msCols.flag());
    158180
     181
    159182  if ((cGetXPol = cPKSMS.isColumn(MSMainEnums::DATA))) {
    160183    if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
    … …  
    177200
    178201  // Number of IFs.
    179   uInt nIF = dataDescCols.nrow();
     202  //uInt nIF = dataDescCols.nrow();
     203  uInt nIF =spWinCols.nrow();
    180204  IFs.resize(nIF);
    181205  IFs = True;
    … …  
    253277        Double &mjd,
    254278        Double &refFreq,
    255         Double &bandwidth)
     279        Double &bandwidth,
     280        String &fluxunit)
    256281{
    257282  if (!cMSopen) {
    … …  
    277302  }
    278303
     304  fluxunit = "";
     305  const TableRecord& keywordSet
     306     = cFloatDataCol.columnDesc().keywordSet();
     307  if(keywordSet.isDefined("UNIT")) {
     308    fluxunit = keywordSet.asString("UNIT");
     309  }
    279310
    280311  // Coordinate equinox.
    … …  
    538569  Int ibeam;
    539570  Int iIF;
     571  Int iDataDesc;
     572
    540573  while (True) {
    541574    ibeam = cBeamNoCol(cIdx);
    542     iIF   = cDataDescIdCol(cIdx);
     575    iDataDesc   = cDataDescIdCol(cIdx);
     576    iIF   =cSpWinIdCol(iDataDesc);
    543577    if (cBeams(ibeam) && cIFs(iIF)) {
    544578      break;
    … …  
    552586
    553587  // Renumerate scan no. Here still is 1-based
    554   scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
    555 
     588  //scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
     589  scanNo = cScanNoCol(cIdx);
     590 
    556591  if (scanNo != cScanNo) {
    557592    // Start of new scan.
    … …  
    576611
    577612  Int srcId = cSrcIdCol(fieldId);
    578   srcName = cSrcNameCol(srcId);
     613  //For source with multiple spectral window setting, this is not
     614  // correct. Source name of srcId may not be at 'srcId'th row of SrcNameCol
     615  //srcName = cSrcNameCol(srcId);
     616  for (uInt irow = 0; irow < cSrcId2Col.nrow(); irow++) {
     617    if (cSrcId2Col(irow) == srcId) {
     618      srcName = cSrcNameCol(irow);
     619    }
     620  }
     621
    579622  srcDir  = cSrcDirCol(srcId);
    580623  srcPM   = cSrcPMCol(srcId);
    … …  
    587630  }
    588631
     632  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
     633  String telescope = antennaCols.name()(0);
     634  Bool cGBT = telescope.contains("GBT");
    589635  // Observation type.
    590636  Int stateId = cStateIdCol(cIdx);
    591   obsMode = cObsModeCol(stateId);
     637  if (stateId == -1) {
     638    //obsMode = " ";
     639  } else {
     640    obsMode = cObsModeCol(stateId);
     641    Bool sigState =cSigStateCol(stateId);
     642    Bool refState =cRefStateCol(stateId);
     643    //DEBUG
     644    //cerr <<"stateid="<<stateId<<" obsmode="<<obsMode<<endl;
     645    if (cGBT) {
     646      // split the obsMode string and append a proper label
     647      // (these are GBT specific)
     648      int epos = obsMode.find_first_of(':');
     649      int nextpos = obsMode.find_first_of(':',epos+1);
     650      string obsMode1 = obsMode.substr(0,epos);
     651      string obsMode2 = obsMode.substr(epos+1,nextpos-epos-1);
     652     
     653      //cerr <<"obsMode2= "<<obsMode2<<endl;
     654      if (!srcName.contains("_ps")
     655          &&!srcName.contains("_psr")
     656          &&!srcName.contains("_nod")
     657          &&!srcName.contains("_fs")
     658          &&!srcName.contains("_fsr")) {
     659        // if Nod mode observation , append '_nod'
     660        if (obsMode1 == "Nod") {
     661          srcName.append("_nod");
     662        } else if (obsMode1 == "OffOn") {
     663        // for GBT position switch observations (OffOn or OnOff)
     664          if (obsMode2 == "PSWITCHON") srcName.append("_ps");
     665          if (obsMode2 == "PSWITCHOFF") srcName.append("_psr");
     666        } else {
     667          if (obsMode2 == "FSWITCH") {
     668          // for GBT frequency switch mode
     669            if (sigState) srcName.append("_fs");
     670            if (refState) srcName.append("_fsr");
     671          }
     672        }
     673      }
     674    }
     675  }
     676  // CAL state
     677  // this should be apply just for GBT data?
     678  Double Cal;
     679  if (stateId==-1) {
     680    Cal = 0;
     681  } else {
     682    Cal = cCalCol(stateId);
     683  }
     684  if (cGBT) {
     685    if (Cal > 0 && !srcName.contains("_calon")) {
     686      srcName.append("_calon");
     687    }
     688  }
    592689
    593690  IFno = iIF + 1;
    594691  Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
    595 
     692 
    596693  // Minimal handling on continuum data.
    597694  Vector<Double> chanFreq = cChanFreqCol(iIF);
    … …  
    607704      freqInc = chanFreq(0) - chanFreq(1);
    608705    }
    609 
    610706    refFreq  = chanFreq(cRefChan(iIF)-1);
    611707    restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
    … …  
    616712  tcal      = 0.0f;
    617713  tcalTime  = "";
    618   azimuth   = 0.0f;
    619   elevation = 0.0f;
     714  //azimuth   = 0.0f;
     715  //elevation = 0.0f;
    620716  parAngle  = 0.0f;
    621717  focusAxi  = 0.0f;
    … …  
    649745  beamNo  = ibeam + 1;
    650746
    651   Matrix<Double> pointingDir = cPointingCol(fieldId);
    652   direction = pointingDir.column(0);
    653   uInt ncols = pointingDir.ncolumn();
     747  //Matrix<Double> pointingDir = cPointingCol(fieldId);
     748  //pointingDir = cPointingCol(fieldId);
     749  //direction = pointingDir.column(0);
     750  //uInt ncols = pointingDir.ncolumn();
     751  //if (ncols == 1) {
     752  //  scanRate = 0.0f;
     753  //} else {
     754  //  scanRate  = pointingDir.column(1);
     755  //}
     756
     757  // Get direction from FIELD table
     758  // here, assume direction to be the field direction not pointing
     759  Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
     760  direction = delayDir.column(0);
     761  uInt ncols = delayDir.ncolumn();
    654762  if (ncols == 1) {
    655763    scanRate = 0.0f;
    656764  } else {
    657     scanRate  = pointingDir.column(1);
    658   }
     765    scanRate  = delayDir.column(1);
     766  }
     767
     768  // caluculate azimuth and elevation
     769  // first, get the reference frame
     770  MVPosition mvpos(antennaCols.position()(0));
     771  MPosition mp(mvpos);
     772  Quantum<Double> qt(time,"s");
     773  MVEpoch mvt(qt);
     774  MEpoch me(mvt);
     775  MeasFrame frame(mp, me);
     776  //
     777  ROMSFieldColumns fldCols(cPKSMS.field());
     778  Vector<MDirection> vmd(1);
     779  MDirection md;
     780  fldCols.delayDirMeasCol().get(fieldId,vmd);
     781  md = vmd[0];
     782  //Vector<Double> dircheck = md.getAngle("rad").getValue();
     783  //cerr<<"dircheck="<<dircheck<<endl;
     784
     785  Vector<Double> azel =
     786        MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
     787                                                frame)
     788                            )().getAngle("rad").getValue();
     789  //cerr<<"azel="<<azel<<endl;
     790  azimuth = azel[0];
     791  elevation = azel[1];
    659792
    660793  // Get Tsys assuming that entries in the SYSCAL table match the main table.
    … …  
    675808  cSigmaCol.get(cIdx, sigma, True);
    676809
     810  //get Tcal if available
     811  if (cHaveTcal) {
     812    Int nTcalColRow = cTcalCol.nrow();
     813    uInt nBeam = cBeams.nelements();
     814    uInt nIF = cIFs.nelements();
     815    uInt nrws = nBeam * nIF;
     816    if (nTcalColRow > 0) { 
     817    // find tcal match with the data with the data time stamp
     818      Double mjds = mjd*(24*3600);
     819      Double dtcalTime;
     820      if ( mjd > lastmjd || cIdx==0 ) {
     821        //Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
     822        tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
     823        //DEBUG
     824        //if (cIdx == 0) {
     825        //  cerr<<"inital table retrieved"<<endl;
     826        //}
     827       
     828      }
     829
     830      if (nBeam == 1) {
     831        tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
     832      } else {
     833        tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
     834                              tmptab.col("FEED_ID") == ibeam , 1);
     835      }
     836      //cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
     837      int syscalrow = tmptab2.nrow();
     838      ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
     839      ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
     840      if (syscalrow==0) {
     841        cerr<<"Cannot find any matching Tcal at/near the data timestamp."
     842           << " Set Tcal=0.0"<<endl;
     843      } else {
     844        tcalCol.get(0, tcal);
     845        tcalTimeCol.get(0,dtcalTime);
     846        tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
     847        //DEBUG
     848        //cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
     849        tmptab.markForDelete();
     850        tmptab2.markForDelete();
     851      }
     852    }
     853    lastmjd = mjd;
     854  }
     855
    677856  // Calibration factors (if available).
    678857  calFctr.resize(cNPol(iIF));
    … …  
    695874    baseSub.resize(0,0);
    696875  }
    697 
    698876
    699877  // Get spectral data.
    … …  
    780958  Int ibeam;
    781959  Int iIF;
     960  Int iDataDesc;
    782961  while (True) {
    783962    ibeam = cBeamNoCol(cIdx);
    784     iIF   = cDataDescIdCol(cIdx);
     963    //iIF   = cDataDescIdCol(cIdx);
     964    iDataDesc   = cDataDescIdCol(cIdx);
     965    iIF   = cSpWinIdCol(iDataDesc);
    785966    if (cBeams(ibeam) && cIFs(iIF)) {
    786967      break;
    … …  
    794975
    795976  IFno = iIF + 1;
    796 
    797977  // Get Tsys assuming that entries in the SYSCAL table match the main table.
    798978  cTsysCol.get(cIdx, tsys, True);
Note: See TracChangeset for help on using the changeset viewer.