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

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

Passage a 6 proc/thread max ds BRVisibilityCalculator , Reza 18/05/2010

File size: 41.9 KB
RevLine 
[3635]1#include "racquproc.h"
2
3#include <stdlib.h>
[3642]4#include <string.h>
[3635]5#include <unistd.h>
6#include <fstream>
7#include <signal.h>
8
9#include "pexceptions.h"
10#include "tvector.h"
[3646]11#include "ntuple.h"
[3647]12#include "datatable.h"
[3652]13#include "histos.h"
[3635]14#include "fioarr.h"
[3655]15#include "matharr.h"
[3635]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
[3683]28//---------------------------------------------------------------------
29// Classe de traitement - calcul de visibilite pour n fibres
30//---------------------------------------------------------------------
[3635]31
[3683]32/* --Methode-- */
[3776]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())
[3692]36 // , dtfos_(outpath+"visdt.fits", Fits_Create), visdt_(dtfos_, 1024, true);
37{
[3776]38 DefineRank(1,0);
[3689]39
[3776]40 uint_4 maxnpairs = (2*memgr_.NbFibres()+1)*memgr_.NbFibres();
41 chanum_.SetSize(maxnpairs);
[3692]42 sa_size_t k=0;
[3694]43 for(size_t i=0; i<2*memgr_.NbFibres(); i++) vpdata_[i]=NULL;
[3692]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 }
[3776]49 SelectPairs();
50
[3692]51 // visdt_.AddFloatColumn("mfc");
[3696]52 visdt_.AddDoubleColumn("mfc");
53 visdt_.AddDoubleColumn("mtt");
[3692]54 visdt_.AddIntegerColumn("jfreq");
55 visdt_.AddIntegerColumn("numch");
56 visdt_.AddFloatColumn("vre");
57 visdt_.AddFloatColumn("vim");
[3689]58
[3692]59 if (nmean_ < 1) nmean_=memgr_.NbPaquets();
60 if (nmean_ < 1) nmean_=1;
61
[3776]62 cout << " BRVisibilityCalculator::/Info nmean=" << nmean_ << endl;
63
[3692]64 totnbpaq_=0;
65 numfile_=0;
[3777]66 nb_flop_=0.;
[3692]67 moyfc_=moytt_=0.;
68
[3726]69 fgallfibok=NULL;
[3705]70 fgcktt_=false;
[3778]71 setNameId("viscalc", 0);
[3683]72}
73
74/* --Methode-- */
75BRVisibilityCalculator::~BRVisibilityCalculator()
76{
[3692]77 cout << " BRVisibilityCalculator - Visibility Datatable : " << endl;
78 cout << visdt_;
[3776]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);
[3692]87 po << visdt_;
[3776]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_;
[3705]101 }
102 }
[3683]103}
104
105/* --Methode-- */
[3776]106void BRVisibilityCalculator::DefineRank(uint_4 nbc, uint_4 cid)
107{
[3780]108 if ((nbc>6)||(cid>=nbc))
109 throw ParmError("BRVisibilityCalculator::DefineRank() NbCalc > 6 !");
[3776]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);
[3780]119 MemZaction mmzas[6]={MemZA_ProcA,MemZA_ProcB,MemZA_ProcC,MemZA_ProcD,MemZA_ProcE,MemZA_ProcF};
[3776]120 SetMemZAction(mmzas[calcid_]);
[3778]121 setNameId("viscalc_grp", calcid_);
[3776]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-- */
[3724]176void BRVisibilityCalculator::run()
177{
[3776]178 cout << " BRVisibilityCalculator[" << calcid_ << "/" << nbcalc_
179 << "]::run() - Starting " << " NFibers=" << memgr_.NbFibres()
180 << " NChan=" << 2*memgr_.NbFibres() << " NPairs=" << nbpairs_ << " First:" << pairst_ << endl;
181
[3726]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 }
[3724]193 paralex_.SetParallelTask(*this);
194 cout << " BRVisibilityCalculator::run()/Info : starting ParallelExecutor with nThreads="
195 << paralex_.nThreads() << " ... " << endl;
196 paralex_.start();
[3726]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;
[3774]205 int mid = memgr_.FindMemZoneId(mmact_); // (MemZA_ProcA);
[3779]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;
[3726]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
[3778]217 cout << "BRVisibilityCalculator[" << calcid_ << "]::run()/ERROR memgr.GetMemZone(" << mid << "," << fib << ") -> NULL" << endl;
[3726]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
[3776]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_);
[3726]237 string flnm=outpath_+nfile;
238 POutPersist po(flnm);
239 po << vismtx_;
[3778]240 cout << numfile_ << "-BRVisibilityCalculator[" << calcid_ << "]::run() NPaqProc="
[3726]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]);
[3776]266 if ((fgcktt_)&&(calcid_==0)) CheckTimeTag();
[3726]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();
[3778]273 if (rcpex!=0) cout << " BRVisibilityCalculator[" << calcid_ << "]::run() / Error Rc[paralex_.execute()]=" << rcpex << endl;
[3726]274
[3774]275 memgr_.FreeMemZone(mid, mmsta_); // (MemZS_ProcA);
[3726]276 } // Fin de boucle sur les zones a traiter
277 //------------------------------------
[3778]278 cout << " --------- END BRVisibilityCalculator[" << calcid_ << "]::run() , TotNbProcPaq=" << totprocnpaq_ << endl;
[3726]279 /*
280 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) vpchk_[fib].Print();
281 cout << " ------------------------------------ " << endl;
282 */
283 delete[] fgallfibok;
[3724]284 }
[3726]285 catch (std::exception& exc) {
[3778]286 cout << " BRVisibilityCalculator[" << calcid_ << "]::run()/catched std::exception " << exc.what() << endl;
[3726]287 setRC(98);
288 return;
289 }
290 catch(...) {
[3778]291 cout << " BRVisibilityCalculator[" << calcid_ << "]::run()/catched unknown ... exception " << endl;
[3726]292 setRC(99);
293 return;
294 }
295
[3724]296}
297
298/* --Methode-- */
[3683]299int BRVisibilityCalculator::Process()
300{
[3696]301
[3689]302 for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) {
[3694]303 vpdata_[2*fib] = vpaq_[fib].Data1C();
304 vpdata_[2*fib+1] = vpaq_[fib].Data2C();
305 }
[3696]306
[3692]307 if (totnbpaq_%nmean_ == 0) {
308 if (totnbpaq_ > 0) {
[3701]309 moyfc_/=nmean_;
310 moytt_/=nmean_;
311 vismtx_.Info()["MeanFC"] = moyfc_;
312 vismtx_.Info()["MeanTT"] = moytt_;
313 vismtx_.Info()["NPAQSUM"] = nmean_;
314
[3692]315 // ATTENTION : Matrice visibilites non moyennee
[3776]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_);
[3692]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_++;
[3689]328 }
[3692]329 vismtx_ = complex<r_4>((r_4)0.,(r_4)0.);
330 moyfc_=moytt_=0.;
[3689]331 }
[3724]332
[3726]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++) {
[3776]336 size_t kpair=i*vpdata_.size()+j;
337 if (kpair<pairst_) continue;
338 if (kpair>=(pairst_+nbpairs_)) break;
[3726]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());
[3724]343 }
[3777]344 nb_flop_ += (7.*(r_8)vis.Size());
[3724]345 }
346 }
[3726]347
[3724]348 moyfc_ += curfc_[0];
349 moytt_ += (vpaq_[0].TimeTag()-ttfirst_[0]);
[3776]350 if ((fgcktt_)&&(calcid_==0)) CheckTimeTag();
[3724]351 totnbpaq_++;
352 return 0;
353}
354
355/* --Methode-- */
356int BRVisibilityCalculator::execute(int tid)
357{
[3726]358 vector<TwoByteComplex*> pvpdata(2*memgr_.NbFibres());
359 size_t paqsz=memgr_.PaqSize();
360 BRPaquet ppaq(paqsz);
361
[3724]362 sa_size_t fdelt = vismtx_.NCols()/nparthr_;
363 sa_size_t fdeb = tid*fdelt;
364 sa_size_t ffin = (tid+1)*fdelt;
[3726]365
[3724]366 if (fdeb<1) fdeb=1;
367 if ((ffin>vismtx_.NCols())||(tid==(nparthr_-1))) ffin=vismtx_.NCols();
[3726]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++) {
[3776]379 size_t kpair=i*vpdata_.size()+j;
380 if (kpair<pairst_) continue;
381 if (kpair>=(pairst_+nbpairs_)) break;
[3726]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 }
[3777]387 nb_flop_ += (7.*(r_8)(ffin-fdeb));
[3692]388 }
389 }
[3726]390
391 } // Fin de boucle sur les paquets
392
[3686]393 return 0;
[3683]394}
395
[3692]396/* --Methode-- */
397int BRVisibilityCalculator::FillVisibTable(double fcm, double ttm)
398{
[3776]399 double xnt[10];
400 xnt[0]=fcm; xnt[1]=ttm/1.25e8;
[3692]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++) {
[3776]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);
[3692]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.;
[3705]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++) {
[3696]421 moyreal+=vismtx_(rv,jjf).real();
422 moyimag+=vismtx_(rv,jjf).imag();
[3692]423 }
[3776]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);
[3692]429 }
430 }
431 }
432 return 0;
433}
434
[3705]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++) {
[3709]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];
[3708]451 /* if ( (ld < mindeltatt_) || (fd<-5) || (fd>5)) { vbadtt_[fib]++; vnsamett_[fib]++; }
[3705]452 else {
453 if (fd!=0) vnsamett_[fib]++;
454 }
[3708]455 */
456 if (ld < mindeltatt_) vbadtt_[fib]++;
[3709]457 else {
458 if (fd != 0) vnsamett_[fib]++;
[3710]459 if ((fd<-5)||(fd>5)) vndiff5tt_[fib]++;
[3709]460 }
[3705]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
[3776]470//-------------------------------------------------------------------------------
471// Classe Groupe (ensemble) de Calculateur de Visibilites, tourant en parallele
472//-------------------------------------------------------------------------------
[3774]473
[3776]474/* --Methode-- */
475BRVisCalcGroup::BRVisCalcGroup(size_t nbcalc, RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, size_t nthr)
[3777]476 : tm_(false)
[3776]477{
[3780]478 if ((nbcalc<1)||(nbcalc>6))
479 throw ParmError("BRVisCalcGroup::BRVisCalcGroup NbCalc > 6 !");
[3776]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();
[3777]505 tm_.SplitQ();
[3776]506}
507/* --Methode-- */
508void BRVisCalcGroup::join()
509{
[3777]510 r_8 totflop=0.;
511 for(size_t i=0; i<viscalcp_.size(); i++) {
[3776]512 viscalcp_[i]->join();
[3777]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;
[3776]523}
524
525
[3774]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
[3635]580//---------------------------------------------------------------
581// Classe thread de traitement donnees ADC avec 2 voies par frame
582//---------------------------------------------------------------
583
[3649]584// Mutex pour eviter le plantage du a FFTW qui ne semble pas thread-safe
585static ZMutex* pmutfftw=NULL;
586
[3645]587/* --Methode-- */
[3683]588BRProcA2C::BRProcA2C(RAcqMemZoneMgr& mem, string& path, bool fgraw, uint_4 nmean,
[3656]589 uint_4 nmax, bool fghist, uint_4 nfsmap, bool fgnotrl, int card)
[3635]590 : memgr(mem)
591{
[3683]592 fgraw_ = fgraw;
[3635]593 nmax_ = nmax;
594 nmean_ = nmean;
[3683]595 if (fgraw_) cout << " BRProcA2C::BRProcA2C() - constructeur RAW data - NMean=" << nmean_ << endl;
596 else cout << " BRProcA2C::BRProcA2C() - constructeur FFT data - NMean=" << nmean_ << endl;
[3656]597 nfsmap_ = nfsmap;
[3635]598 stop_ = false;
599 path_ = path;
[3640]600 fgnotrl_ = fgnotrl;
[3652]601 fghist_ = fghist;
[3645]602 card_ = card;
[3649]603 if (pmutfftw==NULL) pmutfftw=new ZMutex;
[3635]604}
605
[3645]606/* --Methode-- */
[3683]607void BRProcA2C::Stop()
[3635]608{
609 stop_=true;
[3683]610 // cout <<" BRProcA2C::Stop ... > STOP " << endl;
[3635]611}
612
613
614static inline r_4 Zmod2(complex<r_4> z)
615{ return (z.real()*z.real()+z.imag()*z.imag()); }
616
[3645]617static inline string card2name_(int card)
618{
619 if (card==2) return " (Chan3,4) ";
620 else return " (Chan1,2) ";
621}
622/* --Methode-- */
[3683]623void BRProcA2C::run()
[3635]624{
625 setRC(1);
626 try {
[3683]627 Timer tm("BRProcA2C", false);
[3635]628 TimeStamp ts;
[3646]629 BRPaqChecker pcheck(!fgnotrl_); // Verification/comptage des paquets
[3640]630
631 size_t totnbytesout = 0;
632 size_t totnbytesproc = 0;
633
[3683]634 cout << " BRProcA2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
[3645]635 << " NMean=" << nmean_ << card2name_(card_) << endl;
[3683]636 cout << " BRProcA2C::run()... - Output Data Path: " << path_ << endl;
[3635]637 char fname[512];
638// sprintf(fname,"%s/proc.log",path_.c_str());
639// ofstream filog(fname);
[3683]640// filog << " BRProcA2C::run() - starting log file " << ts << endl;
[3635]641// filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl;
642
[3647]643/*----DELETE NTuple
[3646]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;
[3647]649----*/
[3683]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
[3635]664// Initialisation pour calcul FFT
[3640]665 TVector< complex<r_4> > cfour1; // composant TF
[3635]666 uint_4 paqsz = memgr.PaqSize();
667 uint_4 procpaqsz = memgr.ProcPaqSize();
[3646]668
669
[3635]670 BRPaquet pq(NULL, NULL, paqsz);
671 TVector<r_4> vx(pq.DataSize()/2);
[3648]672 int szfour = pq.DataSize()/2/2+1;
673 cfour1.SetSize(szfour);
674/*
[3635]675 vx = (r_4)(0.);
676 FFTPackServer ffts;
[3640]677 ffts.FFTForward(vx, cfour1);
[3648]678 szfour = cfour1.Size();
679*/
680
[3656]681 bool fgtimfreq = false; // true->cartes temps<>frequences
682 if (nfsmap_>0) fgtimfreq=true;
683
[3640]684 TVector< complex<r_4> > cfour2(cfour1.Size());
[3635]685
[3640]686 TVector<r_4> spectreV1(cfour1.Size());
687 TVector<r_4> spectreV2(cfour1.Size());
[3655]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());
[3640]692 TVector< complex<r_4> > visiV12( cfour1.Size() );
[3635]693
[3656]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 }
[3683]699 cout << " *DBG*BRProcA2C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
[3646]700 << " procpaqsz/2=" << procpaqsz/2 << " cfour1.Size()=" << cfour1.Size()
701 << " *8=" << cfour1.Size()*8 << endl;
[3635]702
[3649]703 pmutfftw->lock();
[3640]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);
[3649]708 pmutfftw->unlock();
[3640]709
[3635]710 uint_4 ifile = 0;
[3655]711 uint_4 nzm = 0; // Nb de paquets moyennes pour le calcul de chaque spectre
712 uint_4 nmoyspec = 0; // Nb de spectres moyennes
[3647]713
714 uint_4 curfc=0;
715 uint_8 curtt=0;
716 uint_8 firsttt=0;
717 bool fgfirst=true;
[3658]718 double moysig[2]={0.,0.};
719 double sigsig[2]={0.,0.};
720 uint_8 nbsig[2]={0,0};
721
[3635]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) {
[3683]727 cout << " BRProcA2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl;
[3640]728 break;
[3635]729 }
730 Byte* procbuff = memgr.GetProcMemZone(mid);
731 if (procbuff == NULL) {
[3683]732 cout << " BRProcA2C::run()/ERROR memgr.GetProcMemZone(" << mid << ") -> NULL" << endl;
[3640]733 break;
[3635]734 }
[3647]735//---- DELETE nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=0.;
[3645]736 for(uint_4 i=0; i<memgr.NbPaquets(); i++) {
[3640]737 BRPaquet paq(NULL, buff+i*paqsz, paqsz);
738 if (!pcheck.Check(paq)) continue; // on ne traite que les paquets OK
[3647]739 if (fgfirst) { firsttt=paq.TimeTag(); fgfirst=false; }
740 curfc=paq.FrameCounter();
741 curtt=paq.TimeTag()-firsttt;
[3635]742// Traitement voie 1
[3652]743 if (fghist_) {
744 for(sa_size_t j=0; j<vx.Size(); j++) {
[3683]745 r_4 vts=(fgraw_)?((r_4)(*(paq.Data1()+j))):((r_4)(*(paq.Data1S()+j)));
746 ph1->Add((r_8)vts);
[3658]747 moysig[0] += (double)vts;
748 sigsig[0] += ((double)vts)*((double)vts);
749 nbsig[0]++;
[3652]750 }
[3683]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);
[3658]754 moysig[1] += (double)vts;
755 sigsig[1] += ((double)vts)*((double)vts);
756 nbsig[1]++;
[3652]757 }
758 }
[3683]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 }
[3635]788 for(sa_size_t j=0; j<spectreV2.Size(); j++)
[3640]789 spectreV2(j) += Zmod2(cfour2(j)); // Zmod2(zp2[j]);
790 memcpy(procbuff+i*procpaqsz+procpaqsz/2, cfour2.Data(), sizeof(complex<r_4>)*cfour2.Size());
[3656]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 }
[3635]797
798// Calcul correlation (visibilite V1 * V2)
[3640]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];
[3647]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 }
[3640]810 nzm++;
[3647]811/*----DELETE
[3646]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++;
[3647]820----*/
[3640]821 totnbytesproc += paq.DataSize();
822 totnbytesout += (2*sizeof(complex<r_4>)*cfour1.Size());
[3635]823
[3640]824 } // Fin de boucle sur les paquets d'une zone
[3647]825
826/*---- DELETE
[3646]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 }
[3647]837----*/
[3640]838 if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) {
[3635]839 spectreV1 /= (r_4)(nzm);
840 spectreV2 /= (r_4)(nzm);
841
[3655]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
[3635]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;
[3647]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;
[3652]860 {
[3635]861 sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile);
862 POutPersist po(fname);
[3645]863 string tag1="specV1";
864 string tag2="specV2";
865 string tag12="visiV12";
[3652]866 string tagh1="tshV1";
867 string tagh2="tshV2";
[3656]868 string tagtf1="timfreqV1";
869 string tagtf2="timfreqV2";
[3645]870 if (card_==2) {
871 tag1 = "specV3";
872 tag2 = "specV4";
[3652]873 tagh1 = "tshV1";
874 tagh2 = "tshV2";
[3645]875 tag12="visiV34";
[3656]876 tagtf1="timfreqV3";
877 tagtf2="timfreqV4";
[3645]878 }
879 po << PPFNameTag(tag1) << spectreV1;
880 po << PPFNameTag(tag2) << spectreV2;
881 po << PPFNameTag(tag12) << visiV12;
[3652]882 if (fghist_) {
[3683]883 po << PPFNameTag(tagh1) << (*ph1);
884 po << PPFNameTag(tagh2) << (*ph2);
[3658]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
[3656]912 if (fgtimfreq) {
913 timfreqV1 /= (r_4)nzm;
914 timfreqV2 /= (r_4)nzm;
915 po << PPFNameTag(tagtf1) << timfreqV1;
916 po << PPFNameTag(tagtf2) << timfreqV2;
917 }
[3652]918 }
[3635]919 spectreV1 = (r_4)(0.);
920 spectreV2 = (r_4)(0.);
921 visiV12 = complex<r_4>(0., 0.);
[3652]922 if (fghist_) {
[3683]923 ph1->Zero();
924 ph2->Zero();
[3658]925 moysig[0]=moysig[1]=0.;
926 sigsig[0]=sigsig[1]=0.;
927 nbsig[0]=nbsig[1]=0;
[3652]928 }
[3656]929 if (fgtimfreq) {
930 timfreqV1 = (r_4)(0.);
931 timfreqV2 = (r_4)(0.);
932 }
[3635]933 nzm = 0; ifile++;
934// ts.SetNow();
935// filog << ts << " : proc file " << fname << endl;
[3683]936 cout << " BRProcA2C::run() created file " << fname << card2name_(card_) << endl;
[3640]937 }
[3635]938
939 memgr.FreeMemZone(mid, MemZS_ProcA);
[3640]940 } // Fin de boucle sur les zones a traiter
[3683]941 cout << " ------------ BRProcA2C::run() END " << card2name_(card_)
[3645]942 << " ------------ " << endl;
[3647]943/*---- DELETE
944 {
945 nt.Info()["FirstTT"]=firsttt;
[3646]946 cout << nt;
947 sprintf(fname,"%s_nt.ppf",path_.c_str());
948 POutPersist po(fname);
949 po << PPFNameTag("ntv12") << nt;
[3683]950 cout << " BRProcA2C::run() created NTuple file " << fname << card2name_(card_) << endl;
[3647]951 }
952---- */
[3655]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;
[3683]974 cout << " BRProcA2C::run() created moysigspec file " << fname << card2name_(card_) << endl;
[3655]975 }
976
[3683]977 if (fghist_) {
978 delete ph1;
979 delete ph2;
980 }
[3640]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;
[3683]987 cout << " BRProcA2C::run()/Timing: " << card2name_(card_) << endl;
[3640]988 tm.Print();
989 cout << " ---------------------------------------------------------- " << endl;
990
[3635]991 }
992 catch (PException& exc) {
[3683]993 cout << " BRProcA2C::run()/catched PException " << exc.Msg() << endl;
[3635]994 setRC(3);
995 return;
996 }
997 catch(...) {
[3683]998 cout << " BRProcA2C::run()/catched unknown ... exception " << endl;
[3635]999 setRC(4);
1000 return;
1001 }
1002 setRC(0);
1003 return;
1004}
1005
[3645]1006//---------------------------------------------------------------------
[3683]1007// Classe thread de traitement 2 x 2 voies/frames (Apres BRProcA2C)
[3645]1008//---------------------------------------------------------------------
[3635]1009
[3645]1010/* --Methode-- */
[3683]1011BRProcB4C::BRProcB4C(RAcqMemZoneMgr& mem1, RAcqMemZoneMgr& mem2, string& path,
1012 bool fgraw, uint_4 nmean, uint_4 nmax, bool fgnotrl)
[3645]1013 : memgr1(mem1), memgr2(mem2)
1014{
[3683]1015 fgraw_ = fgraw;
[3645]1016 nmax_ = nmax;
1017 nmean_ = nmean;
[3683]1018 if (fgraw_) cout << " BRProcB4C::BRProcB4C() - constructeur RAW data - NMean= " << nmean_ << endl;
1019 else cout << " BRProcB4C::BRProcB4C() - constructeur FFT data - NMean= " << nmean_ << endl;
[3645]1020 stop_ = false;
1021 path_ = path;
1022 fgnotrl_ = fgnotrl;
1023}
[3635]1024
[3645]1025/* --Methode-- */
[3683]1026void BRProcB4C::Stop()
[3645]1027{
1028 stop_=true;
[3683]1029 // cout <<" BRProcB4C::Stop ... > STOP " << endl;
[3645]1030}
[3635]1031
[3645]1032
1033/* --Methode-- */
[3683]1034void BRProcB4C::run()
[3645]1035{
1036 setRC(1);
1037 try {
[3683]1038 Timer tm("BRProcB4C", false);
[3645]1039 TimeStamp ts;
[3646]1040 BRPaqChecker pcheck1(!fgnotrl_); // Verification/comptage des paquets
1041 BRPaqChecker pcheck2(!fgnotrl_); // Verification/comptage des paquets
[3645]1042
1043 size_t totnbytesout = 0;
1044 size_t totnbytesproc = 0;
1045
[3683]1046 cout << " BRProcB4C::run() - Starting " << ts << " NMaxMemZones=" << nmax_
[3645]1047 << " NMean=" << nmean_ << endl;
[3683]1048 cout << " BRProcB4C::run()... - Output Data Path: " << path_ << endl;
[3645]1049
1050 uint_4 paqsz = memgr1.PaqSize();
1051 uint_4 procpaqsz = memgr1.ProcPaqSize();
1052 if ((paqsz != memgr2.PaqSize())||(procpaqsz!= memgr2.ProcPaqSize())) {
[3683]1053 cout << "BRProcB4C::run()/ERROR : different paquet size -> stop \n ...(PaqSz1="
[3645]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
[3648]1061 BRPaquet pq(NULL, NULL, paqsz);
[3646]1062/*
[3645]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() );
[3646]1072*/
[3648]1073 int szfour = pq.DataSize()/2/2+1;
1074// int szfour = (paqsz-40)/2+1;
[3646]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;
[3683]1080 cout << " *DBG*BRProcB4C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz
[3646]1081 << " procpaqsz/2=" << procpaqsz/2 << " cfour.Size()=" << szfour
1082 << " *8=" << szfour*8 << endl;
[3645]1083
[3647]1084 DataTable dt;
1085 dt.AddLongColumn("fc1");
[3651]1086 dt.AddLongColumn("tt1");
[3647]1087 dt.AddLongColumn("fc2");
1088 dt.AddLongColumn("tt2");
1089 DataTableRow dtr = dt.EmptyRow();
1090
[3645]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;
[3647]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;
[3645]1104 for (uint_4 kmz=0; kmz<nmax_; kmz++) {
1105 uint_4 noksfc = 0;
1106 uint_4 nokpaq = 0;
1107 if (stop_) break;
[3646]1108 // cout << " *DBG*BBBBB" << kmz << endl;
1109
[3645]1110 int mid1 = memgr1.FindMemZoneId(MemZA_ProcB);
1111 Byte* buff1 = memgr1.GetMemZone(mid1);
1112 if (buff1 == NULL) {
[3683]1113 cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid1 << ") -> NULL" << endl;
[3645]1114 break;
1115 }
1116 Byte* procbuff1 = memgr1.GetProcMemZone(mid1);
1117 if (procbuff1 == NULL) {
[3683]1118 cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid1 << ") -> NULL" << endl;
[3645]1119 break;
1120 }
1121 int mid2 = memgr2.FindMemZoneId(MemZA_ProcB);
1122 Byte* buff2 = memgr2.GetMemZone(mid2);
1123 if (buff1 == NULL) {
[3683]1124 cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid2 << ") -> NULL" << endl;
[3645]1125 break;
1126 }
1127 Byte* procbuff2 = memgr2.GetProcMemZone(mid2);
1128 if (procbuff2 == NULL) {
[3683]1129 cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid2 << ") -> NULL" << endl;
[3645]1130 break;
1131 }
1132 uint_4 i1,i2;
1133 i1=i2=0;
[3646]1134// cout << " *DBG*CCCCCC " << kmz << " memgr1.NbPaquets() =" << memgr1.NbPaquets() << endl;
[3645]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
[3646]1150 if ((i1>=memgr1.NbPaquets())||(i2>=memgr1.NbPaquets())) {
1151 cout << " *BUG*BUG i1=" << i1 << " i2=" << i2 << endl;
1152 break;
1153 }
[3645]1154 // Les deux framecounters sont identiques ...
1155 noksfc++;
[3647]1156 curfc=paq1.FrameCounter();
1157 curfc2=paq2.FrameCounter();
1158 if (fgfirst) {
[3651]1159 firsttt=paq1.TimeTag(); firsttt2=paq2.TimeTag();
[3683]1160 cout << " BRProcB4C()/Info First FC="<<curfc<<" , "<<curfc2<<" -> TT="
[3647]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
[3645]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 }
[3647]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 }
[3645]1190 nzm++; i1++; i2++;
1191 totnbytesproc += 2*paq1.DataSize();
1192 } // Fin de boucle sur les paquets d'une zone
[3646]1193 memgr1.FreeMemZone(mid1, MemZS_ProcB);
1194 memgr2.FreeMemZone(mid2, MemZS_ProcB);
1195
[3645]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;
[3647]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;
[3645]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.);
[3646]1226 nzm = 0; ifile++;
[3645]1227// ts.SetNow();
1228// filog << ts << " : proc file " << fname << endl;
[3683]1229 cout << " BRProcB4C::run() created file " << fname << endl;
[3645]1230 }
1231 double okfrac = (nokpaq>1)?((double)noksfc/(double)nokpaq*100.):0.;
[3683]1232 cout << "BRProcB4C ["<<kmz<<"] NOKPaq=" << nokpaq << " NSameFC=" << noksfc
[3645]1233 << " (" << okfrac << " %)" << endl;
1234 totnokpaq += nokpaq;
1235 totnoksfc += noksfc;
1236 } // Fin de boucle sur les zones a traiter
[3683]1237 cout << " ------------------ BRProcB4C::run() END ----------------- " << endl;
[3647]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;
[3683]1246 cout << " BRProcB4C::run() created TimeTag/FrameCounter file " << fname << endl;
[3647]1247 }
[3645]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;
[3683]1257 cout << " BRProcB4C::run()/Timing: \n";
[3645]1258 tm.Print();
1259 cout << " ---------------------------------------------------------- " << endl;
1260}
1261 catch (PException& exc) {
[3683]1262 cout << " BRProcB4C::run()/catched PException " << exc.Msg() << endl;
[3645]1263 setRC(3);
1264 return;
1265 }
1266 catch(...) {
[3683]1267 cout << " BRProcB4C::run()/catched unknown ... exception " << endl;
[3645]1268 setRC(4);
1269 return;
1270 }
1271 setRC(0);
1272 return;
1273}
1274
1275
Note: See TracBrowser for help on using the repository browser.