source: Sophya/trunk/AddOn/TAcq/svv2mtx.cc@ 3702

Last change on this file since 3702 was 3702, checked in by ansari, 16 years ago

Correction bug, Reza 4/12/09

File size: 12.8 KB
RevLine 
[3646]1// Utilisation de SOPHYA pour faciliter les tests ...
2#include "sopnamsp.h"
3#include "machdefs.h"
4
5/* ----------------------------------------------------------
[3701]6 Projet BAORadio - (C) LAL/IRFU 2008-2010
7
8 Programme de lecture des fichiers vecteurs/matrices de
9 visibilites produits par mcrd / vismfib
10 R. Ansari, C. Magneville - LAL/Irfu
[3646]11 V : Mai 2009
12 ---------------------------------------------------------- */
13
14// include standard c/c++
15#include <math.h>
16#include <stdio.h>
17#include <stdlib.h>
18#include <string.h>
19
20#include <iostream>
21#include <string>
22
23#include "pexceptions.h"
24#include "tvector.h"
25#include "fioarr.h"
26// #include "tarrinit.h"
27#include "ntuple.h"
[3701]28#include "datatable.h"
[3646]29#include "histinit.h"
30#include "matharr.h"
31#include "timestamp.h"
[3701]32#include <utilarr.h>
[3646]33
34// include sophya mesure ressource CPU/memoire ...
35#include "resusage.h"
36#include "ctimer.h"
37#include "timing.h"
38
39
[3701]40//--------------------------- Fonctions de ce fichier -------------------
[3646]41int Usage(bool fgshort=true);
42
[3701]43int DecodeProc(int narg, char* arg[]);
44int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq,
45 vector<sa_size_t> tfrlist);
46int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > visdt_, TVector< uint_4> chanum_,
47 sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) ;
48
49int DecodeProcVJun09(int narg, char* arg[]);
50int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card);
51//------------------------------------------------------------------------------------------------------------
52
53/* --Fonction-- */
54int Usage(bool fgshort)
55{
56 cout << " --- svv2mtx.cc : Read PPF files produced by mcrd/visfmib to make time-frequency\n"
57 << " matrices and Visibilites=f(time) datatable " << endl;
58 cout << " Usage: svv2mtx -mcrd InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [card=1] \n"
59 << " OR \n"
60 << " svv2mtx InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq VisMtxRowList" << endl;
61 if (fgshort) {
62 cout << " svv2mtx -h for detailed instructions" << endl;
63 return 1;
64 }
65 cout << " -mcrd : Read mcrd output files -> ProcSVFilesVJun09() \n"
66 << " InOutPath : Input/Output directory name \n"
67 << " Imin,Imax,IStep: Input PPF files sequence number \n"
68 << " FileNames=InOutPath/Ch12_II.fits Imin<=II<=Imax II+=IStep \n"
69 << " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple\n"
70 << " VisMtxRowList : List of visibiliy matrix rows (0,2,...) -> time-freq map \n"
71 << " card=1 Ch12 , card=2 Ch34 (mcrd output files) " << endl;
72 return 1;
73}
74
[3646]75//----------------------------------------------------
76//----------------------------------------------------
77int main(int narg, char* arg[])
78{
79 if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
80 if (narg<4) return Usage(true);
81 HiStatsInitiator _inia;
82 // TArrayInitiator _inia;
83
84 int rc = 0;
85 try {
[3701]86 bool fgvjun09=false;
87 if ((strcmp(arg[1],"-mcrd")==0)) fgvjun09=true;
[3646]88 ResourceUsage resu;
[3701]89 if (fgvjun09) DecodeProcVJun09(narg-2, arg+2);
90 else DecodeProc(narg-1, arg+1);
[3646]91 resu.Update();
92 cout << resu;
93 }
94 catch (PException& exc) {
[3702]95 cerr << " svv2mtx.cc catched PException " << exc.Msg() << endl;
[3646]96 rc = 77;
97 }
98 catch (std::exception& sex) {
99 cerr << "\n svv2mtx.cc std::exception :"
100 << (string)typeid(sex).name() << "\n msg= "
101 << sex.what() << endl;
102 rc = 78;
103 }
104 catch (...) {
105 cerr << " svv2mtx.cc catched unknown (...) exception " << endl;
106 rc = 79;
107 }
108
109 cout << ">>>> svv2mtx.cc ------- END ----------- RC=" << rc << endl;
110 return rc;
111
112}
113
[3701]114//--------------------------------------------------------------------
115// Traitement fichiers produits par vismfib (V Nov09)
116/* --Fonction-- */
117int DecodeProc(int narg, char* arg[])
118{
119 // Decodage des arguments et traitement
120 string inoutpath = arg[0];
121 int imin=0;
122 int imax=0;
123 int istep=1;
124 sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
125 int jf1=0;
126 int jf2=0;
127 int nfreq=0;
128 sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
129 int card=1;
130 vector<sa_size_t> rowlist;
131 if (narg>3) {
132 EnumeratedSequence es;
133 int nbad;
134 es.Append(arg[3], nbad, ",");
135 for(int j=0; j<es.Size(); j++) rowlist.push_back((int)es.Value(j));
136 }
137 if (rowlist.size()<1) rowlist.push_back(0);
[3646]138
[3701]139 cout << " ----- svv2mtx/DecodeProc - Start - InOutPath= " << inoutpath << " IMin,Max,Step="
140 << imin << "," << imax << "," << istep << " Card=" << card << endl;
141 cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
142 cout << " RowList: " ;
143 for(int j=0; j<rowlist.size(); j++) cout << rowlist[j] << " , " ;
144 cout << endl;
145 int rc=ProcSVFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, rowlist);
146 return rc;
147}
[3646]148
149// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
[3701]150/* --Fonction-- */
151int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq,
152 vector<sa_size_t> rowlist)
153{
154 Timer tm("ProcSVFiles");
155
156 vector< TMatrix< complex<r_4> > > vmtf(rowlist.size());
157
158 DataTable visdt;
159 visdt.AddDoubleColumn("mfc");
160 visdt.AddDoubleColumn("mtt");
161 visdt.AddIntegerColumn("jfreq");
162 visdt.AddIntegerColumn("numch");
163 visdt.AddFloatColumn("vre");
164 visdt.AddFloatColumn("vim");
165
166 if (jf1<1) jf1=1;
167 if ((jf2<1)||(jf2<jf1)) jf2=jf1;
168 if (nfreq<1) nfreq=1;
169 int djf=(jf2-jf1)/nfreq;
170 if (djf<1) djf=0;
171
172 char fname[1024];
173
174 TVector< uint_4 > chanum;
175 sprintf(fname, "%s/chanum.ppf",inoutpath.c_str());
176 {
177 PInPersist pic(fname);
178 pic >> chanum;
179 }
180
181 sa_size_t nrows = (imax-imin+1)/istep;\
182 sa_size_t kr=0;
183
184 for(int ifile=imin; ifile<=imax; ifile+=istep) {
185 sprintf(fname, "%s/vismtx%d.ppf",inoutpath.c_str(),ifile);
186 cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
187 PInPersist pin(fname);
188 TMatrix< complex<r_4> > vismtx;
189 pin >> vismtx;
190
191 if (ifile==imin) {
192 sa_size_t ncols = vismtx.NCols();
193 cout << " ProcSVFilesVJun09/Info: Output Time-Frequency matrices NRows=NFiles"
194 << nrows << " NCols=NFreq=" << ncols << endl;
195 for(size_t j=0; j<rowlist.size(); j++) vmtf[j].SetSize(nrows, ncols);
196 }
197 for(size_t j=0; j<rowlist.size(); j++)
198 vmtf[j].Row(kr) = vismtx.Row(rowlist[j]);
199 kr++;
200
201// Calcul moyenne dans des bandes en frequence
202 FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
203 }
204
205 sprintf(fname, "%s/vistfreqmtx.ppf",inoutpath.c_str());
206 cout << "ProcSVFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;
207 POutPersist po(fname);
208 char tagb[64];
209 for(size_t j=0; j<rowlist.size(); j++) {
210 sprintf(tagb,"visft%d", chanum(rowlist[j]));
211 po << PPFNameTag(tagb) << vmtf[j];
212 }
213 cout << visdt;
214 sprintf(fname, "%s/svvdt.ppf",inoutpath.c_str());
215 cout << "ProcSVFile: writing visibility datatable to file " << fname << endl;
216 POutPersist pod(fname);
217 pod << visdt;
218
219 return 0;
220}
221
222/* --Fonction-- */
223int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_4> chanum_,
224 sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_)
225{
226 double xnt_[20];
227 xnt_[0]=(double)vismtx_.Info()["MeanFC"];
228 xnt_[1]=(double)vismtx_.Info()["MeanTT"]/1.25e8;
229
230 uint_4 nmean_=vismtx_.Info()["NPAQSUM"];
231 if (djf_<2) {
232 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
233 for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
234 xnt_[2]=jf;
235 xnt_[3]=chanum_(rv);
236 xnt_[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
237 xnt_[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
238 visdt_.AddRow(xnt_);
239 }
240 }
241 }
242 else {
243 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
244 for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
245 r_4 moyreal=0.;
246 r_4 moyimag=0.;
247 for(sa_size_t jjf=jf; jjf<jf+djf_; jjf++) {
248 moyreal+=vismtx_(rv,jjf).real();
249 moyimag+=vismtx_(rv,jjf).imag();
250 }
251 xnt_[2]=jf+djf_/2;
252 xnt_[3]=chanum_(rv);
253 xnt_[4]=moyreal/(r_4)(nmean_*djf_);
254 xnt_[5]=moyimag/(r_4)(nmean_*djf_);
255 visdt_.AddRow(xnt_);
256 }
257 }
258 }
259 return 0;
260}
261
262
263//--------------------------------------------------------------------
264// Traitement fichiers produits par mcrd (V Jun09)
265/* --Fonction-- */
266int DecodeProcVJun09(int narg, char* arg[])
267{
268 // Decodage des arguments et traitement
269 string inoutpath = arg[0];
270 int imin=0;
271 int imax=0;
272 int istep=1;
273 sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
274 int jf1=0;
275 int jf2=0;
276 int nfreq=0;
277 sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
278 int card=1;
279 if (narg>3) card=atoi(arg[3]);
280 cout << " ----- svv2mtx/DecodeProcVJun09 - Start - InOutPath= " << inoutpath << " IMin,Max,Step="
281 << imin << "," << imax << "," << istep << " Card=" << card << endl;
282 cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
283 int rc=ProcSVFilesVJun09(inoutpath, imin, imax, istep, jf1, jf2, nfreq, card);
284 return rc;
285}
286
287
288// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
289/* --Fonction-- */
[3699]290int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card)
[3646]291{
292 Timer tm("ProcSVFiles");
293 char fname[512];
294// NTuple
[3647]295 const char* nnames[10] = {"fcsm","ttsm","jfreq","s1","s2","s12","s12re","s12im","s12phi","s12mod"};
296 NTuple nt(10, nnames);
297 double xnt[15];
[3646]298 uint_4 nmnt = 0;
[3647]299 double ms1,ms2,ms12,ms12re,ms12im,ms12phi,ms12mod;
[3646]300
301 TMatrix<r_4> s1, s2;
[3647]302 TMatrix<r_4> v12re, v12im, v12phi,v12mod;
[3646]303 sa_size_t ncols = (imax-imin+1)/istep;
304 sa_size_t nrows = 10;
305 sa_size_t kc=0;
306 for(int ifile=imin; ifile<=imax; ifile+=istep) {
307 if (card==2)
308 sprintf(fname, "%s/Ch34_%d.ppf",inoutpath.c_str(),ifile);
309 else
310 sprintf(fname, "%s/Ch12_%d.ppf",inoutpath.c_str(),ifile);
[3701]311 cout << " ProcSVFilesVJun09[" << ifile << "] opening file " << fname << endl;
[3646]312 PInPersist pin(fname);
313 string tag1="specV1";
314 string tag2="specV2";
315 string tag12="visiV12";
316 if (card==2) {
317 tag1 = "specV3";
318 tag2 = "specV4";
319 tag12="visiV34";
320 }
321 TVector<r_4> sv1;
322 TVector<r_4> sv2;
323 TVector< complex<r_4> > vv12;
324 pin >> PPFNameTag(tag1) >> sv1;
325 pin >> PPFNameTag(tag2) >> sv2;
326 pin >> PPFNameTag(tag12) >> vv12;
327 if (ifile==imin) {
328 nrows = sv1.Size();
[3699]329 cout << " ProcSVFilesVJun09/Info: Output s1,s2 matrix size NRows=NFreq="
[3646]330 << nrows << " NCols=NFiles=" << ncols << endl;
331 s1.SetSize(nrows, ncols);
332 s2.SetSize(nrows, ncols);
333 nrows = vv12.Size();
[3699]334 cout << " ProcSVFilesVJun09/Info: Output v12 matrix size NRows=NFreq="
[3646]335 << nrows << " NCols=NFiles=" << ncols << endl;
336 v12re.SetSize(nrows, ncols);
337 v12im.SetSize(nrows, ncols);
338 v12phi.SetSize(nrows, ncols);
[3647]339 v12mod.SetSize(nrows, ncols);
[3646]340 }
341 s1.Column(kc) = sv1;
342 s2.Column(kc) = sv2;
343 v12re.Column(kc) = real(vv12);
344 v12im.Column(kc) = imag(vv12);
345 v12phi.Column(kc) = phase(vv12);
[3647]346 v12mod.Column(kc) = module(vv12);
347
348// Calcul moyenne dans des bandes en frequence
349 int deltajf=(jf2-jf1)/nfreq;
350 if (deltajf<1) deltajf=1;
351 for(int kf=0; kf<nfreq; kf++) {
352 sa_size_t jfstart=jf1+kf*deltajf;
353 sa_size_t jfend=jfstart+deltajf;
354 if (jfend>jf2) break;
355 nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=ms12mod=0.;
356 for(sa_size_t jf=jfstart; jf<jfend; jf++) {
357 ms1 += s1(jf,kc);
358 ms2 += s2(jf,kc);
359 ms12re += v12re(jf,kc);
360 ms12im += v12im(jf,kc);
361 ms12phi += v12phi(jf,kc);
362 ms12mod += v12mod(jf,kc);
363 }
364 nmnt = (jfend-jfstart);
365 if (nmnt>0) {
366 double fnorm = (double)nmnt;
367 xnt[0] = ((int_8)(sv1.Info()["StartFC"])+(int_8)(sv1.Info()["EndFC"]))*0.5;
368 xnt[1] = ((int_8)(sv1.Info()["StartTT"])+(int_8)(sv1.Info()["EndTT"]))*0.5;
369 xnt[2] = kf;
370 xnt[3] = ms1/fnorm;
371 xnt[4] = ms2/fnorm;
372 xnt[5] = ms12/fnorm;
373 xnt[6] = ms12re/fnorm;
374 xnt[7] = ms12im/fnorm;
375 xnt[8] = ms12phi/fnorm;
376 xnt[9] = ms12mod/fnorm;
377 nt.Fill(xnt);
378 }
[3646]379 }
380 kc++;
381
382 }
383 if (card==2)
384 sprintf(fname, "%s/Ch34mtx.ppf",inoutpath.c_str());
385 else
386 sprintf(fname, "%s/Ch12mtx.ppf",inoutpath.c_str());
387
388 cout << nt;
[3699]389 cout << "ProcSVFilesVJun09: Opening file " << fname << " for writing" << endl;
[3646]390 POutPersist po(fname);
391 string tag1="s1";
392 string tag2="s2";
393 string tag12r="v12re";
394 string tag12i="v12im";
395 string tag12p="v12phi";
396 string tagnt="nt12";
397 if (card==2) {
398 tag1="s3";
399 tag2="s4";
400 tag12r="v34re";
401 tag12i="v34im";
402 tag12p="v34phi";
403 tagnt="nt34";
404 }
405 po << PPFNameTag(tag1) << s1;
406 po << PPFNameTag(tag2) << s2;
407 po << PPFNameTag(tag12r) << v12re;
408 po << PPFNameTag(tag12i) << v12im;
409 po << PPFNameTag(tag12p) << v12phi;
410 po << PPFNameTag(tagnt) << nt;
[3699]411 cout << "ProcSVFilesVJun09: Matrices s1, s2, v12re, v12im, v12phi, NTuple nt written to file " << fname << endl;
[3646]412 return 0;
413}
414
415
Note: See TracBrowser for help on using the repository browser.