source: Sophya/trunk/AddOn/TAcq/brproc.cc@ 3646

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

Ajout du programme svv2mtx.cc de lecture des fichiers de mcrd.cc pour faire des tableaux 2D (matrices) temps-frequence + corrections mineures, Reza 09/06/2009
File size: 15.9 KB
Line 
1#include "racquproc.h"
2
3#include <stdlib.h>
4#include <string.h>
5#include <unistd.h>
6#include <fstream>
7#include <signal.h>
8
9#include "pexceptions.h"
10#include "tvector.h"
11#include "ntuple.h"
12#include "fioarr.h"
13#include "timestamp.h"
14#include "ctimer.h"
15#include "fftpserver.h"
16#include "fftwserver.h"
17
18#include "FFTW/fftw3.h"
19
20
21#include "pciewrap.h"
22#include "brpaqu.h"
23#include "brproc.h"
24
25
26//---------------------------------------------------------------
27// Classe thread de traitement donnees ADC avec 2 voies par frame
28//---------------------------------------------------------------
29
30/* --Methode-- */
31BRProcARaw2C::BRProcARaw2C(RAcqMemZoneMgr& mem, string& path, uint_4 nmean,
32 uint_4 nmax, bool fgnotrl, int card)
33 : memgr(mem)
34{
35 nmax_ = nmax;
36 nmean_ = nmean;
37 stop_ = false;
38 path_ = path;
39 fgnotrl_ = fgnotrl;
40 card_ = card;
41}
42
43/* --Methode-- */
44void BRProcARaw2C::Stop()
45{
46 stop_=true;
47 // cout <<" BRProcARaw2C::Stop ... > STOP " << endl;
48}
49
50
51static inline r_4 Zmod2(complex<r_4> z)
52{ return (z.real()*z.real()+z.imag()*z.imag()); }
53
54static inline string card2name_(int card)
55{
56 if (card==2) return " (Chan3,4) ";
57 else return " (Chan1,2) ";
58}
59/* --Methode-- */
60void BRProcARaw2C::run()
61{
62 setRC(1);
63 try {
64 Timer tm("BRProcARaw2C", false);
65 TimeStamp ts;
66 BRPaqChecker pcheck(!fgnotrl_); // Verification/comptage des paquets
67
68 size_t totnbytesout = 0;
69 size_t totnbytesproc = 0;
70
71 cout << " BRProcARaw2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
72 << " NMean=" << nmean_ << card2name_(card_) << endl;
73 cout << " BRProcARaw2C::run()... - Output Data Path: " << path_ << endl;
74 char fname[512];
75// sprintf(fname,"%s/proc.log",path_.c_str());
76// ofstream filog(fname);
77// filog << " BRProcARaw2C::run() - starting log file " << ts << endl;
78// filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl;
79
80// NTuple
81 const char* nnames[8] = {"fcs","tts","s1","s2","s12","s12re","s12im","s12phi"};
82 NTuple nt(8, nnames);
83 double xnt[10];
84 uint_4 nmnt = 0;
85 double ms1,ms2,ms12,ms12re,ms12im,ms12phi;
86// Initialisation pour calcul FFT
87 TVector< complex<r_4> > cfour1; // composant TF
88 uint_4 paqsz = memgr.PaqSize();
89 uint_4 procpaqsz = memgr.ProcPaqSize();
90
91
92 BRPaquet pq(NULL, NULL, paqsz);
93 TVector<r_4> vx(pq.DataSize()/2);
94 vx = (r_4)(0.);
95 FFTPackServer ffts;
96 ffts.FFTForward(vx, cfour1);
97// cfour1.SetSize((paqsz-40)/2+1);
98 TVector< complex<r_4> > cfour2(cfour1.Size());
99
100 TVector<r_4> spectreV1(cfour1.Size());
101 TVector<r_4> spectreV2(cfour1.Size());
102 TVector< complex<r_4> > visiV12( cfour1.Size() );
103
104 cout << " *DBG*BRProcARaw2C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
105 << " procpaqsz/2=" << procpaqsz/2 << " cfour1.Size()=" << cfour1.Size()
106 << " *8=" << cfour1.Size()*8 << endl;
107
108 fftwf_plan plan1 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(),
109 (fftwf_complex*)cfour1.Data(), FFTW_ESTIMATE);
110 fftwf_plan plan2 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(),
111 (fftwf_complex*)cfour2.Data(), FFTW_ESTIMATE);
112
113 uint_4 ifile = 0;
114 uint_4 nzm = 0;
115 for (uint_4 kmz=0; kmz<nmax_; kmz++) {
116 if (stop_) break;
117 int mid = memgr.FindMemZoneId(MemZA_ProcA);
118 Byte* buff = memgr.GetMemZone(mid);
119 if (buff == NULL) {
120 cout << " BRProcARaw2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
121 break;
122 }
123 Byte* procbuff = memgr.GetProcMemZone(mid);
124 if (procbuff == NULL) {
125 cout << " BRProcARaw2C::run()/ERROR memgr.GetProcMemZone(" << mid << ") -> NULL" << endl;
126 break;
127 }
128 nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=0.;
129 for(uint_4 i=0; i<memgr.NbPaquets(); i++) {
130 BRPaquet paq(NULL, buff+i*paqsz, paqsz);
131 if (!pcheck.Check(paq)) continue; // on ne traite que les paquets OK
132
133// Traitement voie 1
134 for(sa_size_t j=0; j<vx.Size(); j++)
135 vx(j) = (r_4)(*(paq.Data1()+j))-127.5;
136// fftwf_complex* coeff1 = (fftwf_complex*)(procbuff+i*procpaqsz);
137 fftwf_execute(plan1);
138// complex<r_4>* zp1 = (complex<r_4>*)(vx.Data());
139// ffts.FFTForward(vx, cfour1);
140 for(sa_size_t j=0; j<spectreV1.Size(); j++)
141 spectreV1(j) += Zmod2(cfour1(j));
142 memcpy(procbuff+i*procpaqsz, cfour1.Data(), sizeof(complex<r_4>)*cfour1.Size());
143// Traitement voie 2
144 for(sa_size_t j=0; j<vx.Size(); j++)
145 vx(j) = (r_4)(*(paq.Data2()+j))-127.5;
146
147 fftwf_execute(plan2);
148 for(sa_size_t j=0; j<spectreV2.Size(); j++)
149 spectreV2(j) += Zmod2(cfour2(j)); // Zmod2(zp2[j]);
150 memcpy(procbuff+i*procpaqsz+procpaqsz/2, cfour2.Data(), sizeof(complex<r_4>)*cfour2.Size());
151
152// Calcul correlation (visibilite V1 * V2)
153 for(sa_size_t j=0; j<visiV12.Size(); j++)
154 visiV12(j)+=cfour1(j)*conj(cfour2(j));
155// for(sa_size_t j=0; j<visiV12.Size(); j++) visiV12(j)+=zp1[j]*zp2[j];
156 nzm++;
157 if (nmnt==0) { xnt[0]=paq.FrameCounter(); xnt[1]=paq.TimeTag(); }
158 for(sa_size_t j=2700; j<2800; j++) {
159 ms1 += Zmod2(cfour1(j)); ms2 += Zmod2(cfour2(j));
160 complex<r_4> zvis = cfour1(j)*conj(cfour2(j));
161 ms12 += Zmod2(zvis); ms12re += zvis.real(); ms12im += zvis.imag();
162 ms12phi+= atan2(zvis.imag(),zvis.real());
163 }
164 nmnt++;
165 totnbytesproc += paq.DataSize();
166 totnbytesout += (2*sizeof(complex<r_4>)*cfour1.Size());
167
168 } // Fin de boucle sur les paquets d'une zone
169 if (nmnt>0) {
170 double fnorm = (double)nmnt*(2800-2700);
171 xnt[2] = ms1 /= fnorm;
172 xnt[3] = ms2 /= fnorm;
173 xnt[4] = ms12 /= fnorm;
174 xnt[5] = ms12re /= fnorm;
175 xnt[6] = ms12im /= fnorm;
176 xnt[7] = ms12phi /= fnorm;
177 nt.Fill(xnt);
178 }
179 if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
180 spectreV1 /= (r_4)(nzm);
181 spectreV2 /= (r_4)(nzm);
182
183 visiV12 /= complex<r_4>((r_4)nzm, 0.);
184
185 spectreV1.Info()["NPaqMoy"] = nzm;
186 spectreV2.Info()["NPaqMoy"] = nzm;
187 visiV12.Info()["NPaqMoy"] = nzm;
188 {
189 sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
190 POutPersist po(fname);
191 string tag1="specV1";
192 string tag2="specV2";
193 string tag12="visiV12";
194 if (card_==2) {
195 tag1 = "specV3";
196 tag2 = "specV4";
197 tag12="visiV34";
198 }
199 po << PPFNameTag(tag1) << spectreV1;
200 po << PPFNameTag(tag2) << spectreV2;
201 po << PPFNameTag(tag12) << visiV12;
202 }
203 spectreV1 = (r_4)(0.);
204 spectreV2 = (r_4)(0.);
205 visiV12 = complex<r_4>(0., 0.);
206 nzm = 0; ifile++;
207// ts.SetNow();
208// filog << ts << " : proc file " << fname << endl;
209 cout << " BRProcARaw2C::run() created file " << fname << card2name_(card_) << endl;
210 }
211
212 memgr.FreeMemZone(mid, MemZS_ProcA);
213 } // Fin de boucle sur les zones a traiter
214 cout << " ------------ BRProcARaw2C::run() END " << card2name_(card_)
215 << " ------------ " << endl;
216 {
217 cout << nt;
218 sprintf(fname,"%s_nt.ppf",path_.c_str());
219 POutPersist po(fname);
220 po << PPFNameTag("ntv12") << nt;
221 cout << " BRProcARaw2C::run() created NTuple file " << fname << card2name_(card_) << endl;
222 }
223 ts.SetNow();
224 tm.SplitQ();
225 cout << " TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= "
226 << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s"
227 << " ProcDataOut=" << totnbytesout/(1024*1024) << " MB" << endl;
228 cout << pcheck;
229 cout << " BRProcARaw2C::run()/Timing: " << card2name_(card_) << endl;
230 tm.Print();
231 cout << " ---------------------------------------------------------- " << endl;
232
233 }
234 catch (PException& exc) {
235 cout << " BRProcARaw2C::run()/catched PException " << exc.Msg() << endl;
236 setRC(3);
237 return;
238 }
239 catch(...) {
240 cout << " BRProcARaw2C::run()/catched unknown ... exception " << endl;
241 setRC(4);
242 return;
243 }
244 setRC(0);
245 return;
246}
247
248//---------------------------------------------------------------------
249// Classe thread de traitement 2 x 2 voies/frames (Apres BRProcARaw2C)
250//---------------------------------------------------------------------
251
252/* --Methode-- */
253BRProcBRaw4C::BRProcBRaw4C(RAcqMemZoneMgr& mem1, RAcqMemZoneMgr& mem2,
254 string& path, uint_4 nmean, uint_4 nmax, bool fgnotrl)
255 : memgr1(mem1), memgr2(mem2)
256{
257 nmax_ = nmax;
258 nmean_ = nmean;
259 stop_ = false;
260 path_ = path;
261 fgnotrl_ = fgnotrl;
262}
263
264/* --Methode-- */
265void BRProcBRaw4C::Stop()
266{
267 stop_=true;
268 // cout <<" BRProcBRaw4C::Stop ... > STOP " << endl;
269}
270
271
272/* --Methode-- */
273void BRProcBRaw4C::run()
274{
275 setRC(1);
276 try {
277 Timer tm("BRProcBRaw4C", false);
278 TimeStamp ts;
279 BRPaqChecker pcheck1(!fgnotrl_); // Verification/comptage des paquets
280 BRPaqChecker pcheck2(!fgnotrl_); // Verification/comptage des paquets
281
282 size_t totnbytesout = 0;
283 size_t totnbytesproc = 0;
284
285 cout << " BRProcBRaw4C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
286 << " NMean=" << nmean_ << endl;
287 cout << " BRProcBRaw4C::run()... - Output Data Path: " << path_ << endl;
288
289 uint_4 paqsz = memgr1.PaqSize();
290 uint_4 procpaqsz = memgr1.ProcPaqSize();
291 if ((paqsz != memgr2.PaqSize())||(procpaqsz!= memgr2.ProcPaqSize())) {
292 cout << "BRProcBRaw4C::run()/ERROR : different paquet size -> stop \n ...(PaqSz1="
293 << paqsz << " Sz2=" << memgr2.PaqSize() << " ProcPaqSz1="
294 << procpaqsz << " Sz2=" << memgr2.ProcPaqSize() << " )" << endl;
295 setRC(9);
296 return;
297 }
298
299 TVector< complex<r_4> > cfour; // composant TF
300/*
301 BRPaquet pq(NULL, NULL, paqsz);
302 TVector<r_4> vx(pq.DataSize()/2);
303 vx = (r_4)(0.);
304 FFTPackServer ffts;
305 ffts.FFTForward(vx, cfour);
306
307 TVector< complex<r_4> > visiV13( cfour.Size() );
308 TVector< complex<r_4> > visiV14( cfour.Size() );
309 TVector< complex<r_4> > visiV23( cfour.Size() );
310 TVector< complex<r_4> > visiV24( cfour.Size() );
311*/
312 int szfour = (paqsz-40)/2+1;
313 TVector< complex<r_4> > visiV13( szfour );
314 TVector< complex<r_4> > visiV14( szfour );
315 TVector< complex<r_4> > visiV23( szfour );
316 TVector< complex<r_4> > visiV24( szfour );
317 // cout << " *DBG*AAAAA ---- Vectors OK" << endl;
318 cout << " *DBG*BRProcBRaw4C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
319 << " procpaqsz/2=" << procpaqsz/2 << " cfour.Size()=" << szfour
320 << " *8=" << szfour*8 << endl;
321
322 uint_4 nzm = 0;
323 uint_4 totnoksfc = 0;
324 uint_4 totnokpaq = 0;
325 uint_4 totnpaq = 0;
326 uint_4 ifile = 0;
327 for (uint_4 kmz=0; kmz<nmax_; kmz++) {
328 uint_4 noksfc = 0;
329 uint_4 nokpaq = 0;
330 if (stop_) break;
331 // cout << " *DBG*BBBBB" << kmz << endl;
332
333 int mid1 = memgr1.FindMemZoneId(MemZA_ProcB);
334 Byte* buff1 = memgr1.GetMemZone(mid1);
335 if (buff1 == NULL) {
336 cout << " BRProcBRaw4C::run()/ERROR memgr.GetMemZone(" << mid1 << ") -> NULL" << endl;
337 break;
338 }
339 Byte* procbuff1 = memgr1.GetProcMemZone(mid1);
340 if (procbuff1 == NULL) {
341 cout << " BRProcBRaw4C::run()/ERROR memgr.GetProcMemZone(" << mid1 << ") -> NULL" << endl;
342 break;
343 }
344 int mid2 = memgr2.FindMemZoneId(MemZA_ProcB);
345 Byte* buff2 = memgr2.GetMemZone(mid2);
346 if (buff1 == NULL) {
347 cout << " BRProcBRaw4C::run()/ERROR memgr.GetMemZone(" << mid2 << ") -> NULL" << endl;
348 break;
349 }
350 Byte* procbuff2 = memgr2.GetProcMemZone(mid2);
351 if (procbuff2 == NULL) {
352 cout << " BRProcBRaw4C::run()/ERROR memgr.GetProcMemZone(" << mid2 << ") -> NULL" << endl;
353 break;
354 }
355 uint_4 i1,i2;
356 i1=i2=0;
357// cout << " *DBG*CCCCCC " << kmz << " memgr1.NbPaquets() =" << memgr1.NbPaquets() << endl;
358 while((i1<memgr1.NbPaquets())&&(i2<memgr2.NbPaquets())) {
359 BRPaquet paq1(NULL, buff1+i1*paqsz, paqsz);
360 BRPaquet paq2(NULL, buff2+i2*paqsz, paqsz);
361 totnpaq++;
362// cout << " DBG["<<kmz<<"] i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
363//<<","<<paq2.FrameCounter()<<endl;
364 // on ne traite que les paquets OK
365 if (!pcheck1.Check(paq1)) { i1++; continue; }
366 if (!pcheck2.Check(paq2)) { i2++; continue; }
367 nokpaq++;
368 if (paq1.FrameCounter()<paq2.FrameCounter()) { i1++; continue; }
369 if (paq2.FrameCounter()<paq1.FrameCounter()) { i2++; continue; }
370// cout << " DBG["<<kmz<<"]OKOK i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
371// <<","<<paq2.FrameCounter()<<endl;
372
373 if ((i1>=memgr1.NbPaquets())||(i2>=memgr1.NbPaquets())) {
374 cout << " *BUG*BUG i1=" << i1 << " i2=" << i2 << endl;
375 break;
376 }
377 // Les deux framecounters sont identiques ...
378 noksfc++;
379 complex<r_4>* zp1 = (complex<r_4>*)(procbuff1+i1*procpaqsz);
380 complex<r_4>* zp2 = (complex<r_4>*)(procbuff1+i1*procpaqsz+procpaqsz/2);
381 complex<r_4>* zp3 = (complex<r_4>*)(procbuff2+i2*procpaqsz);
382 complex<r_4>* zp4 = (complex<r_4>*)(procbuff2+i2*procpaqsz+procpaqsz/2);
383 for(sa_size_t j=0; j<visiV13.Size(); j++) {
384 visiV13(j)+=zp1[j]*conj(zp3[j]);
385 visiV14(j)+=zp1[j]*conj(zp4[j]);
386 visiV23(j)+=zp2[j]*conj(zp3[j]);
387 visiV24(j)+=zp2[j]*conj(zp4[j]);
388 }
389 nzm++; i1++; i2++;
390 totnbytesproc += 2*paq1.DataSize();
391 } // Fin de boucle sur les paquets d'une zone
392 memgr1.FreeMemZone(mid1, MemZS_ProcB);
393 memgr2.FreeMemZone(mid2, MemZS_ProcB);
394
395 if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
396 visiV13 /= complex<r_4>((r_4)nzm, 0.);
397 visiV14 /= complex<r_4>((r_4)nzm, 0.);
398 visiV23 /= complex<r_4>((r_4)nzm, 0.);
399 visiV24 /= complex<r_4>((r_4)nzm, 0.);
400 visiV13.Info()["NPaqMoy"] = nzm;
401 visiV14.Info()["NPaqMoy"] = nzm;
402 visiV23.Info()["NPaqMoy"] = nzm;
403 visiV24.Info()["NPaqMoy"] = nzm;
404 char fname[512];
405 {
406 sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
407 POutPersist po(fname);
408 po << PPFNameTag("visiV13") << visiV13;
409 po << PPFNameTag("visiV14") << visiV14;
410 po << PPFNameTag("visiV23") << visiV23;
411 po << PPFNameTag("visiV24") << visiV24;
412 }
413 visiV13 = complex<r_4>(0., 0.);
414 visiV14 = complex<r_4>(0., 0.);
415 visiV23 = complex<r_4>(0., 0.);
416 visiV24 = complex<r_4>(0., 0.);
417 nzm = 0; ifile++;
418// ts.SetNow();
419// filog << ts << " : proc file " << fname << endl;
420 cout << " BRProcBRaw4C::run() created file " << fname << endl;
421 }
422 double okfrac = (nokpaq>1)?((double)noksfc/(double)nokpaq*100.):0.;
423 cout << "BRProcBRaw4C ["<<kmz<<"] NOKPaq=" << nokpaq << " NSameFC=" << noksfc
424 << " (" << okfrac << " %)" << endl;
425 totnokpaq += nokpaq;
426 totnoksfc += noksfc;
427 } // Fin de boucle sur les zones a traiter
428 cout << " ------------------ BRProcBRaw4C::run() END ----------------- " << endl;
429 ts.SetNow();
430 tm.SplitQ();
431 cout << " TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= "
432 << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s" << endl;
433 double totokfrac = (totnokpaq>1)?((double)totnoksfc/(double)totnokpaq*100.):0.;
434 cout << " NOkPaq1,2=" << totnokpaq << " /TotNPaq=" << totnpaq << " TotNSameFC="
435 << totnoksfc << " (" << totokfrac << " %)" << endl;
436// cout << pcheck1;
437// cout << pcheck2;
438 cout << " BRProcBRaw4C::run()/Timing: \n";
439 tm.Print();
440 cout << " ---------------------------------------------------------- " << endl;
441}
442 catch (PException& exc) {
443 cout << " BRProcBRaw4C::run()/catched PException " << exc.Msg() << endl;
444 setRC(3);
445 return;
446 }
447 catch(...) {
448 cout << " BRProcBRaw4C::run()/catched unknown ... exception " << endl;
449 setRC(4);
450 return;
451 }
452 setRC(0);
453 return;
454}
455
456
Note: See TracBrowser for help on using the repository browser.