source: Sophya/trunk/ArchTOIPipe/TestPipes/quickmap_fio.cc@ 2461

Last change on this file since 2461 was 2328, checked in by aubourg, 23 years ago

pour piolib

File size: 13.2 KB
Line 
1// GPH 424.1 Planck HFI-L2 Simple Map Making
2// Eric Aubourg CEA/DAPNIA/SPP
3// This version is able to produce local maps as well
4// $Id: quickmap_fio.cc,v 1.1 2003-02-24 09:18:01 aubourg Exp $
5
6#include <stdlib.h>
7#include <stdio.h>
8#include "toi.h"
9#include "toiprocessor.h"
10#include "fitstoirdr.h"
11#include "fitstoiwtr.h"
12#include "toimanager.h"
13#include "toisegment.h"
14
15#include "sambainit.h"
16#include "toi2map.h"
17#include "fitsspherehealpix.h"
18#include "localmap.h"
19#include "fitslocalmap.h"
20
21#include <util/PlatformUtils.hpp>
22#include <sax/SAXException.hpp>
23#include <sax/SAXParseException.hpp>
24#include <parsers/DOMParser.hpp>
25#include <dom/DOM_DOMException.hpp>
26
27
28void usage(void);
29
30void usage(void) {
31 cerr << "usage: gph424_1g xml_options_file" << endl;
32 exit(-1);
33}
34
35// Prefs structure
36
37struct Gph424_1_Prefs {
38 struct {
39 struct {
40 string fname;
41 string flagfname;
42 vector<FlagToiDef> flags;
43 string coord1; // TOI name
44 string coord2; // TOI name
45 unsigned long typcoord;
46 } pointing;
47 struct {
48 string fname;
49 string flagfname;
50 vector<FlagToiDef> flags;
51 string toiname;
52 } signal;
53 } input;
54
55 struct {
56 string fname;
57 string wfname; // weight
58 unsigned long typcoord;
59 string typmap; // local, healpix...
60 struct {
61 int nlat;
62 } healpix_opt;
63 struct {
64 int npx, npy;
65 double org1,org2;
66 double anglex, angley;
67 } local_opt;
68 } output;
69
70 int nproc;
71 long snb, sne;
72
73 void dump() {
74 cout << "*** GPH424-1 settings\n";
75 cout << " input\n";
76 cout << " | pointing\n";
77 cout << " | | fname " << input.pointing.fname << "\n";
78 cout << " | | flagfname " << input.pointing.flagfname << "\n";
79 for (int i=0; i<input.pointing.flags.size(); i++) {
80 cout << " | | | flag " << i << " : " << input.pointing.flags[i] << "\n";
81 }
82 cout << " | | coord1 " << input.pointing.coord1 << "\n";
83 cout << " | | coord2 " << input.pointing.coord2 << "\n";
84 cout << " | | typcoord "
85 << ((input.pointing.typcoord & TypCoordGal) ? "G" : "E")
86 << ((input.pointing.typcoord & TypCoord1H) ? "H" :
87 (input.pointing.typcoord & TypCoord1D) ? "D" : "R")
88 << ((input.pointing.typcoord & TypCoord1C) ? "C" : "L")
89 << ((input.pointing.typcoord & TypCoord2H) ? "H" :
90 (input.pointing.typcoord & TypCoord2D) ? "D" : "R")
91 << ((input.pointing.typcoord & TypCoord2C) ? "C" : "L")
92 << "\n";
93 cout << " | signal\n";
94 cout << " | | fname " << input.signal.fname << "\n";
95 cout << " | | flagfname " << input.signal.flagfname << "\n";
96 {for (int i=0; i<input.signal.flags.size(); i++) {
97 cout << " | | | flag " << i << " : " << input.signal.flags[i] << "\n";
98 }}
99 cout << " | | toiname " << input.signal.toiname << "\n";
100 cout << " output\n";
101 cout << " | fname " << output.fname << "\n";
102 cout << " | wfname " << output.wfname << "\n";
103 cout << " | typcoord "
104 << ((output.typcoord & TypCoordGal) ? "G" : "E")
105 << "\n";
106 cout << " | typmap " << output.typmap << "\n";
107 if (output.typmap == "healpix") {
108 cout << " | | nlat " << output.healpix_opt.nlat << "\n";
109 } else if (output.typmap == "local") {
110 cout << " | | npix " << output.local_opt.npx << ", " << output.local_opt.npy << "\n";
111 cout << " | | origin " << output.local_opt.org1 << ", " << output.local_opt.org2 << "\n";
112 cout << " | | extension " << output.local_opt.anglex << ", " << output.local_opt.angley << "\n";
113 }
114 cout << " nproc " << nproc << "\n";
115 if (sne >= 0 && snb >= 0) {
116 cout << " samples " << snb << " " << sne << "\n";
117 }
118 cout << "***" << endl;
119 }
120};
121
122DOM_Node DOMGetNode(DOM_Document& doc, string path) {
123 DOM_Element el = doc.getDocumentElement();
124
125 while (path != "") {
126 if (path[0] == '/') {
127 path = path.substr(1);
128 }
129 if (path == "") break;
130 int x = path.find('/');
131 string pathElem;
132 if (x == -1) {
133 pathElem = path;
134 path = "";
135 } else {
136 pathElem = path.substr(0, x);
137 path = path.substr(x);
138 }
139 DOM_Node n = el.getElementsByTagName(pathElem.c_str()).item(0);
140 if (n==0) return n;
141 el = (DOM_Element&) n;
142 }
143 return el;
144}
145
146string DOMGetString(DOM_Document& doc, string path, string def="") {
147 DOM_Node n = DOMGetNode(doc, path);
148 if (n==0) return def;
149 char* s = n.getFirstChild().getNodeValue().transcode();
150 char* p=s; while ((*p)==' ') p++;
151 char* q=p+strlen(p)-1;
152 while (q>p && *q==' ') {
153 *q = '\0';
154 q--;
155 }
156 string ss = p;
157 free(s);
158 return ss;
159}
160
161bool DOMHasOption(DOM_Document& doc, string path) {
162 DOM_Node n = DOMGetNode(doc, path);
163 return n != 0;
164}
165
166void parseFlags(DOM_Node flagnode, vector<FlagToiDef>& flgs) {
167 DOM_Element elt = (DOM_Element&) flagnode;
168 DOM_NodeList l = elt.getElementsByTagName("flag");
169 // prepare array
170 flgs.clear();
171 flgs.insert(flgs.begin(), l.getLength(), (FlagToiDef)0);
172 for (int i = 0; i<l.getLength(); i++) {
173 DOM_Node node = l.item(i);
174 DOM_Element elt = (DOM_Element&) node;
175 DOMString scol = elt.getAttribute("column");
176 DOMString sval = elt.getAttribute("kind");
177 int col = atoi(scol.transcode());
178 int value = atoi(sval.transcode());
179 flgs[col-1] = (FlagToiDef) value;
180 }
181}
182
183void parsePrefs(string prefFile, Gph424_1_Prefs& prefs) {
184 try {
185 XMLPlatformUtils::Initialize();
186 } catch (const XMLException& toCatch) {
187 cerr << "Error during XML initialization! :\n"
188 << toCatch.getMessage() << endl;
189 exit(-1);
190 }
191
192 DOMParser* parser = new DOMParser;
193
194 try {
195 parser->parse(prefFile.c_str());
196 } catch (const XMLException& toCatch) {
197 cerr << "\nError during parsing: '" << prefFile << "'\n"
198 << "Exception message is: \n"
199 << toCatch.getMessage() << "\n" << endl;
200 exit(-1);
201 } catch (const DOM_DOMException& toCatch) {
202 cerr << "\nDOM Error during parsing: '" << prefFile << "'\n"
203 << "DOMException code is: \n"
204 << toCatch.code << "\n" << endl;
205 exit(-1);
206 } catch (...) {
207 cerr << "\nUnexpected exception during parsing: '" << prefFile << "'\n";
208 exit(-1);
209 }
210
211 DOM_Document doc = parser->getDocument();
212
213 prefs.input.pointing.fname = DOMGetString(doc, "/input/pointing/fname");
214 if (DOMHasOption(doc, "/input/pointing/flags")) {
215 prefs.input.pointing.flagfname =
216 DOMGetString(doc, "/input/pointing/flags/fname");
217 DOM_Node flagnode = DOMGetNode(doc, "/input/pointing/flags");
218 parseFlags(flagnode, prefs.input.pointing.flags);
219 }
220 prefs.input.pointing.coord1 = DOMGetString(doc, "/input/pointing/coord1","phi");
221 prefs.input.pointing.coord2 = DOMGetString(doc, "/input/pointing/coord2","theta");
222 bool inequ = DOMHasOption(doc, "/input/pointing/equatorial");
223 bool ingal = DOMHasOption(doc, "/input/pointing/galactic");
224 prefs.input.pointing.typcoord = ingal ? TypCoordGal : TypCoordEq;
225 bool inhd = DOMHasOption(doc, "/input/pointing/hourdeg");
226 bool indd = DOMHasOption(doc, "/input/pointing/degdeg");
227 bool inrr = DOMHasOption(doc, "/input/pointing/radian");
228 bool iscolat = DOMHasOption(doc, "/input/pointing/colat");
229 if (inhd) {
230 prefs.input.pointing.typcoord |= TypCoord1H|TypCoord2D;
231 } else if (indd) {
232 prefs.input.pointing.typcoord |= TypCoord1D|TypCoord2D;
233 } else if (inrr) {
234 prefs.input.pointing.typcoord |= TypCoord1R|TypCoord2R;
235 } else if (prefs.input.pointing.typcoord & TypCoordEq) {
236 prefs.input.pointing.typcoord = TypCoordEqStd;
237 } else {
238 prefs.input.pointing.typcoord = TypCoordGalStd;
239 }
240 if (iscolat) {
241 prefs.input.pointing.typcoord |= TypCoord2C;
242 }
243
244 prefs.input.signal.fname = DOMGetString(doc, "/input/signal/fname");
245 prefs.input.signal.toiname = DOMGetString(doc, "/input/signal/toiname");
246 if (DOMHasOption(doc, "/input/signal/flags")) {
247 prefs.input.signal.flagfname =
248 DOMGetString(doc, "/input/signal/flags/fname");
249 DOM_Node flagnode = DOMGetNode(doc, "/input/signal/flags");
250 parseFlags(flagnode, prefs.input.signal.flags);
251 }
252 prefs.output.fname = DOMGetString(doc, "/output/fname");
253 prefs.output.wfname = DOMGetString(doc, "/output/wfname");
254 bool outequ = DOMHasOption(doc, "/output/equatorial");
255 bool outgal = DOMHasOption(doc, "/output/galactic");
256 prefs.output.typcoord = outgal ? TypCoordGal : TypCoordEq;
257 bool outhd = DOMHasOption(doc, "/output/hourdeg");
258 bool outdd = DOMHasOption(doc, "/output/degdeg");
259 bool outrr = DOMHasOption(doc, "/output/radian");
260 if (outhd) {
261 prefs.output.typcoord |= TypCoord1H | TypCoord2D;
262 } else if (outdd) {
263 prefs.output.typcoord |= TypCoord1D | TypCoord2D;
264 } else if (outrr) {
265 prefs.output.typcoord |= TypCoord1R | TypCoord2R;
266 } else if (prefs.output.typcoord & TypCoordEq) {
267 prefs.output.typcoord = TypCoordEqStd;
268 } else {
269 prefs.output.typcoord |= TypCoordGalStd;
270 }
271
272 prefs.output.typmap = DOMHasOption(doc, "/output/local") ? "local" :
273 (DOMHasOption(doc, "/output/healpix") ? "healpix" : "healpix");
274 prefs.output.healpix_opt.nlat = atoi(DOMGetString(doc, "/output/healpix/nlat", "128").c_str());
275
276 string s = DOMGetString(doc, "/output/local/npx");
277 prefs.output.local_opt.npx = prefs.output.local_opt.npy = -1;
278 int nc = sscanf(s.c_str(), "%d %d", &prefs.output.local_opt.npx, &prefs.output.local_opt.npy);
279 if (nc == 1) prefs.output.local_opt.npy = prefs.output.local_opt.npx;
280
281 s = DOMGetString(doc, "/output/local/center");
282 prefs.output.local_opt.org1 = prefs.output.local_opt.org2 = 0;
283 sscanf(s.c_str(), "%lg %lg", &prefs.output.local_opt.org1, &prefs.output.local_opt.org2);
284
285 s = DOMGetString(doc, "/output/local/angle");
286 prefs.output.local_opt.anglex = prefs.output.local_opt.angley = 0;
287 nc = sscanf(s.c_str(), "%d %d", &prefs.output.local_opt.anglex, &prefs.output.local_opt.angley);
288 if (nc == 1) prefs.output.local_opt.angley = prefs.output.local_opt.anglex;
289
290 prefs.nproc = atoi(DOMGetString(doc, "/nproc", "1").c_str());
291
292 s = DOMGetString(doc, "/samples", "-1 -1");
293 sscanf(s.c_str(), "%ld %ld", &prefs.snb, &prefs.sne);
294}
295
296int main(int argc, char** argv) {
297
298 // Read prefs
299 Gph424_1_Prefs prefs;
300 parsePrefs(argv[1], prefs);
301 prefs.dump();
302
303 // Assemble pipeline
304
305 TOIManager* mgr = TOIManager::getManager();
306
307 // Two input files
308
309 FITSTOIReader inPoint(prefs.input.pointing.fname);
310 FITSTOIReader inSig (prefs.input.signal.fname);
311
312 inPoint.setImplicitSN();
313 inSig.setImplicitSN();
314
315 // Associate flag files
316
317 if (prefs.input.pointing.flagfname != "") {
318 inPoint.setFlagFile(prefs.input.pointing.flagfname, prefs.input.pointing.flags);
319 }
320
321 if (prefs.input.signal.flagfname != "") {
322 inSig.setFlagFile(prefs.input.signal.flagfname, prefs.input.signal.flags);
323 }
324
325 // Prepare maps for output
326
327 PixelMap<r_8>* map;
328 PixelMap<r_8>* wmap;
329
330 if (prefs.output.typmap == "local") {
331 LocalMap<r_8>* lmap = new LocalMap<r_8>(prefs.output.local_opt.npx, prefs.output.local_opt.npy);
332 lmap->SetOrigin(prefs.output.local_opt.org1, prefs.output.local_opt.org2);
333 lmap->SetSize(prefs.output.local_opt.anglex, prefs.output.local_opt.angley);
334 map = lmap;
335 wmap = new LocalMap<r_8>(*lmap);
336 } else if (prefs.output.typmap == "healpix") {
337 map = new SphereHEALPix<r_8>(prefs.output.healpix_opt.nlat);
338 wmap = new SphereHEALPix<r_8>(prefs.output.healpix_opt.nlat);
339 } else {
340 cout << "typmap " << prefs.output.typmap << " unknown" << endl;
341 exit(-1);
342 }
343
344 TOI2Map toi2m(map,wmap);
345
346 toi2m.SetEquinox(2000.);
347 toi2m.SetCoorIn((TypAstroCoord) prefs.input.pointing.typcoord);
348 toi2m.SetCoorMap((TypAstroCoord) prefs.output.typcoord);
349
350 TOISegmented * toicoord1in = new TOISegmented(prefs.input.pointing.coord1);
351 inPoint.addOutput(prefs.input.pointing.coord1,toicoord1in);
352 toi2m.addInput("Coord1In",toicoord1in);
353
354 TOISegmented * toicoord2in = new TOISegmented(prefs.input.pointing.coord2);
355 inPoint.addOutput(prefs.input.pointing.coord2,toicoord2in);
356 toi2m.addInput("Coord2In",toicoord2in);
357
358 TOISegmented * toibolin = new TOISegmented("toi_bolo_in");
359 inSig.addOutput(prefs.input.signal.toiname,toibolin);
360 toi2m.addInput("BoloIn",toibolin);
361
362 // Start processing
363
364 if (prefs.snb >= 0 && prefs.sne >= 0) {
365 toi2m.setRequestedSample(prefs.snb, prefs.sne);
366 }
367
368 mgr->startAll();
369 mgr->waitForAll();
370
371 // Save maps
372
373 {
374 FitsOutFile sfits(prefs.output.fname,FitsFile::clear);
375 cout<<"Creating map fits file "<<prefs.output.fname<<endl;
376 // Would need some generic mechanism here $CHECK$ $$TBD$$
377 if (map->TypeOfMap() == "RING") {
378 sfits << * (SphereHEALPix<r_8>*)map;
379 } else if (map->TypeOfMap() == "LOCAL") {
380 sfits << * (LocalMap<r_8>*)map;
381 } else {
382 cout << "error, unknown type " << map->TypeOfMap() << endl;
383 }
384 }
385
386 {
387 FitsOutFile swfits(prefs.output.wfname,FitsFile::clear);
388 cout<<"Creating map weight fits file "<<prefs.output.wfname<<endl;
389 // swfits << *wmap;
390 if (wmap->TypeOfMap() == "RING") {
391 swfits << * (SphereHEALPix<r_8>*)wmap;
392 } else if (wmap->TypeOfMap() == "LOCAL") {
393 swfits << * (LocalMap<r_8>*)wmap;
394 } else {
395 cout << "error, unknown type " << wmap->TypeOfMap() << endl;
396 }
397 }
398
399 delete map;
400 delete wmap;
401
402 return 0;
403}
Note: See TracBrowser for help on using the repository browser.