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

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

Correction bugs ProcD/E/F ds memmgr , Reza 18/05/2010

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