Changeset 553 for BAORadio/AmasNancay/trunk/mergeRawOnOff.cc
- Timestamp:
- Oct 6, 2011, 7:32:54 AM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
BAORadio/AmasNancay/trunk/mergeRawOnOff.cc
r552 r553 40 40 const sa_size_t NUMBER_OF_CHANNELS = 2; 41 41 const sa_size_t NUMBER_OF_FREQ = 8192; 42 const sa_size_t TOTAL_NUM_CYCLES = 10000;42 const sa_size_t TOTAL_NUM_CYCLES = 500; 43 43 const r_4 LOWER_FREQUENCY = 1250.0; //MHz 44 44 const r_4 TOTAL_BANDWIDTH = 250.0; //MHz … … 56 56 //Utility functions 57 57 58 sa_size_t freqTochan(r_4 f){ 59 return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ); 60 } 61 //--------- 58 62 class median_of_empty_list_exception:public std::exception{ 59 63 virtual const char* what() const throw() { … … 309 313 listOfSpectra.sort(stringCompare); 310 314 //Loop of spectra matrix 311 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd ; ++iSpec){315 for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < TOTAL_NUM_CYCLES ; ++iSpec){ 312 316 if (para.debuglev_>90){ 313 317 cout << " spactra <" << *iSpec << ">" << endl; … … 316 320 fin.GetObject(aSpec,*iSpec); 317 321 318 if ((sa_size_t)nSpectra < TOTAL_NUM_CYCLES ) { 319 tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec; 320 } else { 321 stringstream tmp; 322 tmp<<nSpectra; 323 throw "FATAL: SIZE ERROR: too many cycles:"+tmp.str(); 324 } 322 tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec; 323 325 324 nSpectra++; 326 325 }//eo loop on spectra in a file 327 326 }//eo loop on files 328 327 329 //Resize330 328 331 329 … … 347 345 reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx); 348 346 347 348 sa_size_t f1320=freqTochan(1320.); 349 sa_size_t f1345=freqTochan(1345.); 350 sa_size_t f1355=freqTochan(1355.); 351 sa_size_t f1380=freqTochan(1380.); 352 //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380] 353 if (para.debuglev_>9){ 354 cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl; 355 } 356 TVector<r_4>meanMed(NUMBER_OF_CHANNELS); 357 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){ 358 double meanMed1; 359 double sigmaMed1; 360 TVector<r_4> band1; 361 band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions(); 362 MeanSigma(band1,meanMed1,sigmaMed1); 363 double meanMed2; 364 double sigmaMed2; 365 TVector<r_4> band2; 366 band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions(); 367 MeanSigma(band2,meanMed2,sigmaMed2); 368 meanMed(iCh) = (meanMed1+meanMed2)*0.5; 369 } 370 meanMed.Print(cout); 371 372 //Compute the sigma in the range 1320MHz-1380MHz 373 if (para.debuglev_>9){ 374 cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl; 375 } 376 TVector<r_4>sigmaMed(NUMBER_OF_CHANNELS); 377 sa_size_t redf1320=(1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_; 378 sa_size_t redf1380=(1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_; 379 380 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){ 381 double meanSigma; 382 double sigmaSigma; 383 TVector<r_4> band; 384 band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions(); 385 MeanSigma(band,meanSigma,sigmaSigma); 386 meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra 387 sigmaMed(iCh) = meanSigma; 388 } 389 sigmaMed.Print(cout); 390 391 392 if (para.debuglev_>9){ 393 cout << "Compute medianOfSpectraNorm" << endl; 394 } 395 TMatrix<r_4> medianOfSpectraNorm(medianOfSpectra,false); //do not share the data... 396 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){ 397 medianOfSpectraNorm.Row(iCh) -= meanMed(iCh); 398 medianOfSpectraNorm.Row(iCh) /= sigmaMed(iCh); 399 } 400 401 402 349 403 {//Save the result 350 404 stringstream tmp; … … 356 410 string tag = "median"; 357 411 fos << PPFNameTag(tag) << medianOfSpectra; 412 413 tag = "medianNorm"; 414 fos << PPFNameTag(tag) << medianOfSpectraNorm; 415 416 358 417 tag = "meanmedred"; 359 418 fos << PPFNameTag(tag) << meanRedMtx;
Note: See TracChangeset
for help on using the changeset viewer.