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

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

Ajout numero d'identification a BRBaseProcesseur, Reza 18/05/2010

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