Changeset 563 for BAORadio/AmasNancay
- Timestamp:
- Oct 10, 2011, 4:31:37 PM (13 years ago)
- Location:
- BAORadio/AmasNancay/trunk
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
BAORadio/AmasNancay/trunk/mergeAnaFiles.cc
r562 r563 36 36 37 37 //----------------------------------------------- 38 //Usage 39 // 40 // ./Objs/mergeAnaFiles -act meanCalibBAODiffOnOff -inPath /sps/baoradio/AmasNancay/JEC/ -src Abell85 -ppfFile dataRaw -debug 1 -mxcycles 500 41 // 42 // 43 // 44 //----------------------------------------------- 38 45 const sa_size_t NUMBER_OF_CHANNELS = 2; 39 46 const sa_size_t NUMBER_OF_FREQ = 8192; 40 47 const r_4 LOWER_FREQUENCY = 1250.0; //MHz 41 48 const r_4 TOTAL_BANDWIDTH = 250.0; //MHz 42 //-----------------------------------------------43 49 //Input parameters 44 50 struct Param { … … 68 74 sa_size_t chHigh, 69 75 TVector<r_4>& vec){ 76 70 77 sa_size_t nr = mtx.NRows(); 78 71 79 for (sa_size_t ir=0; ir<nr; ir++){ 72 80 TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false); … … 242 250 //------------------------------------------------------- 243 251 //Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra 244 //Used like:245 252 // 246 // void meanCalibBAODiffOnOffCycles() throw(string) {247 248 // list<string> listOfFiles;249 // string directoryName;250 // directoryName = para.inPath_ + "/" + para.sourceName_;251 252 // //Make the listing of the directory253 // listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);254 255 // list<string>::const_iterator iFile, iFileEnd, iSpecOff, iSpecOffEnd, iSpecOn, iSpecOnEnd;256 // iFileEnd = listOfFiles.end();257 258 259 // //Loop on files260 // uint_4 nRuns=0;261 // TArray<r_4> tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,para.maxNumberCycles_); //para.maxNumberCycles_ should be large enough...262 263 264 // for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {265 // if (para.debuglev_>90){266 // cout << "load file <" << *iFile << ">" << endl;267 // }268 269 // vector<string> tokens;270 // split(*File,"_",tokens);271 // string dateOfRun = tokens[1];272 // string srcLower = tokens[2];273 // if (para.debuglev_>90){274 // cout << "date <" << dateOfRun << ">" << endl;275 // }276 277 // PInPersist fin(*iFile);278 // vector<string> vec = fin.GetNameTags();279 280 // if (para.typeOfCalib_ == "perRun") {281 // ///////////////////282 // //make the calibration of the mean of all Off and On of the run and perform the difference283 284 // vector<string> modeList;285 // modeList.push_back("On");286 // modeList.push_back("Off");287 // vector<string>::const_iterator iMode;288 // map<string, TMatrix<r_4> > spectreModeCollect;289 290 // for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {291 // ///////////////////292 // //293 // //Compute the mean of the mode294 // //295 // list<string> listOfSpectra;296 // //Keep only required PPF objects297 // string matchstr = "specRaw"+(*iMode)+"[0-9]+";298 // std::remove_copy_if(299 // vec.begin(), vec.end(), back_inserter(listOfSpectra),300 // not1(StringMatch(matchstr))301 // );302 303 // listOfSpectra.sort(stringCompare);304 // iSpecEnd = listOfSpectra.end();305 // TMatrix<r_4> meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);306 // uint_4 nSpectra=0;307 // //Loop of spectra matrix308 // for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd; ++iSpec){309 // if (para.debuglev_>90){310 // cout << " spactra <" << *iSpec << ">" << endl;311 // }312 // TMatrix<r_4> aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);313 // fin.GetObject(aSpec,*iSpec);314 // //How to see if the GetObject is ok?? Ask Reza315 // nSpectra++;316 // meanOfSpectra+=aSpec;317 // }//end loop Off318 // //Mean319 // if(nSpectra>0)meanOfSpectra=(r_4)(nSpectra);320 321 // //BAO Calibrator322 // string calibFileName = directoryName + "/"323 // + "calib_" + dateOfRun + "_" + srcLower + "_"+(*iMode)+"_"324 // + para.calibFreq_ +"MHz-All.txt";325 // if(debuglev_>0) cout << "Read Calib file " << calibFileName << endl;326 // ifstream ifs(calibFileName.c_str());327 // if ( ! ifs.is_open() ) {328 // rc = 999;329 // throw calibFileName + " cannot be opened...";330 // }331 // TVector<r_4> calibBAOfactors;332 // sa_size_t nr,nc; //values read333 // calibBAOfactorsOff.ReadASCII(ifs,nr,nc);334 // if(debuglev_>9){335 // cout << "(nr,nc): "<< nr << "," << nc << endl;336 // calibBAOfactors.Print(cout);337 // }338 339 // for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){340 // meanOfSpectra(Range(iCh), Range::all()) /=calibBAOfactors(iCh);341 // }342 // //storage343 // spectreModeCollect.insert(pair<string, TMatrix<r_4> >(*iMode,TMatrix<r_4>(meanOfSpectra,false))); //do not share data (cf. SOPHYA)344 // }//end of mode345 346 // //Take the difference ON-OFF in current run347 // TMatrix<r_4>diffOnOff = spectreModeCollect["On"]-spectreModeCollect["Off"];348 349 350 // } else if (para.typeOfCalib_ == "perCycle") {351 // //perform the calibration of the OFF and ON per cycle, then make the mean and take the diff352 // } else {353 // string msg="FATAL (meanCalibBAODiffOnOffCycles); unknown calibration mode "354 // + para.typeOfCalib_ ;355 // throw(msg);356 // }357 358 359 // nRuns++;360 361 // }//eo loop on spectra in a file362 // }363 253 void meanCalibBAODiffOnOffCycles() throw(string) { 364 254 … … 378 268 TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //set to 0 379 269 char* onffTupleName[7]={"cycle", 380 "onoffRaw0","onoffRun0","onoffCycle0", 381 "onoffRaw1","onoffRun1","onoffCycle1"}; 270 "onoffRaw0","onoffRaw1", 271 "onoffRun0","onoffRun1", 272 "onoffCycle0","onoffCycle1"}; 382 273 NTuple onoffevolution(7,onffTupleName); 383 274 r_4 xnt[7]; 384 275 385 276 //Lower and Higher freq. bin to perform mean follow up 386 sa_size_t chLow = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5));277 sa_size_t chLow = freqToChan(para.rcalibFreq_ - (para.rcalibBandFreq_*0.5)); 387 278 sa_size_t chHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5)); 279 if (para.debuglev_>90){ 280 cout << "freq. band for follow up [" << chLow << ", "<< chHigh << "]" << endl; 281 } 388 282 389 283 //Loop on files/run … … 501 395 double mean,sigma; 502 396 MeanSigma(coef,mean,sigma); 503 cout << "Mean: " << mean << " sigma " << sigma << endl;504 397 calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean); 505 398 } … … 533 426 double mean,sigma; 534 427 MeanSigma(coef,mean,sigma); 535 cout << "Mean: " << mean << " sigma " << sigma << endl;428 // cout << "Mean: " << mean << " sigma " << sigma << endl; 536 429 calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean); 537 430 } … … 566 459 double mean,sigma; 567 460 MeanSigma(coef,mean,sigma); 568 cout << "Mean: " << mean << " sigma " << sigma << endl;461 // cout << "Mean: " << mean << " sigma " << sigma << endl; 569 462 calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean); 570 463 } … … 597 490 double mean,sigma; 598 491 MeanSigma(coef,mean,sigma); 599 cout << "Mean: " << mean << " sigma " << sigma << endl;492 // cout << "Mean: " << mean << " sigma " << sigma << endl; 600 493 calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean); 601 494 } … … 623 516 nr = calibBAOfactors_Off_Run_Ch0.NRows(); 624 517 for (sa_size_t ir=0; ir<nr; ++ir){ 518 // cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val " 519 // << calibBAOfactors_Off_Run_Ch0(ir,1) << endl; 520 625 521 calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)] 626 522 = calibBAOfactors_Off_Run_Ch0(ir,1); … … 631 527 calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)] 632 528 = calibBAOfactors_Off_Cycle_Ch1(ir,1); 633 } 529 }//eo loop on coef Off 634 530 635 531 nr = calibBAOfactors_On_Run_Ch0.NRows(); … … 643 539 calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)] 644 540 = calibBAOfactors_On_Cycle_Ch1(ir,1); 645 } 541 }//eo loop on coef On 646 542 647 648 649 // for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){650 // meanOfSpectra(Range(iCh), Range::all()) /=calibBAOfactors(iCh);651 // }652 653 543 //Loop on cyles 654 544 for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){ 655 545 656 if(para.debuglev_>9){ 546 //Look if the cycle has been calibrated... 547 bool isCycleCalibrated = 548 ( calibBAO_On_Run_Ch0.count(*ic) * 549 calibBAO_On_Run_Ch1.count(*ic) * 550 calibBAO_Off_Run_Ch0.count(*ic) * 551 calibBAO_Off_Run_Ch1.count(*ic) * 552 calibBAO_On_Cycle_Ch0.count(*ic) * 553 calibBAO_On_Cycle_Ch1.count(*ic) * 554 calibBAO_Off_Cycle_Ch0.count(*ic) * 555 calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false; 556 557 if(para.debuglev_>9) { 657 558 cout << "Calibration coefficients for cycle "<<*ic << endl; 658 cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " " 659 << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n" 660 << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " " 661 << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n" 662 << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " " 663 << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n" 664 << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " " 665 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl; 666 } 559 if (isCycleCalibrated) { 560 cout << "Cycle calibrated " << endl; 561 cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " " 562 << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n" 563 << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " " 564 << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n" 565 << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " " 566 << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n" 567 << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " " 568 << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl; 569 } else { 570 cout << "Cycle NOT calibrated " << endl; 571 } 572 }//debug 573 574 575 if ( ! isCycleCalibrated ) continue; 667 576 668 577 string ppftag; … … 706 615 meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib; 707 616 617 // 618 totalNumberCycles++; 708 619 //Fill NTuple 709 xnt[0] = *ic;620 xnt[0] = totalNumberCycles; 710 621 711 622 TVector<r_4> meanInRange_noCalib(NUMBER_OF_CHANNELS); … … 721 632 TVector<r_4> meanInRange_perCycleCalib(NUMBER_OF_CHANNELS); 722 633 meanInRange(diffOnOff_perCycleCalib,chLow,chHigh,meanInRange_perCycleCalib); 723 xnt[ 3] = meanInRange_perCycleCalib(0);724 xnt[ 4] = meanInRange_perCycleCalib(1);634 xnt[5] = meanInRange_perCycleCalib(0); 635 xnt[6] = meanInRange_perCycleCalib(1); 725 636 726 637 onoffevolution.Fill(xnt); 727 638 728 totalNumberCycles++;639 //Quit if enough 729 640 if (totalNumberCycles >= para.maxNumberCycles_) break; 730 641 … … 735 646 cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl; 736 647 //Normalisation 737 if(totalNumberCycles >0){648 if(totalNumberCycles > 0){ 738 649 meanDiffONOFF_noCalib /= (r_4)totalNumberCycles; 739 650 meanDiffONOFF_perRunCalib /= (r_4)totalNumberCycles; … … 758 669 tmp << totalNumberCycles; 759 670 string fileName = para.outPath_+"/onoffsurvey_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf"; 671 760 672 POutPersist fos(fileName); 673 cout << "Save output in " << fileName << endl; 761 674 762 675 string tag = "meanNoCalib"; -
BAORadio/AmasNancay/trunk/plotDiffOnOff.pic
r556 r563 8 8 set action $4 9 9 10 set defatt "font=helvetica,bold,20 "10 set defatt "font=helvetica,bold,20 fixedfontsize" 11 11 12 12
Note: See TracChangeset
for help on using the changeset viewer.