// Utilisation de SOPHYA pour faciliter les tests ... #include //#include #include #include "sopnamsp.h" #include "machdefs.h" #include #include #include // include standard c/c++ #include #include #include #include #include #include #include #include #include #include // include sophya mesure ressource CPU/memoire ... #include "resusage.h" #include "ctimer.h" #include "timing.h" #include "timestamp.h" #include "strutilxx.h" #include "ntuple.h" #include "fioarr.h" #include "tarrinit.h" #include "histinit.h" #include "fitsioserver.h" #include "fiosinit.h" #include "ppersist.h" //----------------------------------------------- const sa_size_t NUMBER_OF_CHANNELS = 2; const sa_size_t NUMBER_OF_FREQ = 8192; const sa_size_t TOTAL_NUM_CYCLES = 500; const r_4 LOWER_FREQUENCY = 1250.0; //MHz const r_4 TOTAL_BANDWIDTH = 250.0; //MHz //----------------------------------------------- //Input parameters struct Param { int debuglev_; //debug string inPath_; //root directory of the input files string outPath_; //output files are located here string sourceName_; //source name & subdirectory of the input files string ppfFile_; //generic name of the input files int nSliceInFreq_; //used by reduceSpectra() fnc string typeOfCalib_;//type of calibration to be done } para; //-------------------------------------------------------------- //Utility functions sa_size_t freqTochan(r_4 f){ return (sa_size_t)((f-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*NUMBER_OF_FREQ); } //--------- class median_of_empty_list_exception:public std::exception{ virtual const char* what() const throw() { return "Attempt to take the median of an empty list of numbers. " "The median of an empty list is undefined."; } }; template double median(RandAccessIter begin, RandAccessIter end) throw(median_of_empty_list_exception){ if(begin == end){ throw median_of_empty_list_exception(); } std::size_t size = end - begin; std::size_t middleIdx = size/2; RandAccessIter target = begin + middleIdx; std::nth_element(begin, target, end); if(size % 2 != 0){ //Odd number of elements return *target; }else{ //Even number of elements double a = *target; RandAccessIter targetNeighbor= target-1; std::nth_element(begin, targetNeighbor, end); return (a+*targetNeighbor)/2.0; } } //-------------------------------------------------------------- char *regexp (const char *string, const char *patrn, int *begin, int *end) { int i, w=0, len; char *word = NULL; regex_t rgT; regmatch_t match; regcomp(&rgT,patrn,REG_EXTENDED); if ((regexec(&rgT,string,1,&match,0)) == 0) { *begin = (int)match.rm_so; *end = (int)match.rm_eo; len = *end-*begin; word=(char*)malloc(len+1); for (i=*begin; i<*end; i++) { word[w] = string[i]; w++; } word[w]=0; } regfree(&rgT); return word; } //------- sa_size_t round_sa(r_4 r) { return static_cast((r > 0.0) ? (r + 0.5) : (r - 0.5)); } //----- string StringToLower(string strToConvert){ //change each element of the string to lower case for(unsigned int i=0;i tolower( *rit ) ) return false; return false; } //----- list ListOfFileInDir(string dir, string filePettern) throw(string) { list theList; DIR *dip; struct dirent *dit; string msg; string fileName; string fullFileName; size_t found; if ((dip=opendir(dir.c_str())) == NULL ) { msg = "opendir failed on directory "+dir; throw msg; } while ( (dit = readdir(dip)) != NULL ) { fileName = dit->d_name; found=fileName.find(filePettern); if (found != string::npos) { fullFileName = dir + "/"; fullFileName += fileName; theList.push_back(fullFileName); } }//eo while if (closedir(dip) == -1) { msg = "closedir failed on directory "+dir; throw msg; } theList.sort(stringCompare); return theList; } // class StringMatch : public unary_function { public: StringMatch(const string& pattern): pattern_(pattern){} bool operator()(const string& aStr) const { int b,e; regexp(aStr.c_str(),pattern_.c_str(),&b,&e); // cout << "investigate " << aStr << " to find " << pattern_ // << "[" <& specMtxInPut, TMatrix& meanMtx, TMatrix& sigmaMtx) { sa_size_t nSliceFreq = para.nSliceInFreq_; sa_size_t deltaFreq = NUMBER_OF_FREQ/nSliceFreq; for (sa_size_t iSlice=0; iSlice reducedRow; reducedRow = specMtxInPut.SubMatrix(Range(iCh),Range(freqLow,freqHigh)).CompactAllDimensions(); double mean; double sigma; MeanSigma(reducedRow,mean,sigma); meanMtx(iCh,iSlice) = mean; sigmaMtx(iCh,iSlice) = sigma; }//eo loop on channels }//eo loop on slices } //------------------------------------------------------- //Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra //Used like: // void meanCalibBAODiffOnOffCycles() throw(string) { list listOfFiles; string directoryName; directoryName = para.inPath_ + "/" + para.sourceName_; //Make the listing of the directory listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_); list::const_iterator iFile, iFileEnd, iSpecOff, iSpecOffEnd, iSpecOn, iSpecOnEnd; iFileEnd = listOfFiles.end(); //Loop on files for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) { if (para.debuglev_>90){ cout << "load file <" << *iFile << ">" << endl; } PInPersist fin(*iFile); vector vec = fin.GetNameTags(); if (para.typeOfCalib_ == "perRun") { /////////////////// //make the calibration of the mean of all Off and On of the run and perform the difference /////////////////// // //Compute the mean of the Off // list listOfSpectraOff; //Keep only required PPF objects std::remove_copy_if( vec.begin(), vec.end(), back_inserter(listOfSpectra), not1(match) ); listOfSpectraOff.sort(stringCompare); iSpecOffEnd = listOfSpectraOff.end(); TMatrix meanOfSpectraOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); uint_4 nSpectraOff=0; //Loop of spectra matrix for (iSpecOff = listOfSpectraOff.begin(); iSpecOff !=iSpecOffEnd; ++iSpecOff){ if (para.debuglev_>90){ cout << " spactra <" << *iSpecOff << ">" << endl; } TMatrix aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); fin.GetObject(aSpec,*iSpecOff); //How to see if the GetObject is ok?? Ask Reza nSpectraOff++; meanOfSpectraOff+=aSpec; }//end loop Off //Normalisation if(nSpectraOff>0)meanOfSpectraOff/=(r_4)(nSpectraOff); // //Compute the mean of the On // list listOfSpectraOn; //Keep only required PPF objects std::remove_copy_if( vec.begin(), vec.end(), back_inserter(listOfSpectra), not1(match) ); listOfSpectraOn.sort(stringCompare); iSpecOnEnd = listOfSpectraOn.end(); TMatrix meanOfSpectraOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); uint_4 nSpectraOn=0; //Loop of spectra matrix for (iSpecOn = listOfSpectraOn.begin(); iSpecOn !=iSpecOnEnd; ++iSpecOn){ if (para.debuglev_>90){ cout << " spactra <" << *iSpecOn << ">" << endl; } TMatrix aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); fin.GetObject(aSpec,*iSpecOn); //How to see if the GetObject is ok?? Ask Reza nSpectraOn++; meanOfSpectraOn+=aSpec; }//end loop On //Normalisation if(nSpectraOn>0)meanOfSpectraOn/=(r_4)(nSpectraOn); // //Get the calibration file // } else if {para.typeOfCalib_ == "perCycle") { //perform the calibration of the OFF and ON per cycle, then make the mean and take the diff } else { string msg="FATAL (meanCalibBAODiffOnOffCycles); unknown calibration mode " + para.typeOfCalib_ ; throw(msg); } }//eo loop on spectra in a file }//eo loop on files } //------------------------------------------------------- //Compute the mean of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra //Used like: // void meanRawDiffOnOffCycles() throw(string) { list listOfFiles; string directoryName; directoryName = para.inPath_ + "/" + para.sourceName_; //Make the listing of the directory listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_); list::const_iterator iFile, iFileEnd, iSpec, iSpecEnd; iFileEnd = listOfFiles.end(); StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects TMatrix meanOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); uint_4 nSpectra=0; //Loop on files for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) { if (para.debuglev_>90){ cout << "load file <" << *iFile << ">" << endl; } PInPersist fin(*iFile); vector vec = fin.GetNameTags(); list listOfSpectra; //Keep only required PPF objects std::remove_copy_if( vec.begin(), vec.end(), back_inserter(listOfSpectra), not1(match) ); listOfSpectra.sort(stringCompare); iSpecEnd = listOfSpectra.end(); //Loop of spectra matrix for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd; ++iSpec){ if (para.debuglev_>90){ cout << " spactra <" << *iSpec << ">" << endl; } TMatrix aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); fin.GetObject(aSpec,*iSpec); //How to see if the GetObject is ok?? Ask Reza nSpectra++; meanOfSpectra+=aSpec; }//eo loop on spectra in a file }//eo loop on files //Normalisation if(nSpectra>0)meanOfSpectra/=(r_4)(nSpectra); //Compute the reduced version of the mean and sigma TMatrix meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_); TMatrix sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_); reduceSpectra(meanOfSpectra,meanRedMtx,sigmaRedMtx); {//Save the result stringstream tmp; tmp << nSpectra; string fileName = para.outPath_+"/meanDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf"; cout << "Save mean based on " << nSpectra << " cycles " << endl; POutPersist fos(fileName); string tag = "mean"; fos << PPFNameTag(tag) << meanOfSpectra; tag = "meanred"; fos << PPFNameTag(tag) << meanRedMtx; tag = "sigmared"; fos << PPFNameTag(tag) << sigmaRedMtx; } } //------------------------------------------------------- //Compute the median of Diff ON-OFF Raw spectra and also the mean/sigma of rebinned spectra //Used like: // void medianRawDiffOnOffCycles() throw(string) { list listOfFiles; string directoryName; directoryName = para.inPath_ + "/" + para.sourceName_; //Make the listing of the directory listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_); list::const_iterator iFile, iFileEnd, iSpec, iSpecEnd; iFileEnd = listOfFiles.end(); TArray tableOfSpectra(NUMBER_OF_FREQ,NUMBER_OF_CHANNELS,TOTAL_NUM_CYCLES); //TOTAL_NUM_CYCLES should be large enough... StringMatch match("specONOFFRaw[0-9]+"); //Tag of the PPF objects uint_4 nSpectra=0; //Loop on files for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) { if (para.debuglev_>90){ cout << "load file <" << *iFile << ">" << endl; } PInPersist fin(*iFile); vector vec = fin.GetNameTags(); list listOfSpectra; //Keep only required PPF objects std::remove_copy_if( vec.begin(), vec.end(), back_inserter(listOfSpectra), not1(match) ); listOfSpectra.sort(stringCompare); iSpecEnd = listOfSpectra.end(); //Loop of spectra matrix for (iSpec = listOfSpectra.begin(); iSpec !=iSpecEnd && (sa_size_t)nSpectra < TOTAL_NUM_CYCLES ; ++iSpec){ if (para.debuglev_>90){ cout << " spactra <" << *iSpec << ">" << endl; } TMatrix aSpec(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); fin.GetObject(aSpec,*iSpec); tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec; nSpectra++; }//eo loop on spectra in a file }//eo loop on files TMatrix medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //Compute the median for each freq. and Channel for (sa_size_t iCh=0; iCh tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions()); vector tmp; tmp0.FillTo(tmp); medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end()); } } //Compute the reduced version of the mean and sigma TMatrix meanRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_); TMatrix sigmaRedMtx(NUMBER_OF_CHANNELS,para.nSliceInFreq_); reduceSpectra(medianOfSpectra,meanRedMtx,sigmaRedMtx); sa_size_t f1320=freqTochan(1320.); sa_size_t f1345=freqTochan(1345.); sa_size_t f1355=freqTochan(1355.); sa_size_t f1380=freqTochan(1380.); //Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380] if (para.debuglev_>9){ cout << "Compute baseline arround 1350Mhz on [1320-1345] U [1355-1380]" << endl; } TVectormeanMed(NUMBER_OF_CHANNELS); for (sa_size_t iCh=0; iCh band1; band1 = medianOfSpectra(Range(iCh),Range(f1320,f1345)).CompactAllDimensions(); MeanSigma(band1,meanMed1,sigmaMed1); double meanMed2; double sigmaMed2; TVector band2; band2 = medianOfSpectra(Range(iCh),Range(f1355,f1380)).CompactAllDimensions(); MeanSigma(band2,meanMed2,sigmaMed2); meanMed(iCh) = (meanMed1+meanMed2)*0.5; } meanMed.Print(cout); //Compute the sigma in the range 1320MHz-1380MHz if (para.debuglev_>9){ cout << "Compute the sigma in the range 1320MHz-1380MHz" << endl; } TVectorsigmaMed(NUMBER_OF_CHANNELS); sa_size_t redf1320=(sa_size_t)((1320.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_); sa_size_t redf1380=(sa_size_t)((1380.0-LOWER_FREQUENCY)/TOTAL_BANDWIDTH*para.nSliceInFreq_); for (sa_size_t iCh=0; iCh band; band = sigmaRedMtx(Range(iCh),Range(redf1320,redf1380)).CompactAllDimensions(); MeanSigma(band,meanSigma,sigmaSigma); meanSigma *= sqrt(para.nSliceInFreq_); //to scale to orignal spectra sigmaMed(iCh) = meanSigma; } sigmaMed.Print(cout); if (para.debuglev_>9){ cout << "Compute medianOfSpectraNorm" << endl; } TMatrix medianOfSpectraNorm(medianOfSpectra,false); //do not share the data... for (sa_size_t iCh=0; iCh tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1))); fos << PPFNameTag(tag) << tarr; } } //------------------------------------------------------- int main(int narg, char* arg[]) { int rc = 0; //return code string msg; //message used in Exceptions //default value for initial parameters (see Para structure on top of the file) string debuglev = "0"; string action = "meanDiffOnOff"; string inputPath = "."; string outputPath = "."; string sourceName = "Abell85"; string ppfFile; string nSliceInFreq = "32"; string typeOfCalib="perRun"; // bool okarg=false; int ka=1; while (ka\n" << " (ex. /sps/baoradio/AmasNancay/JEC/\n " << " -outPath [.]| \n" << " -nSliceInFreq [32]|\n" << " -ppfFile (ex. diffOnOffRaw)" << " -debug " << endl; return 0; } else if (strcmp(arg[ka],"-debug")==0) { debuglev=arg[ka+1]; ka+=2; } else if (strcmp(arg[ka],"-act")==0) { action=arg[ka+1]; ka+=2; } else if (strcmp(arg[ka],"-calibopt")==0) { typeOfCalib=arg[ka+1]; ka+=2; } else if (strcmp(arg[ka],"-inPath")==0) { inputPath=arg[ka+1]; ka+=2; } else if (strcmp(arg[ka],"-outPath")==0) { outputPath=arg[ka+1]; ka+=2; } else if (strcmp(arg[ka],"-src")==0) { sourceName=arg[ka+1]; ka+=2; } else if (strcmp(arg[ka],"-ppfFile")==0) { ppfFile=arg[ka+1]; ka+=2; } else if (strcmp(arg[ka],"-nSliceInFreq")==0) { nSliceInFreq=arg[ka+1]; ka+=2; } else ka++; }//eo while para.debuglev_ = atoi(debuglev.c_str()); para.inPath_ = inputPath; para.outPath_ = outputPath; para.sourceName_ = sourceName; para.ppfFile_ = ppfFile; para.nSliceInFreq_ = atoi(nSliceInFreq.c_str()); cout << "Dump Initial parameters ............" << endl; cout << " action = " << action << "\n" << " inputPath = " << para.inPath_ << "\n" << " outputPath = " << para.outPath_ << "\n" << " sourceName = " << para.sourceName_ << "\n" << " ppfFile = " << para.ppfFile_ << "\n" << " nSliceInFreq = " << para.nSliceInFreq_ << "\n" << " debuglev = " << para.debuglev_ << "\n"; cout << "...................................." << endl; if ( "" == ppfFile ) { cerr << "mergeAnaFiles.cc: you have forgotten ppfFile option" << endl; return 999; } try { // int b,e; // char *match=regexp("truc0machin","[a-z]+[0-9]*",&b,&e); // printf("->%s<-\n(b=%d e=%d)\n",match,b,e); if ( action == "meanRawDiffOnOff" ) { meanRawDiffOnOffCycles(); } else if (action == "medianRawDiffOnOff") { medianRawDiffOnOffCycles(); } else if (action == "meanCalibBAODiffOnOff") { meanCalibBAODiffOnOffCycles(); } else { msg = "Unknown action " + action; throw msg; } } catch (std::exception& sex) { cerr << "mergeRawOnOff.cc std::exception :" << (string)typeid(sex).name() << "\n msg= " << sex.what() << endl; rc = 78; } catch ( string str ) { cerr << "mergeRawOnOff.cc Exception raised: " << str << endl; } catch (...) { cerr << "mergeRawOnOff.cc catched unknown (...) exception " << endl; rc = 79; } return 0; }