source: Sophya/trunk/SophyaProg/PMixer/extractRS.cc@ 880

Last change on this file since 880 was 880, checked in by ansari, 25 years ago

Adaptation SphereGorski -> SphereHEALPix apres modifs Sophie - Reza 11/4/2000
File size: 13.1 KB
RevLine 
[761]1#include "machdefs.h"
2#include <stdlib.h>
3#include <stdio.h>
4#include <fstream.h>
5#include <iostream.h>
6#include <math.h>
7
8#include <string>
9#include <vector>
10
11#include "tmatrix.h"
12#include "tvector.h"
13#include "vector3d.h"
[880]14#include "srandgen.h"
[761]15#include "timing.h"
16#include "sambainit.h"
17#include "pexceptions.h"
18#include "datacards.h"
19#include "fitsioserver.h"
20
21#include "radspecvector.h"
22#include "blackbody.h"
23#include "nupower.h"
24
25#include "squarefilt.h"
26#include "trianglefilt.h"
27#include "specrespvector.h"
28#include "gaussfilt.h"
29
30// -----------------------------------------------------------------
31// ------------- Function declaration ------------------------------
32int CheckCards(DataCards & dc, string & msg);
33char * BuildFITSFileName(string const & fname);
34SpectralResponse * getSpectralResponse(DataCards & dc);
35template <class T>
36void MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig);
37float ComputeFrequency(DataCards &,string&);
38
39// -----------------------------------------------------------------
40
41// ----- Global (static) variables ------------
42static bool rdmap = false; // true -> Read map first
43static char mapPath[256]; // Path for input maps
[877]44static int hp_nside = 32; // HealPix NSide
[761]45static int nskycomp = 0; // Number of sky components
46static int debuglev = 0; // Debug Level
47static int printlev = 0; // Print Level
48static POutPersist * so = NULL; // Debug PPFOut file
49
50// -------------------------------------------------------------------------
51// main program
52// -------------------------------------------------------------------------
53int main(int narg, char * arg[])
54{
[877]55 if ((narg < 3) || ((narg > 1) && (strcmp(arg[1], "-h") == 0) )) {
56 cout << " Usage: extractRS parameterFile outputfitsnameformaps [outppfname]" << endl;
[761]57 exit(0);
58 }
59
60 InitTim();
61
62 string msg;
63 int rc = 0;
64 RadSpectra * es = NULL;
65 SpectralResponse * sr = NULL;
66 double moy, sig;
67
68 DataCards dc;
69 so = NULL;
70
71 try {
72 string dcard = arg[1];
73 cout << " Decoding parameters from file " << dcard << endl;
74 dc.ReadFile(dcard);
75
76 rc = CheckCards(dc, msg);
77 if (rc) {
78 cerr << " Error condition -> Rc= " << rc << endl;
79 cerr << " Msg= " << msg << endl;
80 return(rc);
81 }
82 }
83 catch (PException exc) {
84 msg = exc.Msg();
85 cerr << " !!!! extractRS/ Readcard - Catched exception - Msg= " << exc.Msg() << endl;
86 return(90);
87 }
88
89
[880]90 cout << " extractRS/Info : NComp = " << nskycomp << " SphereHEALPix_NSide= " << hp_nside << endl;
[761]91 cout << " ... MapPath = " << (string)mapPath << " DebugLev= " << debuglev
92 << " PrintLev= " << printlev << endl;
93
94
95 // We create an output persist file for writing debug objects
96 if (debuglev > 0) so = new POutPersist("extrRSdbg.ppf");
[880]97 SphereHEALPix<float> SourceMap(hp_nside);
98 SphereHEALPix<float> OutputMap(hp_nside);
[877]99 bool okout = true;
[761]100
101 try {
102
103 FitsIoServer fios; // Our FITS IO Server
104 char * flnm, buff[90];
105 string key;
106 string key2;
107
108 double K = 1.;
109
110 // Loop over sky frequencies components
111 int maxOfSkyComp = nskycomp;
112 int nbOfPix;
113 int sk;
114 Vector freq(maxOfSkyComp);
115 TMatrix<float> spectrum;
116 TMatrix<float> fdenu; //par freq et par pixel
117 int nb_source;
118 TMatrix<float> rareSource; //par freq et par numero de source
[877]119
[761]120 for(sk = 0; sk<maxOfSkyComp; sk++)
121 {
122 TMatrix<float> ings;
123 // Opening the file containing the map
124
125 if (printlev > 2) cout << " Processing map's frequency component No " << sk+1 << endl;
126 sprintf(buff, "%d", sk+1);
127 key = (string)"RADSPECMAP" + buff;
128
129 flnm = BuildFITSFileName(dc.SParam(key,0));
130 if (printlev > 2) cout << " Reading Input FITS map " << (string)flnm << endl;
131
132 fios.load(ings, flnm) ;
133 if (printlev > 2)
134 {
[877]135 double min = 1000000000;
136 double max = -10000000000;
137 for(int x=0;x<ings.NRows(); x++)
138 {
139 for(int y=0;y<ings.NCols(); y++)
140 {
141 if(min>ings(x,y)) min = ings(x,y);
142 if(max<ings(x,y)) max = ings(x,y);
143 }
144 }
145 cout << "min and max values for input map" << min << " " << max << endl;
[761]146 }
147 // Computing the frequency according to the name of the file
148 freq(sk) = ComputeFrequency(dc,key);
149 if (printlev > 2) cout << "filling spectrum for freq = " << freq(sk) << endl;
[877]150
[761]151 // Opening the file containing the point source
152 sprintf(buff, "%d", sk+1);
153 key2 = (string)"SOURCEPTMAP" + buff;
154 flnm = BuildFITSFileName(dc.SParam(key2,0));
155 const char* nameOfFile = flnm ;
156 if (printlev > 2) cout << " Reading Input FITS map " << nameOfFile << endl;
157 ifstream from(nameOfFile);
158 char ch;
159 string dum;
160 from >> nb_source >> dum;
161 if(sk==0) rareSource.ReSize(maxOfSkyComp,nb_source);
162 for (int ii=0; ii< nb_source; ii++)
163 {
164 double value = 0;
165 from >> value;
166 rareSource(sk,ii) = value*1.e-26; //for Jansky->W
167 }
168 from.close();
169
170 nbOfPix = ings.NRows()*ings.NCols();
171 spectrum.ReSize(ings.NRows(),ings.NCols());
172 if(sk==0) fdenu.ReSize(maxOfSkyComp,nbOfPix);
173 int xy=0;
174 for(int x=0;x<ings.NRows(); x++)
175 {
[877]176 for(int y=0;y<ings.NCols(); y++)
[761]177 {
178 spectrum(x,y) = ings(x,y)*1.e-26; //for Jansky->W
179 fdenu(sk,xy) = spectrum(x,y);
180 xy += 1;
181 }
182 }
183 if (printlev > 2)
184 {
185 MeanSig(spectrum.DataBlock(), moy, sig );
186 cout << " MeanSig for spectrum map - Mean= " << moy << " Sigma= " << sig << endl;
187 }
188
189
190 K = dc.DParam(key, 1, 1.);
191 if (debuglev > 4) { // Writing the input map to the outppf
[880]192 FIO_TArray<float> fiog;
[877]193 fiog.Write(*so, key);
[761]194 }
195
196 if(printlev>2)
197 {
198 sprintf(buff, "End of Proc. Comp. %d ", sk+1);
199 cout << " --------------------------------------- \n" << endl;
200 PrtTim(buff);
201 }
202 } // End of sky component loop
203
204 // ---------------------------------------------------------------------
205 // Integration in the detector band of the maps
206 // ---------------------------------------------------------------------
207 SpectralResponse* sr = NULL;
208 sr = getSpectralResponse(dc);
209
210 if (printlev > 2) cout << "SpectralResponse!" << *sr << endl;
211 TVector<float> outCoeff(nbOfPix);
212
[877]213 for(int pix=0;pix<nbOfPix;pix++)
[761]214 {
215 Vector smallFDeNu(maxOfSkyComp);
216 for(sk = 0; sk<maxOfSkyComp; sk++)
217 {
218 smallFDeNu(sk) = fdenu(sk,pix);
219 }
220
221 RadSpectraVec radSV(freq,smallFDeNu);
222 if (printlev > 10) cout << "RadiationSpectra" << radSV << endl;
223
224 outCoeff(pix) = radSV.filteredIntegratedFlux(*sr);
225 }
[877]226
227 for(int i=0;i<OutputMap.NbPixels();i++)
228 {
229 float random = frand01();
230 if(random == 1) random = frand01();
231 int randomPix = random*OutputMap.NbPixels();
232 int outpix = random*nbOfPix;
233 OutputMap(i) = outCoeff(outpix)*1.e18;
234 }
[761]235
[877]236
[761]237 // ---------------------------------------------------------------------
238 // Integration in the detector band of the sources
239 // ---------------------------------------------------------------------
240 TVector<float> sourceCoeff(nb_source);
241 for(int nsource=0;nsource<nb_source;nsource++)
242 {
243 Vector SourceFDeNu(maxOfSkyComp);
244 for(sk = 0; sk<maxOfSkyComp; sk++)
245 {
246 SourceFDeNu(sk) = rareSource(sk,nsource);
247 }
248 RadSpectraVec radSV_Source(freq,SourceFDeNu);
249 sourceCoeff(nsource) = radSV_Source.filteredIntegratedFlux(*sr);
250 }
251
[877]252 for(int i=0;i<nb_source;i++)
253 {
254 float random = frand01();
255 int outPix = random*OutputMap.NbPixels();
256 SourceMap(outPix) = sourceCoeff(i);
257 OutputMap(outPix) = OutputMap(i)+ SourceMap(outPix)*1.E18;
258 }
259 }
[761]260 catch (PException exc) {
261 msg = exc.Msg();
262 cerr << " !!!! extractRS - Catched PException - Msg= " << exc.Msg() << endl;
263 rc = 50;
264 }
265 catch (PThrowable exc) {
266 msg = exc.Msg();
267 cerr << " !!!! extractRS - Catched PThrowable - Msg= " << exc.Msg() << endl;
268 rc = 50;
269 }
270 catch (...) {
271 cerr << " Une exception de type inconnue !!! " << endl;
272 exit(99);
273 }
274
275 // Saving the output map in FITS format
[877]276 if (okout)
277 {
278 // {
279 // FitsIoServer fios;
280 // fios.save(SourceMap, arg[3]);
281 // double moy,sig;
282 // MeanSig(SourceMap.DataBlock(),moy,sig);
283 // cout << " MeanSig for Source Map - Mean= " << moy << " Sigma= " << sig << endl;
[880]284 // cout << "Source Map (SphereHEALPix<float>) written to FITS file "
[877]285 // << (string)(arg[3]) << endl;
286 // }
287 {
288 FitsIoServer fios;
289 fios.save(OutputMap, arg[2]);
290 double moy,sig;
291 MeanSig(OutputMap.DataBlock(),moy,sig);
292 cout << " MeanSig for Output Map - Mean= " << moy << " Sigma= " << sig << endl;
[880]293 cout << "Output Map (SphereHEALPix<float>) written to FITS file "
[877]294 << (string)(arg[2]) << endl;
295 }
296 if(printlev>2) PrtTim("End of WriteFITS ");
297 // Saving the output map in PPF format
298 if (narg > 2) {
299 POutPersist s(arg[3]);
[880]300 FIO_SphereHEALPix<float> fiog(OutputMap);
[877]301 fiog.Write(s);
[880]302 cout << "Output Map (SphereHEALPix<float>) written to POutPersist file "
[877]303 << (string)(arg[3]) << endl;
304 if(printlev>2) PrtTim("End of WritePPF ");
305 }
[761]306 }
307 if (so) delete so; // Closing the debug ppf file
308 return(rc);
309}
310
311
312
313/* Nouvelle-Fonction */
314int CheckCards(DataCards & dc, string & msg)
315 // Function to check datacards
316{
317 rdmap = false;
318 mapPath[0] = '\0';
319 nskycomp = 0;
320 debuglev = 0;
321 printlev = 0;
322
323 int rc = 0;
324 string key, key2;
325
326 // Checking datacards
[877]327 if (dc.NbParam("EXT_RS") < 2) {
[761]328 rc = 71;
329 msg = "Invalid parameters - Check @EXT_RS card ";
330 return(rc);
331 }
332 int kc;
333 char buff[256];
334 bool pb = false;
335 string key3;
336 int ncomp = 0;
337 ncomp = dc.IParam("EXT_RS", 0);
[877]338 int nside = 32;
339 nside = dc.IParam("EXT_RS", 1);
[761]340
341 for(kc=0; kc<ncomp; kc++)
342 {
343 key = "RS_EXT_PATH";
344 if (dc.NbParam(key) < 3) strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
345 mapPath[255] = '\0';
346 sprintf(buff, "RADSPECMAP%d", kc+1);
347 key = buff;
348
349 if (dc.NbParam(key) < 1) {
350 msg = "Missing or invalid card : " + key;
351 pb = true; break;
352 }
353 }
354
355 // Checking detection filter specification
356
357 key = "GAUSSFILTER";
358 key2 = "FILTERFITSFILE";
359 if ( (dc.NbParam(key) < 5) && (dc.NbParam(key2) < 3)) {
360 msg = "Missing card or parameters : Check @GAUSSFILTER or @FILTERFITSFILE";
361 rc = 73; return(rc);
362 }
363
364 if (pb) {
365 rc = 75;
366 return(75);
367 }
368
369
370 // Initialiazing parameters
371 rc = 0;
372 msg = "OK";
373 nskycomp = ncomp;
[877]374 hp_nside = nside;
[761]375 // Checking for PATH definition card
376 key = "RS_EXT_PATH";
377 if (dc.NbParam(key) < 3) strncpy(mapPath, dc.SParam(key, 0).c_str(), 255);
378 mapPath[255] = '\0';
379 key = "DEBUGLEVEL";
380 debuglev = dc.IParam(key, 0, 0);
381 key = "PRINTLEVEL";
382 printlev = dc.IParam(key, 0, 0);
383 return(rc);
384}
385
386static char buff_flnm[1024]; // Mal protege !
387/* Nouvelle-Fonction */
388char* BuildFITSFileName(string const & fname)
389{
390if (mapPath[0] != '\0') sprintf(buff_flnm, "%s/%s", mapPath, fname.c_str());
391else sprintf(buff_flnm, "%s", fname.c_str());
392return(buff_flnm);
393}
394
395
396/* Nouvelle-Fonction */
397template <class T>
398void MeanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig)
399
400{
401 gmoy=0.;
402 gsig = 0.;
403 double valok;
404 for(int k=0; k<(int)dbl.Size(); k++) {
405 valok = dbl(k);
406 gmoy += valok; gsig += valok*valok;
407 }
408 gmoy /= (double)dbl.Size();
409 gsig = gsig/(double)dbl.Size() - gmoy*gmoy;
410 if (gsig >= 0.) gsig = sqrt(gsig);
411}
412
413/* Nouvelle-Fonction */
414SpectralResponse * getSpectralResponse(DataCards & dc)
415{
416 SpectralResponse * filt = NULL;
417
418 string key = "FILTERFITSFILE";
419 string key2 = "GAUSSFILTER";
420 string ppfname = "filter";
421
422 if (dc.HasKey(key) ) { // Reading FITS filter file
423 FitsIoServer fios;
424 char ifnm[256];
425 strncpy(ifnm, dc.SParam(key, 1).c_str(), 255);
426 ifnm[255] = '\0';
427 Matrix mtx(2,10);
428 fios.load(mtx, ifnm);
429 double numin = dc.DParam(key, 2, 1.);
430 double numax = dc.DParam(key, 3, 9999.);
431 Vector nu(mtx.NCols());
432 Vector fnu(mtx.NCols());
433 for(int k=0; k<mtx.NCols(); k++) {
434 nu(k) = mtx(0, k);
435 fnu(k) = mtx(1, k);
436 }
437 filt = new SpecRespVec(nu, fnu, numin, numax);
438 ppfname = key;
439 }
440 else if (dc.HasKey(key2) ) { // creating GaussianFilter
441 double nu0 = dc.DParam(key2, 0, 10.);
442 double s = dc.DParam(key2, 1, 1.);
443 double a = dc.DParam(key2, 2, 1.);
444 double numin = dc.DParam(key2, 3, 0.1);
445 double numax = dc.DParam(key2, 4, 9999);
446 filt = new GaussianFilter(nu0, s, a, numin, numax);
447 ppfname = key2;
448 }
449 if (filt == NULL) throw ParmError("datacard error ! No detection filter");
450 cout << " Filter decoded - Created " << endl;
451 cout << *filt << endl;
452
453// for debug
454 // if (debuglev > 1) SpectralResponse2Nt(*filt, *so, ppfname);
455 return(filt);
456}
457
458
459float ComputeFrequency(DataCards & dc,string& key)
460{
461 ofstream filename("temp_filename");
462 filename << dc.SParam(key,0) << '\n';
463 filename.close();
464
465 ifstream fromtempfile("temp_filename");
466 char chtemp;
467 string temp_str;
468 int count = 0;
469 do {
470 fromtempfile.get(chtemp);
471 if(count> 12 && count < 20)
472 {
473 temp_str+= chtemp;
474 }
475 count+=1;
476 } while(chtemp!='\n');
477 fromtempfile.close();
478
479 ofstream filename2("temp_filename");
480 filename2 << temp_str << '\n';
481 filename2.close();
482
483 int frequency;
484 ifstream fromtempfile2("temp_filename");
485 fromtempfile2 >> frequency;
486 fromtempfile2.close();
487 return ((float)frequency/1000.);
488}
Note: See TracBrowser for help on using the repository browser.