Changeset 524 for BAORadio/AmasNancay
- Timestamp:
- Sep 23, 2011, 10:18:17 AM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
BAORadio/AmasNancay/analyse.cc
r523 r524 10 10 #include <stdlib.h> 11 11 #include <dirent.h> 12 12 #include <matharr.h> 13 13 14 14 // include standard c/c++ … … 202 202 virtual int processCmd() throw(string); 203 203 }; 204 205 //JEC 22/9/11 Make ON-OFF analysis WO any calibration START 206 //------------ 207 //Process ON/OFF Raw data 208 //------------ 209 class ProcessONOFFRawData : public ProcessBase { 210 211 public: 212 ProcessONOFFRawData(){} 213 virtual ~ProcessONOFFRawData(){} 214 215 virtual int processCmd() throw(string); 216 }; 217 //JEC 22/9/11 Make ON-OFF analysis WO any calibration END 218 204 219 //------------ 205 220 //Process Gain … … 263 278 //Init process types 264 279 map<string,IProcess*> process; 280 //JEC 22/9/11 Make ON-OFF analysis WO any calibration START 281 process["rawOnOff"] = new ProcessONOFFRawData(); 282 //JEC 22/9/11 Make ON-OFF analysis WO any calibration END 265 283 process["dataOnOff"] = new ProcessONOFFData(); 266 284 process["gain"] = new ProcessGain(); … … 363 381 364 382 //JEC 21/9/11 Give the input parameters START 365 cout << " action " << action << "\n" 383 cout << "Dump Iiitial parameters ............" << endl; 384 cout << " action = " << action << "\n" 366 385 << " inputPath = " << inputPath << "\n" 367 386 << " outputPath = " << outputPath << "\n" … … 376 395 << " numcycle = " << numcycle << "\n" 377 396 << " calibrationOpt = " << calibrationOpt << endl; 397 cout << "...................................." << endl; 378 398 //JEC 21/9/11 Give the input parameters END 379 399 … … 388 408 } 389 409 390 391 392 393 410 394 411 // … … 450 467 throw e.what(); 451 468 } 469 470 //JEC 22/9/11 Make ON-OFF analysis WO any calibration START 471 // try { 472 // ProcessONOFFRawData* procRawdata = dynamic_cast<ProcessONOFFRawData*>(process[action]); 473 // } 474 // catch(exception& e){ 475 // throw e.what(); 476 // } 477 //JEC 22/9/11 Make ON-OFF analysis WO any calibration END 478 452 479 453 480 try { … … 543 570 int Usage(bool flag) { 544 571 cout << "Analyse.cc usage...." << endl; 545 cout << "analyse -act <action_type>: dataOnOff, gain, calib\n"572 cout << "analyse -act <action_type>: dataOnOff, rawOnOff, gain, calib\n" 546 573 << " -inPath <path for input files: default='.'>\n" 547 574 << " -outPath <path for output files: default='.'>\n" … … 620 647 return rc; 621 648 } 649 //JEC 22/9/11 Make ON-OFF analysis WO any calibration START 622 650 //---------------------------------------------- 623 int ProcessONOFF Data::processCmd() throw(string) {651 int ProcessONOFFRawData::processCmd() throw(string) { 624 652 int rc = 0; 625 653 try { … … 629 657 throw s; 630 658 } 631 if(debuglev_>0)cout << "Process Data" << endl;659 if(debuglev_>0)cout << "Process Raw Data ON OFF" << endl; 632 660 vector<string> modeList; 633 661 modeList.push_back("On"); … … 641 669 //Process to get sucessively 642 670 //Raw Spectra, 643 //BAO Calibrated Spectra644 //and RT Calibrated Spectra645 671 //The pocesses are separated to allow intermediate save of results 646 672 … … 665 691 // 666 692 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) { 667 //TOBE FIXED directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/";668 693 directoryName = "./" + mode + "/"; 669 694 stringstream sicycle; … … 697 722 cout << "Save mean raw spectra" << endl; 698 723 string fileName; 724 fileName = "./dataRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf"; 725 726 POutPersist fos(fileName); 727 id=0; 728 iSpectreEnd = spectreCollect.end(); 729 for (iSpectre = spectreCollect.begin(); 730 iSpectre != iSpectreEnd ; ++iSpectre, ++id) { 731 tag = "specRaw"; 732 tag += (iSpectre->first).first; 733 stringstream sid; 734 sid << (iSpectre->first).second; 735 tag += sid.str(); 736 if(debuglev_>9) { 737 cout << "save tag<" << tag << ">" << endl; 738 } 739 fos << PPFNameTag(tag) << iSpectre->second; 740 } 741 }//end of save fits 742 743 744 //------------------------------------------ 745 // Perform ON-OFF 746 //------------------------------------------ 747 748 map< sa_size_t, TMatrix<r_4> > diffCollect; 749 map< sa_size_t, TMatrix<r_4> >::iterator iDiff, iDiffEnd; 750 751 TMatrix<r_4> diffMeanOnOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //init zero 752 r_4 nCycles = 0; 753 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) { 754 nCycles++; 755 TMatrix<r_4> specmtxOn(spectreCollect[make_pair(modeList[0],icycle)],false); //clone the memory 756 TMatrix<r_4> specmtxOff(spectreCollect[make_pair(modeList[1],icycle)],false); //clone the memory 757 TMatrix<r_4> diffOnOff = specmtxOn - specmtxOff; 758 diffCollect.insert(pair< sa_size_t,TMatrix<r_4> >(icycle,TMatrix<r_4>(diffOnOff,false) )); 759 diffMeanOnOff += diffOnOff; 760 }//end loop on cycle 761 if(nCycles>0) diffMeanOnOff/=nCycles; 762 763 //exctract channels and do the mean 764 TVector<r_4> meanOfChan(NUMBER_OF_FREQ); //implicitly init to 0 765 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh) { 766 meanOfChan += diffMeanOnOff.Row(iCh).Transpose(); 767 } 768 meanOfChan /= (r_4)NUMBER_OF_CHANNELS; 769 770 771 772 {//save diff ON-OFF on Raw data 773 if(debuglev_>0)cout << "save ON-OFF RAW spectra" << endl; 774 string fileName; 775 fileName = "./diffOnOffRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf"; 776 POutPersist fos(fileName); 777 iDiffEnd = diffCollect.end(); 778 id = 0; 779 780 //JEC 22/9/11 Mean & Sigma in 32-bins size START 781 sa_size_t nSliceFreq = 32; //TODO: put as an input parameter option ? 782 sa_size_t deltaFreq = NUMBER_OF_FREQ/nSliceFreq; 783 //JEC 22/9/11 Mean & Sigma in 32-bins size END 784 785 for (iDiff = diffCollect.begin();iDiff != iDiffEnd ; ++iDiff, id++) { 786 tag = "specONOFFRaw"; 787 stringstream sid; 788 sid << iDiff->first; 789 tag += sid.str(); 790 fos << PPFNameTag(tag) << iDiff->second; 791 792 //JEC 22/9/11 Mean & Sigma in 32-bins size START 793 if (debuglev_>9) { 794 cout << "Debug slicing: slice/delta " << nSliceFreq << " " << deltaFreq << endl; 795 } 796 TMatrix<r_4> reducedMeanDiffOnOff(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default 797 TMatrix<r_4> reducedSigmaDiffOnOff(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default 798 for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){ 799 sa_size_t freqLow= iSlice*deltaFreq; 800 sa_size_t freqHigh= freqLow + deltaFreq -1; 801 if (debuglev_>9) { 802 cout << "Debug .......... slicing ["<< iSlice << "]: low/high freq" << freqLow << "/" << freqHigh << endl; 803 } 804 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){ 805 if (debuglev_>9) { 806 cout << "Debug .......... Channel " << iCh; 807 } 808 TVector<r_4> reducedRow; 809 reducedRow = (iDiff->second).SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions(); 810 double mean; 811 double sigma; 812 MeanSigma(reducedRow,mean,sigma); 813 if (debuglev_>9) { 814 cout << "mean/signa " << mean << "/" << sigma << endl; 815 } 816 reducedMeanDiffOnOff(iCh,iSlice) = mean; 817 reducedSigmaDiffOnOff(iCh,iSlice) = sigma; 818 }//loop on Channel 819 }//loop on Freq. slice 820 tag = "redMeanONOFFRaw"; 821 tag += sid.str(); 822 fos << PPFNameTag(tag) << reducedMeanDiffOnOff; 823 tag = "redSigmaONOFFRaw"; 824 tag += sid.str(); 825 fos << PPFNameTag(tag) << reducedSigmaDiffOnOff; 826 //JEC 22/9/11 END 827 828 }//loop on ON-OFF spectre 829 //save the mean also 830 fos << PPFNameTag("specONOFFRawMean") << diffMeanOnOff; 831 832 //JEC 22/9/11 START 833 TMatrix<r_4> reducedMeanDiffOnOffAll(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default 834 TMatrix<r_4> reducedSigmaDiffOnOffAll(NUMBER_OF_CHANNELS,nSliceFreq); //init 0 by default 835 for (sa_size_t iSlice=0; iSlice<nSliceFreq; iSlice++){ 836 sa_size_t freqLow= iSlice*deltaFreq; 837 sa_size_t freqHigh= freqLow + deltaFreq -1; 838 if (debuglev_>9) { 839 cout << "Debug .......... slicing ["<< iSlice << "]: low/high freq" << freqLow << "/" << freqHigh << endl; 840 } 841 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; ++iCh){ 842 if (debuglev_>9) { 843 cout << "Debug .......... Channel " << iCh; 844 } 845 TVector<r_4> reducedRow; 846 reducedRow = diffMeanOnOff.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions(); 847 double mean; 848 double sigma; 849 MeanSigma(reducedRow,mean,sigma); 850 if (debuglev_>9) { 851 cout << "mean/signa " << mean << "/" << sigma << endl; 852 } 853 reducedMeanDiffOnOffAll(iCh,iSlice) = mean; 854 reducedSigmaDiffOnOffAll(iCh,iSlice) = sigma; 855 }//loop on Channel 856 }//loop on Freq. slice 857 tag = "redMeanONOFFRawAll"; 858 fos << PPFNameTag(tag) << reducedMeanDiffOnOffAll; 859 tag = "redSigmaONOFFRawAll"; 860 fos << PPFNameTag(tag) << reducedSigmaDiffOnOffAll; 861 //JEC 22/9/11 END 862 863 864 865 fos << PPFNameTag("specONOFFRaw2ChanMean") << meanOfChan; 866 }//end of save fits 867 868 return rc; 869 } //ProcessONOFFRawData::processCmd 870 871 //JEC 22/9/11 Make ON-OFF analysis WO any calibration END 872 //---------------------------------------------- 873 int ProcessONOFFData::processCmd() throw(string) { 874 int rc = 0; 875 try { 876 rc = ProcessBase::processCmd(); 877 } 878 catch (string s) { 879 throw s; 880 } 881 if(debuglev_>0)cout << "Process Data" << endl; 882 vector<string> modeList; 883 modeList.push_back("On"); 884 modeList.push_back("Off"); 885 vector<string>::const_iterator iMode; 886 887 uint_4 id; 888 string tag; 889 890 // 891 //Process to get sucessively 892 //Raw Spectra, 893 //BAO Calibrated Spectra 894 //and RT Calibrated Spectra 895 //The pocesses are separated to allow intermediate save of results 896 897 map< pair<string, sa_size_t>, TMatrix<r_4> > spectreCollect; 898 map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iSpectre, iSpectreEnd; 899 900 for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) { 901 string mode = *iMode; 902 if(debuglev_>0)cout << "Process RAW Mode " << mode << endl; 903 904 //------------------------------------------ 905 //Produce Raw spectra per cycle 906 //------------------------------------------ 907 908 string directoryName; 909 list<string> listOfSpecFiles; 910 list<string>::const_iterator iFile, iFileEnd; 911 912 913 // 914 //loop on cycles 915 // 916 for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) { 917 //TOBE FIXED directoryName = "./" + sourceName_ + "/"+ dateOfRun_ + StringToLower(sourceName_) + "/" +mode + "/"; 918 directoryName = "./" + mode + "/"; 919 stringstream sicycle; 920 sicycle << icycle; 921 directoryName += spectraDirectory_ + sicycle.str() + "/"; 922 923 //read directory 924 listOfSpecFiles = ListOfFileInDir(directoryName,typeOfFile_); 925 926 927 //compute mean of spectra created in a cycle 928 if(debuglev_>0)cout << "compute mean for cycle " << icycle << endl; 929 TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0 930 iFileEnd = listOfSpecFiles.end(); 931 r_4 nSpectres = 0; 932 for (iFile = listOfSpecFiles.begin(); iFile != iFileEnd; ++iFile) { 933 FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO); 934 TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); 935 aSpectrum >> spectre; 936 spectreMean += spectre; 937 nSpectres++; 938 }// end of for files 939 if(nSpectres>0) spectreMean /= nSpectres; 940 941 //save mean spectrum 942 spectreCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreMean,false) )); 943 }//end of for cycles 944 }//end loop on mode for raw preocess 945 946 if(debuglev_>1) {//save mean spectra on file 947 cout << "Save mean raw spectra" << endl; 948 string fileName; 699 949 //TOBE FIXED fileName = "./" + sourceName_ + "/" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_" + "dataRaw" + ".ppf"; 700 950 fileName = "./dataRaw_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf"; … … 734 984 //Read BAO calibration files 735 985 sa_size_t nr,nc; //values read 736 //JEC 20/9/11 use mean calibration coeff upon all cycles START737 // bool first = true;738 // TMatrix<r_4> calibBAOfactors;739 // for (sa_size_t iCh=0;iCh<NUMBER_OF_CHANNELS;++iCh){740 // string calibFileName = inputPath_+ "/"741 // + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_)742 // + "/calib_" + dateOfRun_ + "_" + StringToLower(sourceName_) + "_"743 // + mode + "_" + freqBAOCalibration_ + "MHz";744 // stringstream channel;745 // channel << iCh;746 // calibFileName += "-Ch" + channel.str() + "Cycles.txt";747 // if(debuglev_>0) cout << "Read file " << calibFileName << endl;748 // ifstream ifs(calibFileName.c_str());749 // if ( ! ifs.is_open() ) {750 // rc = 999;751 // throw calibFileName + " cannot be opened...";752 // }753 // TVector<r_4> aCalib;754 // if(debuglev_>9) cout << "Debug 1" << endl;755 // aCalib.ReadASCII(ifs,nr,nc);756 // if(debuglev_>9) cout << "Debug 2" << endl;757 // if(first){758 // first = false;759 // if(debuglev_>9) cout << "Debug 3" << endl;760 // calibBAOfactors.SetSize(NUMBER_OF_CHANNELS,nr);761 // if(debuglev_>9) cout << "Debug 4" << endl;762 // }763 // calibBAOfactors( Range(iCh), Range::all() ) = aCalib.Transpose();764 // if(debuglev_>9) cout << "Debug 5" << endl;765 // ifs.close();766 // }//end of loop on channels767 768 986 769 987 string calibFileName = inputPath_+ "/"
Note: See TracChangeset
for help on using the changeset viewer.