| 1 | //---------------------------------------------------------------- | 
|---|
| 2 | // Projet BAORadio - (C) LAL/IRFU  2008-2010 | 
|---|
| 3 | // Classes de threads de traitement pour BAORadio | 
|---|
| 4 | //---------------------------------------------------------------- | 
|---|
| 5 |  | 
|---|
| 6 | #include <stdlib.h> | 
|---|
| 7 | #include <string.h> | 
|---|
| 8 | #include <unistd.h> | 
|---|
| 9 | #include <fstream> | 
|---|
| 10 | #include <signal.h> | 
|---|
| 11 |  | 
|---|
| 12 | #include "pexceptions.h" | 
|---|
| 13 | #include "tvector.h" | 
|---|
| 14 | #include "ntuple.h" | 
|---|
| 15 | #include "datatable.h" | 
|---|
| 16 | #include "histos.h" | 
|---|
| 17 | #include "fioarr.h" | 
|---|
| 18 | #include "matharr.h" | 
|---|
| 19 | #include "timestamp.h" | 
|---|
| 20 | #include "ctimer.h" | 
|---|
| 21 | #include "fftpserver.h" | 
|---|
| 22 | #include "fitsarrhand.h" | 
|---|
| 23 |  | 
|---|
| 24 | #include "FFTW/fftw3.h" | 
|---|
| 25 |  | 
|---|
| 26 |  | 
|---|
| 27 | #include "pciewrap.h" | 
|---|
| 28 | #include "brpaqu.h" | 
|---|
| 29 | #include "brproc.h" | 
|---|
| 30 |  | 
|---|
| 31 |  | 
|---|
| 32 |  | 
|---|
| 33 | //--------------------------------------------------------------------- | 
|---|
| 34 | // Classe de traitement simple - calcul de spectres moyennes / voie | 
|---|
| 35 | //--------------------------------------------------------------------- | 
|---|
| 36 | /* --Methode-- */ | 
|---|
| 37 | BRMeanSpecCalculator::BRMeanSpecCalculator(RAcqMemZoneMgr& memgr, string outpath, uint_4 nmean, | 
|---|
| 38 | bool fgdatafft, bool fgsinglechan) | 
|---|
| 39 | : BRBaseProcessor(memgr), outpath_(outpath), nmean_(nmean), | 
|---|
| 40 | fgdatafft_(fgdatafft), fgsinglechannel_(fgsinglechan), | 
|---|
| 41 | clnflg_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()), | 
|---|
| 42 | nbadpaq_(fgsinglechan?memgr_.NbFibres():2*memgr_.NbFibres()) | 
|---|
| 43 | { | 
|---|
| 44 | setNameId("meanSpecCalc",1); | 
|---|
| 45 | BRPaquet paq(memgr_.PaqSize()); | 
|---|
| 46 | if (fgsinglechannel_) { | 
|---|
| 47 | mspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2); | 
|---|
| 48 | sigspecmtx_.SetSize(memgr_.NbFibres(), paq.DataSize()/2); | 
|---|
| 49 | sgain_.SetSize(memgr_.NbFibres(), paq.DataSize()/2); | 
|---|
| 50 | } | 
|---|
| 51 | else { | 
|---|
| 52 | mspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4); | 
|---|
| 53 | sigspecmtx_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4); | 
|---|
| 54 | sgain_.SetSize(2*memgr_.NbFibres(), paq.DataSize()/4); | 
|---|
| 55 | } | 
|---|
| 56 | mspecmtx_=(r_4)(0.); | 
|---|
| 57 | sigspecmtx_=(r_4)(0.); | 
|---|
| 58 | sgain_=(r_4)(1.);    // Gain en fonction de la frequence, à 1 par defaut | 
|---|
| 59 |  | 
|---|
| 60 | numfile_=0; | 
|---|
| 61 | nbpaq4mean_=0; | 
|---|
| 62 | totnbpaq_=0; | 
|---|
| 63 |  | 
|---|
| 64 | SetVarianceLimits(); | 
|---|
| 65 |  | 
|---|
| 66 | ofsdtp_=NULL; | 
|---|
| 67 | dtp_=NULL; | 
|---|
| 68 | string dtfile="!"+outpath_+"/dtspec.fits"; | 
|---|
| 69 | ofsdtp_ = new FitsInOutFile(dtfile,FitsInOutFile::Fits_Create); | 
|---|
| 70 | dtp_ = new SwFitsDataTable(*ofsdtp_,1024,true); | 
|---|
| 71 | int nchan=(fgsinglechannel_?memgr_.NbFibres():2*memgr_.NbFibres()); | 
|---|
| 72 | char cnom[32]; | 
|---|
| 73 | for(int i=0; i<nchan; i++) { | 
|---|
| 74 | sprintf(cnom,"var%d",i); | 
|---|
| 75 | dtp_->AddFloatColumn(cnom); | 
|---|
| 76 | } | 
|---|
| 77 | for(int i=0; i<nchan; i++) { | 
|---|
| 78 | sprintf(cnom,"varnorm%d",i); | 
|---|
| 79 | dtp_->AddFloatColumn(cnom); | 
|---|
| 80 | } | 
|---|
| 81 | xnt_=new double[nchan*2]; | 
|---|
| 82 | } | 
|---|
| 83 |  | 
|---|
| 84 | /* --Methode-- */ | 
|---|
| 85 | BRMeanSpecCalculator::~BRMeanSpecCalculator() | 
|---|
| 86 | { | 
|---|
| 87 | if (nbpaq4mean_>1)  SaveSpectra(); | 
|---|
| 88 | cout << " ---------------- BRMeanSpecCalculator()_Finalizing -------------------- " << endl; | 
|---|
| 89 | for(size_t i=0; i<nbadpaq_.size(); i++) { | 
|---|
| 90 | cout << " Channel " << i << "  NBadPaq=" << nbadpaq_[i] << " / TotNbPaq=" << totnbpaq_ << endl; | 
|---|
| 91 | } | 
|---|
| 92 | if (dtp_) { | 
|---|
| 93 | cout << *dtp_; | 
|---|
| 94 | delete dtp_; | 
|---|
| 95 | delete ofsdtp_; | 
|---|
| 96 | delete xnt_; | 
|---|
| 97 | } | 
|---|
| 98 | cout << " ------------------------------------------------------------------------ " << endl; | 
|---|
| 99 | } | 
|---|
| 100 |  | 
|---|
| 101 | /* --Methode-- */ | 
|---|
| 102 | void BRMeanSpecCalculator::ReadGainFitsFile(string filename) | 
|---|
| 103 | { | 
|---|
| 104 | cout << " BRMeanSpecCalculator::ReadGainFitsFile() - reading file " << filename; | 
|---|
| 105 | FitsInOutFile fis(filename, FitsInOutFile::Fits_RO); | 
|---|
| 106 | fis >> sgain_; | 
|---|
| 107 | cout << " MeanGain=" << sgain_.Sum()/sgain_.Size() << endl; | 
|---|
| 108 | } | 
|---|
| 109 |  | 
|---|
| 110 | static inline r_4 Zmod2(complex<r_4> z) | 
|---|
| 111 | { return (z.real()*z.real()+z.imag()*z.imag()); } | 
|---|
| 112 |  | 
|---|
| 113 | /* --Methode-- */ | 
|---|
| 114 | int BRMeanSpecCalculator::Process() | 
|---|
| 115 | { | 
|---|
| 116 | if (nbpaq4mean_==nmean_)  SaveSpectra(); | 
|---|
| 117 | FlagBadPackets(); | 
|---|
| 118 |  | 
|---|
| 119 | if (fgdatafft_) {  // Donnees firmware FFT | 
|---|
| 120 | for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) { | 
|---|
| 121 | if (clnflg_[i])  { nbadpaq_[i]++;  continue; }   // si le paquet a ete flagge mauvais ( clnflg_[i] <> 0 ) | 
|---|
| 122 | TwoByteComplex* zp=NULL; | 
|---|
| 123 | if (fgsinglechannel_) { | 
|---|
| 124 | zp=vpaq_[i].Data1C(); | 
|---|
| 125 | } | 
|---|
| 126 | else { | 
|---|
| 127 | zp=vpaq_[i/2].Data1C(); | 
|---|
| 128 | if (i%2==1)  zp=vpaq_[i/2].Data2C(); | 
|---|
| 129 | } | 
|---|
| 130 | TVector< r_4 > spec = mspecmtx_.Row(i); | 
|---|
| 131 | TVector< r_4 > sspec = sigspecmtx_.Row(i); | 
|---|
| 132 | for(sa_size_t f=1; f<spec.Size(); f++) { | 
|---|
| 133 | r_4 m2zf=zp[f].module2F();; | 
|---|
| 134 | spec(f) += m2zf; | 
|---|
| 135 | sspec(f) += m2zf*m2zf; | 
|---|
| 136 | } | 
|---|
| 137 | } | 
|---|
| 138 | } | 
|---|
| 139 | else {  // Donnees RAW qui ont du etre processe par BRFFTCalculator | 
|---|
| 140 | for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) { | 
|---|
| 141 | if (clnflg_[i])  { nbadpaq_[i]++;  continue; }   // si le paquet a ete flagge mauvais ( clnflg_[i] <> 0 ) | 
|---|
| 142 | complex<ODT>* zp=NULL; | 
|---|
| 143 | if (fgsinglechannel_) { | 
|---|
| 144 | zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i]); | 
|---|
| 145 | } | 
|---|
| 146 | else { | 
|---|
| 147 | zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i/2]); | 
|---|
| 148 | if (i%2==1)  zp= reinterpret_cast< complex<ODT>* >(vprocpaq_[i/2]+memgr_.ProcPaqSize()/2) ; | 
|---|
| 149 | } | 
|---|
| 150 | TVector< r_4 > spec = mspecmtx_.Row(i); | 
|---|
| 151 | TVector< r_4 > sspec = sigspecmtx_.Row(i); | 
|---|
| 152 | for(sa_size_t f=1; f<spec.Size(); f++)  { | 
|---|
| 153 | r_4 m2zf=Zmod2(zp[f]); | 
|---|
| 154 | spec(f) += m2zf; | 
|---|
| 155 | sspec(f) += (m2zf*m2zf); | 
|---|
| 156 | } | 
|---|
| 157 | } | 
|---|
| 158 | } | 
|---|
| 159 | nbpaq4mean_++;   totnbpaq_++; | 
|---|
| 160 | return 0; | 
|---|
| 161 | } | 
|---|
| 162 |  | 
|---|
| 163 | /* --Methode-- */ | 
|---|
| 164 | void BRMeanSpecCalculator::FlagBadPackets() | 
|---|
| 165 | { | 
|---|
| 166 | if (fgdatafft_) {  // Donnees firmware FFT | 
|---|
| 167 | for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) { | 
|---|
| 168 | TwoByteComplex* zp=NULL; | 
|---|
| 169 | if (fgsinglechannel_) { | 
|---|
| 170 | zp=vpaq_[i].Data1C(); | 
|---|
| 171 | } | 
|---|
| 172 | else { | 
|---|
| 173 | zp=vpaq_[i/2].Data1C(); | 
|---|
| 174 | if (i%2==1)  zp=vpaq_[i/2].Data2C(); | 
|---|
| 175 | } | 
|---|
| 176 | double variance=0.; | 
|---|
| 177 | double varnorm=0.; | 
|---|
| 178 | for(sa_size_t f=1; f<mspecmtx_.NCols(); f++) { | 
|---|
| 179 | double modsq=(double)zp[f].module2F(); | 
|---|
| 180 | variance+=modsq; | 
|---|
| 181 | varnorm+=modsq/sgain_(i,f); | 
|---|
| 182 | } | 
|---|
| 183 | xnt_[i]=variance; | 
|---|
| 184 | xnt_[i+mspecmtx_.NRows()]=varnorm; | 
|---|
| 185 | clnflg_[i]=0; | 
|---|
| 186 | if (varnorm<varmin_) clnflg_[i]=1; | 
|---|
| 187 | else if (varnorm>varmax_) clnflg_[i]=2; | 
|---|
| 188 | } | 
|---|
| 189 | } | 
|---|
| 190 | else {  // Donnees RAW qui ont du etre processe par BRFFTCalculator | 
|---|
| 191 | for(sa_size_t i=0; i<(size_t)mspecmtx_.NRows(); i++) { | 
|---|
| 192 | complex<ODT>* zp=NULL; | 
|---|
| 193 | if (fgsinglechannel_) { | 
|---|
| 194 | zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i]); | 
|---|
| 195 | } | 
|---|
| 196 | else { | 
|---|
| 197 | zp=reinterpret_cast< complex<ODT>* > (vprocpaq_[i/2]); | 
|---|
| 198 | if (i%2==1)  zp= reinterpret_cast< complex<ODT>* >(vprocpaq_[i/2]+memgr_.ProcPaqSize()/2) ; | 
|---|
| 199 | } | 
|---|
| 200 | double variance=0.; | 
|---|
| 201 | double varnorm=0.; | 
|---|
| 202 | for(sa_size_t f=1; f<mspecmtx_.NCols(); f++) { | 
|---|
| 203 | double modsq=(double)Zmod2(zp[f]); | 
|---|
| 204 | variance+=modsq; | 
|---|
| 205 | varnorm+=modsq/sgain_(i,f); | 
|---|
| 206 | } | 
|---|
| 207 | xnt_[i]=variance; | 
|---|
| 208 | xnt_[i+mspecmtx_.NRows()]=varnorm; | 
|---|
| 209 | clnflg_[i]=0; | 
|---|
| 210 | if (varnorm<varmin_) clnflg_[i]=1; | 
|---|
| 211 | else if (varnorm>varmax_) clnflg_[i]=2; | 
|---|
| 212 | } | 
|---|
| 213 | } | 
|---|
| 214 | dtp_->AddRow(xnt_); | 
|---|
| 215 | return; | 
|---|
| 216 | } | 
|---|
| 217 |  | 
|---|
| 218 | /* --Methode-- */ | 
|---|
| 219 | void BRMeanSpecCalculator::SaveSpectra() | 
|---|
| 220 | { | 
|---|
| 221 | mspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_; | 
|---|
| 222 | mspecmtx_ /= (double)nbpaq4mean_; | 
|---|
| 223 | sigspecmtx_.Info()["NPAQSUM"] = nbpaq4mean_; | 
|---|
| 224 | sigspecmtx_ /= (double)nbpaq4mean_; | 
|---|
| 225 | sigspecmtx_ -= (mspecmtx_ && mspecmtx_);  // Mean(X^2) - [ Mean(X) ]^2 | 
|---|
| 226 | char nfile[64]; | 
|---|
| 227 | string flnm; | 
|---|
| 228 | /* | 
|---|
| 229 | { | 
|---|
| 230 | sprintf(nfile,"mspecmtx%d.ppf",numfile_); | 
|---|
| 231 | flnm=outpath_+nfile; | 
|---|
| 232 | POutPersist po(flnm); | 
|---|
| 233 | po << mspecmtx_; | 
|---|
| 234 | } | 
|---|
| 235 | { | 
|---|
| 236 | sprintf(nfile,"sigspecmtx%d.ppf",numfile_); | 
|---|
| 237 | flnm=outpath_+nfile; | 
|---|
| 238 | POutPersist po(flnm); | 
|---|
| 239 | po << sigspecmtx_; | 
|---|
| 240 | } | 
|---|
| 241 | */ | 
|---|
| 242 | { | 
|---|
| 243 | sprintf(nfile,"mspecmtx%d.fits",numfile_); | 
|---|
| 244 | flnm="!"+outpath_+nfile; | 
|---|
| 245 | FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create); | 
|---|
| 246 | fos << mspecmtx_; | 
|---|
| 247 | } | 
|---|
| 248 | { | 
|---|
| 249 | sprintf(nfile,"sigspecmtx%d.fits",numfile_); | 
|---|
| 250 | flnm="!"+outpath_+nfile; | 
|---|
| 251 | FitsInOutFile fos(flnm,FitsInOutFile::Fits_Create); | 
|---|
| 252 | fos << sigspecmtx_; | 
|---|
| 253 | } | 
|---|
| 254 |  | 
|---|
| 255 | cout << numfile_ << "-BRMeanSpecCalculator::SaveSpectra() NPaqProc=" | 
|---|
| 256 | << totnbpaq_ << "  -> Mean/Sig spectra Matrix in " << flnm << " /sigspec...ppf" << endl; | 
|---|
| 257 | numfile_++;  nbpaq4mean_=0; | 
|---|
| 258 | mspecmtx_ = (r_4)(0.); | 
|---|
| 259 | sigspecmtx_ = (r_4)(0.); | 
|---|
| 260 | return; | 
|---|
| 261 | } | 
|---|
| 262 |  | 
|---|
| 263 | //--------------------------------------------------------------------- | 
|---|
| 264 | // Classe de thread de calcul de FFT sur donnees RAW | 
|---|
| 265 | //--------------------------------------------------------------------- | 
|---|
| 266 | /* --Methode-- */ | 
|---|
| 267 | BRFFTCalculator::BRFFTCalculator(RAcqMemZoneMgr& memgr, bool fgsinglechannel) | 
|---|
| 268 | : BRBaseProcessor(memgr), fgsinglechannel_(fgsinglechannel), totnbfftpaq_(0) | 
|---|
| 269 | { | 
|---|
| 270 | BRPaquet paq(memgr_.PaqSize()); | 
|---|
| 271 | setNameId("FFTCalc",2); | 
|---|
| 272 | ffts_.SetInDataSize((fgsinglechannel_)?paq.DataSize():paq.DataSize()/2); | 
|---|
| 273 | } | 
|---|
| 274 |  | 
|---|
| 275 | /* --Methode-- */ | 
|---|
| 276 | BRFFTCalculator::~BRFFTCalculator() | 
|---|
| 277 | { | 
|---|
| 278 | } | 
|---|
| 279 |  | 
|---|
| 280 |  | 
|---|
| 281 | /* --Methode-- */ | 
|---|
| 282 | int BRFFTCalculator::Process() | 
|---|
| 283 | { | 
|---|
| 284 | for(size_t fib=0; fib<(size_t)memgr_.NbFibres(); fib++) { | 
|---|
| 285 | ffts_.DoFFT( reinterpret_cast<IDT *>(vpaq_[fib].Data1() ), | 
|---|
| 286 | reinterpret_cast< complex<ODT>* > (vprocpaq_[fib]) ); | 
|---|
| 287 | totnbfftpaq_++; | 
|---|
| 288 | if ( fgsinglechannel_ ) continue; | 
|---|
| 289 | ffts_.DoFFT( reinterpret_cast<IDT *>(vpaq_[fib].Data2() ), | 
|---|
| 290 | reinterpret_cast< complex<ODT>* > (vprocpaq_[fib]+memgr_.ProcPaqSize()/2) ); | 
|---|
| 291 | totnbfftpaq_++; | 
|---|
| 292 | } | 
|---|
| 293 | return 0; | 
|---|
| 294 | } | 
|---|
| 295 |  | 
|---|
| 296 |  | 
|---|
| 297 | //------------------------------------------------------------------------- | 
|---|
| 298 | // Classe WBRFFT : Calcul de TF sur donnees brutes (firmware RAW) | 
|---|
| 299 | //------------------------------------------------------------------------- | 
|---|
| 300 | ZMutex* WBRFFT::mtx_fftwp_ = NULL; | 
|---|
| 301 |  | 
|---|
| 302 | /* --Methode-- */ | 
|---|
| 303 | WBRFFT::WBRFFT(uint_4 sz) | 
|---|
| 304 | : sz_(sz) | 
|---|
| 305 | { | 
|---|
| 306 | if (mtx_fftwp_ == NULL) mtx_fftwp_ = new ZMutex; | 
|---|
| 307 | if (sz>0)  SetInDataSize(sz); | 
|---|
| 308 | } | 
|---|
| 309 |  | 
|---|
| 310 | /* --Methode-- */ | 
|---|
| 311 | WBRFFT::~WBRFFT() | 
|---|
| 312 | { | 
|---|
| 313 | } | 
|---|
| 314 |  | 
|---|
| 315 | /* --Methode-- */ | 
|---|
| 316 | void WBRFFT::SetInDataSize(uint_4 sz) | 
|---|
| 317 | { | 
|---|
| 318 | sz_ = sz; | 
|---|
| 319 | if (sz_<1) return; | 
|---|
| 320 | inp.SetSize(sz); | 
|---|
| 321 | outfc.SetSize(sz/2+1); | 
|---|
| 322 | mtx_fftwp_->lock(); | 
|---|
| 323 | myplan_ = fftwf_plan_dft_r2c_1d(inp.Size(), inp.Data(), | 
|---|
| 324 | (fftwf_complex*)outfc.Data(), FFTW_ESTIMATE); | 
|---|
| 325 | mtx_fftwp_->unlock(); | 
|---|
| 326 | } | 
|---|
| 327 |  | 
|---|
| 328 | /* --Methode-- */ | 
|---|
| 329 | void WBRFFT::DoFFT( IDT *indata, complex<ODT> * ofc) | 
|---|
| 330 | { | 
|---|
| 331 | if (sz_<1) return; | 
|---|
| 332 | for(uint_4 k=0; k<inp.Size(); k++)   inp(k)=(ODT)indata[k]; | 
|---|
| 333 | fftwf_execute(myplan_); | 
|---|
| 334 | for(uint_4 k=0; k<outfc.Size(); k++)   ofc[k]=outfc(k); | 
|---|
| 335 | return; | 
|---|
| 336 | } | 
|---|
| 337 |  | 
|---|
| 338 | /* --Methode-- */ | 
|---|
| 339 | void WBRFFT::PrintData(IDT *indata,  complex<ODT> * ofc, uint_4 sz) | 
|---|
| 340 | { | 
|---|
| 341 | cout << " --- WBRFFT::PrintData() size=" << sz << endl; | 
|---|
| 342 | for(uint_4 k=0; k<sz; k+=8) { | 
|---|
| 343 | IDT* in = indata+k; | 
|---|
| 344 | cout << " Indata[" << k << "..." << k+8 << "]= "; | 
|---|
| 345 | for(uint_4 i=0; i<8; i++)   cout << (IIDT)in[i] << "  "; | 
|---|
| 346 | cout << endl; | 
|---|
| 347 | } | 
|---|
| 348 | cout << endl; | 
|---|
| 349 | for(uint_4 k=0; k<sz/2; k+=4) { | 
|---|
| 350 | complex< ODT>* out = ofc+k; | 
|---|
| 351 | cout << " OutFC[" << k << "..." << k+4 << "]= "; | 
|---|
| 352 | for(uint_4 i=0; i<4; i++)   cout << out[i] << "  "; | 
|---|
| 353 | cout << endl; | 
|---|
| 354 | } | 
|---|
| 355 | return; | 
|---|
| 356 |  | 
|---|
| 357 | } | 
|---|
| 358 |  | 
|---|
| 359 |  | 
|---|
| 360 | //--------------------------------------------------------------- | 
|---|
| 361 | // Classe thread de traitement donnees ADC avec 2 voies par frame | 
|---|
| 362 | //         !!!! OBSOLETE  !!!! | 
|---|
| 363 | //--------------------------------------------------------------- | 
|---|
| 364 |  | 
|---|
| 365 | // Mutex pour eviter le plantage du a FFTW qui ne semble pas thread-safe | 
|---|
| 366 | static ZMutex* pmutfftw=NULL; | 
|---|
| 367 |  | 
|---|
| 368 | /* --Methode-- */ | 
|---|
| 369 | BRProcA2C::BRProcA2C(RAcqMemZoneMgr& mem, string& path, bool fgraw, uint_4 nmean, | 
|---|
| 370 | uint_4 nmax, bool fghist, uint_4 nfsmap, bool fgnotrl, int card) | 
|---|
| 371 | :  memgr(mem) | 
|---|
| 372 | { | 
|---|
| 373 | fgraw_ = fgraw; | 
|---|
| 374 | nmax_ = nmax; | 
|---|
| 375 | nmean_ = nmean; | 
|---|
| 376 | if (fgraw_) cout << " BRProcA2C::BRProcA2C() - constructeur RAW data - NMean=" << nmean_ << endl; | 
|---|
| 377 | else cout << " BRProcA2C::BRProcA2C() - constructeur FFT data - NMean=" << nmean_ << endl; | 
|---|
| 378 | nfsmap_ = nfsmap; | 
|---|
| 379 | stop_ = false; | 
|---|
| 380 | path_ = path; | 
|---|
| 381 | fgnotrl_ = fgnotrl; | 
|---|
| 382 | fghist_ = fghist; | 
|---|
| 383 | card_ = card; | 
|---|
| 384 | if (pmutfftw==NULL) pmutfftw=new ZMutex; | 
|---|
| 385 | } | 
|---|
| 386 |  | 
|---|
| 387 | /* --Methode-- */ | 
|---|
| 388 | void BRProcA2C::Stop() | 
|---|
| 389 | { | 
|---|
| 390 | stop_=true; | 
|---|
| 391 | // cout <<" BRProcA2C::Stop ... > STOP " << endl; | 
|---|
| 392 | } | 
|---|
| 393 |  | 
|---|
| 394 |  | 
|---|
| 395 |  | 
|---|
| 396 |  | 
|---|
| 397 | static inline string card2name_(int card) | 
|---|
| 398 | { | 
|---|
| 399 | if (card==2) return " (Chan3,4) "; | 
|---|
| 400 | else return " (Chan1,2) "; | 
|---|
| 401 | } | 
|---|
| 402 | /* --Methode-- */ | 
|---|
| 403 | void BRProcA2C::run() | 
|---|
| 404 | { | 
|---|
| 405 | setRC(1); | 
|---|
| 406 | try { | 
|---|
| 407 | Timer tm("BRProcA2C", false); | 
|---|
| 408 | TimeStamp ts; | 
|---|
| 409 | BRPaqChecker pcheck(!fgnotrl_);  // Verification/comptage des paquets | 
|---|
| 410 |  | 
|---|
| 411 | size_t totnbytesout = 0; | 
|---|
| 412 | size_t totnbytesproc = 0; | 
|---|
| 413 |  | 
|---|
| 414 | cout << " BRProcA2C::run() - Starting " << ts << " NMaxMemZones=" << nmax_ | 
|---|
| 415 | << " NMean=" << nmean_ << card2name_(card_) << endl; | 
|---|
| 416 | cout << " BRProcA2C::run()... - Output Data Path: " << path_ << endl; | 
|---|
| 417 | char fname[512]; | 
|---|
| 418 | //    sprintf(fname,"%s/proc.log",path_.c_str()); | 
|---|
| 419 | //    ofstream filog(fname); | 
|---|
| 420 | //    filog << " BRProcA2C::run() - starting log file " << ts << endl; | 
|---|
| 421 | //    filog << " ... NMaxMemZones=" << nmax_ << " NMean=" << nmean_ << " Step=" << step_ << endl; | 
|---|
| 422 |  | 
|---|
| 423 | /*----DELETE   NTuple | 
|---|
| 424 | const char* nnames[8] = {"fcs","tts","s1","s2","s12","s12re","s12im","s12phi"}; | 
|---|
| 425 | NTuple nt(8, nnames); | 
|---|
| 426 | double xnt[10]; | 
|---|
| 427 | uint_4 nmnt = 0; | 
|---|
| 428 | double ms1,ms2,ms12,ms12re,ms12im,ms12phi; | 
|---|
| 429 | ----*/ | 
|---|
| 430 | // Time sample (raw data) /FFT coeff histograms | 
|---|
| 431 | Histo* ph1=NULL; | 
|---|
| 432 | Histo* ph2=NULL; | 
|---|
| 433 | if (fghist_) { | 
|---|
| 434 | if (fgraw_) { | 
|---|
| 435 | ph1 = new Histo(-0.5, 255.5, 256); | 
|---|
| 436 | ph2 = new Histo(-0.5, 255.5, 256); | 
|---|
| 437 | } | 
|---|
| 438 | else { | 
|---|
| 439 | ph1 = new Histo(-128.5, 128.5, 257); | 
|---|
| 440 | ph2 = new Histo(-128.5, 128.5, 257); | 
|---|
| 441 | } | 
|---|
| 442 | } | 
|---|
| 443 |  | 
|---|
| 444 | // Initialisation pour calcul FFT | 
|---|
| 445 | TVector< complex<r_4> > cfour1;  // composant TF | 
|---|
| 446 | uint_4 paqsz = memgr.PaqSize(); | 
|---|
| 447 | uint_4 procpaqsz = memgr.ProcPaqSize(); | 
|---|
| 448 |  | 
|---|
| 449 |  | 
|---|
| 450 | BRPaquet pq(NULL, NULL, paqsz); | 
|---|
| 451 | TVector<r_4> vx(pq.DataSize()/2); | 
|---|
| 452 | int szfour = pq.DataSize()/2/2+1; | 
|---|
| 453 | cfour1.SetSize(szfour); | 
|---|
| 454 | /* | 
|---|
| 455 | vx = (r_4)(0.); | 
|---|
| 456 | FFTPackServer ffts; | 
|---|
| 457 | ffts.FFTForward(vx, cfour1); | 
|---|
| 458 | szfour = cfour1.Size(); | 
|---|
| 459 | */ | 
|---|
| 460 |  | 
|---|
| 461 | bool fgtimfreq = false;  // true->cartes temps<>frequences | 
|---|
| 462 | if (nfsmap_>0) fgtimfreq=true; | 
|---|
| 463 |  | 
|---|
| 464 | TVector< complex<r_4> > cfour2(cfour1.Size()); | 
|---|
| 465 |  | 
|---|
| 466 | TVector<r_4> spectreV1(cfour1.Size()); | 
|---|
| 467 | TVector<r_4> spectreV2(cfour1.Size()); | 
|---|
| 468 | TVector<r_4> moyspecV1(cfour1.Size());   // Moyenne des Spectres | 
|---|
| 469 | TVector<r_4> moyspecV2(cfour1.Size()); | 
|---|
| 470 | TVector<r_4> sigspecV1(cfour1.Size());   // Sigma des Spectres | 
|---|
| 471 | TVector<r_4> sigspecV2(cfour1.Size()); | 
|---|
| 472 | TVector< complex<r_4> > visiV12( cfour1.Size() ); | 
|---|
| 473 |  | 
|---|
| 474 | TMatrix<r_4> timfreqV1, timfreqV2;   // Cartes temps<>frequences | 
|---|
| 475 | if (fgtimfreq) { | 
|---|
| 476 | timfreqV1.SetSize(nmean_, spectreV1.Size()/nfsmap_); | 
|---|
| 477 | timfreqV2.SetSize(nmean_, spectreV2.Size()/nfsmap_); | 
|---|
| 478 | } | 
|---|
| 479 | cout << " *DBG*BRProcA2C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz | 
|---|
| 480 | << " procpaqsz/2=" << procpaqsz/2 << " cfour1.Size()=" << cfour1.Size() | 
|---|
| 481 | << " *8="  << cfour1.Size()*8 << endl; | 
|---|
| 482 |  | 
|---|
| 483 | pmutfftw->lock(); | 
|---|
| 484 | fftwf_plan plan1 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(), | 
|---|
| 485 | (fftwf_complex*)cfour1.Data(), FFTW_ESTIMATE); | 
|---|
| 486 | fftwf_plan plan2 = fftwf_plan_dft_r2c_1d(vx.Size(), vx.Data(), | 
|---|
| 487 | (fftwf_complex*)cfour2.Data(), FFTW_ESTIMATE); | 
|---|
| 488 | pmutfftw->unlock(); | 
|---|
| 489 |  | 
|---|
| 490 | uint_4 ifile = 0; | 
|---|
| 491 | uint_4 nzm = 0;  // Nb de paquets moyennes pour le calcul de chaque spectre | 
|---|
| 492 | uint_4 nmoyspec = 0;  // Nb de spectres moyennes | 
|---|
| 493 |  | 
|---|
| 494 | uint_4 curfc=0; | 
|---|
| 495 | uint_8 curtt=0; | 
|---|
| 496 | uint_8 firsttt=0; | 
|---|
| 497 | bool fgfirst=true; | 
|---|
| 498 | double moysig[2]={0.,0.}; | 
|---|
| 499 | double sigsig[2]={0.,0.}; | 
|---|
| 500 | uint_8 nbsig[2]={0,0}; | 
|---|
| 501 |  | 
|---|
| 502 | for (uint_4 kmz=0; kmz<nmax_; kmz++) { | 
|---|
| 503 | if (stop_) break; | 
|---|
| 504 | int mid = memgr.FindMemZoneId(MemZA_ProcA); | 
|---|
| 505 | Byte* buff = memgr.GetMemZone(mid); | 
|---|
| 506 | if (buff == NULL) { | 
|---|
| 507 | cout << " BRProcA2C::run()/ERROR memgr.GetMemZone(" << mid << ") -> NULL" << endl; | 
|---|
| 508 | break; | 
|---|
| 509 | } | 
|---|
| 510 | Byte* procbuff = memgr.GetProcMemZone(mid); | 
|---|
| 511 | if (procbuff == NULL) { | 
|---|
| 512 | cout << " BRProcA2C::run()/ERROR memgr.GetProcMemZone(" << mid << ") -> NULL" << endl; | 
|---|
| 513 | break; | 
|---|
| 514 | } | 
|---|
| 515 | //---- DELETE      nmnt=0;  ms1=ms2=ms12=ms12re=ms12im=ms12phi=0.; | 
|---|
| 516 | for(uint_4 i=0; i<memgr.NbPaquets(); i++) { | 
|---|
| 517 | BRPaquet paq(NULL, buff+i*paqsz, paqsz); | 
|---|
| 518 | if (!pcheck.Check(paq)) continue;   // on ne traite que les paquets OK | 
|---|
| 519 | if (fgfirst) { firsttt=paq.TimeTag(); fgfirst=false; } | 
|---|
| 520 | curfc=paq.FrameCounter(); | 
|---|
| 521 | curtt=paq.TimeTag()-firsttt; | 
|---|
| 522 | // Traitement voie 1 | 
|---|
| 523 | if (fghist_) { | 
|---|
| 524 | for(sa_size_t j=0; j<vx.Size(); j++) { | 
|---|
| 525 | r_4 vts=(fgraw_)?((r_4)(*(paq.Data1()+j))):((r_4)(*(paq.Data1S()+j))); | 
|---|
| 526 | ph1->Add((r_8)vts); | 
|---|
| 527 | moysig[0] += (double)vts; | 
|---|
| 528 | sigsig[0] += ((double)vts)*((double)vts); | 
|---|
| 529 | nbsig[0]++; | 
|---|
| 530 | } | 
|---|
| 531 | for(sa_size_t j=0; j<vx.Size(); j++) { | 
|---|
| 532 | r_4 vts=(fgraw_)?((r_4)(*(paq.Data2()+j))):((r_4)(*(paq.Data2S()+j))); | 
|---|
| 533 | ph2->Add((r_8)vts); | 
|---|
| 534 | moysig[1] += (double)vts; | 
|---|
| 535 | sigsig[1] += ((double)vts)*((double)vts); | 
|---|
| 536 | nbsig[1]++; | 
|---|
| 537 | } | 
|---|
| 538 | } | 
|---|
| 539 | if (fgraw_) { | 
|---|
| 540 | for(sa_size_t j=0; j<vx.Size(); j++) | 
|---|
| 541 | vx(j) = (r_4)(*(paq.Data1()+j))-127.5; | 
|---|
| 542 | //        fftwf_complex* coeff1 = (fftwf_complex*)(procbuff+i*procpaqsz); | 
|---|
| 543 | fftwf_execute(plan1); | 
|---|
| 544 | // Traitement voie 2 | 
|---|
| 545 | for(sa_size_t j=0; j<vx.Size(); j++) | 
|---|
| 546 | vx(j) = (r_4)(*(paq.Data2()+j))-127.5; | 
|---|
| 547 | fftwf_execute(plan2); | 
|---|
| 548 | } | 
|---|
| 549 | else { | 
|---|
| 550 | for(sa_size_t j=1; j<cfour1.Size()-1; j++) { | 
|---|
| 551 | cfour1(j) = complex<r_4>((r_4)paq.Data1C()[j].realB(), (r_4)paq.Data1C()[j].imagB()); | 
|---|
| 552 | cfour2(j) = complex<r_4>((r_4)paq.Data2C()[j].realB(), (r_4)paq.Data2C()[j].imagB()); | 
|---|
| 553 | } | 
|---|
| 554 | cfour1(0) = complex<r_4>((r_4)paq.Data1C()[0].realB(), (r_4)0.); | 
|---|
| 555 | cfour1(cfour1.Size()-1) = complex<r_4>((r_4)paq.Data1C()[0].imagB(), (r_4)0.); | 
|---|
| 556 | cfour2(0) = complex<r_4>((r_4)paq.Data2C()[0].realB(), (r_4)0.); | 
|---|
| 557 | cfour2(cfour2.Size()-1) = complex<r_4>((r_4)paq.Data2C()[0].imagB(), (r_4)0.); | 
|---|
| 558 | } | 
|---|
| 559 | for(sa_size_t j=0; j<spectreV1.Size(); j++) | 
|---|
| 560 | spectreV1(j) += Zmod2(cfour1(j)); | 
|---|
| 561 | memcpy(procbuff+i*procpaqsz, cfour1.Data(), sizeof(complex<r_4>)*cfour1.Size()); | 
|---|
| 562 | if (fgtimfreq) {   // Remplissage tableau temps-frequence | 
|---|
| 563 | for(sa_size_t c=1; c<timfreqV1.NCols(); c++) { | 
|---|
| 564 | for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++) | 
|---|
| 565 | timfreqV1(nzm, c) += Zmod2(cfour1(j)); | 
|---|
| 566 | } | 
|---|
| 567 | } | 
|---|
| 568 | for(sa_size_t j=0; j<spectreV2.Size(); j++) | 
|---|
| 569 | spectreV2(j) += Zmod2(cfour2(j));  // Zmod2(zp2[j]); | 
|---|
| 570 | memcpy(procbuff+i*procpaqsz+procpaqsz/2, cfour2.Data(), sizeof(complex<r_4>)*cfour2.Size()); | 
|---|
| 571 | if (fgtimfreq) {   // Remplissage tableau temps-frequence | 
|---|
| 572 | for(sa_size_t c=1; c<timfreqV2.NCols(); c++) { | 
|---|
| 573 | for(sa_size_t j=c*nfsmap_; j<(c+1)*nfsmap_; j++) | 
|---|
| 574 | timfreqV2(nzm,c) += Zmod2(cfour2(j)); | 
|---|
| 575 | } | 
|---|
| 576 | } | 
|---|
| 577 |  | 
|---|
| 578 | // Calcul correlation (visibilite V1 * V2) | 
|---|
| 579 | for(sa_size_t j=0; j<visiV12.Size(); j++) | 
|---|
| 580 | visiV12(j)+=cfour1(j)*conj(cfour2(j)); | 
|---|
| 581 | //        for(sa_size_t j=0; j<visiV12.Size(); j++) visiV12(j)+=zp1[j]*zp2[j]; | 
|---|
| 582 | if (nzm==0) { | 
|---|
| 583 | spectreV1.Info()["StartFC"] = curfc; | 
|---|
| 584 | spectreV2.Info()["StartFC"] = curfc; | 
|---|
| 585 | visiV12.Info()["StartFC"] = curfc; | 
|---|
| 586 | spectreV1.Info()["StartTT"] = curtt; | 
|---|
| 587 | spectreV2.Info()["StartTT"] = curtt; | 
|---|
| 588 | visiV12.Info()["StartTT"] = curtt; | 
|---|
| 589 | } | 
|---|
| 590 | nzm++; | 
|---|
| 591 | /*----DELETE | 
|---|
| 592 | if (nmnt==0)  { xnt[0]=paq.FrameCounter();  xnt[1]=paq.TimeTag(); } | 
|---|
| 593 | for(sa_size_t j=2700; j<2800; j++) { | 
|---|
| 594 | ms1 += Zmod2(cfour1(j)); ms2 += Zmod2(cfour2(j)); | 
|---|
| 595 | complex<r_4> zvis =  cfour1(j)*conj(cfour2(j)); | 
|---|
| 596 | ms12 += Zmod2(zvis);   ms12re += zvis.real();  ms12im += zvis.imag(); | 
|---|
| 597 | ms12phi+= atan2(zvis.imag(),zvis.real()); | 
|---|
| 598 | } | 
|---|
| 599 | nmnt++; | 
|---|
| 600 | ----*/ | 
|---|
| 601 | totnbytesproc += paq.DataSize(); | 
|---|
| 602 | totnbytesout += (2*sizeof(complex<r_4>)*cfour1.Size()); | 
|---|
| 603 |  | 
|---|
| 604 | } // Fin de boucle sur les paquets d'une zone | 
|---|
| 605 |  | 
|---|
| 606 | /*---- DELETE | 
|---|
| 607 | if (nmnt>0)  { | 
|---|
| 608 | double fnorm = (double)nmnt*(2800-2700); | 
|---|
| 609 | xnt[2] = ms1 /= fnorm; | 
|---|
| 610 | xnt[3] = ms2 /= fnorm; | 
|---|
| 611 | xnt[4] = ms12 /= fnorm; | 
|---|
| 612 | xnt[5] = ms12re /= fnorm; | 
|---|
| 613 | xnt[6] = ms12im /= fnorm; | 
|---|
| 614 | xnt[7] = ms12phi /= fnorm; | 
|---|
| 615 | nt.Fill(xnt); | 
|---|
| 616 | } | 
|---|
| 617 | ----*/ | 
|---|
| 618 | if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) { | 
|---|
| 619 | spectreV1 /= (r_4)(nzm); | 
|---|
| 620 | spectreV2 /= (r_4)(nzm); | 
|---|
| 621 |  | 
|---|
| 622 | // pour le calcul des moyennes et sigmas de ces spectres | 
|---|
| 623 | moyspecV1 += spectreV1; | 
|---|
| 624 | moyspecV2 += spectreV2; | 
|---|
| 625 | sigspecV1 += (spectreV1 && spectreV1); | 
|---|
| 626 | sigspecV2 += (spectreV2 && spectreV2); | 
|---|
| 627 | nmoyspec++; | 
|---|
| 628 |  | 
|---|
| 629 | visiV12 /= complex<r_4>((r_4)nzm, 0.); | 
|---|
| 630 |  | 
|---|
| 631 | spectreV1.Info()["NPaqMoy"] = nzm; | 
|---|
| 632 | spectreV2.Info()["NPaqMoy"] = nzm; | 
|---|
| 633 | visiV12.Info()["NPaqMoy"] = nzm; | 
|---|
| 634 | spectreV1.Info()["EndFC"] = curfc; | 
|---|
| 635 | spectreV2.Info()["EndFC"] = curfc; | 
|---|
| 636 | visiV12.Info()["EndFC"] = curfc; | 
|---|
| 637 | spectreV1.Info()["EndTT"] = curtt; | 
|---|
| 638 | spectreV2.Info()["EndTT"] = curtt; | 
|---|
| 639 | visiV12.Info()["EndTT"] = curtt; | 
|---|
| 640 | { | 
|---|
| 641 | sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile); | 
|---|
| 642 | POutPersist po(fname); | 
|---|
| 643 | string tag1="specV1"; | 
|---|
| 644 | string tag2="specV2"; | 
|---|
| 645 | string tag12="visiV12"; | 
|---|
| 646 | string tagh1="tshV1"; | 
|---|
| 647 | string tagh2="tshV2"; | 
|---|
| 648 | string tagtf1="timfreqV1"; | 
|---|
| 649 | string tagtf2="timfreqV2"; | 
|---|
| 650 | if (card_==2) { | 
|---|
| 651 | tag1 = "specV3"; | 
|---|
| 652 | tag2 = "specV4"; | 
|---|
| 653 | tagh1 = "tshV1"; | 
|---|
| 654 | tagh2 = "tshV2"; | 
|---|
| 655 | tag12="visiV34"; | 
|---|
| 656 | tagtf1="timfreqV3"; | 
|---|
| 657 | tagtf2="timfreqV4"; | 
|---|
| 658 | } | 
|---|
| 659 | po << PPFNameTag(tag1) << spectreV1; | 
|---|
| 660 | po << PPFNameTag(tag2) << spectreV2; | 
|---|
| 661 | po << PPFNameTag(tag12) << visiV12; | 
|---|
| 662 | if (fghist_) { | 
|---|
| 663 | po << PPFNameTag(tagh1) << (*ph1); | 
|---|
| 664 | po << PPFNameTag(tagh2) << (*ph2); | 
|---|
| 665 |  | 
|---|
| 666 | double sspvmax[3] = {0.,0.,0.}; | 
|---|
| 667 | int_4 sspvmaxidx[3] = {-1,-1,-1}; | 
|---|
| 668 | for(int jji=1;jji<visiV12.Size()-1;jji++) { | 
|---|
| 669 | r_4 zmv2 = Zmod2(visiV12(jji)); | 
|---|
| 670 | if (zmv2>sspvmax[2]) { sspvmax[2]=zmv2; sspvmaxidx[2]=jji; } | 
|---|
| 671 | } | 
|---|
| 672 | TVector<r_4>& sspv = spectreV1; | 
|---|
| 673 | for(int ic=0; ic<2; ic++) { | 
|---|
| 674 | if (ic==1) sspv = spectreV2; | 
|---|
| 675 | for(int jji=1;jji<sspv.Size()-1;jji++) | 
|---|
| 676 | if (sspv(jji)>sspvmax[ic]) { sspvmax[ic]=sspv(jji); sspvmaxidx[ic]=jji; } | 
|---|
| 677 | if (nbsig[ic] < 1) { moysig[ic]=sigsig[ic]=-1.;  } | 
|---|
| 678 | else { | 
|---|
| 679 | moysig[ic] /= (double)nbsig[ic]; | 
|---|
| 680 | sigsig[ic] /= (double)nbsig[ic]; | 
|---|
| 681 | sigsig[ic] -= (moysig[ic]*moysig[ic]); | 
|---|
| 682 | sigsig[ic] = sqrt(sigsig[ic]); | 
|---|
| 683 | cout << "===Voie " << ic << " Moy=" << moysig[ic] << " Sig=" << sigsig[ic] | 
|---|
| 684 | << " MaxSpec Amp= " << sqrt(sspvmax[ic])/double(pq.DataSize()/2/2) | 
|---|
| 685 | << " Pos=" << sspvmaxidx[ic] << "  (NPts=" << nbsig[ic] << ")" << endl; | 
|---|
| 686 | } | 
|---|
| 687 | } | 
|---|
| 688 | cout << "=== Voie1x2 MaxSpec Amp= " << sqrt(sqrt(sspvmax[2])/double(pq.DataSize()/2/2)) | 
|---|
| 689 | << " Pos=" << sspvmaxidx[2] << endl; | 
|---|
| 690 | }  // fin if (fghist_) | 
|---|
| 691 |  | 
|---|
| 692 | if (fgtimfreq) { | 
|---|
| 693 | timfreqV1 /= (r_4)nzm; | 
|---|
| 694 | timfreqV2 /= (r_4)nzm; | 
|---|
| 695 | po << PPFNameTag(tagtf1) << timfreqV1; | 
|---|
| 696 | po << PPFNameTag(tagtf2) << timfreqV2; | 
|---|
| 697 | } | 
|---|
| 698 | } | 
|---|
| 699 | spectreV1 = (r_4)(0.); | 
|---|
| 700 | spectreV2 = (r_4)(0.); | 
|---|
| 701 | visiV12 = complex<r_4>(0., 0.); | 
|---|
| 702 | if (fghist_) { | 
|---|
| 703 | ph1->Zero(); | 
|---|
| 704 | ph2->Zero(); | 
|---|
| 705 | moysig[0]=moysig[1]=0.; | 
|---|
| 706 | sigsig[0]=sigsig[1]=0.; | 
|---|
| 707 | nbsig[0]=nbsig[1]=0; | 
|---|
| 708 | } | 
|---|
| 709 | if (fgtimfreq) { | 
|---|
| 710 | timfreqV1 = (r_4)(0.); | 
|---|
| 711 | timfreqV2 = (r_4)(0.); | 
|---|
| 712 | } | 
|---|
| 713 | nzm = 0;  ifile++; | 
|---|
| 714 | //        ts.SetNow(); | 
|---|
| 715 | //        filog << ts << " :  proc file  " << fname << endl; | 
|---|
| 716 | cout << " BRProcA2C::run() created file  " << fname << card2name_(card_) << endl; | 
|---|
| 717 | } | 
|---|
| 718 |  | 
|---|
| 719 | memgr.FreeMemZone(mid, MemZS_ProcA); | 
|---|
| 720 | }  // Fin de boucle sur les zones a traiter | 
|---|
| 721 | cout << " ------------  BRProcA2C::run() END " << card2name_(card_) | 
|---|
| 722 | << " ------------ " << endl; | 
|---|
| 723 | /*---- DELETE | 
|---|
| 724 | { | 
|---|
| 725 | nt.Info()["FirstTT"]=firsttt; | 
|---|
| 726 | cout << nt; | 
|---|
| 727 | sprintf(fname,"%s_nt.ppf",path_.c_str()); | 
|---|
| 728 | POutPersist po(fname); | 
|---|
| 729 | po << PPFNameTag("ntv12") << nt; | 
|---|
| 730 | cout << " BRProcA2C::run() created NTuple file " << fname << card2name_(card_) << endl; | 
|---|
| 731 | } | 
|---|
| 732 | ---- */ | 
|---|
| 733 | if (nmoyspec>0) {  // Calcul des moyennes et sigmas des spectres | 
|---|
| 734 | r_4 fnms = nmoyspec; | 
|---|
| 735 | moyspecV1 /= fnms; | 
|---|
| 736 | moyspecV2 /= fnms; | 
|---|
| 737 | sigspecV1 /= fnms; | 
|---|
| 738 | sigspecV2 /= fnms; | 
|---|
| 739 | sigspecV1 -= (moyspecV1 && moyspecV1); | 
|---|
| 740 | sigspecV2 -= (moyspecV2 && moyspecV2); | 
|---|
| 741 | sigspecV1 = Sqrt(sigspecV1); | 
|---|
| 742 | sigspecV2 = Sqrt(sigspecV2); | 
|---|
| 743 | TVector<r_4> rsbV1, rsbV2;   // Rapport signal/bruit | 
|---|
| 744 | moyspecV1.DivElt(sigspecV1, rsbV1, false, true); | 
|---|
| 745 | moyspecV2.DivElt(sigspecV2, rsbV2, false, true); | 
|---|
| 746 | sprintf(fname,"%s_ms.ppf",path_.c_str()); | 
|---|
| 747 | POutPersist po(fname); | 
|---|
| 748 | po << PPFNameTag("moyspecV1") << moyspecV1; | 
|---|
| 749 | po << PPFNameTag("moyspecV2") << moyspecV2; | 
|---|
| 750 | po << PPFNameTag("sigspecV1") << sigspecV1; | 
|---|
| 751 | po << PPFNameTag("sigspecV2") << sigspecV2; | 
|---|
| 752 | po << PPFNameTag("rsbV1") << rsbV1; | 
|---|
| 753 | po << PPFNameTag("rsbV2") << rsbV2; | 
|---|
| 754 | cout << " BRProcA2C::run() created moysigspec file " << fname << card2name_(card_) << endl; | 
|---|
| 755 | } | 
|---|
| 756 |  | 
|---|
| 757 | if (fghist_) { | 
|---|
| 758 | delete ph1; | 
|---|
| 759 | delete ph2; | 
|---|
| 760 | } | 
|---|
| 761 | ts.SetNow(); | 
|---|
| 762 | tm.SplitQ(); | 
|---|
| 763 | cout << "  TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= " | 
|---|
| 764 | << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s" | 
|---|
| 765 | << " ProcDataOut=" <<  totnbytesout/(1024*1024) << " MB" << endl; | 
|---|
| 766 | cout << pcheck; | 
|---|
| 767 | cout << " BRProcA2C::run()/Timing: " << card2name_(card_) << endl; | 
|---|
| 768 | tm.Print(); | 
|---|
| 769 | cout << " ---------------------------------------------------------- " << endl; | 
|---|
| 770 |  | 
|---|
| 771 | } | 
|---|
| 772 | catch (PException& exc) { | 
|---|
| 773 | cout << " BRProcA2C::run()/catched PException " << exc.Msg() << endl; | 
|---|
| 774 | setRC(3); | 
|---|
| 775 | return; | 
|---|
| 776 | } | 
|---|
| 777 | catch(...) { | 
|---|
| 778 | cout << " BRProcA2C::run()/catched unknown ... exception " << endl; | 
|---|
| 779 | setRC(4); | 
|---|
| 780 | return; | 
|---|
| 781 | } | 
|---|
| 782 | setRC(0); | 
|---|
| 783 | return; | 
|---|
| 784 | } | 
|---|
| 785 |  | 
|---|
| 786 |  | 
|---|
| 787 | //--------------------------------------------------------------------- | 
|---|
| 788 | // Classe thread de traitement 2 x 2 voies/frames (Apres BRProcA2C) | 
|---|
| 789 | //         !!!! OBSOLETE  !!!! | 
|---|
| 790 | //--------------------------------------------------------------------- | 
|---|
| 791 |  | 
|---|
| 792 | /* --Methode-- */ | 
|---|
| 793 | BRProcB4C::BRProcB4C(RAcqMemZoneMgr& mem1, RAcqMemZoneMgr& mem2, string& path, | 
|---|
| 794 | bool fgraw, uint_4 nmean, uint_4 nmax, bool fgnotrl) | 
|---|
| 795 | :  memgr1(mem1), memgr2(mem2) | 
|---|
| 796 | { | 
|---|
| 797 | fgraw_ = fgraw; | 
|---|
| 798 | nmax_ = nmax; | 
|---|
| 799 | nmean_ = nmean; | 
|---|
| 800 | if (fgraw_) cout << " BRProcB4C::BRProcB4C() - constructeur RAW data - NMean= " << nmean_ << endl; | 
|---|
| 801 | else cout << " BRProcB4C::BRProcB4C() - constructeur FFT data - NMean= " << nmean_ << endl; | 
|---|
| 802 | stop_ = false; | 
|---|
| 803 | path_ = path; | 
|---|
| 804 | fgnotrl_ = fgnotrl; | 
|---|
| 805 | } | 
|---|
| 806 |  | 
|---|
| 807 | /* --Methode-- */ | 
|---|
| 808 | void BRProcB4C::Stop() | 
|---|
| 809 | { | 
|---|
| 810 | stop_=true; | 
|---|
| 811 | // cout <<" BRProcB4C::Stop ... > STOP " << endl; | 
|---|
| 812 | } | 
|---|
| 813 |  | 
|---|
| 814 |  | 
|---|
| 815 | /* --Methode-- */ | 
|---|
| 816 | void BRProcB4C::run() | 
|---|
| 817 | { | 
|---|
| 818 | setRC(1); | 
|---|
| 819 | try { | 
|---|
| 820 | Timer tm("BRProcB4C", false); | 
|---|
| 821 | TimeStamp ts; | 
|---|
| 822 | BRPaqChecker pcheck1(!fgnotrl_);  // Verification/comptage des paquets | 
|---|
| 823 | BRPaqChecker pcheck2(!fgnotrl_);  // Verification/comptage des paquets | 
|---|
| 824 |  | 
|---|
| 825 | size_t totnbytesout = 0; | 
|---|
| 826 | size_t totnbytesproc = 0; | 
|---|
| 827 |  | 
|---|
| 828 | cout << " BRProcB4C::run() - Starting " << ts << " NMaxMemZones=" << nmax_ | 
|---|
| 829 | << " NMean=" << nmean_ << endl; | 
|---|
| 830 | cout << " BRProcB4C::run()... - Output Data Path: " << path_ << endl; | 
|---|
| 831 |  | 
|---|
| 832 | uint_4 paqsz = memgr1.PaqSize(); | 
|---|
| 833 | uint_4 procpaqsz = memgr1.ProcPaqSize(); | 
|---|
| 834 | if ((paqsz != memgr2.PaqSize())||(procpaqsz!= memgr2.ProcPaqSize())) { | 
|---|
| 835 | cout << "BRProcB4C::run()/ERROR : different paquet size -> stop \n ...(PaqSz1=" | 
|---|
| 836 | << paqsz << " Sz2=" << memgr2.PaqSize() << " ProcPaqSz1=" | 
|---|
| 837 | << procpaqsz << " Sz2=" << memgr2.ProcPaqSize() << " )" << endl; | 
|---|
| 838 | setRC(9); | 
|---|
| 839 | return; | 
|---|
| 840 | } | 
|---|
| 841 |  | 
|---|
| 842 | TVector< complex<r_4> > cfour;  // composant TF | 
|---|
| 843 | BRPaquet pq(NULL, NULL, paqsz); | 
|---|
| 844 | /* | 
|---|
| 845 | TVector<r_4> vx(pq.DataSize()/2); | 
|---|
| 846 | vx = (r_4)(0.); | 
|---|
| 847 | FFTPackServer ffts; | 
|---|
| 848 | ffts.FFTForward(vx, cfour); | 
|---|
| 849 |  | 
|---|
| 850 | TVector< complex<r_4> > visiV13( cfour.Size() ); | 
|---|
| 851 | TVector< complex<r_4> > visiV14( cfour.Size() ); | 
|---|
| 852 | TVector< complex<r_4> > visiV23( cfour.Size() ); | 
|---|
| 853 | TVector< complex<r_4> > visiV24( cfour.Size() ); | 
|---|
| 854 | */ | 
|---|
| 855 | int szfour = pq.DataSize()/2/2+1; | 
|---|
| 856 | //    int szfour = (paqsz-40)/2+1; | 
|---|
| 857 | TVector< complex<r_4> > visiV13( szfour  ); | 
|---|
| 858 | TVector< complex<r_4> > visiV14( szfour ); | 
|---|
| 859 | TVector< complex<r_4> > visiV23( szfour ); | 
|---|
| 860 | TVector< complex<r_4> > visiV24( szfour ); | 
|---|
| 861 | // cout << " *DBG*AAAAA ---- Vectors OK" << endl; | 
|---|
| 862 | cout << " *DBG*BRProcB4C PaqSz=" << paqsz << " ProcPaqSize=" << procpaqsz | 
|---|
| 863 | << " procpaqsz/2=" << procpaqsz/2 << " cfour.Size()=" << szfour | 
|---|
| 864 | << " *8="  << szfour*8 << endl; | 
|---|
| 865 |  | 
|---|
| 866 | DataTable dt; | 
|---|
| 867 | dt.AddLongColumn("fc1"); | 
|---|
| 868 | dt.AddLongColumn("tt1"); | 
|---|
| 869 | dt.AddLongColumn("fc2"); | 
|---|
| 870 | dt.AddLongColumn("tt2"); | 
|---|
| 871 | DataTableRow dtr = dt.EmptyRow(); | 
|---|
| 872 |  | 
|---|
| 873 | uint_4 nzm = 0; | 
|---|
| 874 | uint_4 totnoksfc = 0; | 
|---|
| 875 | uint_4 totnokpaq = 0; | 
|---|
| 876 | uint_4 totnpaq = 0; | 
|---|
| 877 | uint_4 ifile = 0; | 
|---|
| 878 |  | 
|---|
| 879 | uint_4 curfc=0; | 
|---|
| 880 | uint_8 curtt=0; | 
|---|
| 881 | uint_4 curfc2=0; | 
|---|
| 882 | uint_8 curtt2=0; | 
|---|
| 883 | uint_8 firsttt=0; | 
|---|
| 884 | uint_8 firsttt2=0; | 
|---|
| 885 | bool fgfirst=true; | 
|---|
| 886 | for (uint_4 kmz=0; kmz<nmax_; kmz++) { | 
|---|
| 887 | uint_4 noksfc = 0; | 
|---|
| 888 | uint_4 nokpaq = 0; | 
|---|
| 889 | if (stop_) break; | 
|---|
| 890 | // cout << " *DBG*BBBBB" << kmz << endl; | 
|---|
| 891 |  | 
|---|
| 892 | int mid1 = memgr1.FindMemZoneId(MemZA_ProcB); | 
|---|
| 893 | Byte* buff1 = memgr1.GetMemZone(mid1); | 
|---|
| 894 | if (buff1 == NULL) { | 
|---|
| 895 | cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid1 << ") -> NULL" << endl; | 
|---|
| 896 | break; | 
|---|
| 897 | } | 
|---|
| 898 | Byte* procbuff1 = memgr1.GetProcMemZone(mid1); | 
|---|
| 899 | if (procbuff1 == NULL) { | 
|---|
| 900 | cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid1 << ") -> NULL" << endl; | 
|---|
| 901 | break; | 
|---|
| 902 | } | 
|---|
| 903 | int mid2 = memgr2.FindMemZoneId(MemZA_ProcB); | 
|---|
| 904 | Byte* buff2 = memgr2.GetMemZone(mid2); | 
|---|
| 905 | if (buff1 == NULL) { | 
|---|
| 906 | cout << " BRProcB4C::run()/ERROR memgr.GetMemZone(" << mid2 << ") -> NULL" << endl; | 
|---|
| 907 | break; | 
|---|
| 908 | } | 
|---|
| 909 | Byte* procbuff2 = memgr2.GetProcMemZone(mid2); | 
|---|
| 910 | if (procbuff2 == NULL) { | 
|---|
| 911 | cout << " BRProcB4C::run()/ERROR memgr.GetProcMemZone(" << mid2 << ") -> NULL" << endl; | 
|---|
| 912 | break; | 
|---|
| 913 | } | 
|---|
| 914 | uint_4 i1,i2; | 
|---|
| 915 | i1=i2=0; | 
|---|
| 916 | //      cout << " *DBG*CCCCCC " << kmz << " memgr1.NbPaquets() =" << memgr1.NbPaquets() << endl; | 
|---|
| 917 | while((i1<memgr1.NbPaquets())&&(i2<memgr2.NbPaquets())) { | 
|---|
| 918 | BRPaquet paq1(NULL, buff1+i1*paqsz, paqsz); | 
|---|
| 919 | BRPaquet paq2(NULL, buff2+i2*paqsz, paqsz); | 
|---|
| 920 | totnpaq++; | 
|---|
| 921 | //        cout << " DBG["<<kmz<<"] i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter() | 
|---|
| 922 | //<<","<<paq2.FrameCounter()<<endl; | 
|---|
| 923 | // on ne traite que les paquets OK | 
|---|
| 924 | if (!pcheck1.Check(paq1)) { i1++; continue; } | 
|---|
| 925 | if (!pcheck2.Check(paq2)) { i2++; continue; } | 
|---|
| 926 | nokpaq++; | 
|---|
| 927 | if (paq1.FrameCounter()<paq2.FrameCounter()) { i1++; continue; } | 
|---|
| 928 | if (paq2.FrameCounter()<paq1.FrameCounter()) { i2++; continue; } | 
|---|
| 929 | //        cout << " DBG["<<kmz<<"]OKOK i1,i2=" << i1 <<","<<i2<<" FC1,FC2=" <<paq1.FrameCounter() | 
|---|
| 930 | // <<","<<paq2.FrameCounter()<<endl; | 
|---|
| 931 |  | 
|---|
| 932 | if ((i1>=memgr1.NbPaquets())||(i2>=memgr1.NbPaquets())) { | 
|---|
| 933 | cout << " *BUG*BUG i1=" << i1 << " i2=" << i2 << endl; | 
|---|
| 934 | break; | 
|---|
| 935 | } | 
|---|
| 936 | // Les deux framecounters sont identiques ... | 
|---|
| 937 | noksfc++; | 
|---|
| 938 | curfc=paq1.FrameCounter(); | 
|---|
| 939 | curfc2=paq2.FrameCounter(); | 
|---|
| 940 | if (fgfirst) { | 
|---|
| 941 | firsttt=paq1.TimeTag();  firsttt2=paq2.TimeTag(); | 
|---|
| 942 | cout << " BRProcB4C()/Info First FC="<<curfc<<" , "<<curfc2<<" -> TT=" | 
|---|
| 943 | << firsttt<<" , "<<firsttt2 <<endl; | 
|---|
| 944 | fgfirst=false; | 
|---|
| 945 | } | 
|---|
| 946 | curtt=paq1.TimeTag()-firsttt; | 
|---|
| 947 | curtt2=paq2.TimeTag()-firsttt2; | 
|---|
| 948 | dtr[0]=curfc;  dtr[1]=curtt; | 
|---|
| 949 | dtr[2]=curfc2;  dtr[3]=curtt2; | 
|---|
| 950 | dt.AddRow(dtr); | 
|---|
| 951 |  | 
|---|
| 952 | complex<r_4>* zp1 = (complex<r_4>*)(procbuff1+i1*procpaqsz); | 
|---|
| 953 | complex<r_4>* zp2 = (complex<r_4>*)(procbuff1+i1*procpaqsz+procpaqsz/2); | 
|---|
| 954 | complex<r_4>* zp3 = (complex<r_4>*)(procbuff2+i2*procpaqsz); | 
|---|
| 955 | complex<r_4>* zp4 = (complex<r_4>*)(procbuff2+i2*procpaqsz+procpaqsz/2); | 
|---|
| 956 | for(sa_size_t j=0; j<visiV13.Size(); j++) { | 
|---|
| 957 | visiV13(j)+=zp1[j]*conj(zp3[j]); | 
|---|
| 958 | visiV14(j)+=zp1[j]*conj(zp4[j]); | 
|---|
| 959 | visiV23(j)+=zp2[j]*conj(zp3[j]); | 
|---|
| 960 | visiV24(j)+=zp2[j]*conj(zp4[j]); | 
|---|
| 961 | } | 
|---|
| 962 | if (nzm==0) { | 
|---|
| 963 | visiV13.Info()["StartFC"] = curfc; | 
|---|
| 964 | visiV14.Info()["StartFC"] = curfc; | 
|---|
| 965 | visiV23.Info()["StartFC"] = curfc; | 
|---|
| 966 | visiV24.Info()["StartFC"] = curfc; | 
|---|
| 967 | visiV13.Info()["StartTT"] = curtt; | 
|---|
| 968 | visiV14.Info()["StartTT"] = curtt; | 
|---|
| 969 | visiV23.Info()["StartTT"] = curtt; | 
|---|
| 970 | visiV24.Info()["StartTT"] = curtt; | 
|---|
| 971 | } | 
|---|
| 972 | nzm++;  i1++; i2++; | 
|---|
| 973 | totnbytesproc += 2*paq1.DataSize(); | 
|---|
| 974 | } // Fin de boucle sur les paquets d'une zone | 
|---|
| 975 | memgr1.FreeMemZone(mid1, MemZS_ProcB); | 
|---|
| 976 | memgr2.FreeMemZone(mid2, MemZS_ProcB); | 
|---|
| 977 |  | 
|---|
| 978 | if ((nzm >= nmean_) || ((kmz==(nmax_-1))&&(nzm>1))) { | 
|---|
| 979 | visiV13 /= complex<r_4>((r_4)nzm, 0.); | 
|---|
| 980 | visiV14 /= complex<r_4>((r_4)nzm, 0.); | 
|---|
| 981 | visiV23 /= complex<r_4>((r_4)nzm, 0.); | 
|---|
| 982 | visiV24 /= complex<r_4>((r_4)nzm, 0.); | 
|---|
| 983 | visiV13.Info()["NPaqMoy"] = nzm; | 
|---|
| 984 | visiV14.Info()["NPaqMoy"] = nzm; | 
|---|
| 985 | visiV23.Info()["NPaqMoy"] = nzm; | 
|---|
| 986 | visiV24.Info()["NPaqMoy"] = nzm; | 
|---|
| 987 | visiV13.Info()["EndFC"] = curfc; | 
|---|
| 988 | visiV14.Info()["EndFC"] = curfc; | 
|---|
| 989 | visiV23.Info()["EndFC"] = curfc; | 
|---|
| 990 | visiV24.Info()["EndFC"] = curfc; | 
|---|
| 991 | visiV13.Info()["EndTT"] = curtt; | 
|---|
| 992 | visiV14.Info()["EndTT"] = curtt; | 
|---|
| 993 | visiV23.Info()["EndTT"] = curtt; | 
|---|
| 994 | visiV24.Info()["EndTT"] = curtt; | 
|---|
| 995 | char fname[512]; | 
|---|
| 996 | { | 
|---|
| 997 | sprintf(fname,"%s_%d.ppf",path_.c_str(),(int)ifile); | 
|---|
| 998 | POutPersist po(fname); | 
|---|
| 999 | po << PPFNameTag("visiV13") << visiV13; | 
|---|
| 1000 | po << PPFNameTag("visiV14") << visiV14; | 
|---|
| 1001 | po << PPFNameTag("visiV23") << visiV23; | 
|---|
| 1002 | po << PPFNameTag("visiV24") << visiV24; | 
|---|
| 1003 | } | 
|---|
| 1004 | visiV13 = complex<r_4>(0., 0.); | 
|---|
| 1005 | visiV14 = complex<r_4>(0., 0.); | 
|---|
| 1006 | visiV23 = complex<r_4>(0., 0.); | 
|---|
| 1007 | visiV24 = complex<r_4>(0., 0.); | 
|---|
| 1008 | nzm = 0;  ifile++; | 
|---|
| 1009 | //        ts.SetNow(); | 
|---|
| 1010 | //        filog << ts << " :  proc file  " << fname << endl; | 
|---|
| 1011 | cout << " BRProcB4C::run() created file  " << fname << endl; | 
|---|
| 1012 | } | 
|---|
| 1013 | double okfrac = (nokpaq>1)?((double)noksfc/(double)nokpaq*100.):0.; | 
|---|
| 1014 | cout << "BRProcB4C ["<<kmz<<"] NOKPaq=" << nokpaq << " NSameFC=" << noksfc | 
|---|
| 1015 | << " (" << okfrac << " %)" << endl; | 
|---|
| 1016 | totnokpaq += nokpaq; | 
|---|
| 1017 | totnoksfc += noksfc; | 
|---|
| 1018 | }  // Fin de boucle sur les zones a traiter | 
|---|
| 1019 | cout << " ------------------  BRProcB4C::run() END ----------------- " << endl; | 
|---|
| 1020 | { | 
|---|
| 1021 | dt.Info()["FirstTT1"]=firsttt; | 
|---|
| 1022 | dt.Info()["FirstTT2"]=firsttt2; | 
|---|
| 1023 | cout << dt; | 
|---|
| 1024 | char fname[512]; | 
|---|
| 1025 | sprintf(fname,"%s_fctt.ppf",path_.c_str()); | 
|---|
| 1026 | POutPersist po(fname); | 
|---|
| 1027 | po << PPFNameTag("ttfc") << dt; | 
|---|
| 1028 | cout << " BRProcB4C::run() created TimeTag/FrameCounter file " << fname  << endl; | 
|---|
| 1029 | } | 
|---|
| 1030 | ts.SetNow(); | 
|---|
| 1031 | tm.SplitQ(); | 
|---|
| 1032 | cout << "  TotalProc= " << totnbytesproc/(1024*1024) << " MBytes, rate= " | 
|---|
| 1033 | << (double)(totnbytesproc)/1024./tm.PartialElapsedTimems() << " MB/s" << endl; | 
|---|
| 1034 | double totokfrac = (totnokpaq>1)?((double)totnoksfc/(double)totnokpaq*100.):0.; | 
|---|
| 1035 | cout << " NOkPaq1,2=" << totnokpaq << " /TotNPaq=" << totnpaq << " TotNSameFC=" | 
|---|
| 1036 | << totnoksfc << " (" << totokfrac << " %)" << endl; | 
|---|
| 1037 | //  cout << pcheck1; | 
|---|
| 1038 | //  cout << pcheck2; | 
|---|
| 1039 | cout << " BRProcB4C::run()/Timing: \n"; | 
|---|
| 1040 | tm.Print(); | 
|---|
| 1041 | cout << " ---------------------------------------------------------- " << endl; | 
|---|
| 1042 | } | 
|---|
| 1043 | catch (PException& exc) { | 
|---|
| 1044 | cout << " BRProcB4C::run()/catched PException " << exc.Msg() << endl; | 
|---|
| 1045 | setRC(3); | 
|---|
| 1046 | return; | 
|---|
| 1047 | } | 
|---|
| 1048 | catch(...) { | 
|---|
| 1049 | cout << " BRProcB4C::run()/catched unknown ... exception " << endl; | 
|---|
| 1050 | setRC(4); | 
|---|
| 1051 | return; | 
|---|
| 1052 | } | 
|---|
| 1053 | setRC(0); | 
|---|
| 1054 | return; | 
|---|
| 1055 | } | 
|---|
| 1056 |  | 
|---|
| 1057 |  | 
|---|