source: Sophya/trunk/SophyaProg/PMixer/skymixer.cc@ 1356

Last change on this file since 1356 was 1355, checked in by sophie, 25 years ago

sophie: changes according to Garching requirements
File size: 31.6 KB
Line 
1#include "pmixer.h"
2
3/*!
4 * \defgroup PMixer PMixer module
5 * This module contains programs which:
6 * <UL>
7 * <LI> add several sky components, taking into account their
8 * radiation spectra and convoluting them with a given filter
9 * response : skymixer
10 * <LI> create a map with point source : extractRS
11 * <LI> generate sky components, radiation spectra and spectral
12 * response (small generator of maps) : tgsky and tgrsr
13 * </UL>
14 * A detailed description may be found at:
15 */
16/*!
17 * \ingroup PMixer
18 * \file skymixer.cc
19 *\brief \b PROGRAM \b skyMixer <BR>
20 * add several sky components, taking into account their
21 *radiation spectra and convoluting them with a given filter
22 *response
23 */
24
25// -----------------------------------------------------------------
26// ------------- Function declaration ------------------------------
27int CheckCards(DataCards & dc, string & msg, FitsOutFile& fios);
28char * BuildFITSFileName(string const & fname);
29SpectralResponse * getSpectralResponse(DataCards & dc, FitsOutFile& fios);
30RadSpectra * getEmissionSpectra(DataCards & dc, int nc, FitsOutFile& fios);
31void SKM_MergeFITSKeywords(char * flnm);
32void SKM_MergeFITSKeywords2(DVList & dvl);
33void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name);
34void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name);
35
36// to add different sky components and corresponding tools
37//----------------------------------------------------------
38template <class T>
39void addComponent(SpectralResponse& sr,
40 PixelMap<T>& finalMap,
41 PixelMap<T>& mapToAdd,
42 RadSpectra& rs, double K=1.);
43//
44template <class T>
45void addComponentBeta(SphereHEALPix<T>& finalMap,
46 SphereHEALPix<T>& mapToAdd,SpectralResponse& sr,
47 SphereHEALPix<T>& betaMap, double normFreq, double K);
48//
49template <class T>
50void integratedMap(SpectralResponse& sr,
51 SphereHEALPix<T>& betaMap, double normFreq, SphereHEALPix<T>& intBetaMap);
52
53//
54template <class T>
55void addComponentBeta(SphereHEALPix<T>& finalMap,
56 SphereHEALPix<T>& mapToAdd,
57 SphereHEALPix<T>& intBetaMap, double K);
58//
59template <class T>
60void addDipole(SpectralResponse& sr, PixelMap<T>& finalMap,
61 double theta,double phi,double amp,double temp);
62//
63// -----------------------------------------------------------------
64
65// ----- Global (static) variables ------------
66static bool rdmap = false; // true -> Read map first
67static char mapPath[256]; // Path for input maps
68static int hp_nside = 32; // HealPix NSide
69static int nskycomp = 0; // Number of sky components
70static int debuglev = 0; // Debug Level
71static int printlev = 0; // Print Level
72static POutPersist * so = NULL; // Debug PPFOut file
73static DVList * dvl_fitskw = NULL; // Global DVList for all FITS Keywords
74bool checkpdmtype=false;
75bool dvliston= false;
76// --------- SkyMixer Version --------------
77static double skm_version = 1.5;
78
79// -------------------------------------------------------------------------
80// main program
81// -------------------------------------------------------------------------
82int main(int narg, char * arg[])
83{
84 if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
85 cout << " Usage: skymixer parameterFile outputfitsname [outppfname]" << endl;
86 exit(0);
87 }
88
89 InitTim();
90 FitsOutFile fios(arg[2]);
91
92 cout << " skymixer Version= " << skm_version << " SophyaVersion= "
93 << SophyaVersion() << endl;
94
95 string msg;
96 int rc = 0;
97 RadSpectra * es = NULL;
98 SpectralResponse * sr = NULL;
99 double moy, sig;
100
101 DataCards dc;
102 so = NULL;
103 // DVList for merging all FITS keywords and datacards
104 dvl_fitskw = new DVList;
105
106 try {
107 string dcard = arg[1];
108 if(printlev > 1) cout << " Decoding parameters from file " << dcard << endl;
109 dc.ReadFile(dcard);
110
111 rc = CheckCards(dc, msg,fios);
112 if (rc) {
113 cerr << " Error condition -> Rc= " << rc << endl;
114 cerr << " Msg= " << msg << endl;
115 return(rc);
116 }
117 }
118 catch (PException exc) {
119 msg = exc.Msg();
120 cerr << " !!!! skymixer/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
121 return(90);
122 }
123
124
125 cout << " skymix/Info : NComp = " << nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
126 cout << " ... MapPath = " << (string)mapPath << " DebugLev= " << debuglev
127 << " PrintLev= " << printlev << endl;
128
129 // We create an output persist file for writing debug objects
130 if (debuglev > 0) so = new POutPersist("skymixdbg.ppf");
131 int countHead = 2;
132 SphereHEALPix<float> outgs(hp_nside);
133 try{
134 if (rdmap) { // Reading map from FITS file
135 char ifnm[256];
136 strncpy(ifnm, dc.SParam("READMAP", 0).c_str(), 255);
137 ifnm[255] = '\0';
138 cout << " Reading output HealPix map from FITS file " << (string)ifnm << endl;
139 {
140 FitsInFile fitsin(ifnm);
141 fitsin >> outgs;
142 // Getting FITS keywords in primary header
143 // SKM_MergeFITSKeywords(ifnm);
144 // Getting FITS keywords in data object header
145 if(checkpdmtype)
146 {
147 int status = 1;
148 string whatsup=fitsin.getStringKeyword(1,"PDMTYPE",status);
149 if(!(status==0&&whatsup=="COMPMAP"))
150 {
151 cerr << " !!!! skymixer - file" << ifnm << " is not a COMPMAP " << endl;
152 throw ParmError("Error");
153 }
154 else cout << "read the comment COMPMAP from " << ifnm << endl;
155 }
156 if(dvliston) SKM_MergeFITSKeywords2(outgs.Info());
157 else
158 {
159 fios.appendInputHeader(fitsin, countHead);
160 }
161 }
162 if(printlev>0)
163 cout << " Output HealPIx Map read - NbPixels= " <<
164 outgs.NbPixels() << endl;
165 if (printlev > 0) {
166 MeanSig(outgs.DataBlock(), moy, sig );
167 cout << " MeanSig for outpout map - Mean= " <<
168 moy << " Sigma= " << sig << endl;
169 }
170 }
171 else {
172 if(printlev>0)
173 cout << " Output HealPix Map created - NbPixels= " <<
174 outgs.NbPixels() << endl;
175 outgs.SetPixels(0.);
176 }
177
178 // Decoding detection pass-band filter
179 sr = getSpectralResponse(dc,fios);
180 PrtTim(" After FilterCreation ");
181
182 char * flnm, buff[90];
183 string key;
184
185 double K = 1.;
186 double freqOfMap = 1.;
187 // Loop over sky component
188 int sk;
189 for(sk = 0; sk<nskycomp; sk++) {
190 cout << "------------------------------------" << endl;
191 cout << " Processing sky component No " << sk+1 << endl;
192
193
194 sprintf(buff, "%d", sk+1);
195 key = (string)"DIPOLE" + buff;
196 // check for an eventual dipole
197 if(dc.HasKey(key))
198 {
199 if (es) { delete es; es = NULL; }
200 double temp = -10.;
201 double theta= dc.DParam(key,1,1.);
202 double phi = dc.DParam(key,2,1.);
203 double amp = dc.DParam(key,3,1.);
204 if(dc.NbParam(key)>3)
205 {
206 temp = dc.DParam(key,4,1.);
207 }
208 cout << " creating dipole " << temp << " " << amp << " " << phi << " " << theta << " " << endl;
209 addDipole(*sr, outgs,theta,phi,amp,temp);
210 }
211 else
212 {
213 sprintf(buff, "%d", sk+1);
214 key = (string)"MAPFITSFILE" + buff;
215 flnm = BuildFITSFileName(dc.SParam(key, 0));
216
217 K = dc.DParam(key, 1, 1.);
218
219 cout << " Reading Input FITS map " << (string)flnm << endl;
220 SphereHEALPix<float> ings(hp_nside);
221 {
222 FitsInFile fiosIn(flnm);
223 fiosIn >> ings;
224 // Getting FITS keywords in primary header
225 // SKM_MergeFITSKeywords(flnm);
226 // Getting FITS keywords in data object header
227 //if(fiosIn.hasKeyword(string("PDMTYPE"), 1)) cout << "read the comment COMPMAP from " << flnm << endl;
228 if(checkpdmtype)
229 {
230 int status = 1;
231 string ch = "PDMTYPE";
232 string whatsup=fiosIn.getStringKeyword(1,ch,status);
233 if(!(status==0&&whatsup=="COMPMAP"))
234 {
235 cerr << " !!!! skymixer - file" << flnm << " is not a COMPMAP " << endl;
236 throw ParmError("Error");
237 }
238 else cout << "read the comment COMPMAP from " << flnm << endl;
239 }
240 if(dvliston) SKM_MergeFITSKeywords2(ings.Info());
241 else
242 {
243 fios.appendInputHeader(fiosIn, countHead);
244 }
245 }
246 if (debuglev > 4) { // Writing the input map to the outppf
247 FIO_SphereHEALPix<float> fiog(ings);
248 fiog.Write(*so, key);
249 }
250 if (printlev > 2) {
251 MeanSig(ings.DataBlock(), moy, sig );
252 cout << " MeanSig for input map - Mean= " << moy << " Sigma= " << sig << endl;
253 }
254 bool mapDependentOfFreq = false;
255 key = (string)"BETAFITSFILE"+ buff;
256 if(dc.HasKey(key))
257 {
258 mapDependentOfFreq = true;
259 }
260
261 // getting Emission spectra
262 if(!mapDependentOfFreq)
263 {
264 if (es) { delete es; es = NULL; }
265 es = getEmissionSpectra(dc, sk,fios);
266 addComponent(*sr, outgs, ings, *es, K);
267 }
268 else
269 {
270 key = (string)"BETAFITSFILE"+ buff;
271 //SphereHEALPix<float> betaMap;
272 flnm = BuildFITSFileName(dc.SParam(key, 0));
273 double normFreq = dc.DParam(key, 1, 1.);
274 if (printlev > 4) cout << "....BetaFits... normalization Freq = " << normFreq << endl;
275 int nSideForInt = dc.DParam(key, 2, 1.);
276 if (printlev > 4) cout << "....BetaFits... NSide for Integration map = " << nSideForInt << endl;
277 cout << "....BetaFits... Reading Beta FITS map " << (string)flnm << endl;
278 SphereHEALPix<float> betaMap(hp_nside);
279 {
280 FitsInFile fiosBM(flnm);
281 fiosBM >> betaMap;
282 // Getting FITS keywords in primary header
283 // SKM_MergeFITSKeywords(flnm);
284 // Getting FITS keywords in data object header
285 //if(fiosBM.hasKeyword(string("COMPMAP"), 1)) cout << "read the comment COMPMAP from " << flnm << endl;
286 if(checkpdmtype)
287 {
288 int status = 1;
289 string whatsup=fiosBM.getStringKeyword(1,"PDMTYPE",status);
290 if(!(status==0&&whatsup=="INDEXMAP"))
291 {
292 cerr << " !!!! skymixer - file" << flnm << " is not an INDEXMAP " << endl;
293 throw ParmError("Error");
294 }
295 else cout << "read the comment INDEXMAP from " << flnm << endl;
296 }
297 if(dvliston) SKM_MergeFITSKeywords2(betaMap.Info());
298 else
299 {
300 fios.appendInputHeader(fiosBM, countHead);
301 }
302 }
303 if (printlev > 2) {
304 MeanSig(betaMap.DataBlock(), moy, sig );
305 cout << " MeanSig for Beta map - Mean= " << moy << " Sigma= " << sig << endl;
306 }
307 if (debuglev > 4) { // Writing the input map to the outppf
308 FIO_SphereHEALPix<float> fiogs(betaMap);
309 fiogs.Write(*so, key);
310 }
311
312 if(nSideForInt<0) nSideForInt = sqrt((double)betaMap.NbPixels()/12);
313 bool bydefault = true;
314 if(!bydefault)
315 addComponentBeta(outgs,ings,*sr,betaMap,normFreq, K);
316 else
317 {
318 // integrate the betamap over the SpectralResponse
319 SphereHEALPix<float> intBetaMap(nSideForInt);
320 integratedMap(*sr, betaMap, normFreq, intBetaMap);
321 if (debuglev > 4) { // Writing the input map to the outppf
322 FIO_SphereHEALPix<float> fiogs2(intBetaMap);
323 fiogs2.Write(*so, "INTBETAMAP");
324 }
325
326 betaMap.Resize(8);
327 if (printlev > 4)
328 {
329 MeanSig(intBetaMap.DataBlock(), moy, sig );
330 cout << "....BetaFits... MeanSig for intBetaMap - Mean= "
331 << moy << " Sigma= " << sig << endl;
332 }
333 // add the integrated beta map
334 addComponentBeta(outgs,ings,intBetaMap, K);
335 }
336 }
337
338 MeanSig(outgs.DataBlock(), moy, sig );
339 cout << " MeanSig for Sum map - Mean= " << moy << " Sigma= " << sig << endl;
340 cout << "-------------------------------------------------" << endl;
341
342 sprintf(buff, "End of Processing Component %d ", sk+1);
343 PrtTim(buff);
344 }
345 }
346 }
347 catch(PException exc)
348 {
349 cout << "catched PException" << endl;
350 msg = exc.Msg();
351 cerr << " !!!! skymixer - Catched exception - Msg= " << exc.Msg() << endl;
352 rc = 50;
353 return(50);
354 }
355
356 try {
357 // Saving the output map in FITS format
358 cout << "Output Map (SphereHEALPix<float>) written to FITS file "
359 << (string)(arg[2]) << endl;
360 {
361 // FitsOutFile fios(arg[2]);
362 fios.firstImageOnPrimaryHeader(false); // Use secondary header
363 DVList& dvl = (*dvl_fitskw);
364 dvl["PDMTYPE"] = "COMPMAP";
365 dvl.SetComment("PDMTYPE", "Planck Data Model Type");
366 dvl["SOPHYVER"] = SophyaVersion();
367 dvl.SetComment("SOPHYVER", "Sophya Version number");
368 dvl["SKYMVER"] = skm_version;
369 dvl.SetComment("SKYMVER", "skymixer Version number");
370 // fios.DVListIntoPrimaryHeader(dvl);
371 dvl["CREATOR"] = "SkyMixer";
372 outgs.Info() = dvl;
373 fios << outgs ;
374 }
375 PrtTim("End of WriteFITS ");
376 // Saving the output map in PPF format
377 if (narg > 3) {
378 POutPersist s(arg[3]);
379 FIO_SphereHEALPix<float> fiog(&outgs) ;
380 fiog.Write(s);
381 cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "
382 << (string)(arg[3]) << endl;
383 PrtTim("End of WritePPF ");
384 }
385 }
386 catch(PException exc)
387 {
388 cout << "catched PException (2)" << endl;
389 msg = exc.Msg();
390 cerr << " !!!! skymixer(2) - Catched exception - Msg= " << exc.Msg() << endl;
391 rc = 55;
392 return(55);
393 }
394
395 if (so) delete so; // Closing the debug ppf file
396 return(rc);
397}
398
399/* Nouvelle-Fonction */
400int CheckCards(DataCards & dc, string & msg, FitsOutFile& fios)
401// Function to check datacards
402{
403rdmap = false;
404mapPath[0] = '\0';
405hp_nside = 32;
406nskycomp = 0;
407debuglev = 0;
408printlev = 0;
409
410int rc = 0;
411string key, key2,key3;
412
413 DVList & dvl = *dvl_fitskw;
414 // Cheking datacards
415 if (dc.NbParam("SKYMIX") < 2) {
416 rc = 71;
417 msg = "Invalid parameters - Check @SKYMIX card ";
418 return(rc);
419 }
420 if(dvliston)
421 dvl.SetS("SKYMIX", dc.GetParams("SKYMIX"));
422 else
423 {
424 fios.insertCommentLineOnHeader("--------------------");
425 fios.insertCommentLineOnHeader("SkyMixer DataCard");
426 fios.insertCommentLineOnHeader("--------------------");
427
428 string key = "SKYMIX"+ (string)dc.GetParams("SKYMIX");
429 fios.insertCommentLineOnHeader(key);
430 }
431
432 key = "CHECKPDMTYPE";
433 if (dc.HasKey(key)) {
434 checkpdmtype=true;
435 }
436
437
438 key = "READMAP";
439 if (dc.HasKey(key)) {
440 if (dc.NbParam(key) < 1) {
441 rc = 72;
442 msg = "Invalid parameters - Check @READMAP card ";
443 return(rc);
444 }
445 else rdmap = true;
446
447 if(dvliston) dvl.SetS(key, dc.GetParams(key));
448 else
449 {
450 string keynew = key + (string)dc.GetParams(key);
451 fios.insertCommentLineOnHeader(keynew);
452 }
453 }
454
455// Checking detection filter specification
456 key = "GAUSSFILTER";
457 key2 = "FILTERFITSFILE";
458 key3 = "DIPOLE";
459 if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3) && (dc.NbParam(key3) < 3)) {
460 msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE or @DIPOLE";
461 rc = 73; return(rc);
462 }
463
464 if (dc.HasKey(key))
465 {
466 if(dvliston) dvl.SetS("GAUSFILT", dc.GetParams(key));
467 else
468 {
469 string keynew = key + (string)dc.GetParams(key);
470 fios.insertCommentLineOnHeader(keynew);
471 }
472 }
473 if (dc.HasKey(key2))
474 {
475 if(dvliston) dvl.SetS("FILTFILE", dc.GetParams(key2));
476 else
477 {
478 string keynew = key2 + (string)dc.GetParams(key2);
479 fios.insertCommentLineOnHeader(keynew);
480 }
481 }
482 if (dc.HasKey(key3))
483 {
484 if(dvliston) dvl.SetS("DIPOLE", dc.GetParams(key3));
485 else
486 {
487 string keynew = key3 + (string)dc.GetParams(key3);
488 fios.insertCommentLineOnHeader(keynew);
489 }
490 }
491 // Decoding number of component and pixelisation parameter
492 int mg = 32;
493 int ncomp = 0;
494 ncomp = dc.IParam("SKYMIX", 0, 0);
495 mg = dc.IParam("SKYMIX", 1, 32);
496 if (ncomp < 1) {
497 msg = "Invalid parameters - Check datacards @SKYMIX ";
498 rc = 74;
499 return(rc);
500 }
501
502 // Checking detection filter specification
503 // Checking input FITS file specifications
504 int kc;
505 char buff[256];
506 bool pb = false;
507 string key4;
508 string key5;
509 string key6;
510 for(kc=0; kc<ncomp; kc++) {
511 sprintf(buff, "MAPFITSFILE%d", kc+1);
512 key = buff;
513 sprintf(buff, "DIPOLE%d", kc+1);
514 key3 = buff;
515 if (dc.NbParam(key) < 1 && (dc.NbParam(key3)<1))
516 {
517 msg = "Invalid card or Missing card" + key + " or " + key3 ;
518 pb = true; break;
519 }
520
521 sprintf(buff, "MAPFIL%d", kc+1);
522 if (dc.HasKey(key))
523 {
524 if(dvliston) dvl.SetS(buff, dc.GetParams(key));
525 else
526 {
527 string keynew = buff + (string)dc.GetParams(key);
528 fios.insertCommentLineOnHeader(keynew);
529 }
530 }
531 sprintf(buff, "DIPOLE%d", kc+1);
532 if (dc.HasKey(key3))
533 {
534 if(dvliston) dvl.SetS(buff, dc.GetParams(key3));
535 else
536 {
537 string keynew = buff + (string)dc.GetParams(key3);
538 fios.insertCommentLineOnHeader(keynew);
539 }
540 }
541 sprintf(buff, "SPECTRAFITSFILE%d", kc+1);
542 key = buff;
543 sprintf(buff, "BLACKBODY%d", kc+1);
544 key2 = buff;
545 sprintf(buff, "POWERLAWSPECTRA%d", kc+1);
546 key3 = buff;
547 sprintf(buff, "BETAFITSFILE%d", kc+1);
548 key4 = buff;
549 sprintf(buff, "DIPOLE%d", kc+1);
550 key5 = buff;
551 sprintf(buff, "DERIVBB%d", kc+1);
552 key6 = buff;
553 if ( (dc.NbParam(key) < 3) && (dc.NbParam(key2) < 1) && (dc.NbParam(key3) < 6) && (dc.NbParam(key4)<2)
554 && (dc.NbParam(key6)<1) && (dc.NbParam(key5)<1)) {
555 msg = "Missing card or invalid parameters : " + key + " " + key2 + " " + key3 + " " + key4+ " " + key5;
556 pb = true; break;
557 }
558
559 sprintf(buff, "SPECTF%d", kc+1);
560 if (dc.HasKey(key)) {
561 if(dvliston) dvl.SetS(buff, dc.GetParams(key));
562 else {
563 string keynew = buff + (string)dc.GetParams(key);
564 fios.insertCommentLineOnHeader(keynew);}
565 }
566 sprintf(buff, "BLACKB%d", kc+1);
567 if (dc.HasKey(key2)) {
568 if(dvliston) dvl.SetS(buff, dc.GetParams(key2));
569 else {
570 string keynew = buff + (string)dc.GetParams(key2);
571 fios.insertCommentLineOnHeader(keynew);}
572 }
573
574 sprintf(buff, "PLAWSP%d", kc+1);
575 if (dc.HasKey(key3)) {
576 if(dvliston) dvl.SetS(buff, dc.GetParams(key3));
577 else {
578 string keynew = buff + (string)dc.GetParams(key3);
579 fios.insertCommentLineOnHeader(keynew);}
580 }
581
582 sprintf(buff, "BETAFI%d", kc+1);
583 if (dc.HasKey(key4)) {
584 if(dvliston) dvl.SetS(buff, dc.GetParams(key4));
585 else {
586 string keynew = buff + (string)dc.GetParams(key4);
587 fios.insertCommentLineOnHeader(keynew);}
588 }
589
590 sprintf(buff, "DIPOLE%d", kc+1);
591 if (dc.HasKey(key5)) {
592 if(dvliston) dvl.SetS(buff, dc.GetParams(key5));
593 else {
594 string keynew = buff + (string)dc.GetParams(key5);
595 fios.insertCommentLineOnHeader(keynew);}
596 }
597
598 sprintf(buff, "DERIBB%d", kc+1);
599 if (dc.HasKey(key6)) {
600 if(dvliston) dvl.SetS(buff, dc.GetParams(key6));
601 else {
602 string keynew = buff + (string)dc.GetParams(key6);
603 fios.insertCommentLineOnHeader(keynew);}
604 }
605
606 }
607
608 if (pb) {
609 rc = 75;
610 return(75);
611 }
612
613
614// Initialiazing parameters
615 rc = 0;
616 msg = "OK";
617 nskycomp = ncomp;
618 hp_nside = mg;
619
620
621// Checking for PATH definition card
622 key = "MAPPATH";
623 if (dc.NbParam(key) < 3) strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
624 mapPath[255] = '\0';
625 key = "DEBUGLEVEL";
626 debuglev = dc.IParam(key, 0, 0);
627 key = "PRINTLEVEL";
628 printlev = dc.IParam(key, 0, 0);
629 return(rc);
630}
631
632static char buff_flnm[1024]; // Mal protege !
633/* Nouvelle-Fonction */
634char* BuildFITSFileName(string const & fname)
635{
636if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
637else sprintf(buff_flnm, "%s", fname.c_str());
638return(buff_flnm);
639}
640
641/* Nouvelle-Fonction */
642SpectralResponse * getSpectralResponse(DataCards & dc, FitsOutFile& fios)
643{
644 SpectralResponse * filt = NULL;
645
646 string key = "FILTERFITSFILE";
647 string key2 = "GAUSSFILTER";
648 string ppfname = "filter";
649
650 if (dc.HasKey(key) ) { // Reading FITS filter file
651 char ifnm[256];
652 strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
653 ifnm[255] = '\0';
654 Matrix mtx;
655 FitsInFile fiis(ifnm);
656 fiis.firstImageOnPrimaryHeader(false); // Use secondary header HDU=2
657 fiis >> mtx ;
658 // Getting FITS keywords in primary header
659 // SKM_MergeFITSKeywords(ifnm);
660 // Getting FITS keywords in data object header
661 if(checkpdmtype)
662 {
663 int status = 1;
664 string whatsup=fiis.getStringKeyword(1,"PDMTYPE",status);
665 if(!(status==0&&whatsup=="COMPMAP"))
666 {
667 cerr << " !!!! skymixer - file" << ifnm << " is not a COMPMAP " << endl;
668 throw ParmError("Error");
669 }
670 else cout << "read the comment COMPMAP from " << ifnm << endl;
671 }
672 if(dvliston) SKM_MergeFITSKeywords2(mtx.Info());
673 else
674 {
675 fios.appendInputHeader(fiis, 2);
676 }
677 double numin = dc.DParam(key, 2, 1.);
678 double numax = dc.DParam(key, 3, 9999.);
679 Vector nu(mtx.NCols());
680 Vector fnu(mtx.NCols());
681 for(int k=0; k<mtx.NCols(); k++) {
682 nu(k) = mtx(0, k);
683 fnu(k) = mtx(1, k);
684 }
685 filt = new SpecRespVec(nu, fnu, numin, numax);
686 ppfname = key;
687 }
688 else if (dc.HasKey(key2) ) { // creating GaussianFilter
689 double nu0 = dc.DParam(key2, 0, 10.);
690 double s = dc.DParam(key2, 1, 1.);
691 double a = dc.DParam(key2, 2, 1.);
692 double numin = dc.DParam(key2, 3, 0.1);
693 double numax = dc.DParam(key2, 4, 9999);
694 filt = new GaussianFilter(nu0, s, a, numin, numax);
695 ppfname = key2;
696 }
697 if (filt == NULL) throw ParmError("datacard error ! No detection filter");
698 if(printlev>0)
699 {
700 cout << endl;
701 cout << " Filter decoded - Created " << endl;
702 cout << *filt << endl;
703 }
704 // for debug
705 if (debuglev > 1) SpectralResponse2Nt(*filt, *so, ppfname);
706 return(filt);
707}
708
709/* Nouvelle-Fonction */
710RadSpectra * getEmissionSpectra(DataCards & dc, int nc, FitsOutFile& fios)
711{
712 char numb[16];
713 sprintf(numb, "%d", nc+1);
714
715 string key = (string)"SPECTRAFITSFILE" + numb;
716 string key2 = (string)"BLACKBODY" + numb;
717 string key5 = (string)"DERIVBB" + numb;
718 string key3 = (string)"POWERLAWSPECTRA" + numb;
719 string ppfname = "espectra";
720
721 RadSpectra * rs = NULL;
722 if (dc.HasKey(key) ) { // Reading emission spectra from file
723 char * ifnm = BuildFITSFileName(dc.SParam(key, 0));
724 cout << " Reading Input FITS spectra file " << (string)ifnm << endl;
725 Matrix mtx;
726 FitsInFile fiis(ifnm);
727 fiis.firstImageOnPrimaryHeader(false); // Use secondary header HDU=2
728 fiis >> mtx ;
729 // Getting FITS keywords in primary header
730 // SKM_MergeFITSKeywords(ifnm);
731 // Getting FITS keywords in data object header
732 //if(fiis.hasKeyword(string("PDMTYPE"), 1)) cout << "read the comment COMPMAP from " << ifnm << endl;
733 if(checkpdmtype)
734 {
735 int status = 1;
736 string whatsup=fiis.getStringKeyword(1,"PDMTYPE",status);
737 if(!(status==0&&whatsup=="FGRSPEC"))
738 {
739 cerr << " !!!! skymixer - file" << ifnm << " is not a FGRSPEC" << endl;
740 throw ParmError("Error");
741 }
742 else cout << "read the comment FGRSPEC from " << ifnm << endl;
743 }
744 if(dvliston) SKM_MergeFITSKeywords2(mtx.Info());
745 else
746 {
747 fios.appendInputHeader(fiis, 2);
748 }
749 double numin = dc.DParam(key, 2, 1.);
750 double numax = dc.DParam(key, 3, 9999.);
751 Vector nu(mtx.NCols());
752 Vector tnu(mtx.NCols());
753 for(int k=0; k<mtx.NCols(); k++) {
754 nu(k) = mtx(0, k);
755 tnu(k) = mtx(1, k);
756 }
757 rs = new RadSpectraVec(nu, tnu, numin, numax);
758 ppfname = key;
759 }
760 else if (dc.HasKey(key2) ) { // Creating BlackBody emission spectra
761 rs = new BlackBody(dc.DParam(key2, 0, 2.726));
762 ppfname = key2;
763 }
764 else if (dc.HasKey(key5) ) { // Creating Dipole
765 rs = new DerivBlackBody(dc.DParam(key5, 0, 3.E-3));
766 ppfname = key5;
767 }
768 else if (dc.HasKey(key3) ) { // Creating PowerLaw emission spectra
769 double a = dc.DParam(key3, 0, 1.);
770 double nu0 = dc.DParam(key3, 1, 100.);
771 double dnu = dc.DParam(key3, 2, 10.);
772 double b = dc.DParam(key3, 3, 0.);
773 double numin = dc.DParam(key3, 4, 0.1);
774 double numax = dc.DParam(key3, 5, 9999);
775 rs = new PowerLawSpectra(a, b, nu0, dnu, numin, numax);
776 ppfname = key3;
777 }
778 if (rs == NULL) throw ParmError("datacard error ! missing Emission spectra");
779 cout << " Emission spectra decoded - Created (" << ppfname << ")" << endl;
780 cout << *rs << endl;
781// for debug
782 if (debuglev > 2) RadSpec2Nt(*rs, *so, ppfname);
783 return(rs);
784}
785
786
787
788
789template <class T>
790void addDipole(SpectralResponse& sr, PixelMap<T>& finalMap,
791 double theta,double phi,double amp,double temp)
792{
793 DerivBlackBody dbb;
794 if(temp>0) dbb.setTemperature(temp);
795 double coeff = dbb.filteredIntegratedFlux(sr) * amp;
796 UnitVector vd(theta,phi);
797 UnitVector vc(theta,phi);
798
799 for(int i=0; i<finalMap.NbPixels(); i++)
800 {
801 double thetar,phir;
802 finalMap.PixThetaPhi(i,thetar,phir);
803 vc.SetThetaPhi(thetar, phir);
804 finalMap(i) += vd.Psc(vc)*coeff;
805 }
806 if (debuglev > 4) { // Writing the input map to the outppf
807 SphereHEALPix<float> ings(sqrt((double)finalMap.NbPixels()/12));
808 for(int i=0; i<finalMap.NbPixels(); i++)
809 {
810 double thetar,phir;
811 finalMap.PixThetaPhi(i,thetar,phir);
812 vc.SetThetaPhi(thetar, phir);
813 ings(i) = vd.Psc(vc);
814 }
815 FIO_SphereHEALPix<float> fiog(ings);
816 fiog.Write(*so, "dipole");
817 cout << "Debug the dipole map....saved in debug file !" << endl;
818 }
819}
820/* Nouvelle-Fonction */
821template <class T>
822void addComponent(SpectralResponse& sr, PixelMap<T>& finalMap,
823 PixelMap<T>& mapToAdd, RadSpectra& rs, double K)
824{
825 // finalMap = finalMap + coeff* mapToAdd
826 // coeff = convolution of sr and rs
827 // compute the coefficient corresponding to mapToAdd
828 if (finalMap.NbPixels() != mapToAdd.NbPixels())
829 throw SzMismatchError("addComponent()/Error: Unequal number of Input/Output map pixels");
830 double coeff = rs.filteredIntegratedFlux(sr) * K;
831 if (printlev > 1)
832 cout << " addComponent - Coeff= " << coeff << " (K= " << K << ")" << endl;
833 for(int i=0; i<finalMap.NbPixels(); i++)
834 {
835 finalMap(i) += coeff * mapToAdd(i);
836 }
837}
838/* Nouvelle-Fonction */
839template <class T>
840void addComponentBeta(SphereHEALPix<T>& finalMap,
841 SphereHEALPix<T>& mapToAdd,SpectralResponse& sr,
842 SphereHEALPix<T>& betaMap, double normFreq, double K)
843{
844 // finalMap = finalMap + coeff* mapToAdd
845 // coeff = convolution of sr and rs
846 // compute the coefficient corresponding to mapToAdd
847 // betaMap is the map of (beta(theta,phi))
848
849 int nbpix = finalMap.NbPixels();
850 if (nbpix != mapToAdd.NbPixels())
851 throw SzMismatchError("addComponentBeta()/Error: Unequal number of Input/Output map pixels");
852 if (printlev > 1)
853 {
854 cout << "addComponentBeta - Coeff= " << K << endl;
855 cout << "nb pixels: " << finalMap.NbPixels() << endl;
856 }
857 SphereHEALPix<T> bigBetaMap(sqrt((double)nbpix/12));
858 if(nbpix != betaMap.NbPixels())
859 {
860 Sph2Sph(betaMap,bigBetaMap);
861 }
862 for(int i=0; i<finalMap.NbPixels(); i++)
863 {
864 // coeff = integration of (nu/normFreq)^(-beta(theta,phi)) in each pixels
865 RadSpectra* rs = new PowerLawSpectra
866 (1.,-bigBetaMap(i), 0., normFreq, 0.01, 800.);
867 double coeff = rs->filteredIntegratedFlux(sr);
868 finalMap(i) += coeff*K*mapToAdd(i);
869 }
870}
871
872template <class T>
873void integratedMap(SpectralResponse& sr,
874 SphereHEALPix<T>& betaMap,
875 double normFreq,
876 SphereHEALPix<T>& intBetaMap)
877{
878 PowerLawSpectra rs(1.,-2., 0., normFreq);
879
880 if(betaMap.NbPixels()!=intBetaMap.NbPixels())
881 {
882 Sph2Sph(betaMap,intBetaMap);
883 for(int i=0; i<intBetaMap.NbPixels(); i++)
884 {
885 rs.setExp(-intBetaMap(i));
886 double coeff = rs.filteredIntegratedFlux(sr);
887 intBetaMap(i) = coeff;
888 }
889 }
890 else
891 {
892 for(int i=0; i<intBetaMap.NbPixels(); i++)
893 {
894 rs.setExp(-betaMap(i));
895 double coeff = rs.filteredIntegratedFlux(sr);
896 intBetaMap(i) = coeff;
897 }
898 }
899}
900
901template <class T>
902void addComponentBeta(SphereHEALPix<T>& finalMap,
903 SphereHEALPix<T>& mapToAdd,SphereHEALPix<T>& intBetaMap, double K)
904{
905 // finalMap = finalMap + coeff* mapToAdd
906 // coeff = convolution of sr and rs
907 // compute the coefficient corresponding to mapToAdd
908 // integBetaMap is the map of the integration (nu/normFreq)^(-beta(theta,phi)) over
909 // the spectralResponse
910 // different from addComponentBeta(PixelMap<T>& finalMap,
911 // PixelMap<T>& mapToAdd,SpectralResponse& sr, PixelMap<T>& betaMap, double normFreq, double K)
912 // since it permits to use a intBetaMap with a different number of pixels than
913 // the other maps
914
915 int nbpix = finalMap.NbPixels();
916 if (nbpix != mapToAdd.NbPixels())
917 throw SzMismatchError("addComponentBeta(PixelMap<T>&,PixelMap<T>&,PixelMap<T>&,double)/Error: Unequal number of Input/Output map pixels");
918 double coeff = K;
919
920 if(nbpix != intBetaMap.NbPixels())
921 {
922 for(int i=0; i<finalMap.NbPixels();i++)
923 {
924 double teta,phi,val;
925 finalMap.PixThetaPhi(i, teta, phi);
926 int pixel = intBetaMap.PixIndexSph(teta,phi);
927 val = intBetaMap.PixVal(pixel);
928 finalMap(i) += coeff*mapToAdd(i)*val;
929 }
930 }
931 else
932 {
933 for(int i=0; i<finalMap.NbPixels();i++)
934 {
935 finalMap(i) += coeff*mapToAdd(i)*intBetaMap(i);
936 }
937 }
938 if (printlev > 1)
939 {
940 cout << "addComponentBeta(SG<T>,SG<T>,SG<T>,double) - Coeff= " << K << endl;
941 }
942}
943
944
945
946
947
948/* Nouvelle-Fonction */
949void SKM_MergeFITSKeywords(char * flnm)
950{
951 DVList dvl;
952 FitsFile::FitsExtensionType typeOfExtension;
953 int naxis;
954 vector<int> naxisn;
955 FitsFile::FitsDataType dataType;
956 FitsInFile::GetBlockType(flnm, 1, typeOfExtension, naxis, naxisn, dataType, dvl);
957 // Cleaning the keywords
958#define SZexlst 21
959 char *exlst[SZexlst]=
960 {"SIMPLE","BITPIX" ,"NAXIS" ,"NAXIS#" ,"PCOUNT","GCOUNT",
961 "EXTEND","ORIGIN" ,"DATE*" ,"TFIELDS","TTYPE#","TFORM#",
962 "TUNIT#","EXTNAME","CTYPE#","CRVAL#" ,"CRPIX#","CDELT#",
963 "XTENSION","INSTRUME","TELESCOP"};
964 char kwex[32];
965 int i,l;
966 for (i=0; i<SZexlst; i++) {
967 strncpy(kwex, exlst[i], 32);
968 l = strlen(kwex)-1;
969 if ((kwex[l] != '*') && (kwex[l] != '#')) {
970 dvl.DeleteKey(kwex);
971 }
972 else {
973 bool fgd = (kwex[l] == '#') ? true : false;
974 list<string> lstsup;
975 kwex[l] = '\0';
976 DVList::ValList::const_iterator it;
977 for (it = dvl.Begin(); it != dvl.End(); it++) {
978 if ((*it).first.substr(0,l) != kwex) continue;
979 if (fgd && !isdigit((*it).first[l])) continue;
980 lstsup.push_back((*it).first);
981 }
982 list<string>::iterator it2;
983 for (it2 = lstsup.begin(); it2 != lstsup.end(); it2++)
984 dvl.DeleteKey(*it2);
985 }
986 }
987 dvl_fitskw->Merge(dvl);
988}
989
990/* Nouvelle-Fonction */
991void SKM_MergeFITSKeywords2(DVList & dvl)
992{
993 // Cleaning the keywords
994#define SZexlst 21
995 char *exlst[SZexlst]=
996 {"SIMPLE","BITPIX" ,"NAXIS" ,"NAXIS#" ,"PCOUNT","GCOUNT",
997 "EXTEND","ORIGIN" ,"DATE*" ,"TFIELDS","TTYPE#","TFORM#",
998 "TUNIT#","EXTNAME","CTYPE#","CRVAL#" ,"CRPIX#","CDELT#",
999 "XTENSION","INSTRUME","TELESCOP"};
1000 char kwex[32];
1001 int i,l;
1002 for (i=0; i<SZexlst; i++) {
1003 strncpy(kwex, exlst[i], 32);
1004 l = strlen(kwex)-1;
1005 if ((kwex[l] != '*') && (kwex[l] != '#')) {
1006 dvl.DeleteKey(kwex);
1007 }
1008 else {
1009 bool fgd = (kwex[l] == '#') ? true : false;
1010 list<string> lstsup;
1011 kwex[l] = '\0';
1012 DVList::ValList::const_iterator it;
1013 for (it = dvl.Begin(); it != dvl.End(); it++) {
1014 if ((*it).first.substr(0,l) != kwex) continue;
1015 if (fgd && !isdigit((*it).first[l])) continue;
1016 lstsup.push_back((*it).first);
1017 }
1018 list<string>::iterator it2;
1019 for (it2 = lstsup.begin(); it2 != lstsup.end(); it2++)
1020 dvl.DeleteKey(*it2);
1021 }
1022 }
1023 dvl_fitskw->Merge(dvl);
1024}
1025
1026/* Nouvelle-Fonction */
1027void RadSpec2Nt(RadSpectra & rs, POutPersist & so, string name)
1028{
1029 char *ntn[2] = {"nu","fnu"};
1030 NTuple nt(2,ntn); // Creation NTuple (AVEC new )
1031 float xnt[2];
1032 double nu;
1033 double numin = rs.minFreq();
1034 double numax = rs.maxFreq();
1035 int nmax = 500;
1036 double dnu = (numax-numin)/nmax;
1037 for(int k=0; k<nmax; k++) {
1038 nu = numin+k*dnu;
1039 xnt[0] = nu;
1040 xnt[1] = rs.flux(nu);
1041 nt.Fill(xnt);
1042 }
1043 ObjFileIO<NTuple> oiont(nt);
1044 oiont.Write(so, name);
1045 return;
1046}
1047
1048/* Nouvelle-Fonction */
1049void SpectralResponse2Nt(SpectralResponse& sr, POutPersist & so, string name)
1050{
1051 char *ntn[2] = {"nu","tnu"};
1052 NTuple nt(2,ntn); // Creation NTuple (AVEC new )
1053 float xnt[2];
1054 double nu;
1055 double numin = sr.minFreq();
1056 double numax = sr.maxFreq();
1057 int nmax = 500;
1058 double dnu = (numax-numin)/nmax;
1059 for(int k=0; k<nmax; k++) {
1060 nu = numin+k*dnu;
1061 xnt[0] = nu;
1062 xnt[1] = sr.transmission(nu);
1063 nt.Fill(xnt);
1064 }
1065 ObjFileIO<NTuple> oiont(nt);
1066 oiont.Write(so, name);
1067 return;
1068}
1069
Note: See TracBrowser for help on using the repository browser.