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

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

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