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

Last change on this file since 928 was 928, checked in by ansari, 26 years ago

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