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

Last change on this file since 3903 was 3714, checked in by cmv, 16 years ago

protection contre djf=0, rz+cmv 11/12/2009

File size: 13.0 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 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
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"
28#include "datatable.h"
29#include "histinit.h"
30#include "matharr.h"
31#include "timestamp.h"
32#include <utilarr.h>
33
34// include sophya mesure ressource CPU/memoire ...
35#include "resusage.h"
36#include "ctimer.h"
37#include "timing.h"
38
39
40//--------------------------- Fonctions de ce fichier -------------------
41int Usage(bool fgshort=true);
42
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> > vismtx_, 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
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 {
86 bool fgvjun09=false;
87 if ((strcmp(arg[1],"-mcrd")==0)) fgvjun09=true;
88 ResourceUsage resu;
89 if (fgvjun09) DecodeProcVJun09(narg-2, arg+2);
90 else DecodeProc(narg-1, arg+1);
91 resu.Update();
92 cout << resu;
93 }
94 catch (PException& exc) {
95 cerr << " svv2mtx.cc catched PException " << exc.Msg() << endl;
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
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);
138
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}
148
149// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
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 DataTable visdt;
157 visdt.AddDoubleColumn("mfc");
158 visdt.AddDoubleColumn("mtt");
159 visdt.AddIntegerColumn("jfreq");
160 visdt.AddIntegerColumn("numch");
161 visdt.AddFloatColumn("vre");
162 visdt.AddFloatColumn("vim");
163
164 vector< TMatrix< complex<r_4> > > vmtf;
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=1;
171
172 char fname[1024];
173
174 cout << " ---- ProcSVFiles()-Begin - Reading chanum vector" << endl;
175 TVector< uint_4 > chanum;
176 {
177 sprintf(fname, "%s/chanum.ppf",inoutpath.c_str());
178 PInPersist pic(fname);
179 pic >> chanum;
180 }
181
182 cout << " ---- ProcSVFiles()-Read chanum done " << endl;
183 sa_size_t nrows = (imax-imin+1)/istep;
184 sa_size_t kr=0;
185
186 for(int ifile=imin; ifile<=imax; ifile+=istep) {
187 sprintf(fname, "%s/vismtx%d.ppf",inoutpath.c_str(),ifile);
188 cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
189 PInPersist pin(fname);
190 TMatrix< complex<r_4> > vismtx;
191 pin >> vismtx;
192
193 if (ifile==imin) {
194 sa_size_t ncols = vismtx.NCols();
195 cout << " ProcSVFilesVJun09/Info: Output Time-Frequency matrices NRows=NFiles"
196 << nrows << " NCols=NFreq=" << ncols << endl;
197 for(size_t j=0; j<rowlist.size(); j++)
198 vmtf.push_back(TMatrix< complex<r_4> >(nrows, ncols) );
199 }
200 for(size_t j=0; j<rowlist.size(); j++)
201 vmtf[j].Row(kr) = vismtx.Row(rowlist[j]);
202 kr++;
203
204// Calcul moyenne dans des bandes en frequence
205 FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
206 }
207
208 sprintf(fname, "%s/vistfreqmtx.ppf",inoutpath.c_str());
209 cout << "ProcSVFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;
210 POutPersist po(fname);
211 char tagb[64];
212 for(size_t j=0; j<rowlist.size(); j++) {
213 sprintf(tagb,"visft%d", chanum(rowlist[j]));
214 po << PPFNameTag(tagb) << vmtf[j];
215 }
216 cout << visdt;
217 sprintf(fname, "%s/svvdt.ppf",inoutpath.c_str());
218 cout << "ProcSVFile: writing visibility datatable to file " << fname << endl;
219 POutPersist pod(fname);
220 pod << visdt;
221
222 return 0;
223}
224
225/* --Fonction-- */
226int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_4> chanum_,
227 sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_)
228{
229 double xnt_[20];
230 xnt_[0]=(double)vismtx_.Info()["MeanFC"];
231 xnt_[1]=(double)vismtx_.Info()["MeanTT"]/1.25e8;
232
233 uint_4 nmean_=vismtx_.Info()["NPAQSUM"];
234 if (djf_<2) {
235 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
236 for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
237 xnt_[2]=jf;
238 xnt_[3]=chanum_(rv);
239 xnt_[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
240 xnt_[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
241 visdt_.AddRow(xnt_);
242 }
243 }
244 }
245 else {
246 for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
247 for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
248 r_4 moyreal=0.;
249 r_4 moyimag=0.;
250 sa_size_t jjfmx=jf+djf_;
251 if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
252 for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
253 moyreal+=vismtx_(rv,jjf).real();
254 moyimag+=vismtx_(rv,jjf).imag();
255 }
256 xnt_[2]=jf+djf_/2;
257 xnt_[3]=chanum_(rv);
258 xnt_[4]=moyreal/(r_4)(nmean_*djf_);
259 xnt_[5]=moyimag/(r_4)(nmean_*djf_);
260 visdt_.AddRow(xnt_);
261 }
262 }
263 }
264 return 0;
265}
266
267
268//--------------------------------------------------------------------
269// Traitement fichiers produits par mcrd (V Jun09)
270/* --Fonction-- */
271int DecodeProcVJun09(int narg, char* arg[])
272{
273 // Decodage des arguments et traitement
274 string inoutpath = arg[0];
275 int imin=0;
276 int imax=0;
277 int istep=1;
278 sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
279 int jf1=0;
280 int jf2=0;
281 int nfreq=0;
282 sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
283 int card=1;
284 if (narg>3) card=atoi(arg[3]);
285 cout << " ----- svv2mtx/DecodeProcVJun09 - Start - InOutPath= " << inoutpath << " IMin,Max,Step="
286 << imin << "," << imax << "," << istep << " Card=" << card << endl;
287 cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
288 int rc=ProcSVFilesVJun09(inoutpath, imin, imax, istep, jf1, jf2, nfreq, card);
289 return rc;
290}
291
292
293// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
294/* --Fonction-- */
295int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card)
296{
297 Timer tm("ProcSVFiles");
298 char fname[512];
299// NTuple
300 const char* nnames[10] = {"fcsm","ttsm","jfreq","s1","s2","s12","s12re","s12im","s12phi","s12mod"};
301 NTuple nt(10, nnames);
302 double xnt[15];
303 uint_4 nmnt = 0;
304 double ms1,ms2,ms12,ms12re,ms12im,ms12phi,ms12mod;
305
306 TMatrix<r_4> s1, s2;
307 TMatrix<r_4> v12re, v12im, v12phi,v12mod;
308 sa_size_t ncols = (imax-imin+1)/istep;
309 sa_size_t nrows = 10;
310 sa_size_t kc=0;
311 for(int ifile=imin; ifile<=imax; ifile+=istep) {
312 if (card==2)
313 sprintf(fname, "%s/Ch34_%d.ppf",inoutpath.c_str(),ifile);
314 else
315 sprintf(fname, "%s/Ch12_%d.ppf",inoutpath.c_str(),ifile);
316 cout << " ProcSVFilesVJun09[" << ifile << "] opening file " << fname << endl;
317 PInPersist pin(fname);
318 string tag1="specV1";
319 string tag2="specV2";
320 string tag12="visiV12";
321 if (card==2) {
322 tag1 = "specV3";
323 tag2 = "specV4";
324 tag12="visiV34";
325 }
326 TVector<r_4> sv1;
327 TVector<r_4> sv2;
328 TVector< complex<r_4> > vv12;
329 pin >> PPFNameTag(tag1) >> sv1;
330 pin >> PPFNameTag(tag2) >> sv2;
331 pin >> PPFNameTag(tag12) >> vv12;
332 if (ifile==imin) {
333 nrows = sv1.Size();
334 cout << " ProcSVFilesVJun09/Info: Output s1,s2 matrix size NRows=NFreq="
335 << nrows << " NCols=NFiles=" << ncols << endl;
336 s1.SetSize(nrows, ncols);
337 s2.SetSize(nrows, ncols);
338 nrows = vv12.Size();
339 cout << " ProcSVFilesVJun09/Info: Output v12 matrix size NRows=NFreq="
340 << nrows << " NCols=NFiles=" << ncols << endl;
341 v12re.SetSize(nrows, ncols);
342 v12im.SetSize(nrows, ncols);
343 v12phi.SetSize(nrows, ncols);
344 v12mod.SetSize(nrows, ncols);
345 }
346 s1.Column(kc) = sv1;
347 s2.Column(kc) = sv2;
348 v12re.Column(kc) = real(vv12);
349 v12im.Column(kc) = imag(vv12);
350 v12phi.Column(kc) = phase(vv12);
351 v12mod.Column(kc) = module(vv12);
352
353// Calcul moyenne dans des bandes en frequence
354 int deltajf=(jf2-jf1)/nfreq;
355 if (deltajf<1) deltajf=1;
356 for(int kf=0; kf<nfreq; kf++) {
357 sa_size_t jfstart=jf1+kf*deltajf;
358 sa_size_t jfend=jfstart+deltajf;
359 if (jfend>jf2) break;
360 nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=ms12mod=0.;
361 for(sa_size_t jf=jfstart; jf<jfend; jf++) {
362 ms1 += s1(jf,kc);
363 ms2 += s2(jf,kc);
364 ms12re += v12re(jf,kc);
365 ms12im += v12im(jf,kc);
366 ms12phi += v12phi(jf,kc);
367 ms12mod += v12mod(jf,kc);
368 }
369 nmnt = (jfend-jfstart);
370 if (nmnt>0) {
371 double fnorm = (double)nmnt;
372 xnt[0] = ((int_8)(sv1.Info()["StartFC"])+(int_8)(sv1.Info()["EndFC"]))*0.5;
373 xnt[1] = ((int_8)(sv1.Info()["StartTT"])+(int_8)(sv1.Info()["EndTT"]))*0.5;
374 xnt[2] = kf;
375 xnt[3] = ms1/fnorm;
376 xnt[4] = ms2/fnorm;
377 xnt[5] = ms12/fnorm;
378 xnt[6] = ms12re/fnorm;
379 xnt[7] = ms12im/fnorm;
380 xnt[8] = ms12phi/fnorm;
381 xnt[9] = ms12mod/fnorm;
382 nt.Fill(xnt);
383 }
384 }
385 kc++;
386
387 }
388 if (card==2)
389 sprintf(fname, "%s/Ch34mtx.ppf",inoutpath.c_str());
390 else
391 sprintf(fname, "%s/Ch12mtx.ppf",inoutpath.c_str());
392
393 cout << nt;
394 cout << "ProcSVFilesVJun09: Opening file " << fname << " for writing" << endl;
395 POutPersist po(fname);
396 string tag1="s1";
397 string tag2="s2";
398 string tag12r="v12re";
399 string tag12i="v12im";
400 string tag12p="v12phi";
401 string tagnt="nt12";
402 if (card==2) {
403 tag1="s3";
404 tag2="s4";
405 tag12r="v34re";
406 tag12i="v34im";
407 tag12p="v34phi";
408 tagnt="nt34";
409 }
410 po << PPFNameTag(tag1) << s1;
411 po << PPFNameTag(tag2) << s2;
412 po << PPFNameTag(tag12r) << v12re;
413 po << PPFNameTag(tag12i) << v12im;
414 po << PPFNameTag(tag12p) << v12phi;
415 po << PPFNameTag(tagnt) << nt;
416 cout << "ProcSVFilesVJun09: Matrices s1, s2, v12re, v12im, v12phi, NTuple nt written to file " << fname << endl;
417 return 0;
418}
419
420
Note: See TracBrowser for help on using the repository browser.