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/PKSMS2writer.cc

    r1325 r1393  
    3838#include <casa/BasicSL/Constants.h>
    3939#include <casa/Quanta/QC.h>
     40#include <casa/Logging/LogIO.h>
    4041#include <measures/Measures/Stokes.h>
    4142#include <tables/Tables/ArrColDesc.h>
    … …  
    8182        const Vector<uInt> nPol,
    8283        const Vector<Bool> haveXPol,
    83         const Bool   haveBase)
     84        const Bool   haveBase,
     85        const String fluxUnit)
    8486{
    8587  // Open a MS table.
    … …  
    9294  Int maxNPol = max(cNPol);
    9395
     96  // check if it is GBT data
     97  cGBT = antName.contains("GBT");
    9498
    9599  // Add the non-standard CALFCTR column.
    … …  
    99103  // Add the optional FLOAT_DATA column.
    100104  MS::addColumnToDesc(pksDesc, MS::FLOAT_DATA, 2);
     105  //pksDesc.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().
     106  //              define("UNIT", String("Jy"));
    101107  pksDesc.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().
    102                 define("UNIT", String("Jy"));
     108                define("UNIT", fluxUnit);
    103109  pksDesc.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().
    104110                define("MEASURE_TYPE", "");
    … …  
    350356  } else if (dopplerFrame == "SOURCE") {
    351357    MFrequency::getType(cDopplerFrame, "REST");
     358  } else if (dopplerFrame == "LSRK") {
     359    MFrequency::getType(cDopplerFrame, "LSRK");
    352360  }
    353361
    … …  
    357365  addDopplerEntry();
    358366  addFeedEntry();
    359   addObservationEntry(observer, project);
     367  //addObservationEntry(observer, project);
     368  addObservationEntry(observer, project, antName);
    360369  addProcessorEntry();
    361370
    … …  
    446455
    447456  // SYSCAL subtable.
    448   addSysCalEntry(beamNo, iIF, time, interval, tcal, tsys);
     457  addSysCalEntry(beamNo, iIF, time, interval, tcal, tsys, nPol);
     458
    449459
    450460  // Handle weather information.
    … …  
    499509    }
    500510  }
    501 
    502511  // Transpose spectra.
    503512  Matrix<Float> tmpData(nPol, nChan);
    … …  
    509518    }
    510519  }
    511 
    512520  cCalFctrCol->put(irow, calFctr);
    513521  cMSCols->floatData().put(irow, tmpData);
    … …  
    522530  cMSCols->sigma().put(irow, sigma);
    523531
    524   Vector<Float> weight(1, 1.0f);
     532  //Vector<Float> weight(1, 1.0f);
     533  Vector<Float> weight(nPol, 1.0f);
    525534  cMSCols->weight().put(irow, weight);
    526535
    527536  // Flag information.
    528537  Cube<Bool> flags(nPol, nChan, 1, False);
    529   cMSCols->flag().put(irow, flags.xyPlane(0));
     538  //cMSCols->flag().put(irow, flags.xyPlane(0));
    530539  cMSCols->flagCategory().put(irow, flags);
    531540  cMSCols->flagRow().put(irow, False);
    … …  
    601610  Int n = cAntenna.nrow() - 1;
    602611
     612  // do specific things for GBT
    603613  // Data.
    604614  cAntennaCols->name().put(n, antName);
    605   cAntennaCols->station().put(n, "ATNF_PARKES");
     615  //cAntennaCols->station().put(n, "ATNF_PARKES");
     616  if (cGBT) {
     617    cAntennaCols->station().put(n, "GREENBANK");
     618  }
     619  else {
     620    cAntennaCols->station().put(n, "ATNF_PARKES");
     621  }
    606622  cAntennaCols->type().put(n, "GROUND-BASED");
    607623  cAntennaCols->mount().put(n, "ALT-AZ");
    … …  
    609625  Vector<Double> antOffset(3, 0.0);
    610626  cAntennaCols->offset().put(n, antOffset);
    611   cAntennaCols->dishDiameter().put(n, 64.0);
    612 
     627  //cAntennaCols->dishDiameter().put(n, 64.0);
     628  if (cGBT) {
     629    cAntennaCols->dishDiameter().put(n, 110.0);
     630  }
     631  else {
     632    cAntennaCols->dishDiameter().put(n, 64.0);
     633  }
    613634  // Flags.
    614635  cAntennaCols->flagRow().put(n, False);
    … …  
    723744        const Int srcId)
    724745{
     746
     747  ROScalarColumn<String> fldn(cField, "NAME");
     748  ROScalarColumn<Int> sourceid(cField, "SOURCE_ID");
     749  Int n;
     750  Int nFld = cField.nrow();
     751  for (n = 0; n < nFld; n++) {
     752    if (fldn(n) == fieldName && sourceid(n) == srcId) {
     753      break;
     754    }
     755  }
     756
    725757  // Extend the FIELD subtable.
    726   cField.addRow();
    727   Int n = cField.nrow() - 1;
    728 
    729   // Data.
    730   cFieldCols->name().put(n, fieldName);
    731   cFieldCols->code().put(n, "DRIFT");
    732   cFieldCols->time().put(n, time);
    733 
    734   Matrix<Double> track(2, 2);
    735   track.column(0) = direction;
    736   track.column(1) = scanRate;
    737   cFieldCols->numPoly().put(n, 1);
    738   cFieldCols->delayDir().put(n, track);
    739   cFieldCols->phaseDir().put(n, track);
    740   cFieldCols->referenceDir().put(n, track);
    741   cFieldCols->sourceId().put(n, srcId);
    742 
    743   // Flags.
    744   cFieldCols->flagRow().put(n, False);
     758  if (n == nFld) {
     759    cField.addRow();
     760    //Int n = cField.nrow() - 1;
     761
     762    // Data.
     763    cFieldCols->name().put(n, fieldName);
     764    if (cGBT) {
     765      cFieldCols->code().put(n, " ");
     766    }
     767    else {
     768      cFieldCols->code().put(n, "DRIFT");
     769    }
     770    cFieldCols->time().put(n, time);
     771
     772    //Matrix<Double> track(2, 2);
     773    Matrix<Double> track(2, 1);
     774    track.column(0) = direction;
     775    //track.column(1) = scanRate;
     776    cFieldCols->numPoly().put(n, 1);
     777    cFieldCols->delayDir().put(n, track);
     778    cFieldCols->phaseDir().put(n, track);
     779    cFieldCols->referenceDir().put(n, track);
     780    cFieldCols->sourceId().put(n, srcId);
     781
     782    // Flags.
     783    cFieldCols->flagRow().put(n, False);
     784  }
    745785
    746786  return n;
    … …  
    753793Int PKSMS2writer::addObservationEntry(
    754794        const String observer,
    755         const String project)
     795        const String project,
     796        const String antName)
    756797{
    757798  // Extend the OBSERVATION subtable.
    … …  
    760801
    761802  // Data.
    762   cObservationCols->telescopeName().put(n, "Parkes");
     803  //cObservationCols->telescopeName().put(n, "Parkes");
     804  cObservationCols->telescopeName().put(n, antName);
    763805  Vector<Double> timerange(2, 0.0);
    764806  cObservationCols->timeRange().put(n, timerange);
    … …  
    766808  Vector<String> log(1, "none");
    767809  cObservationCols->log().put(n, log);
    768   cObservationCols->scheduleType().put(n, "ATNF");
     810  //cObservationCols->scheduleType().put(n, "ATNF");
     811  cObservationCols->scheduleType().put(n, "");
    769812  Vector<String> schedule(1, "Not available");
    770813  cObservationCols->schedule().put(n, schedule);
    … …  
    779822
    780823//--------------------------------------------- PKSMS2writer::addPointingEntry
     824
     825// Modified to fill pointing data if the direction is the pointing direction.
     826// So the following comment is no longer true.
    781827
    782828// Add an entry to the POINTING subtable.  This compulsory subtable simply
    … …  
    790836        const Vector<Double> scanRate)
    791837{
    792   // Extend the POINTING subtable.
    793   cPointing.addRow();
    794   Int n = cPointing.nrow() - 1;
    795 
    796   // Keys.
    797   cPointingCols->antennaId().put(n, 0);
    798   cPointingCols->time().put(n, time);
    799   cPointingCols->interval().put(n, interval);
    800 
    801   // Data.
    802   cPointingCols->name().put(n, fieldName);
    803   cPointingCols->numPoly().put(n, 1);
    804   cPointingCols->timeOrigin().put(n, time);
    805 
    806   Matrix<Double> track(2, 2);
    807   track.column(0) = direction;
    808   track.column(1) = scanRate;
    809   cPointingCols->direction().put(n, track);
    810   cPointingCols->target().put(n, track);
    811   cPointingCols->tracking().put(n, True);
    812 
     838
     839  ROScalarColumn<Double> tms(cPointing, "TIME");
     840  Int n;
     841  Int ntm = cPointing.nrow();
     842  for (n = 0; n < ntm; n++) {
     843    if (tms(n) == time) {
     844      break;
     845    }
     846  }
     847
     848  if (n == ntm) {
     849    // Extend the POINTING subtable.
     850    cPointing.addRow();
     851    //Int n = cPointing.nrow() - 1;
     852
     853    // Keys.
     854    cPointingCols->antennaId().put(n, 0);
     855    cPointingCols->time().put(n, time);
     856    cPointingCols->interval().put(n, interval);
     857
     858    // Data.
     859    cPointingCols->name().put(n, fieldName);
     860    cPointingCols->numPoly().put(n, 1);
     861    cPointingCols->timeOrigin().put(n, time);
     862
     863    //Matrix<Double> track(2, 2);
     864    Matrix<Double> track(2, 1);
     865    track.column(0) = direction;
     866    //track.column(1) = scanRate;
     867    cPointingCols->direction().put(n, track);
     868    cPointingCols->target().put(n, track);
     869    cPointingCols->tracking().put(n, True);
     870  }
    813871  return n;
    814872}
    … …  
    833891  // Data.
    834892  Vector<Int> corrType(2);
     893  if (nPol == 1) {
     894  corrType.resize(1);
     895  corrType(0) = Stokes::XX;
     896  }
     897  else {
     898  //Vector<Int> corrType(2);
    835899  corrType(0) = Stokes::XX;
    836900  corrType(1) = Stokes::YY;
     901  }
    837902  cPolarizationCols->corrType().put(n, corrType);
    838903
    839904  Matrix<Int> corrProduct(2,2,1);
     905  if (nPol == 1) {
     906    corrProduct.resize(2,1,1);
     907    corrProduct(1,0) = 0;
     908  }
    840909  if (nPol == 2) {
    841910    corrProduct(1,0) = 0;
    … …  
    9451014
    9461015  // Data.
    947   cSpWindowCols->name().put(n, "L-band");
     1016  //cSpWindowCols->name().put(n, "L-band");
     1017  cSpWindowCols->name().put(n, " ");
    9481018  cSpWindowCols->refFrequency().put(n, refFreq);
    9491019
    … …  
    10271097        const Double interval,
    10281098        const Vector<Float> tcal,
    1029         const Vector<Float> tsys)
    1030 {
     1099        const Vector<Float> tsys,
     1100        const Int nPol)
     1101{
     1102  LogIO os(LogOrigin("PKSMS2writer", "addSysCalEntry()", WHERE));
     1103
    10311104  // Extend the SYSCAL subtable.
    10321105  cSysCal.addRow();
    10331106  Int n = cSysCal.nrow() - 1;
    10341107
     1108  //check fo consistency with n pol
     1109  //here assume size of Tcal vector = npol
     1110  Vector<Float> inTcal(nPol,0);
     1111  Int ndim = tcal.shape()(0);
     1112  Vector<Float> tmpTcal = tcal;
     1113  if (nPol != ndim) {
     1114    os << LogIO::WARN
     1115       << "Found "<< ndim <<" Tcal value(s) for the data with "<<nPol<<" polarization(s)"
     1116       << "(expecting one Tcal per pol)."<<endl
     1117       << "First "<< nPol << " Tcal value(s) will be filled." << LogIO::POST;
     1118    tmpTcal.resize(nPol, True);
     1119    inTcal = tmpTcal;
     1120  }
    10351121  // Keys.
    10361122  cSysCalCols->antennaId().put(n, 0);
    … …  
    10411127
    10421128  // Data.
    1043   cSysCalCols->tcal().put(n, tcal);
     1129  //cSysCalCols->tcal().put(n, tcal);
     1130  cSysCalCols->tcal().put(n, inTcal);
    10441131  cSysCalCols->tsys().put(n, tsys);
    10451132
Note: See TracChangeset for help on using the changeset viewer.