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

Last change on this file since 1819 was 1819, checked in by aubourg, 24 years ago

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