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

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

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