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

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

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