source: Sophya/trunk/ArchTOIPipe/TestPipes/quickmap.cc@ 3407

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

pour piolib

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