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

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

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