Commits

George Moellenbrock authored 8a8f3645981 Merge
Merge branch 'master' into CAS-12168

code/alma/apps/asdm2MS/asdm2MS.cc

Modified
23 23 #include <string>
24 24 #include <regex>
25 25 #include <algorithm>
26 26 #include <vector>
27 27 #include <iomanip>
28 28
29 29 #include <stdcasa/optionparser.h>
30 30 #include <alma/Options/AlmaArg.h>
31 31 using namespace alma;
32 32
33 +#include <alma/Enumerations/CWindowFunction.h>
33 34 #include <alma/ASDM/ASDMAll.h>
34 35 #include <alma/ASDM/Misc.h>
35 36
36 37 #include <alma/ASDMBinaries/SDMBinData.h>
37 38 using namespace sdmbin;
38 39
39 40 #include <exception>
40 41 using namespace asdm;
41 42 #include <alma/ASDM/IllegalAccessException.h>
42 43
1844 1845 error(errstream.str());
1845 1846 } catch ( std::exception & e) {
1846 1847 errstream.str("");
1847 1848 errstream << e.what();
1848 1849 error(errstream.str());
1849 1850 }
1850 1851
1851 1852 LOGEXIT("fillField");
1852 1853 }
1853 1854
1855 +/**
1856 + * Utility function to sort out the numBin value given appropriate for a given row in the SpectralWindow table
1857 + */
1858 +int getNumBin(SpectralWindowRow *spwRow, const string &telescopeName) {
1859 + // assumes that spRow is already a non-null pointer
1860 + if (spwRow->isNumBinExists()) {
1861 + // use any numBin value that already exists
1862 + return spwRow->getNumBin();
1863 + } else {
1864 + // 1 will be returned if numBin <=0 at the end
1865 + int numBin = 0;
1866 +
1867 + // infer it based on the telescope name
1868 + // both methods require that these 2 scalar values exist
1869 + if (spwRow->isChanWidthExists() && spwRow->isResolutionExists()) {
1870 + if (telescopeName == string("EVLA")) {
1871 + numBin = std::nearbyint(spwRow->getChanWidth().get()/spwRow->getResolution().get());
1872 + } else if (telescopeName == string("ALMA")) {
1873 + // only relevant for HANNING
1874 + if (spwRow->getWindowFunction()==WindowFunctionMod::HANNING) {
1875 + // also requires that effectiveBw exists
1876 + if (spwRow->isEffectiveBwExists()) {
1877 + double chanWidth = spwRow->getChanWidth().get();
1878 + if (spwRow->getResolution().get() != chanWidth) {
1879 + // effectiveBw and resolution have been set according to the Hills memo, get implied numBin
1880 + // the expected values of the ratio of effectiveBw/chanWidth are from the Hills memo
1881 + // second table on p. 6
1882 + double ratio = spwRow->getEffectiveBw().get()/std::abs(chanWidth);
1883 + // no other way to do this than one value at a time, starting from numBin=1
1884 + // the Hills memo values are before decimation, ratio above is after decimation
1885 +
1886 + // tolerance here is 5.0e-4, i.e. precision of values in Hills memo
1887 + double tolerance=5.0e-4;
1888 + if (std::abs(2.667-ratio)<=tolerance) {
1889 + numBin = 1;
1890 + } else if (std::abs(3.200-ratio*2.0)<=tolerance) {
1891 + numBin = 2;
1892 + } else if (std::abs(4.923-ratio*4.0)<=tolerance) {
1893 + numBin = 4;
1894 + } else if (std::abs(8.828-ratio*8.0)<=tolerance) {
1895 + numBin = 8;
1896 + } else if (std::abs(16.787-ratio*16.0)<tolerance) {
1897 + numBin=16;
1898 + }
1899 + }
1900 + }
1901 + }
1902 + }
1903 + }
1904 + if (numBin <= 0) {
1905 + // was not set to a valid value
1906 + // none of these cases generate an error or warning
1907 + // scalar columns don't exist - that's expected and so not an error
1908 + // unrecognzed telescope - that's also probably not an error, no rule for inferring numBin
1909 + // bad value for the VLA from that ratio, that might be an reportable error
1910 + // other windowFunction for ALMA - that's expected and so no an error, but might report if that's not UNIFORM
1911 + // no match within tolerance for ALMA - that might be a reportable error
1912 + numBin = 1;
1913 + }
1914 + return numBin;
1915 + }
1916 + // there's no way to get to here
1917 + throw ASDM2MSException("Impossible code line reached in getNumBin");
1918 +}
1919 +
1854 1920 /**
1855 1921 * This function fills the MS Spectral Window table.
1856 1922 * given :
1857 1923 * @parameter ds_p a pointer to the ASDM dataset.
1858 1924 */
1859 -void fillSpectralWindow(ASDM* ds_p, map<unsigned int, double>& effectiveBwPerSpwId_m) {
1925 +void fillSpectralWindow(ASDM* ds_p, map<unsigned int, double>& effectiveBwPerSpwId_m, const string &telescopeName) {
1860 1926 LOGENTER("fillSpectralWindow");
1861 1927
1862 1928 effectiveBwPerSpwId_m.clear();
1863 1929
1864 1930 try {
1865 1931 SpectralWindowTable& spwT = ds_p->getSpectralWindow();
1866 1932 vector<Tag> reorderedSwIds = reorderSwIds(*ds_p); // The vector of Spectral Window Tags in the order they will be inserted in the MS.
1867 1933 //for (vector<Tag>::size_type i = 0; i != reorderedSwIds.size() ; i++) cerr<<" reorderedSwIds["<<i<<"]="<<reorderedSwIds[i].getTagValue()<<endl;
1868 1934
1869 1935 for (vector<Tag>::size_type i = 0; i != reorderedSwIds.size() ; i++) swIdx2Idx[reorderedSwIds[i].getTagValue()] = i;
2085 2151 int numChan = r->getNumChan();
2086 2152 string name = r->isNameExists()?r->getName():"";
2087 2153 double refFreq = r->getRefFreq().get();
2088 2154 int measFreqRef = r->isMeasFreqRefExists()?FrequencyReferenceMapper::value(r->getMeasFreqRef()):MFrequency::TOPO;
2089 2155 double totalBandwidth = r->getTotBandwidth().get();
2090 2156 int netSideband = r->getNetSideband();
2091 2157 int bbcNo = r->getBasebandName();
2092 2158 int ifConvChain = 0;
2093 2159 int freqGroup = r->isFreqGroupExists()?r->getFreqGroup():0;
2094 2160 string freqGroupName = r->isFreqGroupNameExists()?r->getFreqGroupName().c_str():"";
2161 + // windowFunction is a required field
2162 + std::string windowFunction = CWindowFunction::name(r->getWindowFunction());
2163 + int numBin = getNumBin(r, telescopeName);
2164 + if (telescopeName == "EVLA" && (numBin>1) && (!r->isNumBinExists())) {
2165 + // numBin has been inferred for EVLA data, adjust resolution
2166 + resolution1D = chanWidth1D;
2167 + }
2095 2168
2096 2169 for (map<AtmPhaseCorrectionMod::AtmPhaseCorrection, ASDM2MSFiller*>::iterator iter = msFillers.begin();
2097 2170 iter != msFillers.end(); ++iter) {
2098 2171 iter->second->addSpectralWindow(numChan,
2099 2172 name,
2100 2173 refFreq,
2101 2174 chanFreq1D,
2102 2175 chanWidth1D,
2103 2176 measFreqRef,
2104 2177 effectiveBw1D,
2105 2178 resolution1D,
2106 2179 totalBandwidth,
2107 2180 netSideband,
2108 2181 bbcNo,
2109 2182 ifConvChain,
2110 2183 freqGroup,
2111 2184 freqGroupName,
2112 2185 numAssocValues,
2113 2186 assocSpectralWindowId_,
2114 - assocNature_);
2187 + assocNature_,
2188 + windowFunction,
2189 + numBin);
2115 2190 }
2116 2191 }
2117 2192 if (nSpectralWindow) {
2118 2193 infostream.str("");
2119 2194 infostream << "converted in " << msFillers.begin()->second->ms()->spectralWindow().nrow() << " spectral window(s) in the measurement set(s).";
2120 2195 info(infostream.str());
2121 2196 }
2122 2197 } catch (IllegalAccessException& e) {
2123 2198 errstream.str("");
2124 2199 errstream << e.getMessage();
5039 5114 } catch ( std::exception & e) {
5040 5115 errstream.str("");
5041 5116 errstream << e.what();
5042 5117 error(errstream.str());
5043 5118 }
5044 5119
5045 5120 //
5046 5121 // Process the SpectralWindow table.
5047 5122 //
5048 5123 map<unsigned int, double> effectiveBwPerSpwId_m;
5049 - fillSpectralWindow(ds, effectiveBwPerSpwId_m);
5124 + fillSpectralWindow(ds, effectiveBwPerSpwId_m, telName);
5050 5125
5051 5126 //
5052 5127 // Process the Polarization table
5053 5128 //
5054 5129 Stokes::StokesTypes linearCorr[] = { Stokes::XX, Stokes::XY, Stokes::YX, Stokes::YY };
5055 5130 Stokes::StokesTypes circularCorr[] = { Stokes::RR, Stokes::RL, Stokes::LR, Stokes::LL };
5056 5131 int corrProduct1[] = { 0, 0 };
5057 5132 int corrProduct2[] = { 0, 0, 1, 1};
5058 5133 int corrProduct4[] = { 0, 0, 0, 1, 1, 0, 1, 1 };
5059 5134

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut