33 const char *
const fn =
"getError";
38 std::string tele0, tele;
39 std::vector<std::string> keys = {
"BAND",
"TELESCOP"};
40 hdu.keyWord(keys[0]).value(band0);
41 hdu.keyWord(keys[1]).value(tele0);
44 for (
auto &file: files) {
45 std::unique_ptr<CCfits::FITS>
46 ferr{
new CCfits::FITS(file, CCfits::Read,
false, keys)};
47 auto &phdu = ferr->pHDU();
48 phdu.keyWord(keys[0]).value(band);
49 phdu.keyWord(keys[1]).value(tele);
50 if (band == band0 && tele == tele0) {
53 "Error map '%s' size mismatch (%ud pixels, %ud in parent)",
54 file, err.size(), data.size());
56 MU_INFO(fn,
"Using error map '%s' for %ld micron band", file, band0);
62 double relerr = 0.01*cli.
getDouble(
'E', 3.0);
64 for (
auto &elem: err) { elem *= relerr; }
65 MU_INFO(fn,
"No error map found for %ld micron band, assuming %g%% of input",
86 const char *
const fn =
"process";
91 "Minimum of 2 input bands required (%lu given); try --help",
97 std::vector<std::unique_ptr<CCfits::FITS>> img;
98 std::vector<CCfits::PHDU*> bandHDU;
99 std::vector<mapDataType> bandData, bandError;
100 for (
auto &file: files) {
101 img.emplace_back(
new CCfits::FITS(file, CCfits::Read,
true));
102 bandHDU.push_back(&img.back()->pHDU());
106 auto &hd0 = *bandHDU.front(), &hd = *bandHDU.back();
107 int naxes = hd.axes();
108 MU_EXCEPTION_IF(naxes != 2, fn,
"%s has %ld axes (expected 2)", file,naxes);
109 MU_ERROR_IF(hd.axis(0) != hd0.axis(0) || hd.axis(1) != hd0.axis(1), fn,
110 "%s dimensions are %ldx%ld (fiducial is %ldx%ld)",
111 file, hd.axis(0), hd.axis(1), hd0.axis(0), hd0.axis(1));
115 std::string facility;
117 hd.read(bandData.back());
119 MU_INFO(fn,
"%s (%s): %ldx%ld @ %3ld microns (%s)", file,
120 hd.keyWord(
"OBJECT").value(name),
121 hd.axis(0), hd.axis(1), hd.keyWord(
"BAND").value(band),
122 hd.keyWord(
"TELESCOP").value(facility));
126 hd.readKey(
"BUNIT", units);
128 "%s map units are '%s' (need MJy/sr)", file, units);
132 for (
unsigned i=0; i != bandData.size(); i++) {
133 bandError.push_back(
getError(bandData[i], *bandHDU[i], cli));
137 auto mtfile = cli.
getString(
'o',
"mtcore.fits");
138 if (cli.
getOption(
'F')) { mtfile.insert(0, 1,
'!'); }
139 std::unique_ptr<CCfits::FITS>
140 outFITS(
new CCfits::FITS(mtfile, BYTE_IMG, 0,
nullptr));
143 std::vector<mapDataType> outMap;
144 std::vector<CCfits::ExtHDU*> outHDU;
145 solve(outMap, outHDU, outFITS.get(), bandData, bandError, bandHDU, cli);
148 for (
unsigned i=0, n=outMap.size(); i != n; i++) {
149 setKeys(*outHDU.at(i), *bandHDU[0]);
150 outHDU[i]->write(1, outMap[i].size(), outMap[i]);
155 for (
unsigned i=0, n=bandData.size(); i != n; i++) {
157 std::string bunit =
"BUNIT", units,
158 comm = bandHDU[i]->getComments(),
159 hist = bandHDU[i]->getHistory();
160 bandHDU[i]->readKey(bunit, units);
161 bandHDU[i]->keyWord(
"BAND").value(band);
162 std::string name =
"BAND" + std::to_string(band);
164 std::vector<long> axes = {bandHDU[i]->axis(0), bandHDU[i]->axis(1)};
165 auto eHDU = outFITS->addImage(name,
mapDataCode, axes);
166 eHDU->addKey(bunit, units,
"Image flux units");
167 eHDU->copyAllKeys(bandHDU[i]);
168 eHDU->writeComment(comm);
169 eHDU->writeHistory(hist);
170 eHDU->write(1, bandData[i].size(), bandData[i]);
174 eHDU->writeComment(
"ERROR map used for the associated parent image");
175 eHDU->addKey(bunit, units,
"Error flux units");
176 eHDU->copyAllKeys(bandHDU[i]);
177 eHDU->writeComment(comm);
178 eHDU->writeHistory(hist);
179 eHDU->write(1, bandError[i].size(), bandError[i]);
195 void setKeys(CCfits::ExtHDU &hdu,
const CCfits::PHDU &hdu0)
197 const char *
const fn =
"setKeys";
200 hdu.addKey(
"CREATOR", std::string(
"Manticore v") +
version,
201 "FITS file creator");
207 auto name = hdu.name().c_str();
208 if (name[0] ==
'N' || name[1] ==
'N') {
209 hdu.addKey<std::string>(
"BUNIT",
"cm^-2",
"Gas column density units");
211 if (name[0] ==
'T' || name[1] ==
'T') {
212 hdu.addKey<std::string>(
"BUNIT",
"K",
"Gas temperature units");
214 if (name[0] ==
'C') {
215 hdu.addKey<std::string>(
"BUNIT",
"chi2/DOF",
"Goodness of fit units");
219 auto &inKeys = hdu0.keyWord();
220 for (
auto &key: {
"OBJECT",
221 "CTYPE1",
"CRVAL1",
"CRPIX1",
"CDELT1",
222 "CTYPE2",
"CRVAL2",
"CRPIX2",
"CDELT2",
224 auto kp = inKeys.find(key);
225 if (kp != inKeys.end()) {
226 hdu.addKey(kp->second);
228 MU_EXCEPTION(fn,
"Fiducial image missing keyword '%s'", key);
void process(const mutils::CommandLine &cli)
Processes multi-wavelength maps into temperature/column density maps.
void setKeys(CCfits::ExtHDU &hdu, const CCfits::PHDU &hdu0)
Sets header keywords for an output map.
void solve(std::vector< mapDataType > &outMap, std::vector< CCfits::ExtHDU *> &outHDU, CCfits::FITS *outFITS, const std::vector< mapDataType > &bandData, const std::vector< mapDataType > &bandError, const std::vector< CCfits::PHDU *> &bandHDU, const mutils::CommandLine &cli)
Performs least-squares fitting of band SEDs.
unsigned getOption(const Option &option, Arguments *values=nullptr, unsigned nMax=-1) const
Retrieves all copies of a command line option.
const Arguments & arguments() const
Returns remaining arguments (read-only).
Command line options and arguments.
std::string getString(char shortName, const std::string &defValue) const
Retrieves unique string option.
std::vector< std::string > Arguments
Command line arguments type.
#define MU_EXCEPTION_IF(cond, src, msg,...)
Conditionally throws mutils::Exception.
mapDataType getError(const mapDataType &data, CCfits::PHDU &hdu, const mutils::CommandLine &cli)
Generates image error map.
constexpr const char *const version
Version string.
#define MU_EXCEPTION(src, msg,...)
Throws mutils::Exception using uniform messaging.
#define MU_ERROR_IF(cond, src, msg,...)
Conditionally log a message at level Log::ERROR.
#define MU_INFO(src, msg,...)
Log a message at level Log::INFO.
constexpr int mapDataCode
Output FITS maps data code (cf. mapDataType).
double getDouble(char shortName, double defValue) const
Retrieves unique floating point option.
std::valarray< float > mapDataType
Output FITS maps data type (cf. mapDataCode).