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

Last change on this file since 3809 was 2615, checked in by cmv, 21 years ago

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