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

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

Suite devel. BRVisibilityCalculator::Process() - Reza 28/11/2009
File size: 27.3 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 "datatable.h"
13#include "histos.h"
14#include "fioarr.h"
15#include "matharr.h"
16#include "timestamp.h"
17#include "ctimer.h"
18#include "fftpserver.h"
19#include "fftwserver.h"
20
21#include "FFTW/fftw3.h"
22
23
24#include "pciewrap.h"
25#include "brpaqu.h"
26#include "brproc.h"
27
28//---------------------------------------------------------------------
29// Classe de traitement - calcul de visibilite pour n fibres
30//---------------------------------------------------------------------
31
32/* --Methode-- */
33BRVisibilityCalculator::BRVisibilityCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean,
34 uint_4 freq1, uint_4 freq2, uint_4 nbfreq)
35 : BRBaseProcessor(memgr), outpath_(outpath), nmean_(nmean), numfreq1_(freq1), numfreq2_(freq2), nbinfreq_(nbfreq)
36{
37 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
38 vframecount_.push_back(0);
39 vtimetag_.push_back(0);
40 vnpaqm_.push_back(0);
41 }
42
43 BRPaquet paq(memgr_.PaqSize());
44
45 for(size_t i=0; i<2*memgr_.NbFibres(); i++)
46 for(size_t j=i; j<2*memgr_.NbFibres(); j++)
47 vvis_.push_back( TVector< complex<r_4> >(paq.DataSize()/4) );
48
49 cout << " BRVisibilityCalculator/Info NMean= " << nmean_ << endl;
50 totnbpaq_ = 0;
51}
52
53/* --Methode-- */
54BRVisibilityCalculator::~BRVisibilityCalculator()
55{
56}
57
58/* --Methode-- */
59int BRVisibilityCalculator::Process()
60{
61 /*
62 if (totnbpaq_%nmean_ == 0) {
63 cout << " BRVisibilityCalculator::Process() " << totnbpaq_ << " FrameCnt=" ;
64 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)
65 cout << curfc_[fib] << "," ;
66 cout << endl;
67 cout << " TimeTag (sec) : " ;
68 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)
69 cout << vpaq_[fib].TimeTag()/125000000 << "," ;
70 cout << endl;
71 }
72 */
73 vector<TwoByteComplex*> vpdata;
74 vector< TVector< complex<r_4> > > vvis;
75
76 /*
77 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
78 vpdata.push_back(vpaq_[fib].Data1C());
79 vpdata.push_back(vpaq_[fib].Data2C());
80 }
81 */
82 // On s'arrange pour mettre les canaux dans l'ordre pour les observations de Pittsburgh Nov 2009
83 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)
84 vpdata.push_back(vpaq_[fib].Data1C());
85 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++)
86 vpdata.push_back(vpaq_[fib].Data2C());
87
88 if (totnbpaq_%nmean_ == 0)
89 for(size_t k=0; k<vvis_.size(); k++) vvis_[k] = complex<r_4>((r_4)0.,(r_4)0.);
90
91 for(size_t i=0; i<vvis_.size(); i++) {
92 for(sa_size_t f=1; f<vvis_[i].Size(); f++) {
93 vvis_[i](f) += complex<r_4>((r_4)vpdata[0][f].realB(), (r_4)(r_4)vpdata[0][f].imagB());
94 }
95 }
96
97 totnbpaq_++;
98 return 0;
99}
100
101//---------------------------------------------------------------
102// Classe thread de traitement donnees ADC avec 2 voies par frame
103//---------------------------------------------------------------
104
105// Mutex pour eviter le plantage du a FFTW qui ne semble pas thread-safe
106static ZMutex* pmutfftw=NULL;
107
108/* --Methode-- */
109BRProcA2C::BRProcA2C(RAcqMemZoneMgr& mem, string& path, bool fgraw, uint_4 nmean,
110 uint_4 nmax, bool fghist, uint_4 nfsmap, bool fgnotrl, int card)
111 : memgr(mem)
112{
113 fgraw_ = fgraw;
114 nmax_ = nmax;
115 nmean_ = nmean;
116 if (fgraw_) cout << " BRProcA2C::BRProcA2C() - constructeur RAW data - NMean=" << nmean_ << endl;
117 else cout << " BRProcA2C::BRProcA2C() - constructeur FFT data - NMean=" << nmean_ << endl;
118 nfsmap_ = nfsmap;
119 stop_ = false;
120 path_ = path;
121 fgnotrl_ = fgnotrl;
122 fghist_ = fghist;
123 card_ = card;
124 if (pmutfftw==NULL) pmutfftw=new ZMutex;
125}
126
127/* --Methode-- */
128void BRProcA2C::Stop()
129{
130 stop_=true;
131 // cout <<" BRProcA2C::Stop ... > STOP " << endl;
132}
133
134
135static inline r_4 Zmod2(complex<r_4> z)
136{ return (z.real()*z.real()+z.imag()*z.imag()); }
137
138static inline string card2name_(int card)
139{
140 if (card==2) return " (Chan3,4) ";
141 else return " (Chan1,2) ";
142}
143/* --Methode-- */
144void BRProcA2C::run()
145{
146 setRC(1);
147 try {
148 Timer tm("BRProcA2C", false);
149 TimeStamp ts;
150 BRPaqChecker pcheck(!fgnotrl_); // Verification/comptage des paquets
151
152 size_t totnbytesout = 0;
153 size_t totnbytesproc = 0;
154
155 cout << " BRProcA2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
156 << " NMean=" << nmean_ << card2name_(card_) << endl;
157 cout << " BRProcA2C::run()... - Output Data Path: " << path_ << endl;
158 char fname[512];
159// sprintf(fname,"%s/proc.log",path_.c_str());
160// ofstream filog(fname);
161// filog << " BRProcA2C::run() - starting log file " << ts << endl;
162// filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl;
163
164/*----DELETE NTuple
165 const char* nnames[8] = {"fcs","tts","s1","s2","s12","s12re","s12im","s12phi"};
166 NTuple nt(8, nnames);
167 double xnt[10];
168 uint_4 nmnt = 0;
169 double ms1,ms2,ms12,ms12re,ms12im,ms12phi;
170----*/
171// Time sample (raw data) /FFT coeff histograms
172 Histo* ph1=NULL;
173 Histo* ph2=NULL;
174 if (fghist_) {
175 if (fgraw_) {
176 ph1 = new Histo(-0.5, 255.5, 256);
177 ph2 = new Histo(-0.5, 255.5, 256);
178 }
179 else {
180 ph1 = new Histo(-128.5, 128.5, 257);
181 ph2 = new Histo(-128.5, 128.5, 257);
182 }
183 }
184
185// Initialisation pour calcul FFT
186 TVector< complex<r_4> > cfour1; // composant TF
187 uint_4 paqsz = memgr.PaqSize();
188 uint_4 procpaqsz = memgr.ProcPaqSize();
189
190
191 BRPaquet pq(NULL, NULL, paqsz);
192 TVector<r_4> vx(pq.DataSize()/2);
193 int szfour = pq.DataSize()/2/2+1;
194 cfour1.SetSize(szfour);
195/*
196 vx = (r_4)(0.);
197 FFTPackServer ffts;
198 ffts.FFTForward(vx, cfour1);
199 szfour = cfour1.Size();
200*/
201
202 bool fgtimfreq = false; // true->cartes temps<>frequences
203 if (nfsmap_>0) fgtimfreq=true;
204
205 TVector< complex<r_4> > cfour2(cfour1.Size());
206
207 TVector<r_4> spectreV1(cfour1.Size());
208 TVector<r_4> spectreV2(cfour1.Size());
209 TVector<r_4> moyspecV1(cfour1.Size()); // Moyenne des Spectres
210 TVector<r_4> moyspecV2(cfour1.Size());
211 TVector<r_4> sigspecV1(cfour1.Size()); // Sigma des Spectres
212 TVector<r_4> sigspecV2(cfour1.Size());
213 TVector< complex<r_4> > visiV12( cfour1.Size() );
214
215 TMatrix<r_4> timfreqV1, timfreqV2; // Cartes temps<>frequences
216 if (fgtimfreq) {
217 timfreqV1.SetSize(nmean_, spectreV1.Size()/nfsmap_);
218 timfreqV2.SetSize(nmean_, spectreV2.Size()/nfsmap_);
219 }
220 cout << " *DBG*BRProcA2C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
221 << " procpaqsz/2=" << procpaqsz/2 << " cfour1.Size()=" << cfour1.Size()
222 << " *8=" << cfour1.Size()*8 << endl;
223
224 pmutfftw->lock();
225 fftwf_plan plan1 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(),
226 (fftwf_complex*)cfour1.Data(), FFTW_ESTIMATE);
227 fftwf_plan plan2 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(),
228 (fftwf_complex*)cfour2.Data(), FFTW_ESTIMATE);
229 pmutfftw->unlock();
230
231 uint_4 ifile = 0;
232 uint_4 nzm = 0; // Nb de paquets moyennes pour le calcul de chaque spectre
233 uint_4 nmoyspec = 0; // Nb de spectres moyennes
234
235 uint_4 curfc=0;
236 uint_8 curtt=0;
237 uint_8 firsttt=0;
238 bool fgfirst=true;
239 double moysig[2]={0.,0.};
240 double sigsig[2]={0.,0.};
241 uint_8 nbsig[2]={0,0};
242
243 for (uint_4 kmz=0; kmz<nmax_; kmz++) {
244 if (stop_) break;
245 int mid = memgr.FindMemZoneId(MemZA_ProcA);
246 Byte* buff = memgr.GetMemZone(mid);
247 if (buff == NULL) {
248 cout << " BRProcA2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
249 break;
250 }
251 Byte* procbuff = memgr.GetProcMemZone(mid);
252 if (procbuff == NULL) {
253 cout << " BRProcA2C::run()/ERROR memgr.GetProcMemZone(" << mid << ") -> NULL" << endl;
254 break;
255 }
256//---- DELETE nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=0.;
257 for(uint_4 i=0; i<memgr.NbPaquets(); i++) {
258 BRPaquet paq(NULL, buff+i*paqsz, paqsz);
259 if (!pcheck.Check(paq)) continue; // on ne traite que les paquets OK
260 if (fgfirst) { firsttt=paq.TimeTag(); fgfirst=false; }
261 curfc=paq.FrameCounter();
262 curtt=paq.TimeTag()-firsttt;
263// Traitement voie 1
264 if (fghist_) {
265 for(sa_size_t j=0; j<vx.Size(); j++) {
266 r_4 vts=(fgraw_)?((r_4)(*(paq.Data1()+j))):((r_4)(*(paq.Data1S()+j)));
267 ph1->Add((r_8)vts);
268 moysig[0] += (double)vts;
269 sigsig[0] += ((double)vts)*((double)vts);
270 nbsig[0]++;
271 }
272 for(sa_size_t j=0; j<vx.Size(); j++) {
273 r_4 vts=(fgraw_)?((r_4)(*(paq.Data2()+j))):((r_4)(*(paq.Data2S()+j)));
274 ph2->Add((r_8)vts);
275 moysig[1] += (double)vts;
276 sigsig[1] += ((double)vts)*((double)vts);
277 nbsig[1]++;
278 }
279 }
280 if (fgraw_) {
281 for(sa_size_t j=0; j<vx.Size(); j++)
282 vx(j) = (r_4)(*(paq.Data1()+j))-127.5;
283 // fftwf_complex* coeff1 = (fftwf_complex*)(procbuff+i*procpaqsz);
284 fftwf_execute(plan1);
285 // Traitement voie 2
286 for(sa_size_t j=0; j<vx.Size(); j++)
287 vx(j) = (r_4)(*(paq.Data2()+j))-127.5;
288 fftwf_execute(plan2);
289 }
290 else {
291 for(sa_size_t j=1; j<cfour1.Size()-1; j++) {
292 cfour1(j) = complex<r_4>((r_4)paq.Data1C()[j].realB(), (r_4)paq.Data1C()[j].imagB());
293 cfour2(j) = complex<r_4>((r_4)paq.Data2C()[j].realB(), (r_4)paq.Data2C()[j].imagB());
294 }
295 cfour1(0) = complex<r_4>((r_4)paq.Data1C()[0].realB(), (r_4)0.);
296 cfour1(cfour1.Size()-1) = complex<r_4>((r_4)paq.Data1C()[0].imagB(), (r_4)0.);
297 cfour2(0) = complex<r_4>((r_4)paq.Data2C()[0].realB(), (r_4)0.);
298 cfour2(cfour2.Size()-1) = complex<r_4>((r_4)paq.Data2C()[0].imagB(), (r_4)0.);
299 }
300 for(sa_size_t j=0; j<spectreV1.Size(); j++)
301 spectreV1(j) += Zmod2(cfour1(j));
302 memcpy(procbuff+i*procpaqsz, cfour1.Data(), sizeof(complex<r_4>)*cfour1.Size());
303 if (fgtimfreq) { // Remplissage tableau temps-frequence
304 for(sa_size_t c=1; c<timfreqV1.NCols(); c++) {
305 for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++)
306 timfreqV1(nzm, c) += Zmod2(cfour1(j));
307 }
308 }
309 for(sa_size_t j=0; j<spectreV2.Size(); j++)
310 spectreV2(j) += Zmod2(cfour2(j)); // Zmod2(zp2[j]);
311 memcpy(procbuff+i*procpaqsz+procpaqsz/2, cfour2.Data(), sizeof(complex<r_4>)*cfour2.Size());
312 if (fgtimfreq) { // Remplissage tableau temps-frequence
313 for(sa_size_t c=1; c<timfreqV2.NCols(); c++) {
314 for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++)
315 timfreqV2(nzm,c) += Zmod2(cfour2(j));
316 }
317 }
318
319// Calcul correlation (visibilite V1 * V2)
320 for(sa_size_t j=0; j<visiV12.Size(); j++)
321 visiV12(j)+=cfour1(j)*conj(cfour2(j));
322// for(sa_size_t j=0; j<visiV12.Size(); j++) visiV12(j)+=zp1[j]*zp2[j];
323 if (nzm==0) {
324 spectreV1.Info()["StartFC"] = curfc;
325 spectreV2.Info()["StartFC"] = curfc;
326 visiV12.Info()["StartFC"] = curfc;
327 spectreV1.Info()["StartTT"] = curtt;
328 spectreV2.Info()["StartTT"] = curtt;
329 visiV12.Info()["StartTT"] = curtt;
330 }
331 nzm++;
332/*----DELETE
333 if (nmnt==0) { xnt[0]=paq.FrameCounter(); xnt[1]=paq.TimeTag(); }
334 for(sa_size_t j=2700; j<2800; j++) {
335 ms1 += Zmod2(cfour1(j)); ms2 += Zmod2(cfour2(j));
336 complex<r_4> zvis = cfour1(j)*conj(cfour2(j));
337 ms12 += Zmod2(zvis); ms12re += zvis.real(); ms12im += zvis.imag();
338 ms12phi+= atan2(zvis.imag(),zvis.real());
339 }
340 nmnt++;
341----*/
342 totnbytesproc += paq.DataSize();
343 totnbytesout += (2*sizeof(complex<r_4>)*cfour1.Size());
344
345 } // Fin de boucle sur les paquets d'une zone
346
347/*---- DELETE
348 if (nmnt>0) {
349 double fnorm = (double)nmnt*(2800-2700);
350 xnt[2] = ms1 /= fnorm;
351 xnt[3] = ms2 /= fnorm;
352 xnt[4] = ms12 /= fnorm;
353 xnt[5] = ms12re /= fnorm;
354 xnt[6] = ms12im /= fnorm;
355 xnt[7] = ms12phi /= fnorm;
356 nt.Fill(xnt);
357 }
358----*/
359 if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
360 spectreV1 /= (r_4)(nzm);
361 spectreV2 /= (r_4)(nzm);
362
363 // pour le calcul des moyennes et sigmas de ces spectres
364 moyspecV1 += spectreV1;
365 moyspecV2 += spectreV2;
366 sigspecV1 += (spectreV1 && spectreV1);
367 sigspecV2 += (spectreV2 && spectreV2);
368 nmoyspec++;
369
370 visiV12 /= complex<r_4>((r_4)nzm, 0.);
371
372 spectreV1.Info()["NPaqMoy"] = nzm;
373 spectreV2.Info()["NPaqMoy"] = nzm;
374 visiV12.Info()["NPaqMoy"] = nzm;
375 spectreV1.Info()["EndFC"] = curfc;
376 spectreV2.Info()["EndFC"] = curfc;
377 visiV12.Info()["EndFC"] = curfc;
378 spectreV1.Info()["EndTT"] = curtt;
379 spectreV2.Info()["EndTT"] = curtt;
380 visiV12.Info()["EndTT"] = curtt;
381 {
382 sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
383 POutPersist po(fname);
384 string tag1="specV1";
385 string tag2="specV2";
386 string tag12="visiV12";
387 string tagh1="tshV1";
388 string tagh2="tshV2";
389 string tagtf1="timfreqV1";
390 string tagtf2="timfreqV2";
391 if (card_==2) {
392 tag1 = "specV3";
393 tag2 = "specV4";
394 tagh1 = "tshV1";
395 tagh2 = "tshV2";
396 tag12="visiV34";
397 tagtf1="timfreqV3";
398 tagtf2="timfreqV4";
399 }
400 po << PPFNameTag(tag1) << spectreV1;
401 po << PPFNameTag(tag2) << spectreV2;
402 po << PPFNameTag(tag12) << visiV12;
403 if (fghist_) {
404 po << PPFNameTag(tagh1) << (*ph1);
405 po << PPFNameTag(tagh2) << (*ph2);
406
407 double sspvmax[3] = {0.,0.,0.};
408 int_4 sspvmaxidx[3] = {-1,-1,-1};
409 for(int jji=1;jji<visiV12.Size()-1;jji++) {
410 r_4 zmv2 = Zmod2(visiV12(jji));
411 if (zmv2>sspvmax[2]) { sspvmax[2]=zmv2; sspvmaxidx[2]=jji; }
412 }
413 TVector<r_4>& sspv = spectreV1;
414 for(int ic=0; ic<2; ic++) {
415 if (ic==1) sspv = spectreV2;
416 for(int jji=1;jji<sspv.Size()-1;jji++)
417 if (sspv(jji)>sspvmax[ic]) { sspvmax[ic]=sspv(jji); sspvmaxidx[ic]=jji; }
418 if (nbsig[ic] < 1) { moysig[ic]=sigsig[ic]=-1.; }
419 else {
420 moysig[ic] /= (double)nbsig[ic];
421 sigsig[ic] /= (double)nbsig[ic];
422 sigsig[ic] -= (moysig[ic]*moysig[ic]);
423 sigsig[ic] = sqrt(sigsig[ic]);
424 cout << "===Voie " << ic << " Moy=" << moysig[ic] << " Sig=" << sigsig[ic]
425 << " MaxSpec Amp= " << sqrt(sspvmax[ic])/double(pq.DataSize()/2/2)
426 << " Pos=" << sspvmaxidx[ic] << " (NPts=" << nbsig[ic] << ")" << endl;
427 }
428 }
429 cout << "=== Voie1x2 MaxSpec Amp= " << sqrt(sqrt(sspvmax[2])/double(pq.DataSize()/2/2))
430 << " Pos=" << sspvmaxidx[2] << endl;
431 } // fin if (fghist_)
432
433 if (fgtimfreq) {
434 timfreqV1 /= (r_4)nzm;
435 timfreqV2 /= (r_4)nzm;
436 po << PPFNameTag(tagtf1) << timfreqV1;
437 po << PPFNameTag(tagtf2) << timfreqV2;
438 }
439 }
440 spectreV1 = (r_4)(0.);
441 spectreV2 = (r_4)(0.);
442 visiV12 = complex<r_4>(0., 0.);
443 if (fghist_) {
444 ph1->Zero();
445 ph2->Zero();
446 moysig[0]=moysig[1]=0.;
447 sigsig[0]=sigsig[1]=0.;
448 nbsig[0]=nbsig[1]=0;
449 }
450 if (fgtimfreq) {
451 timfreqV1 = (r_4)(0.);
452 timfreqV2 = (r_4)(0.);
453 }
454 nzm = 0; ifile++;
455// ts.SetNow();
456// filog << ts << " : proc file " << fname << endl;
457 cout << " BRProcA2C::run() created file " << fname << card2name_(card_) << endl;
458 }
459
460 memgr.FreeMemZone(mid, MemZS_ProcA);
461 } // Fin de boucle sur les zones a traiter
462 cout << " ------------ BRProcA2C::run() END " << card2name_(card_)
463 << " ------------ " << endl;
464/*---- DELETE
465 {
466 nt.Info()["FirstTT"]=firsttt;
467 cout << nt;
468 sprintf(fname,"%s_nt.ppf",path_.c_str());
469 POutPersist po(fname);
470 po << PPFNameTag("ntv12") << nt;
471 cout << " BRProcA2C::run() created NTuple file " << fname << card2name_(card_) << endl;
472 }
473---- */
474 if (nmoyspec>0) { // Calcul des moyennes et sigmas des spectres
475 r_4 fnms = nmoyspec;
476 moyspecV1 /= fnms;
477 moyspecV2 /= fnms;
478 sigspecV1 /= fnms;
479 sigspecV2 /= fnms;
480 sigspecV1 -= (moyspecV1 && moyspecV1);
481 sigspecV2 -= (moyspecV2 && moyspecV2);
482 sigspecV1 = Sqrt(sigspecV1);
483 sigspecV2 = Sqrt(sigspecV2);
484 TVector<r_4> rsbV1, rsbV2; // Rapport signal/bruit
485 moyspecV1.DivElt(sigspecV1, rsbV1, false, true);
486 moyspecV2.DivElt(sigspecV2, rsbV2, false, true);
487 sprintf(fname,"%s_ms.ppf",path_.c_str());
488 POutPersist po(fname);
489 po << PPFNameTag("moyspecV1") << moyspecV1;
490 po << PPFNameTag("moyspecV2") << moyspecV2;
491 po << PPFNameTag("sigspecV1") << sigspecV1;
492 po << PPFNameTag("sigspecV2") << sigspecV2;
493 po << PPFNameTag("rsbV1") << rsbV1;
494 po << PPFNameTag("rsbV2") << rsbV2;
495 cout << " BRProcA2C::run() created moysigspec file " << fname << card2name_(card_) << endl;
496 }
497
498 if (fghist_) {
499 delete ph1;
500 delete ph2;
501 }
502 ts.SetNow();
503 tm.SplitQ();
504 cout << " TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= "
505 << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s"
506 << " ProcDataOut=" << totnbytesout/(1024*1024) << " MB" << endl;
507 cout << pcheck;
508 cout << " BRProcA2C::run()/Timing: " << card2name_(card_) << endl;
509 tm.Print();
510 cout << " ---------------------------------------------------------- " << endl;
511
512 }
513 catch (PException& exc) {
514 cout << " BRProcA2C::run()/catched PException " << exc.Msg() << endl;
515 setRC(3);
516 return;
517 }
518 catch(...) {
519 cout << " BRProcA2C::run()/catched unknown ... exception " << endl;
520 setRC(4);
521 return;
522 }
523 setRC(0);
524 return;
525}
526
527//---------------------------------------------------------------------
528// Classe thread de traitement 2 x 2 voies/frames (Apres BRProcA2C)
529//---------------------------------------------------------------------
530
531/* --Methode-- */
532BRProcB4C::BRProcB4C(RAcqMemZoneMgr& mem1, RAcqMemZoneMgr& mem2, string& path,
533 bool fgraw, uint_4 nmean, uint_4 nmax, bool fgnotrl)
534 : memgr1(mem1), memgr2(mem2)
535{
536 fgraw_ = fgraw;
537 nmax_ = nmax;
538 nmean_ = nmean;
539 if (fgraw_) cout << " BRProcB4C::BRProcB4C() - constructeur RAW data - NMean= " << nmean_ << endl;
540 else cout << " BRProcB4C::BRProcB4C() - constructeur FFT data - NMean= " << nmean_ << endl;
541 stop_ = false;
542 path_ = path;
543 fgnotrl_ = fgnotrl;
544}
545
546/* --Methode-- */
547void BRProcB4C::Stop()
548{
549 stop_=true;
550 // cout <<" BRProcB4C::Stop ... > STOP " << endl;
551}
552
553
554/* --Methode-- */
555void BRProcB4C::run()
556{
557 setRC(1);
558 try {
559 Timer tm("BRProcB4C", false);
560 TimeStamp ts;
561 BRPaqChecker pcheck1(!fgnotrl_); // Verification/comptage des paquets
562 BRPaqChecker pcheck2(!fgnotrl_); // Verification/comptage des paquets
563
564 size_t totnbytesout = 0;
565 size_t totnbytesproc = 0;
566
567 cout << " BRProcB4C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
568 << " NMean=" << nmean_ << endl;
569 cout << " BRProcB4C::run()... - Output Data Path: " << path_ << endl;
570
571 uint_4 paqsz = memgr1.PaqSize();
572 uint_4 procpaqsz = memgr1.ProcPaqSize();
573 if ((paqsz != memgr2.PaqSize())||(procpaqsz!= memgr2.ProcPaqSize())) {
574 cout << "BRProcB4C::run()/ERROR : different paquet size -> stop \n ...(PaqSz1="
575 << paqsz << " Sz2=" << memgr2.PaqSize() << " ProcPaqSz1="
576 << procpaqsz << " Sz2=" << memgr2.ProcPaqSize() << " )" << endl;
577 setRC(9);
578 return;
579 }
580
581 TVector< complex<r_4> > cfour; // composant TF
582 BRPaquet pq(NULL, NULL, paqsz);
583/*
584 TVector<r_4> vx(pq.DataSize()/2);
585 vx = (r_4)(0.);
586 FFTPackServer ffts;
587 ffts.FFTForward(vx, cfour);
588
589 TVector< complex<r_4> > visiV13( cfour.Size() );
590 TVector< complex<r_4> > visiV14( cfour.Size() );
591 TVector< complex<r_4> > visiV23( cfour.Size() );
592 TVector< complex<r_4> > visiV24( cfour.Size() );
593*/
594 int szfour = pq.DataSize()/2/2+1;
595// int szfour = (paqsz-40)/2+1;
596 TVector< complex<r_4> > visiV13( szfour );
597 TVector< complex<r_4> > visiV14( szfour );
598 TVector< complex<r_4> > visiV23( szfour );
599 TVector< complex<r_4> > visiV24( szfour );
600 // cout << " *DBG*AAAAA ---- Vectors OK" << endl;
601 cout << " *DBG*BRProcB4C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
602 << " procpaqsz/2=" << procpaqsz/2 << " cfour.Size()=" << szfour
603 << " *8=" << szfour*8 << endl;
604
605 DataTable dt;
606 dt.AddLongColumn("fc1");
607 dt.AddLongColumn("tt1");
608 dt.AddLongColumn("fc2");
609 dt.AddLongColumn("tt2");
610 DataTableRow dtr = dt.EmptyRow();
611
612 uint_4 nzm = 0;
613 uint_4 totnoksfc = 0;
614 uint_4 totnokpaq = 0;
615 uint_4 totnpaq = 0;
616 uint_4 ifile = 0;
617
618 uint_4 curfc=0;
619 uint_8 curtt=0;
620 uint_4 curfc2=0;
621 uint_8 curtt2=0;
622 uint_8 firsttt=0;
623 uint_8 firsttt2=0;
624 bool fgfirst=true;
625 for (uint_4 kmz=0; kmz<nmax_; kmz++) {
626 uint_4 noksfc = 0;
627 uint_4 nokpaq = 0;
628 if (stop_) break;
629 // cout << " *DBG*BBBBB" << kmz << endl;
630
631 int mid1 = memgr1.FindMemZoneId(MemZA_ProcB);
632 Byte* buff1 = memgr1.GetMemZone(mid1);
633 if (buff1 == NULL) {
634 cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid1 << ") -> NULL" << endl;
635 break;
636 }
637 Byte* procbuff1 = memgr1.GetProcMemZone(mid1);
638 if (procbuff1 == NULL) {
639 cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid1 << ") -> NULL" << endl;
640 break;
641 }
642 int mid2 = memgr2.FindMemZoneId(MemZA_ProcB);
643 Byte* buff2 = memgr2.GetMemZone(mid2);
644 if (buff1 == NULL) {
645 cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid2 << ") -> NULL" << endl;
646 break;
647 }
648 Byte* procbuff2 = memgr2.GetProcMemZone(mid2);
649 if (procbuff2 == NULL) {
650 cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid2 << ") -> NULL" << endl;
651 break;
652 }
653 uint_4 i1,i2;
654 i1=i2=0;
655// cout << " *DBG*CCCCCC " << kmz << " memgr1.NbPaquets() =" << memgr1.NbPaquets() << endl;
656 while((i1<memgr1.NbPaquets())&&(i2<memgr2.NbPaquets())) {
657 BRPaquet paq1(NULL, buff1+i1*paqsz, paqsz);
658 BRPaquet paq2(NULL, buff2+i2*paqsz, paqsz);
659 totnpaq++;
660// cout << " DBG["<<kmz<<"] i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
661//<<","<<paq2.FrameCounter()<<endl;
662 // on ne traite que les paquets OK
663 if (!pcheck1.Check(paq1)) { i1++; continue; }
664 if (!pcheck2.Check(paq2)) { i2++; continue; }
665 nokpaq++;
666 if (paq1.FrameCounter()<paq2.FrameCounter()) { i1++; continue; }
667 if (paq2.FrameCounter()<paq1.FrameCounter()) { i2++; continue; }
668// cout << " DBG["<<kmz<<"]OKOK i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter()
669// <<","<<paq2.FrameCounter()<<endl;
670
671 if ((i1>=memgr1.NbPaquets())||(i2>=memgr1.NbPaquets())) {
672 cout << " *BUG*BUG i1=" << i1 << " i2=" << i2 << endl;
673 break;
674 }
675 // Les deux framecounters sont identiques ...
676 noksfc++;
677 curfc=paq1.FrameCounter();
678 curfc2=paq2.FrameCounter();
679 if (fgfirst) {
680 firsttt=paq1.TimeTag(); firsttt2=paq2.TimeTag();
681 cout << " BRProcB4C()/Info First FC="<<curfc<<" , "<<curfc2<<" -> TT="
682 << firsttt<<" , "<<firsttt2 <<endl;
683 fgfirst=false;
684 }
685 curtt=paq1.TimeTag()-firsttt;
686 curtt2=paq2.TimeTag()-firsttt2;
687 dtr[0]=curfc; dtr[1]=curtt;
688 dtr[2]=curfc2; dtr[3]=curtt2;
689 dt.AddRow(dtr);
690
691 complex<r_4>* zp1 = (complex<r_4>*)(procbuff1+i1*procpaqsz);
692 complex<r_4>* zp2 = (complex<r_4>*)(procbuff1+i1*procpaqsz+procpaqsz/2);
693 complex<r_4>* zp3 = (complex<r_4>*)(procbuff2+i2*procpaqsz);
694 complex<r_4>* zp4 = (complex<r_4>*)(procbuff2+i2*procpaqsz+procpaqsz/2);
695 for(sa_size_t j=0; j<visiV13.Size(); j++) {
696 visiV13(j)+=zp1[j]*conj(zp3[j]);
697 visiV14(j)+=zp1[j]*conj(zp4[j]);
698 visiV23(j)+=zp2[j]*conj(zp3[j]);
699 visiV24(j)+=zp2[j]*conj(zp4[j]);
700 }
701 if (nzm==0) {
702 visiV13.Info()["StartFC"] = curfc;
703 visiV14.Info()["StartFC"] = curfc;
704 visiV23.Info()["StartFC"] = curfc;
705 visiV24.Info()["StartFC"] = curfc;
706 visiV13.Info()["StartTT"] = curtt;
707 visiV14.Info()["StartTT"] = curtt;
708 visiV23.Info()["StartTT"] = curtt;
709 visiV24.Info()["StartTT"] = curtt;
710 }
711 nzm++; i1++; i2++;
712 totnbytesproc += 2*paq1.DataSize();
713 } // Fin de boucle sur les paquets d'une zone
714 memgr1.FreeMemZone(mid1, MemZS_ProcB);
715 memgr2.FreeMemZone(mid2, MemZS_ProcB);
716
717 if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
718 visiV13 /= complex<r_4>((r_4)nzm, 0.);
719 visiV14 /= complex<r_4>((r_4)nzm, 0.);
720 visiV23 /= complex<r_4>((r_4)nzm, 0.);
721 visiV24 /= complex<r_4>((r_4)nzm, 0.);
722 visiV13.Info()["NPaqMoy"] = nzm;
723 visiV14.Info()["NPaqMoy"] = nzm;
724 visiV23.Info()["NPaqMoy"] = nzm;
725 visiV24.Info()["NPaqMoy"] = nzm;
726 visiV13.Info()["EndFC"] = curfc;
727 visiV14.Info()["EndFC"] = curfc;
728 visiV23.Info()["EndFC"] = curfc;
729 visiV24.Info()["EndFC"] = curfc;
730 visiV13.Info()["EndTT"] = curtt;
731 visiV14.Info()["EndTT"] = curtt;
732 visiV23.Info()["EndTT"] = curtt;
733 visiV24.Info()["EndTT"] = curtt;
734 char fname[512];
735 {
736 sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
737 POutPersist po(fname);
738 po << PPFNameTag("visiV13") << visiV13;
739 po << PPFNameTag("visiV14") << visiV14;
740 po << PPFNameTag("visiV23") << visiV23;
741 po << PPFNameTag("visiV24") << visiV24;
742 }
743 visiV13 = complex<r_4>(0., 0.);
744 visiV14 = complex<r_4>(0., 0.);
745 visiV23 = complex<r_4>(0., 0.);
746 visiV24 = complex<r_4>(0., 0.);
747 nzm = 0; ifile++;
748// ts.SetNow();
749// filog << ts << " : proc file " << fname << endl;
750 cout << " BRProcB4C::run() created file " << fname << endl;
751 }
752 double okfrac = (nokpaq>1)?((double)noksfc/(double)nokpaq*100.):0.;
753 cout << "BRProcB4C ["<<kmz<<"] NOKPaq=" << nokpaq << " NSameFC=" << noksfc
754 << " (" << okfrac << " %)" << endl;
755 totnokpaq += nokpaq;
756 totnoksfc += noksfc;
757 } // Fin de boucle sur les zones a traiter
758 cout << " ------------------ BRProcB4C::run() END ----------------- " << endl;
759 {
760 dt.Info()["FirstTT1"]=firsttt;
761 dt.Info()["FirstTT2"]=firsttt2;
762 cout << dt;
763 char fname[512];
764 sprintf(fname,"%s_fctt.ppf",path_.c_str());
765 POutPersist po(fname);
766 po << PPFNameTag("ttfc") << dt;
767 cout << " BRProcB4C::run() created TimeTag/FrameCounter file " << fname << endl;
768 }
769 ts.SetNow();
770 tm.SplitQ();
771 cout << " TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= "
772 << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s" << endl;
773 double totokfrac = (totnokpaq>1)?((double)totnoksfc/(double)totnokpaq*100.):0.;
774 cout << " NOkPaq1,2=" << totnokpaq << " /TotNPaq=" << totnpaq << " TotNSameFC="
775 << totnoksfc << " (" << totokfrac << " %)" << endl;
776// cout << pcheck1;
777// cout << pcheck2;
778 cout << " BRProcB4C::run()/Timing: \n";
779 tm.Print();
780 cout << " ---------------------------------------------------------- " << endl;
781}
782 catch (PException& exc) {
783 cout << " BRProcB4C::run()/catched PException " << exc.Msg() << endl;
784 setRC(3);
785 return;
786 }
787 catch(...) {
788 cout << " BRProcB4C::run()/catched unknown ... exception " << endl;
789 setRC(4);
790 return;
791 }
792 setRC(0);
793 return;
794}
795
796
Note: See TracBrowser for help on using the repository browser.