| [2328] | 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 | 
 | 
|---|
 | 28 | void usage(void);
 | 
|---|
 | 29 | 
 | 
|---|
 | 30 | void usage(void) {
 | 
|---|
 | 31 |   cerr << "usage: gph424_1g xml_options_file" << endl;
 | 
|---|
 | 32 |   exit(-1);
 | 
|---|
 | 33 | }
 | 
|---|
 | 34 | 
 | 
|---|
 | 35 | // Prefs structure
 | 
|---|
 | 36 | 
 | 
|---|
 | 37 | struct 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 | 
 | 
|---|
 | 122 | DOM_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 | 
 | 
|---|
 | 146 | string 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 | 
 | 
|---|
 | 161 | bool DOMHasOption(DOM_Document& doc, string path) {
 | 
|---|
 | 162 |   DOM_Node n = DOMGetNode(doc, path);
 | 
|---|
 | 163 |   return n != 0;
 | 
|---|
 | 164 | }
 | 
|---|
 | 165 | 
 | 
|---|
 | 166 | void 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 | 
 | 
|---|
 | 183 | void 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 | 
 | 
|---|
 | 296 | int 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 | }
 | 
|---|