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

Last change on this file since 3774 was 3774, checked in by ansari, 15 years ago

Ajout/modif classes (brproc.h,.cc) pour analyse de donnees BAORadio, Reza 17/05/2010

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