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

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

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