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

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

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