Line data Source code
1 : //# MSContinuumSubtractor.cc: Subtract continuum from spectral line data
2 : //# Copyright (C) 2004
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 : //# $Id$
27 : //#
28 :
29 : #include <spectrallines/Splatalogue/ListConverter.h>
30 :
31 : #include <casacore/casa/Quanta/Quantum.h>
32 : #include <casacore/casa/Quanta/Unit.h>
33 : #include <casacore/casa/OS/RegularFile.h>
34 : #include <casacore/casa/IO/RegularFileIO.h>
35 : #include <casacore/casa/Utilities/Regex.h>
36 : #include <casacore/tables/Tables/ScaColDesc.h>
37 : #include <casacore/tables/Tables/ScalarColumn.h>
38 : #include <casacore/tables/Tables/SetupNewTab.h>
39 : #include <casacore/tables/DataMan/StandardStMan.h>
40 : #include <casacore/tables/Tables/TableDesc.h>
41 :
42 : using namespace casacore;
43 : namespace casa {
44 :
45 2 : ListConverter::ListConverter(const Vector<String>& filenames, const String& table) :
46 2 : _log(new LogIO()), _listFiles(Vector<RegularFile>(0)), _tableName(table), _freqUnit("") {
47 2 : RegularFile t(_tableName);
48 2 : if (t.exists()) {
49 0 : *_log << "File " << _tableName << " already exists and will not be overwritten" << LogIO::EXCEPTION;
50 : }
51 2 : File::FileWriteStatus tableStatus = t.getWriteStatus();
52 2 : if(tableStatus == File::NOT_CREATABLE ) {
53 0 : *_log << "Table " << _tableName << " cannot be created" << LogIO::EXCEPTION;
54 : }
55 4 : for (Vector<String>::const_iterator iter = filenames.begin(); iter != filenames.end(); iter++) {
56 2 : RegularFile file(*iter);
57 2 : if (! file.exists()) {
58 0 : *_log << "File " << *iter << " does not exist" << LogIO::EXCEPTION;
59 : }
60 2 : if (! file.isReadable()) {
61 0 : *_log << "File " << *iter << " does not have read permission" << LogIO::EXCEPTION;
62 : }
63 2 : _listFiles.resize(_listFiles.size()+1, true);
64 2 : _listFiles[_listFiles.size() - 1] = file;
65 4 : }
66 2 : }
67 :
68 2 : ListConverter::~ListConverter() {
69 2 : delete _log;
70 2 : }
71 :
72 2 : SplatalogueTable* ListConverter::load() {
73 2 : _parseLists();
74 2 : SplatalogueTable *newTable = _defineTable(_species.size());
75 2 : _addData(newTable);
76 2 : newTable->flush();
77 2 : return newTable;
78 : }
79 :
80 2 : void ListConverter::_parseLists() {
81 2 : _log->origin(LogOrigin("ListConverter", __FUNCTION__));
82 :
83 2 : uInt transitionIndex = 0;
84 2 : Regex blankLine("^[ \n\t\r\v\f]+$");
85 2 : Regex lastLine("^Found [0-9]+ rows");
86 2 : Regex debye2("D<sup>2</sup>");
87 2 : Regex k("(K)");
88 2 : Regex cm1("(cm<sup>-1</sup>)");
89 :
90 2 : _frequency.resize(0);
91 2 : for (
92 4 : Vector<RegularFile>::const_iterator iter = _listFiles.begin();
93 4 : iter != _listFiles.end(); iter++
94 : ) {
95 2 : RegularFileIO fileIO(*iter);
96 2 : Int bufSize = 10000;
97 2 : char *buffer = new char[bufSize];
98 : int nRead;
99 2 : String contents;
100 644 : while ((nRead = fileIO.read(bufSize, buffer, false)) == bufSize) {
101 642 : String chunk(buffer, bufSize);
102 642 : contents += chunk;
103 642 : }
104 : // get the last chunk
105 2 : String chunk(buffer, nRead);
106 2 : contents += chunk;
107 :
108 2 : Vector<String> lines = stringToVector(contents, '\n');
109 2 : uInt newSize = _frequency.size() + lines.size();
110 2 : _species.resize(newSize, true);
111 2 : _chemName.resize(newSize, true);
112 2 : _qns.resize(newSize, true);
113 2 : _lineList.resize(newSize, true);
114 2 : _recommended.resize(newSize, true);
115 2 : _frequency.resize(newSize, true);
116 2 : _intensity.resize(newSize, true);
117 2 : _smu2.resize(newSize, true);
118 2 : _logA.resize(newSize, true);
119 2 : _eL.resize(newSize, true);
120 2 : _eU.resize(newSize, true);
121 2 : uInt lineCount = 0;
122 60004 : for(Vector<String>::iterator liter=lines.begin(); liter!=lines.end(); liter++) {
123 60006 : if (
124 120002 : liter->empty() || liter->firstchar() == '#'
125 120002 : || liter->matches(blankLine) || liter->matches(lastLine)
126 : ) {
127 : // skip comment
128 4 : lineCount++;
129 6 : continue;
130 : }
131 59998 : uInt tabCount = liter->freq('\t');
132 59998 : if (tabCount != 10) {
133 0 : *_log << "bad format for line " << (lineCount+1)
134 0 : << " in file " << iter->path().baseName() << LogIO::EXCEPTION;
135 : }
136 59998 : Vector<String> parts = stringToVector(*liter, '\t');
137 719976 : for (Vector<String>::iterator viter = parts.begin(); viter != parts.end(); viter++) {
138 659978 : viter->trim();
139 59998 : }
140 119996 : String filename = iter->path().dirName() + "/" + iter->path().baseName();
141 59998 : String species = parts[0];
142 59998 : String recommended = parts[1];
143 59998 : String chemicalName = parts[2];
144 59998 : String frequency = parts[3];
145 59998 : String qns = parts[4];
146 59998 : String intensity = parts[5];
147 59998 : String smu2 = parts[6];
148 59998 : String logA = parts[7];
149 59998 : String eL = parts[8];
150 59998 : String eU = parts[9];
151 59998 : String lineList = parts[10];
152 59998 : if (lineCount == 0) {
153 : // header line look for units
154 2 : Vector<String> freqHeaderParts = stringToVector(frequency, ' ');
155 6 : for (
156 4 : Vector<String>::const_iterator viter=freqHeaderParts.begin();
157 8 : viter!=freqHeaderParts.end(); viter++
158 : ) {
159 : try {
160 6 : Unit freqUnit(*viter);
161 4 : Quantity qq(1, freqUnit);
162 4 : if (qq.isConform(Unit("Hz"))) {
163 2 : _freqUnit = freqUnit.getName();
164 : }
165 4 : }
166 2 : catch (AipsError x) {}
167 2 : }
168 2 : Vector<String> smuHeaderParts = stringToVector(smu2, ' ');
169 4 : for (
170 4 : Vector<String>::const_iterator viter=smuHeaderParts.begin();
171 6 : viter!=smuHeaderParts.end(); viter++
172 : ) {
173 : try {
174 4 : String smu2(*viter);
175 4 : if (smu2.matches(debye2)) {
176 0 : _smu2Unit = "Debye*Debye";
177 :
178 : }
179 4 : }
180 0 : catch (AipsError x) {}
181 2 : }
182 2 : Vector<String> elHeaderParts = stringToVector(eL, ' ');
183 4 : for (
184 4 : Vector<String>::const_iterator viter=elHeaderParts.begin();
185 6 : viter!=elHeaderParts.end(); viter++
186 : ) {
187 : try {
188 4 : String el = *viter;
189 4 : if (el.matches(k)) {
190 0 : _elUnit = "K";
191 : }
192 4 : else if (el.matches(cm1)) {
193 0 : _elUnit = "cm-1";
194 : }
195 4 : }
196 0 : catch (AipsError x) {}
197 2 : }
198 2 : Vector<String> euHeaderParts = stringToVector(eU, ' ');
199 4 : for (
200 4 : Vector<String>::const_iterator viter=euHeaderParts.begin();
201 6 : viter!=euHeaderParts.end(); viter++
202 : ) {
203 : try {
204 4 : String eu = *viter;
205 4 : if (eu.matches(k)) {
206 0 : _euUnit = "K";
207 : }
208 4 : else if (eu.matches(cm1)) {
209 0 : _euUnit = "cm-1";
210 : }
211 4 : }
212 0 : catch (AipsError x) {}
213 2 : }
214 2 : lineCount++;
215 2 : continue;
216 2 : }
217 :
218 59996 : if (lineCount % 20000 == 0) {
219 2 : *_log << LogIO::NORMAL << "Parsing line " << lineCount << " of file "
220 2 : << iter->path().baseName() << LogIO::POST;
221 : }
222 59996 : Bool recFreq = recommended.length() > 0;
223 59996 : if (! frequency.matches(RXdouble) ) {
224 0 : *_log << "File " << filename << ", line " << lineCount
225 : << ": frequency value " << frequency << " is not numeric"
226 0 : << LogIO::EXCEPTION;
227 : }
228 59996 : if (! intensity.matches(RXdouble) ) {
229 0 : *_log << "File " << filename << ", line " << lineCount
230 : << ": intensity value " << intensity << " is not numeric"
231 0 : << LogIO::EXCEPTION;
232 : }
233 59996 : if (! smu2.matches(RXdouble) ) {
234 0 : *_log << "File " << filename << ", line " << lineCount
235 : << ": S_ij*mu**2 value " << smu2 << " is not numeric"
236 0 : << LogIO::EXCEPTION;
237 : }
238 59996 : if (! logA.matches(RXdouble) ) {
239 0 : *_log << "File " << filename << ", line " << lineCount
240 : << ": log(A) value " << logA << " is not numeric"
241 0 : << LogIO::EXCEPTION;
242 : }
243 59996 : if (! eL.matches(RXdouble) ) {
244 0 : *_log << "File " << filename << ", line " << lineCount
245 : << ": E_l value " << eL << " is not numeric"
246 0 : << LogIO::EXCEPTION;
247 : }
248 59996 : if (! eU.matches(RXdouble) ) {
249 0 : *_log << "File " << filename << ", line " << lineCount
250 : << ": E_u value " << eU << " is not numeric"
251 0 : << LogIO::EXCEPTION;
252 : }
253 59996 : _species[transitionIndex] = species;
254 59996 : _chemName[transitionIndex] = chemicalName;
255 59996 : _qns[transitionIndex] = qns;
256 59996 : _lineList[transitionIndex] = lineList;
257 59996 : _recommended[transitionIndex] = recFreq;
258 59996 : _frequency[transitionIndex] = String::toDouble(frequency);
259 59996 : _intensity[transitionIndex] = String::toFloat(intensity);
260 59996 : _smu2[transitionIndex] = String::toFloat(smu2);
261 59996 : _logA[transitionIndex] = String::toFloat(logA);
262 59996 : _eL[transitionIndex] = String::toFloat(eL);
263 59996 : _eU[transitionIndex] = String::toFloat(eU);
264 59996 : transitionIndex++;
265 59996 : lineCount++;
266 60024 : }
267 4 : }
268 2 : _species.resize(transitionIndex, true);
269 2 : _chemName.resize(transitionIndex, true);
270 2 : _qns.resize(transitionIndex, true);
271 2 : _lineList.resize(transitionIndex, true);
272 2 : _recommended.resize(transitionIndex, true);
273 2 : _frequency.resize(transitionIndex, true);
274 2 : _intensity.resize(transitionIndex, true);
275 2 : _smu2.resize(transitionIndex, true);
276 2 : _logA.resize(transitionIndex, true);
277 2 : _eL.resize(transitionIndex, true);
278 2 : _eU.resize(transitionIndex, true);
279 2 : }
280 :
281 :
282 2 : SplatalogueTable* ListConverter::_defineTable(const uInt nrows) {
283 2 : String descName = _tableName + "_desc";
284 2 : TableDesc td(descName,TableDesc::Scratch);
285 2 : td.addColumn (ScalarColumnDesc<String>(SplatalogueTable::SPECIES, SplatalogueTable::SPECIES));
286 2 : td.addColumn (ScalarColumnDesc<Bool>(SplatalogueTable::RECOMMENDED, SplatalogueTable::RECOMMENDED));
287 2 : td.addColumn (ScalarColumnDesc<String>(SplatalogueTable::CHEMICAL_NAME, SplatalogueTable::CHEMICAL_NAME));
288 2 : td.addColumn (ScalarColumnDesc<Double>(SplatalogueTable::FREQUENCY, SplatalogueTable::FREQUENCY));
289 2 : td.addColumn (ScalarColumnDesc<String>(SplatalogueTable::QUANTUM_NUMBERS, SplatalogueTable::QUANTUM_NUMBERS));
290 2 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::INTENSITY, SplatalogueTable::INTENSITY));
291 2 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::SMU2, SplatalogueTable::SMU2));
292 2 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::LOGA, SplatalogueTable::LOGA));
293 2 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::EL, SplatalogueTable::EL));
294 2 : td.addColumn (ScalarColumnDesc<Float>(SplatalogueTable::EU, SplatalogueTable::EU));
295 2 : td.addColumn (ScalarColumnDesc<String>(SplatalogueTable::LINELIST, SplatalogueTable::LINELIST));
296 2 : SetupNewTable tableSetup(_tableName, td, Table::New);
297 2 : StandardStMan stm;
298 2 : tableSetup.bindColumn(SplatalogueTable::SPECIES, stm);
299 2 : tableSetup.bindColumn(SplatalogueTable::RECOMMENDED, stm);
300 2 : tableSetup.bindColumn(SplatalogueTable::CHEMICAL_NAME, stm);
301 2 : tableSetup.bindColumn(SplatalogueTable::FREQUENCY, stm);
302 2 : tableSetup.bindColumn(SplatalogueTable::QUANTUM_NUMBERS, stm);
303 2 : tableSetup.bindColumn(SplatalogueTable::INTENSITY, stm);
304 2 : tableSetup.bindColumn(SplatalogueTable::SMU2, stm);
305 2 : tableSetup.bindColumn(SplatalogueTable::LOGA, stm);
306 2 : tableSetup.bindColumn(SplatalogueTable::EL, stm);
307 2 : tableSetup.bindColumn(SplatalogueTable::EU, stm);
308 2 : tableSetup.bindColumn(SplatalogueTable::LINELIST, stm);
309 2 : if (_freqUnit.empty()) {
310 0 : *_log << LogIO::WARN << "Could not determine frequency unit, assuming GHz" << LogIO::POST;
311 0 : _freqUnit = "GHz";
312 : }
313 4 : return new SplatalogueTable(tableSetup, nrows, _freqUnit, _smu2Unit, _elUnit, _euUnit);
314 2 : }
315 :
316 2 : void ListConverter::_addData(const SplatalogueTable* table) const {
317 2 : ScalarColumn<String> col(*table, SplatalogueTable::SPECIES);
318 2 : col.putColumn(_species);
319 2 : ScalarColumn<Bool> rcol(*table, SplatalogueTable::RECOMMENDED);
320 2 : rcol.putColumn(_recommended);
321 2 : ScalarColumn<String> cncol(*table, SplatalogueTable::CHEMICAL_NAME);
322 2 : cncol.putColumn(_chemName);
323 2 : ScalarColumn<Double> fcol(*table,SplatalogueTable::FREQUENCY);
324 2 : fcol.putColumn(_frequency);
325 2 : ScalarColumn<String> qncol(*table, SplatalogueTable::QUANTUM_NUMBERS);
326 2 : qncol.putColumn(_qns);
327 2 : ScalarColumn<Float> icol(*table, SplatalogueTable::INTENSITY);
328 2 : icol.putColumn(_intensity);
329 2 : ScalarColumn<Float> smu2col(*table, SplatalogueTable::SMU2);
330 2 : smu2col.putColumn(_smu2);
331 2 : ScalarColumn<Float> lacol(*table, SplatalogueTable::LOGA);
332 2 : lacol.putColumn(_logA);
333 2 : ScalarColumn<Float> elcol(*table, SplatalogueTable::EL);
334 2 : elcol.putColumn(_eL);
335 2 : ScalarColumn<Float> eucol(*table, SplatalogueTable::EU);
336 2 : eucol.putColumn(_eU);
337 2 : ScalarColumn<String> llcol(*table, SplatalogueTable::LINELIST);
338 2 : llcol.putColumn(_lineList);
339 2 : }
340 :
341 : } //# NAMESPACE CASA - END
342 :
|
