source: Sophya/trunk/AddOn/TAcq/visfits2dt.cc@ 3961

Last change on this file since 3961 was 3961, checked in by ansari, 15 years ago

suite debug visfits2dt.cc, Reza+Cedric, 04/03/2011
File size: 8.4 KB
Line 
1// Utilisation de SOPHYA pour faciliter les tests ...
2#include "sopnamsp.h"
3#include "machdefs.h"
4
5/* ----------------------------------------------------------
6 Projet BAORadio - (C) LAL/IRFU 2011
7
8 Programme de lecture des fichiers matrices de visibilites
9 au format fits pour ecrire une DataTable de visib_moyenne(temps)
10 Fichier de visibilites produits par le prog vismf
11 R. Ansari, C. Magneville - LAL/Irfu
12 V : Mai 2009
13 ---------------------------------------------------------- */
14
15// include standard c/c++
16#include <math.h>
17#include <stdio.h>
18#include <stdlib.h>
19#include <string.h>
20
21#include <iostream>
22#include <string>
23
24#include "pexceptions.h"
25#include "tvector.h"
26#include "fioarr.h"
27// #include "tarrinit.h"
28#include "ntuple.h"
29#include "datatable.h"
30#include "histinit.h"
31#include "matharr.h"
32#include "timestamp.h"
33#include <utilarr.h>
34#include "fitsarrhand.h"
35#include "fitshdtable.h"
36
37// include sophya mesure ressource CPU/memoire ...
38#include "resusage.h"
39#include "ctimer.h"
40#include "timing.h"
41
42
43//--------------------------- Fonctions de ce fichier -------------------
44int Usage(bool fgshort=true);
45
46int DecodeProc(int narg, char* arg[]);
47int ProcVisMtxFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq,
48 vector<sa_size_t> tfrlist);
49int FillVisDTable(BaseDataTable& visdt_, TMatrix< r_4 > vismtx_, TVector< int_4> chanum_,
50 sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) ;
51
52//------------------------------------------------------------------------------------------------------------
53
54/* --Fonction-- */
55int Usage(bool fgshort)
56{
57 cout << " --- visfits2dt : Read fits visibility matrices produced by visfmib to make \n"
58 << " time-frequency matrices and Visibilites=f(time) datatable " << endl;
59 cout << " Usage: visfits2dt InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [VisMtxRowList] \n" << endl;
60 if (fgshort) {
61 cout << " svv2mtx -h for detailed instructions" << endl;
62 return 1;
63 }
64 cout << " InOutPath : Input/Output directory name \n"
65 << " Imin,Imax,IStep: Input FITS files sequence number (vismtxII.fits) \n"
66 << " FileNames=InOutPath/vismtxII.fits Imin<=II<=Imax II+=IStep \n"
67 << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for DataTable\n"
68 << " VisMtxRowList : List of visibiliy matrix rows (0,2,...) -> time-freq map \n" << endl;
69 return 1;
70}
71
72//----------------------------------------------------
73//----------------------------------------------------
74int main(int narg, char* arg[])
75{
76 if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
77 if (narg<4) return Usage(true);
78 HiStatsInitiator _inia;
79 // TArrayInitiator _inia;
80
81 int rc = 0;
82 try {
83 bool fgvjun09=false;
84 ResourceUsage resu;
85 DecodeProc(narg-1, arg+1);
86 resu.Update();
87 cout << resu;
88 }
89 catch (PException& exc) {
90 cerr << " visfits2dt.cc catched PException " << exc.Msg() << endl;
91 rc = 77;
92 }
93 catch (std::exception& sex) {
94 cerr << "\n visfits2dt.cc std::exception :"
95 << (string)typeid(sex).name() << "\n msg= "
96 << sex.what() << endl;
97 rc = 78;
98 }
99 catch (...) {
100 cerr << " visfits2dt.cc catched unknown (...) exception " << endl;
101 rc = 79;
102 }
103
104 cout << ">>>> visfits2dt.cc ------- END ----------- RC=" << rc << endl;
105 return rc;
106
107}
108
109//--------------------------------------------------------------------
110// Traitement fichiers produits par vismfib (V Nov09)
111/* --Fonction-- */
112int DecodeProc(int narg, char* arg[])
113{
114 // Decodage des arguments et traitement
115 string inoutpath = arg[0];
116 int imin=0;
117 int imax=0;
118 int istep=1;
119 sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
120 int jf1=0;
121 int jf2=0;
122 int nfreq=0;
123 sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
124 int card=1;
125 vector<sa_size_t> rowlist;
126 if (narg>3) {
127 EnumeratedSequence es;
128 int nbad;
129 es.Append(arg[3], nbad, ",");
130 for(int j=0; j<es.Size(); j++) rowlist.push_back((int)es.Value(j));
131 }
132 // if (rowlist.size()<1) rowlist.push_back(0);
133
134 cout << " ----- visfits2dt/DecodeProc - Start - InOutPath= " << inoutpath << " IMin,Max,Step="
135 << imin << "," << imax << "," << istep << " Card=" << card << endl;
136 cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
137 cout << " RowList: " ;
138 for(int j=0; j<rowlist.size(); j++) cout << rowlist[j] << " , " ;
139 cout << endl;
140 int rc=ProcVisMtxFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, rowlist);
141 return rc;
142}
143
144// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
145/* --Fonction-- */
146int ProcVisMtxFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq,
147 vector<sa_size_t> rowlist)
148{
149 Timer tm("ProcVisMtxFiles");
150
151 DataTable visdt;
152 visdt.AddDoubleColumn("mfc");
153 visdt.AddDoubleColumn("mtt");
154 visdt.AddIntegerColumn("jfreq");
155 visdt.AddIntegerColumn("numch");
156 visdt.AddFloatColumn("vre");
157 visdt.AddFloatColumn("vim");
158
159 vector< TMatrix< r_4 > > vmtf;
160
161 if (jf1<1) jf1=1;
162 if ((jf2<1)||(jf2<jf1)) jf2=jf1;
163 if (nfreq<1) nfreq=1;
164 int djf=(jf2-jf1)/nfreq;
165 if (djf<1) djf=1;
166
167 cout << " --- ProcVisMtxFiles jf1=" << jf1 << " jf2= " << jf2 << " djf=" << djf << endl;
168 char fname[1024];
169
170 cout << " ---- ProcVisMtxFiles()-Begin - Reading chanum vector" << endl;
171 TVector< int_4 > chanum;
172 {
173 sprintf(fname, "%s/chanum.fits",inoutpath.c_str());
174 FitsInOutFile fic(fname, FitsInOutFile::Fits_RO);
175 fic.MoveAbsToHDU(3);
176 fic >> chanum;
177 }
178
179 cout << " ---- ProcVisMtxFiles()-Read chanum done " << endl;
180 sa_size_t nrows = (imax-imin+1)/istep;
181 sa_size_t kr=0;
182
183 sa_size_t mtf_binfreq=25;
184 sa_size_t mtf_bintime=5;
185
186 sa_size_t ncols=1;
187
188 for(int ifile=imin; ifile<=imax; ifile+=istep) {
189 sprintf(fname, "%s/vismtx%d.fits",inoutpath.c_str(),ifile);
190 cout << " ProcVisMtxFiles[" << ifile << "] opening file " << fname << endl;
191 FitsInOutFile fin(fname, FitsInOutFile::Fits_RO);
192 TMatrix< r_4 > vismtx;
193 fin >> vismtx;
194 // cout << vismtx.Info();
195
196 if ((ifile==imin)&&(rowlist.size()>0)) {
197 ncols = vismtx.NCols();
198 cout << " ProcVisMtxFiles/Info: Output Time-Frequency matrices NRows(time) "
199 << nrows/mtf_bintime+1 << " NCols(freq) =" << ncols/mtf_binfreq+1 << endl;
200 for(size_t j=0; j<rowlist.size(); j++)
201 vmtf.push_back(TMatrix< r_4 >(nrows/mtf_bintime+1, ncols/mtf_binfreq+1));
202 }
203
204 if (rowlist.size()>0) {
205 for(size_t j=0; j<rowlist.size(); j++)
206 for(sa_size_t icf=0; icf<ncols; icf++) {
207 vmtf[j](kr/mtf_bintime,icf/mtf_binfreq)+=vismtx(rowlist[j],icf);
208 }
209 }
210 // cout << " DBG* kr=" << kr << " kr/mtf_bintime=" << kr/mtf_bintime
211 // << " ncols/2/mtf_binfreq=" << ncols/2/mtf_binfreq << endl;
212 kr++;
213
214
215// Calcul moyenne dans des bandes en frequence
216 FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
217 }
218
219 if (rowlist.size()>0) {
220 sprintf(fname, "!%s/vistfreqmtx.fits",inoutpath.c_str());
221 cout << "ProcVisMtxFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;
222 FitsInOutFile fo(fname, FitsInOutFile::Fits_Create);
223 for(size_t j=0; j<rowlist.size(); j++)
224 fo << vmtf[j];
225 }
226 cout << visdt;
227 sprintf(fname, "!%s/visidt.ppf",inoutpath.c_str());
228 cout << "ProcVisMtxFiles: writing visibility datatable to file " << fname << endl;
229 FitsInOutFile fod(fname, FitsInOutFile::Fits_Create);
230 fod << visdt;
231
232 return 0;
233}
234
235/* --Fonction-- */
236int FillVisDTable(BaseDataTable& visdt_, TMatrix< r_4 > vismtx_, TVector< int_4> chanum_,
237 sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_)
238{
239 double xnt_[20];
240 xnt_[0]=(double)vismtx_.Info()["MEANFC"];
241 xnt_[1]=(double)vismtx_.Info()["MEANTT"]/1.25e8;
242
243 uint_4 nmean_=vismtx_.Info()["NPAQSUM"];
244 if (djf_<2) {
245 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
246 for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
247 xnt_[2]=jf;
248 xnt_[3]=chanum_(rv);
249 xnt_[4]=vismtx_(rv,jf*2)/(r_4)(nmean_);
250 xnt_[5]=vismtx_(rv,jf*2+1)/(r_4)(nmean_);
251 visdt_.AddRow(xnt_);
252 }
253 }
254 }
255 else {
256 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
257 for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
258 r_4 moyreal=0.;
259 r_4 moyimag=0.;
260 sa_size_t jjfmx=jf+djf_;
261 if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
262 for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
263 moyreal+=vismtx_(rv,jjf*2);
264 moyimag+=vismtx_(rv,jjf*2+1);
265 }
266 xnt_[2]=jf+djf_/2;
267 xnt_[3]=chanum_(rv);
268 xnt_[4]=moyreal/(r_4)(nmean_*djf_);
269 xnt_[5]=moyimag/(r_4)(nmean_*djf_);
270 visdt_.AddRow(xnt_);
271 }
272 }
273 }
274 return 0;
275}
276
277
278
279
Note: See TracBrowser for help on using the repository browser.